A Single-Launch Deployment Strategy for Lunar Constellations

Abstract

Satellite constellations can provide communication and navigation services to support future lunar missions, and are attracting growing interest from both the scientific community and industry. The deployment of satellites in orbital planes that can have significantly different inclinations and right ascension of the ascending node requires dedicated launches and represents a non-trivial issue for lunar constellations, due to the complexity and low accessibility of launches to the Moon. In this work, a strategy to deploy multiple satellites in different orbital planes around the Moon in a single launch is examined. The launch vehicle moves along a conventional lunar escape trajectory, with parameters selected to take advantage of gravity-braking upon encountering the Moon. A maneuver at the periselenium allows the transfer of the spacecraft along a trajectory converging to the equilibrium region about the Earth–Moon libration point $L_1$, where the satellites are deployed. Providing a small $\Delta V$, each satellite is transferred into a low-energy trajectory with the desired inclination, right ascension of the ascending node, and periselenium radius. A final maneuver, if required, allows the adjustment of the semimajor axis and the eccentricity. The method is verified using numerical integration using high-fidelity orbit propagators. The results indicate that the deployment could be accomplished within one sidereal month with a modest $\Delta V$ budget.

1. Introduction

Lunar exploration has seen renewed interest in recent years, with the Artemis-1 mission validating the new NASA Space Launch System and Orion spacecraft, designed to take humans back to the Moon, and deploying the first CubeSat missions to the cislunar and translunar space [1,2,3,4]. The infrastructure that could support the forthcoming missions, either manned or autonomous, often relies on distributed satellite systems, such as constellations, that can provide navigation, telecommunication, and data relay services to both in-orbit and on-ground facilities [5,6,7,8,9,10].

This deployment represents one of the most critical phases in the implementation of a satellite constellation, and its complexity is due to the need to distribute the satellites on different orbital planes, a task that is usually accomplished with dedicated launches. For constellations operating in low-Earth orbit (LEO), this issue is mitigated by the large number of launch opportunities available on the market. Unfortunately, this is not an option yet for lunar constellations, since the number of launches to the Moon is still limited to no more than a few per year. It should be observed that, when dealing with a small satellite constellation, limitations on the propellant (mass) budget and on the propulsion systems characterizing this class of spacecraft further complicate this issue, since traditional plane change maneuvers are not possible [11,12].

A strategy is proposed here to deploy multiple satellites in different orbital planes based on a single-launch opportunity equivalent to that of the Artemis-1 mission. This result is achieved by taking advantage of low-energy trajectories associated with the Earth–Moon libration point $L_1$ [13].

Libration points and their associated trajectories have been studied for decades [13,14,15,16,17,18,19,20,21,22]. In 1978, the International Sun–Earth Explorer (ISEE)-3 was the first spacecraft to be deployed to a libration point orbit (LPO) [23], followed in 1991 by the Hiten mission [24], the first succeeding in an Earth–Moon transfer using a low-energy trajectory. In the late 1960s, Farquhar was a pioneer in proposing constellations of satellites located in the neighborhood of Earth–Moon collinear libration points $L_1$ and $L_2$ [25], devoted to communication operations. A few decades later, Lee et al. extended this design to Earth communication constellations [26] and, more recently, the use of LPO constellations has been investigated for lunar navigation services [27,28,29].

It should be noted that, when compared to LEO ones, LPO constellations allow the reduction of the fueling effort for station-keeping and the eclipse period but, by virtue of their larger distance from the ground, they are characterized by larger communication delay times and higher power consumption. Therefore, LEO constellations could be preferred for some applications [30,31].

In 2004, Chase et al. observed that a satellite can be transferred from a Lissajous orbit around the Earth–Moon libration point $L_1$ to an Earth orbit whose inclination and right ascension of the ascending node (RAAN) depend only on the epoch at the beginning of the transfer [32]. Based on this fact, they proposed a constellation deployment strategy in which all the satellites are transferred from the same launch trajectory to the same Lissajous orbit where station-keeping is performed independently by each satellite waiting for the right epoch to transfer back to an Earth orbit with the desired inclination and RAAN values. Operating in the dynamic framework of the circular restricted 3-body problem (CR3BP), Nadoushan and Novinzadeh proposed a similar solution, in which the satellites are first transferred to the same halo orbit about the Earth–Moon $L_1$ and are transferred back to LEO at different times, taking advantage of nodal precession to target the arrival orbital plane [33].

The single-launch deployment strategy proposed here is developed based on the characterization of low-energy trajectories which cross the phase space surrounding the libration point $L_1$, named the equilibrium region [34,35]. Operating in the dynamic framework of the elliptic restricted four-body problem (ER4BP), which is significantly more accurate than the CR3BP to model the Sun–Earth–Moon system, it is proved that a satellite can be transferred from a low-energy trajectory to an orbit with desired semimajor axis, eccentricity, inclination, and RAAN by providing a single $\Delta V$ when crossing the equilibrium region. If compared to other solutions available in the literature [32,33], this novel deployment strategy (i) does not require the station-keeping phase around any LPO, and (ii) allows the selection of the target orbital plane by a single maneuver.

These concepts are applied to design the transfer of multiple satellites departing from the same trajectory into different lunar orbits with desired semimajor axis, eccentricity, inclination, and RAAN. The mission consists of two phases. During the first one, the satellites are stowed in the launch vehicle, moving along a lunar escape trajectory with parameters selected such that a powered gravity-braking at the periselenium allows the transfer of the vehicle and its payload along a low-energy trajectory directed towards the equilibrium region. The second phase starts when each satellite maneuvers independently to be injected into a low-energy trajectory approaching the Moon with the desired inclination, RAAN, and periselenium radius. A final maneuver to adjust the semimajor axis concludes the deployment. The method is evaluated using numerical tools simulating the deployment of a lunar constellation of elliptic inclined orbits providing global coverage [30].

The paper is organized as follows. Section 2 collects the theoretical background on the CR3BP, the equilibrium region, the topological characterization of low-energy trajectories, and an extension of these solutions to the Sun–Earth–Moon system. In Section 3, the deployment strategy is developed, and in Section 4, it is verified through numerical simulations, and the results are compared to other solutions available in the literature.

2. Low-Energy Captures in the Sun–Earth–Moon System

2.1. The Circular Restricted 3-Body Problem

The CR3BP represents the most common dynamic framework to study the motion of a spacecraft in the Earth–Moon system. This model is based on the hypothesis that the mass m of the spacecraft is negligible compared to that of the primaries, the Earth ($m_1$) and the Moon ($m_2$), and that in the inertial space, the relative trajectory of $m_2$ about $m_1$ is circular [36].

The dynamic equations of motion for the CR3BP are expressed in a reference frame ${\mathcal{F}}_r=[\widehat{x},\widehat{y},\widehat{z}]$, centered at the center of mass of the system O, with $\widehat{x}$ pointing from $m_1$ to $m_2$, $\widehat{z}$ orthogonal the orbital plane of the primaries and $\widehat{y}$ completing the rectangular reference frame [37]

$$\left\{\begin{array}{c}\dot{x}=v_x\hfill \\ \dot{y}=v_y\hfill \\ \dot{z}=v_z\hfill \\ {\dot{v}}_x=2v_y+\frac{\partial U}{\partial x}\hfill \\ {\dot{v}}_y=-2v_x+\frac{\partial U}{\partial y}\hfill \\ {\dot{v}}_z=\frac{\partial U}{\partial z}\hfill \end{array}\right.$$

(1)

with

$$U=\frac{1}{2}\left(x^2+y^2\right)+\frac{1-\mu }{r_1}+\frac{\mu }{r_2}$$

(2)

$$\mu =\frac{m_2}{m_1+m_2}$$

(3)

where $r_i$ indicates the distance between the spacecraft and the i-th primary. It is worth recalling that Equation (1) are expressed in terms of units of distance $DU=a_M$ and time $TU=1/{\omega }_M$, where $a_M$ and ${\omega }_M$ are, respectively, the mean distance and the orbital angular speed of the Moon with respect to the Earth.

The CR3BP admits five equilibrium points $L_i$, also named the Lagrange or libration points [38], and one constant integral of motion, named the Jacobi constant

$$C=-\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)+2U$$

(4)

In fact, the Jacobi constant is proportional to the opposite of the total energy $\mathcal{E}$ of the spacecraft; therefore, C decreases as $\mathcal{E}$ increases.

Trajectories that can evolve in the neighborhood of both the primaries are characterized by $C<C_1$, where $C_1$ is the Jacobi constant calculated in $L_1={[x_{L_1},0,0,0,0,0]}^T$, the libration point lying on $\widehat{x}$ between $m_1$ and $m_2$ [13,39]. The deployment strategy proposed in this work is based on the use of those low-energy trajectories, which verify the condition $C=C_1-2h$, where $h>0$ is an arbitrarily small constant named the energy level. The characterization of these solutions is discussed in the following section.

2.2. Characterization of Low-Energy Trajectories Crossing the Equilibrium Region

It is known that the ultimate behavior of low-energy trajectories characterized by a low energy level h can be predicted based on their state representation inside the equilibrium region, i.e., the phase space surrounding the libration point $L_1$ [34,40]. The rigorous definition of the equilibrium region and the representation are developed hereafter using the Hamiltonian formalism and the results are later applied to design trajectories that, after crossing the equilibrium region, reach the Moon with desired semimajor axis, eccentricity, inclination, and RAAN.

The Hamilton’s equations are given by

$$\left\{\begin{array}{c}\dot{\mathbf{q}}=\frac{\partial H}{\partial \mathbf{p}}\hfill \\ \dot{\mathbf{p}}=-\frac{\partial H}{\partial \mathbf{q}}\hfill \end{array}\right.$$

(5)

where $\mathbf{q}={[x-x_{L_1},y,z]}^T$ and $\mathbf{p}={[v_x-y,v_y+x-x_{L_1},v_z]}^T$ are the position and conjugate momenta. The Hamiltonian function for the CR3BP is by definition $H=-C/2$. Aiming at investigating the dynamics in the neighborhood of $L_1$, the following expression for H can be conveniently derived from quadratic expansion of Equation (4) about $L_1$ [35]

$$H(\mathbf{q},\mathbf{p})=\frac{1}{2}(p_1^2+p_2^2+p_3^2+2p_1{q}_2-2p_2{q}_1)-\frac{K}{2}(2q_1^2-q_2^2-q_3^2)$$

(6)

with $K=\frac{1-\mu }{|x_{L_1}{+\mu |}^2}+\frac{\mu }{|x_{L_1}{+\mu -1|}^2}$.

The linear system associated with Equations (5) and (6) is characterized by a saddle–center–center type of equilibrium, having one couple of real eigenvalues ($\pm \rho$) and two couples of complex conjugate ones ($\pm j{\lambda }_1$, $\pm j{\lambda }_2$). Therefore, according to Morse’s lemma, the Hamiltonian function can be represented as the sum of three local integrals of motion, each dependent on a different pair of state variables and associated with one of the three above-mentioned eigenspaces. This new form of H can be conveniently derived by applying a canonical transformation ${[\mathbf{x},\mathbf{y}]}^T={\mathbf{T}}_N{[\mathbf{q},\mathbf{p}]}^T$, originally proposed by Siegel and Moser [41], producing

$$H(\mathbf{x},\mathbf{y})=\rho x_1{y}_1+\frac{{\lambda }_1}{2}(x_2^2+y_2^2)+\frac{{\lambda }_2}{2}(x_3^2+y_3^2)=h$$

(7)

with $x_1,y_1\in \mathbb{R}$, $x_2,x_3,y_2,y_3\in \mathbb{C}$ and $h>0$ is the arbitrarily small energy level. Equation (7) is the mathematical representation of the equilibrium region, the phase space in the neighborhood of $L_1$.

Based on a theorem by Moser [42], Conley proved that if a low-energy trajectory is inside the equilibrium region at a given time $t_0$, then its long-term behavior is fully characterized by its topological location in the equilibrium region and provided the following classification [34]

• $x_1\left(t_0\right)y_1\left(t_0\right)<0$, transit trajectories which evolve alternately around one of the two primaries, crossing the equilibrium region multiple times
• $x_1\left(t_0\right)y_1\left(t_0\right)>0$, bouncing trajectories which never cross the equilibrium region, thus evolve only around one of the two primaries
• $x_1\left(t_0\right)\to 0$ or $y_1\left(t_0\right)\to 0$, capture trajectories that cross only once the equilibrium region and then evolve around one of the primaries indefinitely in time
• $x_1\left(t_0\right)=y_1\left(t_0\right)=0$, quasi-periodic orbits which never depart from the equilibrium region.

Within the scope of this research, trajectories transiting from the Earth to the Moon (i.e., transits or captures) are of interest. In particular, low-energy captures are attractive because they do not require any further maneuvering to keep the satellites in the neighborhood of the Moon, thus ensuring relaxed operation times to perform any eventual corrective maneuver.

For this class of trajectories, the osculating orbital elements at capture can be characterized, as for the long-term behavior, by their state representation at time $t_0$ when the trajectory is inside the equilibrium region [40]. Equation (7) can be rearranged as follows [43]

$$\frac{{\lambda }_1}{2}(x_2^2+y_2^2)+\frac{{\lambda }_2}{2}(x_3^2+y_3^2)=h(1-\epsilon )$$

(8)

with $\epsilon >0$ arbitrarily small. The two constant terms of Equation (8) are named energy fractions and indicated as $h_1=\frac{{\lambda }_1}{2}(x_2^2+y_2^2)$ and $h_2=\frac{{\lambda }_2}{2}(x_3^2+y_3^2)$. The value of the energy fractions can be related to those of some osculating orbital elements. For the sake of clarity, the following expressions for the properties of the canonical transformation ${\mathbf{T}}_N$ shall be recalled:

$$\left\{\begin{array}{c}{y}_2=-j{\overline{x}}_2(x,y,v_x,v_y)\hfill \\ y_3=-j{\overline{x}}_3(z,v_z)\hfill \end{array}\right.$$

(9)

From the conservation of the angular momentum, the following expression for the inclination i can be derived [44]

$$cosi=\frac{1}{\sqrt{1+\frac{h_2}{h_1}}}$$

(10)

Therefore, a target inclination at capture can be fixed by selecting any couple ($x_2\left(t_0\right),x_3\left(t_0\right)$) verifying Equation (10). By virtue of Equation (9), Equation (10) represents a constraint for the six state variables in the position space. Moreover, if h is not fixed (i.e., it can be varied within an admissible range) for a given set of in-plane variables $(x,y,v_x,v_y)$, the target inclination can be fixed by selecting only the out-of-plane variables $(z,v_z)$.

Another relationship between the energy fractions can be derived based on Tisserand’s parameter

$$\sqrt{a(1-e^2)}cosi\propto const.$$

(11)

where a and e are the semimajor axes and the eccentricity of the capture orbit. Introducing Equation (10) into (11) leads to

$$a(1-e^2)\propto \frac{h_1+h_2}{h_1}=\frac{h}{h_1}$$

(12)

Equations (10) and (12), which were verified using numerical analysis in a previous study by the author [44], indicate that fixing the value of the inclination and a combination of a and e (or equivalently, the desired pericenter distance $r_p$ at capture) the energy fractions $h_1$ and $h_2$-and therefore the values of $x_2\left(t_0\right)$ and $x_3\left(t_0\right)$-are determined.

In Section 3, the above-mentioned constraints are rearranged in terms of position and velocity in ${\mathcal{F}}_r$ and implemented to design the deployment strategy.

2.3. Extension to the Sun–Earth–Moon System

The design of low-energy captures by the Moon, as well as the deployment strategy derived from it, cannot ignore two conditions that are overlooked in the CR3BP model: neither the eccentricity of the Earth and Moon orbits nor the gravitational field of the Sun are negligible in the real environment [45,46].

The topological characterization presented in Section 2.2 is here extended to the more accurate dynamic framework of the elliptic restricted four-body problem (ER4BP) that is used to model the Sun–Earth–Moon system. The dynamic equations of motion are developed in a reference frame ${\mathcal{F}}_r^{\left(e\right)}=[\widehat{x},\widehat{y},\widehat{z}]$ centered in the center of mass of the system, with $\widehat{x}$ pointing from the Sun ($m_0$) to the center of mass of the Earth–Moon system, $\widehat{z}$ orthogonal to the ecliptic plane and $\widehat{y}$ completing the rectangular frame [47]. Denoting by $r_E$ and ${\theta }_E$ the distance and the true anomaly of the Earth–Moon center of mass O with respect to the Sun, the equations can be expressed as follows [48]

$$\left\{\begin{array}{c}{x}^{\prime }=v_x\hfill \\ y^{\prime }=v_y\hfill \\ z^{\prime }=v_z\hfill \\ v_x^{\prime }=2v_y+\tau \left(\frac{\partial u}{\partial x}+x\right)\hfill \\ v_y^{\prime }=-2v_x+\tau \left(\frac{\partial u}{\partial y}+y\right)\hfill \\ v_z^{\prime }=\tau \left(\frac{\partial u}{\partial z}+z\right)-z\hfill \end{array}\right.$$

(13)

where $\tau ={\left(1+e_{M}cos{\theta }_E\right)}^{-1}$, ${\mu }_i=\frac{m_i}{m_0+m_1+m_2}$, $m_0$ is the mass of the Sun, $u={\displaystyle \sum _{i=1}^3}\frac{{\mu }_i}{r_i}$, and the prime ${}^{\prime }$ indicates the derivative with respect to ${\theta }_E$. Equation (13) are expressed in terms of the units of distance $D{U}^{\left(e\right)}=r_E\left({\theta }_E\right)$ and time $T{U}^{\left(e\right)}=\sqrt{\frac{r_E^3\left({\theta }_E\right)}{G\left(m_0+m_1+m_2\right)}}$ where G is the gravitational constant.

Because of its dependence on time, the ER4BP does not admit equilibrium points. Nevertheless, their instantaneous dynamic substitutes can be computed at any given time considering the corresponding geometric configuration of the primaries [48]. Operating as in Section 2.2, the set of Equation (13) can be linearized about the instantaneous libration point $L_1$ of the Earth–Moon system and expressed using the Hamiltonian formalism using the following Hamiltonian function [49]

$$\begin{array}{c}{H}^{\left(e\right)}=\frac{1}{2}\left(p_1^2+p_2^2+p_3^2\right)+p_1{q}_2-p_2{q}_1-\frac{K}{2}\left(2q_1^2-q_2^2-q_3^2\right)+\\ +e_M{c}_{1}cos{\theta }_E+\frac{{\mu }_0}{r_0}{c}_2\end{array}$$

(14)

where $r_0$ denotes the mean distance between the Sun and the center of mass of the Earth–Moon system and the expressions for the coefficients $c_i$ are reported in Appendix A. It can be observed that the expressions of H and $H^{\left(e\right)}$ are equivalent, except for the last two terms of Equation (14), which represent the effects of the eccentric motion of the primaries ($e_M$) and the solar gravitational perturbation ($\frac{{\mu }_0}{r_0}$). Examining in detail the two coefficients, it is easy to verify that for the Earth–Moon system $e_M<<1$ and in the neighborhood of $L_1$ also $\frac{{\mu }_0}{r_0}<<1$, so they can be regarded as small perturbations acting onto the CR3BP.

Conley and Easton proved that the basic topological properties of the CR3BP in the neighborhood of $L_1$ are preserved in the presence of small perturbations [50]. In fact, a canonical transformation $\mathcal{T}:(\mathbf{q},\mathbf{p})\to (\mathbf{Q},\mathbf{P})$, developed in Appendix A, can be introduced to rearrange Equation (14) to a form equivalent to H plus negligible higher order terms in the perturbations [49]

$${\tilde{H}}^{\left(e\right)}(\mathbf{Q},\mathbf{P})=\frac{1}{2}\left(P_1^2+P_2^2+P_3^2\right)+P_1{Q}_2-P_2{Q}_1-\frac{K}{2}\left(2Q_1^2-Q_2^2-Q_3^2\right)+o\left(e_M,\frac{{\mu }_0}{r_0}\right)$$

(15)

Because ${\tilde{H}}^{\left(e\right)}(\mathbf{Q},\mathbf{P})$ is equivalent to $H(\mathbf{q},\mathbf{p})$, all the results derived from the application of the canonical transformation ${\mathbf{T}}_N$, and in particular the characterization of low-energy captures, are still valid [51]. Capture orbits with the desired orbital elements can be designed in the $(\mathbf{x},\mathbf{y})$ coordinates and then converted to position and velocity coordinates in ${\mathcal{F}}_r^{\left(e\right)}$, applying in order the inverse of the two canonical transformations ${[\mathbf{Q},\mathbf{P}]}^T={\mathbf{T}}_N^{-1}[\mathbf{x},\mathbf{y}]$ and ${\mathcal{T}}^{-1}:(\mathbf{Q},\mathbf{P})\to (\mathbf{q},\mathbf{p})$.

3. Design of the Deployment Strategy

The constellation deployment strategy designed in this section is developed to allow the transfer of satellites from the same launch trajectory to different orbital planes about the Moon. The strategy is organized into two phases. First, the launch vehicle is transferred to a low-energy trajectory approaching the equilibrium region. Here the satellites are released from the launcher and, independently, are transferred to low-energy lunar captures with prescribed orbital elements. The design of the two phases is presented separately and the design parameters discussed hereafter are summarized in

Table 1. Summary of free and constrained parameters during the different phases of the deployment.

Phase:	Free	Constrained
Lunar flyby	${\Omega }_d$ ; ${\omega }_d$ ; $i_d$	$a_d$ ; $e_d$ ; ${\theta }_d$
Injection into equilibrium region	${\mathbf{V}}_{p,a}$	${\mathbf{r}}_{p,a}$
Targeting capture parameters	$v_{x,er}$ ; $v_{y,er}$ ; $v_{z,er}$	$x_{er}$ ; $y_{er}$ ; $z_{er}$
Corrections	${\mathbf{V}}_{p,c}$	${\mathbf{r}}_{p,c}$

and

Table 2. Summary of control and target parameters during the different phases of the deployment.

Phase:	Control	Target
Lunar flyby	${\Omega }_d$ ; ${\omega }_d$ ; $i_d$	$r_{p,a}$ ; $y_a$ ; $i_a$
Injection into equilibrium region	$\Delta V_{x,a}$ ; $\Delta V_{y,a}$ ; $\Delta V_{z,a}$	$x_{er}$ (${\Omega }_c$) ; $y_{er}$ ; $h_{er}$
Targeting capture parameters	$v_{x,er}$ ; $v_{y,er}$ ; $v_{z,er}$	$x_1$ ; $i_c$ ; $r_{p,c}$
Corrections	$\Delta V_{x,c}$ ; $\Delta V_{y,c}$ ; $\Delta V_{z,c}$	$e_c$

. For the sake of clarity, the algorithm implementing the design of the deployment strategy, used for the numerical study in Section 4, is reported at the end of this section (see Algorithm 1).

Algorithm 1 Deployment strategy design.

(1) Determine the orbital elements ${\Omega }_d$, ${\omega }_d$, $i_d$ at the initial epoch $t_0$ by solving a two-point boundary value problem (TPBVP) using the following sets of boundary conditions $a_d$, $e_d$, ${\theta }_d$ at $t_0$ $r_{p,a}$, $y_a$, $i_a$ at the time of periselenium passage ($t_a$) (2) Using the orbital elements determined in Step 1, propagate the trajectory of the launch vehicle to the periselenium and compute the $\Delta V_a$ to be provided for reducing the energy level to the desired value $h_{er}$ (3) After Step 2, the satellites are deployed from the launch vehicle and maneuvered independently; for each satellite select a guess value for $x_{er}$ within a suitable interval propagate the trajectory from the periselenium until it reaches the coordinate $x_{er}$ and compute the remaining state variables $y_{er}$, $z_{er}$, $v_x{}_{er}$, $v_y{}_{er}$ and $v_z{}_{er}$ solve Equation (17) to determine the $\Delta V_{er}$ to be provided to establish the desired $i_c$ and $r_{p,c}$ at capture propagate the trajectory to the periselenium if the difference between the actual value of RAAN and its target value ${\Omega }_c$ is larger than a prescribed error then repeat from Step 3.1 with a different value of $x_{cr}$ else proceed to Step 4 (4) Compute the $\Delta V_c$ eventually required to correct $e_c$

3.1. Injection to the Low-Energy Trajectory

In the first phase of the deployment strategy the launch vehicle, in which the satellites to be deployed are stowed, is transferred from its departure trajectory to the equilibrium region. The initial epoch $t_0$ of the operations is fixed here to the end of the translunar injection (TLI) maneuver. The position of the spacecraft at $t_0$ depends on the performance of the launch vehicle and on the launch profile [52]; therefore, it is reasonable to assume that the semimajor axis $a_d$, eccentricity $e_d$ and true anomaly ${\theta }_d$ at $t_0$ are given. By contrast, the launch date and the value of the other orbital elements at $t_0$ (the RAAN ${\Omega }_d$, the argument of perigee ${\omega }_d$ and the inclination $i_d$) can be selected in the design phase.

The injection into the equilibrium region is performed by taking advantage of lunar gravity-braking; therefore, the orbital elements ${\Omega }_d$, ${\omega }_d$, and $i_d$ are selected such that the departure trajectory reaches the Moon at a minimum distance $r_{p,a}$, inclination $i_a$ and negative value for the position coordinate $y_a$. In particular, the condition $y_a<0$ is necessary for the spacecraft energy to decrease after the flyby (i.e., $y_a>0$ would result in a gravity-assist) while the minimization of $r_{p,a}$ and $i_a$ is pursued to maximize the gravity-braking effect [53].

Once reached, the pericenter ${\mathbf{r}}_{p,a}$, a $\Delta V_{p,a}$ is provided along the direction of the velocity to reduce the energy level to a target value $h_{er}$ (see Section 3.2), injecting the launcher and its payload into a low-energy trajectory heading towards the equilibrium region. Finally, the satellites are released from the launcher and, from now on, they are maneuvered independently.

3.2. Deployment to the Final Orbit

The characterization presented in Section 2.2 indicates that, once the satellites have reached the equilibrium region, a small change in the velocity $\Delta V_{er}$ can be implemented to transfer each one of them to a low-energy capture with desired inclination $i_c$ and periselenium radius $r_{p,c}$. The components $\Delta V_{x,er}$, $\Delta V_{y,er}$ and $\Delta V_{z,er}$ shall be selected to satisfy the following conditions

$$\left\{\begin{array}{c}{x}_1=\epsilon \hfill \\ \frac{h_2}{h_1}=\frac{1}{{cos}^2{i}_c}-1\hfill \\ {tan}^{-1}\left(\frac{Im\left(x_2\right)}{1+Re\left(x_2\right)}\right)=\frac{2r_{p,c}}{r_{SOI}}\hfill \\ y_2=-j{\overline{x}}_2\hfill \\ y_3=-j{\overline{x}}_3\hfill \end{array}\right.$$

(16)

with $\epsilon >0$ arbitrarily small and $r_{SOI}=a_M{\left(\frac{m_2}{m_1}\right)}^{\frac{2}{5}}$ is the radius of the sphere of influence of the Moon. Equation (16) represents five constraints in the six transformed Hamiltonian variables ($\mathbf{x},\mathbf{y}$). The first three represent the capture condition and the inclination and periselenium radius at capture, while the last two result from the properties of the canonical transformation ${\mathbf{T}}_N$ given by Equation (9).

For the sake of clarity, Equation (16) can be transformed by making explicit the dependence on the position space variables

$$\left\{\begin{array}{c}\frac{\left(1+2{\gamma }^2\right)(x-x_{L_1})+\left(1-{\gamma }^2\right)y+\alpha v_x-\alpha \sigma v_y}{2\left[\left(1-{\gamma }^2\right){\sigma }^2+{\rho }^2\right]}=\epsilon \hfill \\ v_z=\sqrt{\frac{2h_1{tan}^2{i}_c}{{\lambda }_2}}\hfill \\ \frac{t_{5,1}x+t_{5,5}{v}_y}{1+t_{2,2}y+t_{2,4}{v}_x}=tan\left(\frac{2r_{p,c}}{r_{SOI}}\right)\hfill \\ \\ h_1=2{\lambda }_1\left(t_{2,2}y+t_{2,4}{v}_x\right)\left(t_{5,1}x+t_{5,5}{v}_y\right)\hfill \\ h_2=2{\lambda }_2\left(t_{3,3}z\right)\left(t_{6,6}{v}_z\right)=h-h_1\hfill \end{array}\right.$$

(17)

where ${\gamma }^2=\left(-\frac{\mu }{|x_{L_1}{+\mu -1|}^3}-\frac{1-\mu }{|x_{L_1}{+\mu |}^3}\right)$, $\sigma =-2\rho /\left({\rho }^2+{\gamma }^2-1\right)$, $\tau =\frac{2({\lambda }_2^2-1)}{{\lambda }_2^4-3{\lambda }_2^2+2}$ and $t_{5,1}=0.1764$, $t_{2,2}=1.0387$, $t_{2,4}=0.1409$, $t_{5,5}=0.5053$, $t_{3,3}=1.0651$, $t_{6,6}=0.4694$ are the $t_{i,j}$ elements of the ${\mathbf{T}}_N$ matrix. Once a suitable combination of $\epsilon$, $i_c$,$r_{p,c}$ and h has been fixed, Equation (17) can be used to select five out of the six position and velocity coordinates.

The last degree of freedom allows the selection of the RAAN at capture ${\Omega }_c$. The value of ${\Omega }_c$ is strictly related to the time $t^{★}$ at which the low-energy capture trajectory crosses the Moon equatorial plane for the first time, which can be delayed or anticipated by changing the time at which the $\Delta V_{er}$ is provided. Recalling now that all the satellites enter the equilibrium region at the same time, a simple way to select different values of ${\Omega }_c$ for different satellites is to provide the $\Delta V_{er}$ at different coordinates x. In any case, because a closed-form solution equivalent to Equation (16) to relate the difference in x to that in ${\Omega }_c$ is not available, the problem must be solved iteratively. Nevertheless, this does not represent an actual issue, since numerical solvers can find a solution in a few iterations.

As mentioned above, providing $\Delta V_{er}$ changes the energy level from the value $h_{er}$, before the maneuver, to $h_c$. Therefore, another iterative process can be implemented to minimize the difference between the energy level $h_{er}$ and $h_c$. Even though this step is not strictly required, it is useful to minimize the total $\Delta V$.

The low-energy capture trajectory obtained from the process above is inherently chaotic; consequently, the values of its orbital elements change in time affecting the performance of the constellation. To overcome this issue, if required, a final corrective maneuver can be performed at the periselenium, providing a $\Delta V_c$ to reduce the semimajor axis.

4. Analysis of the Constellation Deployment

The deployment strategy designed in Section 3 is here examined using numerical methods. Initial conditions equivalent to those of the Artemis-1 mission are considered [54], corresponding to $a_d$ = 202,700 km, $e_d=0.968$ and ${\theta }_d=29$ deg. The satellites are assumed to be deployed in elliptical inclined lunar orbits, whose orbital elements are reported in

Table 3. Orbital elements for the satellites of the constellation.

Element	Sat-1	Sat-2	Sat-3	Sat-4	Sat-5
$a_c$	$6541.4$ km	$9873.0$ km	$7500.0$ km	$7494.7$ km	$6541.4$ km
$e_c$	$0.6$	$0.2$	$0.05$	$0.05$	$0.6$
$i_c$	$56.2$ deg	45 deg	40 deg	40 deg	$56.2$ deg
${\Omega }_c$	50 deg	30 deg	50 deg	10 deg	140 deg

Expressed in the Moon-centered MJ2000Eq reference frame.

, which are suitable for constellations that shall provide polar and global coverage of the Moon [30]. The analysis is performed using the General Mission Analysis Tool R2020a developed by NASA, in which the algorithm presented in Section 3 was implemented. The propagator includes the Sun and the Earth as point masses while the gravitational field of the Moon is described by the LP-165 model. The motion of the celestial bodies is modeled using the NASA JPL Development Ephemerides (DE405).

4.1. Application of the Deployment Strategy

Phase 1, consisting of the transfer to the equilibrium region, is common to all the satellites. Steps 1–2 of the algorithm reported in Section 3 are implemented in GMAT after selecting the following target parameters for the lunar flyby $r_{p,a}=2037.1\pm 50$ km, $y_a<0$ and $i_a=30\pm 2$ deg, while the initial epoch is fixed to 17 October 2023, 00:00. The value of $r_{p,a}$, corresponding to an altitude of $300\pm 50$ km from the lunar surface, is selected because it represents an effective trade-off between the gravity-braking effect and the harmonic perturbations of the lunar gravitational field [55]. The analysis provides the orbital elements of the departure trajectory corresponding to the above-mentioned target parameters, both reported in

Table 4. Orbital elements of the departure trajectory.

Element	Value at Departure ${}^1$	Value at Arrival ${}^2$
a	202,700 km	$-5694$ km
e	$0.968$	$1.369$
i	$33.72$ deg	$28.50$ deg
$\Omega$	$42.97$ deg	$355.04$ deg
$\omega$	$62.67$ deg	$173.02$ deg
$\theta$	$29.00$ deg	0 deg

¹ Expressed in the Earth-centered MJ2000Eq reference frame. ² Expressed in the Moon-centered MJ2000Eq reference frame.

.

At the periselenium, the first maneuver performed by the launch vehicle concludes Phase 1, transferring it to a trajectory heading towards the equilibrium region. Based on previous analyses on capture in the Earth–Moon system [44], for a low-energy capture to incline the range reported in Table 3, an energy value $h_{er}=2\times {10}^{-2}$ is required. To achieve it, a $\Delta V_p=228.4$ m/s is provided in the direction opposite to the orbital velocity. Figure 1 shows the trajectory of Phase 1 represented in the rotating frame ${\mathcal{F}}_r$ with the origin shifted to $x_{L_1}$. The injection into the equilibrium region is achieved on 1 November 2023, at 12:39, about 15.5 days after the departure epoch.

The iterative process reported in Step 3 of the algorithm is implemented to determine, independently and for each satellite, the coordinate $x_{er}$ at which the capture maneuver shall be performed to achieve the desired ${\Omega }_c$ at capture. The following parameters are considered:

• the coordinate $x_{er}$ is selected within the interval $[200,1200]$ km, to ensure that the satellites are maneuvered close enough to $L_1$, where the linear model provides an accurate description of the dynamics;
• after propagation to $x_{er}$, the velocity variation $\Delta V_c$ required to achieve $i_c$ and $r_{p,c}$ reported in Table 3 is calculated, assuming a tolerance of, respectively, 2 deg, 10 km and $\epsilon ={10}^{-10}$;
• the capture trajectory is propagated to the periselenium and the value of the RAAN is verified, the process is repeated until the RAAN converges to ${\Omega }_c\pm 2$ deg;
• a final iteration allows refining the values of $x_{er}$ and $\Delta V_c$ to achieve the target orbital elements within a tolerance of $0.1$ deg and 1 km.

The process described above was applied to the five satellites obtaining the state variables reported in

Table 5. State variables in the equilibrium region before and after the maneuver.

Element	Sat-1	Sat-2	Sat-3	Sat-4	Sat-5
$x_{er}$	909 km	859 km	908 km	908 km	209 km
$y_{er}$	2719 km	2706 km	2719 km	2719 km	2542 km
$z_{er}$	3073 km	3070 km	3072 km	3072 km	3035 km
$v_{x,er}$	$-239.0$ m/s	$-239.2$ m/s	$-238.6$ m/s	$-238.6$ m/s	$-239.4$ m/s
$v_{y,er}$	$-62.2$ m/s	$-61.9$ m/s	$-62.1$ m/s	$-62.1$ m/s	$-58.8$ m/s
$v_{z,er}$	$-12.7$ m/s	$-12.7$ m/s	$-12.7$ m/s	$-12.7$ m/s	$-13.0$ m/s
$\Delta v_{x,er}$	$1.8$ m/s	$263.3$ m/s	$308.3$ m/s	$239.4$ m/s	$144.2$ m/s
$\Delta v_{y,er}$	$252.9$ m/s	$84.8$ m/s	$97.8$ m/s	70 m/s	$482.8$ m/s
$\Delta v_{z,er}$	$-68.2$ m/s	$-87.1$ m/s	$-72.5$ m/s	$-22.6$ m/s	$-143.0$ m/s

.

A final corrective maneuver $\Delta V_c$ at the periselenium also allows the fixing of the value of $a_c$ and $e_c$. The values of all the $\Delta V$, provided in the direction opposite to the orbital velocity, are reported in

Table 6. Required $\Delta V$ by each satellite on the different mission phases.

Element	Sat-1	Sat-2	Sat-3	Sat-4	Sat-5
$\Delta V_a$ [m/s]	228	228	228	228	228
$\Delta V_{er}$ [m/s]	262	290	332	251	524
$\Delta V_c$ [m/s]	158	163	232	236	192
$\Delta V_{TOT}$ [m/s]	648	681	792	715	944

, and the trajectories resulting from the deployment strategy are shown in Figure 2, Figure 3 and Figure 4, represented in the ${\mathcal{F}}_r$ reference frame with origin shifted to $x_{L_1}$. The deployment of the last satellite will be concluded on 11 November 2023, at 23:56, approximately 26 days after the initial epoch.

4.2. Discussion on the Results

The analysis proves the suitability of the single-launch deployment strategy to reach the desired targets within one sidereal month and a total $\Delta V$ less than 1 km/s for each satellite. Recalling that the injection into the equilibrium region is performed by the launch vehicle, the $\Delta V$ that shall be provided by the satellite never exceeds 730 m/s, a value that can be provided by most small satellites. For instance, considering a small satellite of wet mass at launch of 100 kg, a monopropellant propulsion system providing a thrust $T=20$ N with a specific impulse $I_{sp}=230$ s can easily perform the maneuvers requiring a total propellant mass $m_p=32.4$ kg, which is shown hereafter to be lower than that required by other innovative deployment strategies.

It is worth now comparing the deployment strategy to other innovative ones available in the literature. Mahdisoozani et al. and Koblick and Choi proposed interesting solutions to deploy satellites from the same launch trajectory to orbits with different RAAN and inclination, respectively [56,57]. It should be noted that the solution by Mahdisoozani et al. addresses Earth constellations. Compared to these strategies, the one developed here offers the advantages of allowing the selection of both parameters in a single maneuver requiring 50 % or less of $\Delta V$ for this task, leading to relevant savings in the propellant mass budget.

Nadoushan et al. proposed the use of the CR3BP and nodal precession due to gravitational perturbations to deploy multiple satellites in different orbital planes from the same launch trajectory [33]. Additionally, this method is designed for Earth constellations, but can be extended to lunar ones. Compared to it, the strategy proposed in this work is more rapid, because it does not require the nodal precession to select the RAAN, and can be obtained at slightly lower $\Delta V$, because the CR3BP solutions developed by Nadoushan et al. are constructed from libration point halo orbits that are characterized by a higher energy level with respect to that of low-energy captures ($h={10}^{-2}$) considered here.

A different solution is proposed by Pontani et al., who suggest the use of low-thrust propulsion and nonlinear orbit control to deploy a satellite constellation around Mars [58], with consequent substantial savings in the required mass of propellant. Even though it is not trivial to predict the $\Delta V$ required for the deployment in the hypothesis of extending the single-launch deployment to Martian constellations, which is an actual possibility [59], it is reasonable to suppose that it would be greater than the one required using low-thrust propulsion. On the contrary, the deployment time of the low-energy strategy could reasonably be shorter and it requires fewer resources from the satellite platform.

5. Conclusions

A strategy to deploy the satellites of a lunar constellation to multiple orbital planes from the same launch trajectory is investigated in this research.

The strategy is developed by taking advantage of low-energy capture trajectories existing in the Sun–Earth–Moon system, which allows the targeting of the orbital elements of the capture orbit by providing small $\Delta V$ in the neighborhood of libration point $L_1$. The theoretical background on low-energy captures and their characterization in the dynamic framework of the ER4BP are provided, and the strategy is derived from them.

The analysis verifies the strategy: it is applied to a scenario of actual interest, namely the deployment of a constellation that can provide polar and global coverage. Results from the analysis prove its effectiveness and demonstrate its performance in terms of deployment time, which does not exceed one sidereal month (26 days), and $\Delta V$, which is lower than 1 km/s for the worst case examined.

Converting $\Delta V$ to finite-time maneuvers, it can be verified that they can be successfully performed by the same chemical propulsion systems available to small satellites, requiring less than $1/3$ of the mass for a 100 kg spacecraft.