Celest Mech Dyn Astr (2009) 105:3–17

DOI 10.1007/s10569-009-9222-5

ORIGINAL ARTICLE

Computation of weak stability boundaries: Sun–Jupiter

system

Francesco Topputo · Edward Belbruno

Received: 15 October 2008 / Revised: 12 May 2009 / Accepted: 7 July 2009 /

Published online: 15 August 2009
© Springer Science+Business Media B.V. 2009

Abstract We compute the weak stability boundary in the planar circular restricted three-
body problem starting from the algorithmic definition, and its generalization by García and
Gómez. In addition, we consider a new set of primaries, Sun–Jupiter, to replace the case
of Earth–Moon considered in previous studies. Numerical enhancements are described and
compared to previous methods. This includes defining the equations of motion in polar coor-
dinates and a modified numerical scheme for the derivation of both stable sets and their
boundaries. These enhancements decrease the computational time. New results are obtained
by considering the Sun–Jupiter case which we compare to the Earth–Moon case.

Keywords Restricted three-body problem · Weak stability boundary · Ballistic capture ·

Invariant manifolds · Lyapunov orbits · Stable sets

1 Introduction

The concept of weak stability boundary, or equivalently WSB, was first introduced in 1986
to design low energy transfers to the Moon (Belbruno and Miller 1993; Belbruno 2004). The
idea was to to find a location about the Moon, in phase space, where a spacecraft could be
captured ballistically, that is, with no propellant. The utilization of this region has generated
new low energy transfers to the Moon, and other locations in space, as is well known, such
as Japan’s Hiten, for example. The WSB region gives rise to a sensitive dynamics about the
Moon, as it defines a location that lies at the transition between lunar capture and escape.

F. Topputo (B)
Aerospace Engineering Department, Politecnico di Milano, Via La Masa, 34, 20156 Milan, Italy
e-mail: [email protected]

E. Belbruno
Department of Astrophysical Sciences, Princeton University, Peyton Hall, Ivy Lane, Princeton,
NJ 08544, USA
e-mail: [email protected]

123
4 F. Topputo, E. Belbruno

In fact, a hyperbolic network was proven to be associated to this region in the restricted
three-body problem (Belbruno 2004).
The original computation of the WSB was based on a straightforward numerical algo-
rithm. This algorithm measures the change in stability of the motion of zero-mass point, P3 ,
in the restricted three-body problem, as it moves about the secondary mass point, P2 , after
making one cycle about it. The motion of P3 is gravitationally perturbed by the primary
mass point P1 , where it is assumed that the mass of P2 is much smaller than the mass of P1 .
Although this definition worked well for applications, it seemed that it could be generalized
further. For example, instead of measuring the stability of motion after one cycle about P2 ,
why not after two or more cycles?
A satisfactory generalization of this region was published by García and Gómez (2007).
They examined the stability of motion of P3 about P2 after it performed n cycles. A general-
ized numerical algorithm was described to do this. Using this methodology, they were able
to define an n-th weak stability boundary. A generalized weak stability boundary was then
defined as the union of these sets for n = 1, 2, 3, . . . They did this for the case of the mass
ratio of P1 , P2 for the Earth–Moon system.
This paper has two purposes. The first is to reexamine the numerical approach in García
and Gómez (2007) by using certain modifications to enhance the numerics. These modi-
fications include the use of polar coordinates, to compute the stable sets, and a bisection
algorithm, to obtain their boundaries. We will demonstrate that these modifications allow us
to compute the stable sets more efficiently and their boundaries with greater accuracy. This
is demonstrated with numerical examples. The second purpose of this paper is to compute
the weak stability boundary for another basic example not yet considered, the Sun–Jupiter
system. In this way, insights into the mass parameter variability of the stable sets can be
gained. The change of mass ratio has indeed the effect of changing the geometry of the stable
sets, where new regions about P2 are found if the Sun–Jupiter case is compared to the Earth–
Moon system (García and Gómez 2007; Romagnoli and Circi 2009). Finally, we discuss the
possible relationship of the weak stability boundary to the invariant manifolds associated to
the Lyapunov orbits.
The paper is organized as follows. The planar circular restricted three-body problem
and its representation in various coordinate systems, including polar coordinates, is defined
below. In Sect. 2 a general algorithmic definition of the WSB is given and its computation is
described in Sect. 3. The computed stable sets and WSB are presented in Sect. 4 for the Sun–
Jupiter system for a wide class of parameter values (i.e., eccentricity and stability number).
Some comments into the dynamics governing the orbits lying on the WSB are described.
Concluding remarks are given in Sect. 5.
This paper can be viewed as a note to García and Gómez (2007) where an additional
basic example is given, using the Sun–Jupiter system, and some numerical enhancements are
described.

1.1 Planar circular restricted three-body problem

The model that we use to describe the motion of a particle of zero mass is the planar circular
restricted three-body problem (or simply RTBP). In this model, relative to an inertial frame,
two primary bodies P1 , P2 of masses m 1 >> m 2 > 0, respectively, move under the mutual
gravity on circular orbits about their common center of mass. The third body P3 , assumed to
be of zero mass, moves under the gravity of the primaries in their plane of motion. In our case,
P3 represents, for example, a spacecraft, and P1 , P2 represent Sun and Jupiter, respectively.

123
Computation of weak stability boundaries 5

Fig. 1 Rotating reference system

The motion of P3 is studied in a rotating coordinate system with coordinates x, y and

whose origin lies at the center of mass of P1 , P2 . This coordinate system rotates with the
same angular velocity as P1 , P2 about their center of mass. In this system, we can assume
P1 , P2 are fixed on the x−axis. Let µ = m 2 /(m 1 + m 2 ). It is well known that we can nor-
malize the units so that the distance between P1 , P2 and their mutual angular velocity are
both 1. Also, we can assume that P1 is located at (−µ, 0) and P2 is located at (1 − µ, 0),
where the mass of P1 is 1 − µ and the mass of P2 is µ. In the case of the Sun–Jupiter system
µ = 9.538754 × 10−4 , which is assumed for the remainder of the paper (Szebehely 1967).
The motion of P3 relative to the normalized co-rotating frame is described by the following
system of second-order differential equations
∂ ∂
ẍ − 2 ẏ = , ÿ + 2 ẋ = , (1)
∂x ∂y
where  is given by
1 2 1−µ µ 1
(x, y) = (x + y 2 ) + + + µ(1 − µ). (2)
2 r1 r2 2
r1 and r2 represent the distances from P3 to P1 and P2 , respectively, and are given by
r12 = (x + µ)2 + y 2 , r22 = (x + µ − 1)2 + y 2 (Fig. 1).
The motion described by (1) has five equilibrium points, known as the Euler–Lagrange
libration points, labeled L k , k = 1, . . . , 5. Three of these, L 1 , L 2 , L 3 , lie along the x-axis,
while, L 4 , L 5 lie at the vertices of two equilateral triangles with common base extending
from P1 to P2 . The RTBP admits also an integral of motion, the Jacobi integral,
J (x, y, ẋ, ẏ) = 2(x, y) − (ẋ 2 + ẏ 2 ). (3)
Thus, each solution x(t) = (x(t), y(t), ẋ(t), ẏ(t)) of (1) always lies on an energy surface
J (C) = {(x, y, ẋ, ẏ) ∈ R4 |J (x, y, ẋ, ẏ) = C}, (4)
for some energy C. This energy surface is a three-dimensional in the four-dimensional phase
space (x, y, ẋ, ẏ) ∈ R4 . The projection of the energy surface J (C) onto the configuration
space (x, y) is called a Hill’s region, and is given by
H(C) = {(x, y) ∈ R2 |2(x, y) − C ≥ 0}. (5)
The motion of P3 is always confined to the Hill’s region of the corresponding Jacobi energy
C. The Hill’s regions vary with the Jacobi energy C. Their topology changes at the values
of C = Ck corresponding to the libration points L k (see Belbruno 2004; Szebehely 1967 for
more details).

123
6 F. Topputo, E. Belbruno

Computing WSB implies integrating many thousands of orbits, and some of them can
result in collision of P3 with either P1 or P2 . As a result, the numerical integration fails as the
equations of motion (1) become singular for r1,2 → 0. It is therefore necessary to regularize
the equations to avoid such singularities. Levi-Citiva regularization is convenient to use. It
is a local regularization method, and requires two different transformations to eliminate the
singularities at both P1 and P2 (see Szebehely 1967; Celletti 2006 for details). We briefly
summarize:
Let z = x + i y be a vector defined in the complex plane having real part x and imaginary
part y, and let w = u + iv be another vector defined in the regularized complex plane (u, v).
The Levi-Civita transformation that regularizes a collision at P2 can be written as

z = f (w) = w2 + 1 − µ. (6)

It is also necessary to transform the independent variable t to a new variable τ . This is given
by
dt
= g(w) = | f  |2 = 4|w|2 . (7)
dτ
It can be shown that the regularized equations of motion can be written as (Szebehely 1967)

w + 2i| f  |2 w = ∇w | f  |2 U, (8)

where U =  − C/2, ∇w = (∂/∂u, ∂/∂v), and (·) = d/dτ . Analogously, it can be shown
that the transformation z = f (w) = w 2 − µ, together with the time transformation (7),
regularizes collisions with P1 .
The regularized equations of motion (8) are numerically integrated when r1 < r LC or
r2 < r LC , i.e., when P3 enters a circle of radius r LC centered at either P1 or P2 . Outside
these circles, the standard equations of motion (1) or those written in polar coordinates are
used.

1.1.1 Polar coordinates

It is convenient to use polar coordinates relative to P2 , (r2 , θ2 ), in our analysis below (see
Fig. 1). These P2 -centered polar coordinates are related to the Cartesian coordinates through

x = 1 − µ + r2 cos θ2 , y = r2 sin θ2 . (9)

The pair (r2 , θ2 ) defines a polar rotating reference system having the first unit vector aligned
with the P2 P3 line (outward) and the second perpendicular to it (in the direction of increasing
θ2 ). In this reference system the equations of motion for P3 turn out to be



1 1−µ µ
r̈2 − r2 θ̇2 − 2r2 θ̇2 = (1 − µ) cos θ2 1 − 3 + r2 1 −
2
− 2,
r1 r13 r2


1
r2 θ̈2 + 2ṙ2 θ̇2 + 2ṙ2 = (1 − µ) sin θ2 −1 , (10)
r13

where r1 , the distance from P3 to P1 , is expressed in P2 -centered polar coordinates as r1 =


r22 + 2r2 cos θ2 + 1. Analogously, using the coordinates (r1 , θ1 ), the equations of motion
in the P1 -centered polar reference frame are

123
Computation of weak stability boundaries 7



1 µ 1−µ
r̈1 − r1 θ̇12 − 2r1 θ̇1 = µ cos θ1 − 1 + r1 1 − 3 − ,
r23 r2 r12


1
r1 θ̈1 + 2ṙ1 θ̇1 + 2ṙ1 = µ sin θ1 1 − 3 , (11)
r2

with r2 = r12 − 2r1 cos θ1 + 1. P1 -centered polar coordinates map into Cartesian coordi-
nates through

x = −µ + r1 cos θ1 , y = r1 sin θ1 . (12)

Both systems (10) and (11) are equivalent to equations of motion (1). The polar coor-
dinates are used in cases where the angle with respect to one of the two primaries, θ1 or
θ2 , is required as a smooth function of time. Nevertheless, when written as systems of four
first-order differential equations, the vector fields generated by (10) and (11) are computa-
tionally more intensive to evaluate than the vector field associated to (1) due to the presence
of trigonometric functions.

2 Definition of weak stability boundary

Let H2 be the Kepler energy of P3 with respect to the the primary P2 . It is given by

1 2 µ
H2 = v − , (13)
2 2 r2
where v2 is the speed of P3 relative to the P2 -centered inertial reference frame. We consider
the case where the motion of P3 starts from the periapsis of an osculating ellipse around P2 .
In this case we have

µ(1 + e)
r2 = a(1 − e), v2 = , (14)
r2

where a and e are the semi-major axis and the eccentricity of the osculating ellipse. Through-
out the paper we consider the case where P3 initially lies at the periapsis of a prograde ellipse,
i.e., P3 rotates counterclockwise in a P2 -centered inertial frame. With relations (14) the initial
Kepler energy becomes

1 µ(e − 1)
H2 = . (15)
2 r2
Let x(t) be a solution of (1). We say that P3 is ballistically captured by P2 at time t1 ,
if H2 (x(t1 )) < 0, and it is temporary ballistically captured (or weakly captured) by P2 , if
H2 (x(t)) < 0 for t1 ≤ t ≤ t2 and H2 (x(t)) > 0 for t < t1 and t > t2 , for finite times
t1 , t2 , t1 < t2 . Regarding the opposite type of behavior, P3 is ballistically ejected (or ballis-
tically escapes) from P2 at a time t1 if H2 (x(t)) < 0 for t < t1 and H2 (x(t)) ≥ 0 for t ≥ t1 .
The region where weak capture occurs can be used to define the weak stability boundary. This
was first defined in Belbruno and Miller (1993) and Belbruno (2004), and further improved
in García and Gómez (2007). We describe the generalized version of the definition of the
WSB stated in García and Gómez (2007), with some refinements.

123
8 F. Topputo, E. Belbruno

Fig. 2 Example of 1-stable and

unstable trajectories relative to P2

We consider trajectories of P3 with the following initial conditions:

(i) The initial position of the trajectory is on a radial segment l(θ ) departing from P2 and
making an angle θ with the P1 P2 line, relative to the rotating system. The trajectory is
assumed to start at the periapsis of an osculating ellipse around P2 , whose semi-major
axis lies on l(θ ) and whose eccentricity e is held fixed along l(θ ).
(ii) The initial velocity of the trajectory is perpendicular to l(θ ), and the Keplerian energy
H2 of P3 relative to P2 is negative, i.e., H2 < 0. The motion, for fixed values of the
parameters θ and e, and for a choice of direction of the initial velocity vector such that
a prograde osculating ellipse is achieved, depends only on the initial distance r .
(iii) The motion is said to be n-stable if the infinitesimal mass P3 leaves l(θ ), makes n
complete turns about P2 , and returns to l(θ ) at a point with negative Kepler energy H2
with respect to P2 , without making a complete turn around P1 along this trajectory.
The motion is otherwise said to be n-unstable (Fig. 2).
It is worth observing that the motion of P3 is unstable if either P3 performs a full circle about
P2 and returns on l(θ ) on a point where H2 ≥ 0 or P3 moves away from P2 and performs
a full circle about P1 . The former is a ballistic escape, the latter is said primary interchange
escape (Belbruno 2004).
We note that the n-stability condition is an open condition. This is due to the fact that
H2 < 0 is an open condition, and that the n-th return map to l(θ ) is smooth, due to the
smooth dependence of solutions of differential equations on initial conditions. Thus the set
of the n-stable points on l(θ ) is an open subset of l(θ ), hence it is a countable union of open
intervals

∗ ∗
Wn (θ, e) = (r2k−1 , r2k ), (16)
k≥1

with r1∗ = 0. The points of the type r ∗ that are at endpoints of intervals above (except for r1∗ )
are n-unstable.
For the sake of a clear exposition, we define the set Wn (e) obtained by varying the param-
eter θ and, at the same time, by holding fixed the eccentricity e

Wn (e) = Wn (θ, e). (17)
θ ∈[0,2π ]

Thus, the n-stable set is


Wn = Wn (e). (18)
e∈[0,1)

123
Computation of weak stability boundaries 9

This is also an open set since the n-stability of points depends smoothly on e and θ .

Definition 1 The weak stability boundary of order n, denoted by ∂ Wn is the locus of all points
r ∗ (θ, e) along the radial segment l(θ ) for which there is a change of stability of the initial
∗
trajectory, that is, r ∗ (θ, e) is one of the endpoints of an interval (r2k−1 ∗ ) characterized by
, r2k
∗ ∗ ∗
the fact that for all r ∈ (r2k−1 , r2k ) the motion is n-stable, and there exist r   ∈ (r2k−1 ∗ ),
, r2k
∗
arbitrarily close to either r2k−1 ∗ for which the motion is n-unstable. Thus
or r2k
∂ Wn = {r ∗ (θ, e) | θ ∈ [0, 2π], e ∈ [0, 1)}.

Analogously, we can define a subset of the n-th weak stability boundary, ∂ Wn (e), obtained
by fixing the eccentricity of the osculating ellipse,
∂ Wn (e) = {r ∗ (θ, e) | θ ∈ [0, 2π]}. (19)
Both Wn (e) and ∂ Wn (e) have been defined for illustration purposes as they have one dimen-
sion less than Wn and ∂ Wn , respectively.
The generalized definition of the weak stability boundary just given slightly differs from
the definition given in García and Gómez (2007) as the openness of Wn sets has been here
considered.

3 Computation of stable sets and weak stability boundaries

From the definition given in the previous section, both the n-stable sets, Wn , and their weak
stability boundaries, ∂ Wn , have to be determined for various values of n. We compute Wn
first, and then its enclosure ∂ Wn . The derivation of these two sets involve the definition of a
computational grid of initial conditions and the integration of hundred thousands of orbits.
The first set, Wn , is defined over a relatively coarse grid, whereas the second set is achieved
on a locally refined grid in order to determine ∂ Wn with a high level of accuracy.

3.1 Computation of Wn

To accomplish this task it is necessary to integrate a large number of orbits and to keep track
of those that are n-stable. An initial condition is uniquely specified once the triple (r, θ, e) is
given. The initial condition to flow under dynamical system (1) is indeed
x = r cos θ + 1 − µ, y = r sin θ,
ẋ = −v sin θ + r sin θ, ẏ = v cos θ − r cos θ, (20)
√
with v = µ(1 + e)/r . The associated orbits have to be integrated for a large time span in
order to assess whether the orbit is stable or not. In García and Gómez (2007), an interval
of 80 adimensional units has been used. Instead of assuming a fixed final time, in this work
we use a different methodology that allows us to get the final time in an automatic way,
depending on the orbit being integrated. In this section we show that this method requires
the orbits to be integrated for a much shorter time.
We recall that the orbits are taken on a radial line l(θ ) emanating from P2 and making
an angle of θ with the x-axis. Using P2 -centered polar coordinates (r2 , θ2 ) introduced in
Sect. 1.1.1, and the associated equations of motion (10), is a more natural choice to derive
Wn sets. Indeed, integrating the system in polar coordinates allows us to handle the angle θ2
as a smooth function of time, θ2 = θ2 (t) [discontinuous inverse trigonometric functions are

123
10 F. Topputo, E. Belbruno

(a) (b)
5
10

4
10

# stable orbits
3
10

2
10

1
10

0
10
0 10 20 30 40 50 60 70 80
t*

Fig. 3 A preliminary view of the W1 (0.0) set derived over the grid (25). In a the black points stand for the
initial conditions that generate orbits requiring a regularization of the dynamics. The integration times, t ∗ , are
reported in b versus the number of stable orbits. a W1 (0.0): regularized orbits; b W1 (0.0): integration times

involved if θ2 is derived inverting Eq. (9)]. By tracking this angle during the integration, we
are able to establish wether P3 makes a complete turn about P2 or not. We say in fact that P3
completes a full turn about P2 at time t ∗ if

|θ2 (t ∗ ) − θ2 (t0 )| = 2π, (21)

for the smallest t ∗ > t0 such that (21) is verified (t0 is the initial time). Using instead the
P1 -centered polar coordinates (r1 , θ1 ), it is possible to define the condition for a primary inter-
change escape to occur in a totally analogous fashion. We say that P3 performs a primary
interchange escape at time t ∗ if

|θ1 (t ∗ ) − θ1 (t0 )| = 2π. (22)

Thus, the numerical integration terminates at t = t ∗ , i.e., when at least one of the two con-
ditions above is verified for the first time. The time span for integration is therefore [t0 , t ∗ ].
Numerical experiments show that t ∗ 80 for all orbits considered (see Fig. 3b where the
t ∗ values are reported for W1 (0.0)). If condition (21) is verified (and the orbit is such that
H2 (t ∗ ) < 0) the motion is stable, otherwise, if condition (22) is satisfied, the orbit is unstable.
(There are certain cases where the two conditions are satisfied simultaneously; in these cases
the orbits are judged unstable.) Thus, in order to assess the stability of the motion, both θ1
and θ2 have to be handled as a smooth function of time. This is possible by simultaneously
integrating the two systems (10) and (11) with initial conditions

r2 (t0 ) = r, θ2 (t0 ) = θ,

µ(1 + e)
ṙ2 (t0 ) = 0, θ̇2 (t0 ) = − 1, (23)
r3

123
Computation of weak stability boundaries 11

Table 1 Computational time required to construct the stable sets integrating the equations of motion in polar
and Cartesian coordinates
CPU time (polar) (h) CPU time (Cartesian) (h) # Stable orbits # Regularized orbits

W1(0.0) 27.3 132.2 28212 3612

W1(0.2) 27.4 132.8 24035 4119
W1(0.4) 27.6 130.3 18816 4230
W1(0.6) 27.8 127.0 14479 4815
W1(0.8) 29.7 125.7 10719 4257
W1(0.95) 29.8 122.3 9106 3180
Both codes are written in Matlab language. The CPU time is relative to a Dual Xeon 2.33 GHz platform running
Linux. Integration tolerance is set to 10−14

for system (10) and

 
r sin θ
r1 (t0 ) = r 2 + 2r cos θ + 1, θ1 (t0 ) = tan−1 ,
1 + r cos θ
ṙ1 (t0 ) = ṙ2 (t0 ) cos(θ − θ1 (t0 )) − r2 (t0 )θ̇2 (t0 ) sin(θ − θ1 (t0 )), (24)
ṙ2 (t0 ) r2 (t0 )θ̇2 (t0 )
θ̇1 (t0 ) = sin(θ − θ1 (t0 )) + cos(θ − θ1 (t0 )),
r1 (t0 ) r1 (t0 )
for system (11). Integrating the coupled system of equations (10, 11) is much more com-
putational intensive than simply flowing Eq. (1). Nevertheless, as clearly shown in Table 1,
integrating the coupled system (10, 11) with a variable final time, determined by either
condition (21) or (22), is a more efficient choice than integrating (1) for a fixed, long time
span.
The coupled initial conditions (23, 24) still depend on the triple (r, θ, e). These three
parameters are defined on the grid
r = {0, 2 × 10−3 , . . . , 1.5},
θ = {0, 2π/1000, . . . , 2π}, (25)
e = {0, 0.05, . . . , 0.95},
and three nested loops are necessary to describe all their combinations. About 7.5 × 105 inte-
grations are needed to generate each Wn (e), whereas the whole Wn sets can be obtained with
1.5 × 107 integrations. A 7th–8th order Runge–Kutta–Fehlberg integration scheme imple-
mented in Matlab is used to flow system (10, 11). When P3 enters either of the two circles
r1,2 < r LC , the state is transformed into Cartesian coordinates and then, through Eq. (6)
(or the equivalent form for regularizing impacts with P1 ) it is mapped into the regularized
variables (u, v). Regularized system (8) is used in this region until P3 reaches the boundary
r1,2 = r LC . In this paper we have taken r LC = 10−3 .
In Fig. 3a we have reported the set W1 (0.0) generated through the computational grid
defined in (25) with the eccentricity held fixed to zero. The 1-stable set is made up by 2.8×104
points, out of the 7.5 × 105 integrated orbits. Among these orbits, a total of 3.6 × 103 has
required a regularization of the dynamics (black points in Fig. 3a). If either Cartesian or polar
equations of motions were used to flow such points, the numerical integration would have
either failed or required a long computational time. The inclusion of a regularization strategy
is therefore crucial for the derivation of the stable sets as a relative high number of initial
conditions result in either collision or close encounter with P2 .

123
12 F. Topputo, E. Belbruno

In Fig. 3b the integration times, t ∗ , in Eq. (21) have been reported for the orbits making
up W1 (0.0). The orbits with final integration times differing by less than 5 adimensional
units have been grouped in subsets for the sake of exposition. It can be observed that the vast
majority of stable orbits perform a loop about P2 , or equivalently verify Eq. (21), in much
less than 80 time units. More specifically, 95% of the stable orbits have a final integration
time less than 10 units, and even no stable orbits with t ∗ greater than 65 units have been
found.
Table 1 reports a comparison of the CPU time required for the derivation of various W1 (e)
sets with increasing eccentricity. The method described so far (i.e., coupled integration of
systems (10) and (11) together with stopping conditions (21) and (22)) is compared with
the integration of Cartesian equations of motion for 80 time units. The advantages of having
introduced a method based on polar coordinates over the method described in García and
Gómez (2007) are evident.

3.2 Computation of ∂ Wn

Once the stable sets Wn have been computed, for fixed values of (θi , e j ) on the discrete
grid (25), there are one or more boundary points on Wn , r s = r (θi , e j ), such that either
r (θi , e j ) + r or r (θi , e j ) − r is an unstable point (r = 2 × 10−3 in agreement with the
grid defined in (25)). In the following we consider the case where r s is an exterior boundary
of Wn (e j ), i.e., r u = r s + r is unstable; extending the method to interior stable points is
trivial.
The stable sets Wn (e j ) have been obtained on the relatively coarse grid (25). The focus
now is the fine determination of the weak stability boundaries ∂ Wn (e j ) on a dynamically
refined grid. We first assume that there is only one point of the type r ∗ (see Eq. (16)) within the
interval [r s , r u ]. This can be claimed as r in (25) is reasonably small, and it is verified by
numerical experiments. Nevertheless, as the stable sets may possess a Cantor-like structure,
it may occur that more than one r ∗ point exists within the interval [r s , r u ]; this case is not
considered in the derivation of the weak stability boundaries.
We implement a bisection method on the interval [r s , r u ]. The midpoint of the interval,
rk = (r s + r u )/2 is in fact taken as a new starting point, and the following initial condition
is flown
r2 (t0 ) = rk , θ2 (t0 ) = θi ,

µ(1 + e j )
ṙ2 (t0 ) = 0, θ̇2 (t0 ) = − 1. (26)
rk3
If (26) generates a stable orbit, then the new stable point is set to r s = rk , otherwise, if (26)
results in an unstable orbit, the new unstable initial condition is updated, i.e., r u = rk . In
any case the bisection method is repeated on the updated, half-size interval [r s , r u ]. In our
implementation, iterations go on until the last interval has amplitude |r u − r s | ≤ 10−8 . The
midpoint of this last interval is assumed to be the point on the weak stability boundary. The
whole process is repeated for all (θi , e j ) defined on the grid (25). In this way we reach an
elevate accuracy in the determination of ∂ Wn since the maximum error in radial distance
from P2 is of the order of 10−9 .
The bisection method described above is computationally intensive and has a slow (linear)
convergence. Nevertheless, this method is effective in the chaotic regions of the phase space
where the weak stability boundary is defined. Furthermore, the high precision reached in the
determination of ∂ Wn will serve, in future works, to demonstrate the relations between these

123
Computation of weak stability boundaries 13

Fig. 4 Weak stability boundary ∂W1 (0.0) (black points) and 1-stable set W1 (0.0) (gray points). The refined
weak stability boundary is obtained by application of the bisection method at the boundaries of the coarse
grid where W1 (0.0) is defined. In b the structure of this grid can be appreciated. a W1 (0.0) and ∂W1 (0.0);
b W1 (0.0) and ∂W1 (0.0), detail

sets and the invariant manifolds of the Lyapunov orbits about both L 1 and L 2 . The weak
stability boundary ∂ W1 (0.0), associated to the set W1 (0.0) shown in Fig. 3a, is reported in
Fig. 4.

4 Visualization of stable sets and weak stability boundaries

The computed n-stable sets, Wn , and their weak stability boundaries, ∂ Wn , are illustrated in
this section for various values of stability number and eccentricity.

4.1 Visualization of Wn

Figure 5 reports the 1-stable sets W1 (e), e = 0.95, in P2 -centered, co-rotating frame intro-
duced in Sect. 1.1. The picture shows 1-stable initial conditions that are highlighted according
to the associated Jacobi energy C, reported on the right sidebar. As observed in García and
Gómez (2007), the set of 1-stable points recalls a Cantor set; i.e., for a fixed value of e,
along the line l(θ ) there are different transitions from stability to instability. This is more
evident for increasing eccentricities. It can be appreciated how the region about P2 shrinks
and, at the same time, becomes more irregular for e → 1. In Fig. 6 we have reported the
n-stable sets Wn (e), n = {1, 2, . . . , 8}, e = 0. It is straightforward that, for fixed eccentric-
ity, Wn (e) ⊆ Wm (e), n ≥ m. It can be observed that the two branches of W1 (0.0) that are
far from P2 suddenly disappears in the other sets Wn≥2 (0.0), whereas the closer branches
becomes thinner for increasing n. Nevertheless, the core of W1 (0.0), i.e., the quasi-circular
region about P2 (see Fig. 6a–h), is conserved in the subsequent Wn≥2 (0) sets. This indicates
that the n-stable sets contain part of the invariant tori surrounding P2 already observed in
Belbruno et al. (2008) for fixed values of C. This property can be summarized by stating

123
14 F. Topputo, E. Belbruno

Fig. 5 The set W1 (0.95). The 1

sidebar on the right reports the
Jacobi energy C of 1-stable initial 3
conditions. The core of
W1 (0.95) has C 3, whereas
the two far branches reach C 1

0.5
2.5

y 0 2

−0.5 1.5

−1
−0.5 0 0.5
x+1−µ

that, in principle, Wn sets are “quasi-invariant” objects for the dynamics: they are regions
in the phase space where orbits stay for a long period of time before leaving for other
regions.

4.2 Visualization of ∂ Wn

The weak stability boundaries, ∂ Wn (e), corresponding to Wn (e) sets discussed previously
are shown. In Fig. 4 the weak stability boundary ∂ W1 (0.0) of the 1-stable set W1 (0.0) is
reported. ∂ W1 (0.0) has been obtained with the bisection method described in Sect. 3.2. In
Fig. 7a, ∂ W1 (e) sets, e = {0, 0.2, . . . , 0.8, 0.95}, associated to the W1 (e) sets shown in Fig. 5
are presented with e on the z-axis. In analogy with the stable sets W1 (e), their boundaries are
more and more irregular for increasing eccentricity. In Fig. 7b the weak stability boundary
∂ Wn (0.0), n = {1, 2, . . . , 8}, associated to the n-stable sets of Fig. 6 are shown with the
stability number n on the z-axis. As in Fig. 6, the central structure is preserved in ∂ Wn (0.0)
as well. This indicates that, as observed for the n-stable sets, the weak stability boundaries
are quasi-invariant objects for the motion of P3 .

4.3 Comparison with the Earth–Moon system

The stable sets and their weak stability boundaries achieved in the Sun–Jupiter system
(µ = 9.538754 × 10−4 ) can be compared to those defined in the Earth–Moon system

123
Computation of weak stability boundaries 15

Fig. 6 Collection of Wn (e), n = {1, 2, . . . , 8}, e = 0. It can be seen that the central structure of W1 (0.0) is
preserved in the subsequent Wn≥2 (0.0) sets. This indicates that the region of Wn (e) about P2 is made up by
invariant tori at different values of C. a W1 (0.0); b W2 (0.0); c W3 (0.0); d W4 (0.0); e W5 (0.0); f W6 (0.0);
g W7 (0.0); h W8 (0.0)

Fig. 7 a ∂W1 (e), e = {0, 0.2, . . . , 0.8, 0.95}, as boundaries of the 1-stable sets W1 (e), shown in Fig. 5.
b ∂Wn (0.0), n = {1, 2, . . . , 8}, as boundaries of the n-stable sets Wn (0.0), shown in Fig. 6

(µ = 1.21506683 × 10−2 ). It is indeed clear that these sets depend on the mass ratio of
the three-body system, as this is the only parameter characterizing differential equation (1).
Plots reported in Figs. 4, 5, 6 can be compared to the visualizations in García and Gómez
(2007), where the Earth–Moon system has been investigated. As a sample comparison, we
have reported in Fig. 8 the set W1 (0.2) defined in both systems.
Not surprisingly, the scale changes, and the relative size of the two sets is different
as the Earth–Moon mass ratio is about 12 times greater than the Sun–Jupiter one.
Nevertheless, the dependence of the stable sets geometry on the mass ratio is not straight-
forward. By comparing Fig. 8a,b, it can be seen that the island on the upper right in the
Earth–Moon system disappears, and branches occur in the Sun–Jupiter system (right part)
that are not present in the Earth–Moon. Moreover, there are substantially new arcs in the upper

123
16 F. Topputo, E. Belbruno

Fig. 8 The set W1 (0.2) in

Sun–Jupiter (left) and
Earth–Moon (right) systems.
a W1 (0.2), Sun–Jupiter;
b W1 (0.2), Earth–Moon

Fig. 9 A portion of 0.2

∂W1 (e), e = 0, and a set of
sample orbits produced by taking 0.15
random points with
θ ∈ [−π/4, π/4] as initial
conditions (Earth–Moon system). 0.1
There is evidence that these
orbits shadow the stable 0.05
manifolds of the Lyapunov orbits L2
about both L 1 and L 2 (shown
y

0
with ‘+’ markers) L1

−0.05

−0.1

−0.15

−0.2
−0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2
x+1−µ

left part of the Sun–Jupiter W1 (0.2) set. Overall, besides these changes, the two regions are
qualitatively similar, although there is evidence that the smaller mass ratio is, the richer the
stable sets are—which is not trivial.

5 Conclusions

In this paper we have presented a methodology to compute and visualize weak stability
boundaries that is a refinement of García and Gómez (2007), and also for a new fundamental
case, the Sun–Jupiter problem. To make the existing methodology more efficient, a set of
differential equations is written in polar coordinates, together with the formulation of a stop
condition that allows us to significantly reduce the effort required to derive the stable sets.
In addition, a bisection algorithm has been formulated for the fine definition of the weak
stability boundaries. We have also found that the mass ratio variability produces important
changes in the geometry of the stable sets, as the Sun–Jupiter case is different from the
Earth–Moon one.

123
Computation of weak stability boundaries 17

It is well known that the weak stability boundary plays a key role in the design of low
energy transfers as first demonstrated in Belbruno and Miller (1993), Belbruno (2004) and
more recently shown in Romagnoli and Circi (2009). More significantly, however, this object
seems to play a key role in the dynamics of motion of P3 about P2 . It is seen that ∂ Wn is not an
invariant object, which makes its study more difficult. It also allows the dynamics of motion
of P3 about P2 to be studied without directly analyzing the dynamics associated with the
invariant manifolds of the Lyapunov orbits. The weak stability boundary is described by a rel-
atively simple algorithmic definition. Nevertheless, there is evidence that it is a dynamically
rich set that deserves more attention.
As a final example, we show in Fig. 9 portion of ∂ W1 (e), e = 0, in the Earth–Moon
system. We have taken random initial conditions for θ ∈ [−π/4, π/4], and we have numer-
ically integrated these initial points under the RTBP dynamics (1). There is evidence that the
orbits produced by these points shadow that stable manifolds of the Lyapunov orbits about
both L 1 and L 2 . Ongoing work by the same authors is indeed aimed at demonstrating that,
for a significant range of energies such that the Lyapunov orbits exist, a portion of the weak
stability boundaries locally overlap with the invariant manifolds of the Lyapunov orbits about
both L 1 and L 2 .

Acknowledgments The authors would like to acknowledge Marian Gidea for useful discussions and
suggestions. E. B. would like to acknowledge support from NASA/AISR.

References

Belbruno, E.: Capture Dynamics and Chaotic Motions in Celestial Mechanics. Princeton University
Press, Princeton (2004)
Belbruno, E., Miller, J.: Sun-perturbated Earth-to-Moon transfers with Ballistic capture. J. Guid. Control
Dyn. 16, 770–775 (1993)
Belbruno, E., Topputo, F., Gidea, M.: Resonance transitions associated to weak capture in the restricted
three-body problem. Adv. Space Res. 42, 1330–1351 (2008)
Celletti, A.: Basics of Regularization Theory, in Chaotic Worlds: From Order to Disorder in Gravitational
N-Body Dynamical Systems. pp. 203–230. Springer, Netherlands (2006)
García, F., Gómez, G.: A note on weak stability boundaries. Celest. Mech. Dyn. Astron. 97, 87–100 (2007)
Romagnoli, D., Circi, C.: Earth–Moon weak stability boundaries in the restricted three and four body prob-
lem. Celest. Mech. Dyn. Astron. 103, 79–103 (2009)
Szebehely, V.: Theory of Orbits: The Restricted Problem of Three Bodies. Academic Press Inc,
New York (1967)

123