Trajectory Design for the ESA LISA Mission

Abstract

The three Laser Interferometer Space Antenna (LISA) spacecraft are going to be placed in a triangular formation in an Earth-trailing or Earth-leading orbit. They will be launched together on a single rocket and transferred to that science orbit using Solar Electric Propulsion. Since the transfer Δv depends on the chosen science orbit, both transfer and science orbit have been optimised together. For a thrust level of 90 mN, an allocation of 1092 m/s per spacecraft is sufficient for an all-year launch in 2034. For every launch month a dedicated science orbit is designed with a corner angle variation of 60^° ± 1.0^° and an arm length rate of maximum 10 m/s. Moreover, a detailed navigation analysis of the science orbit insertion and the impact on insertion errors on the constellation stability has been conducted. The analysis shows that Range/Doppler measurements together with a series of correction manoeuvres at the beginning of the science orbit phase can reduce insertion dispersions to a level where corner angle variations remain at about 60^° ± 1.1^° at 99% C.L. However, the situation can become significantly worse if the self-gravity accelerations acting during the science orbit phase are not sufficiently characterised prior to science orbit insertion.

Introduction

The LISA mission has been selected as the L3 cornerstone mission by the European Space Agency (ESA) in June 2017 [1, 2]. It will be the first space-based Gravitational Wave detector and will explore a new frequency band and therefore be sensitive to previously undiscovered sources of Gravitational Waves [3]. LISA will be based on laser interferometry between free-flying test masses inside drag-free spacecraft. In order to realise that concept, the three LISA spacecraft will be placed into a heliocentric Earth-trailing or Earth-leading orbit. Together they will form the corners of a nearly equilateral triangular formation with 2.5 × 10⁶ km arm length [4]: the so-called cartwheel formation. Each spacecraft is forced to follow its two test masses along each of the two laser beam axes they define. The orbits have to be selected such that the stability of the constellation for the planned mission lifetime of 10 years (4y nominal + 6y extension) is guaranteed. The main preliminary stability requirements are currently.

• 1. The corner angle variations during 10 years shall not exceed 60^° ± 1.0^°. This constraint comes mainly from the limitation of the optical assembly tracking mechanism (OATM) which ensures that the laser beams point into the correct direction at all times.

• 2. The relative velocity (arm length rate) shall not exceed 10 m/s during 10 years mission duration. This is due to the bandwidth limitation of the interferometric beat note detection of the system [5].

• 3. The arm length shall be constrained within 2.5 × 10⁶ km ± 2.5 × 10⁵ km during 10 years. In practice this constraint is never limiting, however.

These stability requirements can be achieved by maximising the distance of the formation to Earth, and thus minimising its gravitational impact on the formation. On the other hand, communication requirements place a limit on the maximum Earth distance, which is currently assumed to be 65 × 10⁶ km. This value, however, depends on the design of the spacecraft communications system.

It is worth noting that the stated requirement values on corner angle variations and arm length rate are preliminary and are subject to change as the LISA study progresses. The approach taken in this paper is to use the values of 60^° ± 1.0^° and 10 m/s for the nominal mission design. Insertion dispersions and other random perturbations, however, may violate these requirements. The next iteration on system level may either relax the requirements on corner angle variations and arm length rate to fit the existing trajectory design including dispersions, or the trajectory has to be re-designed to accommodate the dispersions in the current requirement.

LISA is going to be launched on a single Ariane 64 rocket and transferred to its operational orbit using Solar Electric Propulsion (SEP). SEP has the advantage of a more efficient use of propellant mass and thus a higher dry mass fraction compared to chemical propulsion. Previous trajectory analyses for LISA [6,7,8,9,10,11,12] have mainly focussed on the optimal cartwheel orbit design and less on the transfer. Moreover, the interdependency between the transfer and the optimal cartwheel orbit has not been treated before. The present paper is going to discuss the cartwheel orbit design at different levels of accuracy: Starting from a review of analytical models to finally presenting fully numerical results for both optimal transfer and cartwheel orbits.

The question of stability of the cartwheel orbit under insertion uncertainties has not been analysed before to the authors’ knowledge. This is going to be addressed in depth in the second part of the present paper.

The document is structured as follows: after a summary of the most important analytical cartwheel models and results in The Cartwheel Formation - Analytic Models, a fully numerical simulation of the optimal cartwheel orbit is presented in The Cartwheel Formation - Fully Numerical Model. Solar Electric Propulsion Transfer will treat the optimisation of the SEP transfer and the interdependency between the transfer and the cartwheel orbit design. A full navigation analysis taking into account all major sources of uncertainty will be presented in Navigation and Insertion Accuracy. Also, the impact of the computed cartwheel insertion accuracy on the stability of the orbit during 10 years of science phase will be addressed there. Finally, Conclusions will summarize the main conclusions.

The Cartwheel Formation - Analytic Models

In the following a series of analytic models of the cartwheel formation with increasing level of fidelity will be discussed.

Linear Two-Body Model

In the two-body problem, the linearised relative motion around a centre on a circular orbit is described by the Clohessy-Wiltshire equations [13]. Imposing no along-track drift, these equations have the solution [14, 15]:

$$ x(t)={\varrho}_x\sin\ \left(\frac{2\pi }{T}t+{\alpha}_x\right), $$

(1)

$$ y(t)={\varrho}_y+2{\varrho}_x\cos\ \left(\frac{2\pi }{T}t+{\alpha}_x\right), $$

(2)

$$ z(t)={\varrho}_z\sin\ \left(\frac{2\pi }{T}t+{\alpha}_z\right), $$

(3)

where x(t), y(t), z(t) are the components of the local orbital frame of the virtual formation centre as a function of time:

• x= radial (oriented positively from central body to spacecraft)

• z= cross-track (along spacecraft angular momentum vector)

• y= completing the right-handed frame (along-track for a circular spacecraft orbit)

This is illustrated in Fig. 1. T is the orbital period of the circular formation centre orbit and α_x, α_z, ϱ_x and ϱ_z are integration constants. The relative trajectory described by eqs. 1-3 is an ellipse in the x − y plane and a closed Lissajous figure in 3D since the phases α_x and α_z are not equal in general. By additionally imposing α_x = α_z ≡ α₀ and ${\varrho}_z=\pm \sqrt{3}{\varrho}_x$ the relative trajectory becomes a circle around the formation centre which is inclined by 60^° w.r.t. the orbital plane of the formation centre. Obviously, there are two solutions depending on the sign of ϱ_z. These will be called the clockwise (+) and counter-clockwise (−) solutions in this document describing the two possible orientations of the cartwheel triangle.¹

Fig. 1

Illustration of the local orbital frame

By placing three spacecraft on this circle with a relative phase of 120^°, an equilateral triangle formation with arm length $d=2\sqrt{3}{\varrho}_x$ is achieved. Using the above equations, the heliocentric orbital elements of the three spacecraftcan be derived. They are summarized in Table 1.

Table 1 Heliocentric, ecliptic orbital elements for LISA in the linear model

The reference frame of these orbital elements is a heliocentric, ecliptic frame, frozen at the initial time, t = 0, and with the axes defined as follows:

• x= from Sun to formation centre (clockwise solution) or from formation centre to Sun (counter-clockwise solution)

• z= along the angular momentum vector of the formation centre

• y= completing the right-handed frame

Note that since the axes are frozen at t = 0, it is an inertial frame.

In order to fix the relative position of the LISA formation with respect to Earth, the Mean Initial Displacement Angle (MIDA) needs to be chosen as well (cf. section 2.3 for details) as a free parameter to fully define the system. For a model with a circular Earth orbit, the MIDA is identical to the instantaneous initial displacement angle. The difference in the case of an elliptic Earth orbit will be clarified in section 2.3. The choice of MIDA defines the rotation of the reference frame above with respect to the standard ecliptic reference frame at t = 0. Overall, there are four free continuous and one discrete free parameter in the system:

• 1. arm length, d

• 2. semi-major axis, a

• 3. right ascension of the ascending node (RAAN) of LISA1, Ω₁

• 4. argument of perihelion (discrete) ω

• 5. MIDA, θ₀

For LISA the semi-major axis would, however, be chosen to be equal to the one of Earth (at least in the linear two-body model) in order to avoid a relative drift.

Keplerian Two-Body Model

The model described in section 2.1 can be improved by taking into account the non-linear effects of the Keplerian motion. Still no Earth perturbations are considered here. This is described in reference [8]. The cited reference also shows that in this approximation the flexing of LISA’s arms can be minimized by deviating slightly from the 60^° triangle inclination. This deviation is parametrised by the angle, δ. However, it is convenient to instead use the dimension less parameter δ₁ defined by δ = αδ_1, with α = d/2a. The optimum value is shown to be δ₁ = 5/8 (this is only true in the Keplerian model). The orbital parameters can be computed from the formulas given in Table 2 where the eccentric anomaly, E_k and the mean anomaly, M_k are related by Kepler’s equation:

Table 2 Heliocentric, ecliptic orbital elements for LISA in the Keplerian model

$$ {E}_k+e\sin\ {E}_k={M}_k=\pi -{\sigma}_k(t),\kern1em k=1,2,3. $$

(4)

And the three clocking angles are defined as:

$$ {\sigma}_k(t)={\sigma}_0+\left(k-1\right)\frac{2\pi }{3}-\sqrt{\frac{\mu_{\mathrm{Sun}}}{a^3}}t,\kern1em k=1,2,3, $$

(5)

where σ₀ is called the initial clocking angle.

The reference frame of the orbital elements in Table 2 is the same as for Table 1 and defined in section 2.1. Overall, there are now five free continuous and one discrete free parameter in the system:

• 1. arm length, d

• 2. semi-major axis, a

• 3. inclination parameter, δ₁

• 4. initial clocking angle, σ₀ = σ₁(t = 0)

• 5. argument of perihelion (discrete) ω

• 6. MIDA, θ₀

These (except ω) are illustrated in Fig. 2.

Fig. 2

Illustration of cartwheel parameters in the Keplerian model

While the corner angles in the case of the linear model described in section 2.1 are constant at 60^°, the non-linear effects of the Keplerian model introduce a “breathing” of the triangle. The resulting evolution of the corner angles, arm length rate and arm length is shown in Fig. 3 for a = 1 AU, d = 2.5 · 10⁶ km and δ₁ = 5/8.

Fig. 3

Evolution of cartwheel corner angles (top), arm length rate (middle) and arm length (bottom) for the Keplerian cartwheel model for LISA 1 (blue), LISA 2 (red) and LISA 3 (green)

Although this model is lacking the effect of the Earth gravity perturbation, which is significant for MIDAs close to 20^°, it is still useful as an initial guess for an optimisation using a full numerical model. Moreover, the Keplerian model can be used for the arrival state during the transfer optimisation. This allows for an easy exploration of a large range of orbit options without previous computationally heavy optimisation of the science orbit.

Mean Initial Displacement Angle (MIDA)

The MIDA, θ₀ defines the location of the centre of the cartwheel with respect to the Earth (see Fig. 2). A negative MIDA denotes a trailing configuration, a positive MIDA a leading configuration. The larger the MIDA magnitude, the lesser the gravitational perturbations of the cartwheel, but the larger the distance to Earth and thus the impact on the communications subsystem. Conversely, the smaller the MIDA, the lower the Δv cost of transfer. A larger MIDA generally also means a higher transfer Δv, due to the larger difference in semi-major axis of the transfer orbit, if the transfer duration is kept fixed.

During the science phase the centre of the LISA formation centre moves on a (very close to) circular orbit around the Sun, while the Earth orbit is eccentric. This leads to a natural variation of the instantaneous Earth displacement angle of about ±2^° over one year. This fact makes the instantaneous displacement angle an inconvenient design parameter when different initial times are compared. In order to compare different launch dates over one year, it is therefore more convenient to use the MIDA, θ₀. The angle averages out the Earth eccentricity by measuring the displacement angle w.r.t. the “Mean Earth”, a virtual body with the same orbital period as the Earth, but a circular orbit around the Sun (see also ref. [12]). The orbital elements of the Mean Earth are summarised in Table 3.

Table 3 Orbital elements of the Mean Earth from True Earth

Earth Perturbations and Choice of Initial Semi-Major Axis

The Keplerian model and cartwheel state definition as described in section 2.2 is suitable for use in transfer optimization and as initial guess for the science orbit optimization. The main perturbation that significantly impacts the LISA science orbit evolution is the Earth’s gravity. It is not possible to analytically solve the model when the Earth’s gravity is included. However, there are a few useful analytical formulas that can be derived to aid the orbit design. These are described in this section.

The main effect of the Earth’s gravity (besides perturbing the triangular formation) is a (near) along-track acceleration on the LISA spacecraft which causes a drift in the LISA semi-major axis. The relative difference in semi-major axis to the one of the Earth orbit leads to a drift of the Earth distance. This will be analysed in the following. The rate of change of the semi-major axis is described by Gauss’ variational eqs. [14]:

$$ \frac{da}{dt}=2\frac{d_R{a}^2e\sin \nu }{\sqrt{\mu_{\mathrm{Sun}}a\left(1-{e}^2\right)}}+2\frac{d_T{a}^3\left(1-{e}^2\right)}{r\sqrt{\mu_{\mathrm{Sun}}a\left(1-{e}^2\right)}}. $$

(6)

Where d_R and d_T are the radial and transversal components of the acceleration vector. Circular orbits are assumed here. θ is assumed small enough such that the radial component can be neglected. The transversal/along-track acceleration from the Earth’s gravity perturbation is given by:

$$ {\displaystyle \begin{array}{c}{d}_T=\frac{\mu_{\mathrm{Earth}}}{r_{\mathrm{Earth}}^2}\cos\ \left(\theta /2\right)\\ {}=\frac{\mu_{\mathrm{Earth}}\cos \left(\theta /2\right)}{4{a}^2{\sin}^2\left(\theta /2\right)}\end{array}} $$

(7)

Where r_Earth is the Earth distance from the constellation centre. The drift rate of the mean semi-major axis from eq. 6 thus becomes:

$$ \frac{da}{dt}=\pm \frac{\mu_{\mathrm{Earth}}\cos \left(\theta /2\right)}{2\sqrt{a{\mu}_{\mathrm{Sun}}}}{\sin}^2\left(\theta /2\right). $$

(8)

The + sign shall be used for the trailing configuration and the – sign for the leading configuration. Evaluating the equation at a = 1AU and mean Earth displacement angle θ = θ₀ gives sufficiently accurate results if θ is not too far from 20^° as a-posteriori verification shows (cf. Figure 4). In that case, θ is still small enough to justify the neglect of radial acceleration component. On the other hand, the Earth’s acceleration is still low enough such that the the evolution of the mean semi-major axis is strictly linear, i.e. $a(t)={a}_0+\dot{a}t$, with $\dot{a}=\mathrm{const}$. The evolution of the mean anomaly is computed by integrating the equation of the mean motion:

$$ {\displaystyle \begin{array}{c}M(t)={M}_0+{\int}_0^t\sqrt{\frac{\mu_{\mathrm{Sun}}}{a{(t)}^3}} dt\\ {}=\sqrt{\mu_{\mathrm{Sun}}}\left(\frac{2}{\dot{a}\sqrt{a_0}}-\frac{2}{\dot{a}\sqrt{a_0+\dot{a}t}}\right)\end{array}}, $$

(9)

Fig. 4

Evolution of mean semi-major axis, a (left), mean Earth displacement angle, θ (middle) and mean Earth range, r_Earth (right), according to approximate equations of section 2.4 (blue) for four different MIDA values: −20^∘, −16^∘, +16^∘, +20^∘ (from top to bottom). The green lines show the comparison with a numerical propagation. The dashed green lines show the true Earth range for the plots on the right

Subtracting this from the linear time-evolution of the Earth’s mean anomaly, yields the time evolution of the mean Earth displacement angle, θ:

$$ {\displaystyle \begin{array}{c}\theta (t)={M}_{\mathrm{Earth}}(t)-M(t)\\ {}={\theta}_0+\sqrt{\mu_{\mathrm{Sun}}}\left(\frac{2}{\dot{a}\sqrt{a_0}}-\frac{2}{\dot{a}\sqrt{a_0+\dot{a}t}}-\frac{t}{\dot{a}\sqrt{{\mathrm{AU}}^3}}\right)\end{array}}, $$

(10)

with $\dot{a}$ given by eq. 8. The mean Earth distance can be then computed from the geometric relation

$$ {r}_{\mathrm{Earth}}(t)=2a\sin \left(\theta (t)/2\right) $$

(11)

The cases where the time evolution of r_Earth, has a turning point are the most interesting ones, because they lead to the least difference between the minimum and maximum Earth distance during a given mission time: e.g. for a trailing configuration the initial semi-major axis is chosen slightly smaller than 1 AU causing an initial decrease of r_Earth. After the turning point is reached at 1 AU semi-major axis, the Earth distance increases again and reaches its maximum value at the end of the mission. In order to write eq. 10 in a form that illustrates this insight clearly, a Taylor expansion can be performed in both the semi-major axis change fraction, $\frac{\overset{.}{a}t}{a_0}$, and the deviation of the initial semi-major axis from 1 AU: ϵ=1 − a₀/AU

$$ \theta (t)={\theta}_0+\frac{2}{3}\sqrt{\frac{\mu_{\mathrm{Sun}}}{{\mathrm{AU}}^3}}\epsilon t-\frac{3}{4}\sqrt{\frac{\mu_{\mathrm{Sun}}}{{\mathrm{AU}}^5}}\overset{\cdotp }{a}{t}^2. $$

(12)

In this form it is clear that the evolution of the mean Earth displacement angle, θ, as a function of time is a parabola. The θ value reached after the science phase duration, Δt, is thus basically determined by the choice of the initial semi-major axis, a₀. Conversely, this allows computing the initial semi-major axis, a₀, from the maximum allowed mean Earth displacement angle, θ(Δt), which occurs at the end of mission, by solving eq. 12 for a₀:

$$ {a}_0=\mathrm{AU}\left(1-\frac{2}{3}\sqrt{\frac{{\mathrm{AU}}^3}{\mu_{\mathrm{Sun}}}}\frac{\theta \left(\varDelta t\right)-{\theta}_0}{\varDelta t}-\frac{\overset{.}{a}\varDelta t}{2\mathrm{AU}}\right). $$

(13)

θ(Δt) is computed from the maximum allowed mean Earth distance. Since for LISA the constraint on the maximum Earth distance is typically on the true Earth distance and not the mean one, the constant maximum difference between the Earth’s mean and true anomaly of almost 1^° has to be taken into account when computing θ(Δt). With some margin, a value of 1.2^° has been found to work well for the considered range of MIDAs.

$$ \theta \left(\Delta t\right)=\pm \left(2{\sin}^2\frac{r_{\mathrm{max}}}{2\mathrm{AU}}-1.{2}^{{}^{\circ}}\right), $$

(14)

where r_max is the maximum allowed true Earth distance. A graphical representation of eqs. 6, 10 and 11 for a number of different MIDA values and r_max = 65 · 10⁶km is shown in Fig. 4. The parabolic shape of the mean Earth displacement angle, θ, and mean Earth range evolution is clearly visible and the time where the turning point occurs depends on the choice of the MIDA. The maximum mean Earth range is always reached at the end of mission. The comparison with an exact numerical propagation of the LISA formation centre is also shown. For the numerical propagation the same final θ value is imposed as for the analytic solution, eq. 14. Note that even for the numeric solution the maximum true Earth distance is not exactly 65 · 10⁶km. This is because the constraint is put on the mean Earth distance in this section whereas the actual constraint is on the true Earth distance. This discrepancy can be refined during the actual numeric science orbit optimisation.

Equation 13 is particularly useful as it provides an initial guess for the initial semi-major axis of the cartwheel orbit and can also be used in conjunction with the Keplerian model from section 2.2 for defining the target orbit in the transfer optimization. For the sake of reference, Table 4 shows the initial semi-major axis values for a number of MIDA values and the comparison with numerically computed (exact) values. Note that for MIDA values below ±11^°, the constraint on the maximum Earth distance cannot be met any more because initial Earth perturbation is too strong.

Table 4 Optimal values for initial mean semi-major axis according to eq. 13 and numerically computed value as a function of MIDA, θ₀. The maximum true Earth distance assumed to compute these values is r_max = 65 · 10⁶km

The Earth’s gravity also has a strong effect on the breathing of the cartwheel and significantly alters the evolution shown for the Keplerian model in Fig. 3. This breathing can only be controlled by a careful adjustment of the initial conditions which will be described in the following section.

The Cartwheel Formation - Fully Numerical Model

In order to design a realistic cartwheel trajectory, a fully numerical model is required. The dynamical model used here takes into account the point mass gravity of the Sun and all relevant planetary bodies (Mercury, Venus, Earth, Moon, Mars, Jupiter, Saturn). Solar radiation pressure (SRP) and other non-gravitational accelerations are not taken into account as they will be compensated by LISA’s Drag Free and Attitude Control System (DFACS) [16,17,18]. This is because the spacecraft trajectories follow the trajectories of the internal test masses which are influenced by gravitational forces only. The effect of the self-gravity acceleration will be discussed later in section 3.2.

Optimisation for 10 Years of Science Phase

The six-parameter Keplerian cartwheel model described in section 2.2 serves as an initial guess for the optimisation. All these parameters, except for the clocking angle, σ₀, are considered “design parameters” here. The assumptions for the values of theses parameters used in this section are shown in Table 5. Note that the scientific requirements allow for a range of MIDAs (including Earth-leading configuration) and also for the counter-clockwise option. These options have been analysed internally at ESA, but will not be presented in this paper for brevity. The clocking angle, σ₀, is considered fixed in the current section as well, but will be used as a free parameter during the transfer optimisation described in section 4.2. It is assumed (and verified a-posteriori) that the success of the cartwheel optimisation does not significantly depend on the initial guess of σ₀.

Table 5 Assumptions for cartwheel design parameters serving as initial guess

The objective of the cartwheel optimisation is to adjust the initial conditions such that the breathing of the operational triangle is controlled within the allowed parameter ranges over 10 years mission time. The optimisation is carried out by allowing all 18 parameters of the initial cartwheel state (3 × 6 Keplerian state parameters) to be adjusted by the optimiser within a user-defined narrow band around the first guess. The initial epoch is fixed and not optimised. There are various ways to formulate the problem. Here the approach of using a constant cost function and adding the stability requirements as constraints has been used. Thus, it is a pure feasibility problem. The assumed constraint definitions are summarised in Table 6.

Table 6 Constraints in the cartwheel optimisation problem

The propagation of the initial state was carried out using a Runge-Kutta (8)7 dense stepper [19, 20]. It also propagates the state transition matrix which is used together with an in-house automatic differentiation software to obtain analytical partials for the optimisation. For the actual optimisation SNOPT [21] is called via the optimisation framework PyGMO [22] from Python. The minimum and maximum values over the mission time used for the constraints evaluation (Table 6) are obtained from sampling the trajectory with a step size of 10 days.

The resulting evolution of the cartwheel geometry is shown in Fig. 5. It can be seen that the corner angles and arm length rates can be constrained within the required windows over the considered mission time of 10 years. The evolution of the Earth range shows a characteristic profile where the maximum distance of 65 · 10⁶ km to the formation centre is reached at the end of mission. It is a result of the initial semi-major axis choice (cf. section 2.4) which leads to an initial drift towards Earth followed by a turning point. Note that the initial semi-major axis determines the Earth distance profile and thus the overall impact of the Earth perturbation on the formation stability. Relaxing the maximum allowed Earth distance decreases the Earth perturbation and thus allows choosing more narrow windows for the corner angles and arm length rates. E.g. by increasing the maximum Earth distance constraint from 65 · 10⁶ km to 75 · 10⁶ km it is possible to reduce the corner angle variations from ±1.0^° to ±0.75^°.

Fig. 5

Evolution of the cartwheel formation over 10 years in a fully numerical model

Figure 6 shows the evolution of the Keplerian elements during the 10 years science phase. The drift of the semi-major axes due to the Earth’s influence is apparent. Moreover, the initial values of the other elements are seen to be close to the initial guess determined by the Keplerian model (cf. section 2.2), but drifting away due to the third body perturbations.

Fig. 6

Evolution of the cartwheel orbital elements in the ecliptic frame over 10 years in a fully numerical model

The initial Keplerian parameters of the cartwheel orbit shown in Figs. 5 and 6 are given in Table 7 for reference.

Table 7 Initial osculating keplerian state in EME2000 for the cartwheel orbit obtained in the fully numerical model. The reference epoch is 2035–08-15 T12:00:00 TDB

Naturally, the chosen MIDA also has a strong impact on the obtainable corner angle variations, because it determines the initial strength of the Earth’s gravitational perturbation. A parametric analysis on different MIDA values is shown in Fig. 7. For MIDA values, ∣θ₀ ∣  < 14^° it was very difficult to achieve convergence at all therefore the solutions don’t appear in the plot.

Fig. 7

Achievable corner angle variations for cartwheel orbits with different MIDA values. The mean arm length and maximum Earth distance are kept at 2.5 · 10⁶ km and 65 · 10⁶ km, respectively

Spacecraft Self-Gravity

The analysis in the previous section neglected one important dynamical effect, namely the spacecraft self-gravity. This is an effective acceleration component acting on the spacecraft resulting from the DFACS following the motion of the two internal test masses which are not located precisely in the centre of mass of the spacecraft. The DFACS is the control system in LISA used to establish drag-free operation: a housing around the test masses shields them from external non-gravitational forces, like SRP, and senses the relative position of test mass and spacecraft. A control system then commands the spacecraft thrusters to follow the free-falling mass. This is required since the test masses are allowed to be influenced by gravitational forces only. More information on the DFACS can be found in [16,17,18].

To understand the effect of self-gravity, first consider a single test mass which is offset from the spacecraft centre of mass. The mass of the spacecraft will exert a small gravitational force on the test mass and will cause it to move towards the spacecraft centre of mass. This motion will be detected by the internal interferometers (between the test mass and the spacecraft) and the DFACS will command the thrusters to accelerate the spacecraft into the opposite direction such that the test mass effectively does not move with respect to the spacecraft. This acceleration produced by the DFACS has an impact on the spacecraft trajectory and needs to be modelled in the trajectory optimisation.

Now LISA does not have one, but two test masses per spacecraft and the situation becomes slightly more complicated: since the two test masses are at different locations with respect to the spacecraft centre of mass, there will also be a relative gravitational (gradient) acceleration between them (in addition to the common-mode acceleration described above). It is thus not possible for the DFACS to follow both test masses at the same time if both are free-falling. The solution is to only have the test masses free-falling along one space dimension (along the laser arm). The self-gravity acceleration on the test masses perpendicular to their laser arm directions is compensated using electrostatic actuators. The reaction forces of these electrostatic actuations on the spacecraft also need to be compensated by the DFACS. This adds another component to the net spacecraft acceleration and needs to be taken into account in the orbit propagation. To summarize, there are two self-gravity acceleration components: a common-mode component, which is also present with only one test mass, and a differential component which results from a relative acceleration between the test masses.

The total self-gravity acceleration depends on the spacecraft configuration which is not known at this point of the project. Theoretically, it can be directed along any of the three spacecraft axes. For simplicity, an analysis has been conducted with a self-gravity acceleration only towards the cartwheel centre (cf. Figure 8). A more detailed Monte-Carlo analysis allowing any direction for the acceleration will be presented in section 5.5.

Fig. 8

Schematic representation of the LISA formation and direction of the self-gravity acceleration. In reality the self-gravity acceleration can be directed in any of the three spacecraft axes

Moreover, due to fuel depletion the magnitude of the acceleration is not constant during the mission time. A detailed model of the acceleration history is yet to be developed and is not subject of the present paper. For simplicity, the current assessment considers acceleration levels changing from −2 nm/s² (start of mission) to +2 nm/s² (end of mission) where a positive value means an acceleration towards the cartwheel centre. Also the range of −4 nm/s² (start of mission) to +4 nm/s² (end of mission) has been looked at for comparison. The state from Table 7 has been propagated taking into account these accelerations and a comparison of the corner angles and arm lengths rate evolution is shown in Figs. 9 and 10, respectively.

Fig. 9

Corner angles evolution (top) and arm length rate evolution (bottom) for added self-gravity acceleration linearly varying from −2 nm/s² to +2 nm/s² in the direction of the cartwheel centre

Fig. 10

Corner angles evolution (top) and arm length rate evolution (bottom) for added self-gravity acceleration linearly varying from −4 nm/s² to +4 nm/s² in the direction of the cartwheel centre

The first case shows a mild violation of the corner angle constraints of ±1.0^° allowed variation. In the second case the violation becomes more severe. Note that no further optimisation has been conducted here. Once a more detailed model of the self-gravity acceleration history is available, it can be implemented as part of the optimisation procedure described in section 3.1. It is expected that accelerations of the levels considered here can be compensated for by a further adjustment of the initial cartwheel state and a stable formation within the requirements can be achieved.

It is considered more problematic if the self-gravity acceleration has a significant unknown contribution (along all three axes). Such a contribution can lead to a violation of the formation constraints, because it cannot be taken into account operationally in the final manoeuvre optimisation at the time of cartwheel orbit insertion. Its effect will be assessed in section 5.5.

Station Keeping

Nominally, no station keeping manoeuvres are foreseen during the 10 years science phase. The operational orbit and insertion conditions are required to be designed such that the requirements on the corner angle variations and arm length rates are fulfilled. Whether this is possible or not mainly depends on the chosen MIDA value as was shown in Fig. 7. The plot shows that if the corner angles are required to stay within a ±1.0^° window, the MIDA cannot be below ±20^°. However, going to lower MIDA values would have the advantage of a reduced transfer Δv (cf. section 4.2). Therefore, if MIDA values lower than ±20^° shall be employed, station keeping is required. Note again that these conclusions hold for a maximum allowed Earth distance of 65 × 10⁶ km. For larger values of that constraint it should be possible to reduce the MIDA below ±20^° without requiring station keeping.

Different station keeping strategies are possible. Manoeuvres imply an interruption of the science operations, and require a re-acquisition of the formation. Therefore, it is desirable to minimize the number of manoeuvres and also perform them simultaneously with all 3 spacecraft, if possible. Gaps in the science operations of longer than one week per year are not permitted. The simplest station keeping strategy that complies with these requirements is to place simultaneous manoeuvres after regular time intervals. The minimum number of manoeuvres required to stabilize the orbit for 10 years as a function of the MIDA value is shown in Fig. 11 (left).

Fig. 11

Number of required station keeping manoeuvres (left) and total station keeping Δv (right) as a function of the MIDA

The right panel in that figure shows the total required station keeping Δv. For each manoeuvres three free parameters were introduced in the cartwheel optimisation process. The corner angles were required to stay in a window of ±1.0^°. No additional objectives like change of the Earth range drift rate were imposed. The manoeuvres were equally spaced in time during the 10 years science duration, i.e. for the case of 2 manoeuvres, they take place 3.33 and 6.66 years after science orbit insertion. The manoeuvre times were not optimised and for each case the minimum number of manoeuvres were used that achieved the objective of maintaining the orbit within the requirements. Note that there is an outlier for the θ₀ =  − 14^° case. This is most likely due to the simple manoeuvre strategy where the manoeuvre times are not optimized. For all the other cases a station keeping budget of 10 m/s per spacecraft is sufficient. In case the Δv budget requires going below a MIDA of ±20^°, it is recommended to not go lower than ±16^°, which still is possible with a single station keeping manoeuvre after 5 years. Note however, this does not yet include perturbations that come from science orbit insertion accuracy. That contribution has been analysed in more detail in section 5.4 and 5.5.

Solar Electric Propulsion Transfer

The current section is going to describe the transfer from launch to cartwheel orbit insertion using Solar electric propulsion (SEP). Since the transfer Δv depends on the cartwheel clocking angle at the time of insertion (cf. section 4.2), both transfer and science phase cannot be strictly separated. Section 4.2 will describe the method that was used to take the interconnection into account.

Direct Escape Launch

A joint launch of all three spacecraft from Kourou into direct escape is the current baseline for ESA. Currently, an all-year launch period is targeted. In order to limit the number of required launcher programs to one, a single separation state (defined in the Earth-fixed) frame has been assumed. Since at the time of this analysis the trajectory analysis from the launch service provider is still under consolidation, a GTO-like escape has been assumed. For simplicity, the assumed orbital parameters at virtual perigee² are defined in the inertial EME2000 reference frame and are summarised in Table 8. The RAAN is a free optimisation parameter and can be fixed by the appropriate choice of the lift-off hour.

Table 8 Assumptions on the escape trajectory in EME2000 reference frame

After separation from the launcher upper stage, the three spacecraft follow independent trajectories, which initially will still be closely grouped.

The SEP is assumed to be used at the earliest 30 days after launch in order to allow for sufficient commissioning time. The actual start time of the first manoeuvre follows from the parametric optimization process of each orbit and can take place much later than 30 days after launch.

Simultaneous Optimisation of Transfer and Science Orbit

The simplest strategy would be to completely decouple the transfer optimisation from the cartwheel orbit optimisation and assume a transfer to a fixed cartwheel orbit. However, the transfer Δv is expected to depend on the chosen cartwheel orbit, in particular the initial clocking angle, σ₀. This dependency for the February 2034 launch will be discussed in more detail in section 4.3 and shown in Fig. 16. Since an end-to-end optimisation in the full dynamical model may be difficult to implement in a robust way, the following three-stage strategy has been employed:

• 1. Assume a Keplerian cartwheel orbit and optimise the transfer with the clocking angle, σ₀ as a free parameter.

• 2. Optimise the cartwheel orbit with full dynamics for the obtained σ₀ value.

• 3. Re-optimise the transfer assuming the obtained cartwheel orbit as fixed.

The change in transfer Δv between stage 1 and stage 3 is expected to be small which will be confirmed later on. Stage 2 has already been discussed in section 3. The following will focus on stage 1.

The individual manoeuvre time-line is optimized for each spacecraft such that its target orbit is met with minimal propellant expenditure. All three spacecraft need to be optimised together since they share the launch vehicle and therefore their (free) launch RAAN needs to be the same. A transfer duration of less than 540 days is targeted. This constraint implies that transfers with more than one revolution around the Sun cannot be used. During SEP operations, the Sun aspect angle (SAA, defined as the angle between the direction of the current thrust acceleration vector and the direction from the spacecraft towards the Sun) must remain within a range of 90^° ± 40^°. This is required in order to guarantee sufficient illumination of the solar arrays, knowing that the solar array normal is perpendicular to the thruster direction.

The starting assumptions for the SEP transfer to the cartwheel orbits are listed in Table 9. In order to account for outages during thrust arcs due to ground communications and contingencies, a duty cycle of 90% has been assumed. This has been modelled simply as a factor on the nominal available thrust level. Hereafter always the effective thrust level (i.e. nominal thrust times duty cycle) will be used to label cases.

Table 9 Assumptions for SEP transfers

The assumed wet mass is compatible with an Ariane 64 launcher performance of 7000 kg into direct escape with an infinite velocity of 300 m/s [23]. The maximum Δv of the three spacecraft is used as an objective function during the optimisation process. This number is relevant for the tank sizing since all three spacecraft are going to be built identically.

The transfer has been transcribed to a multiple shooting problem allowing up to four thrust arcs per spacecraft. The thrust direction of each thrust arc is parametrised by two angles and their rates which are free parameters. The thrust magnitude has been assumed constant. WORHP [24] has been used from PyGMO [22] as an optimisation algorithm. Initial guesses for the optimisation have been generated based on simple Lambert arc solutions with impulsive manoeuvres. Starting from these, a series of optimisation runs have been conducted for each trajectory, successively decreasing the thrust magnitude to the target value of 81 mN. Moreover, the transfer optimisation has first been done for each spacecraft individually keeping the launch RAAN and the clocking angle fixed to the initial guess value. In a second step the launch RAAN and the clocking angle have been freed and all three spacecraft have been re-optimised together to arrive at a globally optimal transfer solution.

February 2034 Launch to MIDA= − 20^° Cartwheel

To illustrate the SEP transfer in detail, the launch epoch of 2034-02-21 12:00 is chosen as an example. In a later analysis of 12 representative launch epochs per year the February launch turns out to be the sizing case for 2034 in terms of Δv. The targeted operational orbit is one with a MIDA of −20^° with parameter assumptions as shown in Table 5. Up to four thrust arcs are allowed for the 1-revolution transfer. The SEP transfer trajectory for an effective thrust of 81 mN is shown in Fig. 12 (left) in the Ecliptic Reference Frame of epoch 2000. The right plot in the Earth local orbital frame illustrates the path that LISA takes w.r.t. Earth. The Earth local orbital frame is defined as follows:

• X-axis along the Sun-Earth line

• Z-axis along the Earth angular momentum vector

• Y-axis completing the right-handed system

Fig. 12

Trajectory plots for the February launch transfer in the Ecliptic Frame (left) and the Earth local orbital frame (right). The Earth orbit is displayed in blue, thrust arcs are red and the arrival points are indicated by red circles

For the Earth-trailing option the transfer will always first increase the semi-major axis in order to achieve the correct phasing. This implies that the spacecraft can pass through eclipses by the Earth in the early transfer phase especially for launch dates close to the equinoxes. The transfer in February is eclipse-free. The evolution of the orbital elements for the March launch are shown in the Fig. 13. The arrival state fulfils the conditions described in section 2.2. The evolution of the geometry w.r.t. to the Sun and Earth is shown in Fig. 14. The Earth eclipse margin is an angle quantifying the eclipse condition. If this angle becomes negative, the spacecraft enters eclipse. The thrust aspect angles evolution (middle right) show that the angular rates required for the spacecraft are very low and thus no operational issues are expected.

Fig. 13

Evolution of the osculating orbital elements for the February launch option

Fig. 14

Evolution of the Earth, Sun and thrust direction geometry for the February launch option

An important aspect for communications is the pointing direction of the High-Gain Antenna (HGA). For this analysis it is assumed that the antenna is mounted on the opposite side from the solar array which is facing the spacecraft +Z axis. Therefore, if the HGA elevation becomes too large, the field of view is obstructed by the spacecraft. A maximum allowed elevation of +25^° is assumed for the time being. The attitude of the spacecraft during transfer is assumed as follows:

• During thrust arcs the acceleration vector is always along the spacecraft +Y axis.

• The Sun incident angle on the spacecraft +Z axis (solar array) is maximized. For coast arcs this means that the Sun direction is always along the +Z axis.

• The spacecraft X axis completes the right-handed system.

During coast arcs this attitude is not completely fixed because a freedom to rotate the spacecraft around the Z axis remains. This implies that the HGA azimuth is not fixed. Figure 15 shows the HGA azimuth and elevation evolution during transfer. It is clearly visible that the prohibited elevation range above +25^° is violated for a large part of the transfer. This is a general feature of transfers to Earth-trailing orbits, because for the initial part of the transfer both Sun and Earth are in the same direction as seen from the spacecraft. This is can be seen in Fig. 12 (right). Various options for solving this problem are discussed, e.g. interruption of thrust arcs for Earth communications or use of Low Gain Antennas. For transfers to an Earth leading orbit the HGA elevation is usually not a problem.

Fig. 15

HGA pointing evolution for the February launch option

Sensitivity on the Initial Cartwheel Clocking Angle

As already mentioned in section 4.2, a key question for the design of the transfer trajectory is whether the operational orbit can be assumed as fixed for all launch dates or whether a custom operational orbit has to be designed as a function of the launch date. The MIDA, arm length and argument of perihelion are considered as (fixed) design parameters. The initial semi-major axis is mainly fixed by the maximum Earth distance constraint (cf. section 2.4) and the inclination parameter, δ is used to adjust the stability of the cartwheel geometry within narrow limits. Therefore, the only significant impact on the transfer Δv is expected from the initial cartwheel clocking angle, σ₀.

In order to assess the impact of this parameter on the transfer, a series of transfers similar to the one described in section 4.3, have been computed, but varying the fixed value for the initial clocking angle. Figure 16 (left) shows the transfer Δv of the three LISA spacecraft as a function of the initial clocking angle. A clear pattern is recognisable, which basically repeats after 120^° due to the symmetry of the formation. The difference between the largest and smallest maximum Δv is about 80 m/s. This shows that the clocking angle has to be optimized together with the transfer, in order to achieve the optimal Δv for the propellant tank sizing. The plot also shows that an indication of the clocking angle being optimal is that the Δv values of two LISA spacecraft is identical while the Δv for the remaining spacecraft is below that value. Conversely, if all three spacecraft have different transfer Δv, it is an indication that the clocking angle has not been chosen optimally. For the transfer described in section 4.3 the clocking angle has been optimised along with the transfer and therefore the maximum transfer Δv is close to 1100 m/s.

Fig. 16

Transfer Δv as function of the arrival orbit initial clocking angle. Individual LISA spacecraft (left) and sum of all LISA spacecraft (right)

Figure 16 (right) shows the sum of the Δv of all three LISA spacecraft as a function of the clocking angle, which is relevant for computing the total wet mass on the launcher. Here the variation is less prominent, but still amounts to about 80 m/s for between the worst and the best clocking angle value. Note that the minima for the Δv sum (right) don’t occur at the same clocking angle as for the maximum Δv value (left).

Full Year 2034 Launch Period

Over the year a seasonal variation in the transfer Δv is expected mainly for two reasons: the true anomaly of the Earth at departure and the relative geometry of the Earth equator to the ecliptic. The total Δv is shown in Fig. 17 (left). The sizing case happens in February which amounts to 1092 m/s without margin (cf. section 4.3). Note that the difference in Δv between optimisation stages 1 and 3 (cf. section 4.2) was found to be below 12 m/s for all months in 2034. For the February launch the difference is 7 m/s. Therefore, even running only stage 1 yields an excellent estimate on the transfer Δv. The transfer Δv depends on how much the initial ecliptic RAAN has to be adjusted during the transfer of the individual spacecraft. Since the RAANs of the three spacecraft in the target orbit are separated by 120^°, there can be significant differences between the three spacecraft. Note that the labelling LISA1/LISA2/LISA3 is arbitrary and can change between two months.

Fig. 17

SEP Transfer Δv (left) and final mass (right) for 2034 MIDA= − 20^∘ clockwise. The crosses indicate actually computed cases. The lines are only drawn to guide the eye

The total transfer duration varies between 440 days and 540 days over the year and between the three spacecraft as shown in Fig. 18 (left). The total thrust-on time fraction, shown in Fig. 18 (right), follows the same functional behaviour as the Δv and has a maximum of about 70% in February. Reducing the thrust level below the 81 mN is therefore expected to make these transfer opportunities more difficult to achieve.

Fig. 18

SEP transfer duration (left) and percentage of thrust-on time (right) for 2034 MIDA= − 20^∘ clockwise

In order to illustrate the difference of the instantaneous initial Earth displacement angle versus the MIDA, these two parameters are plotted for the whole launch period in Fig. 19. As explained in section 2.3, the initial Earth displacement angle (left) experiences a variation of about ±2^° over the year while the MIDA stays constant at −20^°. The slight deviation of the MIDA from −20^° between the different launch opportunities stems from the fact that for the sake of simplicity the value of MIDA= − 20^° was imposed 400 days after launch for all launch dates. Since the transfer duration varies between the different launch opportunities, the actually achieved MIDA has a (negligible) variation.

Fig. 19

Instantaneous initial Earth displacement angle (left) and MIDA (right) for 2034 MIDA= − 20^∘ clockwise

The ranges of Sun aspect angles are always within the allowed range between 50^° and 130^°, as shown in Fig. 20. In most cases either the lower or the upper constraint is active.

Fig. 20

Ranges of Sun aspect angles of the acceleration vector for SEP transfers in 2034 MIDA= − 20^∘ clockwise

In order to assess the dependency of the Δv on the thrust capability of the spacecraft, several cases with a reduced thrust level have been analysed. The results for (effective) thrust levels 81 mN, 72 mN, 64 mN and 54 mN are shown next to each other in Fig. 21. Besides of increasing the required Δv, the effect of the lower thrust is also to reduce the number of launch opportunities. For 54 mN it is quite difficult to achieve convergence at all. Therefore, this option is not very robust and thus not recommended.

Fig. 21

Total Δv for different thrust level values for 2034 MIDA= − 20^∘ clockwise

So far only transfers to an Earth-trailing orbit have been discussed. The analysis for a Earth-leading orbit has also been conducted but will not be presented in detail here for the sake of brevity. The main difference in the transfer Δv is the change in the seasonal variation: The sizing case for a MIDA of +20^° occurs in August instead of February and amounts to 1253.3 m/s. This fact can be exploited in order to reduce the overall sizing Δv by combining launches to trailing and leading orbits depending on the launch month: In the Summer months a launch to a trailing orbit is preferable while in the Winter months a leading orbit yields the lower transfer Δv [12].

Navigation and Insertion Accuracy

This section is presenting all the analyses related to orbit navigation and guidance of the LISA spacecraft during transfer and cartwheel orbit insertion. A key question to be answered here is: how precisely can the three LISA spacecraft be inserted into the cartwheel formation and how does the error impact the stability of the cartwheel? The first part of that question is going to be answered by the presented navigation analysis, the second part by a Monte-Carlo analysis that uses the output dispersion matrix from the navigation as an input. The idea of a navigation analysis is to simulate all ground-based (and potentially inter-spacecraft) measurements as well as orbit correction manoeuvres that will take place during the operations phase. Generally speaking, there are two key outputs of such an analysis that quantify the orbit accuracy that can be reached:

• 1. Knowledge: the difference between the true spacecraft state and the estimated spacecraft state. This is mainly a function of the measurement accuracy and geometry of the ground-based Range and Doppler observations, but is also influenced by noise processes like the SEP acceleration or SRP.

• 2. Dispersion: the difference between the true spacecraft state and the reference spacecraft state. This is mainly a result of manoeuvre execution errors, but is also influenced by the state knowledge: the dispersion can never get lower than the knowledge at a given epoch, because the precision to which manoeuvres are computed by Flight Dynamics can never exceed the state knowledge on which these computations are based.

Both figures are generally described by covariance matrices and are typically shown in the local orbital frame of the spacecraft.

A navigation analysis can be broken down into two parts: the covariance analysis, which mainly deals with the knowledge and the guidance analysis which mainly deals with the dispersion. Since, however, both are inter-linked, iterations between the two are needed until convergence of all covariance matrices and the manoeuvre budget is reached.

All analyses in this chapter are based on the reference trajectory with launch in March 2034. Similar results are expected for any other transfer.

Covariance Analysis Assumptions

The covariance analysis uses a Square-Root-Information Filter (SRIF) to process all available observation types and compute the covariance matrix of all estimated parameters. Some uncertain parameters in the problem are not estimated (i.e. their a-priori covariance is chosen not to be affected by the observations processing), but their uncertainty still needs to be taken into account. These parameters are called consider parameters. A summary of all uncertain parameters is given in Table 10.

Table 10 Summary of all uncertain parameters that appear in the covariance analysis

Observations are used by the filter to improve the knowledge of the estimated parameters. Currently, only ground-based Range and Doppler measurements are assumed. All observation assumptions are summarized in Table 11.

Table 11 Summary of all observation assumptions in the covariance analysis

All other assumptions related to the covariance analysis are summarized in Table 12.

Table 12 Other assumptions related to covariance analysis

In order to simulate the orbit determination process done in operations as closely as possible, a sliding window is used as an orbit determination arc. This means for every solution epoch measurements 14 days before the corresponding data cut-off epoch are processed (4 days for Range and Doppler). The obtained covariance matrix is then mapped forward to the solution epoch. The process is repeated for the considered range of solution epochs and the resulting knowledge evolution is plotted in a graph. This procedure is illustrated in Fig. 22.

Fig 22

Illustration of orbit determination process using a sliding window

Guidance Analysis Assumptions

The purpose of the guidance analysis is to determine to which extent the spacecraft state dispersion can be reduced by manoeuvres using the state estimation from the covariance analysis as input. To this end, guidance manoeuvres are scheduled that target the reference state at a pre-defined epochs. For simplicity, the current implementation of the analysis assumes impulsive guidance manoeuvres although LISA is only equipped with SEP. This approach is justified since the resulting Δv per manoeuvre is expected to be small, which is confirmed a-posteriori. Therefore, these guidance manoeuvres can also be performed by the SEP with roughly the same Δv by simply extending the burn duration. During the nominal trajectory design a duty cycle of 90% on the thrust availability has been taken in order to accommodate (amongst other things) such guidance manoeuvres even if they happen during a deterministic thrust arc (cf. Table 9). If due to extended guidance thrust arcs tracking periods fall into the time of thrusting, a different ground station scheduling can be considered. If that’s not possible, a deterioration of the tracking performance is expected. This is, however, only expected to be a problem for the large correction manoeuvres, TCM1 and TCM3. E.g, a 1 m/s guidance manoeuvre will have a duration of about 6 h when executed with the SEP which can easily be accommodated by an adaptation of the ground station schedule.

Table 13 summarises the general assumptions for the guidance analysis.

Table 13 Assumptions for guidance analysis

For the cartwheel insertion analysis the Trajectory Correction Manoeuvre (TCM) schedule is centred on the cartwheel arrival date, which is defined at the end of the last thrust arc. Since the last thrust arc ends at a different epoch for all three spacecraft, we use the latest date of the three as a reference, which is 2035-09-12 12:00:00 for the reference trajectory used here. In order to accommodate enough time for ground station observations and ground processing between the manoeuvres, they are scheduled every 14 days as shown in Table 14. For simplicity, currently all three spacecraft are assumed to execute the TCMs at the same epochs.

Table 14 TCM schedule for cartwheel insertion

Figure 23 shows a graphical representation of the manoeuvre schedule for the three spacecraft, also showing the SEP thrust arcs. This illustrates that only for LISA2 the first two TCMs happen during an SEP thrust arc. Due to SEP noise these manoeuvres are expected to be much less efficient than for LISA1 and LISA3 where they happen during a coast arc. In that sense LISA2 will represent the sizing/worst case in terms of insertion accuracy. On the other hand, the analysis for LISA1 and LISA3 will show how much can be gained by not doing any guidance manoeuvres during thrust arcs for the cartwheel insertion. The targeting strategy is based on the idea that in order to match both position and velocity of the cartwheel state at the end of the sequence, at least two TCMs are needed: the first one will match the position only and the second one will match the velocity. However, due to manoeuvre execution error, the match will never be exact. Therefore, two manoeuvres are assumed for each target point, a larger and a smaller one. Finally, TCM5 is an additional manoeuvre to clean up any residual velocity dispersion.

Fig. 23

Illustration of cartwheel insertion TCM strategy. The green arrows indicate the targeting point of the TCMs. TCM5 targets the velocity vector at its epoch

Navigation Results

Figure 24 shows the knowledge evolution of both position and velocity in the local orbital frame for all three LISA spacecraft (3 − σ values). The initially large covariance is reduced as soon measurements are taken (indicated by shaded vertical bars). This happens with a latency of 4 days due to the data cut-off. When no measurements are taken, the knowledge slowly degrades owing to the natural dynamics as well as the system noise. The impact of the noise is particularly strong for LISA2 which is the only spacecraft that has an SEP thrust arc in the considered time period. Jumps in the velocity knowledge are also caused by the guidance manoeuvres, TCM1 - TCM5. The along-track component in most cases is the one with the best knowledge. This is because of the measurement geometry: in a trailing/leading orbit ground-based Range and Doppler always measure the along-track component directly. The signature in an extended Doppler measurement arc imparted by the Earth’s rotation also measures the plane-of-sky components. However, the vertical component (= cross-track) is not measured well if the spacecraft is located close to zero declination w.r.t. the Earth’s equator, which is the case here. This is generally known as the zero declination problem. A way to significantly improve the plane-of-sky resolution, is to use ΔDOR measurements. As will be shown in section 5.4, the along-track knowledge is the driving one for the cartwheel stability. Therefore, ΔDOR measurements are not required here.

Fig. 24

Knowledge covariance evolution of all three LISA spacecraft for Range and Doppler (8 h per week per s/c)

The knowledge covariance is used to determine guidance manoeuvres in order to reduce the initially large dispersion. The post-manoeuvre dispersion resulting from the guidance analysis are shown in Fig. 25. One can see how the dispersion is gradually improved as more guidance manoeuvres are executed. The state components with the best knowledge also end up with the least dispersion.

Fig. 25

Post-TCM dispersion of all three LISA spacecraft for Range and Doppler (8 h per week per s/c)

Finally, Tables 15, 16 and 17 present the stochastic manoeuvre budget resulting from the guidance analysis. It is clear that TCM1 and TCM3 are the largest manoeuvres because they are the first manoeuvres for each of the two target points. TCM2, TCM4 and TCM5 are follow-up manoeuvres which merely clean up the mechanisation error of the previous manoeuvres. For the fuel tank sizing the 99% CL value shall be taken which amounts to 18.3 m/s per spacecraft. This value is still preliminary since the initial dispersion used here was guessed and does not originally come from the launcher dispersion and an end-to-end navigation analysis. This end-to-end navigation analysis shall also use continuous SEP guidance manoeuvres. Possibly, part of the corrections can then be absorbed into the deterministic thrust arcs.

Table 15 Stochastic manoeuvre budget for LISA1

Table 16 Stochastic maneuver budget for LISA2

Table 17 Stochastic maneuver budget for LISA3

Cartwheel Stability under Insertion Dispersions

A key aspect for the LISA operational orbit design is the question whether the formation stability can be maintained given the insertion accuracies. The current section uses the dispersion matrices at the end of the insertion sequence that are a product of the navigation analysis, section 5.3, in order to do a Monte-Carlo analysis. The three spacecraft states after the last TCM are sampled and propagated for 10 years. 10,000 samples are used.

The resulting dispersions are presented in Fig. 26. The left plot shows the distribution of maximum arm length rate vs. maximum corner angle deviation from 60^°. The distribution of the maximum arm length deviation from 2.5 · 10⁶ km vs. the time spent in a formation with at least one corner angle exceeding the 60^° ± 1^° range is shown on the right. Apparently, the deviation of the samples is well controlled for the key parameters (arm length rate and corner angle deviation). The duration spent in a formation violating the corner angle requirement is below one year in more than 99% of the cases. The maximum Earth distance is not affected at all by the insertion inaccuracies, as expected.

Fig 26

Formation stability for MIDA= − 20^∘ Range, Doppler case

A more quantitative summary of the analysis is given in Table 18. The 99% C.I. quantities are recommended to be used as reliable benchmarks of the expected variation. These don’t differ very much from the requirements of 10 m/s arm length rate and 1^° corner angle deviation. Note that the computation of time spent in a formation where at least one of the three corner angles is outside the 60^° ± 1^° range is based on a time discretisation step of 2 days. Therefore, the given number has an intrinsic error coming from the discretisation. The order of this error is 2 days times the number of intervals where a corner angle exceeds the allowed range.

Table 18 Formation stability statistics for MIDA= − 20^° Range, Doppler case

Cartwheel Stability Including Unknown Self-Gravity Accelerations

The picture changes if one adds the uncertainty from the acceleration component due to the spacecraft self-gravity. The effect of a known deterministic self-gravity acceleration has been analysed in section 3.2. In reality, the magnitude and direction of the self-gravity acceleration will depend on the spacecraft configuration, component placement tolerances and fuel depletion. This is not known at this point. In any case, the will be a remaining unknown component of the self-gravity acceleration even after the launch of LISA. The current section will analyse the impact of such an unknown component. For the current analysis it has been assumed that this unknown component has a Gaussian distribution with a standard deviation of 1 nm/s² per spacecraft axis. It has been assumed constant in the spacecraft-fixed frame during the mission time. This is an oversimplification and represents a conservative case until it can be better quantified how this acceleration changes over the mission time. The results presented in Fig. 27 and Table 19 and show a substantial increase both in the corner angle deviation and in the arm length rate compared to those in section 5.4.

Fig 27

Formation stability for MIDA= − 20^∘ Range, Doppler case with self-gravity accelerations

Table 19 Formation stability statistics for MIDA= − 20^° Range, Doppler case with added random self-gravity accelerations

As mentioned in the introduction, insertion dispersions currently violate the nominal preliminary requirements on corner angle variations and arm length rates (60^° ± 1.0^° and 10 m/s maximum). During the pending next iteration at system level it is going to be decided whether the spacecraft design allows for a relaxation of these requirements in order to fit the current trajectory design including dispersions. Alternatively, the trajectory design has to be adapted such that the 60^° ± 1.0^° and 10 m/s requirements are met even with insertion dispersions. There are two ways how this could be achieved:

• 3. Increasing the MIDA value will reduce the Earth perturbations on the constellation and thus decrease the corner angle variations and arm length rates (see also Fig. 7). The downside of this is an increase of transfer Δv.

• 4. Increasing the allowed maximum Earth range (e.g. to 75 × 10⁶ km) will also reduce the overall impact of Earth perturbations and allow for a compliant design. This will have an impact on the spacecraft communications system.

A combination of all options is, of course, also possible.

Cartwheel Stability per Component

Finally, in this section the importance of the initial spacecraft dispersion along the individual axes is analysed. To this end, instead of using the dispersion matrices coming from the navigation analysis, as was described in section 5.3, a fixed value along the different axes in the local orbital frame is assumed while the other axes covariances are set to zero. Additionally, a case with a spherical position or velocity dispersion is run to compare the other cases to. The resulting statistics are summarized in Table 20.

Table 20 Statistics of different cases with pre-defined initial dispersion along individual axes

From these results it is clear that the dominant contributions to the de-stabilization of the cartwheel are the radial position dispersion and the along-track velocity dispersion. These are the components which are associated with the semi-major axis. The cross-track component has the smallest contribution. Considering that ground-based Doppler observations mainly measure the along-track velocity component for LISA, this is a very favourable result. Added ΔDOR measurements will therefore offer only a limited benefit.

To understand the dynamics a little better, cases have also been analysed where the initial dispersion was set to zero in all six state components and only self-gravity acceleration along individual axes has been added. The results for the self-gravity accelerations defined in the local orbital frame are shown in Table 21.

Table 21 Statistics of different cases with zero initial dispersion and constant self-gravity accelerations along individual axes in the local orbital frame

It is clear from these results that again the main contribution to the cartwheel de-stabilization comes from the along-track acceleration component. In contrast to an initial along-track velocity dispersion, a constant along track acceleration has a much stronger impact on the formation stability. The radial and in particular the cross-track acceleration have a minor impact. Since, however, the self-gravity accelerations are tied to the spacecraft-fixed frame, the analysis was repeated defining constant accelerations in the frame where:

• X - along the formation centre

• Z - perpendicular to cartwheel plane

• Y - completing the right-handed frame

Results are shown in Table 22. Because the spacecraft frame is rotating with a period of a year, there is no constant along-track acceleration. Therefore, the overall impact is milder than observed in the local orbital frame. The X and Y components show the strongest contribution to the cartwheel de-stabilization.

Table 22 Statistics of different cases with zero initial dispersion and constant self-gravity accelerations along individual axes in the spacecraft-fixed frame

Conclusions

After a review of analytic models for the LISA cartwheel formation, a fully numerical optimisation analysis of both the transfer and science phase has been presented. The interdependency of both mission phases via the cartwheel clocking angle has been taken into account. The SEP transfer Δv varies depending on the launch month in 2034. But for the sizing month an allocation of 1092 m/s per spacecraft is sufficient with the assumptions taken. For the science orbit, with the assumed arm length of $2.5\mathrm{x}{10}^6$ km and MIDA=−20^°, a corner angle variation of close to 60^° ± 1.0^° during 10 years is feasible.

The expected cartwheel insertion accuracy has been estimated in a full navigation analysis where the most important sources of error have been taken into account. An accuracy of the order of $\mathcal{O}\left(10\ \mathrm{km},5\ \mathrm{mm}/\mathrm{s}\right)$ (along-track and radial) and $\mathcal{O}\left(100\ \mathrm{km},50\ \mathrm{mm}/\mathrm{s}\right)$ (cross-track) are achievable with Range/Doppler and several weeks of insertion sequence using 18 m/s per spacecraft. The radial position error and along-track velocity error are the driving ones for the stability of the subsequent science cartwheel orbit.

With the obtained insertion dispersion, a Monte-Carlo analysis has been conducted to analyse the impact on the corner angle variations during 10 years of science phase. If the self-gravity accelerations are perfectly known at the time of cartwheel insertion, the expected deterioration of corner angle variations is about 60^° ± 1.1^° at 99% C.L. violating the current requirement. The arm length rate is hardly impacted. However, if there is a remaining unknown constant component of the self-gravity accelerations of the order of 1 nm/s², the impact on the corner angle variations and arm length rate is significant. This underlines the importance of a good characterisation of the self-gravity accelerations prior to launch. Ways of dealing with the current requirements violation have been discussed.

In the frame of the currently ongoing Phase A, the following future Mission Analysis work is envisioned:

• Implementing the pending system level iteration on the cartwheel stability requirements.

• A better objective function for the cartwheel orbit optimisation: not applying strict constraints on the corner angles and arm length rates, but rather maximising the mission time spent within the allowed bands.

• Cartwheel orbit optimisation including spacecraft self-gravity.

• Relaxation of the maximum Earth distance constraint taking into account the actual achievable data rate as a function of the Earth distance.

• Implementation of a variable thrust model for the transfer optimisation. Instead of a constant thrust value, variations of the thrust as a function of Sun distance and SAA shall be taken into account.