Abstract
Two self-consistent (N-body) non-rotating equilibrium models of elliptical galaxies with smooth central density profiles (called ‘Q’ and ‘C’ models) are constructed, starting from quiet and clumpy cosmological initial conditions, respectively. Both models are triaxial. The Q model has an E7 maximum ellipticity in the inner parts and tends to E6 or E5 maximum ellipticity in the outer parts. The C model has a maximum ellipticity E4 in the inner parts and tends to an E2 or E1 in the outer parts.
For each model, we identify the particles moving in chaotic orbits with the Lyapunov number exceeding a particular threshold (namely, 10−2.8, in units of the inverse radial periods of the particular orbits). At energy levels in the deepest 30 per cent of the potential well, no chaotic orbits were detected in the above limit of chaoticity. In the Q model, the detected chaotic part is 32 per cent of the total mass. This part has a nearly spherical distribution. It imposes limitations on the maximum ellipticity of the system, in spite of the fact that only a part of less than about 8 per cent of the total mass moves in chaotic orbits and is able to develop chaotic diffusion within a Hubble time. In the C model, the detected chaotic part is about 26 per cent of the total mass, but only less than 2 per cent can develop chaotic diffusion within a Hubble time.
These chaotic components produce surface density profiles flatter than the profiles of the rest of the mass, particularly in the Q model. The two profiles intersect at a given distance, where the overall profile forms an observable hump, especially if the surface density profiles are taken along the shortest axis of the projection.
1 Introduction
It is generally believed that chaos in triaxial galactic potentials with smooth central regions is limited. Smooth central regions can be represented by the central regions of integrable models such as the Staekel potentials (e.g. Binney & Tremaine 1987; Contopoulos 1994). A particular model with a smooth central profile is the perfect ellipsoid, which is a special case of the Staekel model family (Kuzmin 1953, 1956; de Zeeuw 1985; de Zeeuw & Lynden-Bell 1985). In such models, the regular orbits can be classified into four main types, i.e. box orbits, inner long axis tubes, outer long axis tubes and short axis tubes (de Zeeuw 1985; Statler 1987). For simplicity, we call the inner long axis tubes ‘box-like orbits’ and the short axis tubes ‘1 : 1 tube orbits’. In weakly non-integrable models, the chaotic regions of phase space are located in between the regions defined by the various types of regular orbits. The triaxial shape is mainly supported by the box or the box-like regular orbits.
Chaotic motion becomes more pronounced under the presence of a central black hole or of a central cusp in the density profile. Both cases have the effect of destroying the regular character of many box or box-like orbits which pass close to the centre (Gerhard & Binney 1985; Merritt & Fridman 1996).
N-body simulations of gravitational dissipationless collapses give relaxed systems that are very similar to elliptical galaxies. The resulting self-consistent potential is non-integrable, although it is in general close to a potential of the Staekel type. There are, however, non-negligible deviations from an integrable potential which produce chaotic motion. While most orbits are not expected to be strongly chaotic, the overall influence of chaotic orbits in the distribution function is important.
We are thus led to consider the question of the relative importance of regular versus chaotic orbits in self-consistent triaxial N-body systems. In particular, we address the following questions: (i) find estimates of the number of particles in chaotic orbits versus the number of particles in regular orbits and (ii) find the regions of phase space occupied by particles of one or the other type of motion. Notice that self-consistency imposes restrictions to the amount of the material allowed to occupy the various regions of phase space, regular or chaotic (Contopoulos, Efthymiopoulos & Voglis 2000; Contopoulos, Voglis & Kalapotharakos 2002).
Two more questions are of interest. First, the dependence of the relative amount of particles in chaotic orbits on the particular shape of a system. Secondly, the identification of possible observable tracers of the chaotic component. In particular, we calculate several observable quantities that depend on the amount of mass in chaotic orbits.
In the present paper we investigate the above questions in the case of non-cuspy self-consistent models. For this purpose, we construct two self-consistent N-body models in equilibrium (N≃ 105) representing elliptical galaxies. Both models are produced from cosmological initial conditions consistent with a power-law spectrum P(k) ∝kn of the density perturbation field.
The two models have many similarities, but also considerable differences due to the fact that one model is produced by quiet initial conditions and the other by clumpy initial conditions.
In the case of quiet initial conditions, the particles start from a high ratio of radial to tangential velocity dispersion. This causes the onset of a radial orbit instability during the collapse phase. The final system is triaxial and it is characterized by the presence of a strong, non-rotating bar. This sort of instability is very well known in stellar dynamics (e.g. Polyachenko & Shukhman 1981; Fridman & Polyachenko 1984; Palmer 1994).
On the other hand, clumpy initial conditions produce hot systems in which the ratio of radial to tangential velocity dispersion is not so high. The resulting systems are still triaxial, but closer to spherical systems.
These results, either in the case of quiet or in the case of clumpy initial conditions have been verified by N-body simulations first by Carpintero & Muzzio (1995); but also by Efthymiopoulos & Voglis (2001) and Contopoulos et al. (2002).
After the initial transient period of relaxation the systems are lead to a quasi-equilibrium state in which the potential is smooth in space, and has small time-fluctuations. Considering the limit of a time-independent smooth potential, the orbits of stars are given, to a very good approximation, by the equations of motion of an autonomous Hamiltonian function (Contopoulos et al. 2000).
The construction of these models and the discussion of their basic features is made in Section 2. A description of the process and of the methods used to distinguish regular from chaotic motions is given in Section 3. The main results are presented in Section 4. Our main conclusions are summarized in Section 5.
2 The Self-Consistent Equilibrium Models
Cosmological N-body initial conditions that produce self-consistent galactic equilibrium models are obtained by a two-step procedure.
At the first step, a spherical part of the Universe, corresponding to the galactic mass scale, is considered at an epoch of high redshift. This spherical part is resolved into a number of about N1≈ 5 × 103 particles, which are given initially Hubble flow motions consistent with the Einstein—de Sitter model.
At the second step, a density perturbation is imposed to the initial system of particles. Two different models, Q and C, are obtained by two different ways of perturbing the initial system.
The C-model is evolved analytically to the same redshift by using the Zeldovich approximation (Zeldovich 1970).
The subsequent evolution is followed by a tree code (Hernquist 1987). When the systems are well relaxed to their equilibrium states (after about 10 dynamical times) every particle is split into 27 smaller particles, in a way that the coarse grained distribution function of each system remains unaffected. This procedure is described in detail in Voglis (1994). The new data of the positions and velocities of the N2= 27 ×N1≈ 1.5 × 105 particles are used as initial conditions for a further self-consistent evolution of the system by a conservative-technique code (Allen, Palmer & Papaloizou 1990). The code provides the self-consistent potential as an expansion in spherical Bessel functions and spherical harmonics up to the quadrupole terms. The coefficients of this expansion are evaluated at every time step.
The equilibrium configurations are triaxial in both models (Q and C). We call the direction of the shortest axis x and the direction of the longest axis z.
In Figs 1(a) and (b) we give the projections, on the x–z plane, of the Q and the C model, respectively, with several equidensity curves. The scaling unit of length is defined for each system to be equal to its half-mass radius. The ellipticity on this plane varies with the radius R (along the major axis). The Q model has an almost constant maximum ellipticity of about E7 in the inner parts, that decreases to about E6 or E5 in the outer parts. The C model starts with a maximum ellipticity of about E4 near the centre which decreases gradually outwards to a value between E2 and E1.
Projection on the x–z plane (shortest–longest axes) of the Q-model (panel a) and the C-model (panel b), with several equidensity contours.
Projection on the x–z plane (shortest–longest axes) of the Q-model (panel a) and the C-model (panel b), with several equidensity contours.
3 Distinguishing Between Ordered And Chaotic Orbits: The Methods
The aim is to estimate the level of chaos in each system and investigate possible implications of chaos to observable quantities. It is obvious that particularly weak chaotic orbits, beyond the fact that their identification is a very demanding task, cannot play a role in the system much different than the role of regular orbits at least as regards observable quantities.
For this reason we restrict our investigation in identifying the most chaotic orbits by defining a particular threshold in the values of their Lyapunov Characteristic Number (LCN), as described bellow.
When dealing with chaos in a galactic system, one unavoidably faces the problem of very different orbital periods of particular orbits that coexist in the same system. For example, the particles moving in the centre may have periods even two or more orders of magnitude smaller than the periods of the particles in outer orbits. Then, if one wishes to calculate Lyapunov times, one has to decide whether these times should be expressed in terms of the particular periods of individual orbits, or of the half mass crossing time which measures the average dynamical period in the system. In the following paragraphs, we describe a methodology which provides a compromise for this problem, and gives a not very much biased estimator of the chaoticity of the orbits, as the latter reflects to the macroscopic properties of the system within a Hubble time.
3.1 The method of the Specific Finite Time Lyapunov Characteristic Number
If we use a Finite Time Lyapunov Characteristic Number (FT-LCN) and integrate the orbits for times comparable to a Hubble time, then we expect to be able to identify only those chaotic orbits with LCNs greater than a limiting value (say about 0.1 in units of the inverse dynamical time of the system). However, such orbits may not exist at all, or there may be only a few, as chaos is expected to be weak. Our aim is to detect chaotic orbits with much smaller LCNs. Although the chaotic diffusion of such orbits may not be dynamically important in a Hubble time, they can still have remarkable implications for the distribution function in phase space, as shown below.
Looking for chaotic orbits with LCNs smaller than 0.1 implies that the integration time must be long enough, i.e. a multiple of the Lyapunov time TLCN= LCN−1 for every orbit and that it must scale proportionally to the individual dynamical times of the orbits.
The quantity Trj is the average radial period (time to go from pericentre to apocentre and back to pericentre) of the orbit of the particle j. This period scales with the energy Ej of the orbit as Trj≈∣Ej∣−3/2 (almost independently of the angular momentum, Voglis 1994).
The above definition of the S-FT-LCN allows us to extend our calculations to different times for every orbit, up to the point when the integration time is long enough to detect LCNs with small values reaching the desirable limit.
We consider the Hubble time to be equal to 100 half-mass crossing times (Thmct) of the system. Let the values of Lj of all the orbits are calculated for the same number of radial periods, e.g. tj/Trj= 1200. The shortest Trj is larger than 0.34 Thmct. Thus, the above number of radial periods, ensures that the particles with the shortest radial periods are integrated for more than 4 Hubble times. This number increases with the decreasing absolute value of the binding energy, reaching to about 300 Hubble times in the outer parts of the system.
This limit of integration time can detect chaotic orbits with Ljs stabilized at values no smaller than 10−3. Of course, the system contains chaotic orbits with much smaller values of Lj, as it is, for example, the case of chaotic orbits that are trapped between invariant tori. Such orbits need exponentially long times to escape by the Arnold diffusion mechanism in larger chaotic areas. Their chaotic behaviour can not be efficient in a Hubble time. For the purposes of this paper, i.e. to reveal the role of chaos, it is better such orbits to be classified as regular, because their features are closer to regular orbits than to chaotic orbits.
We adopt this integration time, tj/Trj= 1200 for each orbit, and we consider the orbits for which Lj has reached a constant value within this time to be classified as the most chaotic orbits of the system. This threshold is rather arbitrary but as we will see below it is satisfactory for our purposes.
The values of Lj can give a good answer as regards the question of the character of an orbit (chaotic, up to the defined level, or not). They can also compare the chaotic characters of two orbits, provided that their Ljs have already been stabilized in the time available. However, as they are measured in different units for every orbit, they cannot be comparable as regards the question of the efficiency of their chaoticity within a Hubble time.
For this reason we convert Ljs in common units by defining the Common Unit Finite Time Lyapunov Characteristic Number (CU-FT-LCN), or simply Lcu, as Lcu=LjThmct/Trj, in units of 1/Thmct. This number compares the chaotic orbits as regards their efficiency within a Hubble time.
3.2 The method of Alignment Index
A second method for distinguishing between ordered and chaotic orbits, that may be used parallel to the S-FT-LCN method, is the method of Alignment Index. This method is based on some known properties of the time evolution of deviation vectors, described in detail in Voglis, Contopoulos & Efthymiopoulos (1998, 1999).
In particular, let us consider the evolution of two arbitrary different initial deviation vectors ξj1 and ξj2 of the same orbit. If the orbit is chaotic, then the two deviation vectors tend to become parallel or anti-parallel to each other (depending on the initial values of the deviation vectors). This is due to the differential stretching in phase space caused by the effect of the variational equations. This forces the deviation vectors to turn along the direction of a local phase space manifold with the largest eigenvalue. In this case, the difference dj−(t) =∣ξj1(t) −ξj2(t)∣ (parallel deviation vectors) or the sum dj+(t) =∣ξj1(t) +ξj2(t)∣ (anti-parallel deviation vectors), tends exponentially to zero.
On the other hand, in the case of motion on a torus (regular motion), the two deviation vectors, after a transient period, tend to become tangent to the surface of the torus and oscillate with respect to each other. Then, dj− or dj+ have an oscillatory behaviour around a roughly constant mean value. This value is usually close to unity, and at any rate it is not less than 10−3.
Thus, by looking at the evolution of the minimum between the two indices dj−, dj+, called the Smaller Alignment Index, (SALI), or simply Alignment Index, (AI), we can distinguish between regular and chaotic orbits.
The main characteristic of the SALI method is its efficiency in distinguishing regular orbits (with Lyapunov numbers equal to zero) from weakly chaotic orbits with even very small Lyapunov Numbers. Numerical examples of the SALI method, applied to multi-dimensional symplectic maps, were given by Skokos (2001).
3.3 Combining the two methods
As an example of combination of the two methods described above, let us consider the orbits of two particles in the Q model under the Hamiltonian equation (3). We calculate Lj and AI of the two orbits for times tj/Trj= 1200. The results are shown in Figs 2(a) and (b) in log—log scales. The value of Lj of the chaotic orbit tends to a constant value of about 10−2.3 (bold line in Fig. 2a), while the value of Lj of the other orbit (thin line in Fig. 2a) continues to decrease with a logarithmic slope −1 as expected for a regular orbit. It reaches a value of about 10−3 at the end of the integration time. In this case it is not clear whether Lj would be stabilized to any particular value smaller than 10−3 at longer integration times (if the orbit were chaotic), or that it would continue to decrease (if the orbit were regular). For the reasons explained in Section 3.1, we use conventionally the value of 10−3 as a threshold, considering regular all the orbits with S-FT-LCN less than 10−3.
(a) Evolution of the S-FT-LCN (Lj) of the orbits of two particles in the Hamiltonian (3), a chaotic (bold line) and a regular (thin line). (b) The AI of the ordered orbit (thin line) remains almost constant in time while the AI of the chaotic orbit (bold line) reaches the limit of accuracy very quickly.
(a) Evolution of the S-FT-LCN (Lj) of the orbits of two particles in the Hamiltonian (3), a chaotic (bold line) and a regular (thin line). (b) The AI of the ordered orbit (thin line) remains almost constant in time while the AI of the chaotic orbit (bold line) reaches the limit of accuracy very quickly.
The distinction between the two sorts of orbits can also be made by using the AI method. In Fig. 2(b) the bold line gives the evolution of the AI of the chaotic orbit. This index is quite sensitive, reaching the limit of the computer accuracy (≈10−12) at the end of the integration time. In contrast, the AI of the regular orbit (thin line in Fig. 3b) remains above the value of 10−3 all along the integration time. Thus, the distinction made by the AI method is clear and fast.
The evolution of the orbits of the particles in the sample on the plane of the log (AI)—log (L)at (a) 20, (b) 100 and (c) 1200 radial periods. In (d) we plot the Finite Time LCN in common units of (1/Thmct), as derived from the data of (c).
The evolution of the orbits of the particles in the sample on the plane of the log (AI)—log (L)at (a) 20, (b) 100 and (c) 1200 radial periods. In (d) we plot the Finite Time LCN in common units of (1/Thmct), as derived from the data of (c).
Notice that, as seen in Fig. 2(b), the AI of the chaotic orbit (bold line) maintains large values, i.e. above 10−2, for a period of about t/Tr= 200, before jumping to very small values. This means that such an integration time is too short to decide about the character of this orbit. Such a transient effect is quite common for most chaotic orbits, and it is expected as a result of the relatively low level of chaos in our systems.
The two methods are then combined in order to obtain a quite reliable distinction between chaotic and regular orbits on the plane of AI and S-FT-LCN.
Given computer time restrictions, we calculate orbits for a representative sample of particles of the whole system. Namely, we take one in every four particles uniformly distributed along the whole set of particles of the system.
Using as initial conditions the three coordinates and velocities of each particle of the sample at the snapshot (A), we integrate every orbit in the Hamiltonian (3) for a maximum time tj= 1200 Trj. The results are presented in the next section.
4 Results
4.1 The Q model
In Fig. 3(a) the values of the S-FT-LCN, i.e. (Lj) of the particles in the sample versus their AIs are plotted at tj/Trj= 20, in log—log scale. At this time, most of the points appear concentrated in a single group having a triangular shape around a mean value of log (Lj) about −1.4 with log (AI) greater than −3. We call this a ‘group of regular orbits’ or a ‘regular group’. In the same figure there are also a number of points forming a lane emanating from the upper end of the regular group towards smaller values of the AI. The points on the lane correspond to orbits that have just started indicating their chaotic character.
At tj/Trj= 100 (Fig. 3b) the main part of the ‘regular group’ is displaced towards lower values of Lj following the t−1 law of Fig. 2(a) (thin line). However, a good number of points have followed a streaming motion along the lane towards smaller values of AI and larger values of Lj and tend to form another group, i.e. the ‘chaotic group’ with very small values of AI (less than 10−10) as displayed in Fig. 3(b).
As time increases, the number of orbits in the chaotic group increases, but more and more slowly. The situation at tj/Trj= 1200 is shown in Fig. 3(c). The chaotic group is well separated from the regular group. The streaming of points along the lane from the regular group to the chaotic group has become slower, but non-negligible.
The points on this lane correspond to more weakly chaotic orbits, or orbits that are still in a phase of ‘stable chaos’ i.e. a phase of almost constant frequencies (Milani, Nobili & Knezevic 1997; Tsiganis, Varvoglis & Hadjidemetriou 2000, 2002).
Our results of Figs 3(a), (b) and (c) have been tested as regards the following two points.
The computation of the Lyapunov Numbers may be affected by the numerical errors during the integration processes. In the results presented above (Figs 3a, b and c) the orbits are integrated with a tolerance in the conservation of their energy δE/E≤ 10−7. The calculations were repeated for δE/E≤ 10−8. We found that only a number less than 2 per cent of the total number of orbits appeared in a different group in the second case, e.g. ordered instead of chaotic. Thus, the general features of our classification are preserved.
As small time-variations in the self-consistent potential have been ignored by using the autonomous Hamiltonian (3), there is a doubt whether the identities and the number of particles would be the same if another snapshot of the potential is used. As is known (e.g. Habib, Kandrup & Mahon 1997; Kandrup 1998), the presence of noise in the potential may change the character of a number of orbits from regular to chaotic and vice versa.
For this reason we have performed the following test. The run of the self-consistent N-body model is continued for another period of 100 half-mass crossing times creating a new snapshot (B) at the end of this time. We update the coefficients of the potential in the Hamiltonian (3) for the snapshot B and we repeat the calculations of Figs 3(a), (b) and (c). A comparison of the two results shows only small differences in both the number of particles in the two classes (regular or chaotic) and the identities of the particles in each class. In Fig. 4(a) we plot in log—log scale the values of the S-FT-LCN, Lj(A), derived from the snapshot (A) versus the values of the S-FT-LCN, Lj(B), derived from the snapshot (B). The particles in ordered motion are near the left-hand bottom corner, while those in chaotic motion are mainly near the right-hand top corner.
Only a small percentage of particles (less than 3 per cent, plotted away from the diagonal) appear with a different character of their motion in the potential of the snapshot B than in the potential of the snapshot A.
Only a small percentage of particles (less than 3 per cent, plotted away from the diagonal) appear with a different character of their motion in the potential of the snapshot B than in the potential of the snapshot A.
Similarly, the values of the AI are plotted in Fig. 4(b). The ordered motion corresponds to large values of AI (right-hand top corner) and the chaotic motion to small values of AI (left-hand bottom corner). The identities of the majority of the particles that belong to one or the other class are preserved, except of a fraction less than 3 per cent of the total number of particles that have changed the character of their motion and they are plotted away from the diagonal in these diagrams.
Thus, within an uncertainty of less than ∼3 per cent we can use the results of Fig. 3(c) for distinguishing the most chaotic orbits in the Q model. If the points along the lane of weakly chaotic particles are counted with the chaotic class, the threshold in the values of Lj is taken 10−2.8. This is the value of Lj at which the lane separates from the regular group in Fig. 3(c). We also consider a threshold of the AI equal to the value of 10−3. (We found that smaller values are quite improbable in regular motion.) As a conclusion, we find that the fraction of the detected chaotic orbits (hereafter, the chaotic component) corresponds to about 32 per cent of the total mass. The rest 68 per cent of the mass is considered in regular motion (the ordered, or regular component).
However, if the points along the lane of weakly chaotic particles are not counted with the chaotic orbits, we find that the ‘strictly chaotic’ orbits (AI ≤ 10−10) are about 26 per cent.
In order to find the values of the rate of the exponential divergence that are dynamically comparable to each other as regards the chaotic diffusion of the orbits within a Hubble time, we convert Lj in units of the inverse half-mass crossing time Thmct of the system, to find the CU-FT-LCN, as described in Section 3.1, i.e. Lcu=LjThmct/Trj. Thus, the results of Fig. 3(c) are converted to the results shown in Fig. 3(d).
The Lcu of the detected chaotic orbits range between 10−4.6 and 10−1.4, with a preference above the value of 10−3. It is obvious that a number of orbits that have been characterized as chaotic (with Lj≥ 10−2.8) have values of Lcus much smaller than the minimum value of Lj, because of their long radial periods. These small values of the Lcu describe the very small diffusion in a Hubble time and they do not give a measure of chaoticity of these orbits.
The distribution of all the particles along the energy axis is shown in Fig. 5 by the solid line. The dashed line gives the distribution of the particles of the detected chaotic component. In the region of small energies (i.e. below the energy level of about −60) no chaotic orbits were detected by this threshold of chaoticity.
The distribution of the total mass (solid line) and of the detected chaotic part (dashed line) along the energy axis.
The distribution of the total mass (solid line) and of the detected chaotic part (dashed line) along the energy axis.
The majority of the chaotic orbits with large binding energies are radial orbits travelling partly in the regime of the bar-like potential and partly in the outer, almost spherical, potential. The particles in these orbits spend most of their time in large radii and they mainly contribute in forming the halo of the galaxy.
4.2 The surface density profiles of the Q model
The projections of the particles in the ordered component on the x–z plane are shown in Fig. 6(a), while the projections on the same plane of the particles of the chaotic component are shown in Fig. 6(b), for comparison. A similar comparison is made in Figs 6(c) and (d) regarding the projections of the same two components on the y–z plane. It is obvious from these figures that the surface density profiles of the ordered and the chaotic components along a slit along the short axis of the projections are quite different.
The Q model. Projection on the x–zplane of the particles in ordered motion (a) and of the particles in chaotic motion (b). Panels (c) and (d) are as in (a) and (b), but on the y–z plane.
The Q model. Projection on the x–zplane of the particles in ordered motion (a) and of the particles in chaotic motion (b). Panels (c) and (d) are as in (a) and (b), but on the y–z plane.
The differences of the two profiles are quantitatively seen in Figs 7(a), (c) and (e). In these figures the log—log diagrams of the surface density profiles σ(r) on the three planes (x–y, x–z and y–z, respectively) are given as functions of the distance (r) along a slit along the short axis in every projection. The solid line gives the profile of the total mass, the dashed line gives the profile of the ordered component and the dotted line gives the profile of the chaotic component of the galaxy.
Left column: surface density profiles along a slit on the short axis of each projection in log—log scale. The plane of projection is indicated in the figures. Right column: the logarithmic slope of the surface density as a function of the distance along the same slit. Solid lines refer to the total mass, the dashed lines to the mass in ordered motion and the dotted lines to the mass in chaotic motion. The different slopes of the two components create a hump at about the half-mass radius of the system.
Left column: surface density profiles along a slit on the short axis of each projection in log—log scale. The plane of projection is indicated in the figures. Right column: the logarithmic slope of the surface density as a function of the distance along the same slit. Solid lines refer to the total mass, the dashed lines to the mass in ordered motion and the dotted lines to the mass in chaotic motion. The different slopes of the two components create a hump at about the half-mass radius of the system.
In these figures we see that in the three projections the surface density profiles of the chaotic component are considerably flatter than the surface density profiles of the ordered component. There is a radius rc close to the half mass radius r= 1 at which these two profiles cross each other. At smaller distances, the surface density of the total mass σ(r) is mainly determined by the ordered component, the ‘regime of ordered motion’. It is easy to check by taking a number of Poincare surfaces of section, as in Contopoulos et al. (2002), that the main types of regular orbits in this regime are box or box-like orbits (i.e. boxes or inner long axis tubes). At distances larger than rc the surface density of the total mass is determined by the chaotic component, the ‘regime of chaotic motion’ (compare the lines of the three types in Figs 7a, c and e).
A direct implication of the different profiles of the ordered and the chaotic components along the minor axis of the projection in this model is that the chaotic component imposes limitations on the ellipticity of the equidensity contours of the total mass. The small asphericity of the chaotic component reduces the elongated distribution formed by the ordered component and creates a system with maximum ellipticity not exceeding the value of E7. The decrease of the maximum ellipticity (from E7 towards E5 beyond a certain radius, Section 2) is due to the effect of the chaotic component. It seems that, in principle, box or box-like ordered orbits alone can produce surface distributions with ellipticities even larger than the E7 limit. The chaotic component is a good reason to explain why ellipticities larger than E7 are not observed in real elliptical galaxies.
Another implication of the transition from the regime of ordered motion to the regime of chaotic motion is a ‘hump’ formed in the total profile of the surface density taken along the short axis of the projected system.
The logarithmic slope s(r) of the profile of the total mass is a sensitive measure of this ‘hump’. It is characterized by a fast decrease to a local minimum at the ‘hump’, i.e. at the distance rc where the profile of the chaotic and the regular components intersect. This minimum corresponds to the distance up to which motion occurs mainly in box or box-like orbits. Then it rises again and falls to another minimum etc., following the variations of the density of the chaotic component (Figs 7b, d and f).
It must be noted that variations of s(r) along the radius can appear in both the ordered and the chaotic components. There are other reasons as well, which could form ‘humps’ in the profile of σ(r). For example, bunches of 1:1 tube orbits. Generally speaking, a hump must be a signature of passing to the regime of new types of orbits.
We have checked, for example, that the pronounced increase of s(r) for the regular orbits in Fig. 7(f) (dashed line) around the value y= 2 is due to orbits having the Jx component of their angular momentum quite larger than Jy and Jz. This is an indication that the corresponding maximum is due to a bunch of short axis tube orbits.
Notice that the distinction between the profiles of ordered and chaotic components is more clear on the projection of the smallest—largest-axes plane (x–z), or on the projection of the smallest—intermediate-axes plane (x–y). The two profiles are less separated on the plane of the intermediate—largest-axes (y–z).
If we consider the surface density profiles along a slit along the major axis of each projection, the profiles of the chaotic component are again always flatter than the profiles of the ordered component, as in Figs 7(a), (c) and (e), but the transition distance, where the two profiles cross each other, is now larger by a factor of 2 or 3.
It is worth noticing that the profiles of σ(r) and s(r) have only small variations in time. We found that the profiles derived from the data of the snapshot B are only slightly different, especially inside the radius of r= 3. Also, these profiles do no change appreciably if the orbits with Lcu < 10−3s are counted with the regular orbits.
4.3 The C model and comparison with the Q model
The asphericity of the C model is considerably less than in the Q model, and therefore it is expected to be closer to an integrable system and chaos is expected to be less.
Figs 8(a) and (b)refer to the results of the C model and they are analogous to Figs 3(c) and (d), respectively. If we separate the ordered from the chaotic component, using the same criteria as in the Q model, we find that at the time of tj/Trj= 1200 the chaotic component is 23, 5 per cent of the total mass, if the points along the lane are included. The strictly chaotic component, i.e. without the points on the lane, is 14.2 per cent of the total mass.
In Figs 9(a) and (b) the time-evolution of the percentage of the orbits that belong to the chaotic component at a given time is shown, for the Q model (Fig. 9a) and for the C model (Fig. 9b). The dashed line gives the percentage of the chaotic part including the orbits on the lane, while the solid line is without these orbits. In these figures we see that the fraction of the orbits on the connecting lane is about 6 per cent in the Q model and about 9.3 per cent in the C model. This means that there is still a considerable flux of orbits along the lane towards the small values of the AI.
(a) The evolution of the percentage of orbits characterized as chaotic in the Q model. The solid line corresponds to the ‘strictly’ chaotic part (AI < 10−10) while the dashed line includes the orbits on the lane. (b) As in (a) but for the C model. The different rates of growth reflects the difference in the level of chaos.
(a) The evolution of the percentage of orbits characterized as chaotic in the Q model. The solid line corresponds to the ‘strictly’ chaotic part (AI < 10−10) while the dashed line includes the orbits on the lane. (b) As in (a) but for the C model. The different rates of growth reflects the difference in the level of chaos.
As we have already said, the fraction of the chaotic orbits has not been exhausted and there must be a number of chaotic orbits with smaller Ljs and Lcus that have been considered as regular for the lifetime of a galaxy. In the Q model their number can not be very large, as the shape of the curves in Fig. 9(a) shows that they tend to a saturation limit. It is clear that in the C model the saturation limit is at much longer times than in the Q model. The differences between the curves of Fig. 9(a) and the curves of Fig. 9(b) reflect the differences in the level of chaos found in the two models.
The most interesting remark, however, is that the fraction of the chaotic orbits with Lcu>10−2 in the Q model is about 26 per cent of the chaotic component, that is, about 8 per cent of the total mass and in the C model about 9 per cent of the chaotic component, that is, about 2 per cent of the total mass (compare Fig. 3d and Fig. 8b). For these fractions of mass, chaotic diffusion throughout the phase space of the galactic models is possible within a Hubble time, but they are only small fractions of the total mass to alter the general features of a distribution function based on third integrals of motion (Contopoulos 1960; Contopoulos et al. 2000), particularly in the C model.
Another remarkable difference between the two models is seen in Figs 10(a) and (b). In these figures the ratio of the number of the detected chaotic orbits to the total number of orbits at every energy bin along the energy axis is plotted. In the Q model (Fig. 10a) this ratio increases almost monotonically. In contrast, in the C model (Fig. 10b) this ratio has a pronounced maximum at the energy level of ≈−40 followed by a minimum around the energy level of ≈−30.
(a) The ratio of the number of chaotic orbits to the total number of orbits at every bin of energy along the energy axis in the Q model. This ratio increases with the energy almost monotonically. (b) The ratio of the number of chaotic orbits to the total number of orbits at every bin of energy in the C model. This ratio is not monotonic as in (a). It has a pronounced maximum around E =−40 and a pronounced minimum around E =−30 due to the fact that the C model contains a good number of 1:1 tube orbits that are responsible for the small asphericity of the model.
(a) The ratio of the number of chaotic orbits to the total number of orbits at every bin of energy along the energy axis in the Q model. This ratio increases with the energy almost monotonically. (b) The ratio of the number of chaotic orbits to the total number of orbits at every bin of energy in the C model. This ratio is not monotonic as in (a). It has a pronounced maximum around E =−40 and a pronounced minimum around E =−30 due to the fact that the C model contains a good number of 1:1 tube orbits that are responsible for the small asphericity of the model.
The explanation of this difference is as follows: In the C model we find a good number of ordered 1:1 tube orbits, at the energy levels around −30. The role of these orbits, which are automatically placed there, together with the chaotic orbits that have an almost spherical distribution, is to balance the effect of the box orbits or box-like orbits, so that the asphericity of the system is reduced to the value required by the self-consistent equilibrium.
In the Q model, at this energy level, there is an extensive stable area in phase space corresponding to the 1:1 tube orbits, but it is almost empty, i.e. it is occupied by only a small number of orbits (Contopoulos et al. 2000, 2002). At the same energy levels the Q model possesses a good number of chaotic orbits. These chaotic orbits are flexible to follow either boxy or circular geometries but they are almost spherically distributed. They can only reduce the strongly anisotropic distribution of matter due to the box or box-like orbits and prevent larger departures from sphericity. Thus the system maintains its self-consistent equilibrium at large values of ellipticity but not larger.
For the C model, the surface density profiles σ(r) and the logarithmic slopes s(r) of the total mass, as well as of the chaotic and the regular components, along the minor axis of the projection on the x–y plane, are given in Figs 11(a) and (b) (analogous to the Figs 7a and b of the Q model). We see that the profile of the chaotic orbits (dotted line) is again flatter than the profile of the regular orbits (dashed line). However, beyond some distance, the difference is small and the profile of chaotic orbits crosses the profile of ordered orbits at a large distance of about x= 2 where the corresponding hump appears. The two components appear more mixed in this case. This is because, as we have seen, chaos in the C model is less, both either as regards the values of the Lyapunov numbers or the total amount of mass in chaotic motion.
As in Figs 7(a) and (b)but for the C model. The difference between the profiles of the chaotic and the regular components is smaller. A hump is seen at ∣x∣≃ 2.
As in Figs 7(a) and (b)but for the C model. The difference between the profiles of the chaotic and the regular components is smaller. A hump is seen at ∣x∣≃ 2.
5 Conclusions
Two self-consistent (N-body) galactic models of elliptical galaxies have been studied as regards their content of particles in chaotic motion. A combination of two methods has been used for this purpose, that is the S-FT-LCN, Lj, and the AI, explained in Sections 3.1 and 3.2. The most chaotic component, i.e. particles moving in orbits with Lj≥ 10−2.8, have been identified, and their distribution is compared with the distribution of the rest of the mass. Our main conclusions are as follows.
In the triaxial model of a non-cuspy elliptical galaxy with large ellipticity (the Q model), the most chaotic component is found to be about 32 per cent of the total mass. This component, although weakly chaotic, has an almost spherical distribution in space, while the distribution of the rest of the mass forms a strongly triaxial bar (in which 1:1 tube orbits are very few). The chaotic component imposes an upper limit on the maximum ellipticity of the system. The existence of such a component is probably responsible for the upper limit E7 observed in real elliptical galaxies. Only a part of the chaotic component, less than about 8 per cent of the total mass, can produce chaotic diffusion within a Hubble time.
In the other modestly triaxial model with small ellipticity (the C model), the mass in chaotic motion is about 23 per cent of the total mass, but these orbits are more weakly chaotic than in the previous case. No more than 2 per cent of the total number of particles can produce chaotic diffusion within a Hubble time. The detected chaotic component in collaboration with the contribution of the regular 1:1 tube orbits counteract the effect of the box or box-like orbits to form a nearly spherical system.
The surface density profile along a particular axis of the chaotic component always appears flatter than the surface density profile of the ordered component. Because of this difference, the two profiles can cross each other and create a hump in the surface density profile of the total mass. Such a hump is more observable when the profiles are taken along the short axis of the projection and appears at distances around the half-mass radius of the system or larger. It is more clearly seen in the model of higher ellipticity.
Acknowledgments
We thank Professor G. Contopoulos and Dr Ch. Efthymiopoulos for their comments and helpful discussions. We also thank Dr L. Hernquist as well as Drs A. Allen, P. Palmer and J. Papaloizou for their codes. IS and CK wish to thank the Greek State Scholarship Foundation (IKY) for financial support.
