% IMPORTANT NOTE:
% - figures available on request from the authors

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \documentstyle[12pt,agums]{article}
 \documentstyle[agupp]{article}
% \documentstyle[jgrga]{article}

% \tighten

\lefthead{M\'ETRIS ET~AL.}
\righthead{ LAGEOS ECCENTRICITY EXCITATIONS}
%\received{September 1, 1995}
%\revised{February 12, 1995}
%\accepted{February 13, 1995}
%\journalid{JGRB}{November 1995}
%\articleid{1}{13}
%\paperid{92JB04601}
%\ccc{0148-0227/93/92JB-04601\$05.00}
   % \cpright{PD}{1993}
   % \cpright{Crown}{1993}
   % (Crown copyrights have no "\ccc{}" information.)
%\cpright{AGU}{1993}

\authoraddr{G. M\'etris, OCA/CERGA, Av N Copernic, F-06130 Grasse,
 France (e-mail: metris@ocar01.obs-azur.fr).}
\authoraddr{D. Vokrouhlick\'y, Institute of Astronomy, Charles
 University,\v Sv\'edsk\'a 8, CZ-15000 Prague 5, Czech Republic
 (e-mail: vokrouhl@earn.cvut.cz).}
\authoraddr{ J.C. Ries, Center for Space Research 60605, The
 University of Texas at Austin, Austin, TX 78712, USA
 (e-mail: ries@moebius.csr.utexas.edu).}
\authoraddr{ R.J. Eanes, Center for Space Research 60605, The
 University of Texas at Austin, Austin, TX 78712, USA
 (e-mail: eanes@csr.utexas.edu).}

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

\begin{document}

\title{Non-gravitational phenomena \\
 and the LAGEOS eccentricity excitations}

\author{Gilles M\'etris and David Vokrouhlick\'y\altaffilmark{1}}
\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, The Czech Republic.}
\author{John C. Ries and Richard J. Eanes}
\affil{ Center for Space Research, The University of Texas at Austin,
 Austin, Texas, USA}

\begin{abstract}
 A model for explaining the anomalous eccentricity excitations of
 the LAGEOS~I satellite orbit is considered. It is suggested
 that two phenomena are responsible for the major part of these residuals:
 (i) about $1.2$ \% mismodeling of the LAGEOS~I surface mean radiation
 pressure coefficient, and (ii) the Yarkovsky-Schach thermal effect
 with properly modeled evolution of the satellite spin axis.  Minor
 influence is also attributed to the asymmetric reflectivity of the
 satellite surface and the asymmetric thermal emissivity of the Earth.
 In all cases we give a suitable and compact formulation for the
 long-term orbit analysis. It is argued that contrary to the case
 of the LAGEOS~I along-track residuals adjustment, which does not
 allow the determination of a particular combination of the
 Yarkovsky-Schach effect parameters (being insensitive to it), the
 combined analysis of the along-track residuals and the eccentricity
 excitations allows for their decorrelation. The
 resulting value of the Yarkovsky-Schach amplitude is estimated to be
 approximately $175\;{\rm pm/s}^2$, about twice as large as
 given in the preceding works.
 It is also shown that the effect of the radiation
 pressure of the Earth reflected sunlight (the albedo effect),
 can be modeled with surprising precision by the redefinition of
 the satellite surface reflectivity coefficient.
 In view of the unexplained part of the eccentricity residuals, we
 discuss the possible importance of several other non-gravitational
 effects.
\end{abstract}

\section{Introduction}
Non-gravitational phenomena perturbing the artificial satellite motion
attracted considerable attention in the past years. Part of the cause
is that the precision of the satellite orbit determination continually
increases because of the improved
tracking technologies and increasing density of high quality
tracking stations. As a result, more and more subtle dynamical
effects come to the evidence. A great part of them is (directly
or indirectly) of the non-gravitational origin. Due to the tiny
values of the perturbing accelerations, such phenomena cannot
be determined via the short-arc considerations. Only the resulting
long-term phenomena can be detected and analyzed. Thus the final,
however equally important, condition which raised interest in
studying non-gravitational phenomena is that we have
sufficiently long-period spanning high-quality satellite
orbit data for several satellites dedicated to geodynamics
applications. Among them, the LAGEOS~I orbit analysis would be
of primary interest for us.

Up to now, a considerable effort has been devoted to the
analysis of the non-gravitational perturbations in connection with
the LAGEOS~I semimajor axis (or equivalently along-track force)
residuals [e.g. \markcite{{\it Afonso et~al.,} 1985}; \markcite{{\it
Barlier et~al.,} 1986}; \markcite{{\it Rubincam,} 1987, 1988};
\markcite{{\it Scharroo et~al.,} 1991}; \markcite{{\it Vokrouhlick\'y
and Farinella,} 1995a,b}; \markcite{{\it Farinella et~al.,} 1995}].
The mean semimajor axis is the most precisely defined
orbital element and, thus, it allows studying the influence of tiny
dynamical effects. Nonetheless, at the beginning of
nineties, arguments for simultaneous (global) solution
of all orbital element residuals has been strengthened.
It appears that the major part of the semimajor axis
residuals are currently believed to be understood and the most
unexplained dynamical features act on the other orbital elements --
notably components of the eccentricity vector $(e\cos\omega,
e\sin\omega)$ [e.g. \markcite{{\it Tapley et~al.,} 1993}]. Here,
$e$ is the eccentricity and $\omega$ the argument of perigee measured
from the osculating node.

Only a minor attention has been dedicated to the analysis of the
non-gravitational perturbations of the other mean
orbital elements. Most of it has been pointed to the nodal
residuals motivated by the LAGEOS~III mission proposal and
analysis of the potential bias of the Lense-Thirring
effect detection [e.g. \markcite{{\it Farinella et~al.,} 1990};
\markcite{{\it Watkins et~al.,} 1990}; \markcite{{\it
Ciufolini et~al.,} 1993}; \markcite{{\it Vokrouhlick\'y
et~al.,} 1993}]. Broader context has been considered by
\markcite{{\it Farinella et~al.} [1990]}, who studied the
influence of the thermal
effects on the inclination and node excitations (although
they did not apply their model directly to the orbital data
analysis). Interestingly, they did not pursue the work to the
eccentricity excitations analysis which might have much earlier
shifted general attention to the non-gravitational effects in
discussion of the eccentricity residuals.

Another exception among the previous studies which focused on
the subject of the LAGEOS~I eccentricity residuals in context
of the non-gravitational influence is a recent work by
\markcite{{\it Martin and Rubincam} [1995]}, who studied
perturbations due to the sunlight reflected on the Earth surface
(hereafter called the albedo effect\footnote{ In this paper,
we reserve the meaning of the ``albedo effect'' to the diffusely
reflected part of the sunlight on the Earth surface. We
verified, that the specularly reflected radiation on the
oceans, which is dominantly responsible for the
albedo perturbations of the LAGEOS~I semimajor axis
[\markcite{{\it Vokrouhlick\'y and Farinella,} 1995a,b}],
does not contribute importantly to the eccentricity residuals.}).
Based on detailed data of the ERBE project, {\it Martin and Rubincam}
were able to model very intricate features of the complex phenomenon
enclosed generally by the albedo concept. To their surprise,
they found a strong anticorrelation of the results with the
observed eccentricity excitations concluding that such output
just enhances the urgency to focus on the subject of the
LAGEOS~I eccentricity residuals/excitations. Inspired by these
facts we demonstrate below two important theorems related to
the albedo effect influence on the LAGEOS~I eccentricity
excitations: (i) despite of complexity of the albedo distribution
on the Earth surface, it is merely the mean albedo value (with
subdominant participation of the quadrupole term in the albedo
latitude development) which determines the eccentricity signal,
and (ii) the eccentricity excitation due to the constant
albedo Earth model is equivalent, with surprising degree of
precision, to a small decrease of the radiation pressure coefficient.
Thanks to these results we are able to give simple analytical
formulation of the long-term eccentricity albedo signal.

The aim of this paper is to point out several non-gravitational
effects which may be linked with the observed eccentricity
excitation of the LAGEOS~I satellite. In view of the previous
notes, we have to satisfy one elementary rule: the proposed
effects should not create a conflict with the fitted model
of the LAGEOS~I along-track residuals or, at least, in a carefully
controlled manner. A trivial way to
fulfil this requirement is to consider the non-gravitational
phenomena contributing to the along-track signal itself (if not done
before). Otherwise, one should carefully respect ``the razor''
yielded by the successful model for the along-track residuals.

Considering the guidelines mentioned above, we first investigate the
influence of mismodeling in the LAGEOS~I radiation pressure coefficient
$C_R$. Since the principal action of the direct radiation pressure
is focused on the eccentricity itself, we face a problem of
explaining why the (long-arc)-adjusted solution of $C_R$ for
LAGEOS~I does leave a residual eccentricity signal. Below, we propose
some routes which might indicate such kind of bias: the albedo
effect mismodeling, an intricate interplay between the radiation
pressure effect and other non-gravitational phenomena during the
phases when the LAGEOS~I orbit enters the Earth shadow, etc.. In frame
of our suggestion, about $1.2$~\% increase of the mean (effective)
radiation pressure coefficient may cure an important part of the
eccentricity residuals.

Secondly, we discuss in some detail the possible role of the thermal
effects and the LAGEOS~I surface asymmetric reflectivity effect in
context of the observed eccentricity residuals (thus the natural
candidates to be studied as such phenomena are standardly employed
for interpreting the along-track residuals of the LAGEOS~I orbit).
It might seem to be surprising that no published attempt
to elaborate such study has been performed until now (to the
authors knowledge), despite of the fact that those phenomena
have been seriously considered in the studies of the along-track
residuals since the end of the eighties. Probably one
important circumstance which inhibited extension of the
non-gravitational force influence on all orbital elements of
LAGEOS~I is a lack of a firm knowledge of the satellite spin axis
evolution. Several important effects -- notably the thermal phenomena
and the optical anisotropy effect -- depend critically on the
spin axis orientation. A lack of the theoretical
understanding of this evolution has posed (and still poses) generally
an important obstacle in the non-gravitational force modeling
progress. Alternatively, several authors has attempted to
solve inversely for the LAGEOS~I spin axis evolution via analysis
of the thermal phenomena perturbation of the satellite
mean semimajor axis or, equivalently, the orbit along-track residuals
[\markcite{{\it Ries et~al.,} 1993}; \markcite{{\it Robbins et~al.,}
1994}; \markcite{{\it Kolenkiewicz et~al.,} 1995}]. Such procedure
is in danger due to the mismodeling of the
phenomena acting on the semimajor-axis residuals, a phenomenon
which can seriously pollute the resulting axis evolution [e.g.
\markcite{{\it Robbins et~al.,} 1994}]. Obviously, even greater danger
of mismodeling would arise if one attempted to fit the spin axis
evolution using residuals of the other orbital element
residuals.

A conceptually better way of tackling
the problem is to involve a theoretical model for the satellite
spin axis evolution. As commented by \markcite{{\it Farinella et~al.}
[1995]}, curious misunderstanding has gradually led to a critical
attitude of the specialists towards the model developed by
\markcite{{\it Bertotti and Iess} [1991]}. Only recently,
\markcite{{\it Farinella et~al.} [1995]}
have clarified that the original model is able to recover the
essential features of the spin axis evolution, even though the
fine details of the evolution remain difficult to be explained
with relatively simple model applied both by \markcite{{\it Bertotti
and Iess} [1991]} and \markcite{{\it Farinella et~al.} [1995]}.
Namely, in the last
period, when the LAGEOS~I rotation period has been considerably
decreased by the action of the eddy current dissipation,
one still observes moderate differences between the theoretical
modeling and ``empirical'' axis position [\markcite{{\it Farinella
et~al.,} 1995}]. Hopefully, refined theoretical analysis will possess
improvements of this point in the future.

In any case, availability of the independent spin axis evolution
modeling allows us to decrease the danger of bias when
considering the influence of the thermal phenomena and the optical
anisotropy effect on {\it all} orbital elements. Hereafter,
we shall focus on the analysis of the eccentricity excitations.

In addition to previous phenomena, we also comment on the
influence of the north/south asymmetry of the Earth thermal
emissivity. Though difficult to model, this phenomenon is
naturally present in Earth radiation data. On contrary to the
monopole (constant) and quadrupole parts of the Earth thermal
emissivity development which are ``north/south symmetric'',
it involves the ``north/south asymmetry'' which results in
significant excitation of the satellite's eccentricity.

The paper is organized as follows: in the next section we comment
on modeling of the non-gravitational phenomena considered in this
paper and derive their long-term excitation of the eccentricity
vector components. Next, we proceed with applying previous considerations
to the case of LAGEOS~I and discuss the comparison of the
resulting signal with the observed orbit excitations. Finally, we also
pay some attention to other non-gravitational effects
and their possible influence on the LAGEOS~I eccentricity
residuals, together with a brief conclusion. Technical details of
the (thermal) radiative force calculation are relegated to
Appendix~A.

\section{ Modeling and development of the perturbations}
We reserve this section for discussion of the phenomena which
we basically propose as a good candidates for explaining of the
LAGEOS~I eccentricity excitations (discussion of the less
favorable candidates is postponed to the last section). They are
namely:
\begin{itemize}
\item the effect of small change of the global (long-arc-adjusted)
 value of the radiation pressure coefficient;
\item the Yarkovsky-Schach thermal effect;
\item the effect of optical asymmetry of the two satellite
 hemispheres;
\item the diffuse part of the albedo effect;
\item the north/south asymmetry of the Earth thermal
 emissivity.
\end{itemize}

\noindent In all cases we comment on the basic formulation
of the phenomenon and then derive its long-term excitation of
the eccentricity vector components.

\subsection{ Force modeling}
The first proposed phenomenon consist of redefining the adjusted
effective value of LAGEOS~I surface reflectivity coefficient $C_R^0$
by a small fraction $\Delta C_R$. An obvious output of such
proposal is a supplemental perturbation force oriented oppositely
to the solar direction (specified by unit vector ${\bf n}$)
\begin{equation}
 {\bf a} = -{\pi I_0 R_s^2\over mc} \Delta C_R \,
 \Gamma(f)\, {\bf n} \; , \label{eq1}
\end{equation}
where $I_0$ is the solar constant (effectively entailing ``its
geometric variations'' due to the Earth orbit eccentricity), $R_s$
satellite radius, $m$ satellite mass, $c$ the light velocity.

The factor $\Gamma(f)$ ($f$ being the true anomaly) in Eq.~(\ref{eq1})
represents the shadow function
well-known from the analysis of the direct solar radiation pressure
perturbations [e.g. \markcite{{\it L\'ala and Sehnal,} 1969};
\markcite{{\it Ferraz-Mello,} 1972}]. Its value reads $0$ inside
the geometrical Earth shadow and $1$ outside.
Developing this function in a Fourier series $\Gamma(f)=
\sum_{k\in{\bf Z}} c_k \exp(-\imath kf)$, we obtain\footnote{
Notice, that the following formulation rigorously applies to the
case when the shadow position (as well as its duration) within one
revolution of the satellite is constant. Only then the conditions
of periodicity necessary for application of the Fourier development
are satisfied. However, satellite node and perigee motion as well
as the solar motion continually change the mutual Sun-orbit
configuration. Especially at the beginning/end of the shadow
crossing periods our approach may thus result in some mismodeling.
Discussion of such details is out of the scope of this paper.}
\begin{mathletters}
 \label{eq1bis}
\begin{eqnarray}
 c_0 &=& \;1 - {f_2-f_1 \over 2\pi} \; , \\
 c_k &=&\; {1 \over 2\pi \imath k} \left[ \exp\left(\imath k f_1
  \right)-\exp\left(\imath k f_2\right)\right] \; ,
\end{eqnarray}
\end{mathletters}
($\imath$ is the complex unit) where the second row applies to all
$k\neq 0$. The angles $f_1$ and $f_2$ in the previous formulas denote
the true anomaly of the Earth shadow entry and exit, respectively,
along one satellite revolution around the Earth.

Next, let us focus on the effect of the asymmetric reflectivity of
the satellite surface. Even
that it rests on the level of hypothesis as no such asymmetry
has been detected during the prelaunch tests ({\it L. Anselmo},
private communication), it appears that its existence has been
convincingly confirmed from the ``practical point of view'' via the
analysis of the LAGEOS~I along-track residuals [e.g.
\markcite{{\it Scharroo et~al.,} 1991}; \markcite{{\it Farinella
et~al.,} 1995}]. It has been also recently
suggested, that the origin
of the asymmetry might be during the last stage of the satellite
boosting due to the surface blemishing by the launcher. This
idea might also well correspond to the fact that the two satellite
hemispheres, which asymmetry is discussed, are defined by the
orientation of the spin axis. The instantaneous perturbing
acceleration associated with this effect is given by
\begin{equation}
 {\bf a} = \gamma \Gamma(f) \left[\zeta\left(\theta_r\right){\bf s}
 + 2\cos\theta_r {\bf n}\right] \; , \label{eq2}
\end{equation}
where
\begin{mathletters}
 \label{eq3}
\begin{eqnarray}
 \gamma &=& -{\pi I_0 R_s^2 \over 4 mc} \delta \; , \label{eq3a} \\
 \zeta(\theta_r) &=& \sin^2\theta_r \; , \label{eq3b} \\
 \delta &=& 2\, {C^+-C^-\over C^+ +C^-}\; .  \label{eq3c}
\end{eqnarray}
\end{mathletters}
Here, ${\bf s}$ represents unit vector in the direction of
the satellite spin axis and ${\bf n}$ unit vector in the solar
direction as before, $\theta_r$ mutual angle of the two vectors (i.e.
$\cos\theta_r = {\bf s}.{\bf n}$), and $C^+$ and $C^-$ are the
reflectivity coefficients of the ``northern'' and ``southern''
satellite hemispheres, respectively. It is appropriate to make
several comments related to the modeling of the anisotropic
surface reflection effect since,
to our knowledge, the effect has not been discussed in sufficient
detail in the literature so far. First, it would be mistake
to present the effect as an effective increase of the radiation
reflectivity coefficient on one of the satellite hemispheres while
keeping its nominal value on the other. By doing so, one
introduces apart of the ``true anisotropy effect'' an effective
increase of the mean value of the radiation pressure coefficient.
Hence, speaking about the anisotropic reflectivity effect one
is allowed to change the satellite surface reflectivity coefficient
in an asymmetric way but keeping its mean value unchanged.
One can for instance imagine increasing the ``northern satellite
hemisphere reflectivity'' by a fraction $+{\delta/2}$ of the
mean reflectivity value $(C^+ + C^-)/2$, while decreasing
the ``southern satellite hemisphere reflectivity'' by a fraction
$-{\delta/2}$. Such formulation fits the spirit of the model
presented by \markcite{{\it Scharroo et~al.,} [1991]} and results in
the effective perturbation force as given in Eqs.~(\ref{eq2}) and
(\ref{eq3}). Then, $C^+$ and $C^-$ are the optical coefficients of
reflectivity of the northern and the southern LAGEOS~I hemispheres
(defined by the spin axis orientation) and the parameter
$\delta$ is the degree of optical asymmetry of the two
hemispheres. As mentioned previously, this parameter is not known a
priori, but has been adjusted by the along-track residual analysis to the
value of about $1.4$~\% [e.g. \markcite{{\it Scharroo et~al.,} 1991};
\markcite{{\it Farinella et~al.,} 1995}].

A second important note to be mentioned is that
the anisotropic reflectivity force (\ref{eq2}) obeys two
components: one in the direction of the satellite spin axis ${\bf s}$
the second in the direction of the Sun ${\bf n}$. The latter
component has been totally omitted by \markcite{{\it Scharroo et~al.}
[1991]} as, indeed, it does not contribute to the long term evolution
of the semimajor axis treated by those authors. However, it would be a
mistake to neglect this part when considering the excitations of the
other satellite elements, particularly that of eccentricity.

In view of the previously mentioned hypothesis about the satellite
surface anisotropy origin mentioned before, we tried to refine
our model for the anisotropic reflectivity by supposing another
``latitudinal'' distributions of the reflectivity coefficient on
the LAGEOS~I surface\footnote{ In this way, we can also test the
robustness of the anisotropic effect formulation.}. Indeed, it
does seem artificial to suppose that the optical
properties of the two LAGEOS~I hemispheres are precisely
uniform. We thus consider models which account for distribution of
the reflectivity coefficient proportional to an arbitrary {\it odd} power
($n$, say) of the ``the satellite surface colatitude'' (thus zero on
the LAGEOS~I ``equator'' and maximal on its ``poles''). Surprisingly,
the resulting formulas for the perturbative acceleration rest
simple and consist merely in an appropriate redefinition of the
$\gamma$ factor and $\zeta(\theta_r)$ function in (\ref{eq2}) by
\begin{mathletters}
 \label{eq4}
\begin{eqnarray}
 \gamma &=& -{4\over (n+2)(n+4)} {\pi I_0 R_s^2 \over mc} \delta \; ,
  \label{eq4a} \\
 \zeta(\theta_r) &=& 1+(n-1)\cos^2\theta_r \; . \label{eq4b}
\end{eqnarray}
\end{mathletters}
Parameter $\delta$ loses its simple interpretation (\ref{eq3c}) and
is to be considered as an free parameter of the problem. Naturally,
formulas (\ref{eq3}) might be obtained as a cumulative sum of
particularly weight contributions of (\ref{eq4}). [Remember that
the step function having
the values of $\pm 1$ on the northern and southern hemispheres,
respectively, can be decomposed into an infinite series of the
odd powers of cosine of the colatitude.] Generally, results obtained either
with the simple formulas (\ref{eq3}) or a given value of $n$ in formulas
(\ref{eq4}) are very similar. We shall thus essentially employ formulation
of the anisotropic reflectivity given by Eqs.~(\ref{eq2}) and
(\ref{eq3}), and only briefly demonstrate closeness of the
formulation following formulas (\ref{eq4}).

For the Yarkovsky-Schach effect, related to the recoil force
associated to anisotropic thermal radiation emitted from the satellite
surface asymmetrically heated by the solar radiation, we shall closely
follow the formulation by \markcite{{\it Afonso et~al.} [1989]}. The
instantaneous perturbing acceleration reads
\begin{equation}
 {\bf a} = {\tilde \gamma} \Gamma_{\rm ys}(f) \cos\theta_r
 {\bf s} \; , \label{eq5}
\end{equation}
where $\Gamma_{\rm ys}(f)$ stands for the ``modified
shadow function'' [also the ``physical shadow function''; e.g.
\markcite{{\it Afonso et~al.,} 1989}] reading
\begin{mathletters}
 \label{eq6}
\begin{eqnarray}
 \Gamma_{\rm ys}(f) &=& C \exp\left[\left(f_1 - f\right)
 /f_0\right]\; , \label{eq6a} \\
  &=& C C' \left\{1 - \exp\left[\left(f - f_2\right)/f_0\right]\right\}
   \nonumber \\ & & + C \exp\left[\left(f_2 - f_1\right)/f_0\right]
   . \label{eq6b}  \\
  C &\equiv& {1 - \exp\left[\left(f_2 - f_1 - 2\pi\right)/f_0\right]
   \over 1 - \exp\left[-2\pi/f_0\right]} \; ,\label{eq6c} \\
  C' &\equiv& {1 - \exp\left[\left(f_1 - f_2\right)/f_0\right]\over
   1 - \exp\left[\left(f_2 - f_1 - 2\pi\right)/f_0\right]}\; ,
   \label{eq6d}
\end{eqnarray}
\end{mathletters}
The first line -- (\ref{eq6a}) -- applies for the phase when the
satellite is in the geometrical shadow [$f \in (f_1,f_2)$], while
the second one -- (\ref{eq6b}) -- applies out of the geometrical shadow
[$f \notin (f_1,f_2)$]. The free parameter $f_0$ results from the
thermal lag between the solar heating and equilibrium temperature
partition in the satellite body [see \markcite{{\it Rubincam} 1987,
1988}; \markcite{{\it Afonso et~al.,} 1989}]. Though it is possible
to derive a theoretical model to estimate the value of this thermal
lag parameter, it is usually taken as an empirical and adjusted
coefficient due to the uncertainty of such models [e.g. \markcite{{\it
Afonso et~al.,} 1989}, \markcite{{\it Scharroo et~al.,} 1991};
\markcite{{\it Farinella et~al.,} 1995}]. Notice that in comparison
with \markcite{{\it Scharroo et~al.} [1991]} we do not neglect the
factor $\exp(-2\pi/f_0)$ in the description of the physical shadow
function $\Gamma_{\rm ys}(f)$ which allows us to consider any value
of the thermal lag parameter. Similarly, the
amplitude ${\tilde \gamma}$ of the perturbing force (\ref{eq5}) might
be theoretically estimated from the satellite body properties
[see \markcite{{\it Rubincam} 1987, 1988}; \markcite{{\it Afonso
et~al.,} 1989}]. However, the uncertainty in the equilibrium temperature of
the satellite and the temperature difference of its two
hemispheres estimates as well as the complications due to the
inhomogeneous and anisotropic structure of LAGEOS~I result in
the common practise of adjusting the value from the along-track
residuals. It is important to note that the along-track
adjustment does not allow decorrelation of the two principal
parameters of the Yakovsky-Schach effect as it is nearly insensitive
to their particular (nearly linear) combination [\markcite{{\it
Farinella et~al.,} 1995}]. In the following, we shall see that the
combined analysis of the along-track residuals and the
eccentricity excitations might allow decorrelation
of the two parameters.

Performing again the Fourier series analysis of the modified
shadow function $\Gamma_{\rm ys}(f)= \sum_{k\in{\bf Z}} {\tilde c}_k
\exp(-\imath kf)$ we obtain
\begin{mathletters}
 \label{eq7}
\begin{eqnarray}
 {\tilde c}_0 &=& 1 - {f_2-f_1 \over 2\pi}  \; , \label{eq7a} \\
 {\tilde c}_k &=& {1 \over 2\pi k(\imath +kf_0)} \left[
  \exp\left(\imath k f_1 \right)-\exp\left(\imath k f_2\right)
  \right] \; . \nonumber \\ & & \label{eq7b}
\end{eqnarray}
\end{mathletters}
It is a remarkable that the complex factors $C$ and $C'$ from the
definition (\ref{eq6}) disappeared exactly from the formulas
(\ref{eq7}) giving the Fourier coefficients $({\tilde c}_0,
{\tilde c}_k)$.

In the late stage of working on this paper we have learned
about the detailed numerical thermal model of the LAGEOS
satellite by \markcite{{\it Slabinski} [1995]}. The output
of this model is a unified formulation of the thermal
dynamical effects, both due to the Earth IR heating (the
Rubincam effect) and the solar heating (the Yarkovsky-Schach
effect), developing the resulting force applied to LAGEOS~I
satellite in Fourier series of the satellite longitude
in orbit. Numerical coefficients of these series emerge from
the theoretical considerations, however, in order to
improve the fit of the observed signal one is allowed to
adjust some of them (note that \markcite{{\it Slabinski} [1995]}
proceeds in this way when applying his theory to the LAGEOS
along-track residuals).
We shall not repeat here extensive formulation of this model,
nevertheless we shall compare results based on using
{\it Slabinski's} model with our previous, more standard,
formulation. This comparison will yield some hints in
support of our empirically determined values for the
Yarkovsky-Schach effect parameters.

Motivated by the conclusions of \markcite{{\it Martin and
Rubincam} [1995]} we decided to revisit in more detail the
albedo effect influence on the LAGEOS~I eccentricity excitations.
Our strategy is as follows: (i) first, we develop highly simplified
model of the albedo effect, which allows to calculate
{\it analytically} the long-term eccentricity signal, (ii)
second, we justify the mathematical validity of these results by performing
complete numerical integration of general albedo problem
(allowing thus for any latitude and longitude development
of the albedo coefficient), and (iii) last, we check the pertinence of the
physical model by comparing our results to those of \markcite{{\it Martin and
Rubincam} [1995]} based on observed albedo data. In the following
paragraphs we formulate the simplified albedo problem -- (i) --
relegating the numerical -- (ii) -- and observational -- (iii) --
confirmations to the next section.

\markcite{{\it Borderies and Longaretti} [1990]} succeeded to
reduce the double (2-D) integration necessary for evaluation of the
albedo radiation pressure to a single (1-D) one for a special
case of the constant albedo model, giving the components of
the satellite acceleration projected onto a local orbital
frame $({\bf r}, {\bf t}, {\bf w})$ [${\bf r}$ being the
unit vector in the radial direction, ${\bf t}$ the unit vector
in the transverse direction, and ${\bf w}$ the unit vector normal
to the orbit] as
\begin{mathletters}
 \label{eq1a}
\begin{eqnarray}
 {\bf a}_{\rm alb} &=& R\,{\bf r} + T\,{\bf t} + W\,{\bf w}\; ,
  \label{eq1ala}  \\
 R &=& {R_s^2 I_0\over mc} C_R \,{\bar{\cal A}} \alpha^2 J_{x'}\left(
 \sigma,\alpha\right) \; , \label{eq1alb} \\
 T &=& - {R_s^2 I_0\over mc} C_R \,{({\bf n}.{\bf t})\over \sin
  \sigma}\, {\bar{\cal A}} \alpha^3 J_{y'}\left(\sigma,\alpha\right)
  \; , \label{eq1alc} \\
 W &=& - {R_s^2 I_0\over mc} C_R\, {({\bf n}.{\bf w})\over \sin
  \sigma}\, {\bar{\cal A}} \alpha^3 J_{y'}\left(\sigma,\alpha\right)
  \; , \label{eq1ald}
\end{eqnarray}
\end{mathletters}
where ${\bar{\cal A}}$ denotes the value of the Earth mean albedo,
$\alpha = R_e/a$ with $R_e$ denoting the Earth radius and $a$
the satellite semimajor axis, $\sigma$ is the geocentric angular
distance of the Sun and the satellite. The functions $J_{x'}(\sigma,
\alpha)$ and $J_{y'}(\sigma,\alpha)$ are given for all possible
Sun-satellite mutual configurations by \markcite{{\it Borderies
and Longaretti} [1990]} in form of 1-D integrals. Starting with
these formulas, we adopt the following simplifications: (i) in all
components $(R,T,W)$ we retain only the dominant power
of the satellite ``parallax'' $\alpha$, and (ii) we
disregard the transitional phase of the satellite orbit
when the Earth surface region visible from the satellite
contains the terminator of the solar illumination. Under these
simplifying assumptions the albedo acceleration of the satellite
becomes
\begin{eqnarray}
 {\bf a}_{\rm alb} &=& {\pi R_s^2 I_0\over 4mc} C_R\, {\bar{\cal A}}
 \alpha^2\, \Gamma_{\rm alb}\left(f\right)\times \nonumber \\
 & & \;\; \left[\cos\sigma \left(\case83 +3\alpha\right)\,
  {\bf r} - \alpha\, {\bf n}\right] \; , \label{eq2a}
\end{eqnarray}
where $\Gamma_{\rm alb}(f)$ is ``the albedo shadow function''
equal to $1$ for $\sigma \leq \pi/2$ and $0$ for $\sigma >
\pi/2$. Similarly to the previously considered shadow functions,
$\Gamma_{\rm alb}(f)$ can be represented by a Fourier series
of the true anomaly $f$ with coefficients equal to those
listed in (\ref{eq1bis}) provided the ``shadow'' limit
angles $f_1$ and $f_2$ are appropriately redefined. Because of
a particular symmetry of the albedo shadow definition
the even terms of its Fourier development vanish. This
property importantly simplifies the evaluation of the long-term
eccentricity excitations.

Our approximation of the albedo perturbation thus
accounts for two components, the first being radial while the
second constantly aligned with the direct solar radiation
pressure. Both of them are nonzero only over one half of the
satellite orbit [apart of the singular case of the orbital
plane perpendicular to the solar direction]. Notice, that
the mean albedo value ${\bar{\cal A}}$ of our approximation
(\ref{eq2a}) may depend on time provided the time scale of
its changes are longer than one satellite revolution. Here,
we have essentially in mind seasonal variations clearly
present in the photometric albedo data [e.g. \markcite{{
\it Stephens et~al.,} 1981}]. Although very simple, we shall
demonstrate that this formulation accounts with surprising
precision for the albedo contribution to the eccentricity
excitations.

Finally, we shall treat in some detail the effect of the north/south
thermal (IR) emissivity of the Earth. Obviously, the monopole and
quadrupole parts of the Earth thermal emissivity are dominant as
concerns the magnitude, but one
easily sees that the constant emissivity term contributes
mainly as an effective modification of the gravitational constant, while
the quadrupole term does not yield significant eccentricity
perturbations due to presence of the satellite eccentricity factor
[\markcite{{\it Sehnal,} 1981}]. It is only the asymmetry term
which could
raise perturbations of the zeroth order in eccentricity. Hence,
representing the Earth thermal emissivity ${\cal E}$ in a simple
colatitude $\vartheta$ dependent model
\begin{equation}
 {\cal E} = {\cal E}_0 + {\cal E}_1 \cos\vartheta + {\cal E}_2
 (3\cos^2\vartheta-1)/2 + \ldots \; , \label{eq8}
\end{equation}
we shall focus on the influence of the ${\cal E}_1$--associated
term. Simple algebra shows that the values of the
radial ($R$) and transverse ($T$) accelerations due to the discussed
effect are (see also Appendix~A)
\begin{mathletters}
 \label{eq9}
\begin{eqnarray}
 R &= & \,2{\hat \gamma} \psi_1(a) \sin I \sin\lambda \; ,\label{eq9a}
 \\
 T &= & -{\hat \gamma} \left[\psi_1\left(a\right) - \psi_2\left(a\right)
 \right] \sin I \cos\lambda \; ,\label{eq9b}
\end{eqnarray}
\end{mathletters}
where
\begin{equation}
 {\hat \gamma} = {\pi R_s^2 \over 4mc} {a \over R_e} {\cal C}_R
 {\cal E}_1 \label{eq10}
\end{equation}
and functions $\psi_1(a)$, $\psi_2(a)$ are given by
\begin{mathletters}
 \label{eq11}
\begin{eqnarray}
 \psi_1(a) &=& \alpha^4 + {1 \over 8} \beta^4 \Bigl\{
 \sinh\left[ 4\, {\rm arccosh}\left(1/\beta\right)\right] \nonumber \\
 & & \quad\quad\quad - 4\, {\rm arccosh}\left(1/\beta\right)\Bigr\}
 \; , \label{eq11a}\\
 \psi_2(a) &=& \beta^2 \Bigl\{\sinh\left[ 2\, {\rm arccosh}\left(
 1/\beta\right)\right] \nonumber \\
 & & \qquad\quad - 2\, {\rm arccosh}\left(1/\beta\right)\Bigr\}
 \; . \label{eq11b}
\end{eqnarray}
\end{mathletters}
Here we abbreviate $\beta = \sqrt{1-\alpha^2}$ ($\alpha$ as before).
The angle $\lambda$ in Eqs.~(\ref{eq9}) denotes the mean
satellite longitude in orbit measured from the node. The force components
(\ref{eq9}) are given up to the terms depending on the satellite
eccentricity.

\noindent We can also approximate
\begin{equation}
 \psi_2(a) \simeq {4\over 3} \left(1-\alpha^2\right)
 {\rm arccosh}^3\left[1/\sqrt{1-\alpha^2}\right] \; , \label{eq12}
\end{equation}
with a few percent error at the LAGEOS~I altitude.

\subsection{ Long-term perturbations of the eccentricity vector}
As indicated in the introduction, we are merely interested in the
long-term phenomena associated with the specified  perturbations. Thus,
in the following discussion, we shall involve the averaging
scheme instead of treating the instantaneous influence of the
perturbing forces. In terms of the classical tools of
celestial mechanics, we shall derive a first order
perturbation theory for the satellite mean elements which
is well justified considering the small magnitude
of the effects. Technically, we shall simply average the right hand
sides of the perturbation equations over the unique fast
variable in our formulation, the satellite's mean anomaly.

Perturbation equations for the components of the eccentricity vector
$(h = e\cos\omega,k=e\sin\omega)$ are given by [e.g. \markcite{{\it
Milani et~al.,} 1987}]
%\begin{mathletters}
% \label{eq13}
%\begin{eqnarray}
% {\Delta\dot h} &=& {1 \over na} \left[ R \sin\left(
% \omega + f\right) + 2T \cos\left(\omega + f\right) \right]
% \nonumber \\ & & \; + {\cal O}(e) \; ,  \label{eq13a} \\
% {\Delta\dot k} &=& {1 \over na} \left[ -R \cos\left(
% \omega + f\right) + 2T \sin\left(\omega + f\right) \right]
% \nonumber \\ & & \; + {\cal O}(e) \; ,  \label{eq13b}
%\end{eqnarray}
%\end{mathletters}
\begin{equation}
 {\Delta\dot h} + \imath {\Delta\dot k} = {1 \over na} \left[2T-
 \imath R\right]\exp{\imath\left(\omega + f\right)}+ {\cal O}(e)\; ,
 \label{eq13}
\end{equation}
where we have retained terms only to the
zeroth order in the satellite eccentricity $e$. Because of
a small value of the LAGEOS~I eccentricity (about $0.0044$), the
approximation (\ref{eq13}) is sufficient for our purposes. The
quantities $R$ and $T$, represent, as before, the radial and
transverse components of the perturbing force per unit mass.

We adopt the convention of defining the eccentricity excitation
as the components of the vector ($\Delta\dot h$,
$\Delta\dot k$). The first component -- $\Delta\dot h$ -- is to
be called ``the real part of the excitation'', while the second
component -- $\Delta\dot k$ -- ``the imaginary part of the excitation''.
Note that \markcite{{\it Ries et~al.} [1993]} and
\markcite{{\it Tapley et~al.} [1993]} prefer to take $(-\Delta\dot k)$
for the imaginary part of the excitation (the negative of our
definition).

We shall introduce auxiliary quantities
($k\geq 0$ and $j=1,2$)
\begin{mathletters}
 \label{eq14}
\begin{eqnarray}
 A^{(k)}_j &= & s_h \cos kf_j + s_k \sin kf_j \; ,  \label{eq14a} \\
 B^{(k)}_j &= & -s_h \sin kf_j + s_k \cos kf_j \; , \label{eq14b} \\
 \alpha^{(k)}_j &= & n_h \cos kf_j + n_k \sin kf_j \; , \label{eq14c} \\
 \beta^{(k)}_j &= & -n_h \sin kf_j + n_k \cos kf_j \; , \label{eq14d}
\end{eqnarray}
\end{mathletters}
where $(f_1,f_2)$ are the true anomalies of the satellite entry/exit
from the shadow as mentioned above. The variables $(s_k,s_h)$,
resp. $(n_k,n_h)$, in Eqs.~(\ref{eq14}) denote projections
of the satellite spin axis, resp. solar direction, on the
directions towards the orbital perigee and direction in the
orbital plane normal to the perigee direction.
One easily verifies that our quantities $A^{(1)}_{1}$, $A^{(1)}_{2}$,
$B^{(1)}_{1}$ and $B^{(1)}_{2}$ correspond to $A_{1}$, $A_{2}$,
$B_{1}$ and $B_{2}$
introduced by \markcite{{\it Scharroo et~al.} [1991]}. Simple
geometrical considerations show that $\alpha^{(1)}_1 -\alpha^{(1)}_2
= {\cal O}(e)$ and $\beta^{(1)}_1 +\beta^{(1)}_2= {\cal O}(e)$.
We also introduce special unit vectors: ${\hat {\bf a}}$ in the
direction of the ascending node and ${\hat {\bf b}}$ lying in
the satellite orbit plane and perpendicular to ${\hat {\bf a}}$.
Then, one obtains
\begin{eqnarray}
 s_h &=& ({\bf s}.{\hat {\bf a}}) \cos\omega + ({\bf s}.{\hat {\bf
  b}}) \sin\omega \; , \nonumber \\
 s_k &=& -({\bf s}.{\hat {\bf a}}) \sin\omega + ({\bf s}.{\hat {\bf
  b}})\cos\omega \; \nonumber
\end{eqnarray}
and similarly for the solar direction ${\bf n}$.

A straightforward analysis shows that, during the periods when the
satellite orbit crosses the Earth shadow, the effect of {\it the mean
reflectivity coefficient change} denoted by $\Delta C_R$ results in
the averaged eccentricity excitations
\begin{mathletters}
 \label{eq15bis}
\begin{eqnarray}
 {\Delta\dot h} &=& \eta \biggl[ \cos\omega \left(12\pi c_0
 \beta^{(0)} + \alpha^{(2)}_1 - \alpha^{(2)}_2\right) \nonumber \\
 & &  +\sin\omega \left(12\pi c_0 \alpha^{(0)} + \beta^{(2)}_1 -
 \beta^{(2)}_2\right)\biggr] \nonumber \\ & & + {\cal O}(e)
  \; , \label{eq15bisa} \\
 {\Delta\dot k} &=& \eta \biggl[-\cos\omega \left(12\pi c_0
 \alpha^{(0)} + \beta^{(2)}_1 - \beta^{(2)}_2\right) \nonumber \\
 & &  +\sin\omega \left(12\pi c_0 \beta^{(0)} +
 \alpha^{(2)}_1 - \alpha^{(2)}_2\right)\biggr] \nonumber \\ & &
  + {\cal O}(e) \; , \label{eq15bisb} \\
 \eta &\equiv& -{I_0 R_s^2 \over 8 mc}{\Delta C_R\over na}\; .
  \label{eq15bisc}
\end{eqnarray}
\end{mathletters}
Analogous results for the period free of shadow crossings
can be obtained by setting $f_1 = f_2$. The Eqs.~(\ref{eq15bis})
then simplify to
\begin{mathletters}
 \label{eq16bis}
\begin{eqnarray}
 {\Delta\dot h} &= & -{3I_0\pi R^2_s\over 2mc}{\Delta C_R\over
  na} \left({\bf n}.{\hat {\bf b}}\right) + {\cal O}(e)
  \; ,\nonumber \\ & & \label{eq16bisa} \\
 {\Delta\dot k} &= & {3I_0\pi R^2_s\over 2mc}{\Delta C_R\over
  na} \left({\bf n}.{\hat {\bf a}}\right) + {\cal O}(e)
  \; . \label{eq16bisb}
\end{eqnarray}
\end{mathletters}
An analogous analysis for the {\it the asymmetric reflectivity effect}
results in
\begin{mathletters}
 \label{eq15}
\begin{eqnarray}
 {\Delta\dot h} &=& \kappa \biggl[ \zeta\left(\theta_r\right)
 \cos\omega \left(12\pi c_0 B^{(0)} + A^{(2)}_1 -
 A^{(2)}_2\right) \nonumber \\
 & & + \zeta\left(\theta_r\right) \sin\omega \left(12\pi c_0 A^{(0)}
 + B^{(2)}_1 - B^{(2)}_2\right) \nonumber \\
 & & +2\cos\theta_r \cos\omega \left(12\pi c_0
 \beta^{(0)} + \alpha^{(2)}_1 - \alpha^{(2)}_2\right) \nonumber \\
 & &  +2\cos\theta_r \sin\omega \left(12\pi c_0 \alpha^{(0)} +
 \beta^{(2)}_1 -
 \beta^{(2)}_2\right)\biggr] \nonumber \\ & & + {\cal O}(e)
  \; , \label{eq15a} \\
 {\Delta\dot k} &=& \kappa \biggl[- \zeta\left(\theta_r\right)
  \cos\omega \left(12\pi c_0 A^{(0)} + B^{(2)}_1 -
  B^{(2)}_2\right) \nonumber \\
 & &  + \zeta\left(\theta_r\right) \sin\omega \left(12\pi c_0 B^{(0)}
 + A^{(2)}_1 - A^{(2)}_2\right) \nonumber \\
 & &  -2\cos\theta_r  \cos\omega \left(12\pi c_0
 \alpha^{(0)} + \beta^{(2)}_1 - \beta^{(2)}_2\right) \nonumber \\
 & &  +2\cos\theta_r \sin\omega \left(12\pi c_0 \beta^{(0)} +
 \alpha^{(2)}_1 - \alpha^{(2)}_2\right)\biggr] \nonumber \\ & &
  + {\cal O}(e) \; , \label{eq15b} \\
 \kappa &=& {\gamma \over 8\pi na}\; , \label{eq15c}
\end{eqnarray}
\end{mathletters}
for the eccentricity excitations during the shadow crossing phase.
Eqs.~(\ref{eq15}) hold for all models of the asymmetry of the
reflectivity coefficient introduced in the previous section,
provided one chooses a particular form of the $\zeta(\theta_r)$
function and the amplitude $\gamma$. When there is no shadow crossing,
the Eqs.~(\ref{eq15}) simplify to
\begin{mathletters}
 \label{eq16}
\begin{eqnarray}
 {\Delta\dot h} &= & {3 \gamma \over 2na} \left[ \zeta\left(\theta_r\right)
 \left({\bf s} . {\hat {\bf b}}\right)+ 2\cos\theta_r \left(
 {\bf n}.{\hat {\bf b}}\right)\right]\nonumber \\ & &\,
  + {\cal O}(e)  \; , \label{eq16a} \\
 {\Delta\dot k} &= & -{3 \gamma \over 2na} \left[ \zeta\left(\theta_r\right)
 \left({\bf s} . {\hat {\bf a}}\right)+ 2\cos\theta_r \left({\bf n}.
 {\hat {\bf a}}\right) \right] \nonumber \\ & &\, + {\cal O}(e)
  \; . \label{eq16b}
\end{eqnarray}
\end{mathletters}
Next, we proceed with the averaged eccentricity excitations yielded by
{\it the Yarkovsky-Schach effect}. During the shadow crossing phase
one obtains
\begin{mathletters}
 \label{eq17}
\begin{eqnarray}
 {\Delta\dot h} &= & {\tilde \kappa} \biggl\{ \cos\omega
  \Bigl[ A^{(2)}_1 - A^{(2)}_2 + 2f_0 \left(B^{(2)}_1 -
  B^{(2)}_2\right) \nonumber \\ & &  + 12\pi
  {\tilde c}_0 (1+4f_0^2) B^{(0)}\Bigr]  \nonumber \\
  & & +\sin\omega \Bigl[ B^{(2)}_1 - B^{(2)}_2 - 2f_0
  \left(A^{(2)}_1 - A^{(2)}_2\right) \nonumber \\ & & \quad
  + 12\pi{\tilde c}_0 (1+4f_0^2)
  A^{(0)}\Bigr] \biggr\}\nonumber \\ & & + {\cal O}(e)
  \; ,  \label{eq17a} \\
 {\Delta\dot k} &= & {\tilde \kappa}
  \biggl\{ -\cos\omega \Bigl[ B^{(2)}_1 -
  B^{(2)}_2 - 2f_0 \left(A^{(2)}_1 - A^{(2)}_2\right) \nonumber \\
  & &  + 12\pi
  {\tilde c}_0 (1+4f_0^2) A^{(0)}\Bigr] \nonumber \\
  & &  +\sin\omega \Bigl[ A^{(2)}_1 - A^{(2)}_2 + 2f_0
  \left(B^{(2)}_1 - B^{(2)}_2\right) \nonumber \\ & &
  + 12\pi{\tilde c}_0 (1+4f_0^2)
  B^{(0)}\Bigr] \biggr\} \nonumber \\ &  & + {\cal O}(e) \; , \\
 {\tilde \kappa} &=& {{\tilde \gamma} \cos\theta_r \over
 8\pi na (1+4f_0^2)} \; ,  \label{eq17b}
\end{eqnarray}
\end{mathletters}
while
\begin{mathletters}
 \label{eq18}
\begin{eqnarray}
 {\Delta\dot h} &= & {3{\tilde \gamma} \cos\theta_r \over
 2na} ({\bf s} . {\hat {\bf b}}) + {\cal O}(e) \; , \label{eq18a} \\
 {\Delta\dot k} &= & -{3{\tilde \gamma} \cos\theta_r \over
 2na} ({\bf s} . {\hat {\bf a}}) + {\cal O}(e) \; , \label{eq18b}
\end{eqnarray}
\end{mathletters}
out of the shadow crossing period.

As pointed in the previous section, the simplified model (\ref{eq2a})
for {\it the albedo force} allows the analytical treatment of
the long-term eccentricity perturbation. Detailed
computations show that the radial part of the albedo
acceleration (\ref{eq2a}) results in a signal
completely equivalent to the direct solar radiation pressure
perturbation (\ref{eq16bis}) (without shadow crossing)
provided the effective decrease of the radiation pressure coefficient
is calculated as $\Delta C_R \equiv -{1\over 24}\alpha^2 {\bar{\cal A}}
({8\over 3}+3\alpha) C_R$. Analogously, the solar oriented
piece of the albedo
acceleration (\ref{eq2a}) results in a signal equivalent to the
direct solar radiation pressure perturbation (\ref{eq16bis})
(without shadow crossing) with an effective increase of the
radiation pressure coefficient given by $\Delta C_R \equiv {1\over 8}
\alpha^3 {\bar{\cal A}}$. The total albedo eccentricity perturbation thus
becomes
\begin{mathletters}
 \label{eq3al}
\begin{eqnarray}
 {\Delta\dot h} &= & -{3I_0\pi R^2_s\over 2mc}{\Delta C_R^{\rm alb}
  \over na} \left({\bf n}.{\hat {\bf b}}\right) + {\cal O}(e)
  \; , \label{eq3aa} \\
 {\Delta\dot k} &= & {3I_0\pi R^2_s\over 2mc}{\Delta C_R^{\rm alb}
  \over na} \left({\bf n}.{\hat {\bf a}}\right) + {\cal O}(e)
  \; , \label{eq3ab} \\
 \Delta C_R^{\rm alb} &=& -\case19 \alpha^2 {\bar{\cal A}}\; .
  \label{eq3ac}
\end{eqnarray}
\end{mathletters}
If we consider ${\bar{\cal A}}\simeq 0.3$ for the Earth surface
[\markcite{{\it Stephens et~al.,} 1981}] and
$\alpha \simeq 0.5$ for the LAGEOS~I orbit, we preliminarily conclude
that the albedo effect influence on the eccentricity's
long-term excitation is equivalent to an approximately $0.9$~\%
decrease in the radiation pressure reflectivity coefficient.
In the next section, we shall demonstrate that the minus sign
in (\ref{eq3ac}) is responsible for the anticorrelation
of the observed LAGEOS~I excitations and the modeled
albedo signal by \markcite{{\it Martin and Rubincam} [1995]}.

Employing Eqs.~(\ref{eq9}) we also obtain the long-term excitations
of the
eccentricity vector components due to the north/south {\it asymmetry of
the Earth thermal emissivity}
\begin{mathletters}
 \label{eq19}
\begin{eqnarray}
 {\Delta\dot h} &=&\, {{\hat \gamma} \over na} \psi_2(a)
 \sin I + {\cal O}(e) \; ,  \label{eq19a} \\
 {\Delta\dot k} &=&\, {\cal O}(e) \; . \label{eq19b}
\end{eqnarray}
\end{mathletters}
Remember that the ${\hat \gamma}$ amplitude contains the
first order emissivity coefficient ${\cal E}_1$ of Eq.~(\ref{eq8}).
It is highly probable that this coefficient shows seasonal
variations (or more generally long-term variations; \markcite{{\it
Knocke et~al.,} 1988}) which than
yield corresponding long-term variability of the real part
of the eccentricity excitation ${\Delta\dot h}$.

In the next section, we shall employ the formulas (\ref{eq15bis}) --
(\ref{eq19}) for the analysis of the LAGEOS~I eccentricity vector
residuals.

\section{ Application to LAGEOS~I}
Fig.~1 shows both components of the LAGEOS~I eccentricity
excitation: crosses denote the observed residuals given by the
UT/CSR\footnote{ We adopted a data weighting procedure which
closely follows the effect of gradual increase of data precision
as described by \markcite{{\it Tapley et~al.} [1993]}.}
while the solid curve shows the modeled signal accounting for:
(i) perturbation due to the modified value of the mean reflectivity
coefficient $\Delta C_R$ and (ii) the Yarkovsky-Schach effect.
The adjusted value $\Delta C_R = 0.013$ (with a formal $1\sigma$-error
of $0.15\times 10^{-4}$) corresponds to about
$1.2$~\% of the nominal LAGEOS~I reflectivity value $C_R^0 = 1.14$.
Below, we comment independently on the albedo effect
``contribution'' to biasing the satellite reflectivity
coefficient, leaving more detailed discussion of this fact to
the next section.

It has been noted by \markcite{{\it Farinella et~al.,} [1995]} that
the along-track residual analysis does not provide decorrelation of
the two principal parameters of the Yarkovsky-Schach effect: (i) the
amplitude ${\tilde \gamma}$, and (ii) the thermal lag $f_0$. On
contrary, we have checked that the Yarkovsky-Schach excitation
of the eccentricity vector components depends on the value of the
thermal lag parameter $f_0$ only very weakly, being essentially
determined by its amplitude. This opens a remarkable possibility
to decorrelate the two parameters. It should also be noted
that the Yarkovsky-Schach effect acts mainly on the real part of
eccentricity excitation and is responsible for a significant
increase in the observed signal between the epochs $4500$ and
$5600$ (in days after Jan 1, 1977; see Figs.~1 and 2). This
phenomenon is intimately connected with the intricate evolution
of the LAGEOS~I spin axis and the strong dependence of the
Yarkovsky-Schach effect on the spin axis orientation. To obtain
a good fit of this part of the observed signal, we determined the
value of the Yarkovsky-Schach effect amplitude ${\tilde \gamma}$
to be $-175\;{\rm pm/s}^2$ ($1\sigma$-error of $0.14\,{\rm pm/s}^2$),
nearly twice as large as the value approximately adopted by
\markcite{{\it Scharroo et~al.} [1991]}. Combined with this
result, the adjustment of the LAGEOS~I along-track residuals yields
$f_0 = 143^\circ$ ($1\sigma$-error of $8^\circ$) for the thermal
phase lag [or $\tau_{RR} \simeq 5360$ s for the thermal time
lag; see \markcite{{\it Afonso et~al.,} 1989}]. In this study, we
avoid using a close correlation of the phase-lags $\delta_R$, resp.
$f_0$, corresponding to the Yarkovsky, resp. Yarkovsky-Schach,
effects usually considered in the literature and given by [see
\markcite{{\it Rubincam,} 1987}; \markcite{{\it Scharroo et~al.,}
1991}]
\begin{equation}
 \delta_R \equiv \arctan\left(f_0\right) \; .
 \label{eq19bis}
\end{equation}
\markcite{{\it Slabinski} [1995]} has recently shown, by using a
detailed numerical model, that the usual argument stating that the
two phases relate to the same thermal cooling processes is
incorrect. Diverse constituents of the LAGEOS thermal budget
(retroreflectors, mounting rings, aluminium mantle etc.)
interact differently with the visible and infrared
radiation. Complicated thermal transfer processes then
imply different effective phase-lags of the individual
LAGEOS constituents. These LAGEOS structure effects change
even more significantly the temperature variations over
its surface [see Figures 6 and 7 of \markcite{{\it Slabinski,}
1995}]. As a consequence, the simple correlation between the
amplitudes of the Yarkovsky and the Yarkovsky-Schach effects
due to correlation of the corresponding temperature
differences $\Delta T$ does not hold exactly. We have postponed
a few comments on the comparison of our previous results with that
obtained by using the \markcite{{\it Slabinski} [1995]}
model to the next section. Finally, we
mention that we have employed a numerically integrated evolution
of the LAGEOS~I spin axis given by \markcite{{\it Farinella et~al.}
[1995]}; however, nearly the same results have been obtained
when employing the empirical model for the spin axis evolution
given by \markcite{{\it Ries et~al.} [1993]}.

The RMS of the previous combined fit of both the real and imaginary
parts of the eccentricity excitations is about $25.4\; {\rm mas/yr}$.
Individually, the model for the imaginary part of the shows much
better performance, resulting in an RMS on the order of $13.5\; {\rm
mas/yr}$ only. The much poorer agreement with the real part of the
observed residuals can be understood by the following arguments.
First, the modeling of the LAGEOS~I spin axis evolution in the last
period of the surveyed time interval is an intricate task. Small
defects of the model can result in mismodeling of the real part
of the eccentricity excitations though the significant influence on
the Yarkovsky-Schach effect. Second, assuming that the small
perturbation phenomena missing so far in the dynamical model are
related to the satellite spin axis we conclude by a brief inspection
of Eqs.~(\ref{eq20}) that they act dominantly on the real part
of the excitations (as in the case of the Yarkovsky-Schach effect).
Close alignment of the spin axis with the south pole is responsible
for this result.

To see the robustness of the previous result to adding a new
perturbation, we performed a series of tests. The most important of
them was inclusion of the asymmetric reflectivity effect.
The corresponding  model is shown in Fig.~2. We adopted $1.4$~\% for
the degree of the reflectivity asymmetry $\delta$ defined in
(\ref{eq3c}), the value obtained by the LAGEOS~I along-track residual
analysis. We observe that the contribution of the asymmetric
reflectivity effect to the eccentricity excitations is small,
not exceeding $20\; {\rm mas/yr}$, but not totally negligible.
The general properties of the fit
by using only the mean reflectivity coefficient redefinition and the
Yarkovsky-Schach effect (Fig.~1) are not significantly changed,
although the RMS is slightly increased. None of the other
investigated phenomena yielded significant change of the previous
results (see also the discussion in the next section).

In Fig.~3, we plot the real part of the modeled eccentricity
excitation (normalized to unity) due to the asymmetric reflectivity
effect. The absolute scale of the effect is about $20\; {\rm
mas/yr}$ if we adopted the asymmetry parameters corresponding to the
LAGEOS~I along-track residual analysis. Fig.~3 allows the
comparison of different models of the asymmetric reflectivity:
(i) ``the step-like model'' accounting for homogeneous but different
values of the reflectivity on both LAGEOS~I hemispheres
[Eqs.~(\ref{eq2}) and (\ref{eq3})], and (ii) two models with
smooth latitude dependence of the reflectivity excess of the
LAGEOS~I surface proportional to $\cos^n \theta$ with
$n=1$, resp. $n=3$ [here $\theta$ denotes ``the LAGEOS~I colatitude'';
Eqs.~(\ref{eq2}) and (\ref{eq4})]. None of the models significantly
differ from the others, showing an insensitivity of the orbit
perturbation with respect to changes in the reflectivity distribution
on the satellite surface.

We devote a separate discussion to the albedo effect with
particular care to justifying the approximations (\ref{eq2a})
and (\ref{eq3al}). The first test we performed consisted of
comparing the eccentricity excitation (\ref{eq3al}) with results
of entirely numerical models. In the latter models we
computed numerically both the instantaneous albedo acceleration
components $(R,T,W)$ and the values of the eccentricity
vector excitation $({\Delta {\dot h}},{\Delta {\dot k}})$ averaged
over one satellite revolution.

In the case of a constant Earth albedo, we cross-checked the
results by employing both the reduced 1-D formulation
(\ref{eq1a}) of \markcite{{\it Borderies and Longaretti}
[1990]} and the complete 2-D integration over the
visible Earth surface. Results of both numerical models
for the constant albedo Earth show agreement better then
one percent if compared to the results given by (\ref{eq3al}).
Next, we proceeded by testing how the simple formula
(\ref{eq3al}) can include models with variable Earth albedo.
The most important latitude albedo variability
corresponds to the quadrupole zonal term, accounting
for higher values of the Earth albedo in the polar regions. We
note that this term is usually included in the UTOPIA program of
UT/CSR and other high precision orbit determination programs.
Assuming first the reduced latitudinal
albedo model with ${\cal A}_2(\vartheta) ={\cal A}_0+{\cal A}_2
P_2(\cos \vartheta)$ [$\vartheta$ being the colatitude as
before and $P_2(\cos \vartheta)$ the Legendre polynomial],
we tuned the values of the coefficients ${\cal A}_0=0.25$
and ${\cal A}_2=0.2$ so that the mean albedo remains $0.3
(={\bar{\cal A}})$, the value we adopted in our constant albedo
models. Comparison of the constant albedo ${\bar{\cal A}}$ and the
``pole-quadrupole'' ${\cal A}_2(\vartheta)$
eccentricity excitation is given in Fig.~4. We observe, that
the real part of the excitation is nearly insensitive to
the model change, while the imaginary part signal is
reduced by some $20$~\% in the latitude dependent model.
However, redefining the constant albedo value ${\bar{\cal A}}$ to
about $0.283$ we can model the eccentricity excitation due
to the reduced latitudinal model with RMS level of about
$3.4\, {\rm mas/yr}$. Finally, we performed robustness tests
of the previous conclusions to adding other zonal and tesseral
terms in the general spherical harmonics development ${\cal A}
(\vartheta,\varphi)$ of the Earth albedo. The short-dashed
curve in Fig.~4 shows the results of one of these tests. Here
we adopted the coefficients of the albedo development in
the spherical function given by \markcite{{\it Sehnal} [1979]}
preserving, however, the values of $({\cal A}_0,{\cal A}_2)$
as previously. A set of similar tests with randomly changed
values of the zonal and tesseral coefficients of the albedo
development has also been performed. Common features of all these
tests are: (i) the real part of the eccentricity excitation
is very weakly dependent on the changes of the albedo model,
(ii) an appropriate redefinition of the effective mean
value ${\bar{\cal A}}$ can ``absorb'' the complex features
of the albedo model. As an explicit example, we may quote
that after adjusting the effective mean albedo ${\bar{\cal
A}}_{\rm MR} \simeq 0.36$ the constant albedo model (\ref{eq3al})
fits the \markcite{{\it Martin and Rubincam} [1995]} results with
$7.1\, {\rm mas/yr}$ residuals.

We thus conclude that, contrary to the case of the LAGEOS~I
along-track residuals for which the albedo effects contributes
a highly intricate component [\markcite{{\it Vokrouhlick\'y
et~al.,} 1993a,b}; \markcite{{\it Vokrouhlick\'y and Farinella,}
1995a,b}], the albedo effect on the eccentricity vector
can be, to a high degree of precision, accommodated very efficiently
by a slight decrease of the radiation pressure coefficient. As a
consequence, any mismodeling of the averaged, or effectively
averaged in the sense of our previous discussion, Earth albedo value
by the orbit determination program generates a constant bias of the
LAGEOS reflectivity  coefficient. However, because more than $20$~\%
for such a mismodeling does not seem reasonable, we expect
that the albedo effect does not influence the reflectivity
coefficient by more than $0.2$ \%.

Finally, let us comment briefly on the role of the north/south
asymmetry of the Earth's thermal field. Such a phenomenon is likely
to follow seasonal changes in the solar irradiation of the
Earth surface. Earth radiation budget measurement [e.g.
\markcite{{\it Stephens et~al.,} 1981}] and temperature distribution data
[e.g. \markcite{{\it Comiso,} 1994}] suggest this conclusion, allowing
the value up to $50\, {\rm W/m^2}$ for the coefficient ${\cal E}_1$
of the development (\ref{eq8}). Admitting this value for the
quantitative characterization of the north/south asymmetry of
the Earth thermal emissivity, we arrive at the excitation of the
real part of the eccentricity vector of the magnitude
approximately $13.5\, {\rm mas/yr}$, slightly more then $10$~\%
of the observed phenomena. The CSR/UT model accounts for the
essential annual variations of this coefficient; however, mismodeling
of the actual amplitude of the effect for a given epoch
up to $10$ -- $20$~\% is not excluded. Appropriately smaller
mismodeling in the eccentricity excitation is then to be
expected.

\section{ Discussion}
It is a remarkable property of the model proposed
in this paper that it involves a minimum number of perturbation
phenomena and, what is more important, a minimum number of a free
parameters. The only one we need to adjust freely is the value
of the LAGEOS~I surface mean optical reflectivity increase $\Delta
C_R$, while in case of the amplitude of the Yarkovsky-Schach
effect we employ a degree of freedom left by the analysis of the
along-track residuals of LAGEOS~I [see \markcite{{\it Farinella
et~al.,} 1995}]. Because the reflection coefficient was not known
a priori, it must be estimated as part of the orbit determination
process. While the long spans of data available for LAGEOS~I
allow the determination of $C_R$ with a high degree of precision,
the value will be biased by other systematic effects, such as
those discussed here, if there is some degree of correlation.
In particular, the Yarkovsky-Schach effect may
result in a biased estimate of $C_R$.
A weak aliasing in the real part of the signal with the S1
tide is also possible [\markcite{{\it Tapley et~al.,} 1993}].
However, the theoretical models predict only a very
small amplitude for this tidal component.

The value of $-175\, {\rm pm/s}^2$ for the effective amplitude of
the Yarkovsky-Schach effect rests an empirical finding in this
paper. Recent numerical model by \markcite{{\it Slabinski}
[1995]} does show a tendency for increasing the amplitude
of the thermal effects in comparison with previous results.
However, the suggested value for the amplitude of the
Yarkovsky-Schach effect, $H_{\rm c,0}$ in the {\it Slabinski}
terminology, is about $1.7$ times smaller than our empirical
result. We do not have an explanation for this fact, however,
even \markcite{{\it Slabinski} [1995]} adopts the practise
of modifying the theoretical values of the coefficients of his
theory in order to reach better agreement with the
observations. Although {\it Slabinski's} model is an
important step forward, teaching us mainly that the
thermal processes on the LAGEOS satellite are much more
complex in comparison with previous modeling, it is not
without drawbacks. We will not speculate on possible
weaknesses of {\it Slabinski's} theory, but, as an example
including albedo radiation as a source for the satellite
heating may result in a desirable increase of the amplitude
of those effects\footnote{ Notice that the Earth albedo
radiation flux may exceed occasionally the Earth IR flux
at the LAGEOS altitude by a factor of $2$. Contrary to the
IR flux, the albedo flux is not constant and shows rather complicated
variations as well as the shadow periods.}.
These topics demand specialized study.

Figures~1a and 2a show that the main imperfection of our empirical
model is to predict still too small peaks in the last part of the
data (despite of already spurious factor $1.7$ in the Yarkovsky-Schach
effect amplitude). As
those are clearly related to the spin axis orientation (see
Eq.~(\ref{eq18})), we might expect to deliver a part of this observation
to the errors in the spin axis evolution modeling. However, detailed
analysis of various factors entering the modeled signal ($|
{\bf s}.\hat{\bf b}|$, $\cos\theta_r$ etc.) shows that moderate
change of the spin axis orientation, by about $10$ to $20$ degrees,
is not sufficient to explain the high magnitude of of the observed
extrema. As the spin axis evolution scenarios proposed
by \markcite{{\it Ries et~al.} [1993]} and \markcite{{\it Farinella
et~al.} [1995]} gained credibility by allowing to explain the
along-track LAGEOS~I residuals, we prefer to hypothesize the
increase of the Yarkovsky-Schach effect magnitude ${\tilde
\gamma}$ rather then totally abandoning the spin axis evolution
models.

As we do not feel we can improve significantly on the modeling
of the previously discussed phenomena, partly because of lack of the
detailed information (e.g. concerning the possible geometry of the
optical anisotropy of the satellite) and partly because of
difficulties in the physical (often empirically adjusted) modeling,
and in view of the remaining (post-fit) residuals in the
eccentricity excitations we turn our investigation towards still
another non-gravitational effects.

We recall the general rule mentioned in the Introduction:
the tricky, but interesting,
circumstance in considering the eccentricity excitations is that
the work is already largely constrained by a good fit to the
LAGEOS~I semimajor axis. Put in other words, one is not allowed
to introduce effects which could spoil the semimajor axis analysis.
It appears, that the best candidates are forces
not depending of the orbital motion. Indeed it is easy to verify that
for such a perturbation ${\bf \gamma}$, the long term effects on the
semimajor axis and eccentricity vector are given by
\begin{mathletters}
 \label{eq20}
\begin{eqnarray}
 {\Delta\dot a} &=& {\cal O}(e) \; , \label{eq20a} \\
 {\Delta\dot h} &=& \phantom{-}{3 \over 2na} {\bf \gamma.{\hat b}}+
  {\cal O}(e) \; ,  \label{eq20b} \\
 {\Delta\dot k} &=& -{3 \over 2na} {\bf \gamma.{\hat a}}+ {\cal O}(e)
 \; . \label{eq20c}
\end{eqnarray}
\end{mathletters}
As conjectured above, the forces of interest should favorably depend on
the spin axis orientation. Armed with these principles, we consider
below a set of non-gravitational effects which are possible candidates
for helping to cure the previous difficulties of explaining the LAGEOS~I
eccentricity residuals.

One easily verifies that the thermal effect dominantly
contributing to the along-track residual signal -- the
Rubincam version of the Yarkovsky thermal effect
characterized by the Earth IR radiation thermal heating --
yields only minor action on the eccentricity residuals.
An upper limit is estimated to be about $0.05\,
{\rm mas/yr}$ for both the real and imaginary components of the
excitations. Similar result holds for the charged or
neutral particle drag observed via the along-track residuals.

The next group of potentially interesting non-gravi\-ta\-ti\-onal effects
perturbing the LAGEOS~I motion are the electromagnetic phenomena
embodying: (i) charging of the satellite due to its interaction
with the ionospheric vicinity [\markcite{{\it Al'pert,} 1983};
\markcite{{\it Afonso et~al.,} 1985}], (ii)
induction of the magnetic moment of the satellite due to its
rapid rotation in the Earth's magnetic field [\markcite{{\it Bertotti
and Iess,} 1991}; \markcite{{\it Farinella et~al.,} 1995}],
(iii) electromagnetic
force exerted on the charged and magnetized satellite moving
in the Earth's magnetic field [partly studied by \markcite{{\it Sehnal,}
1969}; \markcite{{\it Vokrouhlick\'y,} 1989}].

We shall limit ourselves to the
inclined dipole approximation for the Earth magnetic field
(denoting $D \simeq 11^0$ its axis deviation from the Earth
rotation axis). The satellite equilibrium charge may be
considered in the first approximation constant, but generally
short- and long-term variations are likely to be expected.
With the constant charge approach one easily shows
that the Lorentz force induced by the Earth's magnetic field
does not yield any long-term excitations of the eccentricity
vector components. The same analysis also shows that only
improbable short-term variations of the satellite charge may
result in the long-term eccentricity signals. However, it should
be pointed that the details of the charging mechanism of a rapidly
rotating conducting satellite has not been studied yet in context
of the satellite dynamics (\markcite{{\it Al'pert} [1983]} or
\markcite{{\it Afonso et~al.} [1985]} restrict themselves to the
case of a non-rotating spherical satellite).

Similarly, interaction of the Earth's magnetic field with the
induced magnetic moment of the satellite results in a perturbing
acceleration. Contrary to the case of the LAGEOS~I despinning
scenario, the more important component here is the potential part of
the satellite magnetic moment [first term in Eq.~(4) of
\markcite{{\it Bertotti and Iess,} 1991}].
Denoting $\alpha'(\omega_r)$ the real part of the volumic average
polarizability of the LAGEOS~I body, $V$ its volume one obtains
for the magnetic force
\begin{equation}
 {\bf a} = 2 {V \alpha'(\omega_r) \over m} \left({\bf s}.{\bf B}\right)
  \left({\bf s}.\nabla {\bf B}\right)\; , \label{eq28}
\end{equation}
where ${\bf s}$ is the vector orientation of the satellite
spin axis and ${\bf B}$ is the induction of the Earth magnetic field.
Clearly, this force depends not only on the satellite axis
orientation but also on the law by which the satellite angular
velocity $\omega_r$ decays. Careful analysis again shows that
the (inclined) dipole approximation for the Earth's magnetic field
${\bf B}$ in (\ref{eq28}) yields eccentricity excitations only in the
first order in the satellite eccentricity itself making the values
negligibly small. Double occurrence of the magnetic field induction
in the formula (\ref{eq28}) still leaves the possibility of interpreting
one of them as the dipole part and the second as a quadrupole
part of the Earth's field. Indeed, such cross-contributions
yield eccentricity excitation which does not vanish even for
the circular (zero eccentricity) orbit. Its value has been, however,
estimated to be smaller than $0.01\,{\rm mas/cy}$.

We thus conclude that it is not likely that the
electromagnetic phenomena contribute significantly to the
observed LAGEOS~I eccentricity excitations.

Most of the software packages employed for the first level orbit
treatment neglect the intricate structure of the transition
into/from the Earth's shadow in context of the perturbation
raised by the direct solar radiation pressure (the dominant
non-gravitational effect for high altitude satellites). The
error linked with those phenomena may be twofold: (i) ``numerical'',
related to the bad behaviour of particular integrators when the
forces show discontinuities, and (ii) ``theoretical'', related to
the proper dynamical effect of the penumbra phase radiation pressure.
As for the first item, we note that
since this phenomenon is observed with different numerical
integrators and a variety of integrator stepsizes, it is
unlikely to be associated with the numerical integration process.
%It should be noted, however, that \markcite{{\it Lundberg} [1995]}
%has recently shown that particularly ill-chosen values of the
%integration stepsize can create an aliasing of the shadowing
%function that can resemble ??? force effects.
 On side of the
second issue mentioned above, we remind that the
natural mechanism of the asymmetric conditions of the the
shadow entry/exit phases led \markcite{{\it Vokrouhlick\'y et~al.}
[1993c, 1994]} to investigate in more detail implications of
those intricate phenomena. The surprising result, that despite
of the extreme localization in time ``the penumbra effect'' can
provoke a long-term along-track signal (though of a minor amplitude)
made us to investigate the same phenomenon in context of
the eccentricity excitations. Employing the model similar to
\markcite{{\it Vokrouhlick\'y et~al.} [1994]} we conclude that
the penumbra excitations occasionally exceed a few
milliarseconds/year, but generally stay smaller than $1$~mas/year,
small enough with respect to the residuals to leave
the penumbra effect on the edge of primary interest for the
explanation of the puzzling LAGEOS~I eccentricity residuals.

\section{ Conclusion}
The main results of this paper can be summarized as follows:
\begin{enumerate}
\item we propose a model for explaining the observed eccentricity
 excitations of the LAGEOS~I orbit by accounting for two effects:
 (i) a $1.2$~\% increase of the mean value of the LAGEOS~I surface
 reflectivity coefficient (with respect to the nominal value
 of $1.14$), and (ii) the Yarkovsky-Schach thermal effect;
\item the combined analysis of the LAGEOS~I along-track residuals and
 the eccentricity excitations allow one to decorrelate the parameters
 of the Yarkovsky-Schach thermal effect leading to the following
 empirical values: (i) ${\tilde \gamma} \simeq -175\; {\rm pm/s}^2$
 for its amplitude, and (ii) $f_0 \simeq 143^\circ$ for the value
 of the thermal lag;
\item the albedo effect influence of the long-term excitations
 of LAGEOS~I eccentricity is equivalent, with unexpected precision,
 to an effective decrease of the radiation pressure coefficient
 by about one percent; however, as the UTOPIA program accounts
 for the main features of the Earth albedo distribution, we
 do not expect that the albedo effect biases the value of LAGEOS
 reflectivity coefficient by more than $0.2$ percent;
\item none of the other investigated phenomena (notably those related
 to the satellite spin axis orientation) contribute significantly
 to the eccentricity excitations.
\end{enumerate}

Though discussed in some detail above, the necessity of redefining the
adjusted value of the LAGEOS~I surface reflectivity coefficient still
needs a careful study. Some routes, especially the remarkable
reduction of the albedo long-term eccentricity effect, for biasing
the adjusted value of the $C_R$ coefficient has been mentioned in this
paper. Another outcome of the current study which
desires theoretical understanding in future is the adjusted amplitude
of the Yarkovsky-Schach thermal effect which is in conflict with
all previous studies. The recent work by \markcite{{\it Slabinski}
[1995]} shows that previous estimates of the LAGEOS
thermal parameters (the amplitude of the thermal effects and
the corresponding phase-lags) were oversimplified and that
the intricate heat transfer between the different parts
of the satellite body tends to increase the amplitude of the
effects. However, the theory suggested value is still
by a factor of $1.7$ smaller than our empirical value
$-175\; {\rm pm/s}^2$. Following {\it Slabinski's} work, it
is thus necessary to revisit the theoretical modeling of the
LAGEOS body thermal effects in the future.

Other contributors to the LAGEOS~I eccentricity excitations observed
by UT/CSR model among a tiny effects of both the gravitational or
the non-gravitational origin are not excluded. However, we
were unable to find really suitable candidates. As for
the former, we can note for instance the S1 tide component
has a correct characteristics.

Finally, let us mention that apart from a better understanding of
how the diverse phenomena usually appearing in the dynamical
models for the artificial satellite motion (such as the
radiation pressure, thermal effects, tides etc.) act in the Earth
environment, the solution of the eccentricity residual problem
potentially hides answers on the important tests of physics.
Determination of the time variability of the $J_3$ zonal
coefficient is among the most urgent tasks [see e.g.
\markcite{{\it Tapley et~al.,} 1993}; \markcite{{\it Eanes and
Bettadpur,} 1995}]. Another interesting use of the ``unpolluted''
LAGEOS~I eccentricity residuals has been recently suggested by
\markcite{{\it Damour and Esposito-Far\`ese} [1994]} in respect
to the testing for the preferred frame effects of gravitational
theories [see also \markcite{{\it Peterson and Eanes,} 1995}].

\appendix
\section{Appendix}
In the following, we shall briefly outline the main ideas needed for
calculation of the radiative force exerted on a spherical satellite
employed in the main text. We shall rely on the formalism developed
by \markcite{{\it Vokrouhlick\'y et~al.} [1993a, 1994]}.
Our method consists of two different steps: (i) computation of
the radiative flux vector (or more generally components of the
radiative field stress-tensor), (ii) relating previous quantities
to the searched radiative force exerted on the body. For the case
of a spherical satellite the second step (ii) is already trivial
[see Eq.~(A4)].

Important technical characteristics of the mentioned method
consists of usage of the local satellite frames
for performing all integrations. Often, the results are more
simple and conceptually more transparent than integration ``over the
Earth surface elements'' (see Fig.~5 for geometrical quantities
used below). Considering uniquely the north/south asymmetric part of
the Earth thermal emissivity (\ref{eq8}) one obtains for the radial
component -- $F_{z'}$ -- and meridional component -- $F_{x'}$
-- of the radiative flux (the longitudinal component --
$F_{y'}$ -- vanishes from obvious symmetry reasons)
\begin{mathletters}
\begin{eqnarray}
 F_{x'} &= & {{\cal E}_1 \over \pi} \int d\nu d\phi\; \mu(\nu,
 \phi)\, \sqrt{1-\nu^2} \cos\phi \; , \nonumber \\ & & \\
 F_{z'} &= & {{\cal E}_1 \over \pi} \int d\nu d\phi\; \mu(\nu,\phi)\,
 \nu \; ,
\end{eqnarray}
\end{mathletters}
where
\begin{eqnarray}
 \alpha \mu(\nu,\phi) &=& \cos\vartheta_s - \left[\nu - \sqrt{\nu^2-
 \beta^2}\right] \\
 & & \; \times \left[\nu \cos\vartheta_s - \sqrt{1-\nu^2} \cos\phi
 \sin\vartheta_s\right] , \nonumber
\end{eqnarray}
The integration (A1) is performed over the domain
$\nu \times \phi = (\beta,1)\times(0,2\pi)$ ($\beta$ as previously).
Both integrals (A1) can be worked out analytically
\begin{mathletters}
\begin{eqnarray}
 F_{x'} &= & {{\cal E}_1 \over 2} {a \over R_e} \psi_1(a) \cos
  \vartheta_s \; ,\\
 F_{z'} &= & {{\cal E}_1 \over 4} {a \over R_e} \left[\psi_1(a)-
  \psi_2(a) \right] \sin\vartheta_s \; ,
\end{eqnarray}
\end{mathletters}
with the functions $\psi_1(a)$ and $\psi_2(a)$ given by
Eqs.~(\ref{eq11}) of the main text.

Once all components of the radiative flux ${\bf F}$ in the satellite
local frame are computed, one can relate them to the components of
the perturbing radiative acceleration for a simple shaped satellites.
Notably, for a spherical satellite one proves [e.g. \markcite{{\it
Vokrouhlick\'y et~al.,} 1993a,c}]
\begin{equation}
 {\bf a} = {\pi R_s^2 \over mc} C_R {\bf F} \; ,
\end{equation}
where $R_s$ is the satellite's radius, $m$ its mass, $c$
the light velocity and $C_R$ the radiation pressure coefficient
which depends on the optical properties of the satellite surface.

\acknowledgements
 We thank L. Anselmo, F. Barlier, P. Exer\-tier and P. Fari\-ne\-lla
 for helpful discussions. D.V. performed the study when staying at
 the OCA/CERGA, Grasse (France) thanks to the H. Poincar\'e
 research fellowship.

\newpage

\begin{references}
\reference
Al'pert, Y.L., {\it The Near-Earth and Interplanetary Plasma},
 Vol.~2, Cambridge Univ. Press, Cambridge, 1983.

\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, \jgr {\it , 90,} 9381-9398, 1985.

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

\reference
Barlier, F., M. Carpino, P. Farinella, et~al., Non-gravitational
 perturbations on the semimajor axis of LAGEOS, {\it Annales
 Geophysicae,} {\it 4A,} 193-210, 1986.

\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
Borderies, N., and P.Y. Longaretti, A new treatment of the albedo
 radiation pressure in the case of a uniform albedo and of a
 spherical satellite, {\it Celest. Mech.,} {\it 49,} 69-98, 1990.

\reference
Ciufolini, I., P. Farinella, A.M. Nobili, D. Lucchesi, and L.
 Anselmo, Results of a joint ASI-NASA study on the LAGEOS gravitomagnetic
 experiment and the nodal perturbations due to radiation pressure
 and particle drag effects, {\it Nuovo Cimento,} {\it B108,} 151-162,
 1993.

\reference
Comiso, J.C., Surface temperatures in the polar regions from
 Nimbus~7 temperature humidity infrared radiometer, \jgr {\it , 99,}
 5181-5200, 1994.

\reference
Damour, T., and G. Esposito-Far\`ese, Testing for the preferred-frame
 effects in gravity with artificial Earth satellites, {\it Phys.
 Rev. D,} {\it 49,} 1693-1706, 1994.

\reference
Eanes, R.J., and S. Bettadpur, Temporal variability of Earth's
 gravitational field from satellite laser ranging, to appear in
 Proceedings of the International Association of Geophysics, Symp.
 G3, {\it Global Gravity Field and its Temporal Variations},
 R.H. Rapp, ed., Springer-Verlag, 1995.

\reference
Farinella, P., D. Vokrouhlick\'y, and F. Barlier, A consistent
 theoretical model for the rotation of LAGEOS, \jgr , in press, 1995.

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

\reference
Ferraz-Mello, S., Analytical study of the Earth's shadowing effects
 on satellite orbits, {\it Celest. Mech.,} {\it 5,} 80-101, 1972.

\reference
Knocke, P., J.C. Ries, and B.D. Tapley, Earth radiation pressure effects
 on satellites, in {\it Proceedings of the AIAA/AAS Astrodynamics
 Conference}, Am. Inst. Aeron. Astronaut., Washington, D.C.,
 577-587, 1988.

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

\reference
L\'ala, P., and L. Sehnal, The Earth's shadowing effects in the
 short-periodic perturbations of satellite orbits, {\it Bull. Astron.
 Inst. Czechosl.,} {\it 20,} 327, 1969.

%\reference
%Lundberg, J.B., Mitigation of satellite orbit errors resulting from
% numerical integration across shadow boundaries, to appear in
% {\it Proceedings of the AIAA/AAS Astrodynamics Specialist Conference},
% American Astronaut. Society, San Diego, CA, 1995.

\reference
Martin, Ch.F., and D.P. Rubincam, Effects of Earth albedo on the
 LAGEOS~I satellite, \jgr , submitted, 1995.

\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
Peterson, G., and R.J. Eanes, Testing the preferred-frame effect
 using the LAGEOS satellites, {\it Technical Memorandum
 CSR-TM-95-01}, University of Texas at Austin, 1995.

\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
Scharroo, R., K.F. Wakker, B.A.C. Ambrosius, and R. Noomen, On the
 along-track acceleration of the LAGEOS satellite, \jgr
 {\it , 96,} 729-740, 1991.

\reference
Sehnal, L., The theory of a charged satellite in the Earth's magnetic
 field, SAO Special Report No.~271, Cambridge, 1969.

\reference
Sehnal, L., The Earth albedo model in spherical harmonics, {\it Bull.
 Astron. Inst. Czech.,} {\it 30,} 199-204, 1979.

\reference
Sehnal, L., Effects of the terrestrial IR radiation pressure on the
 motion of an artificial satellite, {\it Celest. Mech.,} {\it 25,}
 169-179, 1981.

\reference
Slabinski, V.J., A numerical solution for LAGEOS thermal thrust:
 the rapid-spin case, {\it Celest. Mech.,} submitted, 1995.

\reference
Stephens, G.L., G.G. Campbell, and T.H. Vondner Haar, Earth radiation
 budgets, \jgr {\it , 86,} 9739-9760, 1981.

\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., The geomagnetic effects on the motion of an
 electrically charged artificial satellite, {\it Celest. Mech.,}
 {\it 46,} 85-104, 1989.

\reference
Vokrouhlick\'y, D., and P. Farinella, Specular reflection of sunlight from
 wavy ocean surfaces and the albedo effect on satellite orbits, I.
 A statistical model, {\it Astron. Astroph.,} {\it 298,} 307-322, 1995a.

\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. Astroph.,}
{\it 301,} 282-289, 1995b.

\reference
Vokrouhlick\'y, D., P. Farinella, and D. Lucchesi, Al\-be\-do perturbation
 models: general formalism and application to LAGEOS, {\it Celest. Mech.,}
 {\it 57,} 225-244, 1993a.

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

\reference
Vokrouhlick\'y, D., P. Farinella, and F. Mignard, Solar radiation pressure
 perturbations for Earth satellites, I. A complete theory including
 penumbra transitions, {\it Astron. Astroph.,} {\it 280,} 295-312,
 1993c.

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

\reference
Watkins, M.M., B.D. Tapley, R.J. Eanes, and J.C. Ries, Detailed
 simulation of the LAGEOS 3  measurement of the Lense-Thirring precession,
 in {\it Proceedings of the 1$^{\rm st}$ W. Fairbank Conference on
 Relativistic Gravitational Experiments in Space}, Rome, Italy,
 September 1990.
\end{references}

%\end{document}
%\newpage
%.
\newpage

\begin{figure}
\figurewidth{35pc}
\caption{ Real -- {\bf a)} -- and imaginary -- {\bf b)} -- parts of
 the LAGEOS~I eccentricity excitation vs. time. The crosses
 show the observed residuals, while the solid line corresponds
 to the modeled signal accounting for: (i) direct solar radiation
 pressure related to the $1.2$~\% redefinition of the mean
 reflectivity coefficient of LAGEOS, and (ii) the Yarkovsky-Schach
 effect. Parameters of the latter effect are: (i) the amplitude
 ${\tilde \gamma} = -175\; {\rm pm/s}^2$, and (ii) thermal lag
 corresponding to $\tau_{RR} = 3755$ s (or $f_0 \equiv n\tau_{RR}
 = 143^\circ$). These values suit well both the along-track and
 the eccentricity residuals adjustments (see comments in the text).}
\end{figure}

\begin{figure}
\figurewidth{35pc}
\caption{ The same as Fig.~1 but now the asymmetric reflection
 effect, necessary for successful adjustment of the LAGEOS~I
 along-track residuals, included. Though slightly worse, the modeling
 of the eccentricity residuals qualitatively and quantitatively
 suit explanation of the observed values.}
\end{figure}

\begin{figure}
\figurewidth{20pc}
\caption{ The imaginary part of the eccentricity excitation (normalized
 to unity) due to the asymmetric reflectivity effect vs. time.
  Different models of the ``LAGEOS~I latitude''
 distributions of the reflectivity coefficients compared: (i) the
 step model -- solid line, (ii) the model with the $\cos\theta$
 dependence of the reflectivity parameter ($n=1$) -- dashed line,
 and (iii) the model with the $\cos^3\theta$ dependence of the
 reflectivity parameter ($n=3$) -- dash-dotted line. No
 significant dependence of the results on the adopted model is
 observed.}
\end{figure}

\begin{figure}
\figurewidth{35pc}
\caption{ Real -- {\bf a)} -- and imaginary  -- {\bf b)} -- parts
 of the eccentricity excitation due to the albedo models vs.
 time. The covered  period, March 1985 to June 1989, is the same as
 in \markcite{{\it Martin and Rubincam} [1995]}. The solid curve
 accounts corresponds to the model with reduced latitudinally
 dependent albedo distribution ${\cal A}_2(\vartheta)$, the
 long-dashed curve to the constant albedo model ${\bar{\cal A}}$
 and the short-dashed curve to a particular model with
 general albedo distribution ${\cal A}(\vartheta,\varphi)$
 introduced in the text.}
\end{figure}

\begin{figure}
\figurewidth{20pc}
\caption{ Reference frames and angular variables
 introduced in the text. Notice orientation of the local reference
 frame $(x,y,z)$ of the satellite, the Earth surface visible cap
 and an example of the thermal photon which contributes to the
 radiative effect discussed in the text.}
\end{figure}

\end{document}
