We show that the Lyapunov exponents of a periodic orbit can be easily obtained from the eigenvalues of the monodromy matrix. It turns out that the Lyapunov exponents of simply stable periodic orbits are all zero, simply unstable periodic orbits have only one positive Lyapunov exponent, doubly unstable periodic orbits have two different positive Lyapunov exponents and the two positive Lyapunov exponents of complex unstable periodic orbits are equal. We present a numerical example for periodic orbits in a realistic galactic potential. Moreover, the center manifold theorem allowed us to show that stable, simply unstable and doubly unstable periodic orbits are the mothers of families of, respectively, regular, partially and fully chaotic orbits in their neighbourhood.

We have recently presented strong evidence that chaotic orbits that obey one isolating integral besides energy exist in a toy Hamiltonian model with three degrees of freedom and are bounded by regular orbits that isolate them from the Arnold web. The interval covered by those numerical experiments was equivalent to about one million Hubble times in a galactic context. Here we use a four dimensional map to confirm our previous results and to extend that interval fifty times. We show that, at least within that interval, features found in lower dimension Hamiltonian systems and maps are also present in our study, e.g., within the phase space occupied by a chaotic orbit that obeys one integral there are subspaces where that orbit does not enter and are, instead, occupied by regular orbits that, if tori, bound other chaotic orbits obeying one integral and, if cantori, produce stickiness. We argue that the validity of our results might exceed the time intervals covered by the numerical experiments.

Three types of orbits are theoretically possible in autonomous Hamiltonian systems with three degrees of freedom: fully chaotic (they only obey the energy integral), partially chaotic (they obey an additional isolating integral besides energy) and regular (they obey two isolating integrals besides energy). The existence of partially chaotic orbits has been denied by several authors, however, arguing either that there is a sudden transition from regularity to full chaoticity, or that a long enough follow up of a supposedly partially chaotic orbit would reveal a fully chaotic nature. This situation needs clarification, because partially chaotic orbits might play a significant role in the process of chaotic diffusion. Here we use numerically computed Lyapunov exponents to explore the phase space of a perturbed three dimensional cubic force toy model, and a generalization of the Poincar\'e maps to show that partially chaotic orbits are actually present in that model. They turn out to be double orbits joined by a bifurcation zone, which is the most likely source of their chaos, and they are encapsulated in regions of phase space bounded by regular orbits similar to each one of the components of the double orbit.

We built two self-consistent models of triaxial, cuspy, rotating stellar systems adding rotation to non-rotating models presented in previous papers of this series. The final angular velocity of the material is not constant and varies with the distance to the center and with the height over the equator of the systems, but the figure rotation is very uniform in both cases. Even though the addition of rotation to the models modifies their original semiaxes ratios, the final rotating models are considerably flattened and triaxial. An analysis of the orbital content of the models shows that about two thirds of their orbits are chaotic yet the models are very stable over intervals of the order of one Hubble time. The bulk of regular orbits are short axis tubes, while long axis tubes are replaced by tubes whose axes lie on the short-long axes plane, but do not coincide with the major axis. Other types of regular orbits that do not appear in non-rotating systems, like horseshoes and orbits that cross themselves, are also found in the present models. Finally, our frequency maps show empty regions where studies of orbits on fixed potentials found orbits, a likely consequence of the self-consistency of our models that excludes them.

In several previous investigations we presented models of triaxial stellar systems, both cuspy and non cuspy, that were highly stable and harboured large fractions of chaotic orbits. All our models had been obtained through cold collapses of initially spherical $N$--body systems, a method that necessarily results in models with strongly radial velocity distributions. Here we investigate a different method that was reported to yield cuspy triaxial models with virtually no chaos. We show that such result was probably due to the use of an inadequate chaos detection technique and that, in fact, models with significant fractions of chaotic orbits result also from that method. Besides, starting with one of the models from the first paper in this series, we obtained three different models by rendering its velocity distribution much less radially biased (i.e., more isotropic) and by modifying its axial ratios through adiabatic compression. All three models yielded much higher fractions of regular orbits than most of those from our previous work. We conclude that it is possible to obtain stable cuspy triaxial models of stellar systems whose velocity distribution is more isotropic than that of the models obtained from cold collapses. Those models still harbour large fractions of chaotic orbits and, although it is difficult to compare the results from different models, we can tentatively conclude that chaoticity is reduced by velocity isotropy.

In the first paper of this series we used the N--body method to build a dozen cuspy (gamma ~ 1) triaxial models of stellar systems, and we showed that they were highly stable over time intervals of the order of a Hubble time, even though they had very large fractions of chaotic orbits (more than 85 per cent in some cases). The models were grouped in four sets, each one comprising models morphologically resembling E2, E3, E4 and E5 galaxies, respectively. The three models within each set, although different, had the same global properties and were statistically equivalent. In the present paper we use frequency analysis to classify the regular orbits of those models. The bulk of those orbits are short axis tubes (SATs), with a significant fraction of long axis tubes (LATs) in the E2 models that decreases in the E3 and E4 models to become negligibly small in the E5 models. Most of the LATs in the E2 and E3 models are outer LATs, but the situation reverses in the E4 and E5 models where the few LATs are mainly inner LATs. As could be expected for cuspy models, most of the boxes are resonant orbits, i.e., boxlets. Nevertheless, only the (x, y) fishes of models E3 and E4 amount to about 10 per cent of the regular orbits, with most of the fractions of the other boxlets being of the order of 1 per cent or less.

We used the N-body code of Hernquist and Ostriker (1992) to build a dozen cuspy ({\gamma}\approx 1) triaxial models of stellar systems through dissipationless collapses of initially spherical distributions of 10^6 particles. We chose four sets of initial conditions that resulted in models morphologically resembling E2, E3, E4 and E5 galaxies, respectively. Within each set, three different seed numbers were selected for the random number generator used to create the initial conditions, so that the three models of each set are statistically equivalent. We checked the stability of our models using the values of their central densities and of their moments of inertia, which turned out to be very constant indeed. The changes of those values were all less than 3 per cent over one Hubble time and, moreover, we show that the most likely cause of those changes are relaxation effects in the numerical code. We computed the six Lyapunov exponents of nearly 5,000 orbits in each model in order to recognize regular, partially and fully chaotic orbits. All the models turned out to be highly chaotic, with less than 25 per cent of their orbits being regular. We conclude that it is quite possible to obtain cuspy triaxial stellar models that contain large fractions of chaotic orbits and are highly stable. The difficulty to build such models with the method of Schwarzschild (1979) should be attributed to the method itself and not to physical causes.

Here I present a review of the work done on the presence and effects of chaos in elliptical galaxies plus some recent results we obtained on this subject. The fact that important fractions of the orbits that arise in potentials adequate to represent elliptical galaxies are chaotic is nowadays undeniable. Alternatively, it has been difficult to build selfconsistent models of elliptical galaxies that include significant fractions of chaotic orbits and, at the same time, are stable. That is specially true for cuspy models of elliptical galaxies which seem to best represent real galaxies. I argue here that there is no physical impediment to build such models and that the difficulty lies in the method of Schwarzschild, widely used to obtain such models. Actually, I show that there is no problem in obtaining selfconsistent models of elliptical galaxies, even cuspy ones, that contain very high fractions of chaotic orbits and are, nevertheless, highly stable over time intervals of the order of a Hubble time.

Orbits in the principal planes of triaxial potentials are known to be prone to unstable motion normal to those planes, so that three dimensional investigations of those orbits are needed even though they are two dimensional. We present here an investigation of such orbits in the well known logarithmic potential which shows that the third dimension must be taken into account when studying them and that the instability worsens for lower values of the forces normal to the plane. Partially chaotic orbits are present around resonances, but also in other regions. The action normal to the plane seems to be related to the isolating integral that distinguishes regular form partially chaotic orbits, but not to the integral that distinguishes partially from fully chaotic orbits.

Numerical experiments conducted by Fellhauer et al. (MNRAS, 372, 338, 2006) suggest that a supercluster may capture up to about 40 per cent of its mass from the galaxy where it belongs. Nevertheless, in those experiments the cluster was created making appear its mass out of nothing, rather than from mass already present in the galaxy. Here we use a thought experiment, plus a few simple computations, to show that the difference between the dynamical effects of these two scenarios (i.e., mass creation vs. mass concentration) is actually very important. We also present the results of new numerical experiments, simulating the formation of the cluster through mass concentration, that show that trapping depends critically on the process of cluster formation and that the amounts of gained mass are substantially smaller than those obtained from mass creation.

Here I review recent work, by other authors and by myself, on some particular topics related to the regular and chaotic motion in elliptical galaxies. I show that it is quite possible to build highly stable triaxial stellar systems that include large fractions of chaotic orbits and that partially and fully chaotic orbits fill different regions of space, so that it is important not to group them together under the single denomination of chaotic orbits. Partially chaotic orbits should not be confused with weakly fully chaotic orbits either, and their spatial distributions are also different. Slow figure rotation (i.e., rotation in systems with zero angular momentum) seems to be always present in highly flattened models that result from cold collapses, with the rotational velocity diminishing or becoming negligibly small for less flattened models. Finally, I comment on the usefulness and limitation of the classification of regular orbits via frequency analysis.

As a first step towards a comprehensive investigation of stellar motions within globular clusters, we present here the results of a study of stellar orbits in a mildly triaxial globular cluster that follows a circular orbit inside a galaxy. The stellar orbits were classified using the frequency analysis code of Carpintero and Aguilar and, as a check, the Liapunov characteristic exponents were also computed in some cases. The orbit families were obtained using different start spaces. Chaotic orbits turn out to be very common and while, as could be expected, they are particularly abundant in the outer parts of the cluster, they are still significant in the innermost regions. Their relevance for the structure of the cluster is discussed.