Numerical Exploration of the Limit Ring Problem

The aim of this work is to provide an insight of an idealized model of a planetary ring. The model is a limit case of the planar circular restricted 1 + n body problem, where an infinitesimal particle moves under the gravitational influence of a large central body and n smaller bodies located on the vertices of a regular n-gon. When considering n tending to infinity, a model depending on one parameter is obtained. We study the main important structures of the problem depending on this parameter (equilibria, Hill’s regions, linear stability, …). We use Poincaré maps, for different values of the parameter, in order to predict the width of the ring and the richness of the dynamics that occur is discussed. This work is a continuation of the work presented in Barrabés by (SIAM J Appl Dyn Syst 9:634–658, 2010).


Introduction
In the lovely book, "Periodic Solutions of the N-Body Problem" [6], Ken Meyer shows the great power of exploiting the differences in relative masses and relative distances between bodies in order to form interesting and useful limiting problems in Celestial Mechanics. This paper is one more application of this idea.
The study of the dynamics of planetary rings dates back at least to Maxwell [5] who showed, among other things, that the relative equilibrium formed by an n-gon of small equal masses in orbit around a single large mass is spectrally stable (see also Moeckel [8]). Meyer and Schmidt in [7] use also Maxwell's model to study braided rings. While this is an appealing caricature of a planetary ring, observed rings are not so simple, containing numerous particles of various sizes. Adding a swarm of infinitesimal particles orbiting in the vicinity of the n-gon makes a slightly more realistic model. Since each of these particles is assumed to be tiny, one may assume that they do not interact and hence study the restricted problem in which the motion of a test particle is determined by its attraction to the large central mass and the small masses of an orbiting n-gon. This is a restricted (1 + n) + 1-body problem because of the relative sizes of the bodies. One goal of our work is to study the dynamics of particles in such a ring in isolation, without the constraining forces provided by orbiting moons.
This problem has been studied numerically for small n by Kalvouridis (see the review article [4] and the references therein). Since there is no natural choice for the number of bodies in the n-gon and observed rings are not observed to have imbedded n-gons with (relatively) large ring bodies, it is reasonable to look at the limit problem when n is taken very large and the mass of each of the n-gon bodies is very small. We studied this limit in [1], and this paper is a continuation of that study.
There are several possible limits, depending on the relative sizes of the central mass and n-gon bodies, the relative distance between the n-gon bodies with each other and the central mass and the angular velocity of the n-gon. The most interesting limit keeps terms for attraction to the central mass and to the n-gon and contains a parameter ( 0 ) which relates these sizes and distances. In [1], we focused on relatively large values of 0 since those give stability of the n-gon relative equilibrium. In this paper, we study how the dynamics of the limit restricted (1 + n) + 1-body problem evolves with the parameter 0 , particularly the bifurcations and dynamics of the equilibria, their stable and unstable manifolds and the effect on the global dynamics. Section 2 outlines the derivation of the equations of motion of the limit ring problem. More details can be found in [1]. Section 3 discusses the existence and stability of the equilibria of the system as the parameter 0 is varied. Analogies to the restricted three-body problem, particularly the bifurcations in the number and stability of the equilibria and in the "Hill's regions" are discussed.
Symmetries in the limit problem make it possible to set up natural surfaces of section for the computation of Poincaré return maps. Unfortunately, the nature of the equations as infinite sums also immediately introduces a source of computation error that can quickly accumulate when computing an orbit over a long time interval. In Sect. 4 we consider techniques to control this error. The Poincaré maps are studied for various 0 values in Sect. 5 for 0 values not considered in [1] and the stable Fig. 1 The n-gon relative equilibrium configuration of radius 1/ρ around a large central body of mass m 0 , with the distance between two consecutive peripherals equal to 1 P 0 P -1 P 1 1/ρ 1 P k 1 P -k 0 and unstable manifolds are computed numerically. As one might expect, the stable and unstable manifolds give considerable insight into the possible motions of the infinitesimal particle and the richness of the dynamics.

Equations of Motion
In this section we summarize the results stated in [1], Sect. 2, where the limit process and basic properties of the limit system were discussed.
We start considering an infinitesimal particle in rotating coordinates (x, y) with angular velocity ω, moving under the influence of a large central body, called the primary, and n smaller bodies, called peripherals or ring primaries, arranged in a regular n-gon about the central body (see, for example, [3]). The primary is located at the origin with mass m 0 , and the n peripherals have each one the same mass m. The equations of motion of the infinitesimal particle arë where G is the gravitational constant and (x i , y i ) denote the position of the P i ring primary for i = −k 0 . . . k 1 , being k 1 = [n/2] for all n ∈ N, and k 0 = k 1 if n is odd or k 0 = k 1 − 1 if n is even (see Fig. 1).
Since our goal is to create a model for the motion of infinitesimal particles near a regular n-gon relative equilibrium for large n, we choose a limit process that keeps the individual ring primaries distinct and keeps both the force from the central planet and the ring primaries from disappearing.
Choosing the length scale so that the distance between two consecutive peripherals equals one, the radius of the n-gon becomes 1/ρ, ρ = 2 sin(π/n), and moving the origin to the position of P 0 the equations of motion (2.1) becomë where the components of the field (F n,1 , F n,2 ) are The constant β = m 0 /m is the ratio between the mass of the planet and the mass of one ring particle and = ω 2 /(Gm). The relative size of the quantities β, and ρ is important when n grows and are related by the equation We can now begin the process of finding the appropriate limit system. Since the distance from the ring to the planet grows as n → ∞, we must let β = m 0 /m grow so that the force from the planet does not drop out of the limit problem. The next lemma (see [1]) relates the number of ring primaries n to the size of the planet through the parameter β and the radius 1/ρ of the ring.

Lemma 2.1
Assume β is such that βρ 3 has a finite non-zero limit as n tends to infinity and let be the quantity given by Eq. (2.2). Then, has the same finite non-zero limit as βρ 3 .
Throughout the rest of this paper, we assume the hypothesis of Lemma 2.1 and we denote by This means that for a large value of the number of ring primaries n, where m R = n m is the total mass of the ring. We can interpret the physical meaning of 0 in the following way: for a fixed big value of n, the bigger the ratio between the mass of the planet and the mass of the ring, the bigger the value of 0 . We will consider values of 0 ≥ 1, so we avoid the case of massive rings in comparison with the central planet.
Next we want to let n tend to infinity in order to obtain the limit system. For n → ∞, the ring primaries form a very large circle as the radius of the n-gon grows to infinity. Those ring primaries near the origin lie almost along the y axis and they tend toward the points (0, k), k ∈ Z, as n increases (see Fig. 2). However, there are always ring primaries "on the other side" of the planet. These bodies are getting so far away that their effect on the motion of the test particle near the origin tends to zero. Our limit  Fig. 2 Location of the ring primaries as n grows. On the right, a close look around the origin process exploits this subdivision of the ring primaries, and it is explained in detail in [1].
The resulting system, which we call the "limit ring system", LRP, is defined by the equationsẍ Notice that for any solution (x(t), y(t)) of the Eq. (2.3), replacing y(t) by y(t) ± k for any integer k, yields another solution. Hence, any portion of any solution can be translated by an appropriate integer so that −1/2 < y(t) ≤ 1/2. Any solution which decreases across y = −1/2 can be studied by following the translation of this solution "up" by one unit in y, and similarly for solutions which increase through y = 1/2. Hence, it suffices to study solutions in the strip {(x, y)| x ∈ R, −1/2 < y ≤ 1/2}. Topologically, we are creating a new phase space by identifying points which differ by an integer translation in the y coordinate. We can also think of this as identifying the lines y = −1/2 and y = 1/2, giving a phase space of a cylinder for position coordinates. Since the reflection on the y axis also leaves the vector field unchanged, it suffices to study the strip −1/2 < y ≤ 1/2 with x ≥ 0.
Also notice that the Eq. (2.3) has a first integral, called the Jacobi constant given by (2.5) Finally, we will write the LRP as a system of first differential equationṡ where z = (x, y, u, v), being u =ẋ and v =ẏ, and (2.6)

Equilibrium Points and Hill's Regions
The equilibrium points of LRP (2.3) are the solutions of the following algebraic system We will see that bifurcations in the number and type of equilibria occur as the parameter 0 is increased. Recall that, due to the symmetry it suffices to study the zeros of the above system in the strip −1/2 < y ≤ 1/2 and x ≥ 0. For each equilibrium point at (x, y) with x > 0, there is another one located at (±x, y ± k), for all integer k.

Computation of the Equilibrium Points
The existence and the analysis of the equilibria is a straightforward computation. Uniqueness is also straightforward computation, with the exception of the following statement, that ensures that all the equilibrium points are either on the y = 0 axis, or on the boundary of the strip |y| = 1 2 .
Lemma 3.1 Let x be any real number and y ∈ [−1/2, 1/2]. Then the equation is satisfied if and only if y = 0 or y = ± 1 2 .
Clearly, if y = 0 or y = ± 1 2 , the above equation is fulfilled. The uniqueness of these solutions has been only verified numerically.

Equilibrium Points on y = 0
The first equation in (3.1) can be written for y = 0 as since x = 0 is not an admissible solution. The next result states that the above equation has only one solution for x > 0, and also gives some properties of it. 1.
Proof Let us denote by It is easy to see that the function f 1 (x) is positive and decreases strictly on (0, ∞). Moreover, lim x→0 + f 1 (x) = +∞ and lim x→+∞ f 1 (x) = 0. Then, using the classic Theorems of Bolzano and Rolle, given any value 0 , there exist one and only one solution of f 1 (x) = 3 0 for x > 0, denoted by x 1 .

Equilibrium Points on |y| = 1/2
It suffices to consider y = 1/2. Then, the first equation of (3.1) becomes In this case x = 0 is an admissible solution, so we have an equilibrium point at (0, 1/2) for all values of 0 , that we call L 2 .
The next result states that, depending on the value of 0 , there exist another equilibrium point on y = 1/2 and x > 0. , for x ≥ 0, so the solutions of f 3 (x) = 3 0 /16 will be equilibrium points. Clearly, the function f 3 (x) decreases on (0, ∞) and has its maximum value at x = 0, which is given by Thus, using again the classic Theorems of Bolzano and Rolle, given any value 0 , the equation f 3 (x) = 3 0 /16 has a unique positive solution if and only if 0 < 16S 3 /3.

Linear Stability of the Equilibrium Points
In this section we present the results about the linear stability of the equilibrium points of the LRP. The differential of the vector field (2.6) at the equilibrium points As the LRP has a Hamiltonian formulation, the eigenvalues of M are ±λ 1 , ±λ 2 , where each λ k satisfy the equation We have obtained the following results: 1. Equilibrium point L 1 : It can be shown that the eigenvalues of the matrix M evaluated at L 1 are the solutions of In conclusion, L 1 is of type center × saddle. 2. Equilibrium point L 2 : In this case is easy to see that where S 3 is defined in Lemma 3.3. Recall that at value * 0 = 16S 3 /3 the point L 2 bifurcates and a new equilibrium point appear for smaller values of 0 . We obtained three scenarios: (a) For 0 < * 0 , two real and two complex eigenvalues are obtained. Thus, L 2 is of type center×saddle. (b) For 0 = * 0 , the eigenvalues are λ = ± √ 2 and λ = 0 (double). (c) For 0 > * 0 , the sign of the discriminant that appears when calculating the eigenvalues is not constant, so again, we can find different scenarios depending on 0 . By careful calculation, one can prove the following: , there are four complex eigenvalues, so L 2 is of type complex saddle. (iii) For 0 ≥ 32S 3 /σ 2 , there are four pure complex eigenvalues, so L 2 is of type center×center. In Fig. 4, left, some of these cases are shown for a certain range of values of 0 . 3. Equilibrium point L 3 : in this case we have explored the eigenvalues numerically, the results are plotted in Fig. 4, right. We can see that there are two intervals of values of 0 , say (0, r 1 ] and [r 2 , 0 * ), for which the L 3 point is of type saddle×saddle, while for 0 ∈ (r 1 , r 2 ) it is a complex saddle.

Jacobi Constant Value at the Equilibrium Points
Let us denote by C i = C ∞ (x i , y i , 0, 0), i = 1, 2 the value of the Jacobi constant (2.4) at the equilibrium points L 1 and L 2 . Next two Lemmas state the behavior of C i , i = 1, 2 as a function of 0 .

Lemma 3.4
The value of C 1 as a function of 0 satisfies Proof Using (2.4), the value of the Jacobi constant at L 1 can be written as From the proof of Lemma 3.2, we have that lim 0 →0 0 x 2 1 = 2/3, and from (2.5) which proves the first statement. Also, from the proof of Lemma 3.2, we have that lim 0 →+∞ 0 x 3 1 = 1/3, and the second statement follows from the inequality Finally, it can be shown using the definition of x 1 that which implies that C 1 , as a function of 0 , has only one critical point that must be a maximum.
The plot of C 1 for certain values 0 > 1 is shown in Fig. 5. The maximum of the function C 1 is reached for a value 0 < 1.
There exists only one value of 0 such that C 12 = 0.
Proof Using (2.4), the value of the Jacobi constant at L 2 can be written as .
The first two statements come from this expression and Lemma 3.4. Let s be , so that C 2 = 4s/ 0 . Note that C 12 has exactly one critical point, which is a maximum, which implies that the function crosses the horizontal axis only one time.
On one hand, using Lemma 3.2, we have that 1/x 1 is an increasing function of 0 , going from 0 to +∞. On the other hand, 0) is an increasing function. Thus, C 12 only has one critical point, which has to be a maximum. As the energy decreases, Lemma 3.5 implies that, for values of 0 less than a certain value, the equilibrium point L 2 appears before than L 1 , while for a bigger values of 0 , L 1 appears before than L 2 , see Fig. 5, where C 1 and C 2 as functions of 0 ∈ [1, 10] are shown. This fact has to be taken into account when visualizing the zero velocity curves and the Hill's regions.

Zero Velocity Curves
As in the restricted three-body problem the Jacobi constant given by (2.4) can be used to restrict the possible locations of the infinitesimal body.
From (2.4), the regions of admissible motion in the configuration space must have U ∞ ≥ C ∞ /2. The zero velocity curves separate the (x, y) plane into regions that are accessible and inaccessible for solutions with Jacobi constant equal to C ∞ , which are given by The zero velocity curves only exists for values of the Jacobi constant greater than a certain level that depends on 0 . Using that (the last inequality has been used in the proof of Lemma 3.4) we have that the right hand side of Eq. (3.3) has a lower boundary (which depends on 0 ). Thus, as the Jacobi constant decreases, the regions of forbidden motion shrink and finally disappear.
As we have seen before, there is a value for 0 such that, for greater values, as the Jacobi constant decreases, the equilibrium point L 1 appears before L 2 . Furthermore, for 0 > * 0 , the L 3 point does not exist. Thus, the zero velocity curves will look different depending on the values of 0 .
For all the values of 0 and big values of C ∞ , the region of admissible motion is essentially as where B k is a neighborhood of the k-th ring primary and R = {(x, y); |x| > M}, for a certain M. That is, the motion can only take place in any neighborhood of a single body, or far away from them. As the Jacobi constant decreases, the Hill's region grows, the equilibrium points appear and finally, the regions of forbidden motion disappear.
In the following plots, some zero velocity curves are shown for different values of 0 . For each energy level, the forbidden region is, with respect the x direction, the region bounded by the zero velocity curves of that energy. Notice that these regions are unbounded in the y direction.
The plots in Fig. 6 are representative of the following situations: • 0 = 2: for big values of C ∞ , the motion is only allowed around a single ring primary or far from them. As C ∞ decreases, the equilibrium point L 2 appears, and the allowed regions of the motion on either side of the peripherals connect. Next, L 1 appears, and as C ∞ continues decreasing, L 3 appears, after which there is no forbidden region for the motion anymore. • 0 = 50: in this case the problem only has two equilibrium points, and as the energy decreases, L 1 appears before than L 2 .

Approximations and Errors
In order to obtain accurate numerics of the LRP the main problem that must be tackled is the fact that the infinite sum must be replaced by a finite one. The easiest way to do so is using a truncated system, that is, replacing into the equations, so only n = 2N + 1 bodies are considered. In doing so, the error created is of order O(( 0 n) −1 ), but even with large values of N , the error grows as time increases. Furthermore, the truncated problem is a poor approximation when a solution approaches y = ±N .
A new difficulty arises when considering the cylinder {(x, y) | − 1/2 ≤ y < 1/2} with the lines y = 1/2 and y = −1/2 identified as a configuration space of the LRP. Using the symmetries of the problem, any solution which cross one of the two lines is connected via a jump to the other line, so the entire solution, in configuration space, "lives" inside the strip between the two lines. Although using the strip for the numerical computations allow us to avoid approaching the lines y = ±N , when deleting the influence of the bodies located at |y| > N , the symmetry is lost and an error is introduced at each jump. The solutions, in general, perform multiple jumps, so the error is multiplied by the total number of jumps performed.
One way to reduce the error when using an approximated problem for the numerical integrations consists in replacing the terms corresponding to |k| > N in the infinite sums of the vector field of Eq. (2.3) by the approximated integrals. This choice reduce the error to order O( −1 0 n −2 ). In the following subsections, we describe the approximations used and we give an estimate of the errors involved in each one.

The Linear Approximation
The simplest approximation that can be done consist in removing all the influences of the peripherals and keep only the first term of the Eq. (2.3), that is, Writing the above equations as a system of first order with variables (x, y, u =ẋ, v = y), the solution with initial conditions (x 0 , y 0 , u 0 , v 0 ) is given by where A = 3x 0 + 2v 0 , B = u 0 , c 1 = −6x 0 − 3v 0 , and c 2 = y 0 − 2u 0 . The linear system has a "Jacobi constant" given by Notice that the solutions are periodic in the x component, while the y component essentially decreases or increases: c 1 < 0 (or c 1 > 0) for any initial condition with x 0 > 0 (respectively x 0 < 0). That is, for x 0 > 0, the orbits move "downwards" (except for local loops) in configuration space, crossing the lines y = ±1/2−k, k ∈ N. These orbits are in general quasi-periodic orbits in the cylinder {(x, y)| x ∈ R, −1/2 < y ≤ 1/2} (depending only on x 0 they can be also periodic), so all the solutions live inside invariant tori. Going back to the original reference system for a finite value of the number of peripherals, the orbits turn around the ring in clockwise sense without approaching it.
We can expect the same behavior for solutions of the LRP with initial big values of the initial condition x 0 . As we show in [1] (for a specific value of 0 = 110), far from the peripherals, most of the invariant tori survive (see the plots of Poincaré sections in the Sect. 5), while close the ring primaries they are destroyed and the dynamics becomes more chaotic.
It can be proven that the terms removed are of order O(( 0 x) −1 ), so for big values of x, the linear approximation is good enough, but as we get closer to the ring primaries we need to improve the approximation in order to obtain reliable numerics.

First Approximation: Truncated Problem
Consider the LRP given by Eq. (2.3). The easiest thing to do (after the linear approach) in order to give an approximation is to consider only a finite number of terms of the potential U ∞ given by (2.5). That is, fixed an odd number of bodies n = 2N + 1, we write where U n contains n terms corresponding to the bodies k = −N , . . . N , 2 is the distance between the infinitesimal particle and the k-th ring primarie, and W n is the remainder. A first approximation of the LRP can be obtained by removing W n from the potential (and thus, the influence of the farthest ring primaries). Considering again z = (x, y, u, v), we define the truncated LRP as the one given by the equationṡ where F n (z) is defined as F(z) in (2.6), replacing U ∞ by U n . Clearly, the problem has a first integral, that we also call "Jacobi constant" given by We consider the strip without a neighborhood of the origin for a fixed small δ > 0. Applying the Gronwall's Lemma, the next Proposition can be proven.
for all t such that the projections of z(t) and z n (t) in the configuration space are in D, and where K satisfies that While the orbits z(t) and z n (t) remain inside the strip, we have that both solutions differ a quantity of order O = ((n 0 ) −1 ). Each time the solution z(t) performs a jump (from y = −1/2 to y = 1/2 or vice versa), when continuing the integration with the truncated LRP we introduce an error of the same order. Thus, we are approximating an orbit z(t) of the LRP by a sequence of pieces of solutions of the truncated system. We have to take into account the total number of jumps m that the orbit performs to stay inside the strip when considering the total error at the final point, which is O m n 0 .

Second Approximation: Extended Truncated Problem
In order to improve the error, we substitute the terms of the reminder, W n , that we eliminate in the first approximation by its continuous version. Thus, the terms are replaced into the equations of the LRP by N (x, y) .
We define the extended truncated LRP as the problem given by the equationṡ where F n (z) is the same field of the truncated LRP (4.3) and The next Lemma ensures that the vector field F(z) is well defined in the domain D and gives the potential associated to the extended truncated problem.
where U n (x, y) is given in (4.2) and W n (x, y) is defined as Then, 1. the function W n (x, y) is at least C 1 in D and ∇W n = G n . Furthermore, the second order derivatives of W n are bounded functions in D.

U (x, y) is the potential associated to the vector field F(z) and
uniformly in any compact subset of D.
Proof Using Taylor expansions we have that for y 0 = 0 which ensures that W n is continuous.
With respect the derivatives, it is straightforward to check that, and lim (x,y)→(0,y 0 ) so the first part of the statement is proven. Furthermore, it is also easy to see that uniformly on each compact subset of D, which ends the proof.
Again applying the Gronwall's Lemma, we obtain the next Proposition. Therefore, the error involved using this second approximation of a solution of the limit ring problem will be of order being m the number of jumps from y = −1/2 to y = 1/2. Clearly the error involved using the extended truncated problem is much better than the one using the truncated problem. Nevertheless, we have to be careful when evaluating the vector field F close to the x = 0 axis, due to the definition of the function g 1 (x, y). But, close to x = 0, the dynamics is governed by the influence of the ring primary located at the origin, so an infinitesimal particle approaching the axis x = 0 is caught around the peripheral, or revolves around it before leaving a neighborhood of the origin (the presence of the invariant manifolds associated to the equilibrium points are the key to describing the dynamics close and within the ring). In any case, the orbit stays within the strip {(x, y) | |y| < 1/2} without making any "jump". In this case (when the orbit goes close to the x = 0 axis), as there are no errors due to jumps, the approximation given by the truncated problem is good enough.
Therefore, depending on the region where we want to follow the solutions of the LRP, we will use the most suitable of the two approximations explained. Whereas the integration is done far from the x = 0 axis, the approximation used is the one given by the extended truncated problem. This is the case of the Poincaré section plots showed in the next Section. If an integration crosses the x = 0 axis, then, in a neighborhood of the axis the equations used are those of the truncated problem, going back to the extended truncated problem as the orbit moves away from the axis.

Numerical Explorations
In this Section, we show some numerical explorations varying the value of the parameter 0 , in order to compare the similarities and the differences. The main explorations done correspond to Poincaré section maps. The orbits can display two different types of behavior, depending on the distance to the ring primaries. Near the peripherals, the dynamics are governed essentially by the stable/unstable manifolds associated to the equilibria and Lyapunov periodic orbits, so they are rather complicated. Far from the ring primaries, the system is well approximated by the linear approximation, so most of the invariant tori persist. The boundary between these regions is formed by the KAM circles closest to the ring primaries and we are particularly interested in approximating the location of the "innermost" invariant tori that separates the "regular" from the chaotic region. Its distance to the origin can be considered as the "width" of the ring around the peripherals. Of course, it depends not only on the value of the parameter 0 but also on the Jacobi constant, C ∞ . We will show that the width is bigger as the value of 0 and C ∞ are lower.
Furthermore, we will show that for any value of 0 the presence of the invariant manifolds associated to the Lyapunov orbits around the equilibrium point L 1 allows Table 1 For each value of 0 in the first column, the table includes the x coordinate of L 1 , the value of the Jacobi constant C ∞ at the equilibrium points L 1 and L 2 , and the values of C ∞ chosen to show the Poincaré sections 0 The data is computed using the truncated problem with 600001 bodies the existence of heteroclinic connections, as well as transit and non transit orbits connecting neighborhoods of different ring primaries. The values of 0 chosen have been 50, 110 and 550. We have not considered smaller values of it in order to keep the errors reasonably bounded, so the smallest one considered is 50. 0 = 110 is close to the first value of 0 for which the finite 1 + n ring problem is stable. The bigger value considered 550 is chosen so that the equilibrium point L 2 has a different linear stability than in the case of 0 = 110. For each one of these values, we have considered two fixed values of the Jacobi constant C ∞ = C a , C b . The first one, C a is always close to and small than the value of the Jacobi constant at L 1 , C 1 . That means that the region of motion is connected by a bottleneck around the location of L 1 , so motion from away the ring primaries toward them is possible (see Fig. 6, plots in the middle). The second value considered, C b correspond to a value of the Jacobi constant such that the zero velocity curve does not exists, and the entire plane (x, y) is accessible for the motion.
In Table 1, we show, for each chosen values of 0 , the x coordinate of L 1 , the value of the Jacobi constant at L 1 and L 2 and the two values of the Jacobi constant C a and C b chosen to perform the Poincaré sections.

Poincaré Section Maps
We have computed some Poincaré section maps in order to plot them and to show how the invariant tori that persists far from the ring primaries, are destroyed as we move close to the origin. These plots show a region in the configuration space where there is a transition between regular motion (trapped inside invariant tori, which occurs far from the ring) and chaotic motion (the orbits escape towards the ring).
In order to perform the integrations, we have used the extended truncated problem. We follow the Poincaré iterates while the orbit does not approach the x = 0 axis. If an orbit approaches the ring primaries before the total number of Poincaré sections are performed, the integration is stopped.
We have seen that the error accumulated using one of the approximated problems depends not only on the number of bodies n considered, but on the number of jumps from y = −1/2 to y = 1/2 performed during the integration. We would like to have an idea of how big this error can be, for example, when considering a certain number p of iterates of a Poincaré map. For solutions that do not approach the ring, we know that their behavior is quite regular and similar to the solutions of the linear problem, so the number of crossings performed between two consecutive iterates of the Poincaré map can be considered the same as the linear approximate solution. In that case, the number of crossings between two consecutive iterates can be calculated and is almost constant, m p . Then, the total number of crossings with |y| = 1/2 that an orbit can perform is approximately m p p, and the error involved in the approximation will be of order Fixed a value of the Jacobi constant C 0 , we consider the Poincaré section given by where C(x, y, u, v) is the Jacobi constant first integral of the problem, and the Poincaré return map so, m p is the number of crossings with |y| = 1/2 between (x 0 , y 0 ) and (x 1 , y 1 ).
In order to have an idea of how big m p can be, we use the linear approximation. We consider initial conditions z 0 = (x 0 , y 0 , 0, v 0 ), such that 3x 0 + 2v 0 < 0 (so z 0 ∈ 0 ). As we only consider x 0 > 0 (due to the symmetry of the problem with respect the y axis), this implies that the initial velocity v 0 must be negative, so from the Jacobi first integral For the linear problem, the Poincaré return map correspond to follow the flow 2π units of time. Thus, from (4.1), we have that the Poincaré map for the linear problem writes Using the explicit expression of the Poincaré map we have that where It is worth noticing that the bigger the value of x 0 , the bigger the number of crossings m p . In fact, for a fixed value of the Jacobi constant, the bigger the value of the initial x 0 , the bigger the initial velocity, so the orbits move faster downwards.
For fixed values of n, 0 and a precision ε, from the expression of the total error performed we can obtain an upper bound for the number of Poincaré iterates starting at initial condition x 0 that are admissible in order to maintain the desired precision: But m p grows as x 0 increases, so if we want to keep m p bounded, x 0 must be bounded also. In Table 2 the maximum value x 0 that can be considered in order to maintain a given error when using the extended truncated LRP is given. For bigger values of x 0 , the number of jumps between two consecutive Poincaré iterates is too big, but we are interested in exploring the "intermediate" regions where the linear approximation is not good enough and the invariant tori begin to be destroyed. Next, we show some of the Poincaré map plots for different values of 0 and C ∞ (see Table 1). As we have said, orbits of the Poincaré map display two very different types of behavior depending on the distance to the ring primaries: far of them, the orbits are well approximated by the solutions of the linear system, so the motion is "quite" regular, whereas near the primaries, the dynamics are strongly affected by the presence of the stable/unstable manifolds associated to the equilibria and Lyapunov periodic orbits. The boundary between these regions is formed by the KAM circles closest to the ring primaries. We are interested in the last KAM circle separating the  Table 1 "regular" from the more chaotic region, what we call the innermost invariant curve of the Poincaré map. We can plot the Poincaré map iterates in the configuration space, so the innermost will be the closest invariant curve that intersects |y| = 1/2 and its distance (in configuration space) to the origin will give an approximate idea of the width of the ring: the region where the motion is more unpredictable and chaotic, and we can find orbits which cross into the region near the ring primaries.
Clearly, the location of the innermost will depend on 0 and on the value of the Jacobi constant. The bigger the value of 0 and C ∞ , the closer to the origin the innermost will be. This can be seen in Figs. 7, 8, and 9 where 1,000 iterates for each initial conditions of a grid for the Poincaré map P are plotted. The initial conditions have been taken equally spaced in a fixed region in the configuration space. In each of the plots, we show the Poincaré section map for the two values C a , C b shown in Table 1. First, we notice that the bigger the value of 0 , the smaller the distance from the innermost to the origin (compare the range of values of x 0 shown in each plot). As we have expected, the width of the ring will be bigger for smaller values of 0 . Second, the bigger the value of the Jacobi constant, the smaller the bottleneck around the equilibrium point L 1 , so the innermost invariant curve is closer to the origin.  Table 1 This is more difficult to see for bigger values of 0 (compare Figs. 7 and 9), as the linear approximation gives a good approximation for smaller values of x 0 .
The approximate location of the innermost can be done as follows: fixed a value of the Jacobi constant, consider the initial condition (x 0 , 0, 0, v 0 ), and compute a certain number of Poincaré iterates. If they belong to an invariant curve, then it will be possible to compute its rotation number (following the procedure explained in [9], for example). Varying x 0 (decreasingly), the "last" invariant curve should be detected. It is beyond our scope to compute it with a higher precision, so we have just detected, for each 0 and C J , which is, approximately, the value of x 0 corresponding to the invariant curve. These approximate values are given in Table 3. As we have expected, the value of x 0 decreases as 0 or C J increases.

Dynamics Around the Ring Primaries
Near the ring primaries, the existence of the stable/unstable manifolds associated to the equilibria and the Lyapunov periodic orbits play a key role in the dynamics.  Table 1   Table 3 For each value of 0 in the first column, the approximate value of the x 0 coordinate of the initial condition of the innermost invariant curve is given, for each value of the Jacobi constant C a and C b of As we have seen, the equilibrium point L 1 is of type center × saddle for all the values of the parameter 0 . Thus, a family of periodic orbits (the so called Lyapunov orbits) is born from it, as the Jacobi constant decreases, and the orbits inherit the hyperbolicity (at least, for values of the Jacobi constant close to the value at L 1 ). For each Lyapunov periodic orbit we can follow the unstable and stable branches of the invariant manifolds. For each invariant manifold, one branch goes towards the right hand side (from the Lyapunov orbit point of view), whereas the other branch goes towards the ring primary. In Fig. 10, left, the two branches of the unstable manifold W u associated to a Lyapunov periodic orbit are plotted in the configuration space, up to a certain Following the successive intersections of the right branch of W u and W s with y = −1/2, many heteroclinic connections between the periodic orbit and its translates in the y direction can be found. Hence, orbits on the branch of the unstable manifold moving away from the ring particles, can visit the neighborhood of the translates of these ring particles, see Fig. 10 right, were a heteroclinic connection between two different Lyapunov periodic orbits is plotted. This happens for different values of the Jacobi constant (the heteroclinic connections belong to families of orbits), and for the different values of 0 explored. In all the cases, we can see that the orbits of the invariant manifolds expand up to the location of the innermost invariant curve that separates the regular from the chaotic region. Also, for those values of 0 for which the equilibrium point L 2 is hyperbolic, the closure of its 2-dimensional invariant manifold contains the first KAM invariant circle, as we expect (see [1]).
With respect the left branches of the invariant manifolds, following them it is possible to obtain heteroclinic connections with the Lyapunov orbits on the x < 0 half plane. A rigorous study using symbolic dynamics would allow one to show the existence of orbits with prescribed itineraries, visiting the different ring primaries (see for example [2]) in the case of the planar restricted three-body problem). Therefore, we can say that around the ring primaries and until the innermost invariant curve, the dynamics is very rich principally due to the presence of the invariant manifolds of the Lyapunov orbits around L 1 , independently from the value of 0 . The only qualitative difference that can be observed, as we observed in the Poincaré section maps, is the location of this last KAM invariant tori, or equivalently, the width of the chaotic region around the ring primaries that we can consider as the "ring".

Conclusions
The numerical and analytical analysis presented above for the Limiting Ring Problem (LRP) indicate that it contains a diversity of different dynamical behavior. Like many Celestial Mechanics problems, the LRP provides fertile ground for testing different analytical and numerical procedures. It is clear that there is a great deal more structure and dynamics present, including periodic orbits and symbolic dynamics, that could be described. These might be dynamically interesting orbits that "tour" the ring in interesting ways.
As noted, the LRP is just a cartoon model of an actual ring. It is interesting that the dynamics alone provides a way to describe a natural boundary to a swarm of infinitesimal particles forming a diffuse ring engulfing an n-gon of particles, distinguishing those whose orbits are dominated by the n-gon from those away from the ring. That this width is very much less than that observed for thin planetary rings (see [1]), is another indication of importance of external shepherds and resonances in the maintenance of these rings. The inclusion of these external bodies or structure to the ring (ellipticity or thickness) seems quite challenging using our techniques.
One can imagine that the observation and study of the orbits of individual ring particles within the swarm debris making up observed planetary rings is not that far off. The actual orbits of these bodies will surely be at least as complicated and at least as interesting as the dynamics of the LRP.