# Violation of the Leggett-Garg Inequality for Dynamics of a Bose-Einstein
Condensate in a Double-Well Potential

Tsubasa Sakamoto<EMAIL_ADDRESS>Department of Physics, Chuo
University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo 112-8551, Japan  Ryosuke Yoshii
<EMAIL_ADDRESS>Center for Liberal Arts and Sciences, Sanyo-Onoda City
University, 1-1-1 Daigaku-Dori, Sanyo-Onoda, Yamaguchi 756-0884, Japan
International Institute for Sustainability with Knotted Chiral Meta Matter
(WPI-SKCM2), Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526,
Japan  Shunji Tsuchiya<EMAIL_ADDRESS>Department of Physics, Chuo
University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo 112-8551, Japan

###### Abstract

The Leggett-Garg inequality (LGI) serves as a criterion to determine the
adherence of macroscopic system dynamics to macrorealism, as postulated by
Leggett and Garg. A violation of this inequality implies either the absence of
a realistic description of the system or the impracticality of noninvasive
measurements. In this work, we investigate the violation of the LGI for the
system of bosons in a double-well potential. Specifically, we explore the
violation of the LGI in the dynamics of bosons in a double-well potential in
the Bose-Einstein-Condensation (BEC) regime, where the system can be
considered as two weakly coupled Bose condensates, and in the single-particle
regime to establish the conditions under which the violation of the LGI
occurs. Our analysis reveals that the LGI is violated due to Josephson
oscillations, while it remains unviolated in the strong coupling regime,
attributed to the self-trapping phenomena. Notably, we observe that the
violation of the LGI becomes increasingly significant as the particle number
increases. These findings provide valuable insights into the macrorealistic
behavior of Bose condensates and highlight the effect of measurements on the
dynamics of a macroscopic system.

Introduction.—How classical behaviors of objects, which obey our intuition
about how the macroscopic world behaves, can be distinguished from non-
classical behaviors such as what are described by quantum theory? The Leggett-
Garg inequality (LGI) was proposed to answer this question [1, 2, 3]. It
serves as a test of a concept called macrorealism, which consists of two
conditions: macrorealism per se (MRPS) and noninvasive measurability (NIM).
MRPS assumes that a macroscopic system with two or more macroscopically
distinct states available to it is always in one or the other of these states,
and NIM assumes that it is possible to determine the state of the system with
arbitrary small perturbations on its subsequent dynamics. A violation of the
LGI implies a departure from macrorealism. The exploration of the LGI was
inspired by the question of whether macroscopic coherence, as illustrated by
Schrödinger’s cat gedanken experiment, can be realized in a laboratory [4].

The violation of the LGI has been studied in several systems, experimentally
and theoretically [3, 7, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15]. Recently,
experimental LG tests have been achieved in a superconducting flux qubit with
genuine macroscopicity, refuting its classical realistic description [16].

The system of Bose-Einstein condensates (BECs) in a double-well potential is
suitable for testing the LGI to clarify the boundary between classical and
quantum regimes. The dynamics of a BEC in a double-well potential has been
observed in several experiments [19, 21, 20, 22, 23, 17, 18, 25, 24]. In
particular, in the experiment of the Heidelberg group [19], Josephson
oscillations and macroscopic quantum self-trapping were successfully observed,
and the experimental result exhibits agreement with the prediction by the
Gross-Pitaevskii (GP) equation [28, 27, 26]. It implies that a population
imbalance and a phase difference between the two wells have definite values
independent of measurements in the presence of sufficiently large number of
particles in the double-well potential [29].

In this letter, we discuss the violation of the LGI for bosons in a double-
well potential. We find the violation of the LGI in the BEC regime, in which
the system can be considered as two weakly coupled BECs, despite the fact that
their dynamics are well described by coupled semiclassical equations for
definite phase and population differences, which are presumed to be unaffected
by measurements. We reveal that the violation of the LGI results from the
collapse of the state due to projective measurements in the intermediate time,
which necessarily affects the subsequent dynamics. We find that the violation
of the LGI becomes increasingly significant as the particle number increases.
Furthermore, we find that the critical value of interaction strength for the
onset of the violation of the LGI coincides with the one for self-trapping. In
the single-particle regime, in which bosons do not interact each other, the
bosons tunneling backwards induced by the measurement are found to be crucial
for the violation of the LGI.

Model.—The system of bosons in a double-well potential is well described by
the Bose-Josephson junction (BJJ) model [28, 33, 31, 35, 34, 37, 41, 36, 38,
30, 39, 40, 32]. Using the Schwinger boson representation
$\hat{S}_{x}=(\hat{a}_{\rm L}^{\dagger}\hat{a}_{\rm R}+\hat{a}_{\rm
R}^{\dagger}\hat{a}_{\rm L})/2$ and $\hat{S}_{z}=(\hat{n}_{\rm L}-\hat{n}_{\rm
R})/2$, where $\hat{a}_{\alpha}$ ($\hat{a}_{\alpha}^{\dagger}$) and
$\hat{n}_{\alpha}=\hat{a}^{\dagger}_{\alpha}\hat{a}_{\alpha}$ are annihilation
(creation) and number operators of a boson in the well $\alpha={\rm L},{\rm
R}$, the Hamiltonian of the BJJ model can be written as [36, 38, 30]

$\displaystyle\hat{H}_{\rm BJJ}=-J\hat{S}_{x}+\frac{U}{N}\hat{S}_{z}^{2},$ (1)

where $J$ is the hopping strength, $U$ is on-site interaction and $N$ is the
total number of particles. Hereafter, we set $\Lambda=U/J$ and $\hbar=1$.

Method.—The LGI is formulated for the correlation functions $C_{ij}=\langle
Q_{i}Q_{j}\rangle$ of dichotomic variables $Q_{i}=\pm 1$, which is given as
the result of the measurement at time $t_{i}$, where $i,j=1,2,3$. MRPS ensures
the existence of the value of the observable $Q_{i}$ at $t_{i}$, regardless of
whether the measurement is performed or not. Thus we can define the joint
probabilities $P_{ij}(Q_{1},Q_{2},Q_{3})$, where measurements are performed at
$t_{i}$ and $t_{j}$, whereas the system is unmeasured at time $t_{k}$ ($\neq
t_{i},t_{j}$). Under only MRPS, the joint probability may depend on when
measurements are made, since measurements at different times can affect the
time-evolution differently. By adding the NIM condition, however, the
subscripts $i$ and $j$ can be dropped since any measurement does not affect
the system under NIM. Then the joint probability $P(Q_{1},Q_{2},Q_{3})$ yields
all correlation functions and we obtain the LGI [1, 3]

$LG\equiv C_{12}+C_{23}-C_{13}\leq 1.$ (2)

We follow the dynamics of the system numerically by exact diagonalization of
the Hamiltonian (1). In evaluation of the LGI, in this work, we define the
dichotomic observable as $\hat{Q}={\rm sgn}(\hat{S}_{z})$, where ${\rm
sgn}(x)=1$ if $x\geq 0$, otherwise ${\rm sgn}(x)=-1$. For the sake of
simplicity, we set an initial state at the time $t_{1}$. Throughout this
letter, we focus on the identical time interval and define
$\tau=t_{2}-t_{1}=t_{3}-t_{2}$. Then all three two-time correlation functions
depend only on the time interval $\tau$. Based on the above setting, we
calculate $LG$ defined as $LG=C_{12}+C_{23}-C_{13}$ (also see the Supplemental
Material (SM) for details on the calculations of the correlation functions).

BEC regime.—First, we examine violation of the LGI in the BEC regime
$(\Lambda\neq 0$ and $N\gg 1)$, in which the system can be considered as two
weakly coupled Bose condensates. Such Bose condensates have been studied by a
semiclassical approach based on two coupled GP equations, which predict
Josephson oscillation or self-trapping depending on the value of $\Lambda$ and
population imbalance at initial time [19, 21, 28, 30, 43, 42].

Figure 1: (a)-(e) Time evolutions of the occupation probabilities of the left
well $P(n_{\rm L})$ for $\Lambda=0.0$ (a), $1$ (b), $1.5$ (c), $1.9$ (d), and
$2.5$ (e). Here we set $N=2000$. (f)-(h) Colormaps of $LG$ for $N=100$ (f),
$200$ (g) and $2000$ (h). The two black dashed lines in (f)-(h) denote the
critical value of self-trapping $|\Lambda_{\rm c}|=2$. (i) Violation region of
$\Lambda$. The light-blue region denotes the violation region, while the gray
region non-violation region. The red dashed line denotes the boundary between
them, which coincides with the onset of self-trapping.

Population dynamics in the BEC regime.—Figs. 1-1 present the time-evolutions
of the occupation probability of the left well $P(n_{\rm L})$ in the case of
$N=2000$ with the initial state $|N\rangle_{\rm L}|0\rangle_{\rm R}$. The peak
of the occupation probability oscillates between $n_{\rm L}=0$ and $n_{\rm
L}=N$ due to single-particle tunneling effect when $\Lambda=0$ (Fig. 1). A
similar coherent oscillation occurs in the presence of on-site interaction,
but it is rather induced by Josephson effect when $0<\Lambda<2$ (Figs. 1-1).
For strong interaction strength, the Josephson oscillation ceases and self-
trapping occurs when $\Lambda\geq 2$, where the peak of the occupation
probability oscillates keeping the population imbalance (Fig. 1). The critical
value of interaction strength for the onset of self-trapping $\Lambda_{\rm
c}=2$ agrees with the value derived in the semiclassical analysis based on the
GP equation[28].

If the dynamics obeys coupled GP equations, oscillation of the population
imbalance continues without damping [28, 31, 42]. The damping of the
oscillation that occurs in the presence of on-site interaction in Figs. 1-1
can be understood as being due to the intrinsic quantum fluctuation in the
initial state [30]. That is, since the initial population imbalance is fixed,
the initial phase difference between the Bose condensates should be
fluctuating due to the uncertainty relation between population imbalance and
phase difference [30]. As $n_{\rm L}$ and $n_{\rm R}$ becomes larger, the
quantum fluctuation ($\simeq 1/\sqrt{n_{{\rm L},{\rm R}}}$ [44]) becomes more
negligible [38, 28]. This is the reason why, with the increase of $N$, the
coherent Josephson oscillation persists for longer (see Fig. S3 in the SM).

BECs in a double-well potential are considered to possess macroscopic
observables for sufficiently large $N$. In fact, the damping is suppressed for
larger $N$, and the population dynamics is consistent with the semiclassical
description. The experimental data in Ref. [19] have also been shown to be
consistent with the prediction of the semiclassical analysis. Therefore, both
the theoretical and experimental studies support the existence of macroscopic
wave functions of two BECs with definite phase difference and population
imbalance for sufficiently large $N$. Upon such considerations, it seems
reasonable to expect non-violation of the LGI.

Violation of the LGI in the BEC regime.—Fig. 1 shows the result of the LG test
in the BEC regime for $N=2000$: $LG$ is plotted as a function of $\tau$ and
$\Lambda$. $LG$ exhibits periodic oscillations as a function of $\tau$ that
reflects the oscillation of the occupation probabilities in Figs. 1-1.
Contrary to the expectation, the LGI is violated near the peaks of the
oscillations for the interaction $-2<\Lambda<2$, where Josephson oscillations
take place for $\Lambda\neq 0$, while single-particle oscillations for
$\Lambda=0$. It implies that macrorealism is excluded even for such a large
$N$, and thus contradicts what is expected by the GP description; the LGI is
never broken for the GP description [29].

Figure 2: Effect of collapse of the state due to measurement at time $t_{2}$.
(a) Plots of $LG$ in Eq. (2) as a function of $\tau J$. The upper blue region
represents a violation of the LGI in Eq. (2). The blue (red) line denotes the
plot of $LG$ where the state is (not) collapsed at $t_{2}$. (b) $|\langle
Q_{3}\rangle_{{\rm M}_{2}}-\langle Q_{3}\rangle|$ as a function of $\tau J$.
The parameters used in (a) and (b) are $N=100$ and $\Lambda=1$.

As shown in Fig. 1, when $\Lambda$ exceeds the critical value indicated by the
dashed line, the LGI becomes non-violated regardless of $\tau$. It is
noteworthy that this critical value coincides with the one for the onset of
self-trapping $\Lambda_{\rm c}=2$ (see Fig. 1). The self-trapping prevents the
system from coherent oscillation and yields, $C_{12}=C_{13}=1$ for any
interval $\tau$ (see SM). Then $C_{12}$ and $C_{13}$ cancel out and
$LG=C_{23}\leq 1$. Thus, the LGI is non-violated when self-trapping occurs.

Fig. 1 shows that the period of oscillation of $LG$ increases and its
amplitude decreases as $\Lambda\to\Lambda_{\rm c}^{-}$. It is due to the
softening of the Josephson oscillation associated with the transition from the
Josephson regime to the self-trapping regime [38, 45].

The LGI becomes non-violated for $\tau$ greater than a characteristic value
when $0<\Lambda<\Lambda_{\rm c}$, which becomes smaller as the interaction
strength increases but it becomes longer as $N$ increases, as shown in Figs.
1-1. This behavior is due to the damping of Josephson oscillation, which
destroys the correlation in time. The Josephson oscillation indeed decays
after a characteristic time, which becomes smaller as the interaction strength
increases, as shown in Figs. 1-1, but it becomes longer as $N$ increases due
to suppression of fluctuation of phase difference of two Bose condensates.

The LGI is found to be violated symmetrically between repulsive ($\Lambda>0$)
and attractive ($\Lambda<0$) on-site interaction in Figs. 1-1 since the
Hamiltonian with repulsive interaction can be mapped to the one with
attractive interaction as $\hat{R}_{z}(\pi)\hat{H}_{\rm
BJJ}(\Lambda)\hat{R}_{z}^{\dagger}(\pi)=-\hat{H}_{\rm BJJ}(-\Lambda)$, where
$\hat{R}_{z}(\theta)=\exp\quantity(-i\theta\hat{S}_{z})$. From this
relationship, we obtain

$\displaystyle P(n_{\rm L},\Lambda,t)=P(n_{\rm L},-\Lambda,t),$ (3)

where $P(n_{\rm L},\Lambda,t)$ denotes the probability of $n_{\rm L}$ bosons
on the left-well at time $t$ for $\Lambda$ (see SM). This symmetric
relationship between repulsive and attractive on-site interaction results in
the symmetric evolution of $LG$ between them in Figs. 1-1.

Origin of the violation of the LGI.—The violation of the LGI in Fig. 1 results
from the collapse of the state due to projective measurements made at time
$t_{2}$. As shown in Fig. 2, the LGI holds if the state is assumed not to be
collapsed by measurements at $t_{2}$, while if the state is collapsed at
$t_{2}$ the LGI is violated. The effect of the collapse appears explicitly in
Fig. 2 in the comparison between $\langle Q_{3}\rangle_{{\rm M}_{2}}$, where
the state is assumed to be collapsed due to measurements at $t_{2}$, and
$\langle Q_{3}\rangle$, where measurements are not made at $t_{2}$. The
condition $|\langle Q_{3}\rangle_{{\rm M}_{2}}-\langle Q_{3}\rangle|\neq 0$
implies that the dynamics after time $t_{2}$ evolves differently depending on
whether a measurement is made at $t_{2}$ or not. Comparison between Figs. 2
and 2 shows that the LGI is violated for $\tau$ at which $|\langle
Q_{3}\rangle_{{\rm M}_{2}}-\langle Q_{3}\rangle|$ becomes large. Therefore
$Q_{3}$ of $C_{13}$ and that of $C_{23}$ evolve differently, and this fact
leads to violations of the LGI. As shown in Figs. 2 and 2, the effect of
collapse of the state on $LG$ and $\langle Q_{3}\rangle$ due to measurements
at $t_{2}$ disappears as the oscillation of population decays.

Figure 3: (a) Time evolutions of $P(n_{\rm L})$ for $z(0)=0.28$. (b) Violation
of the LGI for $z(0)=0.2$ (the red line), 0.28 (the blue line), and 0.4 (the
green line). These parameters are well below the critical value of macroscopic
quantum self-trapping. In (a) and (b), the parameters are $N=1000$ and
$\Lambda=15$.

Observability of violation of the LGI.—We numerically test the LGI for the
nearly same parameters for which Josephson oscillations were observed in Ref.
[19]. Here, we set $\Lambda=15$ and $N=1000$. The LGI is shown to be broken
for realistic parameters in Fig. 3, in which the evolution of $LG$ is plotted
as a function of $\tau$ for several initial population imbalance $z(0)=(n_{\rm
L}(0)-n_{\rm R}(0))/N$. $z(0)$ needs to be smaller than the critical value for
self-trapping $z_{\rm c}(0)\simeq 0.5$ [31]. $z(0)=0.28$ is expected to be
more desirable for the demonstration of violation of the LGI. Fig. 3 shows the
time evolution of $P(n_{\rm L})$ for $z(0)=0.28$. It exhibits tunneling of
atoms between the two wells, which indeed leads to the violation of the LGI.

Figure 4: (a) Violation of the LGI. (b) Three correlation functions $C_{12}$
(the red dashed line), $C_{23}$ (the black solid line) and $C_{13}$ (the green
dotted line) as functions of $\tau$. (c)/(d) Time evolution of $P(n_{\rm L})$
in the absence/presence of the measurement at $tJ/2\pi=0.2$ (represented by
the gray dotted line). $n_{\rm L}=70$ at $tJ/2\pi=0.2$ in (d). (e) Colormap of
$LG$ for all different initial occupation number in the left well $n_{\rm
L}(0)$. The violation of the LGI is symmetric under the reflection with
respect to the line $n_{\rm L}(0)=50$. The blue color in (e) shows the part
where the LGI is violated. In (a)-(e), the parameters are $N=100$ and
$\Lambda=0$.

Single particle regime.—Next, we study violation of the LGI in the absence of
on-site interaction. In this case the Hamiltonian is simply given by
$\hat{H}_{\rm BJJ}=-J\hat{S}_{x}$. For this model, the Schrödinger equation is
analytically solvable and the resultant wave function at time $t$ evolved from
the initial state $|N-l\rangle_{\rm L}|l\rangle_{\rm R}$ is given by

$\displaystyle\ket{\Psi(t)}=$
$\displaystyle\frac{1}{\sqrt{(N-l)!}}\frac{1}{\sqrt{l!}}\left(\hat{a}^{\dagger}_{\rm
L}\cos\frac{J}{2}t-i\hat{a}^{\dagger}_{\rm R}\sin\frac{J}{2}t\right)^{N-l}$
$\displaystyle\times\left(\hat{a}^{\dagger}_{\rm
R}\cos\frac{J}{2}t-i\hat{a}^{\dagger}_{\rm
L}\sin\frac{J}{2}t\right)^{l}|\mathrm{vac}\rangle.$ (4)

The first/second line of Eq. (Violation of the Leggett-Garg Inequality for
Dynamics of a Bose-Einstein Condensate in a Double-Well Potential) represents
the coherent oscillation of the bosons initially located in the left/right
well. These bosons flow in opposite directions. The former represents bosons
initially tunneling from the left well to the right well, as shown in Fig. 4,
while the latter represents those initially tunneling from the right well to
the left well. As shown in Fig. 4, the peak of the occupation probability
coherently oscillates between $n_{\rm L}=0$ and $n_{\rm L}=N$ without damping.
Due to this coherent oscillation the LGI is violated, as shown in Fig. 4.

Focusing on the first violation in Fig. 4, this behavior can be explained as
follows: Until $C_{13}$ starts to decrease, the LGI holds, since
$C_{12}=C_{13}=C_{23}=1$. The LGI begins to be first violated when $C_{13}$
starts decreasing at $\tau J/2\pi\simeq 0.125$. Since $C_{13}$ abruptly
becomes $-1$ whilst $C_{12}=1$ and $C_{23}$ is still positive, the LGI becomes
violated quite abruptly. The value of the peak of $LG$ becomes larger with the
increase of $N$, due to more abrupt decrease of $C_{13}$ (see SM).

The persistence of this violation is attributed to the bosons tunneling
backwards in Fig. 4 induced by the measurement, which corresponds to the
second line of Eq. (Violation of the Leggett-Garg Inequality for Dynamics of a
Bose-Einstein Condensate in a Double-Well Potential). In the case of $C_{23}$
the bosons tunneling backwards impede the decrease of $C_{23}$, whilst such
effect is absent in the case of $C_{13}$, where no measurement is performed at
$t_{2}$ (see Fig. 4). The LGI again begins to hold approximately when $C_{12}$
becomes zero at $\tau J/2\pi\simeq 0.25$ (see Fig. 4).

We also test the LGI for all different initial occupation numbers of the left
well $n_{\rm L}(0)$ with $N=100$. In Fig. 4, it can be seen that the peak
value of $LG$ increases with a larger initial population imbalance (see SM).

Conclusion and discussion.—In this work, we study violations of the LGI for
bosons in a double-well potential. We find that the LGI is violated due to
Josephson oscillations in the weak coupling regime, while it remains
unviolated in the strong coupling regime, attributed to the self-trapping
phenomena. Remarkably, we observe that the LGI remains violated even in the
large $N$ system. We also find that, in this system, on-site interaction
induces a dephasing effect that destroys temporal correlations, resulting in
the LGI no longer being violated afterwards. We have confirmed that the LGI is
broken below the critical value of interaction strength for self-trapping. In
contrast, it has been shown that the LGI is no longer violated when self-
trapping takes place. Furthermore, we have revealed that the LGI is broken for
nearly the same parameters for which the Josephson oscillations are observed
in the experimental work of Ref. [19], thus indicating the presence of
macroscopic quantum coherence in the semiclassical regime.

Our result suggests that either MRPS or NIM is violated, or possibly both, in
the system of bosons in a double-well potential with a large number of
particles. In our LG tests, we perform projection measurements, which
necessarily affect the subsequent dynamics. Therefore, if measurements can be
made without collapsing the state in any way [46], the result can be different
from ours.

Acknowledgments.—We thank G. Kimura, M. Kunimi and T. Nikuni for fruitful
discussions. ST thanks F. Dalfovo, G. Ferrari, and A. Recati for useful
discussions. RY is partially supported by a Grant-in-Aid of MEXT for
Scientific Research (KAKENHI Grant Nos. 19K14616 and 20H01838). ST is
supported by the Japan Society for the Promotion of Science Grant-in-Aid for
Scientific Research (KAKENHI Grant No. 19K03691).

## References

  * [1] A. J. Leggett and A. Garg, Phys. Rev. Lett. 54, 857 (1985).
  * [2] A. J. Leggett, J. Phys.: Condens. Matter 14, R415 (2002); Rep. Prog. Phys. 71, 022001 (2008).
  * [3] C. Emary, N. Lambert, and F. Nori, Rep. Prog. Phys. 77, 016001 (2014).
  * [4] E. Schr${\rm\ddot{o}}$dinger, Naturwiss. 23, 807 (1935).
  * [5] M. Barbieri, Phys. Rev. A 80, 034102 (2009).
  * [6] D. Avis, P. Hayden, and M. M. Wilde, Phys. Rev. A 82, 030102(R) (2010).
  * [7] C. Emary, N. Lambert, and F. Nori, Phys. Rev. B 86, 235447 (2012).
  * [8] J. Kofler and Č. Brukner, Phys. Rev. A 87, 052115 (2013).
  * [9] C. Budroni and C. Emary, Phys. Rev. Lett. 113, 050401 (2014).
  * [10] S. V. Moreira, A. Keller, T. Coudreau, and P. Milman, Phys. Rev. A 92, 062132 (2015).
  * [11] J. A. Formaggio, D. I. Kaiser, M. M. Murskyj, and T. E. Weiss, Phys. Rev. Lett. 117, 050402 (2016).
  * [12] J. Martin and V. Vennin, Phys. Rev. A 94, 052135 (2016).
  * [13] L. Rosales-Zárate, B. Opanchuk, Q. Y. He, and M. D. Reid, Phys. Rev. A 97, 042114 (2018).
  * [14] A. Matsumura, Y. Nambu, and K. Yamamoto, Phys. Rev. A 106, 012214 (2022).
  * [15] C. Mawby and J. J. Halliwell, Phys. Rev. A 107, 032216 (2023).
  * [16] G. C. Knee, K. Kakuyanagi, M.-C. Yeh, Y. Matsuzaki, H. Toida, H. Yamaguchi, S. Saito, A. J. Leggett, and W. J. Munro, Nature Commun. 7, 13253 (2016).
  * [17] Y. Shin, M. Saba, T. A. Pasquini, W. Ketterle, D. E. Pritchard, and A. E. Leanhardt, Phys. Rev. Lett. 92, 050405 (2004).
  * [18] M. Saba, T. A. Pasquini, C. Sanner, Y. Shin, W. Ketterle, and D. E. Pritchard, Science 307, 1945 (2005).
  * [19] M. Albeiz, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005).
  * [20] R. Gati, B. Hemmerling, J. Fölling, M. Albeiz, and M. K. Oberthaler, Phys. Rev. Lett. 96, 130404 (2006).
  * [21] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Nature (London) 449, 579 (2007).
  * [22] T. Berrada, S. van Frank, T. Schumm, J. F. Schaff, and J. Schmiedmayer, Nat. Commun. 4, 2077 (2013).
  * [23] C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G. Boshier, Phys. Rev. Lett. 111, 205301 (2013).
  * [24] M. Vretenar, B. Kassenberg, S. Bissesar, C. Toebes, and J. Klaers, Phys. Rev. Res. 3, 023167 (2021).
  * [25] A. Mukhopadhyay, X.-W. Luo, C. Schimelfenig, M. K. H. Ome, S. Mossman, C. Zhang, and P. Engels, Phys. Rev. Lett. 132, 233403 (2024).
  * [26] L. P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961).
  * [27] E. P. Gross, Nuovo Cimento 20, 454 (1961); J. Math. Phys. 4, 195 (1963).
  * [28] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Phys. Rev. Lett. 79, 4950 (1997).
  * [29] In the derivation of the GP equation, the field operator $\hat{\Psi}({\bf r},t)$ is replaced with its average $\langle\hat{\Psi}({\bf r},t)\rangle$, which corresponds to the macroscopic wave function of a Bose condensate, where the fluctuation is neglected [26, 27]. It implies that both MRPS and NIM hold for the GP equation, thus satisfying the LGI.
  * [30] G. J. Milburn, J. Corny, E. M. Wright, and D. F. Walls, Phys. Rev. A 55, 4318 (1997).
  * [31] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, Phys. Rev. A 59, 620 (1999).
  * [32] A. Vardi and J. R. Anglin, Phys. Rev. Lett. 86, 568 (2001).
  * [33] I. Zapata, F. Sols, and A. J. Leggett, Phys. Rev. A 57, R28(R) (1998); 67, 021603(R) (2003).
  * [34] R. Hipolito and A. Polkovnikov, Phys. Rev. A 81, 013621 (2010).
  * [35] G. Ferrini, D. Spehner, A. Minguzzi, and F. W. J. Hekking, Phys. Rev. A 82, 033621 (2010); 84, 043628 (2011).
  * [36] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler, Phys. Rev. Lett. 105, 204101 (2010).
  * [37] X.-W. Xu, A.-X. Chen, and Y.-x. Liu, Phys. Rev. A 96, 023832 (2017).
  * [38] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Rev. Mod. Phys. 90, 035005 (2018).
  * [39] Sudip Sinha and S. Sinha, Phys. Rev. E. 100, 032115 (2019).
  * [40] Sudip Sinha and S. Sinha, Phys. Rev. Lett. 125, 134101 (2020).
  * [41] A. Vardi, Phys. Rev. E 106, 064210 (2022).
  * [42] D. Ananikian and T. Bergeman, Phys. Rev. A 73, 013604 (2006).
  * [43] B. Juli${\rm\acute{a}}$-D${\rm\acute{i}}$az, D. Dagnino, M. Lewenstein, J. Martorell, and A. Polls, Phys. Rev. A 81, 023615 (2010).
  * [44] M. Gajda and K. Rza̧żewski, Phys. Rev. Lett. 78, 2686 (1997).
  * [45] A. Farolfi, A. Zenesini, D. Trypogeorgos, C. Mordini, A. Gallemì, A. Roy, A. Recati, G. Lamporesi, and G. Ferrari Nat. Phys. 17, 1359 (2021).
  * [46] Phase contrast imaging is one of the candidates for such nondestructive subsequent measurements for BECs [47, 48, 49, 50].
  * [47] M. R. Andrews, M.-O. Mewes, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Science 273, 84 (1996).
  * [48] C. C. Bradley, C. A. Sackett, and R. G. Hulet, Phys. Rev. A 55, 3951 (1997).
  * [49] J. M. Higbie, L. E. Sadler, S. Inouye, A. P. Chikkatur, S. R. Leslie, K. L. Moore, V. Savalli, and D. M. Stamper-Kurn, Phys. Rev. Lett. 95, 050401 (2005).
  * [50] E. O. Ilo-Okeke and T. Byrnes, Phys. Rev. Lett. 112, 233602 (2014).


# Transfers between moons with escape and capture patterns via Lyapunov
exponent maps

David Canales 111Ph.D. Candidate, Purdue University; Currently, Assistant
Professor in the department of Aerospace Engineering in Embry–Riddle
Aeronautical University, AIAA Member1 Kathleen C. Howell222Hsu Lo
Distinguished Professor of Aeronautics and Astronautics, School of Aeronautics
and Astronautics, AIAA Fellow2 Elena Fantino333Associate Professor, Aerospace
Engineering Department, AIAA Member3 Annika J. Gilliam444Undergraduate
Research Assistant, Aerospace Engineering Department

###### Abstract

This contribution focuses on the design of low–energy transfers between
planetary moons and presents an efficient technique to compute trajectories
characterized by desirable behaviors in the vicinities of the departure and
destination bodies. The method utilizes finite–time Lyapunov exponent maps in
combination with the Moon–to–Moon Analytical Transfer (MMAT) method previously
proposed by the authors. The integration of these two components facilitates
the design of direct transfers between moons within the context of the
circular restricted three–body problem, and allows the inclusion of a variety
of trajectory patterns, such as captures, landings, transits and takeoffs, at
the two ends of a transfer. The foundations and properties of the technique
are illustrated through an application based on impulsive direct transfers
between Ganymede and Europa. However, the methodology can be employed to
assist in the design of more complex mission scenarios, such as moon tours.

Keywords: Multi–Body Dynamics, Circular Restricted Three–Body Problem,
Trajectory Design, Moon–To–Moon Transfers, Dynamical Systems Theory, Galilean
Moons

## Nomenclature

$\Delta{v}$ = | magnitude of velocity variation
---|---
$t_{tot}$ = | total time–of–flight for a transfer
$\delta$ = | perturbation
x = | six–dimensional vector accounting for three position and three velocity components
$\phi(t_{f},t_{0})$ = | state transition matrix (STM) between $t_{f}$ (final time) and $t_{0}$ (initial time)
$n$ = | integer used to enumerate, e.g., $n$ = 1, 2
$m_{p}$ = | mass of the planet
$m_{m}$ = | mass of the moon
$\textbf{r}_{p-s/c}$ = | position vector of the spacecraft relative to the planet
${r}_{p-s/c}$ = | distance between the planet and the spacecraft
$\textbf{r}_{m-s/c}$ = | position vector of the spacecraft relative to the moon
${r}_{m-s/c}$ = | distance between the moon and the spacecraft
$\textbf{r}_{p}$ = | position vector of the planet relative to the origin of the reference frame
$\textbf{r}_{m}$ = | position vector of the moon relative to the origin of the reference frame
$\theta$ = | true anomaly
$t$ = | time
CR3BP = | circular restricted three–body problem
$\mu$ = | mass ratio in the CR3BP
$e$ = | orbital eccentricity
$a$ = | orbital semi–major axis
$i$ = | orbital inclination
$\Omega$ = | right ascension of the ascending node of an orbit
$\omega$ = | argument of periapsis of an orbit
${L_{a}}^{*}$ = | reference length for normalization
r = | position of the spacecraft in the planet–moon barycentric rotating reference frame
$\dot{\textbf{r}}$ = | velocity of the spacecraft in the planet–moon barycentric rotating reference frame
$U^{*}$ = | pseudo–potential function in the CR3BP
$J$ = | Jacobi constant
$L_{1},L_{2},\dots,L_{5}$ = | five equilibrium points of the CR3BP
$x$ = | coordinate on the $\hat{x}$–axis of the rotating reference frame
$y$ = | coordinate on the $\hat{y}$–axis of the rotating reference frame
$z$ = | coordinate on the $\hat{z}$–axis of the rotating reference frame
$\lambda$ = | eigenvalue
$\Sigma$ = | hyperplane associated to a Poincaré section
$\theta_{0_{m}}$ = | true anomaly of the orbit of a moon measured from the ascending node at the initial epoch
$R_{SoI}$ = | radius of the sphere of influence of a moon
$d_{SoI}$ = | ratio between the gravitational acceleration of the moon and that of the planet
$\theta_{Int}$ = | true anomaly of the point of intersection between two confocal conic sections
$\iota$ = | singular values of a matrix
$\Upsilon$ = | diagonal matrix incorporating all singular values
$V$ = | matrix that incorporates the direction of stretching at $t_{0}$
$U$ = | matrix that incorporates the direction of stretching at $t_{f}$
$C$ = | Cauchy–Green Strain Tensor matrix

##  1 Introduction

The recently released Decadal Strategy for Planetary Science and Astrobiology
2023–2032 [1] prioritizes new missions to the gas giants and their moons
(e.g., Enceladus multiple flyby and lander, Saturn probe, Titan orbiter,
Europa lander, Neptune–Triton probe, Uranus orbiter) and firmly recommends the
realization of NASA’s Europa Clipper [2]. The scientific community believes
that the open questions regarding the Sun’s planetary system can only be
addressed through a systematic exploration of the icy worlds, with particular
emphasis on the in situ observation of planetary moons. In this scenario, the
development of efficient tools to design trajectories enabling the execution
of transfers between moons and tours of planetary systems is crucial.

Conventional trajectory design methods based on patched conics and
multi–gravity assist have been extensively applied to real mission scenarios
(e.g., JUpiter ICy moons Explorer – JUICE [3], JIMO [4]) and have been the
focus of numerous studies (see, e.g., [5, 6, 7]). Over the past two decades,
investigations in terms of dynamical systems have demonstrated that it is
possible to fly a spacecraft on low–energy trajectories departing from and
leading to the vicinity of the libration points in the circular restricted
three–body problems (CR3BPs) composed of a planet and individual moons. The
Petit Grand Tour (PGT) has been the first concept of a low–energy tour of the
Jovian system [8, 9]. Here, hyperbolic invariant manifolds of libration point
orbits (LPOs) in the CR3BPs associated with Jupiter and distinct moons are
propagated in the space between the moons, and their intersections are used to
design direct impulsive transfers; a trajectory between LPOs in the vicinity
of Ganymede and Europa costs 1.214 km/s and takes 25 days. In the Multi-Moon
Orbiter (MMO) concept [10, 11, 12], the spacecraft executes several resonant
gravity assists with the moons, reducing drastically the propellant
consumption to tens of m/s at the expense of increased times of flight
(several years). Grover and Ross [13] employed the Keplerian map and mitigated
the long transfer times of the MMO through the introduction of ad hoc small
impulsive maneuvers (summing to $\Delta v$s of 100 m/s). The investigation by
Campagnola and Russell [14] on $V_{\infty}$-Leveraging maneuvers (VILMs) led
to the identification of Ganymede–to–Europa transfers including endgames (the
departure and arrival conditions are low circular orbits) with a total cost of
1.71 km/s (of which 1.41 km/s resulted from achieving escape and capture) and
a minimum time of flight of 151 days. The study of moon tours culminated with
the blending of resonance hopping transfers and multi-body dynamics with
different types of begin and end games, such as circular orbits about the
departure and arrival moons ($\Delta v$ = 1.25 km/s and time of flight of 300
days for the Ganymede–to-Europa transfer) [15], general low-energy initial and
final states (59.5 m/s and 158.5 days) [16], and halo orbits near the
collinear libration points in the vicinity of the two moons (55 m/s and 205
days) [17].

Following up on the direct transfers developed within the PGT, Fantino and
Castelli [18] introduced a patched two–body/three-body (2BP-CR3BP) model to
facilitate the design of minimum–cost single-impulse moon–to–moon trajectories
in the Jovian system using invariant manifolds of planar Lyapunov orbits in
two dimensions (2D). A preliminary extension to trajectories between three-
dimensional (3D) halo orbits (Fantino et al. [19]) was completed by Canales et
al. through the development of an analytical method, termed the Moon–to–Moon
Analytical Transfer (MMAT) technique, to construct impulsive transfers between
2D and 3D LPOs of planet–moon CR3BPs [20, 21, 22]. Invariant manifold
trajectories emanating from a departure and a destination LPO are propagated
in the respective CR3BPs to the limit of the sphere of influence for the
respective moon, where the states of the spacecraft are expressed in a
planet–centered inertial frame and used to produce orbital elements of
osculating Keplerian orbits. Thus, the problem of connecting trajectories
originating from or leading to distinct moons translates into the analytical
computation of the intersection between confocal ellipses, the derivation of
the conditions under which such intersections exist, and the evaluation of the
transfer performance in terms of cost and time of flight. The method
incorporates the inclination of the moon orbits, and, in addition to
single–impulse transfers, can solve problems with intermediate arcs (two– and
three-impulse scenarios) and plane–change maneuvers in a variety of systems,
including trajectories between the Martian moons [23, 24].

For single-impulse direct trajectories between planar Lyapunov orbits at
Ganymede and Europa, the MMAT approach yields a $\Delta v$ of 0.94 km/s [21],
consistent with the results available in the open literature for these types
of trajectories. The predicted time of flight is 9.5 days. In its original
formulation, MMAT deals only with the escape and capture phases of a transfer,
and does not resolve the initial and final portions for departing and
inserting into the science orbit around each moon, i.e., the so-called begin
game/end game problem mentioned previously. The objective of the present
contribution is to lay the foundations of a strategy to link the inter-moon
transfer and the end game design through some desired trajectory patterns. The
methodology is applicable within the context of impulsive direct trajectories
(MMAT scenarios) as well as more efficient and practical moon–tour design
methods (such as resonance hoppings [15, 16, 17]).

With the aid of chaos indicators, the trajectories that depart or approach the
vicinity of a moon can be classified in terms of the motion patterns that they
exhibit in close proximity to the target. In particular, it is possible to
discriminate among temporary captures (with one or more revolutions around the
moon), escapes, takeoffs and impacts. The first use of chaos indicators for
spacecraft trajectory design in multi-body environments is due to Lara et al.
[25] and Villac [26], who employed the Fast Lyapunov Indicator, well-known in
dynamical astronomy. In this work, the properties of finite–time Lyapunov
exponents (FTLEs) [27] and their associated scalar fields that measure the
largest stretching direction in the flow associated with the dynamical
differential equations are utilized to distinguish phase–space regions
corresponding to distinct motion patterns. Low–energy transfers with specific
departure and arrival behaviors are designed by coupling FTLE maps and the
MMAT method. The theoretical foundations, the properties and the advantages of
the technique are illustrated through the classical case of direct single-
impulse transfers between Ganymede and Europa.

The article is organized as follows. Section 2 presents the dynamical model
and summarizes the background on MMAT and FTLEs. Section 3 delves into the
basic components of the method, i.e., the moon–to–moon access map, that
facilitates the selection of trajectory patterns at the beginning and end of a
transfer. Sections 5 to 7 include an analysis of the dependence of the outcome
on parameters such as the Jacobi constant and the departure epoch, whereas
Section 4 describes the choice of specific motion patterns on the basis of
cost, time–of–flight and departure dates through appropriate inspection maps.
Section 8 compares inward and outward transfers, and Section 9 delivers the
concluding remarks.

A preliminary version of the study has been presented by Canales et al. [28],
whereas the application of the technique to a transfer involving three moons
(Io, Europa and Ganymede) has been illustrated in [29].

##  2  Background and methodology

After defining the dynamical model, this section provides a summary of the
MMAT method and an overview of FTLEs and FTLE maps. The discussion employs
direct single-impulse transfers from Ganymede to Europa as applications to
illustrate the methodology and its features.

###  2.1 Dynamical model

Figure 1: Sketch representing the spatial 2BP–CR3BP patched model.

In a transfer between moons, the spacecraft is subject to the gravitational
attraction of multiple bodies. In the CR3BP [30], only two masses (in this
case, the planet and one moon) affect the motion of the spacecraft. This
approximation is not sufficient to obtain accurate moon–to–moon trajectories,
but the results can be refined in a high–fidelity full-ephemeris model. In
this investigation, a spatial 2BP–CR3BP patched model (Fig. 1) is adopted: in
the vicinity of each moon, the motion is modeled in the CR3BP of the
corresponding planet–moon system; far from the moon, the dynamics of the
spacecraft are approximated with the planet–s/c 2BP problem, where the real
inclinations of the orbital planes of the moons are incorporated.

The CR3BP provides an appropriate and convenient framework to study the
spacecraft motion in the vicinity of a moon. In this model, the planet (mass
$m_{p}$) and the moon (mass $m_{m}$) move in circular orbits around the center
of mass of the system. Although they exert gravitational attraction upon the
spacecraft, the latter does not affect their motion. Additionally, the
primaries are spherical and their mass distribution is homogeneous, i.e.,
gravity field irregularities are disregarded. A suitable normalization of
distances, masses and angular velocities and the adoption of an appropriate
barycentric rotating reference frame leads to the following set of
dimensionless differential equations for the motion of the spacecraft [31]:

$\displaystyle\ddot{x}-2\dot{y}=\frac{\delta{U^{*}}}{\delta{x}};\;\;\;\ddot{y}+2\dot{x}=\frac{\delta{U^{*}}}{\delta{y}};\;\;\;\ddot{z}=\frac{\delta{U^{*}}}{\delta{z}}.$
(1)

Here, $\mu=\frac{m_{m}}{(m_{m}+m_{p})}$ is the mass ratio of the system,
whereas the term
$U^{*}=\frac{1-\mu}{r_{p-s/c}}+\frac{\mu}{r_{m-s/c}}+\frac{1}{2}(x^{2}+y^{2})$
represents the pseudo–potential function, $r_{p-s/c}$ and $r_{m-s/c}$ being
the distances of the spacecraft to the planet and the moon, respectively. The
$\hat{x}$–axis of this rotating reference frame contains both primaries (the
planet at $\textbf{r}_{p}=[-\mu,0,0]^{T}$ and the moon at
$\textbf{r}_{m}=[1-\mu,0,0]^{T}$), whereas the $\hat{z}$–axis is aligned with
their orbital angular momentum. In Eq. (1), $\textbf{r}_{rot}=[x,y,z]^{T}$ and
$\dot{\textbf{r}}_{rot}=[\dot{x},\dot{y},\dot{z}]^{T}$ are the position and
the velocity of the spacecraft, respectively. The CR3BP admits five
equilibrium positions (denoted as libration points and labelled $L_{1}$,
$L_{2}$, …, $L_{5}$) i.e., points where the acceleration is zero if the third
body is at rest. As the libration points exhibit linear stability properties,
orbits around them are categorized by families of periodic and quasi–periodic
orbits [32, 33]. The planar Lyapunov orbit family is of interest in this work.
In the vicinity of a moon, hyperbolic invariant manifolds extending from
periodic orbits operate as pathways, including connections to other periodic
orbits within the same system [34]. Additionally, stable manifolds arrive near
a periodic orbit while unstable manifolds depart from its vicinity. Transit
orbits, as defined for this investigation, reach the moon vicinity through the
$L_{1}$ and $L_{2}$ gateways defined by periodic orbits, revolve around the
moon and either depart again or collide with it. The Jacobi constant $J$, that
represents the energy of the system, is a fundamental quantity in the CR3BP.
It is defined as

$\displaystyle J=2U^{*}-(\dot{x}^{2}+\dot{y}^{2}+\dot{z}^{2}).$ (2)

The model in the CR3BP is well-known and serves as the foundational framework
in multi-body regimes.

###  2.2 The MMAT method

Figure 2: Scheme illustrating the MMAT method for constructing direct
transfers between moons.

In the patched 2BP-CR3BP, within the Sphere of Influence (SoI) of the moon the
motion is modeled in the planet-moon CR3BP. When the trajectories cross the
surface of the SoI, they are approximated in the planet–s/c 2BP as conic
sections with one focus at the planet. Therefore, the motion outside the SoI
is considered Keplerian and completely determined by six osculating orbital
elements computed at the surface of the SoI. In this investigation, the
semi–major axis ($a$), the eccentricity ($e$) and the true anomaly ($\theta$)
are key parameters for moon–to–moon transfer design. In this patched 2BP–CR3BP
model, connections between conics that depart from and arrive at distinct
moons are explored through analytical methods. The SoI of a moon is defined as
a spherical region centered at the moon with radius $R_{SoI}$ equal to the
distance from the moon along the $x$–axis at which the ratio ($d_{SoI}$)
between the gravitational acceleration due to the moon and that caused by the
planet equals a certain small quantity, which is a free parameter. The value
adopted for this investigation for the SoI of both Ganymede and Europa is
$d_{SoI}=5\times 10^{-4}$. The selection of this value affects the design of
the moon–to–moon transfer and is examined for this particular application by
Canales et al. [21].

This study incorporates the inclinations of the orbital planes of the moons
and, therefore, the analysis is developed in three-dimensions (3D). Although
there exist many possible departure and arrival paths, this 3D nature makes
the problem dependent on the relative orbital phase between the moons at the
departure epoch. The MMAT method (see Fig. 2) previously presented by the
authors [21] utilizes the patched 2BP–CR3BP: the trajectories arriving or
departing a moon vicinity are represented using conic arcs outside the moons’
SoI. In a transfer from Ganymede to Europa, departure conics approximate the
departure Jupiter–Ganymede (J–G) CR3BP orbits, while arrival conics
approximate the arrival Jupiter–Europa (J–E) CR3BP orbits. The approximation
through planet–centered Keplerian orbits neglects the gravitational attraction
of the moons, but the associated error is small (see [18]) and the solutions
can be efficiently transitioned to a high–fidelity ephemeris model.

The necessary condition for an arrival conic to intersect spatially with a
departure conic can be expressed analytically as:

$\displaystyle
a_{a}(1-e_{a})\leq\frac{a_{d}(1-{e_{d}}^{2})}{1+{e_{d}}\cos(\theta_{d_{Int}}+n\pi)}\leq{a_{a}(1+e_{a})},\;\;\text{with
}n=0,1.$ (3)

Here, $a_{a}$ and $a_{d}$ are the semi–major axes of the arrival and departure
conics, respectively. Similarly, $e_{a}$ and $e_{d}$ are the respective
eccentricities. The true anomalies $\theta_{d_{Int}}$ and
$\theta_{d_{Int}}+\pi$ indicate the two geometrical configurations for which
intersections between the planes of the departure and arrival conics exist.
For a given departure epoch, if the condition in Eq. (3) is satisfied, a
suitable orbital phase for the target arrival moon can be determined. The
correct phasing for the two moons yields a transfer between the conic arcs
through a single impulsive $\Delta{v}$. Each possible transfer is
characterized by a total time–of–flight $t_{tot}$ and requires an impulsive
maneuver of magnitude $\Delta{v}$. In summary, given the departure epoch from
one moon, the MMAT technique identifies possible arrival conditions and
moon–to–moon transfers with one impulse and different performance
characteristics.

###  2.3 Finite–time Lyapunov exponents

The concept of FTLEs is related to that of Lagrangian Coherent Structures
(LCS) [35, 36, 37, 38, 39], i.e., regions bounding different behaviors in a
dynamical flow. Dynamical systems theory leverages different techniques with
two–dimensional Poincaré maps to investigate the long–term dynamics in the
CR3BP. The Cauchy–Green Strain Tensor (CGST) [40] effectively describes the
time evolution of the flow resulting from a perturbation. To introduce the
CGST, the state transition matrix (STM) is required. The STM provides a
relationship between a variation of the initial state $\delta{\textbf{x}_{0}}$
and the resulting deviation of the final state $\delta{\textbf{x}_{f}}$ via
the linear mapping:

$\displaystyle\delta\textbf{x}_{f}=\phi(t_{f},t_{0})\delta\textbf{x}_{0}.$ (4)

The variables $t_{0}$ and $t_{f}$ represent the initial and final times along
the trajectory, respectively, whereas $\phi(t_{f},t_{0})$ is the STM. Through
a singular–value decomposition (SVD) [41] of the STM, it is possible to
characterize a perturbation behavior as divergence or convergence. Such
characteristics of the local phase space are encoded in singular values
($\iota_{i}$) defined by the direction of stretching or contraction
$\textbf{V}_{i}$ at $t_{0}$. The SVD is expressed as:

$\displaystyle\phi(t_{f},t_{0})=\boldsymbol{U}\Upsilon{\boldsymbol{V}^{T}}.$
(5)

Here boldface letters represent matrices. Then, $\boldsymbol{U}$ and
$\boldsymbol{V}$ are mutually orthogonal: $\boldsymbol{V}$ identifies the
directions of stretching at $t_{0}$ for every singular value, whereas the
elements of the columns of $\boldsymbol{U}$ define the directions of
stretching (or contraction) at $t_{f}$. Additionally, $\Upsilon$ is a diagonal
matrix representing the stretching magnitude, written with the singular values
in descending order ($\iota_{1}>\iota_{2}>...>\iota_{n}$):

$\displaystyle\Upsilon={\begin{bmatrix}\iota_{1}&0&\dots&0\\\
0&\iota_{2}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\dots&\iota_{n}\\\
\end{bmatrix}}.$ (6)

Figure 3 sketches the SVD in a simplified two–dimensional space: in this
scheme, $\iota_{1}$ and $\iota_{2}$ indicate the smallest and largest
stretching, respectively, revealing which directions are more or less
sensitive to perturbations, with larger magnitudes indicating higher
sensitivity. Once the STM is defined, the CGST describes the deformation of
the flow as the product of the STM and its transpose [42]:

$\displaystyle\boldsymbol{C}(t_{f},t_{0})=\phi^{T}(t_{f},t_{0})\phi(t_{f},t_{0}).$
(7)

The eigenvalues ($\lambda$) of $\boldsymbol{C}$ are obtained from the
eigen–decomposition of the CGST, and are related to $\Upsilon$ through
$\lambda_{i}={\iota_{i}}^{2}$. Additionally, $\boldsymbol{C}$ and
$\boldsymbol{U}$ possess the same eigenvectors. Given that the largest
singular value corresponds to the largest perturbation growth, the FTLEs are
defined based on $\iota_{1}$ and the propagation time along the trajectory:

$\displaystyle\text{FTLE}=\frac{\iota_{1}}{|\Delta{t}|}.$ (8)

Here $\Delta{t}=t_{f}-t_{0}$. Hence, all stretching directions produce an
FTLE, but only the direction of maximum stretching is employed in Eq. (8).
This value is the one of interest for the definition of LCS. In summary, FTLEs
measure the relative phase–space element growth and contraction over a time
interval related to the system flow.

Figure 3:  Cauchy–Green Strain Tensor associated eigenvector stretching.

###  2.4 FTLE maps

The benefits of FTLEs emerge when paired with Poincaré section
representations. Finite–time Lyapunov exponent maps have been proven useful by
various authors (see, e.g., [37]) to provide quantitative information about
the propagation of a trajectory and the type of motion that occurs at a given
energy level after departure or before arrival at a moon vicinity. Possible
behaviors that are discerned include capture orbits, collision paths (or
landings), departure trajectories (or takeoffs) and close passages by the
moons (or transits). Additionally, FTLE maps illustrate separated flow
patterns entering the region of the moon through the corresponding gateway
(either $L_{1}$ or $L_{2}$). They have been recently employed to design
transfers that approach Oberon and Titania in the Uranus system and are
characterized by a desired behavior [43]. These low–energy trajectories are
produced by coupling Uranus–Oberon and Uranus–Titania CR3BPs and assuming that
the moons revolve in coplanar orbits.

(a) Departure trajectories from Ganymede vicinity through $\boldsymbol{L_{1}}$
in the J–G rotating frame.

(b) Departure FTLE map at a $\boldsymbol{J_{d}=3.007538}$ with a propagation
time $\boldsymbol{t_{d}\approx-11.4}$ days.

(c) Europa arrival trajectories through $\boldsymbol{L_{2}}$ in the J–E
rotating frame.

(d) Arrival FTLE map at a $\boldsymbol{J_{a}=3.00240}$ and with propagation
time $\boldsymbol{t_{a}\approx 5.65}$ days.

Figure 4:  (a) Departure trajectories from Ganymede and (b) departure FTLE
map; (c) Trajectories arriving at Europa and (d) arrival FTLE map.

The computation of FTLE maps involves a number of parameters, e.g., the
location of the Poincaré section, the propagation time and the value of the
Jacobi constant [37]. To illustrate this methodology, trajectories departing
Ganymede and arriving at Europa are employed (see Table 2 for relevant system
data).

Table 2: Orbital data for Ganymede and Europa relative to the Ecliptic and Equinox of J2000.0 Jupiter–centered frame [44]. Moon | Semi–major | Orbital | CR3BP | Eccentricity | Inclination | Longitude
---|---|---|---|---|---|---
| axis | period | mass ratio |  |  | asc. node
|  [$10^{5}$ km] | [day] |  [$10^{-5}$] |  [$10^{-3}$] | [degree] | [degree]
Europa | $6.713$ | 3.554 | 2.528 | 9.170 | 2.1 | 331.4
Ganymede | $10.706$ | 7.158 | 7.804 | 2.542 | 2.2 | 340.3

The Jacobi constant values are ${J}_{d}=3.00754$ and ${J}_{a}=3.00240$,
respectively in the Ganymede departure and Europa arrival CR3BPs. Separate
FTLE maps are created for departure and arrival trajectories (Fig. 4). The
departure map is based on a departure from Ganymede via a Poincaré section
near the Jupiter–Ganymede $L_{1}$ gateway. In Fig. 4(a), the section is
denoted $\Sigma_{departure}$ and is mapped over an interval of normalized time
through the propagation of all the states onto the section. The departure
section is created at $x=0.965$ in the J–G rotating frame, where all departure
states are generated. The $y$ coordinate varies from $-0.006$ to 0.015, while
the range for $\dot{y}$ extends from -0.01 to 0.02. Finally, $\dot{x}$ is
computed for each state by means of Eq. (2). Additionally, in this case, the
interval $t_{d}$ equals $-11.4$ days, corresponding to $-10$ normalized J–G
CR3BP time units. Note that $t_{tot}$ is a crucial parameter for computing the
FTLEs. A numerical analysis has been accomplished to select a value of
$t_{tot}$ that provides different regions of interest within the FTLE maps
and, thus, different trajectory behaviors over the given time interval. Recall
that $t_{tot}$ and the $x$-coordinate of the section are parameters to be
selected when computing FTLE maps.

Similarly, the arrival map (Fig. 4(d)) is constructed on a Poincaré section
($\Sigma_{arrival}$) near the J–E $L_{2}$ gateway. In Fig. 4(c),
$\Sigma_{arrival}$ is the plane $x=1.028$ in the J–E rotating frame. A grid is
generated by varying the $y$ coordinate from $-0.018$ to 0.05 and $\dot{y}$
from $-0.04$ to 0.01. In this study, the step size adopted to produce the
departure and arrival maps is 0.0001 for both coordinates. Using Eq. (2),
$\dot{x}$ is computed for each state. The propagation time interval for
arrival is $t_{a}=10$ normalized Jupiter–Europa CR3BP time units, equivalent
to $5.65$ days.

The selected propagation times for departure and arrival ensure a sufficient
emergence of LCS that include different types of behaviors over a considerable
time. The FTLE maps depend on the $x$ coordinate of the sections and on the
propagation time [45]. Note that, in this work, the normalized units in the
FTLE maps are those of the planet–moon CR3BP for which the map is being
constructed. The trajectories that transit throughout the vicinity of the
moons are contained in the hyperbolic invariant manifolds associated with the
planar Lyapunov orbits at the given Jacobi constant levels near $L_{1}$ and
$L_{2}$, respectively (i.e., unstable for J–G CR3BP, stable for the J–E
CR3BP).

As shown in Fig. 5, strainlines bound regions (or “lobes”) of qualitatively
similar motion. This property allows the designer to select initial conditions
that follow a desired behavior since all the initial conditions in any such
“lobe” lead to similar trajectory patterns. As observed in Fig. 6, FTLE maps
separate initial conditions corresponding to trajectories that enter the
moon’s vicinity from those that do not approach the moon. Finally, within the
maps themselves, strainlines represent various behaviors and outcomes (see
Fig. 7). Hence, they inform in the selection of initial conditions for
specified departure and arrival characteristics. Every “lobe” in Fig. 7 yields
a specific motion pattern, such as captures, tours and collisions. Moreover,
the separatrices between “lobes” identify collision trajectories because when
an impact occurs, the propagation time of the trajectory is relatively short,
and this produces large values of the FTLEs (e.g., 1.352 for the collision
trajectory in Fig. 7), resulting in distinctive patterns. Finally, low FTLE
values are associated with the centers of the “lobes”. This fact is justified
since smaller FTLEs (e.g., 0.315 for the selected tour trajectory in Fig. 7)
involve a lower sensitivity to perturbations. It is important to emphasize
that the only purpose of the raw FTLE values is to identify distinct “lobes”
that generate different trajectory behaviors.

Figure 5: Three departure trajectories with the same departure pattern and the
same associated isoline upon the FTLE map. Figure 6: Different isolines in the
departure FTLE map corresponding to transit and non–transit trajectories.
Figure 7: Three departure trajectories with different patterns and their
distinct isolines on the departure FTLE map ($J_{d}=3.00754$).

##  3 Designing moon–to–moon transfers with different trajectory patterns: the
moon–to–moon access maps methodology

Merging the MMAT method with FTLE maps facilitates the design of moon–to–moon
transfers with selected behaviors in the vicinity of the moons. Consider the
framework for transfers from Ganymede to Europa in Fig. 4. Recall that
moon–to–moon transfers are dependent on the departure epoch because the moons
revolve in their true orbital planes [21]. Equation (3) must be satisfied for
the Ganymede and Europa vicinities to be connected with a single impulsive
maneuver. Therefore, the selection of a suitable departure epoch is critical
for the identification of effective transfers. Assume that all initial
conditions depart Ganymede at epoch $\theta_{0_{Gan}}=82.506^{\circ}$111Here
time is expressed as the true anomaly of the moon, measured from the ascending
node of its orbit.. Experiments demonstrate that this value leads to the
largest number of connections between Ganymede and Europa for the specified
departure and arrival energy levels (see below). Then, since outside the SoIs
the trajectories are approximated by conics in the Jupiter–centered 2BP, the
central expression in Eq. (3) is evaluated for every departure conic generated
by each initial condition on the Poincaré section in the vicinity of the
$L_{1}$ gateway. Note that $\Sigma_{departure}$ lies within the SoI. The map
depicted in Fig. 8(a) is the moon–to–moon tides map. Its color gradients
reflect the value of the expression
$\frac{1}{{L_{a}}^{*}}\frac{a_{d}(1-{e_{d}}^{2})}{1+{e_{d}}\cos(\theta_{d_{Int}})}$,
i.e., the central term of Eq. (3), as each initial condition is propagated
from the Poincaré section to the SoI of Ganymede. This expression is labelled
the “$\theta_{d_{Int}}$ configuration” and includes the reference length
${L_{a}}^{*}$ for normalization in the J–E CR3BP. Observe that the
moon–to–moon tides map and the departure FTLE map have the same shape;
however, the former uses the propagation of initial conditions towards the SoI
and represents the “$\theta_{d_{Int}}$ configuration”, whereas the latter
employs the propagation from the departure moon and illustrates the FTLE. The
Europa arrival states are propagated backwards in time from $\Sigma_{arrival}$
towards the SoI of Europa. Then, apoapsis and periapsis arrival maps (Figs.
8(b) and (c), respectively) are produced by evaluating the expressions on the
left and right sides of Eq. (3) ($a_{a}(1+e_{a})/L_{a}^{*}$ and
$a_{a}(1-e_{a})/L_{a}^{*}$, respectively) for each arrival condition at the
SoI. Note that the reference length for the J–E CR3BP is used to normalize all
three maps. The moon–to–moon tides map and the apoapsis and periapsis arrival
maps are then matched to identify initial conditions that satisfy Eq. (3). For
one–impulse transfers to be feasible, the middle term in this equation must be
simultaneously $\leq a_{a}(1+e_{a})$ (upper constraint) and $\geq
a_{a}(1-e_{a})$ (lower constraint). This methodology is denoted the
moon–to–moon access maps approach and is sketched in Fig. 9: the gray region
corresponds to the Ganymede–to–Europa transfers that satisfy Eq. (3).

Figure 8: Moon–to–moon tides map (a) and apoapsis (b) and periapsis (c)
arrival maps matched with 1.242 isoline for transfers from Ganymede towards
Europa. Figure 9: Schematic of the moon-to-moon access maps method.

A specific isoline is selected to accomplish the matching in Fig. 8:

$\displaystyle\frac{{a}_{d}(1-{{e}_{d}^{2})}}{1+{e}_{d}{\cos({\theta}_{{d}_{Int}}})}={a}_{a}(1+{e}_{a})=1.242\;{L}_{a}^{*}.$
(9)

For transfers from an outer to an inner moon, the moon–to–moon tides map and
the apoapsis arrival map are matched with respect to the requirement
$\frac{a_{d}(1-{e_{d}}^{2})}{1+{e_{d}}\cos(\theta_{d_{Int}}+n\pi)}\leq
a_{a}(1+e_{a})$. The uniform color of the periapsis arrival map in Fig. 8
means that the periapsis radii of all arrival trajectories are similar. Note
that since these are smaller than the middle term in Eq. (3), the lower
constraint is always satisfied for this particular application. Similarly, the
upper constraint is always satisfied for a transfer from an inner to an outer
moon, and the middle term of Eq. (3) must be larger than the arrival conic
periapsis radius. Within the MMAT access maps, red arrows indicate accessible
regions, i.e., sets of initial conditions for which Eq. (3) is fulfilled and
the construction of one–impulse trajectories from Ganymede to Europa with
different trajectory patterns is possible; in other words, an arrival epoch
exists and the transfer is feasible. Projecting the isolines onto the
departure and arrival FTLE maps (Fig. 10) allows inspection for distinct
transfers between Ganymede and Europa that are available leveraging
strainlines associated with different patterns. Consequently, it is possible
to select a feasible transfer between trajectories exhibiting desired
behaviors in the vicinity of each moon. For example, the selected initial
conditions for the departure and arrival FTLE maps in Fig. 10 lead to the
sample transfer plotted in Fig. 11 in which the spacecraft, after completing
two revolutions of Ganymede and implementing one impulsive maneuver, reaches
the interior region of the J–E CR3BP and transits through the vicinity of
Europa. Henceforth, the following notation is adopted to illustrate moon-to-
moon trajectories in the Jupiter-centered frame (e.g., see Fig. 11): instant
t0 denotes the beginning of the transfer; label 0 corresponds to the crossing
of the departure section; label 1 indicates the time at which the departure
arc intersects the surface of the departure moon SoI; 2 represents the
intersection between departure and arrival conics; 3 is the moment in which
the arrival conic crosses the surface of the arrival moon SoI; 4 is the time
of arrival of the spacecraft at the arrival section; finally, tf denotes the
end of the transfer. Moreover, the following color scheme is used to link each
time label with a specific body: black indicates the spacecraft, orange refers
to the departure moon and green corresponds to the arrival moon.

Figure 10:  Available trajectories from Ganymede to Europa corresponding to
the 1.242 isoline shown through the departure and arrival FTLE maps. Figure
11:  Transfer from Ganymede to Europa with a desired behavior, as selected in
the FTLE maps.

##  4 Selection of trajectory patterns through inspection maps

The design of moon–to–moon transfers is greatly aided by combining MMAT access
maps with the information provided by FTLE maps, as demonstrated in Sect. 3.
Once the Jacobi constant values in the arrival and departure CR3BPs are
selected, trajectories that yield a moon–to–moon transfer are identified,
assuming that Eq. (3) is satisfied. Moreover, transfers with specific mission
objectives near each moon are produced. For example, it is possible to
determine all the departure conditions that yield a given arrival trajectory
(Design 1). Alternatively, one can identify all the feasible arrival
conditions for a given departure trajectory (Design 2). For illustration
purposes, transfers of the two types with Jacobi constant values
${J}_{d}=3.00754$ and ${J}_{a}=3.00240$, respectively at Ganymede and Europa,
and a departure epoch $\theta_{0_{Gan}}=50^{\circ}$ are here presented. Both
design types produce complete transfers dependent on the departure or arrival
trajectory, and the other portion of the route (arrival or departure,
respectively) is determined.

###  4.1 Design 1: Identifying departure conditions corresponding to a
specified arrival trajectory

Assume that the desired arrival trajectory is a close flyby of Europa with
approach and departure through the $L_{2}$ gateway, as depicted in Fig. 12.

Figure 12: Design 1: selected Europa arrival condition.

(a) Ganymede departure FTLE map.

(b) Total $\boldsymbol{\Delta{v}}$ budget.

(c) Total time–of–flight.

(d) Phase of Europa at departure.

Figure 13: Design 1: Ganymede departure FTLE map (a), $\Delta v$ magnitude
(b), time–of–flight (c) and Europa’s phase (d) inspection maps for the
Ganymede departure epoch $\boldsymbol{\theta_{0_{Gan}}=50^{\circ}}$.

The $\Delta{v}$ budget, the total time–of–flight and the relative phase
requirements between the moons (Fig. 13) are identified based on the
application of the MMAT method to the FTLE map corresponding to departures
from Ganymede. Equation (3) is evaluated for all possible transfers with
departure conditions from Ganymede matching the selected trajectory arc
approaching Europa. Cost and phase inspection maps are employed to obtain
cost–effective initial conditions, as apparent in Fig. 13, where the black
line corresponds to the unstable manifold and bounds departure options on the
maps. A departure trajectory with two revolutions around Ganymede (temporary
capture) is selected from among all these opportunities. Eventually, combining
the selected departure and arrival trajectories yields the full transfer
plotted in Fig. 16(a).

###  4.2 Design 2: Identifying arrival conditions corresponding to a specified
departure trajectory

In the selected departure trajectory, the spacecraft completes a few
revolutions around Ganymede before departure (temporary capture, see Fig. 14).
By leveraging the MMAT method, all the accessible arrival trajectories near
Europa and their associated costs ($\Delta v$) are determined. The initial
conditions for arrival are matched to the specified departure arc using Eq.
(3), and the inspection maps for $\Delta{v}$ budget, total time–of–flight and
relative phase requirements between moons (Fig. 15) are produced. For example,
the selected solution leads to a landing on the surface of Europa (see Fig.
16(b)).

Figure 14: Design 2: Selected Ganymede departure condition.

(a) Europa arrival availability.

(b) Total $\boldsymbol{\Delta{v}}$ budget.

(c) Total $\boldsymbol{t_{tot}}$ budget.

(d) Europa phase at departure epoch.

Figure 15: Design 2: Europa arrival FTLE map (a), $\Delta v$ magnitude (b),
time–of–flight (c) and Europa’s orbital phase (d) inspection maps for the
selected Ganymede departure epoch
($\boldsymbol{\theta_{0_{Gan}}=50^{\circ}}$).

(a) Design 1: from temporary capture around Ganymede to a flyby of Europa.

(b) Design 2: from temporary capture around Ganymede to landing at Europa.

Figure 16: Designs 1 (a) and 2 (b) for planning arrivals to Europa with
different endgames.

##  5 Dependence on the Jacobi constant

The methodology presented in the previous sections generates different results
depending on the Jacobi constant of the departure and arrival trajectories. If
the Jacobi constant value for the departure trajectories is decreased ($J_{d}$
= 3.0061, i.e., the departure energy is higher than in the previous example)
while the trajectories approaching Europa remain at $J_{a}$ = 3.00240, the
number of initial conditions leading to transfers from Ganymede to Europa
increases, as is apparent by comparing Fig. 10 with Fig. 17. Note that an
isoline of 1.25 is selected in Fig. 17 to illustrate the availability of a
larger number of initial conditions. Moreover, for the same $J_{a}$, the lower
the departure Jacobi constant ($J_{d}$) the lower the resulting $\Delta v$ for
most available sets of initial conditions and for the same departure epoch.
This fact is demonstrated by developing the inspection maps discussed in Sect.
4 for two departure Jacobi constant values and the same departure epoch.

Now, increase the Jacobi constant level for the arrival trajectories to
$J_{a}$ = 3.0030, i.e., the arrival energy is lower than in the previous two
cases, while the trajectories departing Ganymede remain at $J_{d}$ = 3.00754.
If the initial condition with the maximum value of $a_{a}(1+e_{a})/L_{a}^{*}$
is selected from the apoapsis arrival map, there are no transit orbits
departing the Ganymede vicinity that reach the vicinity of Europa (Fig. 18).
Hence, the selection of the departure and arrival Jacobi constant values is
crucial for the moon–to–moon design process.

Figure 17:  Available transfers connecting Ganymede and Europa in the
departure and arrival FTLE maps; the 1.25 isoline is used as a reference value
for selecting initial conditions. Figure 18:  Apoapsis arrival map (a) and
departure FTLE map (b) for which Ganymede and Europa cannot be connected for
the given combination of energy levels.

##  6  Dependence on the departure epoch

Moon–to–moon transfers are epoch–dependent, i.e., the number of connecting
trajectories available on the MMAT maps depends on the departure epoch. In the
previous examples, the epoch is represented as $\theta_{0_{Gan}}$ (true
anomaly of Ganymede at departure). In the following, two strategies for
selecting a departure epoch are discussed: one leverages the arrival
hyperbolic stable manifold and the moon–to–moon tides map, whereas the other
is based on the departure unstable manifold. For the sake of illustration, the
Ganymede–to–Europa scenario is continued.

###  6.1 Method A: arrival stable manifold and moon–to–moon tides map

The apoapsis and periapsis arrival maps respectively represent the apoapsis
and periapsis radii of the arrival conics originating from all the initial
conditions that cross the surface of the arrival moon SoI in the arrival FTLE
map. Transit orbits near the arrival moon are enclosed by the stable
hyperbolic invariant manifolds of planar Lyapunov orbits.

To locate the best departure epoch using this method, the stable manifold
trajectories of arrival are employed. For an inward journey, the stable
manifold trajectory that leads to the maximum
$\frac{a_{a}(1+e_{a})}{L_{a}^{*}}$ value is chosen. For the opposite journey,
a trajectory that provides the minimum value for
$\frac{a_{a}(1-e_{a})}{{L_{a}}^{*}}$ is selected. From this perspective, the
selected trajectory is termed the “minimum access arrival trajectory”, i.e.,
the first trajectory that grants access to the arrival moon. As noted, to
inspect the available connections from Ganymede to Europa (inward journey),
the minimum access arrival trajectory is the one that leads to the maximum
value of $\frac{a_{a}(1+e_{a})}{L_{a}^{*}}$ (first available trajectory to
Europa for the selected $J_{a}$). In this case, the minimum access arrival
trajectory located at isoline 1.2750 is plotted in Fig. 19 on the
corresponding apoapsis arrival map.

(a) Apoapsis arrival map with 1.2750 isoline overlaid.

(b) Stable manifold trajectory in the J–E CR3BP resulting in an
$\boldsymbol{L_{2}}$ Lyapunov orbit.

Figure 19: Minimum access arrival trajectory corresponding to the first stable
manifold with Ganymede access.

To determine whether transit orbits for a departure from Ganymede and arrival
near Europa exist, such an isoline is overlaid onto the moon–to–moon tides map
varying $\theta_{0_{Gan}}$. Figure 20 illustrates, for different values of
$\theta_{0_{Gan}}$, the trajectories departing the Ganymede vicinity that
first deliver access towards Europa.

Figure 20: Europa access variation in moon–to–moon tides map dependent on
$\theta_{0_{Gan}}$.

This procedure demonstrates the variety of possibilities within the range of
transfers based on departure epoch. At $\theta_{0_{Gan}}\approx 29^{\circ}$,
there is no access to Europa, though departure epochs greater than
approximately $35^{\circ}$ do yield transfer options. The number of
opportunities increases continuously until a maximum value corresponding to
$\theta_{0_{Gan}}\approx 82.5^{\circ}$. Then, the number decreases until
$\theta_{0_{Gan}}\approx 130^{\circ}$, where transfers are no longer possible.

###  6.2 Method B: departure unstable manifold

Though more precise, the moon–to–moon tides map approach requires many
computations for FTLE maps to discern the various permutations in the
departure epochs and identify the date with the largest amount of available
transfers. To simplify the process, an alternative strategy utilizing the
departure unstable manifold can be employed. The middle term in Eq. (3) is
computed for every unstable manifold trajectory propagated towards the SoI of
the departure moon over the span of departure epochs from $0^{\circ}$ to
$360^{\circ}$. The apoapsis and periapsis for the minimum access arrival
trajectory are then compared against the middle term in Eq. (3) for each
unstable manifold trajectory. Figure 21 represent such a comparison based on a
departure epoch of ${\theta_{0_{Gan}}=82.5^{\circ}}$: the “departure unstable
manifold angle” on the horizontal axis corresponds to the location of the
departure/arrival arc along the manifold associated with the periodic orbit
and measured clockwise from the $\hat{x}$–axis of the J–G rotating frame. The
results exhibit significant similarities to those produced from Method A.
Utilizing this technique for all $\theta_{0_{Gan}}$ values allows for
selection of the unstable manifold trajectories with minimum and maximum
values for the middle term in Eq. (3) (inward or outward transfer,
respectively). Finally, Fig. 22 illustrates the range of initial conditions
upon the arrival FTLE map from the first departure trajectory from Ganymede
that grants access to Europa, and the variations depending on
$\theta_{0_{Gan}}$. Employing the MMAT technique, a feasibility analysis
concerning the budget for $\Delta{v}$ and $t_{tot}$ may, thus, be accomplished
for this trajectory. By selecting the lowest $\Delta{v}$, the MMAT access maps
only need to be constructed for the given epoch.

Figure 21:  Evolution of Eq. 3 for the entire departure unstable manifold and
compared against the minimum access arrival trajectory.

##  7  Access dependence between two moons

The analysis of the MMAT access maps reveals a relationship concerning the
volume of available departure and arrival trajectories that link two moons.
For inward transfers, recall that the isolines in the MMAT maps are produced
by matching the values of
$\frac{{a}_{d}(1-{{e}_{d}^{2})}}{1+{e}_{d}\cos({\theta}_{{d}_{Int}})}$ and
$a_{a}(1+e_{a})$. For outward transfers, the quantity
$\frac{{a}_{d}(1-{{e}_{d}^{2})}}{1+{e}_{d}\cos({\theta}_{{d}_{Int}})}$ is
instead set equal to $a_{a}(1-e_{a})$. Consider again the Ganymede–to–Europa
case. Extensive experiments with the moon–to–moon tides map and the apoapsis
arrival map demonstrate that, for fixed Jacobi constant values at departure
and arrival, increasing the isoline value yields a higher number of options in
the moon–to–moon tides map, whereas the opposite occurs in the apoapsis
arrival map (see Fig. 23(b)). In contrast, decreasing the isoline value yields
more transfer options in the apoapsis arrival maps and fewer in the
moon–to–moon tides maps (Fig. 23(a)):

  1. 1.

An increase in the number of transit and unstable manifold trajectories at
departure reduces the amount of arrival options;

  2. 2.

An increase in the number of transit and stable manifold trajectories at
arrival reduces the available departure options.

In conclusion, a larger pool of trajectories at departure reduces the set of
possible arrival trajectories and vice versa. With the moons revolving in
their true orbital planes, the above effects are particularly important for
mission planning to deliver transfers with single impulsive maneuvers.

Figure 22:  Apoapsis arrival map depicting accessibility from the first
departure trajectory from Ganymede that grants access to Europa.

(a) Moon–to–moon tides map overlapped with apoapsis arrival map using 1.242
isoline.

(b) Moon–to–moon tides map overlapped with apoapsis arrival map using 1.26
isoline

Figure 23: MMAT maps application to understand access dependence by varying
isolines.

##  8  Outward versus inward transfers

The MMAT maps may also be employed in applications to deliver a spacecraft
from an inner to an outer moon. As an example, a transfer is designed from
Europa to Ganymede with the same Jacobi constant values as in the example of
Sect. 3, i.e., ${J}_{d}=3.00240$ (J–E CR3BP) and ${J}_{a}=3.00754$ (J–G
CR3BP). The transfer options are analyzed using the moon–to–moon tides map in
conjunction with the periapsis arrival map. The reference isoline employed in
this case is 0.779:

$\displaystyle\frac{{a}_{d}(1-{{e}_{d}^{2})}}{1+{e}_{d}\cos({\theta}_{{d}_{Int}})}={a}_{a}(1-{e}_{a})=0.779\;{L}_{a}^{*}.$
(10)

The matching of the two maps is illustrated in Fig. 24. The red arrows
indicate the initial conditions that lead to transfers from Europa to
Ganymede. The direction of the arrows is characteristic of the transfer
direction (inwards or outwards): for a path from an outer to an inner moon,
the arrows are oriented towards the right of the isoline in both FTLE maps
(i.e., see Fig. 8 whereas for transfers from an inner to an outer moon, the
arrows are directed to the left of the isoline (see Fig. 24).

Figure 24:  An overlay between moon–to–moon tides map and periapsis arrival
maps for an Europa to Ganymede transfer.

##  9  Conclusions

This contribution introduces the use of finite–time Lyapunov exponent (FTLE)
maps for the selection of desired motion patterns in the vicinity of the
departure and destination bodies within moon–to–moon transfers. The options
that can be identified include gravitational captures, lunar surface impacts
(or landings), transits through the vicinity of the moons, and insertion into
libration point periodic orbits. The FTLE maps have been blended with the
Moon–to–Moon Analytical Transfer (MMAT) method and applied to the design of
direct 3D single–impulse transfers between two Jovian moons, i.e., Ganymede
and Europa. The technique offers a very simple visualization of different
types of motion through so-called access maps. Their computation, constructed
using a combination of MATLAB and Java scripts, is relatively fast, taking
approximately 10 minutes on a Macbook Pro (2.3 GHz Intel Core i9, 16 GB RAM)
per map at a typical resolution of $10^{-4}$ in $y$ and $\dot{y}$ in
normalized CR3BP units. This performance can be improved by 25-50% if only the
calculation of the inner portion of the map (i.e., internal to the manifold)
is carried out.

The number of solutions, their type and performance (i.e., $\Delta v$ and time
of flight) depend on the launch date and the Jacobi constant of the departure
and destination orbits, hence the technique is suitable for trade off and
optimization studies. Moreover, if the spacecraft performs multiple
revolutions around the planet on the selected Keplerian orbits, more phasing
options become available between departure and destination bodies for the same
$\Delta v$.

Although the case analysed throughout the paper is a simple impulsive
connection between planet–centered ellipses originating from libration point
periodic orbits at the departure and arrival moons, it constitutes a solid
foundation for the development of the methodology. Firstly, the direct
transfer can provide the initial guess for the design of low-thrust
trajectories, in this way alleviating the mass budget corresponding to the
large associated $\Delta v$ values. Secondly, blending FTLE maps with the MMAT
technique has served the purpose of illustrating the method and maintaining
the focus of the discussion on its properties and benefits, rather than on the
inter–moon transfer design methodology. As a matter of fact, FTLE maps can
assist in the design of the begin game and end game of a moon–to–moon tour
constructed on intermediate gravity assists and resonant orbits: the
information provided by FTLE maps allows an assessment of the specific type of
trajectory connecting an initial condition in a low moon orbit with the first
resonance, a valuable possibility when shaping the science phase of a mission.
It is interesting to observe that among the different motion types illustrated
in this paper, takeoffs and landings are begin games and end games themselves
if the moon is in tidal lock with the planet (it is the case of the Galilean
moons as well as the Inner Large Moons of Saturn) because, in this case, the
positions of the launch and landing sites in the planet–moon synodic
barycentric reference frame can be straightforwardly mapped to geographical
coordinates. In conclusion, the proposed FTLE maps are useful to characterize
the space around the moon and study the available options within the envelope
of a loose capture, which is the intermediate stage of any low–energy transfer
involving approach and escape.

## Acknowledgments

Assistance from colleagues in the Multi–Body Dynamics Research group at Purdue
University is acknowledged as is the support from the Purdue University School
of Aeronautics and Astronautics and College of Engineering including access to
the Rune and Barbara Eliasen Visualization Laboratory. D. Canales is currently
supported by the Department of Aerospace Engineering in Embry–Riddle
Aeronautical University. E. Fantino acknowledges Khalifa University’s internal
grant 8474000413/CIRA–2021–065 and projects PID2020–112576GB–C21 and
PID2021–123968NB–100 awarded by the Spanish Ministry of Science and
Innovation.

## References

  * National Academies of Sciences, Engineering, and Medicine [2022] National Academies of Sciences, Engineering, and Medicine, _Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023-2032_ , The National Academies Press, Washington, DC, 2022.  10.17226/26522.
  * Phillips and Pappalardo [2014] Phillips, C. B., and Pappalardo, R. T., “Europa Clipper Mission Concept,” _Eos, Transactions, American Geophysical Union_ , Vol. 95, No. 20, 2014, pp. 165–167.  10.1002/2014EO200002.
  * ESA [2014] ESA, _JUICE - JUpiter ICy moons Explorer, Vol. 1_ , ESA, 2014\.  10.1016/j.pss.2012.12.002.
  * Sims [2006] Sims, J. A., “Jupiter Icy Moons Orbiter mission design overview,” _Proceedings of the AAS/AIAA SpaceFlight Mechanics Meeting_ , Tampa, Florida, 2006.
  * Ross et al. [2003] Ross, S. D., Koon, W. S., Lo, M. W., and Marsden, J. E., “Design of a multi-moon orbiter,” _Proceedings of the 13th AAS/AIAA Space Flight Mechanics Meeting_ , Ponce, Puerto Rico, 2003.  Paper AAS 03-143.
  * Izzo et al. [2013] Izzo, D., Simoes, L. F., Martens, M., de Croon, G. C. H. E., Heritier, A., and Yam, C., “Search for a grand tour of the Jupiter Galilean moons,” _Proceedings of the 15th Annual Conference on Genetic and Evolutionary Computation_ , 2013, pp. 1301–1308.  10.1145/2463372.2463524.
  * Colasurdo et al. [2014] Colasurdo, G., Zavoli, A., Longo, A., Casalino, L., and Simeoni, F., “Tour of Jupiter Galilean moons: Winning solution of GTOC6,” _Acta Astronautica_ , Vol. 102, 2014, pp. 190–199.  10.1016/j.actaastro.2014.06.003.
  * Gómez et al. [2001] Gómez, G., Koon, W. S., Lo, M. W., Marsden, J. E., Masdemont, J. J., and Ross, S. D., “Invariant manifolds, the spatial three-body problem and space mission design,” _Advances in the Astronautical Sciences_ , Vol. 109, 2001, pp. 3–22.  10.1142/9789812704849_0025.
  * Gómez et al. [2004] Gómez, G., Koon, W. S., Lo, M. W., Marsden, J. E., Masdemont, J. J., and Ross, S. D., “Invariant manifolds, the spatial three-body problem and space mission design,” _Nonlinearity_ , Vol. 17, No. 5, 2004, p. 1571.  10.1142/9789812704849_0025.
  * Koon et al. [2000] Koon, W. S., Lo, M. W., Marsden, J. E., and Ross, S. D., “Heteroclinic connections between periodic orbits and resonance transitions in celestial mechanics,” _Chaos_ , Vol. 10, No. 2, 2000, pp. 427–469.  10.1063/1.166509.
  * Koon et al. [2002] Koon, W. S., Lo, M. W., Marsden, J. E., and Ross, S. D., “Constructing a low energy transfer between Jovian moons,” _Celestial Mechanics, Dedicated to Donald Saari for his 60th Birthday_ , 2002, p. 129.  10.1090/conm/292/04919.
  * Koon, W. S. and Lo, M. W. and Marsden, J. E. and Ross, S. D. [2011] Koon, W. S. and Lo, M. W. and Marsden, J. E. and Ross, S. D., _Dynamical Systems, the Three-Body Problem and Space Mission Design_ , Marsden Books, 2011.  10.1142/9789812792617_0222.
  * Grover and Ross [2009] Grover, P., and Ross, S., “Designing Trajectories in a Planet-Moon Environment Using the Controlled Keplerian Map,” _Journal of Guidance, Control, and Dynamics_ , Vol. 32, No. 2, 2009, pp. 437–444.  10.2514/1.38320.
  * Campagnola and Russell [2010] Campagnola, S., and Russell, R. P., “Endgame Problem Part 1: $V_{\infty}$-Leveraging Technique and the Leveraging Graph,” _Journal of Guidance, Control, and Dynamics_ , Vol. 33, No. 2, 2010, pp. 463–475.  https://doi.org/10.2514/1.44258.
  * Campagnola and Russell [2010] Campagnola, S., and Russell, R. P., “Endgame Problem Part 2: Multibody Technique and the Tisserand-Poincaré Graph,” _Journal of Guidance, Control, and Dynamics_ , Vol. 33, No. 2, 2010, pp. 476–486.  10.2514/1.44290.
  * Lantoine et al. [2011] Lantoine, G., Russell, R. P., and Campagnola, S., “Optimization of low-energy resonant hopping transfers between planetary moons,” _Acta Astronautica_ , Vol. 68, No. 7, 2011, pp. 1361–1378.  10.1016/j.actaastro.2010.09.021.
  * Lantoine and Russell [2011] Lantoine, G., and Russell, R. P., “Near Ballistic Halo-to-Halo Transfers between Planetary Moons,” _Journal of the Astronautical Sciences_ , Vol. 58, No. 3, 2011, pp. 335–363.  10.1007/BF03321174.
  * Fantino and Castelli [2017] Fantino, E., and Castelli, R., “Efficient design of direct low-energy transfers in multi-moon systems,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 127, 2017, pp. 429–450.  10.1007/s10569-016-9733-9.
  * Fantino et al. [2018] Fantino, E., Flores, R. M., and Al-Khateeb, A. N., “Efficient Two-Body Approximations of Impulsive Transfers between Halo Orbits,” _Proceedings of the 69th International Astronautical Congress_ , Bremen, Germany, 2018.
  * Canales et al. [2020] Canales, D., Howell, K. C., and Fantino, E., “Moon-to-Moon Transfer Methodology for Multi-Moon Systems in the Coupled Spatial Circular Restricted Three-Body Problem,” _Proceedings of the AAS/AIAA Astrodynamics Specialist Conference_ , South Lake Tahoe, California, 2020.
  * Canales et al. [2021a] Canales, D., Howell, K. C., and Fantino, E., “Transfer design between neighborhoods of planetary moons in the circular restricted three-body problem: the moon-to-moon analytical transfer method,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 133, No. 8, 2021a, p. 36.  10.1007/s10569-021-10031-x.
  * Canales et al. [2022a] Canales, D., Howell, K. C., and Fantino, E., “A versatile moon-to-moon transfer design method for applications involving libration point orbits,” _Acta Astronautica_ , Vol. 198, 2022a, pp. 388–402.  10.1016/j.actaastro.2022.06.010.
  * Canales et al. [2021b] Canales, D., Gupta, M., Park, B., and Howell, K. C., “Exploration of Deimos and Phobos Leveraging Resonant Orbits,” _Proceedings of the 31st AAS/AIAA Space Flight Mechanics Meeting_ , Charlotte, North Carolina, 2021b.
  * Canales et al. [2022b] Canales, D., Gupta, M., Park, B., and Howell, K. C., “A transfer trajectory framework for the exploration of Phobos and Deimos leveraging resonant orbits,” _Acta Astronautica_ , Vol. 194, 2022b, pp. 263–276.  10.1016/j.actaastro.2022.02.001.
  * Lara et al. [2007] Lara, M., Russell, R. P., and Villac, B. F., “Fast estimation of stable regions in real models,” _Meccanica_ , Vol. 42, 2007, pp. 511–515.  https://doi.org/ 10.1007/s11012-007-9060-z.
  * Villac [2008] Villac, B. F., “Using FLI maps for preliminary spacecraft trajectory design in multi-body environments,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 102, 2008, pp. 29–48.  https://doi.org/10.1007/s10569-008-9158-1.
  * Haller and Sapsis [2011] Haller, G., and Sapsis, T., “Lagrangian coherent structures and the smallest finite-time Lyapunov exponent,” _Chaos_ , Vol. 21, 2011, p. 023115\.  10.1063/1.3579597.
  * Canales et al. [2021c] Canales, D., Howell, K. C., and Fantino, E., “Using Finite-Time Lyapunov Exponent maps for planetary moon-tour design,” _Proceedings of the 2021 AAS/AIAA Astrodynamics Specialist Conference_ , Big Sky, Virtual, 2021c.  Paper AAS 21-625.
  * Canales et al. [2021d] Canales, D., Howell, K. C., and Fantino, E., “Using Finite-Time Lyapunov Exponent maps for planetary moon-tour design,” _Proceedings of ASCEND 2021_ , Las Vegas, Nevada & Virtual, 2021d.  10.2514/6.2021-4154.
  * Poincaré [1892] Poincaré, H., _Les méthodes nouvelles de la mécanique céleste_ , Gauthier-Villars et Fils, Paris, France, 1892.  10.1007/BF02742713.
  * Szebehely [1974] Szebehely, V., _The General and Restricted Problems of Three Bodies_ , Springer, Vienna, 1974.  10.1007/978-3-7091-2916-6.
  * Broucke [1968] Broucke, R. A., “Periodic Orbits in the Restricted Three-body Problem with Earth-Moon Masses,” Tech. rep., Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, 1968.  TR 32-1168.
  * Howell and Campbell [1999] Howell, K. C., and Campbell, E. T., “Three-dimensional periodic solutions that bifurcate from halo families in the circular restricted three-body problem,” _Advances in the Astronautical Sciences_ , Vol. 102, 1999, pp. 891–910.
  * Haapala and Howell [2015] Haapala, A., and Howell, K. C., “A Framework for Constructing Transfers Linking Periodic Libration Point Orbits in the Spatial Circular Restricted Three-Body Problem,” _International Journal of Bifurcation and Chaos_ , Vol. 26, 2015, p. 1630013.  10.1142/S0218127416300135.
  * Wiggins and Guckenheimer [1992] Wiggins, S., and Guckenheimer, J., “Chaotic Transport in Dynamic Systems,” _Physics Today_ , Vol. 45, No. 7, 1992, pp. 68–69.  10.1063/1.2809741.
  * Gawlik et al. [2009] Gawlik, E. S., Marsden, J. E., Du Toit, P. C., and Campagnola, S., “Lagrangian coherent structures in the planar elliptic restricted three-body problem,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 103, No. 3, 2009, pp. 227–249.  10.1007/s10569-008-9180-3.
  * Short and Howell [2014] Short, C., and Howell, K. C., “Lagrangian coherent structures in various map representations for application to multi-body gravitational regimes,” _Acta Astronautica_ , Vol. 94, 2014, pp. 592–607.  10.1016/j.actaastro.2013.08.020.
  * Perez-Palau et al. [2012] Perez-Palau, D., Gómez, G., and Masdemont, J., “Detecting Invariant Manifolds Using Hyperbolic Lagrangian Coherent Structures,” _Proceedings of the 1st IAA conference on dynamics and control of spacecrafts_ , Porto, Portugal, 2012.
  * Haller [2011] Haller, G., “A variational theory of hyperbolic Lagrangian Coherent Structures,” _Physica D: Nonlinear Phenomena_ , Vol. 240, No. 7, 2011, pp. 574–598.  10.1016/j.physd.2010.11.010.
  * CGS [1997] “1 - Kinematics of Fluid Flow,” _Fluid Mechanics of Viscoelasticity_ , Rheology Series, Vol. 6, edited by R. Huilgol and N. Phan-Thien, Elsevier, 1997, pp. 1–83.  https://doi.org/10.1016/S0169-3107(97)80003-2.
  * Golub and Kahan [1965] Golub, G., and Kahan, W., “Calculating the Singular Values and Pseudo-Inverse of a Matrix,” _Journal of the Society for Industrial and Applied Mathematics Series B Numerical Analysis_ , Vol. 2, No. 2, 1965, pp. 205–224.  10.1137/0702016, URL https://doi.org/10.1137/0702016.
  * Smith [1993] Smith, D. R., _An Introduction to Continuum Mechanics - after Truesdell and Noll_ , Springer, Dordrecht, The Netherlands, 1993.  10.1007/978-94-017-0713-8.
  * Short et al. [2015] Short, C. R., Blazevski, D., Howell, K. C., and Haller, G., “Stretching in phase space and applications in general non-autonomous multi-body problems,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 122, No. 3, 2015, pp. 213–238.  10.1007/s10569-015-9617-4.
  * SPI [2023] “SPICE - An Observation Geometry System for Space Science Missions,” https://naif.jpl.nasa.gov/naif, 2023.  Accessed: 2023-05-12.
  * Canales [2021] Canales, D., “Transfer design methodology between neighborhoods of planetary moons in the circular restricted three-body problem,” Ph.D. thesis, Purdue University, December 2021.