G Pucacco
email: pucacco@roma2.infn.it

M Bassan

Autonomous perturbations of LISA orbits

We investigate autonomous perturbations on the orbits of LISA, namely the effects produced by gravitational fields that can be expressed only in terms of the position, but not of time in the Hill frame. This first step in the study of the LISA orbits has been the subject of recent papers which implement analytical techniques based on a "post-epicyclic" approximation in the Hill frame to find optimal unperturbed orbits. The natural step forward is to analyze the perturbations to purely Keplerian orbits. In the present work a particular emphasis is put on the tidal field of the Earth assumed to be stationary in the Hill frame. Other relevant classes of autonomous perturbations are those given by the corrections to the Solar field responsible for a slow precession and a global stationary field, associated to sources like the interplanetary dust or a local dark matter component. The inclusion of simple linear contributions in the expansion of these fields produces secular solutions that can be compared with the measurements and possibly used to evaluate some morphological property of the perturbing components.

Introduction

Recent papers have revived the analytical study of the orbits of the three spacecrafts of the LISA 'constellation' [1,2,[START_REF] Sweetser | Epicycles and oscillations: the dynamics of the LISA orbits[END_REF]. Although the topic has been explored and analyzed in some detail with the usual numerical methods and many perturbing effects were already considered and dealt with [START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF][START_REF] Povoleri | LISA Orbits[END_REF], we think it can be useful to pursue an analytical approach, that can be of help in gaining insight on the relevance and hierarchy of the various perturbing effects.

The main results obtained so far include the possibility of realizing first-order stable solutions of arbitrary shape in both 'standard' planes at angles near ±60 • with respect to the ecliptic and the optimization of the configuration which produces a peak-to-peak variation in arm length of about 40000 km and a maximum rate of change of arm length of about 4 m/s. These figures represent the least possible amount of 'flexing' of the arms of the nominal equilateral triangular configuration and refer to a purely Keplerian motion in the field of the Sun only, on suitably phase-shifted Earth-like orbits. They can be obtained either by a trial/error procedure based on the direct computation of the exact Keplerian orbits [START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF] or by a perturbative analytical study. It can be shown [2] that, due to the small eccentricity and inclination of the orbits, a series expansion up to second order in a suitable perturbation parameter differs from the exact Keplerian solution by less than 0.03% making the method of analytical series expansion a viable route to study more complex aspects of the motion of LISA.

There are in fact several sources of perturbations affecting the motion of the three spacecrafts on the unperturbed elliptic orbits. It is easy to convince ourselves that all of them are weak enough to produce very small effects over the time of the mission. However, in view of the need of an accurate knowledge of distances and Doppler shifts in the implementation of the Time Delay Interferometry, knowledge, analysis and, possibly, mitigation of all perturbative effects are a necessary task. A preliminary study can be profitably done with those analytical techniques mentioned above, in parallel with more accurate but painstaking numerical simulations [START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF][START_REF] Povoleri | LISA Orbits[END_REF].

The most useful setting in which to perform the analytical study is that based on the Hill system, recalled in sect. [START_REF] Sweetser | Epicycles and oscillations: the dynamics of the LISA orbits[END_REF]. This is a frame rotating with a constant angular velocity set by the mean motion along a reference orbit. The dynamics in this frame are determined by a tidal approximation of the gravitational potential of the Sun and other sources. The far tide approximation [START_REF] Boccaletti | Theory of Orbits[END_REF] consists in taking only the first order term in the field expansion around the origin of the frame. In order to be consistent with the above mentioned series expansion of the Keplerian orbits, it is sufficient to truncate the expansion to the second ('octupolar') order of the field of the Sun. The same procedure can then be applied to the Earth, planets, etc.

These perturbations can be classified as autonomous and non-autonomous. In the first class fall the fields that can be expressed only in terms of the position (and, possibly, of velocity), but not of time in the Hill frame. In this paper we shall focus on the autonomous perturbation effects on LISA orbits. The most relevant example in this class is the field of the Earth, that in the Hill frame is practically at rest. The equations of motion can be cast in the form of a system of oscillators subject to a linear vector potential (due to the Coriolis acceleration) and coupled through higherorder terms. The linear system can be solved exactly, the nonlinear couplings are then included perturbatively. In sect.4 we investigate the configuration, i.e. the phase of each LISA spacecraft orbit, at the closest approach to the Earth, in order to see the effects on all the flexing and Doppler shift indicators, in particular when the time base-line of the mission is longer than usually assumed.

The effects of global stationary fields can also be studied by their series expansion in the Hill frame. The linear part already provides useful information. A first interesting example, that we address in sect.5 is given by the corrections to the gravitational potential of the Sun, treated as a classical spherical mass. The general relativistic correction in the weak field limit, responsible for the classical relativistic perihelion advance, can be obtained by adding a term scaling as the inverse cube of the distance and parametrized by the Schwarzschild radius. Another interesting correction term is provided by the intrinsic quadrupole moment of the Sun. With the current values of the parameters, the effects produced by both these terms are probably within the detection capabilities of LISA. In sect.6 we discuss the tiny effects of the interplanetary dust or a local dark matter component that can in principle be traced in the spectrum of the relative motion between the spacecrafts and can be used to evaluate some morphological property of the perturbing component.

In general, the perturbative effects due to all other bodies in the Solar System have to be expressed by fields explicitly dependent on time and therefore fall in the non-autonomous class. In sect.7 we conclude with some perspective on the study of this second class of perturbations.

Unperturbed Keplerian orbits

We briefly recall the main results, obtained by other authors, that are useful for our work. We begin by writing the unperturbed motion of the three spacecrafts SC k , k = 1, 2, 3, on Keplerian orbits in the field of the Sun alone. In solar barycentric coordinates, let us call X k , Y k , Z k the positions of the SC k along the three orbits. They are given by [1]:

X k = X 1 (ψ k ) cos 2π 3 (k -1) -Y 1 (ψ k ) sin 2π 3 (k -1), Y k = X 1 (ψ k ) sin 2π 3 (k -1) + Y 1 (ψ k ) cos 2π 3 (k -1), (1)
Z k = Z 1 (ψ k ),
where ψ k is the eccentric anomaly of SC k (see below, Eq.6). We then introduce the natural series expansion parameter

α = 2R 1 60 , (2)
where = 5 × 10 6 km is the nominal length of an arm (the side of the LISA triangle) and R = 1.5 × 10 8 km = 1UA is the semi-major axis of the Earth orbit. The authors of ref [1] choose to fix the eccentricity e (the same for all SC's orbits)

e = 1 + 4 3 α 2 + 4 √ 3 α cos θ -1 ( 3 )
by the choice of α and θ, the inclination of the plane of triangle with respect to the ecliptic,

θ = ± π 3 + δ = ± π 3 + αδ 1 , (4)
where δ (alternatively δ 1 ) denotes the additional tilt angle with respect to the canonical value ± π 3 . Finally, the inclination i of the orbital planes with respect to the ecliptic is given by

tan i = 2 √ 3 α sin θ 1 + 2 √ 3 α cos θ . (5)
We can now write the orbit of one 'reference' spacecraft, e.g. SC 1 , in its parametric form:

X 1 = R(cos ψ 1 + e) cos i, Y 1 = R 1 -e 2 sin ψ 1 , Z 1 = R(cos ψ 1 + e) sin i,
The orbits of SC 2 ed SC 3 , are obtained starting from that of SC 1 and performing two rotations around the Z axis, of an angle 2π/3 and 4π/3 respectively, leading to eqs. (1). The eccentric anomaly as a function of time is given by the solution of the Kepler equation. Actually, we now have three anomalies, given as solutions of the three equations

ψ k + e sin ψ k = Ωt - 2π 3 (k -1) . = φ k , k = 1, 2, 3 ( 6 )
Figure 1: Logarithm of the modulus of the Fourier transform of the difference L 12 -L 13 : standard tilt (θ = π/3, δ 1 = 0) (continuous line); δ 1 = 5/8 (dashed line).

where the mean motion is Ω = 2π/T = GM R 3 . If we measure time in years, then T = 1 and Ω = 2π. It is also useful to introduce the phase shifts

σ k = 2π 3 (k -1), (7)
so that φ k = Ωt -σ k . The distances among spacecrafts are clearly given by the three lengths

L jk = (X j -X k ) 2 + (Y j -Y k ) 2 + (Z j -Z k ) 2 , j,k = 1, 2, 3. (8)
These distances are not constant in time, causing what is usually referred to as 'flexing', even for purely Keplerian orbits. A simple way to see this, is to consider that the instantaneous distance from the Sun are different from one spacecraft to the other. The flexing is quite sensitive to the inclination of the plane of the triangle. The value δ 1 = 5/8, giving the angle θ 60 • 35 45 is proved to be optimal by Rajesh Nayak et al. [2]. The maximum flexing of the optimal configuration is about 40000 km and the maximum rate of change of arm lengths is ±4 m/s. The reduction of the relative velocity and related Doppler phase shift with respect to the case considered in [2] (δ 1 = 0) amounts to a factor 5.5. However, it should be noted [START_REF] Bender | [END_REF] that flexing figures comparable to the optimal ones can be obtained simply with slightly different choices of the Keplerian orbital parameters [START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF][START_REF] Sweetser | Epicycles and oscillations: the dynamics of the LISA orbits[END_REF]. The large amount of flexing considered in [2] with δ 1 = 0 is determined by the choice implicit in definitions [START_REF] Sweetser | Epicycles and oscillations: the dynamics of the LISA orbits[END_REF][START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF][START_REF] Povoleri | LISA Orbits[END_REF].

A nice way to appreciate the improvement provided by the optimal choice of the tilt angle is that of comparing the spectral amplitudes of the differential relative motion. Following ref. [11], we can compute the spectra of the difference

∆L kj = L ij -L ik , i,j,k = 1, 2, 3. (9)
In Fig. 1 the two curves represent the modulus of the Fourier transform of ∆L 32 tapered with a Blackman window. The spectrum is computed on a time span of 30 years chosen just to get a good resolution in frequency. The optimal curve has the second harmonic (frequency 2 yr -1 ) seven orders of magnitudes below that of the standard solution with δ 1 = 0. The higher harmonics are much smaller or even unnoticeable. In both solutions the harmonics 3, 6, 9, ... are absent in view of the triangular symmetry of the configuration.

The Hill system

A perturbative analytical approach is convenient when we have to analyze in the simplest way the geometry of orbits and the most relevant perturbations. The Hill reference system (see, e.g. [START_REF] Boccaletti | Theory of Orbits[END_REF][START_REF] Bik | [END_REF], also called as Clohessy-Wiltshire system [14]) allows us to reduce the equations of motion to those of a set of coupled oscillators. The Hill equations properly said consist of the linear part of this system which is readily solved. Non-linear and non-autonomous terms can be added and solved perturbatively. In the specific instance, the system is at rest with the barycenter of the constellation and so an observer sees the Earth at rest in it. The new reference is introduced according to the following steps:

1. Start by considering a circular orbit of radius R around the Sun and pass from the inertial barycentric system to a new system with its origin rotating with the reference orbit and with axes x = x 1 , y = x 2 , z = x 3 . x 1 is directed radially opposed to the center, x 2 is in the direction tangent to the motion and x 3 is perpendicular to their plane. The system is non-inertial and the equations of motion are

ẍ = -∇Φ -2Ω ∧ ẋ -Ω ∧ (Ω ∧ (x + R)) . (10)
Φ is the total gravitational potential and Ω = (0, 0, Ω), R = (R, 0, 0). 2. The potential is approximated according to the expansion

Φ = Φ| x=0 + ∂Φ ∂x λ x=0 x λ + 1 2 ∂ 2 Φ ∂x µ ∂x λ x=0 x µ x λ + O(|x| 3 ), ( 11
)
where the sum is implicit over the repeated indexes λ, µ = 1, 2, 3. The "far tide" (or Hill) approximation [START_REF] Boccaletti | Theory of Orbits[END_REF] consists of retaining only the first order term into the expansion of the gradient of the potential. Taking

Φ = - GM (x + R) 2 + y 2 + z 2 (12)
the Solar far tide approximation gives

∂Φ ∂x x=0 = GM R 2 = Ω 2 R, ∂Φ ∂y x=0 = ∂Φ ∂z x=0 = 0 (13)
and

∂ 2 Φ ∂x 2 x=0 = - 2GM R 3 = -2Ω 2 , ∂ 2 Φ ∂y 2 x=0 = ∂ 2 Φ ∂z 2 x=0 = GM R 3 = Ω 2 , (14)
all others zero. The next order of approximation (octupole) is:

∂ 3 Φ ∂x 3 x=0 = 6GM R 4 = 6Ω 2 R , (15)
∂ 3 Φ ∂x∂y 2 x=0 = ∂ 3 Φ ∂x∂z 2 x=0 = - 3GM R 4 = - 3Ω 2 R , (16)
all others zero. By using these in (10) we get

ẍ -2Ω ẏ -3Ω 2 x + 3Ω 2 2R (2x 2 -y 2 -z 2 ) = 0, ÿ + 2Ωẋ - 3Ω 2 R xy = 0, ( 17
) z + Ω 2 z - 3Ω 2 R xz = 0.
These are the equations for a second order approximation of the purely Keplerian motion. Any other perturbing force can be derived by an analogous series expansion of the corresponding potential and added to this system. Each component, coherently with a perturbative approach, will be a power series in the coordinates with suitably 'small' coefficients ε = (ε x , ε y , ε z ). Recalling the definition of the parameter α in (2) and introducing explicitly only additional autonomous perturbations represented by the perturbation vector f (x, ε), the system (10) at second order can be written in the form

ẍ -2Ω ẏ -3Ω 2 x + 3Ω 2 α (2x 2 -y 2 -z 2 ) -f x (x, ε) = 0, ÿ + 2Ωẋ - 6Ω 2 α xy -f y (x, ε) = 0, ( 18
) z + Ω 2 z - 6Ω 2 α xz -f z (x, ε) = 0.
Since α is small, the nonlinear terms in the multipole expansion of the Solar field can be included in the perturbations. This is an important remark, since in general the presence of quadratic or higher-order terms in the system (18) makes it non-integrable. The most effective approach is that of finding the exact solution of the linear part and using it recursively to find approximate solution for the full non-linear system. The perturbation vector f(x, ε) could even start simply with linear terms of the form εx + . . . giving explicit small changes to the unperturbed zero-order solution. However, it is more natural to consider them of magnitude comparable to terms of order α 2 or higher. These and possibly higher-order perturbing terms propagate in the full solutions through the non-linear coupling with the other components of the field expansion. The solution of the system can therefore be written in the form of a series expansion

x k (t) = x (0) k + αx (1) k + ε x x (2) k + . . . , y k (t) = y (0) k + αy (1) k + ε y y (2) k + . . . , (19)
z k (t) = z (0) k + αz (1) k + ε z z (2) k + . . . ,
where the value of the upper index denotes the order and with the index k = 1, 2, 3 we denote as usual the three spacecrafts. The 'zero' order solution is obtained by putting α = ε = 0 in (18). It amounts to be

x 1 (t) = 2 2x 1 (0) + 1 Ω ẏ1 (0) + 1 Ω ẋ1 (0) sin Ωt -3x 1 (0) + 2 Ω ẏ1 (0) cos Ωt, y 1 (t) = y 1 (0) - 2 Ω ẋ1 (0) (20) -3 (2Ωx 1 (0) + ẏ1 (0)) t + 2 Ω ẋ1 (0) cos Ωt + 2 3x 1 (0) + 2 Ω ẏ1 (0) sin Ωt, z 1 (t) = z 1 (0) cos Ωt + ż1 (0) Ω sin Ωt.
The solutions for the other two spacecrafts are obtained by substituting Ωt with φ k as in definition [START_REF] Bender | [END_REF]. However, the solution actually employed in the perturbation expansion (the epicyclic motion) is obtained by choosing initial conditions in order to remove constant and secular terms [1]:

x (0) k = 2 √ 3 cos φ k , y (0) k = -√ 3 sin φ k , z (0) k = 2 cos φ k . (21)
It corresponds to harmonic oscillations around the 'reference' x = 0 circular orbit. The order 'one' (in α) solution with the integration constants chosen in order to comply with the LISA constraints is [2] x

(1) k = - 5 24 + 2 1 2 -δ 1 cos φ k - 24 cos 2φ k , y (1) k
= - 1 2 -δ 1 sin φ k + 6 sin 2φ k , (22)
z (1) k = √ 3 4 - 2 √ 3 (1 -δ 1 ) cos φ k - 4 √ 3 cos 2φ k .
The solution

x k (t) = x (0) k + αx (1) k (23)
with these explicit forms is a very good approximation of the exact Keplerian orbits, with an error less than 0.03% on the prediction of the flexing of the arms [2,[START_REF] Sweetser | Epicycles and oscillations: the dynamics of the LISA orbits[END_REF]. In other terms, the plots of Fig. 1 can be reproduced using solution (23) with a remarkable accuracy. This result proves that this solution is a reliable starting point for the study of additional sources of gravitational disturbances. From the analytical point of view it is worth mentioning that also system (18) with just f (x, ε) ≡ 0 (of which ( 23) is a first order approximate solution) is not integrable unlike the original Kepler system which is clearly integrable. Is a common lesson of modern analytical mechanics that truncated solutions of non-integrable approximations can be more useful than formally exact solutions of their integrable counterparts.

The Earth effect

As mentioned above the Earth is practically at rest in the Hill system. The Earth coordinates are

x ⊕ = -R(1 -cos 20 • ) -9.0 × 10 6 km, y ⊕ = R sin 20 • 5.13 × 10 7 km, z ⊕ = 0.
We can therefore consider its effect as a constant + linear perturbation to add to the equations of motion. If

d ⊕ = x 2 ⊕ + y 2 ⊕
5.2 × 10 7 km, is the distance from the origin, we can add the perturbative parameter

ε ⊕ = M ⊕ M R d ⊕ 3 7.3 × 10 -5 , (24)
where we used the ratio

M ⊕ M = 1 328900
.

The equations of motion can now be written in the form (18) with the perturbation vector f(x, ε) given by

f ⊕ = -ε ⊕ Ω 2 (x -x ⊕ ). ( 25
)
Proceeding as assumed in eqs. [START_REF] Weinberg | Gravitation and Cosmology[END_REF], we may look for a solution that can be linearly superposed to (23) so that we have

x k (t) = x (0) k + αx (1) k + ε ⊕ x (2)
k .

(

) 26
The explicit solution is [START_REF] Dhurandhar | [END_REF]:

x (2) k = 2 A k Ω + x ⊕ + 2y ⊕ φ k + B k cos φ k + C k sin φ k - 5 4 √ 3 φ k sin φ k , y (2) k
= -(3 A k Ω + 2x ⊕ )φ k - 3 2 y ⊕ φ 2 k + ( 2 7
)

- 5 2
√ 3 φ k cos φ k + √ 3 2 sin φ k + 2(C k cos φ k -B k sin φ k ) + D k , z (2) k
= E k cos φ k + F k sin φ k - 4 φ k sin φ k ,
where the A k , ..., F k are constants determining the initial conditions. Recalling that the mean anomalies φ k are just phase shifted rescaling of time, in these solutions we see the presence of additional secular terms. In particular, the y coordinate has a term quadratic in time that cannot be eliminated by any choice of the initial conditions. Dhurandhar et al. [START_REF] Dhurandhar | [END_REF] propose to inject LISA on orbits with initial conditions devised to reduce to a minimum the Earth effect. In view of the secular terms, it seems a good strategy that of minimizing the disturbance at 'mid mission': they assume a mission duration of three years, so that this 'optimal time' is at 1.5 years after starting the operation activity. The authors of ref. [START_REF] Dhurandhar | [END_REF] introduce in the mean anomalies φ k a phase shift to accommodate for different starting times and compute the resulting effect on the chosen period. However, the solution is invariant under time translation in view of the autonomous nature of the system, therefore, expressions (27) can be profitably used at arbitrary initial time. If this, for simplicity, is put equal to zero, the values of the constants appearing in (27) are

A k = Ωy ⊕ σ k + Ω √ 3 cos σ k , B k = - 13 8 √ 3 -x ⊕ cos σ k -2y ⊕ sin σ k - √ 3 8 cos 2σ k , C k = x ⊕ sin σ k -2y ⊕ cos σ k + 3 sin 2σ k -10σ k 8 √ 3 , (28)
D k = y ⊕ 4 - 3 2 σ 2 k -2x ⊕ σ k - √ 3 σ k cos σ k + 4 √ 3 sin σ k , E k = - 4 sin 2 σ k , F k = - 8 (2σ k -sin 2σ k ) ,
where the σ k are the phase shifts introduced in ( 7). Experience gathered on many space experiments of the past allows some optimism on expecting the LISA mission to last longer than the nominal duration of three years. A longer lifetime also allows us to better investigate the secular (non periodic) behavior of the orbits. In order to clarify the effect of the Earth tidal field, we have therefore chosen to consider an "optimistic duration" and to double the expected mission time: in Fig. 2 the three arm lengths time are plotted over a period of 6 years with the solutions (26) and the initial conditions above. In Fig. 3 they are plotted, on the same six years period, the relative velocities v jk = Ljk . Now we can see that, if we limit the observation to the three years centered around zero, both the flexing and the relative velocity have an almost constant envelope. However, in the long run the envelope tends to linearly grow for both quantities and is therefore appropriate a comparison with the 'optimal' solution illustrated in section 2. In a three-year (±1.5) mission we have a flexing of 60000 km which increases to 80000 km in six (±3) years: on this longer time span we arrive at a factor 2 of increase of the flexing with respect to the purely Keplerian solution. The maximum rate of change of arm lengths is ±5.5 m/s in the ±1.5 years time span but increases to ±10 m/s in the ±3 years time span, a factor 2.5 of increase with respect to the Keplerian solution.

This behavior is due to the global dynamics in the coupled fields of Sun and Earth as described in the rotating frame. As can be seen from the solution (27) with the constants determined by (28), the secular terms always produce a drift away from Earth: this is linear in x and, as remarked above, quadratic in y. The secular drift in x of the three SC k amounts to ±70000 km on the ±1.5 years time span, whereas the motion in y is of initial approach, a state of rest and a subsequent departure: on the same time span, the drift in y is of 5 × 10 5 km. So, the best that we can do is a choice of initial conditions such that LISA is a little further out at start, approaches the Earth, is at rest ad mid mission and departs after that. The physical reason of this behavior is evident when we consider the nature of the effective potential. LISA is well outside the Roche lobe associated to the Earth; therefore it cannot 'fall' (as a stone or a projectile [START_REF] Dhurandhar | [END_REF]) towards it. Rather, it moves on a slow merry-go-round and feels a small Coriolis apparent force. The combined parabolic-drift motions have therefore quite a relevant role in increasing the flexing on a certain time base.

However, the possibility exists of an appreciable reduction of the flexing due to the Earth tidal field over a limited time span by a suitable shift in the relative phase between the epicyclic Figure 5: Spectrum of the difference ∆L 32 = L 12 -L 13 for the case γ = 0. The dashed line is the Keplerian configuration and the continuous line represents the addition of the Earth perturbation. motion and the linear perturbation. The scenario depicted above can in fact be changed if the constellation gets to the point of minimum distance from the Earth with a different orientation. This situation can be described either by modifying the phase shift in the initial conditions of the linear perturbation or by putting this extra shift directly in the definition [START_REF] Dhurandhar | [END_REF]. Let us modify this in

σk = 2π 3 (k -1) -γ k , k = 1, 2, 3 ( 2 9 )
using the angles γ k to change the orientation of the triangle at a given time. We limit the analysis to identical angles for all the SC k , however a non-rigid rotation could also be studied. The standard zero order solution used above, with γ = 0, of (21) corresponds to a configuration with the spacecraft 1 at the 'highest' point (z 1 = /2) over the ecliptic at time zero (let us say at mid mission). The spacecraft 1 being at z 1 = 0 with either y 1 = + / √ 3 or y 1 = -/ √ 3 correspond to the two extremes in which it is closest or farthest from the Earth at time zero. It is interesting to examine these two situations compared with that above: they are given by putting γ = ±π/2 in (29). In Fig. 4 the first case (γ = -π/2) is illustrated by plotting the three arm lengths time over the period of 6 years: in this case in which SC 1 is the one closest to the Earth at time zero, the two arms 12 and 13 exhibit a variation not larger than the optimal one (40000 km) for the first half-mission period of 3 years. It can be shown that in the second case (SC 1 farthest from the Earth at time zero) the plot is the time reversal of that in Fig. 4, so that the variation of L 12 , L 13 is kept to a minimum in the second 3 years. The price for this reduction is that the flexing of the third (23) arm is higher, but only slightly with respect to the level found with γ = 0 (80000 km) on the same time span.

Another way to examine the effect of the Earth is by comparing the differential relative motion computing the difference ∆L 32 = L 12 -L 13 and the corresponding spectral amplitudes. In Fig. 5 the modulus of the Fourier transform of ∆L 32 is shown together with the optimal Keplerian case of Fig. 1 (no Earth contribution). The continuous line includes the contribution of the Earth in the case of γ = 0. As before, to get a good resolution in frequency, the spectrum is computed on a longer time span (30 years). The effect of the Earth is evident from the appearance of the harmonics 3, 6, 9, ... due to the fact that the triangular symmetry of the configuration is now broken. The configurations with the phases γ = ±π/2 display the same spectral lines with slightly lower values at higher frequencies.

The analysis of these and other possible configurations could be of a certain practical relevance should it be necessary to privilege the proper working of only a pair of arms [16].

The relativistic precession and the effect of the Solar quadrupole

In the list of small perturbative effects given by fixed fields naturally fall those connected with the not exactly Keplerian interaction with the Sun. Two important examples are the relativistic correction of the classical field and the perturbation due to a possible quadrupole moment of the Solar gravitational field. The first produces the slow precession responsible for the advance of the perihelion; the second produces a more complex variation of orbital elements but of (supposedly) much smaller amount.

It is tempting to see if LISA is capable to detect small effects due to the violation of the spherical symmetry of the Solar field, such as we have assumed so far. The small eccentricity and the distance from the Sun of Earth-like orbits makes the amount of the effects quite small, but worth investigating in view of the high sensitivity of LISA to small differential variations of the arm lengths. On this ground there is the recent proposal by Polnarev et al. [15] of detecting with LISA the low frequency g-mode Solar oscillations.

The relativistic correction can be mimicked by the presence of a small additional term in the potential of the form Φ S = -r S L 2 r 3 , where r S = GM c 2 is the Schwarzschild radius of the Sun and L is the angular momentum of the orbit related to the other parameters by L 2 = GM R(1 -e 2 ). The Solar quadrupole is associated with a small term of the form

Φ Q = GM J 2 R 2 r 3 • 2Z 2 -X 2 -Y 2 2r 2 , ( 30
)
where J 2 = Q × 10 -7 characterizes the mass quadrupole of the Sun, Q being a parameter of order 1 according to current estimates [18]. R 7 × 10 5 km is the radius of the Sun and the radial distance is expressed in terms of the coordinates as

r = X 2 + Y 2 + Z 2 = (R + x) 2 + y 2 + z 2 .
An approximate evaluation of the two effects can be attempted in the Hill frame with the usual hypothesis of small perturbations of almost circular orbits. The potentials with the perturbing terms introduced above are expanded in series about the origin of the Hill frame and the first and second order terms are retained and added to the Keplerian terms. We remark that a small constant term appears in the x component because the centrifugal force in the rotating frame exactly balances only the Keplerian monopole at r = R: the extra term is responsible for the classical precession effect common to all spacecrafts. The linear terms are given by the different quadrupole terms arising in the expansions: these give significant contributions in the form of differential effects. The equations of motion are of the same form for both effects and can be written in the form (18) with the perturbation vector f(x, ε) given by

f x = -ε x Ω 2 (R -4x) , f y = -ε y Ω 2 y, f z = -ε z Ω 2 z. (31)
The Schwarzschild correction is obtained if the coupling constants are given by

ε x = ε y = ε z = ε S ≡ 3 r S R 3 × 10 -8 (32)
whereas the Solar quadrupole is included when the coupling constants are given by

ε x = ε y = 1 3 ε z = ε Q ≡ 3 2 J 2 R R 2 3.3 Q × 10 -12 . (33)
Figure 6: Difference ∆L pert jk -∆L kepl jk time for a period of 5 years as a consequence of the differential effect due to the relativistic precession: the continuous line is for jk = 23, the long dashed line for jk = 13 and the short dashed line for jk = 12.

Figure 7: Difference ∆L pert 32 -∆L kepl 32 in units of Q time for a period of 5 years as a consequence of the differential effect due to the Solar quadrupole: the continuous line is for jk = 23, the long dashed line for jk = 13 and the short dashed line for jk = 12.

The structure of the system is the same as in the case of the Earth (cfr. Eqs. (18) with the perturbation vector given by ( 25)) and the strength of the perturbing terms can be evaluated by comparing the above constants to ε ⊕ 7.3 × 10 -5 . For the consistency of the approach we preliminarily check the prediction concerning the overall precession. The constant force appearing in the x component is a factor R/ = (2α) -1 = 30 stronger than that linearly depending on the coordinates. Therefore we can assign to the solution the following structure

x k (t) = x (0) k + αx (1) k + ε x 2α x (2) + x (3) k (ε). ( 34
)
In the second-order term there is no index k because the solution satisfying proper initial conditions on the epicycle is the same for each spacecraft. The structure of this term has to comply with the natural condition that the radial oscillation be symmetric with respect to the circular reference orbit whose radius is slightly reduced by the constant component of the force. The choice of initial conditions

x k (0) = 2 √ 3 cos σ k , y k (0) = √ 3 sin σ k , z k (0) = 2 cos σ k ; ẋk (0) = Ω 2 √ 3 sin σ k , ẏk (0) = - Ω √ 3 cos σ k - ε x 2 √ 3α , żk (0) = Ω 2 sin σ k ,
gives the correct solution

x (2) (t) = - 3 (1 -cos Ωt) , y (2) (t) = Ωt - 2 3 sin Ωt , z (2) (t) = 0, (35)
since the constant part in x(t), -ε x (6α) -1 = -ε x R/3, accounts for the change of radius of the reference orbit. An estimate of the advance of the perihelion after a revolution is given by ∆ϕ ε y 2α

y (2) ( 2π Ω ) R = 2πε y .
Using Eqs.( 32) and (33) we get respectively

∆ϕ S 6πr S R [rad/orbit]
for the relativistic advance of the perihelion and

∆ϕ Q 3πJ 2 R R 2 [rad/orbit]
X -Y plane. Another apparently very crude assumption is that of assuming a linear gradient field as a second order perturbation in our system of equations of motion. This hypothesis is however viable since almost circular orbits are probing the extra field in a small radial range in which it is slowly changing and what can be of practical interest is the possibility of detecting an anisotropy in the field, possibly due to the ellipticity of the spheroidal distribution. The hypothesis is sustained by numerical computations with the interplanetary dust [11] showing that its perturbing effect is weakly affected by the steepness of the radial component of the density function.

The system of equations that we use is of the form (18) with the perturbation vector f(x, ε) given by

f x = -ε D Ω 2 x, f y = -ε Ω 2 y, f z = -βε D Ω 2 z. ( 38
)
Here ε D is the small parameter associated to the diffuse component and we have assumed a spheroidal distribution parametrized by a shape parameter β. In the very simple setting sketched above, we may assume

ε D = 4πGρ D 3Ω 2 , ( 39
)
where ρ D is a suitable average over the density distribution. The solution is still of the form (26) where in the second order terms we add the contribution due to the solution

x (2) k D = 2 A k Ω + B k cos φ k + C k sin φ k - 5 4 √ 3 φ k sin φ k , y (2)
k D = -3 A k Ω φ k + 2(C k cos φ k -B k sin φ k ) + D k - 5 2 √ 3 φ k cos φ k - 3 5 sin φ k , (40) z (2)
k D = E k cos φ k + F k sin φ k - 4 φ k sin φ k
and the arbitrary constants have a prime to distinguish them from those of the solution for the Earth. With the natural condition of putting to zero both the displacement and its time rate at 'zero' (half-mission) time, we get

A k = Ω √ 3 cos σ k , B k = - 13 8 √ 3 - √ 3 8 cos 2σ k , C k = 3 sin 2σ k -10σ k 8 √ 3 , D k = - √ 3 σ k cos σ k + 4 √ 3 sin σ k , (41)
E k = - 4 sin 2 σ k , F k = - 8 (2σ k -sin 2σ k ) .
The components in (40), multiplied by ε D and added to the solution (23), provide the evolution of the arm-lengths under this perturbing effect. To show the small effects of the diffuse component, Cerdonio et al. [11] (whose approach is purely numerical) plot the following quantities

δL jk = ∆L D jk -L kepl jk , (42)
where the arm-lengths L jk D are computed using the full perturbed solution and the L jk are those obtained with the optimal Keplerian solution. In Fig. 9 these relative distances are plotted in the case in which the average density is taken to be ρ D = 9.6 × 10 -20 kg /m 3 as estimated by models of the meteoritic complex [17] and a shape parameter β = √ 3 (spheroidal distribution): we can compare these plots with those for the α = 0 case of reference [11], corresponding to their most homogeneous distribution. Our solution is in good agreement with that obtained with a numerical integration of the Gauss equations [11]. In Fig. 10 we have plotted the Fourier spectrum of two different realizations of the difference δL 32 of Eq.(42). In ref. [11] an analogous plot has been obtained by the discrete Fourier analysis of the numerical solution and we can again appreciate the good accuracy of this simple analytical approach. We can see how the solutions suggest the possibility of assessing the shape of the matter distribution through the spectral analysis of the arm lengths variation. However, we have seen in the previous sections that other autonomous perturbations effects leave on these observables marks at the same frequencies but with higher amplitude.

Conclusions and future work

We have reported on a systematic study of the autonomous perturbations of LISA orbits, namely those stationary with respect to the center of the LISA constellation, i.e. constant in the Hill frame. This is the simplest step in the study of the perturbations on the LISA orbits.

As in other recent papers [1,2,[START_REF] Sweetser | Epicycles and oscillations: the dynamics of the LISA orbits[END_REF], we implement analytical techniques based on a 'post-epicyclic' approximation in the Hill frame. This work is intended as a preliminary study aimed at gathering an overall picture and more physical insight to the issue. In fact, the analysis of changes of the ideal Keplerian orbits may suggest in principle possible changes in the choice of initial conditions. Accurate, dedicated numerical investigations [START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF][START_REF] Povoleri | LISA Orbits[END_REF][START_REF] Bik | [END_REF] are then required in the assessment of the exact amount of each effect. However, a numerical optimization depends, in the case of independent orbits, on eighteen parameters and is yet to be fully explored [START_REF] Hughes | Preliminary optimal orbit design for LISA[END_REF]. A systematic comparison of our analytical approach with these solutions is a further step that must be undertaken.

A particular emphasis was put on the tidal field of the Earth, assumed to be stationary in the Hill frame. The solution of the equations of motion under the perturbation gives a combination of a reversed parabolic motion in the tangential direction and a drift in the radial direction. This sheds light on possible implications on the choice of initial conditions when the time base-line of the mission is longer than that assumed in previous papers. We propose a more systematic investigation of the phase choice affecting the orientation of the constellation at the moment of closest approach to the Earth that could be of a certain practical relevance in the case it should be necessary to privilege the proper working of only a pair of arms.

Another relevant class of autonomous perturbations is that given by a global stationary field, associated either to the classical or relativistic corrections to the spherically symmetric Solar potential or to sources like the interplanetary dust or a local dark matter component. The inclusion of a simple linear contribution representative of these fields produces a secular solution that can be used to evaluate physical parameters of the perturbing force. The relativistic precession and the more complex motion due to a Solar quadrupole contribution are most probably within the detection capabilities of LISA. More problematic, with the figures obtained, is a detection of a possible flatness of the perturbing spheroidal component.

The natural step further on this research line is to attempt the analysis of non-autonomous perturbations. The time-varying contributions due to the Moon, the planets and other celestial bodies must be included in the list of relevant perturbing effects. On the basis of the present work, we are confident that the study of a properly modified version of the equations of motion in the form of a forced system of coupled nonlinear oscillators (e.g. that recently presented in ref. [12]) will provide additional solutions with a degree of accuracy comparable to that presented here.

Figure 2 :

 2 Figure 2: Distances L ij vs time including the Earth influence for a period of 6 years: L 12 , continuous line; L 23 , long dashes; L 13 short dashes.

Figure 3 :

 3 Figure 3: Rates of change v ij = Lij vs time including the Earth influence for a period of 6 years: v 12 , continuous line; v 23 , long dashes; v 13 short dashes.

Figure 4 :

 4 Figure 4: Distances L ij time including the Earth influence for a period of 6 years with an extra γ = -π/2 phase: L 12 , continuous line; L 23 , long dashes; L 13 short dashes.

Figure 9 :√ 3 )

 93 Figure 9: Relative distances δL jk time due to a spheroidal (β = √3) distribution of interplanetary dust, over a time-span of 5 years: the continuous line is for jk = 23, the long dashed line for jk = 13 and the short dashed line for jk = 12.

Figure 10 :

 10 Figure 10: Spectrum of the relative distance δL jk for two values of the shape parameter β: β = √ 3 (continuous line) and β = 1 (dashed line).

m years

Acknowledgments

We thank Pete Bender, Alfio Bonanno, Roberto De Pietri, Oliver Jennrich, David Lucchesi, Ruggero Stanga and Michele Vallisneri for fruitful discussions. GP acknowledges partial support from ICRAnet.

for the corresponding advance due to the J 2 term. These results are in agreement with standard findings in the limit of small eccentricities [START_REF] Weinberg | Gravitation and Cosmology[END_REF] making us confident about the reliability of the solution. The higher order terms are

To better show the effects produced by these perturbations we can use the quantities

where ∆L kj = L ij -L ik and the arm-lengths L jk are computed using the full perturbed solution (34) in ∆L pert jk and including only the Keplerian solutions in ∆L kepl jk . In Fig. 6 we show these quantities time for a period of 5 years as a consequence of the differential effect due to the relativistic precession: the amplitude of the arm length variations can be as large as 10 km. In Fig. 7 we show the same quantities due to the differential effect of the Solar quadrupole: the amplitude of the arm length variations reaches 1.5 m times Q. The modulus of the Fourier transform of ∆L pert 32 -∆L kepl 32 is shown in Fig. 8. In principle, a precise evaluation of all other effects (Schwarzschild precession included) would allow a measurement of Q (and J 2 ) from the residuals. Indeed, Fig. 8 shows that the relative content of higher harmonics with respect to the fundamental is larger for the Solar quadrupole precession than for the relativistic one. However, this measurement is hindered by the absolute magnitude of the first effect being orders of magnitude smaller than the latter.

A global stationary potential

The presence of a stationary global potential in the Solar System can be due to different kinds of diffuse matter distributions. Main candidates are the interplanetary dust and a local dark matter component [10,11]. Our concern here is to show that this kind of gravitational source is readily described as an autonomous perturbation and can be analyzed with the methods illustrated above. We consider a very simple model in which a spheroidal distribution can be assumed with rotational symmetry around the (barycentric) Z-axis and mirror reflection symmetry with respect to the