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.