The LArase Satellites Spin mOdel Solutions (LASSOS): a comprehensive model for the spin evolution of the LAGEOS and LARES satellites
Massimo Visco David M. Lucchesi
gr-qc




	The LArase Satellites Spin mOdel Solutions (LASSOS):  a comprehensive model for the spin evolution of the LAGEOS and LARES satellites 

	M. Visco 1,2,
	D. M. Lucchesi1,2,3

	1IAPS - INAF , Via del Fosso del Cavaliere, 100, 00133 Roma, Italy 

		2 INFN, Sez. Tor Vergata, Via della Ricerca Scientifica 1, 00133  Roma, Italy 
 
		3 ISTI - CNR,  via G. Moruzzi 1, 56124 Pisa, Italy 
 
		[3em]

         


The two LAGEOS and LARES are laser-ranged satellites tracked with the best accuracy ever achieved. Using their range measurements many geophysical parameters were calculated and some General Relativity effects were directly observed. To obtain precise and refined measurements of the effects due to the predictions of General Relativity on the orbit of these satellites,  it is mandatory to model with high precision and accuracy all other forces, reducing the free parameters introduced in the orbit determination. A main category of non-gravitational forces to be considered are those of thermal origin, whose fine modeling strongly depends on the knowledge of the evolution of the spin vector.  
We present a complete model, named LASSOS, to describe the evolution of the spin of the LAGEOS and LARES satellites. In particular, we solved Euler equations of motion considering not-averaged torques. This is the most general case, and the  predictions of the model well fit the available observations of the satellites spin. 
We also present the predictions of our model in the fast-spin limit, based on the application of averaged equations. The results are in good agreement with those already published, but with our approach we have been able to highlight small errors within these previous works. 
LASSOS was developed within the LARASE research program. LARASE aims to improve the dynamical model of the two LAGEOS and LARES satellites to provide very precise and accurate measurements of relativistic effects on their orbit, and also to bring benefits to geophysics and space geodesy. 





§ INTRODUCTION

The two LAGEOS (LAser GEOdynamics Satellite) and LARES (LAser RElativity Satellite) are passive Earth orbiting satellites. They are almost spherical in shape, with a diameter of about 60 cm and mass of about 407 kg and 405 kg, respectively for LAGEOS and LAGEOS II, and a diameter of 36.4 cm and a mass of about 387 kg in the case of LARES.   A large number of retroreflectors,  the so-called cube-corner retro-reflectors (CCRs), is distributed in the form of rings over the satellites' surface. These CCRs allow a precise tracking of the orbit of the satellites by means of a network of Earth bound stations coordinated by the International Laser Ranging Service (ILRS), see <cit.>. In particular, the two LAGEOS have embedded in their surface 426 CCRs, 422 made of fused silica and 4 of germanium, while the smaller LARES has 92 CCRs, all made of fused silica.

The older LAGEOS was launched by NASA in 1976, on May 14. LAGEOS has an orbit inclination of about 109.9^∘ over the Earth's equator, a semi-major axis of about 12,270 km and an eccentricity of about 0.004.
LAGEOS II was launched by NASA/ASI in 1992, on October 22. LAGEOS II has an orbit inclination of about 52.7^∘ over the Earth's equator, a semi-major axis of about 12,162 km and an eccentricity of about 0.014. The orbital periods of the two satellites are respectively of 225.4 min and  222.4 min. LARES was launched by ASI in 2012, on February 13. LARES has an inclination of 69.5^∘ over the Earth's equator, a semi-major axis of about 7,820 km, an eccentricity close to zero and an orbital period of 114.9 min.

The ground stations of the ILRS assure, by means of the Satellite Laser Ranging (SLR) technique, almost continuous measurements of the round trip time of narrow laser pulses and, consequently, of the satellites' distance, which is usually called range. High quality range data, called Normal Points (NPs), are produced averaging the tracking measurements over appropriate time spans, see <cit.>.

The orbital parameters, obtained from the NPs with least squares methods, can be compared with those calculated using parametric dynamical models, allowing on one side an accurate evaluation of many geophysical parameters (see for instance <cit.>,
<cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>) and, on the other side, a direct measurement of the relativistic gravitational effects acting on the satellite orbit, see for instance <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. 

The dynamical model used to calculate the orbit of these geodetic satellites must include also forces related with gravitational and non-gravitational perturbations (NGP). A main disturbing role is played by the surface forces strictly connected with the evolution of the satellites' spin vector (orientation and rate), we refer to <cit.>, <cit.>, <cit.>, <cit.>, <cit.> for details. These are the Earth-Yarkovsky and Yarkovsky-Schach thermal effects (see <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>). 

The two LAGEOS satellites and LARES were injected in their orbit with an initial rotation that slows down in time under the magnetic torque, and with an initial orientation of their spin that were subject to an evolution from the precession due to the combined action of three main forces: magnetic, gravitational and surface forces. Several spin models were introduced in the past. Among these models, those that better describe the experimental data are the ones valid in the so called fast-spin approximation.
In this approximation, the solutions of the models are valid as long as the satellite rotation period is much shorter than the orbital one. 
Nowadays, while LARES spins with a period of about 1000 s, one order of magnitude smaller then the orbital period, the two LAGEOS' spin periods are, respectively, about 3· 10^5 s for the older satellite and about 3· 10^4 s for the younger one. Both satellites spin less than one revolution during each orbit, therefore a more general model is needed beyond the fast-spin approximation.

The aim of the present work is to introduce our spin model. This model has been developed within the activities of the LAser RAnged  Satellites Experiment (LARASE) <cit.>. The main new characteristic of our model is its capacity to simulate completely the spin trend of LAGEOS and LARES satellites in the general case, with no restriction to the mentioned fast-spin approximation. In other words, the model is valid for any rotational period of the satellite. After a brief summary about the models developed in the past (see section <ref>), in section <ref> we focus on the equations we used and the different torques we considered.  Starting from these formulas we have been also able to calculate the  solution valid for averaged equations in the limit of the fast-spin approximation. In section <ref> we describe the numerical integration methods used to solve the equations in the general case and we compare the obtained results with the available experimental observations.
In section <ref> we provide a discussion on the importance of a refined modeling for the spin evolution, and we also focus on the order of magnitude of other torques that are however negligible with respect to those taken into account in previous sections.
Finally, in section <ref> our conclusions and recommendations are provided together with the work needed to further improve the free parameters that characterize our model for the spin of LAGEOS-like satellites.





§ PREVIOUS MODELS FOR THE SPIN OF THE TWO LAGEOS SATELLITES
 

The problem of the comprehension of the rotational behavior of the older of the two LAGEOS satellites was first studied by <cit.>.  In that paper the authors, in order to model the spin of LAGEOS, considered two torques acting over the satellite. The first torque is the one produced by the interaction of the magnetic moment  of the satellite with the Earth's Geomagnetic field. The magnetic moment is the one produced by the eddy currents (or Foucault currents) induced in the conductive body of the satellite from its rotation and motion in the same field. The second torque arises from the action of the Earth gravitational field, because of the non spherical mass distribution of the satellite. 

<cit.> solved the problem in the fast-spin approximation by averaging the torques on the satellite orbital period and on the Earth rotational period.  
In that work for the first time the measured decay of the spin of LAGEOS  was explained with a good accuracy and, at the same time, the complexity and variety of the results in different regions of the space of parameters were shown.
An additional significant result of <cit.> model was that, with respect to the inversion of the initial sense of rotation of the satellite, the evolution of both the direction and rate of the spin are not invariant.

The analysis of the spin model of LAGEOS — as well as the results of this first successfully work — were successively extended by <cit.> and <cit.>.  In <cit.>, some small errors present in the formulas published by <cit.> were corrected and, most importantly, these authors extended the model in order to take in account of a possible misalignment between the spin and symmetry axes. They also showed that the <cit.> spin model leads to a successful fit of the along–track residuals of the  two LAGEOS related with the orbital perturbation produced by the main thermal effects, once the force mechanisms related with the direction of the satellite spin vector were properly modeled. 

Finally, <cit.> considered other possible torques connected to atmospheric skin forces and radiation pressure: one first torque connected to the difference of reflectivity between the two LAGEOS's hemisphere <cit.> and the second arising from the non coincidence between the geometric center and the center of mass of the satellite. However, they wrongly concluded that likely these torques would not play an important role in the evolution of the spin of the two LAGEOS satellites, 
because negligible and that, at first order, they should be ignored. We refer to <cit.> for details. 
This conclusion was probably dictated by the fact that these authors have not tuned their spin model with the available (at that time) observations of the spin of the two LAGEOS satellites, but, conversely, with the residuals in their semi-major axis.

These two (much smaller) torques were inserted later in the LOSSAM  (LageOS Spin Axis Model) model (see <cit.>), that till now was the most complete and reliable model to describe the spin behavior of LAGEOS satellites in the fast-spin approximation (<cit.>) and recently it was also applied with success to LARES, see <cit.>.
Indeed, the LOSSAM model has represented a huge improvement in  modelling the evolution of the spin of the two LAGEOS satellites with respect to previous models, even if the bulk of the characteristics of the torques considered are exactly the one previously developed by <cit.> — for the magnetic and gravitational torques — and by <cit.> for what concerns the additional (minor) torques related with the radiation pressure asymmetry and the possible offset between the satellites geometric center and their center of mass.


The key aspects that have led to the development of a so successful model can be summarized in the following main points:

	
  * 
		the collection of all available observations (up to 2007) of the spin of the two LAGEOS satellites, both in their rate and orientation, see <cit.>, <cit.>, <cit.> and <cit.>;
	
  * a careful analysis of all the parameters that enter in the models, i.e., in the mathematical expressions of the various torques;
	
  * the tuning of a suitable subset of these parameters by means of a least-squares fit of the LOSSAM model to all the available observations.



It is important to underline that, in order to analyze the spin evolution under the most general conditions, the torques expressions used in the equations should not be averaged. In this case it is handy to work in the rotating frame of the satellite and use the Euler's equations. <cit.> were the first to formulate the problem in this way. Their aim was to qualitatively analyze the time evolution of the satellites rotation, and they did not care if their solution did not fit to the experimental data. Anyway, their model was too much simplified to replicate the measured data because of several assumptions. In fact: i) the expression adopted for the polarizability of the satellite is too far from reality, ii) the expression of the magnetic torque is the same of <cit.> and <cit.>, that is valid only in the fast-spin approximation, iii) the Earth magnetic field is simplified as a pure dipole field along the Earth rotation axis, finally iv) the equations of motion are written in a non inertial reference frame, neglecting the precession of the orbital plane with respect to an inertial reference frame related to the Earth.

However, despite of these excessive oversimplifications, <cit.> have correctly introduced the general approach to be followed and based on the resolution of the Euler dynamical equations, and they also introduced the correct general expression for the rate of precession in the case of the gravitational torque, correcting the expression given by <cit.> (that was anyway used later by <cit.>).

The  <cit.> approach based on the numerical solution of the full set of the Euler equations was resumed a few years later by <cit.>, removing the simplifications of the previous model but, at the end, still with an unsuccessful fit to the available observational data. 
Later on, this more general model was also considered by <cit.> and <cit.>, that developed the so called ”real-time” LOSSAM formulation.  We could not find in literature any example of practical application of these models. In fact, <cit.> wrote in his PhD thesis work that 
 “real-time” LOSSAM had not given a better result than the model developed by <cit.>, giving implicitly a negative judgment on the results obtained by both models.   



§ THE LASSOS MODEL FOR THE EVOLUTION OF THE SPIN


In the following we present our new model for the spin of LAGEOS and LARES satellites, that aims to be general and complete. All the main known torques were considered and their mathematical expressions were written in a general (not averaged) way in the rotating frame of the satellite adopting the Euler's rotational equations. 





 §.§ Reference frames

For our analysis, in the following, we will adopt two different reference frames: 

	
  * the Earth Mean Equator and Equinox of Date Frame (J2000)
	
  * the satellite Body Frame (BF).


The J2000 is identified by three axes (x̂^E, ŷ^E, ẑ^E). The origin of this Cartesian reference frames coincides with the Earth's center, the (x̂^E,ŷ^E) plane coincides with the Earth's mean equatorial plane. The x̂^E  axis is directed to the vernal equinox , intersection between the mean equator and the mean equinox of date. 
The J2000 represents a quasi-inertial frame.

The BF (also known as Body-Fixed reference frame) is centered in the body's center of mass with the axes (x̂^b, ŷ^b ,ẑ^b) aligned along the principal axes of inertia. 
The position of the BF can be conveniently expressed with respect to the J2000 reference frame as function of the three Euler's angles θ, ϕ and ψ. 
The nutation angle θ is the angle between ẑ^E and ẑ^b, the spin angle ψ is the angle between the nodal line and x̂^b axis, the precession angle ϕ is the angle between   x̂^E axis and the nodal line, where the nodal line is the intersection between the two planes (x̂^E,ŷ^E) and (x̂^b,ŷ^b). 
The transformations of a vector  V from J2000 to BF reference frame  V^b =ℛ  V^E  is given by the rotation matrix ℛ, that in terms of Euler's angles is (<cit.>):


  
	
    ℛ=(
    	[   cosϕ cosψ - cosθ sinϕ sinψ   sinϕ cosψ + cosθ cosϕ sinψ                    sinθ sinψ; - cosϕ sinψ - cosθ sinϕ cosψ   cosθ cosϕ cosψ - sinϕ sinψ                    sinθ cosψ;                    sinθ sinϕ                  - cosϕ sinθ                         cosθ ]).

 

The satellite motion along its orbit, supposed to be quasi-circular with radius a, will be identified by the Keplerian orbital parameters measured in the J2000 frame: Ω, the right ascension of the ascending node, I, the orbit inclination over the Earth's equator, ω, the argument of pericenter, and M_0, the mean anomaly.



 §.§ Equations of motion

The spin evolution can be conveniently described by the Euler's equations


    I_x ω̇_sx^b-ω_sy^b ω_sz^b (I_y-I_z)=M_x 
    
    I_y ω̇_sy^b-ω_sx^b ω_sz^b (I_z-I_x)=M_y
    
    I_z ω̇_sz^b-ω_sx^b ω_sy^b (I_x-I_y)=M_z
 

where I_x,I_y,I_z  are the moments of inertia with respect to the principal axes of inertia of the satellite body, while M_x, M_y and M_z are the components of total torque M along the same axes. 
The components of the angular velocity measured with respect to the body axes, ω_sx^b, ω_sy^b and ω_sz^b can be substituted by their expressions in term of the Euler angles in the J2000 reference frame <cit.>:

    ω_sx^b   =   ϕ̇ sinθsinψ+ θ̇ cosψ
    ω_sy^b   =   ϕ̇ sinθcosψ- θ̇ sinψ
    ω_sz^b   =   ϕ̇ cosθ + ψ̇.


We finally obtain the Euler's equations in term of Euler's angles measured in the J2000 reference frame:





    θ̈   =   M_x/ I_xcosψ - M_y/ I_ysinψ - I_z/I_yϕ̇ψ̇sinθ  +I_y-I_z/I_xϕ̇^2 sin(2 θ)/2
           +I_x-I_y/I_xΛ/ I_y[θ̇(ψ̇+ϕ̇cosθ)  sin(2 ψ)/2 +ϕ̇^2 sin(2 θ)/2sin^2 ψ - ϕ̇ψ̇sinθ(  I_y -  I_z/Λ - sin^2 ψ) ] 
    
    ϕ̈   =   M_y/I_ycosψ/sinθ+ M_x/I_xsinψ/sinθ +I_z/I_yψ̇θ̇/sinθ - Λ/I_xcosθ/sinθϕ̇θ̇  + 
          I_x-I_y/I_yΛ/I_x[ 1/sinθ( sin^2 ψ - I_x/Λ)ψ̇θ̇
    	-  sin(2 ψ)/2( cosθϕ̇ + ψ̇)  ϕ̇ 
    	- cosθ/sinθϕ̇θ̇cos^2 ψ]  
    
    ψ̈   =   M_z/I_z - cos(θ)/sin(θ) (M_y/I_ycos(ψ) + M_x/I_xsin(ψ)) + ϕ̇θ̇1/sinθ( I_y -  I_z/ I_x cos^2 θ + 1) -  I_z/ I_y ψ̇θ̇cosθ/sinθ+ 
          ( I_x -  I_y ) [1/ I_zϕ̇θ̇1/sinθ(sin^2 θcos(2 ψ)+ Λ/ I_x  I_ycos^2 ψcos^2 θ)- θ̇^2 1/2I_zsin(2ψ)-ϕ̇^21/2I_zsin(2 ψ) . 
          .(Λ I_z/ I_x  I_ycos^2 θ -sin^2 θ)  - ψ̇θ̇1/ I_ycosθ/sin(θ)(Λ/I_xsin^2 ψ - 1 )+ Λ/2 I_x  I_yϕ̇ψ̇cosθsin(2 ψ) ],

	where Λ= I_x+I_y -  I_z.



These equations were written without hypothesizing any peculiar symmetry for the satellite, that is we assumed I_x ≠ I_y ≠ I_z differently from <cit.>, who instead have adopted an axial symmetry for the satellite. 
Coherently,  our expressions (<ref>), (<ref>) and (<ref>) converge, in the limit of I_x=I_y (cylindrical symmetry), to expressions (9)-(14) of <cit.>. 
We underline that even if in this limit the formal expressions are equal, our Euler angles are expressed in the J2000 reference, while <cit.> adopted as main reference frame the orbital one.


In our analysis, for the torques to be included in equations (<ref>)-(<ref>) we considered four contributions: 


	
  * M_mag, i.e. the torque from the Earth magnetic field (section <ref>)
	
  * M_grav, i.e. the torque from the Earth gravitational field  (section <ref>)
	
  * M_off, i.e. the torque  arising from radiation pressure acting on the geometric center of the satellite,  if it does not coincide with its center of mass (section <ref>)
	
  * M_asy, i.e. the torque from radiation pressure acting on the two opposite hemispheres of the two LAGEOS satellites, which have a different reflectivity (section <ref>)


therefore 
M=  M_mag+  M_grav+  M_off+  M_asy.



 §.§ Torques in time domain
 
In the following sub-sections we will calculate the expressions of the main torques, to be inserted in previous equations (<ref>), (<ref>) and (<ref>), as function of time. 



  §.§.§ Torque from the Earth magnetic field

Both LAGEOS satellites (<cit.>) and LARES (<cit.>) are made of conductive material, and even if they are supposed to be uncharged, they assume a magnetic moment while they are spinning inside the Earth magnetic field, the value and the direction of which change along their orbit. The induced magnetic moments m, because of the Earth magnetic field B, produce a torque on the satellites (M_mag=m×B).

In the previously quoted papers, to describe such a phenomena, the LAGEOS satellite was modeled (<cit.>) as a conductive sphere rotating in a static magnetic field. The value of this constant magnetic field was calculated averaging the magnetic field of the Earth over the entire orbit of the satellite. The solution to this problem is well known in the literature, see for instance <cit.>.   

This solution, which is completely valid in a quasi-stationary field, can be suitably used as long as the the rotation period of the satellite is much shorter than its orbital period as well as of the Earth's rotation period, but it could produce wrong results when is used in slow-spin conditions, as in <cit.>. 
These authors have probably adopted this simplification due to the lack of an handy model for a sphere rotating in an alternating magnetic field. The only solution apparently available till now in the literature is the one by <cit.> that unfortunately it is not easy to deal with.

In order to obtain a more general expression of the magnetic torque we faced the problem to find an easily integrable expression for the torque acting on a conducting sphere rotating in an alternating magnetic field. 
We found an expression for the torque that was applied to the satellites (modeled like perfect conducting spheres) while rotating and moving along their orbit:

 
    [                                                                                                                   M_mag^E=                                  V ∑^8_i=0|  B_i |^2/2 | ω_s|{-A”_i[1+cos(2 ω_i t+2 φ_i) ]+D'_i sin(2 ω_i t+2 φ_i) }ω_s^E+;                                                                                                                           ;                                                                                                                            V  ∑^8_i=0 B_i·ω_s/2 | ω_s |^2{[α'(ω_i) -A'_i] [1+cos(2 ω_i t+2 φ_i) ]-[D”_i+α”(ω_i) ] sin(2 ω_i t+2 φ_i) }( ω_s^E× B_i) +;                                                                                                                           ;                                                                                                                                                              V  ∑^8_i=0 B_i·ω_s/2 | ω_s |{A”_i [1+cos(2 ω_i t+2 φ_i) ]-D'_i sin(2 ω_i t+2 φ_i) } B_i, ]

	
    [                                                                                                                     ; A'_i=α'(ω_s^E -ω_i) +α'(ω_s^E + ω_i)/2                                         D'_i=α'(ω_s^E -ω_i)-α'(ω_s^E + ω_i)/2;                                       ; A”_i=α”(ω_s^E -ω_i) +α”(ω_s^E + ω_i)/2                                        D”_i=α”(ω_s^E -ω_i)-α”(ω_s^E + ω_i)/2.;                                        ]

 


In the above expressions, 
α(ω)=α'+jα” is the complex fourier transform of the magnetic polarizabilty per unity of volume of the satellite, V is the satellite volume,  ω_i are the angular velocities of the harmonic components of the magnetic field, while ω_s^E is the satellite spin angular velocity in J2000 frame, whose components can be expressed in term of the Euler angles <cit.>:


    ω_s x^E   =   θ̇ cosϕ+ψ̇ sinθsinϕ
    ω_s   y^E   =   θ̇ sinϕ-ψ̇ sinθcosϕ
    ω_s   z^E   =   ψ̇ cosθ + ϕ̇.
 In appendix A we calculated  the magnetic field experienced by the satellite as function of time along its orbit using the dipole approximation for the Earth's magnetic field (see appendix A for explicit expressions of B_i):


    B   =   ∑^8_i=0 B_icos(ω_i t+φ_i),

where

    ω_0   =   0 
    ω_1   =   ω_2=ω_⊕-2 n  
    ω_3   =   ω_4=ω_⊕+2 n 
    ω_5   =   ω_6=2 n 
    ω_7   =   ω_8=ω_⊕,

and

    φ_i   =   {[ -π/2     ;    0    . ].



In the previous expressions ω_⊕ represents the Earth's rotational angular velocity, while n represents the satellite mean motion in the hypotheses of a circular orbit having radius a:


    n=√(G M_⊕/a^3)(sin I sinΩ, -sin IcosΩ, cos I),
 

where M_⊕ and G are, respectively, the Earth's mass and the gravitational constant.


The torque that originated from Earth magnetic field
(<ref>) shows components at the magnetic field harmonics (<ref>) and their multiples. The torque is expressed as sum of several contributions: i) along the satellite angular speed direction (ω_s^E), ii) along the magnetic field directions B_i  and iii) along directions orthogonal to the previous ones ( ω_s^E× B_i). 

The expression (<ref>) is calculated in J2000 frame, the matrix ℛ transforms it into the equivalent expression in body frame, to be used in the equations of motion (<ref>): M_mag^b=ℛ  M_mag^E.  We have not determined an analytic expression for M_mag^b, because it is calculated numerically while solving the Euler's equations. 
The complex fourier transform  of the magnetic polarizability of the satellite α(ω) can be approximated with the expression valid for a perfect sphere: we adopted the one included in <cit.> calculated for a sphere having radius R, electrical conductivity σ, in a magnetic uniform field changing with an angular frequency ω:

 
	
    α(ω)   =   α'+j α”= 
       =   3/8 π{2 μ_r [1-k ·(k)]+[1-k^2-k ·(k)]/μ_r [1-k ·(k)]-[1-k^2 -k ·(k) ]},

	where
	
    k(ω)=R/δ(ω)(1+j), 
    δ(ω)=c/√(2 πωσμ_r ).

This expression contains the dependence from  the relative magnetic permeability μ_r, that was  ignored  by all the other authors (<cit.>). Even if μ_r can be very close to 1 (for instance for aluminum μ_r-1≅2.2· 10^-5), its contribution could be non negligible at very low frequencies, where it represents the residual magnetic polarizability of a non-spinning satellite.  

Past experience with the two LAGEOS  (<cit.>) has shown that the calculated values of the time evolution of the spin of LAGEOS satellites fit well the experimental data only if a low frequency approximation is chosen for  (<ref>):

 

    α(ω)   =   α'+j α”= 
    
       ≃   [3/4 πμ_r-1/μ_r+2 -9/350 πμ_r (μ_r+9)/(μ_r+2)^3(R/δ)^4 ]β'  +j 9/20 πμ_r/(μ_r+2)^2(R/δ)^2 β”,



β' and β” are dimensionless constants to be determined experimentally in order to take in account the differences of the satellite's shape  with respect to an ideal sphere. 

There is no trivial explanation to justify such a simplified empirical model. It was hypothesized that the composite structure of LAGEOS, an inner cylindrical  core made of brass inserted in two hollow aluminum hemispheres, would invalidate the assumption of a uniform sphere as a model for LAGEOS. We will devote a forthcoming paper to discuss in detail this case that, for now, we adopt as a matter of fact. In the case of LARES, made of a unique block of metal, we believe that it is more correct to apply the expression (<ref>).

The first term of the real part of (<ref>) is introduced for the first time in this work, and it is an important term, because it represents the residual polarizability at zero frequency. This term dominates the real part of the polarizability at frequencies lower than ω_c:


    ω_c=5 c^2 /12 π R^2 σ(μ_r+2) /μ_r√(42(μ_r-1)/μ_r (μ_r+9 )).


This bound frequency corresponds for the two LAGEOS  to a period of about  T_c=2 π/ω_c∼400 s, that was already reached by both satellites. 


We underline that there is a big difference between the formula we have used for the torque of magnetic origin (<ref>), and that valid for static magnetic fields used till now by the others authors (see for instance (2), (5) and (6) of <cit.>).  
In formula (<ref>), the torque along the satellite angular velocity ω_s could change sign when 
ω_s= ω_i, depending on  A_i” that switch from negative to positive. In this case the angular velocity ω_s could be maintained equal to ω_i by the magnetic force. 
This resonant condition, however, appears only if the ratio R/δ (ω_i) is greater than a threshold value that is not reached anyway by LAGEOS or LARES. A more detailed discussion around these very important issues will be included in a forthcoming work.

We point out that two possible exponents can be used in the Fourier transform  (i.e.  jω or -j ω) and, consequently, we can have two different expressions with different complex conjugate for α. One must be careful to maintain the same convention in the inversion formula in order to anti-transform in the time domain.     





  §.§.§ Torque from the Earth Gravitational field

In this section we  calculate the components of the torque of gravitational origin M_grav, also known as gravity gradient torque, to be included in equations (<ref>), (<ref>) and (<ref>). In fact, because the satellites do not have a perfect spherical symmetry in their mass distribution, the Earth gravitational field produces a torque on them. We restrict our analysis to the monopole component of this field, whose general expression is given by (<cit.>):


 

    M_grav^b=3 n^2 { ŝ^b×[ I_x ( ŝ^b·x̂^b) x̂^b  +I_y (  ŝ^b·ŷ^b) ŷ^b +I_z ( ŝ^b·ẑ^b) ẑ^b] } =3 n^2
    	[
    	[   (I_y-I_z) s_y^b s_z^b;   (I_z-I_x) s_x^b s_z^b; (I_x-I_y )  s_x^b s_y^b;                         ]],

 

where ŝ^b is the unit vector Earth-to-satellite, having components s_x^b, s_y^b, s_z^b in the BF reference frame. 

To be conveniently included in the equations of motion (<ref>-<ref>), (<ref>) must be expressed as a function of the Euler angles with respect to the J2000 reference frame.

The expression of the three components of the Earth-to-satellite unit vector can be  easily calculated in the J2000 reference frame (see for instance <cit.>) and  can be further converted into the BF reference frame using (<ref>) for the rotation matrix ℛ, finally obtaining:

 

    s_x^b   =   cos(ω+M_0+n · t) [cos( ϕ- Ω) cos( ψ) -  cos(θ) sin( ϕ- Ω) sin( ψ) ] + 
          sin(ω+M_0+n · t) {sin( ϕ- Ω)   cos(I) cos( ψ)+[ sin(θ) sin(I)+
    	cos(θ) cos( ϕ- Ω) cos(I)  ]  sin(ψ)  }
    
    
    	s_y^b   =    - cos(ω+M_0+n · t)[  cos(θ) sin( ϕ- Ω) cos( ψ) +cos( ϕ- Ω) sin( ψ) ] + 
          sin(ω+M_0+n · t) {[  sin(θ)  sin(I)+ cos(θ) cos( ϕ- Ω)cos(I) ]  cos( ψ)-sin( ϕ- Ω) cos(I) sin( ψ)   }
    
    
    	s_z^b   =   cos(ω+M_0+n · t) sin(θ) sin( ϕ- Ω)+ 
          sin(ω+M_0+n · t) [cos(θ) sin(I) -
    	sin(θ) cos( ϕ- Ω) cos(I) ].

 


Obviously, from equation (<ref>) follows that the gravitational torque is zero if the satellite has a perfect spherical symmetry (i.e. I_x=I_y=I_z). Moreover, the unit vectors of Eqs. (<ref>) are periodic with the satellite mean motion frequency n, therefore their product in the expression (<ref>) for the gravitational torque produces a constant component plus one periodic component with a frequency 2n.




  §.§.§ Offset torque
 
The radiation pressure on the satellite surface acts, at the end, on the geometric center of the satellite. If the geometric center does not coincide with the center of mass of the satellite, the radiation pressure produces a torque (<cit.>). If the vector h^b identify the geometric center position respect to the center of mass in the body reference frame, the torque expression is:


    M_off^b   =   νπ R^2 Φ_⊙/c C_R (  h^b ×ŝ_⊙^b).


In this equation R is the radius of the satellite, C_R represents its radiation coefficient (related with the reflectivity of its surface), Φ_⊙ is the solar flux at the Earth-Sun distance, ŝ_⊙^b is the satellite-to-Sun unity vector, finally ν (with 0≤ν≤ 1) represents the shadow function. This function measures the percentage of solar flux that reaches the satellite during the eclipses by Earth (with ν=1 when there is not an eclipses, and ν=0 when the eclipse is complete; see for instance <cit.>).  
The expression (<ref>) is more general of the one by <cit.>, that hypothesized the coincidence between the rotation and symmetry axes.

To include the torque coming from the offset (<ref>) into the equations of motion  (<ref>-<ref>), the unit vector ŝ_⊙^b
must be expressed in the BF reference frame in term of the Euler angles with respect to the J2000 reference frame. 

If the motion of the Sun in the J2000 reference frame is considered circular with constant angular velocity ω_⊙ along the ecliptic, by denoting with  ϵ the angle between the plane of the ecliptic and the Earth's equator, and with λ_⊙ the Sun's ecliptic longitude at t=0,  the three components s_⊙ x^b, s_⊙ y^b and  s__⊙ z^b of the  Earth-to-Sun unit vector in the BF reference frame are:

 

    s_⊙ x^b   =   sin(λ_⊙ + ω_⊙· t) [ cos(ϵ) (sin( ϕ) cos( ψ)   +  cos( θ)cos( ϕ)    sin( ψ)) +  sin(ϵ)    sin( θ)sin( ψ) ] + 
          cos(λ_⊙ + ω_⊙· t) [ cos( ϕ)  cos( ψ) -  cos( θ)  sin( ϕ)  sin( ψ) ] 
    
    	s_⊙ y^b   =      -  sin(λ_⊙ + ω_⊙· t) [ cos(ϵ) ( sin( ϕ)  sin( ψ) -  cos( ϕ)  cos( ψ)  cos( θ)) -  cos( ψ)  sin(ϵ)  sin( θ) ] - 
          cos(λ_⊙ + ω_⊙· t) [ cos( ϕ)  sin( ψ) +  cos( ψ)  cos( θ)  sin( ϕ) ] 
    
    	s_⊙ z^b   =   sin(λ_⊙ + ω_⊙· t) [ sin(ϵ)  cos( θ) -  cos(ϵ)  cos( ϕ)  sin( θ) ] +   cos(λ_⊙ + ω_⊙· t)  sin( θ)sin( ϕ).

 

This torque, slowly changing with a period of one solar year, is modulated by the shadow function ν.





  §.§.§ Anisotropic  reflection torque
 
A difference in the reflectivity of the two hemisphere of LAGEOS was hypothesized by <cit.> in order to explain the observed along-track residuals in acceleration of the satellite. This empirical result was reasonably explained by <cit.> as due by a possible asymmetry in the satellite's reflection of the visible solar radiation introduced by the four germanium CCRs. If the relative difference between the reflectivity of North and South hemispheres is Δρ=(C_R^N-C_R^S)/C̅_R, with  
C̅_R=(C_R^N+C_R^S)/2, the following torque is produced (<cit.>):

 

    M_ar^b   =   ν 2/3R^3 Φ/c Δρ C_R   ( ẑ^b ×ŝ_⊙^b )| ẑ^b ×ŝ_⊙^b |=ν 2/3R^3 Φ/cΔρ C_R √((s_⊙ x^b)^2+(s_⊙ y^b)^2)[
    	[ -s_⊙ y^b;  s_⊙ x^b;        0 ]] 
    
    	.

 

The difference in reflectivity between the two hemispheres has as consequence the non-coincidence of the center of mass of the satellite with the point of application of the solar radiation pressure. 
Therefore, we could have described the effect within the torque (<ref>) analyzed in previous section (<ref>).  However, we preferred to maintain the expressions for two distinct torques, as done by previous authors, in order to better connect the torque to the measurable parameters, as Δρ for this second torque.

As we can see from Eq. (<ref>), also this torque is slowly changing with a period of one solar year and it is modulated by the shadow function ν. 






 §.§ The solution with averaged equations


The numerical solution of the equations (<ref> - <ref>) can be more quickly found, if averaged expressions are adopted for the torques.  This simplification does not modify the result until the satellite spins with a time scale shorter (or comparable) than the orbital period of the satellites's and of the Earth rotational period. The averaged values for the torques are calculated by integration on these longer periods.
In this case, both  the magnetic  (<ref>)  and the gravitational (<ref>) torques can be averaged since these torques vary with time scales that are related with the pulsation of the satellite mean motion n and of the Earth's angular speed ω_⊕, as well as with their combinations. 



  §.§.§ Averaged Magnetic Torque


The averaged magnetic torque can be easily calculated starting from expression (<ref>), and integrating each one of the harmonic components over its characteristic period. Therefore, in this limit, the variables introduced in (<ref>) become: D'=D”=0, A'_i= α' and  A”_i= α”. 
For the magnetic (averaged) torque expression we obtain:

 

    < M_mag>   =    -V <|  B|^2> α”(ω_s) ω_s/| ω_s| + V α'(0) -α'(ω_s)/| ω_s|^2 <(  B·ω_s) (  ω_s× B) > +  V α”(ω_s)/| ω_s|<(  B·ω_s)    B>= 
       =   -V  α”(ω_s) ω_s/| ω_s|(<|  B|^2>-<  B  B^T > )+ V α'(0) -α'(ω_s)/| ω_s|^2 ω_s×(<  B  B^T>  ω_s),

 

where we used the two relations:


    <(  B·ω_s) (  ω_s× B) > = ω_s×(<  B  B^T>  ω_s),



    <(  B·ω_s)    B> =<  B  B^T >  ω_s.
    

By averaging on the nine frequency characteristic of the magnetic field, see Eq. (<ref>), we get the different terms to be inserted in (<ref>), that is:

 
	
    ⟨ B_i B_j ⟩(a^6/d^2) = α_i,j   =   1/8·[9 sin^2I+1/2sin^2θ_p( 20-27 sin^2I)]δ_i,j
          -1/8(1-3cos^2θ_p )E_iE_j
          +3/8cos I (1-3cos^2θ_p)(E_in̂_j+E_jn̂_i ) 
          -3/8n̂_i n̂_j[3(1-3cos^2I)-1/2sin^2θ_p(5-27cos^2I)],

	
	
    ⟨ B^2 δ_i,j - B_i B_j ⟩(a^6/d^2)   =   β_i,j={1+1/8·[3 sin^2I+1/2sin^2θ_p( 4-9 sin^2I)]}δ_i,j
          +1/8(1-3cos^2θ_p ))E_iE_j  
          -3/8cos I (1-3cos^2θ_p ) (E_i n̂_j+E_j n̂_i ) 
          +3/8n̂_i n̂_j[3(1-3cos^2I)-1/2sin^2θ_p(5-27cos^2I)],

	
	
    <|  B|^2>=<∑^8_i=0|  B_i |^2>   =   ∑_1^3α_i,i=-1/4[ 10-6 cos^2I-sin^2θ_p (3-9cos^2I)],

 

where
E=[0,0,1] is the unit vector along the Earth rotation axis, δ_i,j is kronecker's delta and n̂ is the  mean motion unit vector  normal to the orbital plane of the satellite.

We introduced α_i,j and β_i,j to use  the same notation of expressions (7a) and (7b) of <cit.>. 
The first line of Eq. (<ref>) differs from the analogous expression (7a) and also the last line of Eq. (<ref>) differs from the correspondent line of (7b). 
We are confident in our calculations, that follow a very different path with respect to those of  <cit.>, and we believe that the expressions (7a) and (7b) by <cit.> contain some small errors, as they do not satisfy some obvious conditions. In fact, β_i,j+α_i,j has to be diagonal and equal to ∑_1^3α_i,i. Conversely, these conditions are satisfied by our expressions.



  §.§.§ Averaged Gravitational Torque


The gravitational torque (<ref>) has one constant component and one varying at a frequency twice of the satellite mean motion n. Integrating on the corresponding period is therefore possible to calculate the averaged torque, we obtain:  

 
	
    <  M_grav>   =   3/32 n^2
    	| 
    	[   (I_y-I_z)( M_1cos( ψ) - M_2sin( ψ));    (I_z-I_x)(M_1sin( ψ) + M_2cos( ψ)); (I_x-I_y)(-M_3sin(2 ψ) - M_4cos(2 ψ)) ]|,

	where	
	
    M_1   =   2 {3 cos(2 I)+cos(2 ϕ-2 Ω) [cos(2 I)-1]+1}sin(2 θ)-8 sin(2 I) cos(ϕ-Ω) cos(2 θ) 
    
    	M_2   =   4 [cos(2 I)-1] sin(2 ϕ-2 Ω) sin(θ)-8 sin(2 I) sin(ϕ-Ω) cos(θ)
    
    	M_3   =   [1+3 cos(2 I)] [cos(2 θ)-1] +cos(2 ϕ-2 Ω) [cos(2 I)-1] [cos(2 θ)+3]+4 sin(2 I) cos(ϕ-Ω) sin(2 θ)
    
    	M_4   =   4 sin(2 ϕ-2 Ω) [cos(2 I)-1] cos(θ)+8  sin(ϕ-Ω) sin(2 I) sin(θ).

 




§ SOLUTIONS OF THE EQUATIONS FOR THE TORQUES


The equations (<ref>-<ref>) cannot be integrated analytically. We built a code based on MATLAB routines in order to solve for the differential equations. 
We used as independent variables the three Euler angles θ, ϕ,  ψ and their time derivatives  θ̇,  ϕ̇ and ψ̇, in such a way to transform the solution into that of a system of six first-order differential equations.



 §.§ Physical quantities and initial conditions


A very critical point in solving the equations is in choosing the most likely initial conditions and in defining, with the highest possible precision, the parameters that appear in them. 

The moments of inertia I_x, I_y,I_z of the two LAGEOS were not directly measured before their launch,  but were estimated with a careful analysis by <cit.>. Unfortunately, we could not found any document containing the values for the moments of inertia of LARES.  In this regard, <cit.> computed the moment of inertia of LARES in the hypothesis that the satellite is an homogeneous and symmetrical sphere (with no oblateness). They obtained I = 5.125 kg · m^2. 
With an approach a similar to the one adopted for the two LAGEOS in <cit.>, we built also for LARES a 3D-model and  we calculated the values reported in Table <ref>. The reported error is calculated considering that the structure of LARES was designed with the possibility to add small masses to tune the moments of inertia. In the simulations we opted for a  very small oblateness, but a spherical symmetry is nevertheless possible within the errors. 






In Table <ref> we report the values of the electromagnetic quantities we used in our analyses and simulations. We started from the range of possible changes of the values reported in <cit.>, modifying them within their errors with the aim to minimize the residuals between the predictions of our model and the available observations. 





In Table <ref>, the overall optical properties of the surface of each satellite are  reported. These values are the ones that we found to be the most likely in the literature, and that we confirmed with an independent analysis.





The data of Table <ref> relative to the initial conditions for the spin vector of the two LAGEOS satellites are taken from <cit.>, and modified within their variability interval to fit the experimental data.
Finally, in Table <ref>, the values of the keplerian parameters of the orbit of the three satellites are shown. These values were calculated averaging  the output values from a precise orbit determination (POD) made using the GEODYN II software (<cit.>).












 §.§ Numerical solution and experimental data

In this section we present our results for the time evolution of the spin of the two LAGEOS and LARES satellites. This evolution has been calculated using the equations of motion previously introduced, namely Eqs. (<ref>), (<ref>) and (<ref>).

The expressions that we used for the torques were of two types: the more general one, that we have discussed in section <ref>, and the averaged one (see section <ref>). This last solution differs from the general one for the magnetic and gravitational torques. For the other two torques that we considered, that vary with the periodicity of one year and are locally modulated by the Earth's shadow function, we adopted the general expressions (<ref>) and (<ref>). 

The parameters and the initial conditions we used are those reported in Tables 
<ref>–<ref>. 
For each satellite it has been calculated and drawn the period of rotation and the orientation (right ascension and declination) of the spin as function of time. 
In Figures (<ref>) and (<ref>) we show the results for the older LAGEOS, while in  Figures (<ref>), (<ref>) and in Figures (<ref>), (<ref>) we show, respectively, the results in the case of LAGEOS II and for the newly LARES.  

















In all these plots the results for both the general model (red continuous line) and the averaged model (blue dashed line) are shown. These results are compared with all the available measurements. 
All these comparisons validate very well our models. 
As specified at the beginning of section <ref>, we refer to our new general model as to LASSOS, that is the acronym of LArase Satellites Spin mOdel Solutions.

The bulk of the experimental measurements are of two kind for the two LAGEOS: those obtained using Sun's light reflected by the CCR mirrors (see <cit.>), and those obtained observing the modulation of the laser light reflected back by the satellite's CCR (see <cit.>). The last measurements are indeed available for all three satellites (see <cit.>).
The latter method has increased a lot in efficiency after the introduction of high frequency (kHz) repetition laser in several SLR stations, the first of these stations has been in 2003 Graz.

In the comparison between models and observations, a disagreement is present for LAGEOS's declination in the years 1983-1987. The reason for this disagreement with the data analyzed by <cit.> could come from the still low quality of the laser used by the ILRS stations in that years. 
Another disagreement is present in the case of LAGEOS in the  measurements after 1996, and for LAGEOS II after 2008. These periods correspond to the decrease of the spin frequency of the satellites, close to the limit of the possible measurements with kHz laser. 
We notice that the same disagreements between experimental measurements and model are present when the LOSSAM model is applied, see Figure 1  of <cit.>.

Using our model LASSOS we have also computed the time evolution of the absolute value and of the components of the different torques acting on the satellites in the J2000  reference frame.  
For instance, in the case of LAGEOS in Figure (<ref>) we show the modulus of the torques due to the Earth's gravity gradients and magnetic field, both when they are averaged and when they are not averaged, while in Figure (<ref>) the values of the two torques due to radiation pressure are shown for the same period. The amplitude of these torques is modulated by the eclipses, that bring these torques to zero value.
During the first period the dominating torque is the magnetic one, later it is overcome by the  torque from gravity gradients. In the last years, the effect of the torques that arise from  radiation pressure becomes comparable with that of the other two torques.








Finally, in Figure (<ref>) and (<ref>),  the three components of the considered torques on LAGEOS during the first 20 years of its life, and in the body reference frame, are shown.
The torque due to the Earth magnetic field is the only one acting on the ẑ^b axis with a despinning effect. If the geometric center, where the radiation pressure force is applied,  has coordinates h_x or h_y different from zero, a further component of the torque along the  ẑ^b axis will appear; in this case the torque could also have an effect in increasing the spin rate of the satellite.

Similar considerations to those obtained in the case of LAGEOS are also valid for LAGEOS II and LARES.










§ DISCUSSION ABOUT THE RESULTS


In this work we have faced the problem  of the knowledge of the spin evolution of LAGEOS-like satellites. Our first goal in this field was to significantly improve the  time evolution of the spin vector of the passive LAGEOS, LAGEOS II and LARES satellites.
The knowledge of the right spin evolution for these satellites, both in rate and orientation, is of fundamental importance in order to correctly model the thermal thrust forces acting on their surface, as well as the asymmetric reflectivity observed in the case of the two LAGEOS.
The ability to manage these non-gravitational perturbations through reliable models is, consequently, of utmost importance. This not only for studies focused on measurements of fundamental physics by means of the analysis of the orbit of these satellites, but also for the potential benefits in the fields of space geodesy and geophysics.
 
For instance, in this regard, the ILRS determinations of the tracking station coordinates and of the Earth Orientation Parameters (EOP), represent the starting point for the determination of the Earth's geocenter and the definition of the Earth's International Terrestrial Frame (ITRF).
These, and other geophysical products, require a careful orbit determination of the two LAGEOS and (now) LARES (and other) satellites. Indeed, a POD should be based on the continuous refinement of the perturbations included in the dynamical model of each satellite.



In Section <ref> we have introduced the main torques acting on the two LAGEOS and LARES satellites. Of course, the four effects there considered are not the only physical effects able to produce a torque on a spacecraft in orbit around our planet.
However, these torques are the most important to take into account at the current level of the orbit modelling of the LAGEOS and LARES satellites and also of the their corresponding POD, that depends not only by the models implemented in the software but also from the tracking observations and the Reference Frames accuracy. In fact, the other possible torques that we can take into consideration are negligible, at least by two or three order of magnitude, with respect to the torque due to the asymmetric reflectivity and the one due to the offset between the satellite center of mass and the center of pressure (i.e. the geometric center), that is  the current smallest torques that have a maximum magnitude of about 10^-9 N m, see Figure <ref>.

For these reasons we have not performed a complete analysis to model torques that 𝑎 𝑝𝑟𝑖𝑜𝑟𝑖 should be orders of magnitude smaller than those considered, but we simply estimated the order of magnitude of the effects not included in our model. In particular, we evaluated the order of magnitude of five additional torques originated by:

	
  * Neutral drag 
	
  * Charged drag 
	
  * Earth Yarkovsky 
	
  * Yarkovsky-Schach  
	
  * Inner eddy currents




We show in Table <ref> the order of magnitude of the acceleration and of the corresponding torque produced by these physical effects in the case of LAGEOS. These orders of magnitude are still valid in the case of LAGEOS II. The first four torques are in principle null for a satellite with a perfectly spherical in shape mass distribution. Deviations from this distribution are due to an offset between the geometrical center of the satellite and its center of mass. In the estimation of the torques that arise from these effects a very conservative value of 1 mm for the offset of the three satellites was adopted. The margin adopted in the estimation can be evaluated comparing the value of the offset of  1 mm with those used in the current work that were estimated on the basis of our fit of the available spin observations of the satellites, see previous Table <ref>.   

The accelerations acting on the satellite  in the case of neutral drag  have been obtained from <cit.>, and from <cit.> — as a very pessimistic upper bound — in the case of the drag produced by the Coulomb interaction of the charged spacecraft with the surrounding charged particles.
The thermal drag acceleration produced by the Earth-Yarkovsky effect has been obtained as an upper limit of the acceleration estimated by <cit.>, while in the case of the solar Yarkovsky-Schach effect the value is the one estimated by <cit.>.

These last two torques have to be considered as torques related to a thermoelastic deformation of the satellites, and are operative only if the temperature distribution across the satellite surface is effective in changing the spacecraft mass distribution.
Finally, the last effect is related with the electrodynamical force produced by the eddy currents generated by the photoelectric emission of the satellite, see also <cit.> and <cit.>.


In Table <ref> the order of magnitude of possible additional torques are shown in the case of LARES. Among the non-gravitational perturbations acting on this satellite, the effects of the neutral atmosphere drag need a special attention, since the relatively low orbit of the satellite, at an height of about 1,450 km  where the residual Earth atmosphere cannot be neglected (see <cit.>).





The neutral drag produces on LARES an acceleration of about 1 × 10^-11 [m/s^2] and a torque, in the presence of the cited offset, of about 4 × 10^-12 [N m], a factor ≈ 40 larger than that obtained in the case of the two LAGEOS, but still negligible with respect to the smaller torque of those considered within the LASSOS model.

The last line of Table <ref> accounts for the possible cumulative contribution to the torque in the presence of an offset of about 1 mm from the other smaller perturbing effects.




§ CONCLUSIONS AND FUTURE WORK


The work here presented falls in the LARASE research program, whose ultimate goal is to provide refined measurements of relativistic physics on the orbit of laser-ranged satellites, as the ones considered in the present paper: the two LAGEOS and LARES.
Indeed, a very important prerequisite in order to reach refined measurements of the tiny gravitational effects on the orbit of a satellite, as those predicted by Einstein's theory of General Relativity, is the inclusion of accurate models for the handling of both gravitational and non-gravitational perturbations in the data reduction of the orbit of a satellite.
In particular, the non-gravitational perturbations are responsible of very subtle effects on the orbit, really quite complex to model in a reliable fashion.
Their modeling will first impact the orbit determination of the satellites considered for the relativistic measurements of interest and, finally, will have a deep impact in the robustness of the final error budget of the measurements.

In the present work we have described the effects of the torques due to: 

	
  * Earth's magnetic field
	
  * Earth's gravitational field	
	
  * not coincidence between geometrical center and center of mass of the satellite
	
  * asymmetries in the reflectivity of the satellite surface.
    

The model for the spin presented in this paper, i.e. LASSOS, is, to our knowledge, the only model — among all the models developed so far — that has succeeded at the same time in answering to all the following aspects:


	
  * it implements not-averaged torques
	
  * it can be used for any orientation of the satellite rotation axis
	
  * it is valid for any spin rate
	
  * it has produced results compatible with the available observations.



The good agreement with the experimental measurements is also consequence of the independent evaluation of the moments of inertia of the satellites (<cit.>). 
A further innovation is the introduction of a description of the torque from Earth's magnetic field using its local values along the satellites' orbit, and not simply by means of a static average value for the magnetic field. Indeed, this represents a very important result.

A step forward of LASSOS with respect to previous models is also represented by the expression that we used for the magnetic polarizability per unit of volume α of the satellite.
In particular, we considered the dependence of α from the relative magnetic permeability μ_r, that provides the residual magnetic polarizability in the case of a non-spinning satellite.
This aspect is quite important for the considered satellites. Indeed, the rotational period of LAGEOS is quite long, probably of the order of about 1 day, while for LAGEOS II it is probably of the order of its orbital period. In the case of LARES the rotational period is probably of a few thousands of seconds, but its slowing down is much faster than that of the two LAGEOS.
An average of our solution on the day and the orbital period, in such a way to reproduce the limit of the spin evolution valid in the so-called fast-spin approximation, has allowed us to point out also some errors present in previous works, which are based on simplified averaged models.

In order to further improve (if possible) the knowledge of the parameters that enter in our spin model, in the near future we plan to develop a more sophisticated tool regarding the fit procedures that we used for the comparison of the available observations (i.e. the measurements of the orientation and rate) of the spin and  the corresponding predictions by the LASSOS model.
We retain that this tool will be fundamental in order to increase the robustness for the prediction of the spin evolution in the periods not covered by measurements, and during which the spin rate is so low that only our spin model, and not one model based on averaged equations, is able to provide a reliable foresight.
Of course, the possibility of comparing the predictions of LASSOS with new observations of the spin would be very interesting. Considering the techniques applied till now for the determination of the spin of the satellites, this is not an easy task in the case of the two LAGEOS, since their very low rotational rate, but it is still possible for LARES with a spectral analysis of the range data obtained by means of kHz lasers. 

A further and very important step to be performed in the near future is to include the contribution of the predictions of the spin  provided by LASSOS in a thermal model of the satellites and compute the changes in their orbit. These changes should than be compared with the corresponding orbital residuals of the satellites obtained within a dedicated POD, in order to see how well these residuals will be reduced by modeling the thermal effects.

The development of a new thermal model for the considered satellites is one of the main targets of LARASE, and we are confident to present soon the model and the results that can be obtained applying the model to a POD.
Nevertheless, we preliminary tested the LASSOS model in this direction with a POD performed for the two LAGEOS with GEODYN in the case of the Yarkovsky thermal drag.
We obtained a slightly reduction in the RMS of the satellites range residuals over a timespan of several years: from 2.8 cm to 2.5 cm in the case of LAGEOS and from 2.5 cm to 2.2 cm for LAGEOS II.
It is important to underline that the thermal model included in Geodyn is a very old model not up-to-date, and it is valid only in the case of the fast spin approximation. Moreover, the starting epoch of the data reduction was that of LAGEOS II launch, consequently LAGEOS was not rotating fast, while LAGEOS II  was in the condition of fast rotation only during the first years of the analysis.

However, this preliminary result represents an argument in favour of the fact that we are on the right path with the LASSOS model.







§ ACKNOWLEDGMENTS

The authors acknowledge the ILRS for providing high quality laser ranging data of the two LAGEOS satellites and of LARES. Special thanks to D. Kucharski for providing the spin measurements of LARES.
We thank our colleagues of the LARASE collaboration and, in particular, R. Peron for useful discussion. This work has been in part supported by the Commissione Scientifica Nazionale II (CSNII) on astroparticle physics experiments of the italian Istituto Nazionale di Fisica Nucleare (INFN). 



§ THE EARTH MAGNETIC FIELD CALCULATED ALONG THE ORBIT OF A SATELLITE

If the Earth magnetic field is approximated by a dipole d_E (with d>0) tilted with respect to the Earth rotational axis in such a way that α_p and θ_p are the longitude and colatitude of the magnetic north pole, the dipole in the J2000 reference frame evolves in time as:


    d_E(t)=d[[ cos(α_p +ω_⊕  t) sin(θ_p); sin(α_p +ω_⊕  t) sin(θ_p);                  cos(θ_p) ]].


The satellite motion is approximated on a circular orbit with Keplerian parameters measured in J2000 frame: M_0, the mean anomaly at t=0; Ω, the longitude of ascending node,  I, the orbit inclination, ω, the argument of pericenter and a, the semimajor axis. Than the satellite motion in the J2000 frame is:



    s=acos(M_0 + n t) [[ cosΩ cosω - cos I  sinΩ sinω; sinΩ cosω + cos I  cosΩ sinω;                  sin I  sinω;                              ]]+ 
    
    - a sin(M_0 +  n t) [
    [ cosΩ sinω + cos I sinΩ cosω; sinΩ sinω - cos I cosΩ cosω;               - sin I  cosω ]].


The magnetic field as function of time over the satellite orbit is therefore given by:

 
	
    B_E=3  s (  d_E· s)/a^5- d_E/a^3= ∑^8_i=0 B_i cos(ω_i t+ φ_i),

	where
	
    ω_0   =   0 
    ω_1   =   ω_2=ω_⊕-2 n  
    ω_3   =   ω_4=ω_⊕+2 n 
    ω_5   =   ω_6=2 n 
    ω_7   =   ω_8=ω_⊕
    φ_i   =   {[ -π/2     ;    0      ].


    B_0   =     - 3  d cos θ_p /4   a^3[ 
    	[  -sin(2I) sinΩ;   sin(2I) cosΩ; -1/3 + cos(2I) ]]
    	
    
     B_1   =   - 3  d sin θ_p (1+cos I )/8   a^3[
    	[ (1-cos I) sin(2ω+2M_0-ϕ_p)+(1+cos I) sin(2Ω+2ω+2M_0-α_p);   (1-cos I)cos(2ω+2M_0-ϕ_p)-(1+cos I)cos(2Ω+2ω+2M_0-α_p);                            -2 sin I ·cos(Ω +2ω+2M_0-α_p) ]] 
    
    
     B_2   =   - 3  d sin θ_p (1+cos I )/8   a^3[
    	[  (1-cos I)cos(2ω+2M_0-ϕ_p)+(1+cos I) cos(2Ω+2ω+2M_0-α_p); -(1-cos I)sin(2ω+2M_0-ϕ_p)+(1+cos I) sin(2Ω+2ω+2M_0-α_p);                            2 sin I ·sin(Ω + 2ω+2M_0-α_p) ]] 
    
    
     B_3   =   - 3  d sin θ_p (1-cos I )/8   a^3[
    	[ -(1+cos I)sin(2ω+2M_0+ϕ_p)+(1-cos I)sin(2Ω-2ω-2M_0-α_p);  (1+cos I)cos(2ω+2M_0+ϕ_p)-(1-cos I)cos(2Ω-2ω-2M_0-α_p);                             2 sin I ·cos(Ω-2ω-2M_0-α_p) ]] 
    
    
     B_4   =   - 3  d sin θ_p (1-cos I )/8   a^3[
    	[ (1+cos I)cos(2ω+2M_0+ϕ_p)+(1-cos I)cos(2Ω-2ω-2M_0-α_p); (1+cos I)sin(2ω+2M_0+ϕ_p)+(1-cos I)sin(2Ω-2ω-2M_0-α_p);                        -2  sin I ·sin(Ω - 2ω-2M_0-α_p) ]] 
    
     B_5   =    - 3  d cos θ_p sin I /2   a^3[
    	[ cosΩcos(2ω+2M_0)- cos I sin(Ω) sin(2ω+2M_0);  sinΩcos(2ω+2M_0)+cos I ·cos(Ω)sin(2ω+2M_0);                          sin I sin(2ω+2M_0) ]]
    
     B_6   =   - 3  d cos θ_p sin I /2   a^3[ 
    	[ cosΩsin(2ω+2M_0)+cos I sin(Ω)cos(2ω+2M_0); sinΩsin(2ω+2M_0)-cos I cos(Ω)cos(2ω+2M_0);                       -sin I cos(2ω+2M_0) ]]
    
     B_7   =    - 3  d sin θ_p /8   a^3[ 
    	[  (1-cos(2I)) sin(2 Ω-ϕ_p)-sin(ϕ_p) (1/3+cos(2I) ); -(1-cos(2I)) cos(2 Ω-ϕ_p)+cos(ϕ_p) (1/3+cos(2I) );                              2 sin(2I) cos(Ω-ϕ_p) ]]
    
     B_8   =   - 3  d sin θ_p /8   a^3[ 
    	[ (1-cos(2I)) cos(2 Ω-ϕ_p)+cos(ϕ_p) (1/3+cos(2I) ); (1-cos(2I)) sin(2 Ω-ϕ_p)+sin(ϕ_p) (1/3+cos(2I) );                             -2 sin(2I)sin(Ω-ϕ_p);                                                  ]].
 

plainnat


