Motion around a monopole + ring system

Motion around a Monopole Ring system:
I. Stability of Equatorial Circular Orbits vs
Regularity of Three-dimensional Motion.

Javier Ramos-Caro, Juan F. Pedraza and Patricio S. Letelier
Departamento de Matemática Aplicada, IMECC, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil
Department of Physics, University of Texas, 1 University Station C1608, Austin, TX 78712, USA
E-mail: javier@ime.unicamp.brE-mail: jpedraza@physics.utexas.eduE-mail:

We study the motion of test particles around a center of attraction represented by a monopole (with and without spheroidal deformation) surrounded by a ring, given as a superposition of Morgan & Morgan discs. We deal with two kinds of bounded orbits: (i) Equatorial circular orbits and (ii) general three-dimensional orbits. The first case provides a method to perform a linear stability analysis of these structures by studying the behavior of vertical and epicyclic frequencies as functions of the mass ratio, the size of the ring and/or the quadrupolar deformation. In the second case, we study the influence of these parameters in the regularity or chaoticity of motion. We find that there is a close connection between linear stability (or unstability) of equatorial circular orbits and regularity (or chaoticity) of the three-dimensional motion.

Planets and satellites: rings, celestial mechanics, chaos.

1 Introduction

Light particles interacting with a central massive body are frequently encountered in diverse fields of physics. Electrons in rotating molecules and particular versions of the many-body problem in celestial mechanics are the most well-known examples, but some nuclear models fall into the same category. Increasing amounts of information about narrow planetary rings suggest that such rings are often associated with the so-called shepherd satellites (Goldreich & Tremaine 1979 ), and may exist due to mechanisms somewhat more complicated than the well known broad rings (Greenberg & Brahic 1984 ; Smith et al. 1981 ; Smith et al. 1982 ; Murray et al. 2005 ).

Recent progress in the subject suggest a generic mechanism, that does not depend on Kepler orbits, to explain the formation of rings around a rotating object which also holds for systems quite different from planetary rings. In rotating scattering systems, the generic saddle-center scenario leads to stable islands in phase space. Non-interacting particles whose initial conditions are defined in such islands will be trapped and form rotating rings. This result is generic and guarantees that the orbits supporting the ring structure are rather insensitive to small perturbations and thus may play a role in different situations of the type mentioned above (Benet & Seligman 2000 ). So, although we are interested in studying planetary rings, there are many other systems that exhibit the same structure and therefore one can perform similar treatments by knowing the nature of the forces involved.

Now then, based on the arguments mentioned above, one may conjecture that once the ring structure is formed around the rotating body, it remains stable, as revealed by the structure of its phase space. In this set of papers, we propose an analytical study about the dynamical aspects concerning the stability of ring structures and the relation with the existence of isolated islands in the phase space. In this first paper we focus in the effects related to the physical parameters like the mass of the central body, the geometry of the central body, the mass of the ring and the size of the ring. In next papers, we plan to examine the influence of the angular momentum of the system in its stability, as well as, large perturbations like the interaction with external satellites and their responses to resonance.

Planetary rings consist in thin discs of cosmic dust and other small colliding particles revolving around a central planet in a flat disc-shaped region. The most spectacular example of ring structures are those around Saturn (Porco et al. 2005 ; Cuzzi et al. 2010 ), but they are a common feature of the other three gas giants of the solar system; Jupiter, Uranus and Neptune possess ring systems of their own. Recent reports have suggested that the Saturnian moon Rhea may have its own tenuous ring system, which would make it the only moon known to possess a ring system (Jones et al. 2008 ).

There are many possible mechanisms to explain the existence of planetary rings, but essentially three of them are the most relevant: from material of the protoplanetary disc that was within the Roche limit of the planet and thus could not coalesce to form moons; from the debris of a moon that was disrupted by a large impact; or from the debris of a moon that was disrupted by tidal stresses when it passed within the planet’s Roche limit. This last one allows us to predict that Phobos, a moon of Mars, will break up and form into a planetary ring in about 50 million years due to its low orbit (Holsapple 2001 ).

Additionally, we should take into account the increasing amount of data we get from extrasolar systems. The discovery of extrasolar planets (which up to date are more than 300) by radial velocity measurements has provided the first dynamical characteristics of planets: orbital elements and mass. The next step will be to investigate physical characteristics: albedo, temperature, radius, etc, and their surroundings. Among the latter are planetary rings. The emitted thermal infrared light from the planet should show no phase effect assuming the planet is in thermal equilibrium. But the reflected visible light will vary with phase angle, as should be shown by a broad-band photometric follow-up of the planet during its orbital motion. In particular, it has been studied from different perspectives how the presence of a ring around a planet would influence its brightness as a function of its orbital position and based on that there have been multiple theoretical predictions for photometric and spectroscopic signatures of rings around transiting extrasolar planets (Barnes & Fortney 2004 ; Arnold & Schneider 2006 ; Ohta, Taruya and Suto 2009 ). So, understanding the dynamics of composed planetary systems would not constitute just an astrophysical curiosity, it would have some impact on the strategy for their detection (see Arnold and Schneider 2004 and the references therein).

On the other hand, it has been shown that collisions play an important role in the dynamical aspects of the rings (Longaretti 1992 ; Sicardy 2006 ). Under some physical assumptions, one can see that as this process dissipates mechanical energy, while conserving angular momentum of the ensemble, tend to flatten the disc perpendicular to its total angular momentum, and to circularize the particles orbits. If the planet is not spherically symmetric, as it is the case for all the oblate giant planets, apart of being flat the configuration becomes a equatorial ring. Of course, some physics is still missing in the description of the rings and an interesting question here would be if these effects that one does not take into account as a first approximation are relevant in the dynamics of the ring. A classical approach is to consider what happens to the system when small perturbations are applied.

Commonly, disc-like systems such as planetary rings are described by a set of coupled equations: two hydrodynamic equations (Euler’s and continuity equation), Poisson’s equation and a equation of state (Toomre 1964 ; Binney & Tremaine 2008 ). For small enough disturbances, these equations can be linearized and solved under further simplifications as there are still missing analytical self-consistent models of planetary rings. Nevertheless, it is possible to obtain useful information of this by using simplified models which illustrates quantitatively how pressure and rotation tend to stabilize the disc against self-gravity.

Although most rings were thought to be unstable and to dissipate over the course of tens or hundreds of millions of years, recent evidence coming from NASA’s Cassini-Huygens Mission suggest that Saturn’s rings might be quite old, dating to the early days of the Solar system. Then, the exact mechanism to explain which physical factors account for the stability of these systems is still an open question to be answered.

Our approach to study the stability of rings will be quite different of the mentioned above. The method we will use is insensible to the hydrodynamical aspects of the rings but with the advantage that we will have under control the geometrical aspects of the models, namely the size and shape of the ring, the mass quotient between the planet and the ring and deviations from the spherical geometry of the planet. Also, similar techniques can be applied to study the effects caused by the angular momentum of the planet and even the influence of external satellites, but these topics are out of the aim of this paper and will be left for future works. Another benefit of this approach is that it becomes easy to elucidate a non-trivial connection between the stability in equatorial circular orbits and the regularity in three-dimensional motion. In particular, we will see how the existence of isolated islands will lead to stable ring configurations and how relevant are the physical parameters mentioned above.

Now then, the exact gravitational potential of a ring of zero thickness and constant linear density is given in its exact form by an elliptic integral that is seldom used for practical purposes. This potential is usually approximated by truncated series of spherical harmonics, i.e., a multipolar expansion (see for example the deep analysis of orbits performed in Tresaco & Ferrer 2010 ). As far as we know there are not simple expressions for the exact gravitational potential of a finite flat ring with inner and outer edges and any surface density.

Meanwhile, simple potential-density pairs for thin discs are known. The simpler is the Plummer-Kuzmin disc (Kuzmin 1956 ) that represents a simple model with a concentration of mass in its center and density that decays as on the plane of the disc. This structure has no boundary even though for practical purposes one can put a cutoff radius wherein the main part of the mass is inside, say 98% of the mass. There are many other models in the literature (Lynden-Bell 1962 ; Mestel 1963 ; Toomre 1963 ; Hunter & Toomre 1969 ; Kalnajs 1972 ; Jiang 2000 ; Jiang & Moss 2002 ; González & Reina 2006 ; Jiang & Ossipkov 2007 ; Pedraza, Ramos-Caro & González 2008 ), some of them of infinite extension but others with an outer edge. Among the last ones, a family of simple models are the Morgan and Morgan discs (Morgan & Morgan 1969 ), which has been inverted (Lemos & Letelier 1994 ) in order to obtain infinite discs with a central hole of the same radius of the original disc. We can also put a cutoff in these structures and therefore the inverted Morgan and Morgan discs can be considered as representing a flat rings. These models have been used by the authors in order to study the superposition of an annular disk with a central black hole in the context of the General Theory of Relativity (see also Semerák & Suková 2010 ; Lora-Clavijo, Ospina-Henao & Pedraza 2010 , for recent related works).

Another approach to construct flat ring structures is by the superposition of different kind of discs. In Letelier 2007 the author superposed Morgan and Morgan discs (Morgan & Morgan 1969 ) while in Vogt & Letelier 2009 the Kuzmin-Toomre discs (Kuzmin 1956 ; Toomre 1963 ) were used. In this paper we will use the first of the above models which are relatively simple and superposing the potential of the planet will allow us to study the effects of those parameters we have been talking about.

2 Models of rings around a spherically symmetric object

In this section we shall focus on simple analytical models that represent sources conformed by an axisymmetric thin ring and a central object that we will assume with spherical symmetry, i.e. a body characterized only by its monopolar moment. The effects concerning with oblate or prolate deformation of the center of attraction, i.e. the inclusion of a quadrupolar moment, will be studied in the last section.

In a recent work, Letelier 2007 , shows how to construct a simple family of potential-density pairs for flat rings by means of the superposition of Morgan and Morgan discs (Morgan & Morgan 1969 ). These discs are finite in extension and have a well-behaved surface mass density with a maximum in the center and monotonically decreasing up to the edge. The structures obtained by the superposition of discs with different densities have a finite outer radius and zero density on their centers, i.e. discs with a hole in their centers, or in other words flat rings. Although the models do not have an inner edge, for practical purposes one can put a cutoff radius and neglect the residual density (which becomes smaller for higher members of the family).

The Morgan & Morgan discs are obtained by solving the Laplace equation in the natural coordinates to represent the gravitational potential of a disc-like structure, i.e. oblate coordinates (defined in the ranges and ) that are related to the usual cylindrical coordinates by


where is a positive constant defining the disc radius. The inverse relations are given by


Note that, according to (1), on the equatorial plane one has to distinguish between two regions: (i) the points inside the disc, with coordinates and ; (ii) the points outside the disc, with coordinates and .

The mass surface density of each disc (labeled with the positive integer ) is given by


where is the total mass and the disc radius. Such mass distribution generates an axially symmetric gravitational potential, that can be written as


Here, and are the usual Legendre polynomials and the Legendre functions of the second kind respectively, and are constants given by


where G is the gravitational constant. If we consider discs of the same radius and decreasing mass,


where is a constant that will be taken equal for all discs of the Morgan and Morgan family, we obtain discs with surface density


In order to obtain a mass density in agreement with a flat ring distribution, Letelier 2007 considered the superposition

Figure 1: The mass surface density for the first ten members of the ring family constructed as a superposition of Morgan and Morgan discs.

We have that the above superpositions give rings of radius and a central residual density that becomes smaller for larger . For practical purposes one can put a cutoff inner radius , which we will take as the radius such that the density falls below the 1% of its maximum. In table 1 we show the values of for the first ten members of this family and, in figure 1, their corresponding mass surface densities. The total mass of each ring is


Thus by increasing , the mass of the flat ring decreases along with the size of the region where it is distributed (near the outer radius). Note that this is a feature associated to the particular family of models we are considering, not with the corresponding physical applications. The values of and are inferred from basic data corresponding to a special case. For example, to describe a ring with mass Kg, inner radius Km and outer radius Km, we have to set and , leading to to a maximum mass density of about (these values agree roughly with the measurements of the so-called Ring A of Saturn (Dougherty et al. 2009 )).

The potentials associated to these structures can be obtained using a superposition with the same coefficients as the ones used for the densities, i.e.,


where is the same as with the masses given by equation (6).

1 0.06210 2 0.23292 3 0.36991 4 0.43940 5 0.54302 6 0.59222 7 0.64335 8 0.67884 9 0.70798 10 0.73231
Table 1: Values of the ratio for the first ten members of the family of rings given by the mass surface (8). Here represents the inner cut-off radius and is the outer radius. For higher members of the family the mass concentration tends to be located near the outer edge.

Finally we add a monopolar term, which represents the exterior field of a central spherically symmetric object, to the ring potential described above. Thus, the total gravitational potential reads


where is the mass of the central object. Although our 4-parametric toy model is quite simple, we believe that this is a good starting point to study the effects on the dynamics caused by the mass ratio and the radius ratio which is one of the purposes of this paper. The effects caused by the rotation of the central body (its own angular momentum would be an additional parameter), will be studied in a next paper through the post-Newtonian scheme.

3 Motion of Test Particles

Let us deal with the problem of motion of test particles around the models described above. Since is static and axially symmetric, the specific energy and the specific axial angular momentum are conserved along the particle motion, which is restricted to a three dimensional subspace of the phase space. Each orbit is determined by the set of equations (Binney & Tremaine 2008 )


where is the effective potential, defined by


in terms of which the particle’s energy can be written as


Relations (12a)-(12b) define an autonomous system whose equilibrium points are and , where must satisfy the equation


that is the condition for a circular orbit in the plane . In other words, the equilibrium points of the system occur when the test particle describes equatorial circular orbits of radius , specific energy , and specific axial angular momentum given by


The description of circular orbits in the equatorial plane is a first step to understand the linear stability of the system, as well as, the regularity or chaoticity of three dimensional orbits. On the one hand, if one assumes as a first approximation that the structures are built from particles moving only in circular orbits, the epicyclic and vertical frequencies of quasi circular orbits provide us a criterion for the system’s stability (see section 4). On the other hand, the analysis of such frequencies also leads to determine the existence of saddle points, which are preliminary indicators of irregular motion.

In order to deal with the problem of correlation between regularity of three dimensional motion and the stability of circular orbits, we have to distinguish between exterior and interior motions of test particles, separately. As we pointed out in the last section, the relation between oblate and cylindrical coordinates is given by (2) or, in a more explicit form




for points located outside the disc zone, while for points inside the disc region, we have


The ambiguity in the sign in the last equation is due to the singular behavior of the coordinate on crossing the disc. Rigorously speaking we have at (upper limit) and at (lower limit).

For later formulae, it is important to point out that the piecewise form of this transformation rule leads also to a piecewise form in the first and second derivatives, when evaluated at the equatorial plane. After some calculations, we obtain the following relations for the first derivatives of and , at :




The ambiguity in the sign in (21) has the same meaning as in (17). These relations are used to compute the second derivatives in the equatorial plane and the result is




The latter two equations will be important in the calculation of epicyclic and vertical frequencies (see section 4).

By introducing (21) in the equations of motion (12b), when evaluated in the inner zone, we can derive the equations of motion for a test particle in the case in wich , . For we have


Note that, in the derivation of this equation, there is no difference between the choice along with , and the other option, along with . This is due to the fact that is an even function of . In contrast, the equation for has a change of sign on crossing the disc:


However, since has symmetry of reflection with respect to the plane , its -derivative must vanish exactly in the equatorial plane and we have


ensuring the existence of circular orbits inside the disc. Now then, according to (25), one can verify that equation (16), for inner equilibrium points, can be cast as


On the other hand, for outer circular orbits (equilibrium points outside the disc) equation (16) becomes


Equations (28) and (29) are relevant in the derivation of the quadratic epicyclic and vertical frequencies, in order to provide a criteria for the linear stability of the structures studied here.

4 Linear Stability of Structures

As it was mentioned above, we assume the simplified model that the structures are built from particles moving in concentric circles. In this first approximation, the stability analysis of circular orbits associated to test particles provides a stability criterion for the structure, assumed to be a rotating ring of fluid (Lord Rayleigh 1916 ; Landau & Lifshitz 1987 ; Letelier 2007 ). For this reason, we now examine the behavior of the epicyclic and vertical frequencies associated to quasi circular orbits. These quantities describe the response of test particles to radial and vertical (-direction) perturbations, when describing a circular motion. The epicycle frequency and the vertical frequency , can be calculated from the effective potential through the following relations (Binney & Tremaine 2008 ):


If relation (16) is introduced in the second derivatives of , we can obtain and as functions of . Thus, values of such that and (or) corresponds to stable quasi circular orbits under small radial and (or) vertical perturbations, respectively. Otherwise we find unstable circular orbits. Since the equatorial plane is composed by two regions, i.e. inside and outside the disc, we must define each of these quantities as piecewise functions. The quadratic epicyclic frequency for equilibrium points in the inner zone can be written as


whereas that for equilibrium points in the outer zone we have


Similar expressions can be derived for quadratic vertical frequency, which inside the disc takes the form


and outside the disc is


Equations (31) and (33) are used to determine the conditions to be satisfied by the parameters of the system, so that the stability condition is fulfilled (the relations (32) and (34) are studied in the next section). By keeping and constants, we can search the set of values that make stable configurations. From here on we shall use without loss of generality. In Figure 2, at the left side, we can see the behavior of (inside de disc) for different values of the planet’s mass, in the case . Since it is a positive concave function with a critical point in the range , we can find the minimum value of (or the maximum rate ) for which is positive in this range. We find such value by solving the simultaneous equations

for the variables and , representing the minimum value of the planet’s mass and the corresponding critical radius, respectively. We can see that by increasing , starting from , we obtain increasing values for .

Figure 2: Behavior of quadratic epicyclic frequency for the model , for . For values larger than (bottom curve on the left and right sides), is negative in the equatorial plane. For smaller values, is greater and the gap of discontinuity at is smaller.

Figure 3 shows the behavior of the quadratic vertical frequency, also for the model , and we note that it is a monotonically decreasing function, in the range , with a minimum at . In this case, to find the minimum value of (which we shall denote as ) for which is positive, it is enough to solve the equation

for the variable . As in the previous case, we note that by increasing , starting from , the values of increase.

Figure 3: Behavior of quadratic vertical frequency for the model , for . For values larger than (bottom curve on the left side and top curve on the right side), is negative in the equatorial plane. For smaller values, is greater and the gap of discontinuity is smaller.

Another feature showed in figures 2 and 3 is that, for large values of , the behavior of and in the range is very different from the behavior in . In general, we see that

(both are finite values), but the difference between these two limits is attenuated by decreasing the ratio . We say that there is a gap of discontinuity at . The same statements hold for and models with . It is clear that the gap disappear as .

The behavior of epicyclic and vertical frequencies sketched above is the same for the other models , so that we can compute the maximum values of the rate leading to or . Such values are listed in Table 2, for the first ten members of the family (11). Since in all cases , we can establish that is the parameter that determines the boundary between stability and instability in configurations characterized by gravitational potentials of the form (11).

1 2.83102 0.04850 2 0.86297 0.02141 3 0.43398 0.01205 4 0.26554 0.00772 5 0.18068 0.00537 6 0.13148 0.00395 7 0.10024 0.00303 8 0.07909 0.00240 9 0.06407 0.00194 10 0.05299 0.00161
Table 2: Ratios and for the first ten members of the family of configurations represented by (11). Here and are the minimum value of the central body’s mass such that and , respectively. In all cases .

It is important to note that, in these models with fixed exterior radius , the smaller the size of the -th ring model, the larger is . The reason is that, for higher members of the family (11), the ring’s mass concentration tends to be located near the outer edge and, therefore, away from the central monopole. Thus, it will be required increasingly central body mass to provide stability to the ring structure, as grows. This fact can be glimpsed through a comparison between tables 2 and 1, from which one might infer that there some kind of correlation between and the ring’s size of the family of models. In figure 4 we show this correlation by plotting the values of tables 2 and 1 and interpolating the corresponding points, for the first ten members of the family. Thus we can sketch a boundary between stability and instability of self-gravitating configurations under study. In figure 4, points located in the grey zone corresponds to parameters leading to stable configurations, while the points in the white region are associated to unstable configurations. We find that the points belonging to the sepparatrix of fig. 4 can be fitted by the relation

Figure 4: Correlation between the rate and according to the behavior of quadratic vertical frequency and quadratic epicyclic frequency . The interpolation line among points is a separatrix between a stability grey region, where and , and an instability white region, where (and may be positive or negative). The points plotted here correspond to the first ten models, from left to right.

5 The Phase Space Structure

In this section, we shall make a description of three-dimensional motion through the study of phase space structure associated to orbits of test particles. We are principally interested in the influence of the mass ratio, in the regularity of three-dimensional bounded orbits. The influence of the size of the ring can be inferred from it, because in the models studied here, is automatically determined by . This influence has already been analyzed in relation with the stability of equatorial circular orbits and now we extend such study to the case of more general orbits (remember that the motion restricted to the equatorial plane, which is completely integrable, can be exclusively classified as stable or unstable). In order to show how the nature of bounded motion is conditionated by the linear stability of the self-gravitating structures, we will use parameters close to the critical values we have defined previously (table 2). The influence of the quadrupolar deformation is considered in section 6.

We shall use cylindrical coordinates and plot surfaces of section corresponding to the equations of motion (12a)-(12b), for different values of the conserved quantities and . It is worth to point out that we need to make explicit distinction between the orbits that cross the plane at and the ones that cross it only at . The former are called disc-crossing orbits (DCO) and the latter we shall denote as non disc-crossing orbits (NDCO). The reason is that, for each of the above situations, the origin of irregular motion is different. For the case of DCO the discontinuity in the -component of the gravitational field (equations (26a)-(26b)) can produce a fairly abrupt change in the curvature, leading to irregular motion (Saa 1999 ; Saa 2000 ; Hunter 2005 ; Ramos-Caro, López-Suspez & González 2008 ). Somehow, this is a problem analogous to the case of the chaotic behavior of Chua’s circuit (Matsumoto, Chua & Komuro 1985 ), which is described by an autonomous system of the type , where is a piecewise function of class (continuous but not differentiable). It is the first example where the existence of such class of function leads to the existence of a chaotic attractor in a dynamical system (Madan 1993 ).

In contrast, in the case of other three-dimensional orbits, chaotic motion is due to the existence of saddle points in the effective gravitational potential. Such saddle points, in the particular case of a potential as , will be located at equatorial plane, outside the disc (note that it is not possible to define saddle points inside the disc, due to the discontinuity in the potential’s -derivative).

Equations (32) and (34) help us to investigate the existence of saddle points in the potential , through the evaluation of the quantity , which is negative when evaluated at saddle points. Since is symmetric with respect to , the term vanishes at equatorial plane and the condition for existence of saddle points is reduced to . In the right side of figures 2 and 3 we can observe the behavior of and , respectively, and the product between them is showed in figure 5, for the case of model . For different values of we note that there is a region of saddle points closed to the outer edge, even for the maximum value leading to a stable configuration. When the ratio decreases and we have even more stable configurations, the range of saddle points decreases as well as the gap between the values of near the outer edge. This suggest that more and more stable configurations, produce less and less irregular orbits.

Figure 5: Behavior of the product for the values (curves from top to bottom). For , the maximum value leading to a stable configuration, there is a small region of saddle points near to the outer edge (right side of the dashed line). This region becomes smaller as decreases.

In order to illustrate the above statements, we plot the surfaces of section corresponding to motion around the first two models, by choosing several values for the parameter . We solve the equations of motion (12a)-(12b) by the Runge-Kutta 4th-method with variable time step and incorporating the Hénon algorithm (Hénon 1982 ), in the LCP (Laboratório de Computação Paralela) at IMECC. In the algorithm we take into account explicitly the discontinuity in the field force by implementing the set of equations (26a)-(27), as well as, the piecewise transformation (17)-(20). We choose initial conditions at , and several values for (the component is given by (14) and does not vanish). Orbits were integrated with a precision characterized by a relative error of or less, in the energy conservation. Integrations were carried out on times of the order of .

In figure 6 we have set , which determine a radially stable configuration, but vertically unstable. There is a variety of central KAM curves and small resonant islands outside the disc zone as well as in the inner region. They alternate with two chaotic regions, one of them is due to DCO (the most prominent) and the other, more dense, is the result of two orbits near the saddle point ( for this case). Here we have a situation where an unstable model admits a variety of irregular orbits.

Figure 6: Surface of section for some orbits around the model , with . Such a rate determine an unstable situation, where (marginally) but . We find two chaotic regions: (i) A prominent zone due to DCO and (ii) a smaller zone corresponding to two orbits with initial conditions next to the saddle point. In this case, we have chosen and .

If we choose a more stable configuration, i.e. characterized by parameters near the sepparatrix of figure 4, we obtain surfaces of section as the shown in figure 7 where the chaotic region is restricted to the small central annulus enclosing the three KAM curves on the rigth. Note the presence of two chains of resonant islands enclosing the stochastic zone and other two small chains within it. In this case we find no chaotic motion associated to disc crossings, but only to exterior saddle points.

Figure 7: Surface of section corresponding to three-dimensional orbits around the model , setting . The model represents an unstable situation (although very near to the separatrix of figure 4) and there are some chaotic orbits. The irregularity associated to these motions is due to the proximity to saddle point near to the disc edge. The orbits plotted have values and .

When we turn to the other side of the border, toward the region of stable configurations in figure 4, regularity is the most common feature of three dimensional orbits. Such is the case illustrated in figure 8, where we find only smooth KAM curves along the surface of section, even in the case of DCO orbits. In this case we can see a chain of 29 small resonant islands enclosing another chain of 3 central islands.

Figure 8: By maintaining the same values for and as in figure 7, but decreasing to , we find a situation when the model is stable (also very near to the separatrix of figure 4) and we note only regular orbits.

This transition from chaos to regularity, characterized by the apparition of increasingly number of islands when we pass from unstable to stable structures, can be seen in figures 9 and 10, which correspond to the model . In the former there is a sequence of three prominent central islands, one of them enclosing a small chaotic region with a chain of 17 small resonant islands. In figure 10 the stochastic region is absent and we see only a dashed KAM curve enclosing three large islands.

It is clear that we expect a similar behavior for the remaining members of the infinite family: the increment of the mass of the spherical central body favors the existence of isolated islands in the phase space.

Figure 9: Surface of section corresponding to three-dimensional motion around the model with , for orbits with and . Since we are dealing with a configuration belonging to the unstable region of figure 4, there is a stochastic zone between a variety of regular KAM curves.
Figure 10: By maintaining the same values for and as in figure 9, and choosing , we obtain a situation when the model is stable and we note only regular orbits.

6 Addition of Quadrupolar Term to the Central Body

It is known that in the universe there exist many centers of attraction with a certain deviation from spherical symmetry, so we are interested in including a quadrupolar term in the description of the central object that makes up the structures studied above and examine the effects of deformation (preserving the axial symmetry). Now the gravitational potential of the structure is


where is the quadrupolar moment that quantifies the oblate () or prolate () deformation of the central body. In general, it is related to the mass density (axially symmetric) through the equation (Binney & Tremaine 2008 )


where and are the spherical coordinates.

The equations of motion for test particles around such configurations are the same (12a)-(12b) but replacing by . Then, the relation for inside the ring becomes


while the relation for , in the inner, remains the same as in section 3. As a consequence, the relation that determines , for inner circular orbits, changes to


and, for outer circular orbits


Therefore, the relations for the quadratic epicyclic frequency turns to




while for the vertical frequency we have now



Figure 11: Behavior of for and (curves from top to bottom). Maintaining the central monopole fixed, the epicyclic frequency grows as the quadrupolar moment decreases.
Figure 12: Behavior of for and (curves from top to bottom). Maintaining the quadrupolar moment fixed, the epicyclic frequency grows as the central mass increases.

We note that the epicyclic frequency behaves in a very similar fashion as in the previous case of a spherically symmetric central body. This can be seen in figures 11 and 12. In the former we have chosen a certain fixed value for and plotted for some positive values of quadrupolar moment (assuming prolate deformation), while in the latter we have fixed and plotted for different values of . These figures reveal that the quadratic epicyclic frequency is a positive concave function with a critical point in the range . Needless to clarify that it is a monotonically decreasing function if we use negative values for . Also, we have to point out that this quantity has the same behavior for the remaining models (the same statement holds for the vertical frequency).

Figure 13: Behavior of for and (curves from top to bottom). Maintaining the central monopole fixed, the range where the quadratic vertical frequency is positive, grows as the quadrupolar moment decreases.
Figure 14: Behavior of for and (curves from top to bottom). Maintaining the quadrupolar moment fixed, the range where the quadratic vertical frequency is positive, grows as the central mass increases.

In contrast, the addition of a quadrupolar term introduces new features in the behavior of the quadratic vertical frequency. It can be observed from figures 13 and 14 that is now a negative concave function with a maximum in the range , for positive values of quadrupolar moment (we have used the same values for the parameters as in figures 11 and 12). It is worth clarifying that, by using negative values for , the quadratic vertical frequency becomes a monotonically decreasing function, as in the case of a spherically symmetric central body.

According to the above statements, the task to find the limiting values of and , leading to stable configurations, is reduced to formulate the simultaneous equations

for the variables and , i.e. the minimum value of the central body’s mass and the maximum value of the quadrupolar moment such that the quadratic vertical frequency (and evidently ) is positive in the range between the cut-off radius and the outer radius. In table 3 we have listed the corresponding values of these quantities for the first ten members of the family. Consequently, we show the correlation between the logarithm of and , plotting the separatrix between the stability and instability (fig. 15). Here, represents the quadrupolar moment of the -th ring model and it is computed using equation (37). We find that the variables and can be fitted by the relation


providing an approximate expression for the separatrix of figure 15.

1 -16.1076 0.04804 2 -0.55925 0.02024 3 -0.12408 0.01040 4 -0.04713 0.00602 5 -0.02301 0.00378 6 -0.01300 0.00253 7 -0.00807 0.00177 8 -0.00537 0.00129 9 -0.00375 0.00097 10 -0.00273 0.00074
Table 3: Ratios and for the first ten members of the family given by (36). Here and represents the minimum value of mass and maximum value of quadrupolar moment such that the structure is linearly stable.
Figure 15: Correlation between the rate and . The grey region represents the set of values leading to stable structures, while the white region corresponds to the stable ones.
Figure 16: Detail of the central part of figure 6 corresponding to , and . The inner stochastic region is generated by two orbits with initial conditions near the saddle point. The outer stochastic region is due to a DCO.
Figure 17: We have used the same parameters and initial conditions as in figure 16, but switching on the quadrupolar moment to . The stochastic regions associated with this unstable configuration are distorted, with respect to the figure 16. The outer chaotic zone of DCO is absent, while the inner chaotic zone is more prominent.

As it was shown in Guéron & Letelier 2001 , for the case of monopole-quadrupole configurations, the chaoticity is induced by prolate (rather than oblate) deformation. In the case studied here, we have also a ring structure with oblate shape and one would expect two situations: (i) an “attenuation” in the chaoticity, for the case of a central body with prolate deformation; (ii) only regular motions, if the central body has also oblate shape.

The surfaces of section corresponding to structures located at stability region of figure 15 are not very different from those shown in plots 8 or 10. However, some differences appear when we deal with unstable configurations. For example, the outer stochastic region due to DCO, in figure 6, disappears when we turn on the quadrupolar moment, but the inner chaotic region is now most prominent and shifted slightly to the left (details are shown in figures 16 and 17).

7 Concluding Remarks

Simple models as the ones described by relations (8) and (10), provide us an useful tool to analyze and understand the dynamics of astrophysical objects that can be modeled as a spherical (or quasi-spherical) mass surrounded by a flat ring structure. One important and fundamental aspect in the dynamics is the stability against proper modes, if we consider in a first approximation that the ring structure is formed principally by particles moving in concentric circular orbits. It is known, in previous studies (Letelier 2007 ), that the structure of one ring standing alone has limited stability but, by adding a central monopole, it can improve. In this study we were able to verify such statement and find the set of values for the parameters leading to linearly stable configurations. Interestingly, we find monopole-ring configurations belonging to the stability region of figure 4, when dealing with parameters of the order of physical measurements performed in the solar system. For example, a structure with the dimensions of Saturn has