The Aditya-L1 probe just (19th of September 2023) performed its last maneuver around Earth to escape its orbit, and is heading towards the Lagrange 1 point. When it gets there, the 1,500 kg (3,300 lb) satellite will use its seven science payloads to measure coronal heating, solar wind acceleration, coronal magnetometry, origin and monitoring of near-UV solar radiation (which drives Earth's upper atmospheric dynamics and global climate), coupling of the solar photosphere to the chromosphere and corona, energetic particle fluxes and magnetic fields of the solar wind.
I could not help myself but to try and reproduce the trajectory of the probe along its long journey from a highly elliptical orbit around Earth to its final destination (a Halo orbit around the Earth-Sun L1), 1.5 million kilometers away.
Optimal DeltaV trajectory for the Aditya-L1 mission.
The mission profile
Starting off with a launch on top of a PSLV-XL rocket on the 2nd of September 2023, the probe was placed on an Earth orbit with a perigee altitude of approximately 235 km and an apogee altitude of 19,500 km. According to the ISRO's Twitter account, the probe performed four (4) earth-bound maneuvers within a period of two weeks in order to increase its apogee, with the final orbit having a perigree altitude of 256 km and an apogee altitude of 122,000 km. From this highly elliptic orbit, it fired its engines once more and placed itself onto a trans-Lagrangian Point 1 trajectory.
Image credit: ISRO
This transfer is going to last approximately 110 days after which Aditya-L1 will start orbiting around the L1 point in a halo orbit.
How does a halo orbit around the L1 look like?
There is a lot of literature available on the Lagrange (also called libration) points, so I will not describe their properties in detail. In short, Lagrange points arise as equilibrium locations, when dealing with a circular restricted celesital three-body problem, i.e. two primaries that move in circles (e.g. Sun and Earth) and a third body with negligible mass (e.g. a spacecraft) moving in the gravitational field of the primaries without affecting them.
At the location of these points, the third body moves in a circular orbit with the same frequency as the primaries, so it appears as stationary in the rotating frame. To obtain the location of these points one has to find the equilibrium points of the circular restricted three-body problem set of ODEs, which describe the motion of the third body in the gravitational field of the primaries.
Once you go through the trouble of doing the math, you will see that there are three collinear points on the x-axis (the axis connecting the two primary bodies), which were discovered by Euler around 1750 and are denoted L1, L2 and L3 and there are two equilateral points discovered by Lagrange around 1760 and are denoted L4, L5.
Lagrange points in the Sun-Earth system (primary bodies not to scale)
Of special interest are the two points L1 and L2 closest to the secondary body (the Earth in this case), which are dynamically center-saddle points as can be shown with a linearized analysis. According to the Lyapunov theorem, there is a family of periodic orbits surrounding each of these points. There are planar periodic orbits called Lyapunov orbits, while their counterparts in the 3D problem are called Halo and Lissajous orbits.
For a detailed overview of the three-body problem that gives rise to the libration points, I can recommend the books Solar System Dynamics by Murray, The Three-Body Problem by Valtonen and the classic Theory of Orbits: The Restricted Problem of Three Bodies by Szebehely.
Halo orbits come in various sizes, as can be seen in the figure below for the Sun-Earth system. The Earth is located at (0,0,0) (illustrated by a tiny dot showing how small our planet is compared to the Halo orbit size), whereas the L1 point is shown with a small 'x'. As the size of the orbits is reduced, they start moving towards the Lagrange point and further away from the Earth, while at the same time, they start "tilting" towards the plane of rotation of the Earth.
This "tilting" can be expessed as a reduction in the vertical (out-of-plane) amplitude of motion Az compared to the horizontal amplitude Ax. This can be quantified in the figure below. At the same time, it can be also seen that the period of the orbit drops with increasing amplitude. The points in the graphs below correspond to the orbits shown above with the same color.
Any halo orbit can be characterized completely by specifying a particular out-of-ecliptic plane amplitude Az. Also, the Ax and Az amplitudes are constrained by a non-linear algebraic releshionship as shown in the left sub-figure below. From the same figure, one can find the minimum permissible value for Ax in order to have a halo orbit (Az > 0). In the case of a periodic solution about L1 in the Sun-Earth system, the min value of Az is about 200,000 km.
As a reference, the Aditya-L1 halo orbit has a period of 177.86 days and Az = 120,000 km according to this source.
For a detailed analysis of the Halo orbits, I would recommend going through Farquhar et al. (1973), Breakwell et al. (1979) and Howell et al. (1983).
There are both analytical and numerical ways available in order to construct a Halo orbit of a given amplitude/period. One of the simplest ways is to use an analytic method as derived by Richardson (1980). This is a third-order analytical solution, obtained by an application of successive approximations using the Lindstedt-Poincare approach.
Of course, this is only an approximate method and is insufficient for serious study of accurate motion near L1 or L2, but can give some qualitative physical insight into its properties.
If you want the "real deal", then one option is to use analytical approximations as a first guess and then use a differential correction process to generate a halo orbit accurate enough for mission design.
This numerical calculation starts with an initial condition and a target final state. Since halo orbits are symmetric about the XZ-plane (y = 0), and they intersect this plane perpendicularly (xdot = zdot = 0), a good candidate for an initial state is the vector X0 = (x0, 0, z0, 0, ydot0, 0). The initial values in this vector can be taken from the Richardson solution at the y=0 location.
The equations of motion are then integrated until the trajectory crosses the XZ-plane (y=0) again. A necessary condition for a periodic orbit is that this second crossing is also perpendicular, and so the desired target state vector needs to also satisfy the xdot_f=zdot_f=0 condition. Due to the inaccuracy of the third-order numerical method, the actual values for xdot_f and zdot_f will probably not be zero at the first crossing. In order to correct that, the three non-zero initial conditions (x0, z0, and ydot0) can be manipulated in an attempt to drive these final velocities xdot_f and zdot_f to zero.
To do so, one can use the transition matrix, which correlates a deviation in the final state vector (where we wanted to go vs where we ended up) to a change that needs to be applied in the initial state vector (how should we steer to get there). In essence, it is a form of sensitivity matrix of the final conditions (after half a period) to the initial conditions. With this correction applied, the revised initial conditions X0 are used to begin a second iteration and the process is continued until xdot_f=zdot_f=0. Usually, convergence to a solution is achieved within a few iterations.
For the Halo orbit of the Aditya-L1 mission, this corrective procedure looks like that:
Differential correction method for the halo orbit of Aditya-L1. Earth is the pale blue dot in the plot and the L1 point is shown with the 'x' marker. The converged periodic orbit (PO) is shown in black.
This is great, but how do I get there? Invariant manifolds (or "highways")
Describing and calculating the halo orbit is the first step. Now, which information from the dynamical system of the three-body problem can we exploit in order to get onto these orbits?
The answer is: stable invariant manifolds!
The stable invariant manifold of any periodic orbit is the set of all points of the phase space which tend asymptotically (as time reaches +inf) to that periodic orbit. In the case of the unstable manifold, it is the same argumentation but the periodic orbit is reached asymtoptically as time approaches -inf.
In a nutshell, the L1 point (as well as the other two collinear points) linearly behave like the product saddle × center × center. The two centers are associated to the in-plane and out-of-plane motions, and give rise to the families of planar and vertical periodic Lyapunov orbits. The Halo orbit is the special case in which the in-plane and out-of-plane frequencies match and hence a periodic orbit emerges. The saddle part is resposible for the presence of the two-dimensional stable and unstable manifolds that emanate from the periodic orbits. More info on periodic orbits, equilibrium points and manifolds in dynamical systems are also given in Nonlinear Dynamics and Chaos by Strogatz.
In the following figure, the stable (green) and unstable (red) invariant manifolds of the 177.86 day period halo orbit which is the target of Aditya-L1 are shown. The halo orbit is shown in black, the L1 point with the 'x' and the Earth is located at (0,0). The Sun is around 100 times the distance between L1 and Earth towards the left of the image so it is omitted for obvious reasons. Two branches are seen for both the stable and unstable manifolds, one towards the right of the halo orbit (the "Earth realm") and one towards the left (towards the "Sun realm").
So... to reach the periodic (halo) orbit for free, one can place the spacecraft on the stable manifold associated to the final (target) orbit. The one coming from the Earth realm is isolated in the following figure. As you can see it twists and turns around the Earth and some of its paths come very close to the vicinity of typical Earth-bound orbits (like the one onto which Aditya-L1 was placed before its injection burn).
Before I show how to jump onto one of these manifolds (or highways as I like to call them), I want to briefly describe how they are calculated.
When using a fully numerical method, a point of the periodic orbit (in this case the halo orbit) is perturbed in the local direction of the manifold and then integrated in time. To construct the whole manifold, this process is repeated for all the points in which the orbit is discretized. To do that, one first needs a linear approximation of the manifold in the neighborhood of the orbit.
If an orbit is periodic, then the eigenvectors of the variational matrix computed over one period of the orbit (the monodromy matrix), give the directions of the tangent to the invariant manifolds at the starting point. Hence the initial conditions on the orbit are perturbed in the direction of the tangent x(t1) ± ε vs(t1) (where vs(t1) is the stable eigenvector of the monodromy matrix). The small displacement ε perturbs x(t1) in the stable direction and it is small enough to preserve the local validity of the linear approximation, but also large enough to prevent from long integration times needed to propagate the manifold in time.
This is shown schematically in the figure below. More info on how to obtain the manifolds can be found in Topputo (2016) and in Dynamical Systems, the Three-Body Problem and Space Mission Design by Koon et al.
Method used to compute the invariant manifolds from Topputo (2016).
Ok enough with the theory, how does this process look like? Taking the same 177.86 day orbit as an example, I first disturbed along the eigenvactor and then propagated 20 points in time along the stable manifold. The resulting "tube" is the same as the one in the previous figure (but with fewer pathways).
BE PATIENT with the video, it takes around 15 seconds before anything happens. This is because the dynamics close to the periodic orbit are so slow that the trajectory along the manifold practically coincides with the trajectory along the periodic orbit. Once we reach the vicinity of the Earth however, things speed up a lot as you can see.
Stable manifold propagation from the halo orbit. The good part of the video starts ~15seconds in, before that the dynamics time-scales are too slow
Getting onto the "highway" from an earth-bound orbit
Ok if you still are following, then it is time to be rewarded for your patience and see how we can put all this together.
We saw how the halo orbits look like.
We saw how we can calculate them accurately.
We saw how the highways that lead to them look like.
The remaining question is how do we get onto the highway from an earth-bound orbit? Rather, to put it better, we don't want to only get onto the highway, but also pay no toll (or at least as little as we are allowed to).
To find the optimal solution, we can perform a trajectory optimization, where the initial condition is a point on an Earth-bound orbit and the target is a point on the manifold. Not just any point but the optimal manifold insertion point is searched for. To make sure that the jump from the Earth-bound orbit to the manifold is physical, we need to enforce a numerical constraint in the form of a boundary condition, involving the matching of the position on the manifold and the initial orbit.
The parameter that is minimized is the velocity difference between the two trajectories (earth orbit and manifold) at this intersection point, as this would need to be applied via the on-board propulsion system in form of a maneuver.
In this particular formulation of the problem, I assumed the apogee and perigee altitudes of the earth-bound orbit to be constant and equal to the values reported by the ISRO (122,000km and 256 km). The orientation of the orbit however was defined as a set of optimization variables (inclination, right ascencion, argument of perigee and true anomaly).
Different combinations of the "open" variables (Kepler elements of the Earth-bound orbit, ID of the manifold and time along the manifold from the PO) are varied at random in the next animation to make the geometric configuration clearer. The final solution shown in the animation is also one that satisfies the "crossing" constraint.
Variation of the optimization parameters.
And by solving this optimization problem...
Ta daaaaa! We have an optimal (DeltaV-minimal) solution for the trip to L1!
Only one impulsive maneuver is applied to embark onto the manifold and after that, the solution asymptotically reaches the periodic halo orbit. No more fuel needs to be expended (in this ideal scenario). For this configuration, we reach a trip duration of ~110 days until we reach the vicinity of the halo orbit, which is very similar to what is reported by the ISRO. The total DeltaV for the maneuver is around ~930 m/s according to this solution.
The solution was done in Matlab using fmincon for the optimization problem. ODE113 was used for all numerical integrations of the dynamic equations.
Note that the analysis is done in the circular restricted three-body framework, so several effects like eccentricity of the earth orbit around the Sun are omitted. Also, the maneuver is considered as impulsive and the analysis does not account for any maneuvers that have to be performed during the lifetime of the mission for stationkeeping, but this is a whole 'nother discussion.
Optimal solution to the halo orbit!
I hope that this guide was helpful and that it sparked some interest in the Aditya mission, dynamical systems and libration point orbits and their corresponding transfer trajectories! I am really excited to read more news about the next phases of the Aditya-L1 mission and wish the probe a successful remaining 108 days to its fascinating halo orbit!