Computation and Analysis of Jupiter-Europa and Jupiter-Ganymede Resonant Orbits in the Planar Concentric Circular Restricted 4-Body Problem

We know that the $m_{1}$-$m_{3}$ frame $x$-axis is the line between those two bodies, which is represented in the $m_{1}$-$m_{2}$ coordinate frame as a line passing through the point $(-\mu,0)$, making an angle $\theta_{3}$ with the $m_{1}$-$m_{2}$ frame $x$-axis, and revolving at an angular rate of $\Omega_{3}-1$. Thus, to transform a state from the $m_{1}$-$m_{2}$ frame into the $m_{1}$-$m_{3}$ frame described earlier, one must follow several steps:

1. 1.

   Shift the origin of the position from the $m_{1}$-$m_{2}$ barycenter to $m_{1}$

2. 2.

   Rotate the $m_{1}$-$m_{2}$ frame position and velocity vectors by the angle $-\theta_{3}$

3. 3.

   Apply the transport theorem to get the actual $m_{1}$-$m_{3}$ frame apparent velocity

4. 4.

   Rescale the length, mass, and time units, and set $\theta_{2}=-\theta_{3}$

5. 5.

   Shift the origin of the position from $m_{1}$ to the $m_{1}$-$m_{3}$ barycenter

Let $(x,y,p_{x},p_{y},\theta_{3})$ be a state in the $m_{1}$-$m_{2}$ frame, and let $(\bar{x},\bar{y},\bar{p}_{x},\bar{p}_{y},\theta_{2})$ denote the equivalent state in the $m_{1}$-$m_{3}$ frame. Denote $R_{\rho}$ as the $2\times 2$ rotation matrix such that for $\mathbb{x}\in\mathbb{R}^{2}$, the vector $R_{\rho}\mathbb{x}$ is rotated by $\rho$ radians counterclockwise. Then, remembering that $\dot{x}=p_{x}+y$ and $\dot{y}=p_{y}-x$, steps 1 and 2 from above for finding $(\bar{x},\bar{y},\bar{p}_{x},\bar{p}_{y})$ require computing

$$\begin{bmatrix}x^{\prime}\\ y^{\prime}\end{bmatrix}=R_{-\theta_{3}}\begin{bmatrix}x+\mu\\ y\end{bmatrix}\quad\quad\begin{bmatrix}\dot{x}^{\prime}\\ \dot{y}^{\prime}\end{bmatrix}=R_{-\theta_{3}}\begin{bmatrix}\dot{x}\\ \dot{y}\end{bmatrix}$$

(5)

After this, we can write expressions for $(\bar{x},\bar{y},\dot{\bar{x}},\dot{\bar{y}})$ in the $m_{1}$-$m_{3}$ frame as

$$\bar{x}=\frac{x^{\prime}}{r_{13}}-\bar{\mu}\quad\quad\bar{y}=\frac{y^{\prime}}{r_{13}}\quad\quad\dot{\bar{x}}=\frac{\dot{x}^{\prime}+(\Omega_{3}-1)y^{\prime}}{r_{13}\Omega_{3}}\quad\quad\dot{\bar{y}}=\frac{\dot{y}^{\prime}-(\Omega_{3}-1)x^{\prime}}{r_{13}\Omega_{3}}$$

(6)

These expressions encompass steps 3, 4, and 5 from the previous list. Last of all, we can compute $\bar{p}_{x}=\dot{\bar{x}}-\bar{y}$ and $\bar{p}_{y}=\dot{\bar{y}}+\bar{x}$; substituting Eq. (6) into these expressions for $\bar{p}_{x}$ and $\bar{p}_{y}$ actually gives an even simpler alternative to Eq. (6). In particular, we find that

$$\bar{p}_{x}=\frac{\dot{x}^{\prime}-y^{\prime}}{r_{13}\Omega_{3}}=\frac{p_{x}^{\prime}}{r_{13}\Omega_{3}}\quad\quad\quad\bar{p}_{y}=\frac{\dot{y}^{\prime}+x^{\prime}}{r_{13}\Omega_{3}}-\bar{\mu}=\frac{p_{y}^{\prime}}{r_{13}\Omega_{3}}-\bar{\mu}$$

(7)

where $\begin{bmatrix}p_{x}^{\prime}&p_{y}^{\prime}\end{bmatrix}^{T}=R_{-\theta_{3}}\begin{bmatrix}p_{x}&p_{y}+\mu\end{bmatrix}^{T}$. Finally, once $(\bar{x},\bar{y},\bar{p}_{x},\bar{p}_{y})$ has been found, one should also keep in mind that $\theta_{2}=-\theta_{3}$, which completes the full extended state required to propagate the $m_{1}$-$m_{3}$ equations of motion given in Eq. (4). We will refer to this transformation as $\Phi_{\theta_{3}}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$; that is, using the earlier notation of this section, $(\bar{x},\bar{y},\bar{p}_{x},\bar{p}_{y})=\Phi_{\theta_{3}}(x,y,p_{x},p_{y})$. Note that $\Phi_{\theta_{3}}$ is just a composition of affine and linear transformations, so its derivative $D\Phi_{\theta_{3}}$ is simple to compute.

The quasi-periodic orbits of interest in the planar CCR4BP lie on 2D unstable invariant tori in the 5D extended phase space $(x,y,p_{x},p_{y},\theta_{p})$, where $p=3$ or $2$ depending on the frame being used. These invariant tori can be parameterized as the image of a function of two angles $K_{2}:\mathbb{T}^{2}\rightarrow\mathbb{R}^{4}\times\mathbb{T}$. Any quasi-periodic trajectory $\mathbb{x}(t)$ lying on this torus can be expressed as

$$\mathbb{x}(t)=K_{2}(\theta,\theta_{p})\quad\quad\quad\theta=\theta_{0}+\Omega_{1}t,\quad\theta_{p}=\theta_{p,0}+(\Omega_{p}-1)t$$

(8)

where $\mathbb{x}(0)$ determines $\theta_{0}$ and $\theta_{p,0}$. $\theta_{p}$ and $\Omega_{p}-1$ are the perturbation phase angle and frequency described earlier, respectively. Defining the stroboscopic map $F:\mathbb{R}^{4}\times\mathbb{T}\rightarrow\mathbb{R}^{4}\times\mathbb{T}$ as the time-$\frac{2\pi}{|\Omega_{p}-1|}$ mapping of extended phase space points by the equations of motion, we have

$$F(K_{2}(\theta,\theta_{p}))=K_{2}(\theta+\omega,\theta_{p}),\text{ where }\omega=2\pi\Omega_{1}/|\Omega_{p}-1|$$

(9)

since the angle $\theta_{p}$ increases or decreases by $2\pi$ in the time $\frac{2\pi}{|\Omega_{p}-1|}$. Since the value of $\theta_{p}$ is invariant under the map $F$, we can fix $\theta_{p}$ and define $K(\theta)=K_{2}(\theta,\theta_{p})$. Without loss of generality, we choose $\theta_{p}=0$ in this study; this has the added benefit that both $\theta_{2}=\theta_{3}=0$, so the $m_{1}$-$m_{2}$ and $m_{1}$-$m_{3}$ frames are aligned with each other at this phase. With $K$ thus defined, Equation (9) becomes

$$F(K(\theta))=K(\theta+\omega)$$

(10)

Ignoring the constant $\theta_{p}$ component of the extended phase space and making a slight abuse of notation, we have $F:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ and $K:\mathbb{T}\rightarrow\mathbb{R}^{4}$. Equation (10) implies that $K$ is an invariant 1D torus (circle) of $F$. Hence, basing our computations on the stroboscopic map $F$ is more efficient than solving for tori invariant under the flow of the equations of motion, since we reduce the phase space dimension from 5D to 4D and the dimension of the unknown invariant tori from 2D to 1D. Once the 1D stroboscopic map torus is computed, one can visualize the full 2D torus of the flow by numerically integrating points from the invariant circle by the equations of motion Eq. (2) or (4).

With the stroboscopic map $F$ defined by the CCR4BP equations of motion, we now wish to find solutions $K(\theta)$ of Equation (10). In this section, we summarize the results of the parameterization methods developed in Kumar, Anderson, and de la Llave[11] for computing invariant tori in periodically-perturbed planar CRTBP models such as the planar CCR4BP; these methods were themselves inspired by methods described in Haro et al. [13] and Zhang and de la Llave[14]. The rotation number $\omega=\frac{2\pi\Omega_{1}}{|\Omega_{p}-1|}$ is generally known; for instance, if the CCR4BP torus being solved for comes from a known PCRTBP periodic orbit, then $\Omega_{1}$ is the frequency of that periodic orbit.

In our previous work[11] , we developed an efficient quasi-Newton method for solving equation (10) given a sufficiently accurate initial guess. Our quasi-Newton method adds an extra equation to be solved in addition to Eq. (10). In particular, as well as $K(\theta)$, we simultaneously solve for matrix-valued periodic functions $P(\theta)$, $\Lambda(\theta):\mathbb{T}\rightarrow\mathbb{R}^{4\times 4}$ satisfying

$$DF(K(\theta))P(\theta)=P(\theta+\omega)\Lambda(\theta)$$

(11)

$P(\theta)$ and $\Lambda(\theta)$ are the matrices of bundles and of Floquet stability, respectively. For each $\theta\in\mathbb{T}$, the columns of $P(\theta)$ are comprised of the (linearly independent) tangent, symplectic conjugate center, stable, and unstable directions of the torus at the point $K(\theta)$, in that order, while $\Lambda$ has the form

$$\Lambda(\theta)=\begin{bmatrix}1&T(\theta)&0&0\\ 0&1&0&0\\ 0&0&\lambda_{s}&0\\ 0&0&0&\lambda_{u}\end{bmatrix}$$

(12)

where $T:\mathbb{T}\rightarrow\mathbb{R}$ and $\lambda_{s},\lambda_{u}\in\mathbb{R}$ are constants with $\lambda_{s}<1$ and $\lambda_{u}>1$. $\lambda_{s}$ and $\lambda_{u}$ are the stable and unstable multipliers for the invariant torus, which we will also find useful in this study.

As it turns out, solving simultaneously for $K$, $P$, and $\Lambda$ not only gives stability information, but actually has lower computational complexity than more commonly-used methods which solve for $K$ alone. When discretizing functions of $\theta$ at $N$ evenly spaced angle values, our method requires only $O(N)$ storage and $O(N\log N)$ operations since the nearly-diagonal $\Lambda$ matrix allows decoupling of the scalar functional equations which are solved in each quasi-Newton step[11]. This is much more efficient than $K$-only methods which require solving an un-decoupled system of functional equations in each step, which after discretization requires $O(N^{3})$ operations due to Gaussian elimination.

In this work, we use the quasi-Newton method of the previous section as part of a numerical continuation scheme for computing tori. Note that the stroboscopic map $F$ defined earlier is obtained by integrating points by the equations of motion Eq. (2) or (4); hence, $F$ depends on the value of $\mu_{3}$ or $\bar{\mu}_{2}$, depending on the frame being used. Due to this, in this section we write $F=F_{\mu_{3}}$ or $F=F_{\bar{\mu}_{2}}$ to signify this dependence on parameters which we will vary during the continuation.

To compute an invariant torus for some desired value of $\mu_{3}$, we can start with a periodic orbit in the $m_{1}$-$m_{2}$ PCRTBP, and then continue this by $\mu_{3}$ until the required $\mu_{3}$ value is reached. A periodic orbit is an invariant circle for $F_{\mu_{3}=0}$, so we can take this periodic orbit as a starting point for the continuation by $\mu_{3}$. We can similarly continue $m_{1}$-$m_{3}$ PCRTBP periodic orbits by $\bar{\mu}_{2}$. One thing to note is that due to Kepler’s third law (see Eq. (1)), fixing the angular velocity $\Omega_{3}$ of $m_{3}$ and varying $\mu_{3}$ during the continuation results in $r_{13}$ changing as well; similarly, $r_{12}$ changes during the continuation by $\bar{\mu}_{2}$ in the case of an $m_{1}$-$m_{3}$ frame continuation. It is desirable to fix $\Omega_{3}$ or $\Omega_{2}$ to its physical value during the entire continuation, so that the invariant circle rotation number $\omega=\frac{2\pi\Omega_{1}}{|\Omega_{3}-1|}$ or $\frac{2\pi\Omega_{1}}{|\Omega_{2}-1|}$ remains constant. Once a single torus has been continued by $\mu_{3}$ or $\bar{\mu}_{2}$ to the desired physical mass value, the same quasi-Newton method can also be used to continue the resulting torus by $\omega$ (with fixed mass values) to find other CCR4BP tori in the same family.

While the quasi-Newton method gives a way to correct an approximate solution of Eq. (10)-(11), we improved continuation performance by adding a predictor step based on the Poincaré-Lindstedt method to complement the corrector. We described the Lindstedt method for $\omega$ continuation in our previous work[11]. Here we derive a similar method for predicting the next torus during a continuation by $\mu_{3}$ (the predictor for $\bar{\mu}_{2}$ continuation works the exact same way).

Assume we have a solution $K_{\mu_{3}}(\theta),P_{\mu_{3}}(\theta),\Lambda_{\mu_{3}}(\theta)$ to Eq.(10)-(11) for some value of $\mu_{3}$. Then, we have that $F_{\mu_{3}}(K_{\mu_{3}}(\theta))=K_{\mu_{3}}(\theta+\omega)$. Differentiating this with respect to $\mu_{3}$, we have

$$\frac{dF_{\mu_{3}}}{d\mu_{3}}(K_{\mu_{3}}(\theta))+DF_{\mu_{3}}(K_{\mu_{3}}(\theta))\frac{dK_{\mu_{3}}}{d\mu_{3}}(\theta)=\frac{dK_{\mu_{3}}}{d\mu_{3}}(\theta+\omega)$$

(13)

where $DF_{\mu_{3}}$ is the matrix-valued derivative of $F_{\mu_{3}}$ with respect to the state $(x,y,p_{x},p_{y})$, and $\frac{dF_{\mu_{3}}}{d\mu_{3}}$ is the vector-valued derivative of $F_{\mu_{3}}$ with respect to the parameter $\mu_{3}$. To solve for $\frac{dK_{\mu_{3}}}{d\mu_{3}}(\theta)$, substitute $\frac{dK_{\mu_{3}}}{d\mu_{3}}(\theta)=P_{\mu_{3}}(\theta)\xi(\theta)$ into Eq. (13). Left multiplying by $P^{-1}_{\mu_{3}}(\theta+\omega)$ and recalling from Eq. (11) that $P_{\mu_{3}}^{-1}(\theta+\omega)DF_{\mu_{3}}(K_{\mu_{3}}(\theta))P_{\mu_{3}}(\theta)=\Lambda_{\mu_{3}}(\theta)$, we get the equation

$$\Lambda_{\mu_{3}}(\theta)\xi(\theta)-\xi(\theta+\omega)=-P_{\mu_{3}}^{-1}(\theta+\omega)\frac{dF_{\mu_{3}}}{d\mu_{3}}(K_{\mu_{3}}(\theta))$$

(14)

Recalling that $F_{\mu_{3}}$ just represents the integration of an ODE by some fixed time, $\frac{dF_{\mu_{3}}}{d\mu_{3}}(K_{\mu_{3}}(\theta))$ can be computed using the usual variational equations for a parameter dependent system. Hence, the RHS of Eq. (14) is known; equations of the form of Eq. (14) are almost diagonal due to Eq. (12), and $\xi$ can easily be solved for using Fourier series as is described in our previous work[11]. Once $\xi$ and hence $\frac{dK_{\mu_{3}}}{d\mu_{3}}(\theta)$ are known, we can use this to predict the torus for $\mu_{3}+\Delta\mu_{3}$ as

$$K_{\mu_{3}+\Delta\mu_{3}}(\theta)\approx K_{\mu_{3}}(\theta)+\Delta\mu_{3}\frac{dK_{\mu_{3}}}{d\mu_{3}}(\theta)$$

(15)

The result of this computation is then used as the initial guess (along with $P_{\mu_{3}}(\theta)$ and $\Lambda_{\mu_{3}}(\theta)$) for the quasi-Newton method for finding an invariant circle of $F_{\mu_{3}+\Delta\mu_{3}}$.

Assume we have a solution $K,P,\Lambda$ of Eq. (10)-(11) with $F=F_{\mu_{3}}$, the stroboscopic map for some fixed $\theta_{3}$ of the $m_{1}$-$m_{2}$ frame equations of motion Eq. (2); we would like to transform this to the corresponding solution $\bar{K},\bar{P},\bar{\Lambda}$ of Eq. (10)-(11) for $F=F_{\bar{\mu}_{2}}$, the $\theta_{2}=-\theta_{3}$ based stroboscopic map of the $m_{1}$-$m_{3}$ frame equations of motion Eq. (4). Recall the function $\Phi_{\theta_{3}}$ defined earlier taking states from the $m_{1}$-$m_{2}$ frame to the $m_{1}$-$m_{3}$ frame. The relationship between $F_{\mu_{3}}$, $F_{\bar{\mu}_{2}}$, and $\Phi_{\theta_{3}}$ is

$$F_{\bar{\mu}_{2}}\circ\Phi_{\theta_{3}}=\Phi_{\theta_{3}}\circ F_{\mu_{3}}$$

(16)

Eq. (16) simply means that the same final state is attained whether one first transforms to the $m_{1}$-$m_{3}$ frame and then applies the $m_{1}$-$m_{3}$ frame stroboscopic map $F_{\bar{\mu}_{2}}$, or one first applies the $m_{1}$-$m_{2}$ frame stroboscopic map $F_{\mu_{3}}$ and then only afterwards transforms to the $m_{1}$-$m_{3}$ frame. This can be justified by noting that substituting $(x,y,p_{x},p_{y})=\Phi_{\theta_{3}}^{-1}(\bar{x},\bar{y},\bar{p}_{x},\bar{p}_{y})$ and $\theta_{3}=-\theta_{2}$ into Eq. (2), followed by a rescaling of time to $m_{1}$-$m_{3}$ frame time units $\Omega_{3}\,dt=d\bar{t}$, yields Eq. (4); the derivation of this is simple but lengthy, so we omit it for the purpose of space. Eq. (16) is simply the discrete-time version of the equivalency of Eq. (2) and Eq. (4) through $\Phi_{\theta_{3}}$.

Since we know that $F_{\mu_{3}}(K(\theta))=K(\theta+\omega)$, Eq. (16) implies that

$$F_{\bar{\mu}_{2}}\left(\Phi_{\theta_{3}}(K(\theta))\right)=\Phi_{\theta_{3}}\left(F_{\mu_{3}}(K(\theta))\right)=\Phi_{\theta_{3}}(K(\theta+\omega))$$

(17)

which means that $\bar{K}(\theta)=\Phi_{\theta_{3}}(K(\theta))$ satisfies Eq. (10) for $F=F_{\bar{\mu}_{2}}$. Similarly, we know that $DF_{\mu_{3}}(K(\theta))P(\theta)=P(\theta+\omega)\Lambda(\theta)$; differentiating Eq. (16) and right-multiplying by $P$ thus gives

$$\begin{split}DF_{\bar{\mu}_{2}}\left(\Phi_{\theta_{3}}(K(\theta))\right)D\Phi_{\theta_{3}}(K(\theta))P(\theta)&=D\Phi_{\theta_{3}}\left(F_{\mu_{3}}(K(\theta))\right)DF_{\mu_{3}}(K(\theta))P(\theta)\\ &=D\Phi_{\theta_{3}}\left(K(\theta+\omega)\right)P(\theta+\omega)\Lambda(\theta)\end{split}$$

(18)

which implies that $\bar{K}(\theta)=\Phi_{\theta_{3}}(K(\theta))$, $\bar{P}(\theta)=D\Phi_{\theta_{3}}(K(\theta))P(\theta)$, and $\bar{\Lambda}(\theta)=\Lambda(\theta)$ satisfies Eq. (11) for $F=F_{\bar{\mu}_{2}}$. Thus, through these equations, we can transform any solution of Eq. (10)-(11) for $F=F_{\mu_{3}}$ to the corresponding solution for $F=F_{\bar{\mu}_{2}}$.

We first compute the quasi-periodic equivalents of Jupiter-Europa 3:4 unstable resonant periodic orbits in the CCR4BP. Starting from a periodic orbit with Jacobi constant value 3.0041 and $\omega=3.097849$ (arbitrarily chosen), the Jupiter-Europa frame planar CCR4BP stroboscopic map invariant circles computed at various $\mu_{3}$ continuation steps are displayed in Figure 1; the continuation step size used was $\Delta\mu_{3}=8\times 10^{-6}$ until the final value of $\mu_{3}=7.804102777055038\times 10^{-5}$ for Ganymede. We also show the orbit of Ganymede on the plot as a red dashed-line circle, although recall that our stroboscopic map is defined using a phase of $\theta_{3}=0$, so Ganymede is on the $x$-axis aligned with Europa at the phase of these computed invariant circles.

Zooming into the neighborhood of Europa in the $\mu_{3}$ continuation plot, it is visible that the invariant circle moves closer to Europa as $\mu_{3}$ increases, potentially to offset some of the gravitational pull of Ganymede in the other direction. In physical units, the closest approach to Europa decreases from 22052 km to 18721 km. Nevertheless, plotting the unstable Floquet multiplier $\lambda_{u}$ (recall Eq. (12)) of the torus versus the parameter $\mu_{3}$ in Figure 2, we can see that the instability of the torus actually decreases slightly as $\mu_{3}$ increases and the invariant circle moves towards Europa.

After computing a single Jupiter-Europa 3:4 unstable resonant CCR4BP torus for the physical value of Ganymede’s mass ratio $\mu_{3}$, we continued this torus by $\omega$ to get a family of invariant circles for the stroboscopic map. Figure 3 shows the resulting tori computed in the Jupiter-Europa frame; the overall structure of the invariant circle family is similar to that of the Jupiter-Europa 3:4 unstable resonant PCRTBP periodic orbit family.

With Jupiter-Europa 3:4 unstable resonant CCR4BP invariant circles computed for the Jupiter-Europa frame stroboscopic map, we can transform the invariant circle states to the Jupiter-Ganymede rotating frame using the function $\Phi_{\theta_{3}}$ defined previously; note that one of the benefits of using $\theta_{3}=0$ as the stroboscopic map phase is that no rotation is required during this transformation between frames. As demonstrated earlier, the resulting set of states will be an invariant circle of the Jupiter-Ganymede frame CCR4BP stroboscopic map based at $\theta_{2}=0$. One such invariant circle is shown in both frames in Figure 4; on the left of the figure is the invariant circle in the Jupiter-Europa frame, while on the right is the same circle in the Jupiter-Ganymede frame. As is clearly visible, the shape of the invariant circle in position space is the same in either frame.

We can also convert this 1D CCR4BP stroboscopic map invariant circle to a 2D invariant torus of the continuous time equations of motion in either frame; this is done by integrating the circle states by their respective CCR4BP equations of motion. The resulting 2D torus corresponding to the invariant circle of Figure 4 is plotted in both frames in Figure 5;

it is clear that the 2D flow-invariant torus looks very different depending on the frame used, even though the map-invariant circle was similarly shaped in both frames. Indeed, in the Jupiter-Europa frame, the continuous-time trajectories closely follow the shape of the invariant circle, whereas in the Jupiter-Ganymede frame, a much larger portion of the position space is explored.

To better visualize what is happening, we numerically integrated a single state from the invariant circle to get a continuous-time trajectory lying on the torus of Figure 5; this was then plotted in each frame, as is shown in Figure 6. The invariant circle is plotted in blue, while the integrated trajectory is in green; the trajectory initial condition was chosen to be the intersection of the invariant circle with the negative $x$-axis. The integration time in the Jupiter-Europa frame was $\frac{6\pi}{|\Omega_{3}-1|}$, while in the Jupiter-Ganymede frame it was $\frac{6\pi}{|\Omega_{2}-1|}$; these are the same physical amount of time and are related to each other through the change of time unit between frames. Comparing the trajectory plotted in the two frames, it is visible that the trajectory’s apoapsis passes seem to rotate clockwise in the Jupiter-Ganymede frame over time; by contrast, in the Jupiter-Europa frame, the apoapsis passes remain in the vicinity of the same three locations over time.

In Figure 6, plotting the 3:4 Jupiter-Europa resonant CCR4BP torus trajectory in the Jupiter-Ganymede frame gave a trajectory somewhat reminiscent of a 3:2 Jupiter-Ganymede PCRTBP resonant periodic orbit. Hence, we next computed 3:2 Jupiter-Ganymede unstable resonant CCR4BP tori. During our first attempts at continuing 3:2 Jupiter-Ganymede unstable PCRTBP periodic orbits into invariant tori of the CCR4BP (we used the Jupiter-Ganymede frame for this), we found our quasi-Newton method failing to converge for even extremely small $\bar{\mu}_{2}$ values such as $10^{-12}$.

Upon further investigation, the reason for this was found to be that the addition of Europa to the model introduces singularities in the equations of motion at points belonging to many of the 3:2 Jupiter-Ganymede periodic orbits. In Figure 7, we plot part of the family of unstable 3:2 Jupiter-Ganymede PCRTBP resonant orbits, as well as a red circle representing the orbit of Europa in the CCR4BP. It is clearly visible that the Europa orbit intersects all of the 3:2 Jupiter-Ganymede orbits which pass closest to Ganymede. Hence, no matter how small we take $\mu_{2}$, unless $\bar{\mu}_{2}=0$, we are introducing an infinite perturbation of the equations of motion at states which are included in those periodic orbits; we suspect this destroys any potential invariant tori which would have existed near those periodic orbits for a finite perturbation. We have not investigated whether such tori might exist if one uses regularized equations of motion to remove the singularity at Europa, however.

While all of the most unstable 3:2 Jupiter-Ganymede PCRTBP unstable resonant periodic orbits intersect Europa’s orbit, some of the less unstable ones do not; some of these are visible as the outermost orbits in Figure 7, having higher perijoves and lower apojoves. We were able to continue one such orbit by $\bar{\mu}_{2}$ to the physical value of $\bar{\mu}_{2}=2.5265115494603433\times 10^{-5}$ for Europa in the Jupiter-Ganymede frame planar CCR4BP. The continuation step size used was $\Delta\bar{\mu}_{2}=2.5\times 10^{-7}$. Some of the computed intermediate stroboscopic map invariant circles are shown in Figure 8,

with a zoomed-in plot near Ganymede given on the right.

Just as was seen in the Jupiter-Europa 3:4 resonant CCR4BP torus computation, we see the invariant circles moving closer to Ganymede as $\bar{\mu}_{2}$ is increased. However, this time, plotting the torus’ unstable Floquet multiplier $\lambda_{u}$ versus the parameter $\bar{\mu}_{2}$ in this case shows the opposite trend compared to the Jupiter-Europa 3:4 resonance. In Figure 9, it can be seen that as $\bar{\mu}_{2}$ increases and the circle’s closest approach moves closer and closer to Ganymede, the torus grows more unstable. However, it should be noted that the value of $\lambda_{u}$ is much smaller for this 3:2 Jupiter-Ganymede resonant torus than it was for the 3:4 Jupiter-Europa torus.

Another resonance of interest for Jupiter-Europa-Ganymede tour design is the Jupiter-Ganymede 7:5 resonance[9]. Hence, we computed the corresponding tori in the CCR4BP model. Figure 10 shows the continuation of one such invariant circle by $\bar{\mu}_{2}$ to the physical value for Europa in the Jupiter-Ganymede frame. As $\bar{\mu}_{2}$ went from zero to its physical value, the number of Fourier modes required to sufficiently accurately represent the torus increased significantly, from 512 to 2048.

Using a combination of continuation by $\bar{\mu}_{2}$ and $\omega$, we next computed a family of 7:5 Jupiter-Ganymede resonant tori in the physical Jupiter-Europa-Ganymede CCR4BP; this is plotted in Figure 11. Zooming into the family as shown in Figure 12, we see some interesting phenomena. First of all, it is clear that certain invariant circles see sharp bends at their apojoves; as indicated in the left plot, these apojoves belong to the same invariant circles which make closer passes of Europa’s orbit (note that the 7:5 invariant circles loop twice around Jupiter in position space before closing).

Second of all, as indicated on the right plot of Figure 12, we see a significant gap between the invariant circles we were able to compute. The reason for this is that in this region, the rotation number of topologically similar tori (had they existed) would have been nearly rational; this was verified by computing the continued fraction expansion of the average rotation number of the two tori on either side of the gap. Equivalently, the period of the perturbation from Europa is nearly resonant with the period of the corresponding Jupiter-Ganymede PCRTBP periodic orbits. This causes torus breakdown. Indeed, the “bulge” shape of the gap is immediately reminiscent of the gaps between invariant circles which appear at rational rotation numbers in the commonly-studied Chirikov standard map[15], as the perturbation strength is increased.

To characterize the type of Jupiter-Ganymede PCRTBP dynamical structure we may be approaching as $\bar{\mu}_{2}$ decreases to zero, we integrated trajectories in the Jupiter-Ganymede frame from the $\omega=3.116809$ and $\bar{\mu}_{2}=2.506370\times 10^{-6}$ invariant circle. As mentioned earlier, this gives us the full flow-invariant 2D torus corresponding to the stroboscopic map-invariant circle. We plotted this 2D torus in both $(x,y)$ space as well as in $(x,y,p_{x})$ space; the results are given in Figure 17. Clearly, the resulting structure is still similar to that of the Jupiter-Europa 3:4 resonant CCR4BP tori computed for full $\bar{\mu}_{2}$, such as the one from Figure 5. Our initial hypothesis before carrying out the computations was that by continuing $\bar{\mu}_{2}$ to zero, we might approach a Jupiter-Ganymede PCRTBP stable 3:2 resonant periodic orbit. However, we see no indications of this. Instead, given the results, what we now believe is occurring is that we are approaching a Jupiter-Ganymede PCRTBP stable 2D non-resonant KAM torus as $\bar{\mu}_{2}\rightarrow 0$.

To see this, we converted the torus states shown in Figure 17 to Jupiter-Ganymede synodic Delaunay coordinates. For full details of this coordinate transformation see Celletti[16]; we summarize the basics here. For any state in the Jupiter-Ganymede rotating cartesian coordinate frame, we can compute the instantaneous semi-major axis $a$, eccentricity $e$, longitude of periapse $g$ relative to the Jupiter-Ganymede rotating $x$-axis, and the mean anomaly $\ell$. The Delaunay coordinates are comprised of $L=\sqrt{\mathcal{G}m_{1}a}$, $G=L\sqrt{1-e^{2}}$, $\ell$, and $g$. The most important thing about synodic Delaunay coordinates is that they are action-angle variables for the rotating frame Kepler problem (the PCRTBP with $\mu=0$). Hence, in the Kepler problem, the Hamiltonian becomes a function of only $L$ and $G$; the actions $L$ and $G$ remain constant along trajectories, while $\ell$ and $g$ advance at constant frequencies which are a function of $L$ only ($\dot{g}=-1$, as that is the rotation rate of the coordinate frame). For any $(L,G)$ such that $\dot{\ell}$ and $\dot{g}$ are non-resonant, this results in an invariant 2D torus with $(\ell,g)$ values in $\mathbb{T}^{2}$ is densely filled over time by any trajectory starting on the torus.

When the rotating-frame Kepler problem is perturbed by $\mu>0$ in the CRTBP, the problem is no longer globally integrable; that is, there are no global action-angle coordinates. Nevertheless, many of the 2D stable non-resonant tori described just earlier persist with just small changes. These tori are called KAM tori, and are discussed further by Celletti[16]; indeed, Celletti and Chierchia[17] rigorously proved the existence of such tori inside an energy surface in the Sun-Jupiter CRTBP. Since KAM tori do not change much compared to the tori of the same frequencies in the Kepler problem, the values of the actions $L$ and $G$ do not change significantly on the torus either. However, the full range of possible $(\ell,g)\in\mathbb{T}^{2}$ should still be densely filled in such a KAM torus. There can exist other tori in the perturbed system which do not visit the whole range of all possible $(\ell,g)$ values; these are referred to as secondary tori, and are topologically different from KAM tori.

Converting the torus states shown in Figure 17 to $(L,G,\ell,g)$ coordinates and plotting the resulting torus in these Delaunay coordinates, it is clear that the torus topology is similar to that expected for a KAM torus. The Delaunay coordinate plots in $(\ell,g,L)$ space and $(\ell,g,G)$ space are shown on the left and right respectively of Figure 18. We indeed observe that the whole range of $(\ell,g)$ values in $\mathbb{T}^{2}$ is densely filled, while the actions $(L,G)$ do not significantly vary on the torus. Hence, we propose that the limiting dynamical structure of the Jupiter-Europa 3:4 resonant CCR4BP torus in the Jupiter-Ganymede PCRTBP is a non-resonant KAM torus.

We believe that the reason for this is the slight incommensurability of the periods of Europa and Ganymede around Jupiter; Europa and Ganymede are in an approximate, not exact, 2:1 resonance with each other. Hence, writing $\Omega_{2}$ for the rotational frequency of Europa around Jupiter in the Jupiter-Ganymede frame time units, the frequencies involved in a Jupiter-Europa 3:4 resonant orbit are $3\Omega_{2}/4$ (mean motion), $-1$ (due to $\dot{g}$ caused by the frame rotating with Ganymede), and $\Omega_{2}-1$ (due to Europa). This creates a single resonance as $4(3\Omega_{2}/4)+3(-1)-3(\Omega_{2}-1)=0$; $\Omega_{2}$ is not rational so there are no other resonances. If Europa’s mass is set to zero, the $\Omega_{2}-1$ frequency is decoupled from the other two, so the resonance is broken. The only coupled frequencies left are the $3\Omega_{2}/4$ and $-1$, which are incommensurate, leading to a KAM torus.