Trajectory Design for the ESA LISA Mission

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 $\Delta 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 close to $60^\circ \pm 1.0^\circ$ 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^\circ \pm 1.1^\circ$ 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.

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 dragfree 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 6 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 (6 years nominal + 4 years extension) is guaranteed. The main 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 6 km ±2.5×10 5 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 6 km. This value, however, depends on the design of the spacecraft communications system. 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 chapter 2, a fully numerical simulation of the optimal cartwheel orbit is presented in chapter 3. Chapter 4 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 chapter 5. 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, chapter 6 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]: 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 Figure 1. T is the orbital period of the circular formation centre orbit and α i and i are integration constants. The relative trajectory described by Equation 1 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 ≡ α 0 and z = ± √ 3 x the relative trajectory becomes a circle around the formation centre which in 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 1 .
By placing three spacecraft on this circle with a relative phase of 120 • , an equilateral triangle formation with arm length d = 2 √ 3 x is achieved. Using the above equations, the heliocentric orbital elements of the three spacecraft can be derived. They are summarized in Table 1 . 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. Overall, there are five free continuous and one discrete free parameter in the system: semi-major axis, a 3. right ascension of the ascending node (RAAN) of LISA1, Ω 1 4. initial mean anomaly of LISA1, M 1 (t = 0) 5. argument of perihelion (discrete) ω 6. MIDA, θ 0 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 nonlinear 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 dimensionless parameter δ 1 defined by δ = αδ 1 , with α = d/2a. The optimum value is shown to be δ 1 = 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 i and the mean anomaly, 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, θ 0 defines the location of the centre of the cartwheel with respect to the Earth (see Figure 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, θ 0 . 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.

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. Assuming circular orbits and neglecting any radial and cross-track contributions, the following formula for the drift rate of the mean semi-major axis can be derived from the Gauss equations: 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 θ = θ 0 gives sufficiently accurate results if θ is not too far from 20 • . In that case it can also be assumed that the evolution of the mean semi-major axis is strictly linear, i.e. a(t) = a 0 +ȧt. Inserting this evolution of the mean semi-major axis into the (linear) time evolution of the mean anomaly and integrating over time, the time evolution of the mean Earth displacement angle, θ, can be derived: withȧ given by equation 6. The mean Earth distance can be then computed from 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 Equation 7 in a form that illustrates this insight clearly, it can be Taylor expanded in both the semi-major axis change fraction, at a 0 , and the deviation of the initial semi-major axis from 1 AU: : 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 0 . This allows computing the initial initial semi-major axis, a 0 , from the maximum allowed mean Earth displacement angle, θ(∆t), which occurs at the end of mission: 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 an 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.
where r max is the maximum allowed true Earth distance. A graphical representation of the equations in this section for a number of different MIDA values and r max = 65 · 10 6 km is shown in Figure 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, Equation 11. Note that even for the numeric solution the maximum true Earth distance is not exactly 65 · 10 6 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 10 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.
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 Figure 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]. 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.  Table 5 Assumptions for cartwheel design parameters serving as initial guess.

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, σ 0 , 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, σ 0 , 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 σ 0 .
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. There are various ways to formulate the problem. Here the approach of using a constant cost function and 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.
The propagation of the initial state was carried out using a Runge-Kutta (8)7 dense stepper [18,19]. 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 [20] is called via the optimisation framework PyGMO [21] 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 Figure 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 6 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 semimajor 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 6 km to 75 · 10 6 km it is possible to reduce the corner angle variations from ±1.0 • to ±0.75 • . 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.
The initial Keplerian parameters of the cartwheel orbit shown in Figures 5 and 6 are given in Table 7 for reference.
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 Figure 7. For MIDA values below ±14 • it was very difficult to achieve convergence at all therefore the solutions don't appear in the plot.

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. To understand this effect, 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 optical bench) and will command the DFACS to accelerate the spacecraft into the opposite direction such that the test mass effectively doesn't move with respect to the spacecraft. 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 acceleration between them (in addition to the commonmode 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  The mean arm length and maximum Earth distance are kept at 2.5 · 10 6 km and 65 · 10 6 km, respectively. 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.
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 2 (start of mission) to +2 nm/s 2 (end of mission) where a positive value means an acceleration towards the cartwheel centre. Also the range of −4 nm/s 2 (start of mission) to +4 nm/s 2 (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 Figures 9  and 10, respectively.
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 Figure 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 6 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 number of manoeuvres required to stabilize the orbit for 10 years as a function of the MIDA value is shown in Figure 11 (left).
The right panel in that figure shows the total required station keeping ∆v. 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. Note that there is an outlier for the θ 0 = −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 window 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 2 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.
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 after launch than 30 days.

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, σ 0 . This dependency for the February 2034 launch will be discussed in more detail in section 4.3 and shown in Figure 16. Since a 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, σ 0 as a free parameter. 2. Optimise the cartwheel orbit with full dynamics for the obtained σ 0 value. 3. Re-optimise the transfer assuming the obtained cartwheel orbit as fixed.
The change in transfer ∆v between state 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.
The assumed wet mass is compatible with an Ariane 64 launcher performance of 7000 kg into direct escape with an infinite velocity of 300m/s. 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 [22] has been used from PyGMO [21] 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 1revolution transfer. The SEP transfer trajectory for an effective thrust of 81 mN is shown in Figure 12 (left) in the Ecliptic Reference Frame of epoch 2000. The Earth orbit is displayed in blue, thrust arcs are red and the arrival points are indicated by red circles. 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 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 in eclipse-free. The evolution of the orbital elements for the March launch are shown in the Figure 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 Figure 14. The Earth eclipse margin is an angle quantifying the eclipse condition. If this angle becomes negative, the spacecraft enters eclipse. 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 Figure 12 (right). Various options for solving this problem are discussed, e.g. interruption of thrust arcs for Earth communications. For transfers to an Earth leading orbit the HGA elevation is usually not a problem.

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, σ 0 .
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. 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 window
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 Figure 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.
The total transfer duration varies between 440 days and 540 days over the year and between the three spacecraft as shown in Figure 18 (left). The total thrust-on time fraction, shown in Figure 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.
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 window in Figure 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.
The ranges of Sun aspect angles are always within the allowed range between 50 • and 130 • , as shown in Figure 20. In most cases either the lower or the upper constraint is active. 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 of 81 mN, 72 mN, 64 mN and 54 mN are shown next to each other in Figure 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.
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 [23].

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.
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.
All other assumptions related to the covariance analysis are summarized in Table 12.
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 Figure 22.

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 trajectory design a margin of 10% on the available thrust magnitude has been taken in order to accommodate (amongst other things) such guidance manoeuvres even if they happen during a deterministic thrust arc. Table 13 summarises the general assumptions for the 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.  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. 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.

Navigation results
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 Figure 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.
Finally, Table 15, 16 and 17 present the stochastic manoeuvre budget resulting from the guidance analysis. It is clear that TCM1 and TCM3 are the largest ma- noeuvres 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.

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 Figure 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 6 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 very well controlled for the key parameters (arm length rate and corner angle deviation). The duration spent in a formation fulfilling the corner angle requirement is below one year in most cases. The maximum Earth distance is not affected at all by the insertion inaccuracies, as expected.  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.

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 2 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 Figure 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.

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 is with a spherical position or velocity dispersion is run to compare the other cases to. The resulting statistics are summarized in Table 20.
From these results it is clear that the dominant contributions to the destabilization of the cartwheel are the radial position dispersion and the along-track velocity dispersion. These are the components which are associated with the semimajor axis. The cross-track component has the smallest contribution. Considering that ground-based Doppler observations mainly measure the along-track velocity To understand the dynamics a little better, also cases have 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.
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.

Conclusions
After a review of analytic models for the LISA cartwheel formation, a fully numerical optimisation analysis of both the transfer and science phase have 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 × 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 O(10 km, 5 mm/s) (along-track and radial) and O(100 km, 50 mm/s) (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.. 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 2 , 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.
In the frame of the currently ongoing Phase A, the following future Mission Analysis work is envisioned: -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.