2 Orbits in non-axisymmetric 2D potentials

  
2.1 Third integral and the pendulum model

In any 2D non-axisymmetric potential the main families of orbits are the so-called loop and box (BT87). Which family dominates the orbital structure depends mainly on the relative value of the rotational kinetic energy with respect to the degree of flatness of the potential (see below). To describe the problem in a more general context, let us consider a spatial axisymmetric galaxy, where we assume that the potential depends on the position through $m_{\rm q}(R,z)=R^2+z^2/q^2$, $(R,\varphi,z)$ being cylindrical coordinates and q<1 the semiaxis ratio of isopotential curves in the Rz-plane. The potential is then $\phi(\vec{r})=\Phi(m_{\rm q}(R,z))$ where $\Phi$ is a smooth function of its argument.

In any case, 3D motion reduces to 2D motion in Cartesian coordinates if we introduce the effective potential $\Phi(m_{\rm q}(R,z))+p^2_{\varphi}/2R^2$, where $p_{\varphi}$ is a global integral. As the second term is the same for any $\Phi$, we focus the attention on the motion of a star in the 2D potential $\phi(x,y)=\Phi(m_{\rm q}(x,y))$ where x,y are coordinates in some meridian plane by setting $p_{\varphi}=0$. Alternatively, $\phi(x,y)$ could represent the motion in the equatorial plane (z=0) of a barlike galaxy, being then x,y coordinates in the latter plane.

The equations of motion, in these variables, are

$$\begin{eqnarray*}\dot{p}_x=-2\Phi'x,\qquad \dot{p}_y=-2\Phi'y/q^2, \end{eqnarray*}$$

where $p_x=\dot{x}$, $p_y=\dot{y}$ and $\Phi'\equiv{\rm d}\Phi/{\rm d}m_{\rm q}$is assumed analytic everywhere. To be $-\vec{\nabla}\phi$ a well defined gravitational field it is necessary to impose the conditions $\Phi'>0$ and $\Phi''<0$. To understand the differences between both families of orbits one can follow different approaches. The "rigorous'' one as follows. For box orbits, one should restrict the flow to the invariant planes p_x=x=0or p_y=y=0 and to investigate the 1D Hamiltonians H_y and H_x, respectively. Take one of them and consider the other as a small perturbation. The next step is to analyse the stability of the periodic orbits in the unperturbed 1D Hamiltonians, at a given energy level, by a linearization of the equations of motion of the full Hamiltonian around these orbits (for instance by a Floquet analysis, see BT87 and Sect. 5). Similar considerations apply to loop orbits: just take values of q very close to 1 so that the field is nearly spherical and can be written as a near integrable one. The stability of the 1:1 (circular) periodic orbit is then analysed. This procedure is followed, for example, by MS89 and PL96 for the 2D logarithmic potential to conclude that, for the energies and values of q they studied, the short-axis periodic orbit (y-axis) is, in general, unstable while the long-axis orbit (x-axis) is, in general, stable for low-to-moderate energies (see Sect. 5). The 1:1 periodic orbit (that bifurcates from the y-axis orbit) turns out to be always stable for any physical value of q. Therefore, box orbits can be thought as perturbations to the x-axis periodic orbits while loop orbits arise from perturbations to the 1:1 (circular) periodic orbit in the spherical system.

A "physical'' interpretation is the following. The angular momentum (or the rotational kinetic energy) plays a crucial role in the existence of both families of orbits. Indeed, take polar coordinates in the xy-plane: $x=r\cos\theta$, $y=r\sin\theta$, so $m_{\rm q}\rightarrow m_{\alpha}= r^2(1+\alpha\sin^2\theta)$, where $\alpha=(1-q^2)/q^2$. Due to the lack of central symmetry, a test star will be acted by a torque $N= -\partial\phi/\partial\theta= -\Phi'\partial m_{\alpha}/\partial\theta = -\alpha r^2\Phi'\sin 2\theta$. If $\Phi'\ne 0$ then for any r>0the torque is null at $\theta=0$, $\pi/2$ (and $\pi$, $-\pi/2$), that is, on the x and y axis. Recalling that $\dot{p}_{\theta}=N$, where $p_{\theta}$ is the angular momentum of the star, we conclude that an orbit with $p_{\theta}=0$ will follow a rectilinear orbit along the x or y axis. A simple inspection of the expression for N shows that the torque is negative in the first and third quadrant, being positive in the others. Therefore, we see why the x-axis periodic orbit is stable while the y-axis one is unstable. The torque confines near the x-axis and pulls away near the y-axis. On the other hand, a simple epicycle approximation (see BT87, Lees & Schwarzschild 1992) shows that the 1:1 (circular) periodic orbit is naturally stable for $\alpha$ not too large.

Let us recall that this description is true provided that $\Phi$ is a smooth function of $m_{\alpha}$. If the potential has a singularity or a cusp at the origin, then the analysis may be different. Therefore the discussion given above is suitable for potentials that are not "hard'' at the origin, that is, those for which the deflection angle $\Delta\theta$ is close to $\pi$ when $p_{\theta}\rightarrow 0$ (see BT87, MS89 for details).

Let us write the full Hamiltonian in polar coordinates

$$\begin{eqnarray*}H(p_{\rm r},p_{\theta},r,\theta)={p^2_{\rm r}\over 2}+{p^2_{\theta}\over 2r^2}+ \Phi\left(m_{\alpha}(r,\theta)\right), \end{eqnarray*}$$

where $p_{\rm r}=\dot{r}$, $p_{\theta}=r^2\dot{\theta}$. Assume that $\alpha$ is small, that is, , so we can expand $\Phi\left(m_{\alpha}(r,\theta)\right)$ in powers of $\alpha$ and we can separate the part independent of $\theta$

 

$$% H(p_{\rm r},p_{\theta},r,\theta) =p^2_{\rm r}/2+p^2_{\theta}/2r^2+\phi_{\alpha}(r)- {\alpha\over 2}f_1(r)\cos 2\theta$$

$$-{\alpha^2\over 4}f_2(r)\left(\cos 2\theta-{1\over 4}\cos 4\theta\right) +\cdots,$$ (1)

where

 

$$% \begin{array}{l} \phi_{\alpha}(r)=\Phi(r^2)+{\alpha\over 2}f_1(... ...\Phi'(r^2)r^2\ge 0,\qquad f_2(r)=\Phi''(r^2)r^4\le 0.\hspace{6.5mm} \end{array}$$ (2)

From (1) and (2) the Hamiltonian can be written as

 

$$% \begin{array}{l} H(p_{\rm r},p_{\theta},r,\theta)=H_0(p_{\rm r}... ..._{\rm m}\cos 2m\theta,\ f_n(r)= \Phi^{(n)}(r^2)r^{2n}. \end{array}\hspace{-8mm}$$ (3)

H₀ is an integrable Hamiltonian being H₀=h₀ itself and $p_{\theta}=p_{\theta}^{{\rm o}}$ the unperturbed integrals and $\alpha^n V_n$ are small perturbations (see the remark at the end of this subsection). So, from now on, when we refer to unperturbed motion, we mean orbits in H₀ even though it depends on $\alpha$.

The unperturbed system is just a central field. So r oscillates between two boundaries, $r_{\rm m}(h_0,p_{\theta}^{{\rm o}})\le r^{{\rm o}}(t)\le r_{\rm M}(h_0,p_{\theta}^{{\rm o}})$, with frequency $\omega_{\rm r}$, while $\theta^{{\rm o}}$ varies on the circle $\S^1$. The frequency in the tangential direction is $\omega_{\theta}=\kappa\,\omega_{\rm r}$where $\kappa=\Delta\theta/2\pi<1$ is, in general, irrational. The time evolution of $\theta$ can be written as $\theta^{{\rm o}}(t)= \omega_{\theta}t+\Theta(t)$ where $\Theta$ is a $2\pi/\omega_{\rm r}$-periodic function of time.

Let us focus the attention, in the perturbed system, on the dynamics in the tangential direction. Keeping terms up to first order in $\alpha$ in (3) we get

$$% \dot{p}_{\theta}=-{\partial H\over\partial\theta}\approx -... ...{\partial V_1\over\partial\theta}= -\alpha f_1(r)\sin 2\theta.$$ (4)

From (4) a simple manipulation shows that the latter can be written as

$$% {{\rm d}{\cal K}\over{\rm d}t}+\alpha g_1(r(t))\sin 2\theta... ...al K}\equiv{p^2_{\theta}\over 2}-{\alpha\over 2}g\cos 2\theta,$$ (5)

where g and g₁(r(t)) are the average and oscillating parts of f₁(r(t))r²(t) respectively:

 

$$% \begin{array}{l} g=\langle f_1(r)r^2 \rangle=\langle \Phi'(r^2)r^4 \rangle>0, \\ g_1 \left(r(t)\right)=f_1\left(r(t)\right)r^2(t)-g. \end{array}$$ (6)

To keep order $\alpha$ in the perturbation, in the second term in the first of (5) we can replace the actual values of $r,\,\theta$ by their unperturbed values $r^{{\rm o}}(t)$, $\theta^{{\rm o}}(t)$. Since the unperturbed motion is completely regular, we can expand $g_1\left(r^{{\rm o}}(t)\right)\sin 2\theta^{{\rm o}}(t)$ in Fourier series, with basic frequencies $\omega_{\rm r}$ and $\omega_{\theta}$

$$% g_1\left(r^{{\rm o}}(t)\right)\sin 2\theta^{{\rm o}}(t)=\Re... ...e{g}_k{\rm e}^{i(k\omega_{\rm r}+2\omega_{\theta})t} \right\},$$ (7)

where $\tilde{g}_k$ are certain complex coefficients. Assuming quasiperiodicity (which is the more abundant behaviour if $\alpha$ is small), i.e. $\kappa$ irrational, we easily see from (7) that $\langle g_1 \sin 2\theta^{{\rm o}} \rangle\approx 0$. Hence if we average the first in (5) over several radial periods we see that (see below)

$$% {\cal K}={p^2_{\theta}\over 2}-\Omega^2\cos 2\theta, \qquad\Omega^2={\alpha\over 2}g> 0,$$ (8)

is an approximate invariant. ${\cal K}$ plays the role of the total energy in a simple pendulum model where $\Omega$ is the small oscillation frequency. Therefore two critical values of ${\cal K}$ exist: $-\Omega^2$ and $\Omega ^2$. For ${\cal K}=-\Omega^2$, $(\theta,\,p_{\theta})=(0,0)$ is a stable equilibrium point: the motion is stable along the x axis. On the other hand, for ${\cal K}=\Omega^2$ we have the separatrix and the unstable equilibrium points are $(\theta,p_{\theta})= (\pm\pi/2,0)$: the motion along the y axis is unstable. The domain of box orbits, that oscillate about the long-axis, corresponds to $\vert{\cal K}\vert<\Omega^2$ and the domain of loop orbits, that rotate about the origin, to ${\cal K}>\Omega^2$. The separatrix, $p_{\theta}^{\rm s}=\pm 2\Omega\cos\theta^{\rm s}$, separates then different kinds of motion: oscillations and rotations; i.e. box and loop orbits. For ${\cal K}\gg \Omega^2$, ${\cal K}\approx p^2_{\theta}/2$: the kinetic energy in the tangential direction. The largest value of ${\cal K}$ corresponds to the largest $p_{\theta}$, which appears for the 1:1 periodic orbit. For $V_1\ne 0$ this periodic orbit should not be circular but elliptic with small eccentricity (see below).

Since $\Omega$ is a measure of the amplitude of the torque, we conclude that the relevant parameter that defines the orbital family is the relative value of the rotational energy with respect to the strength of the torque, which in turn depends on the degree of flatness of the potential.

From the above discussion it turns out that a limit angle, $\theta_{\rm l}$, could exist

$$\begin{eqnarray*}\cos 2\theta_{\rm l}\approx -{{\cal K}\over\Omega^2}, \end{eqnarray*}$$

which is another way to conclude that $\vert{\cal K}\vert<\Omega^2$ for boxes. However it is important to remark that this bound for $\theta$ appears for r bounded away from 0. When $p_{\theta}$ is small, which is the case for boxes, the analysis of the motion in a neighbourhood ${\cal R}$ of r=0 should be done in a different way since the origin is a singular point in this description. As we assume that the potential is regular at r=0, we can approximate $\phi(x,y)$ by a harmonic oscillator in ${\cal R}$. The approximate invariants of motion are then the energy in each degree of freedom $h_x,\, h_y$. But the Lissajous-like orbits in a harmonic oscillator with incommensurable frequencies are dense in ${\cal R}$ whenever $h_x,\,h_y\ne 0$, so no bound for $\theta$ exists while the star is in ${\cal R}$.

The frequency $\Omega$, depends on the average $\langle f_1(r) r^2\rangle$ over the unperturbed motion (Eqs. (6) and (8)). As defined in (2), f₁(r) can be put in terms of the circular speed, $v_{\rm c}$: $f_1(r)=v_{\rm c}^2(r)/2$. So from (8) follows that

$$\begin{eqnarray*}\Omega^2={\alpha\over 4}\langle v_{\rm c}^2(r)r^2 \rangle. \end{eqnarray*}$$

For the realistic case of flat rotation curves at large radii we get

$$% \Omega^2\approx{\alpha\over 4} v^2_{\rm c}\langle r^2 \rangle\sim {\alpha\over 12} v_{\rm c}^2r_{\rm M}^2(1+\beta+\beta^2),$$ (9)

where $0\le \beta =r_{\rm m}/r_{\rm M}\le1,$ and $r_{\rm M}>r_{\rm m}$ are the two roots of the equation (see Eq. (3))

$$\begin{eqnarray*}(p_{\theta}^{{\rm o}})^2-2r^2\left(h_0-\Phi(r^2)\right)+ {1\over 2}\alpha\Phi'(r^2)r^4=0, \end{eqnarray*}$$

which, for $r_{\rm M}$, can be approximated by

$$% (p_{\theta}^{{\rm o}})^2-2r_{\rm M}^2\left(h_0-\Phi(r_{\rm M}^2)\right)+ {1\over 2}\alpha v_{\rm c}^2r_{\rm M}^2\approx 0.$$ (10)

For the estimate in (9), where a factor 2 should be added if $\beta=0$, we approximate the time-average of r² by the r-average over the allowed interval. This is not true in general but it provides a rough estimate of the average that will help us later.

The invariant ${\cal K}$ is in some sense local, since unperturbed orbits with different angular momentum will have different values of the frequency: $\Omega= \Omega(h_0,p^{{\rm o}}_{\theta})$. From (9) and (10) it is not difficult to conclude that the largest $\Omega$ is expected for minimum $\vert p_{\theta}^{{\rm o}}\vert$; $p_{\theta}^{{\rm o}}=0$, i.e., for radial orbits, while the smallest one for maximum $\vert p_{\theta}^{{\rm o}}\vert$, i.e., for circular orbits (see also Sect. 2.2).

For the case of the 1:1 periodic orbit we can write,

$$% \Omega_{1:1}^2\approx{\alpha\over 4}v_{\rm c}^2a^2,$$ (11)

where a is the circular radius defined by

$$% \Phi(a^2)\approx h_0-\left(1+{\alpha\over 2}\right){v_{\rm c}^2\over 2}\cdot$$ (12)

Then, the maximum value of ${\cal K}$ lies somewhere between

 

The whole picture given above is true for small $\alpha$ (q close to 1). Indeed, this approach makes sense when the x-axis periodic orbit is stable. It is well known, that for large values of $\alpha$ ( $\alpha\sim 1$, $q\sim 0.7$) the latter orbit could become unstable bifurcating to another periodic orbit. A sub-family associated to this new orbit appears. It is expected also that the x-axis periodic orbit lies now in a narrow stochastic layer around the separatrix of the resonance (see below). Other high-order resonances would occupy some region of the phase space and many zones of stochastic motion would also appear. So it is hard to speak then only about box or loop orbits when the perturbation is large (in fact, the term boxlets is often used in this case, see for example, MS89). We refer to Sect. 5, where a global study of the logarithmic potential reveals that this "very chaotic'' panorama does not show up even for large values of h and $\alpha$. Nevertheless, in general, bounds to the value of q would appear: $0<q_0\le q\le 1$. This bound comes from the Poisson equation, $\nabla^2\phi=4\pi G\rho$ with $\rho>0$.

One should remark that ${\cal K}$ given by (8) was obtained neglecting high order terms, assuming quasiperiodicity in the unperturbed motion and averaging to zero the oscillating part. Thus the pendulum model is a rough first approximation to the dynamics and other perturbing terms should be present. However, the main effect of perturbations to the pendulum is to distort somehow the invariant curves and to give rise to a stochastic layer around the separatrix. That is, box and loop should be actually separated by a stochastic layer instead of a separatrix. The larger the strength of the perturbation, the larger is the width of the layer (for details about the dynamics of the pendulum model and perturbation effects, see Chirikov 1979).

The derivation given above for ${\cal K}$ is a justification of the invariant introduced ad-hoc to compute certain models of elliptical galaxies that respect a third integral. Indeed, if the potential has the form

$$% \phi(\vec{r})=\psi(r)+{\chi(\theta)\over r^2},$$ (14)

with $\psi$ and $\chi$ arbitrary functions, then a third integral exists

$$\begin{eqnarray*}I_3={p_{\theta}^2\over 2}-\chi(\theta). \end{eqnarray*}$$

The form (14) is a particular case of a more general type of potentials introduced almost one century ago by Eddington (1915) to study oblate distributions where the ellipsoidal velocity law is exactly satisfied (he showed, however, that the latter condition does not hold if in (14) $\chi\ne 0$). Later on, this model was adopted, for instance, by Lynden-Bell (1967) in his investigations on statistical mechanics of violent relaxation in rotating elliptical systems. As was pointed out by Eddington, Lynden-Bell and others, (14) is unsuitable for any galactic potential so, in general, the third integral for a more realistic model is supposed to be

$$\begin{eqnarray*}I_3={p_{\theta}^2\over 2}-\xi(r,\theta), \end{eqnarray*}$$

where $\xi(r,\theta)$ is such that I₃ should satisfy approximately the collisionless Boltzmann equation (see Petrou 1983a,b for more details). No other explicit integral can be expected for a general potential of the form $\phi(r,\theta)=\psi(r)+Q(r)\chi(\theta)$.

Note, however, that for a bar-like galaxy the multipolar expansion of $\phi(\vec{r})$ has as dominant terms

$$\begin{eqnarray*}\phi(\vec{r})\approx\psi(r)+Q(r)P_2(\cos\theta)\equiv\psi_1(r)+Q_1(r)\cos 2\theta, \end{eqnarray*}$$

where $P_2(\mu)$ is the Legendre polynomial of degree 2. If $\psi$ and Qare regular at r=0 we recover the Hamiltonian (3) taken as a model for the above discussion.

  
2.2 Example using the 2D logarithmic potential

The logarithmic potential has the form:

$$% \phi(x,y)={1\over 2}p_0^2\ln(x^2+y^2/q^2+r_{\rm c}^2)$$ (15)

where p²₀ is a constant and $r_{\rm c}$, the core radius, plays the role of a softening parameter to avoid the singularity at the origin. Hence $\phi$ is regular at r=0 and, in ${\cal R}$ defined as $r\ll r_{\rm c}$, $\phi(x,y)\sim x^2+y^2/q^2$. The circular speed is $v_{\rm c}^2\sim p_0^2$ for $r\gg r_{\rm c}$ and in what follows we set p₀²=2, $r_{\rm c}=0.1$ (see, however, next section). The semiaxis ratio, q, can take any value within the range $q_0\le q\le 1$ where $q_0^2=1/2-r_{\rm c}^2/\exp(2h/p_0^2)$, being h the energy level. Note that q is the semiaxis ratio of isopotential curves. The semiaxis ratio for isodensity curves, $\tilde{q}$, is related with q by $\tilde{q}\approx q\sqrt{2q^2-1}$ for $r\gg r_{\rm c}$ (see BT87), so isodensities are always much flatter than isopotentials. Typical orbits and their corresponding surfaces of section for this potential may by found elsewhere, for instance, BT87, MS89, PL96. In particular, the level h=-0.4059 is adopted hereafter following PL96 (see, however, Sect. 5). For the latter value of h, $q_0\approx 0.696$ and the frequency of the long-axis and 1:1 (circular if q=1) periodic orbits are about 2 and 3 respectively (see PL96).

As we showed in the previous section, one of the main differences between loop and box orbits is that the latter pass through the origin (projected on the configuration space, the (x,y)-plane) while the former do not. Since for boxes $\vert{\cal K}\vert\le\Omega^2$, $p_{\theta}=0$ is a suitable choice for the initial angular momentum for that orbits, while ${\cal K}>\Omega^2\rightarrow \vert p_{\theta}\vert>2\Omega$ is enough to ensure a loop orbit (see Eq. (8) and discussion below). In Cartesian variables, since $p_{\theta_0}=x_0p_{y_0}- y_0p_{x_0}$, in any case we can set, for instance, y₀=0 and then p_x0=0for loops and x₀=0 for boxes. Therefore, for the energy level h, x₀ and p_x0 parameterize, respectively, loop and box orbits. The same applies if we exchange $x_0,p_{x_0}\leftrightarrow y_0,p_{y_0}$. This choice of a single Cartesian variable to label orbits in each family will help us later (Sect. 4.2), but to investigate the connection between the actual motion and the pendulum model, we will proceed with polar coordinates.

Figure 1: Surfaces of section $(\theta ,\,p_{\theta })$, $r_{\rm s}=0.5$, $p_{\rm r}>0$ for $\theta$ restricted to the interval $(-\pi /2,\pi /2)$ and for several orbits in the logarithmic potential (15) with $r_{\rm c}=0.1$, p02=2, h=-0.4059 and q=0.9, 0.8 and 0.7 (see text for details)

A suitable map that should reveal the dynamics in the tangential direction is a $(\theta,p_{\theta})$ section for constant r and, for instance, with $p_{\rm r}>0$. Nevertheless one should ensure that the trajectory actually intersects this surface. We expect that box orbits always cross any surface $r=r_{\rm s}$ provided that $r_{\rm m}<r_{\rm s}<r_{\rm M}$. However, this is not the case for all loop orbits. For instance, in the limiting case when $q\rightarrow 1$, the 1:1 periodic orbit becomes circular and it is fully contained in r=a. Another restriction is that $r_{\rm s}>r_{\rm c}$ because, as mentioned, for $r\ll r_{\rm c}$, ${\cal K}$ is not well defined.

In Fig. 1 we show $(\theta,p_{\theta})$ surfaces of section for $r_{\rm s}=0.5$ and three different values of q, 0.9, 0.8 and 0.7. Note that for these values of q, the corresponding values of $\alpha$ are not small, 0.234, 0.562 and 1.04 respectively. Nevertheless the global motion resembles a pendulum model. The origin corresponds to the long-axis periodic orbit while the unstable points $(\pm\pi/2,0)$ to the short-axis periodic orbit. The outermost curves correspond to loop orbits relatively close to the 1:1 periodic orbit (for both senses of rotation). Indeed, from (12) one immediately finds that for this orbit $a\approx 0.45$, 0.42, 0.38 for q=0.9, 0.8, 0.7 respectively. Because $V_1\ne 0$, a gives an estimation of the mean value of the semiaxis of this elliptic orbit (see Sect. 4.2 for a more accurate determination).

From Fig. 1 we see that box and loop families are separated by a relatively narrow stochastic layer and, for large values of the angular momentum, the latter is nearly an integral of motion. As q decreases from 0.9 the departures of the actual motion from the pendulum model become significant. For instance, for q=0.7 we observe that the x-axis periodic orbit becomes unstable and a bifurcation to a 2-periodic orbit appears (see footnote at the end of next section). This sub-family of box family is conformed by the well-known banana orbits (see, for instance, MS89). Note the thin chaotic layer that separates bananas and boxes. We distinguish several resonances and a comparatively large stochastic layer separates boxes and loops. Some small islands are also present in the loop domain. Nevertheless, only a few loop orbits appear in this figure (q=0.7) because almost all the intersections with the surface $r_{\rm s}=0.5$ correspond to box orbits. Anyway the domain of box family increases against that of loops as the flatness of the potential increases.

Figure 2: $\Omega ^2$ vs. $p_{\theta _0}$ for several orbits with r0=0.5, $\theta _0=0$ and $r_{\rm c}=0.1$, h=-0.4059. The small oscillation frequency was computed as $\Omega^2=0.5\alpha \langle\Phi'(r^2)r^4\rangle=0.5\alpha \langle r^4/(r^2+r_{\rm c}^2)\rangle$ over $t\approx 500$ periods and using the actual r values (see text)

In Fig. 2 we show the computed values of $\Omega$ (actually $\Omega ^2$) for several (about 1500) orbits with initial conditions along the $p_{\theta}$axis. We average $\Phi'(r^2)r^4=r^4/(r^2+r_{\rm c}^2)$ over 1500 units of time which is $\sim 500$ periods of the x-axis periodic orbit, using the values of r obtained by solving the equations of motion corresponding to (15) with q=0.9, 0.8, 0.7. We checked that the oscillating part, g₁(r(t)) - see Eq. (6) - averages to zero. For q=0.9 and 0.8 almost all orbits are either box or loop while for q=0.7, those for belong to the banana sub-family (see Fig. 1). In any case $\Omega ^2$ does not change too much over the whole range, less than by a factor 2. This variation of $\Omega$ can be computed using (9), (10), (11) and (12). It is immediate that the change in $\Omega ^2$ from radial ( $p_{\theta}=0$) to circular (maximum $p_{\theta}$) orbits is

$$\begin{eqnarray*}{\Delta\Omega^2\over\Omega_{1:1}^2}\sim{2\over 3}{r_{\rm M}^2(0... ...0)+r_{\rm c}^2\over a^2+r_{\rm c}^2}\approx {{\rm e}^{p_0^2/2}}, \end{eqnarray*}$$

which for $r_{\rm M}(0), a\gg r_{\rm c}$ and p₀²=2 leads to $r_{\rm M}^2/a^2\sim{\rm e}$ and then $\Delta\Omega^2/\Omega^2_{1:1}\sim 0.8$.

Figure 3: ${\cal K}/\Omega ^2$ vs. $p_{\theta _0}$ for the same ensemble of orbits than in Fig. 2. The line ${\cal K}/\Omega ^2=1$ separates the domains of box and loop orbits. The theoretical invariant for the whole ensemble, defined by $\theta _0=0$, was computed as ${\cal K}_{\rm T}=p_{\theta _0}^2/2-\Omega ^2$. The numerical invariant was averaged, for each orbit, over $t\approx 500$ periods (see text)

From this figure, for q=0.9 and 0.8 we see a similar, smooth behaviour of $\Omega$, except in those narrow intervals where a small island seems to be present. Since the width of the stochastic layer for these two values of q is rather small (particularly due to the fact that we are crossing the layer through the thinest part), almost no evidence of its existence could be inferred from the figure. However for q=0.7, $\Omega$ looks noisy in the vicinity of the long-axis periodic orbit as well as on the separatrix, revealing the existence of chaotic zones. Much more discontinuities can be also observed, certainly due to the existence of many resonances. For the banana family, as expected, the behaviour is different. In the latter case $\Omega$ is almost constant in all this range. This can be understood looking back to Fig. 1 and recalling that we are taking initial conditions along the $p_{\theta}$ axis; in all cases the banana orbits are close to the marginal one.

Having the small oscillation frequency we compute ${\cal K}$. First we took a set of representative orbits and we studied the time evolution of ${\cal K}$. We observed a nearly constant value of ${\cal K}$ but, for box orbits, a periodical peak structure was observed. This was not the case for loop orbits. These peaks are due to the fact that ${\cal K}$ is not an invariant when $r<r_{\rm c}$; all of them appeared at minimum r. To reduce their effect on the instantaneous value of the invariant, we averaged ${\cal K}$ over the whole time interval. This procedure, however, leads to slightly larger values of ${\cal K}$ due to the cumulative effect.

In Fig. 3 we show the computed values of ${\cal K}/\Omega ^2$ for the same ensemble of orbits than in Fig. 2. For the long-axis periodic orbit ${\cal K}=-\Omega^2$ (a zoom around $p_{\theta_0}=0$ reveals that the initial value of ${\cal K}/\Omega^2=-1$ also for q=0.7). The line ${\cal K}=\Omega^2$ separates box and loop domains and a simple comparison with Fig. 1 shows that the separatrix (to be precise, the stochastic layer) is very close to that line. The maximum value observed for ${\cal K}$ is, in all cases, within the interval derived in (13).

Figure 4: Theoretical surfaces of section $(\theta ,\,p_{\theta })$ for q=0.9, 0.8 and 0.7 computed using ${\cal K}_{\rm T}$ and $\Omega$ obtained numerically. Compare with Fig. 1

Figure 3 is illustrative to show which family dominates the dynamics at given level of energy and to identify other sub-families of the principal ones. For instance, for q=0.9 box and loop orbits populate the phase space and no other sub-family seems to be significant. On the other hand for q=0.7 almost all the phase space is filled by box orbits but the banana family occupies a fraction of the phase space (see Sect. 4.2). Loops are scanty and some chaotic domains seems to become important. Recall that the actual volume of the phase space occupied by any of these families would not be measured by the size of the corresponding ${\cal K}$ intervals. We also include in Fig. 3 the theoretical value of ${\cal K}$ which, for the ensemble considered, is ${\cal K}_{\rm T}=p^2_{\theta_0}/2-\Omega^2$, where we take for $\Omega$ the computed values given in Fig. 2. Note the good agreement between ${\cal K}$ and ${\cal K}_{\rm T}$ for loops and boxes but, of course, they differ for the banana sub-family.

In Fig. 4 we show the surfaces of section, $(\theta,p_{\theta})$, computed using ${\cal K}_{\rm T}$ and $\Omega$ obtained numerically. A comparison with Fig. 1 reveals that the actual motion is well approximated by a pendulum model for q=0.9 and for q=0.8, but the agreement is not so good, as expected, for q=0.7.