{'fraction_non_alphanumeric': 0.027118852711269775, 'fraction_numerical': 0.012531544285098787, 'mean_word_length': 4.554218514973335, 'pattern_counts': {'":': 1, '<': 6, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Orbital Evolution Induced by Baryonic Condensation in Triaxial Halos

Using spectral methods, we analyse the orbital structure of prolate/triaxial dark matter (DM) halos in N-body simulations to understand the processes that drive the evolution of shapes of DM halos and elliptical galaxies in which central masses are grown. A longstanding issue is whether the change in the shapes of DM halos is the result of chaotic scattering of box orbits, or whether they change shape adiabatically in response to the evolving galactic potential. We use orbital frequencies to classify orbits, to quantify orbital shapes, and to identify resonant orbits and chaotic orbits. The frequency-based method overcomes the limitations of Lyapunov exponents which are sensitive to numerical discreteness effects. Regardless of the distribution of the baryonic component, the shape of a DM halo changes primarily due to changes in the shapes of individual orbits within a given family. Orbits with small pericentric radii are more likely to change both their orbital type and shape than orbits with large pericentric radii. Whether the evolution is regular (and reversible) or chaotic (and irreversible), depends primarily on the radial distribution of the baryonic component. The growth of an extended baryonic component of any shape results in a regular rather than chaotic change in orbital populations. In contrast the growth of a massive and compact central component results in chaotic scattering of a significant fraction of both box and long-axis tube orbits. The growth of a disk causes a significant fraction of halo particles to become trapped by major global orbital resonances. Despite the fact that shape of a DM halo is always quite oblate following the growth of a central baryonic component, a significant fraction of its orbit population has characteristics of its triaxial or prolate progenitor (ABRIDGED).

What is the physical mechanism driving shape change? Options suggested in the literature include two possibilities. The first is that the presence of a central mass concentration scatters box orbits that serve as the backbone of a triaxial potential, rendering them chaotic (Gerhard & Binney 1985;Merritt & Valluri 1996). Chaotic orbits in a stationary potential do not conserve any integrals of motion other than the energy E and consequently are free to uniformly fill their allowed equipotential surface. Since the potential is, in general, rounder than the density distribution chaotic diffusion results in evolution to a more oblate or even a spherical shape (Merritt & Quinlan 1998;Kalapotharakos 2008). The second possibility is that the change of the central potential occurs because the growth of the baryonic component causes orbits of collisionless particles in the halo to respond by changing their shapes in a regular (and therefore reversible) manner (Holley-Bockelmann et al. 2002).
Time dependence in a potential is also believed to result in chaotic mixing (Terzić & Kandrup 2004;Kandrup & Novotny 2004) and has been invoked as the mechanism that drives violent relaxation. However, a more recent analysis of mixing during a major merger showed that the rate and degree of mixing in energy and angular momentum are not consistent with chaotic mixing, but rather that particles retain strong memory of their initial energies and angular momenta even in strongly time dependent potentials (Valluri et al. 2007).
One of the principal features of chaotic evolution is irreversibility. This irreversibility arises from two properties of chaotic orbits. First, chaotic orbits are exponentially sensitive to small changes in initial conditions even in a collisionless system. Second, chaotic systems display the property of chaotic mixing (Lichtenberg & Lieberman 1992). Using this principle of irreversibility, D08 argued that if chaotic evolution is the primary driver of shape change, then if, subsequently, the central mass concentration is artificially "evaporated", the system would not be able to revert to its original triaxial distribution. D08 showed that growing baryonic components inside prolate/triaxial halos led to a large change in the shape of the halo. Despite these large changes, by artificially evaporating the baryons, they showed that the underlying halo phase space distribution is not grossly altered unless the baryonic component is too massive or centrally concentrated, or transfers significant angular momentum to the halo. This led them to argue that chaotic evolution alone cannot explain the shape change since such a process is irreversible. They speculated that at most only slowly diffusive chaos occurred in their simulations. Using test particle orbit integrations they also showed that box orbits largely become deformed, possibly changing into tube orbits, during disk growth, but do not become strongly chaotic. D08 employed irreversibility as a convenient proxy for the presence of chaos. In this paper we undertake an orbital analysis of some of the models studied by D08 to better understand the mechanism that drives shape change. Our goal is to understand whether chaotic orbits are an important driver of shape change and if so under what conditions they are important. We also wish to understand how the orbital populations in halos change when a centrally concentrated baryonic component grows inside a triaxial dark matter halo. Finally we would like to understand under what circumstances orbits change their classification.
This paper is organised as follows. In § 2 we describe the simulations used in this paper and briefly describe three models from D08 as well as two additional simulations. In § 3 we describe the principal technique: Numerical Analysis of Fundamental Frequencies (NAFF) that we use to obtain frequency spectra and fundamental frequencies and describe how these frequencies are used to characterise orbits. In § 4 we describe the results of our analysis of five different simulations. In § 5 we summarise our results and discuss their implications.

NUMERICAL SIMULATIONS
We formed prolate/triaxial halos via mergers of systems, as described in Moore et al. (2004). The initially spherical NFW (Navarro et al. 1996) halos were generated from a distribution function using the method described in  with each halo composed of two mass species arranged on shells. The outer shell has more massive particles than the inner one, similar to the method described by Zemp et al. (2008), which allows for higher mass resolution at small radii. Our model halo A was generated by the head-on merger of two prolate halos, themselves the product of a binary merger of spherical systems. The first merger placed the concentration c = 10 halos 800 kpc apart approaching each other at 50 km s −1 , while the second merger starts with the remnant at rest, 400 kpc from an identical copy. The resulting halo is highly prolate with a mild triaxiality. Halo model B was produced by the merger of two spherical halos starting at rest, 800 kpc apart and is prolate, with b/a = c/a ≃ 0.58. Halo A has b/a ≃ 0.45 and c/a ≃ 0.35 while halo B has b/a = c/a ≃ 0.58 (see Figure 3 of D08 for more details). Both halos A and B consist of 4 × 10 6 particles. The outer particles are ∼ 18 times more massive in halo A and ∼ 5 times more massive in halo B. A large part of the segregation by particle mass persists after the mergers and the small radius regions are dominated by low mass particles (cf. Dehnen 2005). We used a softening parameter ǫ = 0.1 kpc for all halo particles. The radius, r200, at which the halo density is 200 times the mean density of the Universe and the total mass within this radius, M200, are given in Table 1.
Once we produced the prolate/triaxial halos, we inserted a baryonic component, either a disk of particles that remains rigid throughout the experiments or softened point particles. The Table 1. The simulations in this paper. M b is the mass in baryons and f b is the baryonic mass fraction. For the particle simulations (P f B2, P l B3, P l A3, P l A4), R b refers to the softening of the spherical baryonic distribution particle(s). For simulation SA1, R b refers to the scale length of the baryonic disk.
given in Table1. In four of the models (P l A3, P l A4, P f B2 and P l B3) the baryonic component is simply a softened point mass with softening scale length given by R b . In model SA1, the density distribution of the disk was exponential with scale length of the baryonic component R b and Gaussian scale-height z b /R b = 0.05. The disk was placed with its symmetry axis along the triaxial halo's short axis in model SA1 (additional orientations of the disk relative to the principal axes were also simulated but their discussion is deferred to a future paper). Initially, the disk has negligible mass, but it grows adiabatically and linearly with time to a mass M b during a time tg. After this time, we slowly evaporated it during a time te. We stress that this evaporation is a numerical convenience for testing the effect of chaos on the system, and should not be mistaken for a physical evolution. The disk is composed of 300K equal-mass particles each with a softening ǫ = 100 pc. From t = 0 to tg + te the halo particles are free to move and achieve equilibrium with the baryons as their mass changes, but all disk particles are frozen in place. The masses of models with single softened particles are also grown in the same way; these are models P f B2, P l B3 from D08, and P l A3 and P l A4 which are new to this paper. D08's naming convention for these experiments used "P" subscripted by "f" for particles frozen in place and by "l" for live particles free to move.
Three different baryonic components are grown in the triaxial halo A: in model SA1 the baryons are in the form of a disk grown perpendicular to the short axis and the model is referred to as Triax+Disk; in model P l A3 the baryonic component is a softened central point mass resembling a bulge and the model is referred to as Triax+Bulg; finally in model P l A4 the baryonic component is a hard central point mass with a softening of 0.1 kpc and the model is referred to as Triax+hardpt. Two different baryonic components are grown in the prolate halo B: in model P f B2 the baryonic component loosely resembles an elliptical galaxy so this run is referred to as Prolt+Ellip; in model P l B3 the baryonic component is a hard central point mass with a softening of 0.1 kpc and is referred to as Prolt+hardpt. For model P f B2, D08 showed that there is no significant difference in the evolution if the central particle is live instead of frozen, all other things being equal. P l A3 was constructed specifically for this paper in order to have a triaxial halo model with a moderately soft spherical baryonic distribution which can be contrasted with the prolate halo model P f B2, while P l A4 is a triaxial halo model which can be contrasted with the prolate halo model P l B3.
For each model there were 5 phases in evolution. The initial triaxial or prolate halo without the baryonic component is referred to as phase a. There is then a phase (of duration tg) during which the halo's shape is evolving as the baryonic component is grown adiabatically. We do not study orbits in this phase since the potential is evolving with time. When the baryonic component has finished growing to full strength, and the halo has settled to a new equilibrium the model is referred to as being in phase b. The baryonic component is "adiabatically evaporated" over a timescale of duration te listed in Table 1. Again we do not study orbits during this period when the potential is evolving with time. After the baryonic component has been adiabatically evaporated completely and the halo has returned to an equilibrium configuration the halo is referred to as being in phase c. We only study the orbits in the halo during the three phases when the halo is in equilibrium and is stationary (not evolving with time).
The growth of the baryonic component induces several changes in the distributions of the DM halo particles: first is an increase in central density relative to the original NFW halo due to increase in the depth of the central potential (an effect commonly referred to as "baryonic compression"); second, the halos become more oblate especially within 0.3rvir. The details of the changes in the density and velocity distributions of DM particles differ slightly depending on the nature of the baryonic component (D08).
All the simulations in this paper, which are listed in Table 1, were evolved with PKDGRAV an efficient, multi-stepping, parallel tree code (Stadel 2001). We used cell opening angle for the tree code of θ = 0.7 throughout 1 . Additional details of the simulations can be found in D08.

Computing Orbits
In each of the halos studied we selected a subsample of between 1000-6000 particles and followed their orbits in each of the three stationary phases of the evolution described in the previous section. The particles were randomly chosen in the halos at t = 0 such that they were inside a fixed outer radius (either 100 or 200 kpc). Since the particles were selected at random from the distribution function, they have the same overall distribution as the entire distribution function within the outer radius selected. We integrated the motion of each a test particle while holding all the other particles fixed in place. We used a fixed timestep of 0.1 Myr and integrated for 50 Gyr, storing the phase space coordinates of each test particle every 1 Myr. We used such long integration times to ensure we are able to obtain accurate measurements of frequencies (as described in the next section). We carried out this operation for the same subset of particles at phases a, b and c. In model SA1 (Triax+Disk) we integrated the orbits of 6000 particles which in phase a were within r = 200 kpc. In model P l A3 (Triax+Bulg) and P l A4 (Triax+hardpt) we considered a subsample of 5000 particles starting within r = 100 kpc. In models P f B2 and P l B3 we considered orbits of 1000 particles within r = 200 kpc. We integrated their orbits as above but we used a smaller timestep δt = 10 4 years in the case of P l A4 and P l B3, which had harder central point masses. The orbit code computes forces in a frozen potential using an integration scheme that uses forces calculated from the PKD-GRAV tree; we used the orbit integration parameters identical to those used for the evolution of the self-consistent models.

FREQUENCY ANALYSIS
In a 3-dimensional galactic potential that is close to integrable, all orbits are quasi-periodic. If an orbit is quasi-periodic (or regular), then any of its coordinates can be described explicitly as a series, where the ω k 's are the oscillation frequencies and the A k 's are the corresponding amplitudes. In a three dimensional potential, each ω k can be written as an integer linear combination of three fundamental frequencies ω1, ω2, ω3 (one for each degree of freedom). If each component of the motion of a particle in the system (e.g. x(t)) is followed for several (∼ 100) dynamical times, a Fourier transform of the trajectory yields a spectrum with discrete peaks. The locations of the peaks in the spectrum correspond to the frequencies ω k and their amplitudes A k can be used to compute the linearly independent fundamental frequencies (Boozer 1982;Kuo-Petravic et al. 1983;Binney & Tremaine 2008). Binney & Spergel (1982, 1984 applied this method to galactic potentials and obtained the frequency spectra using a least squares technique to measure the frequencies ω k . Laskar (1990Laskar ( , 1993) developed a significantly improved numerical technique (Numerical Analysis of Fundamental Frequencies, hereafter NAFF) to decompose a complex time series of the phase space trajectory of an orbit of the form x(t) + ivx(t), (where vx is the velocity along the x coordinate). Valluri & Merritt (1998) developed their own implementation of this algorithm that uses integer programming to obtain the fundamental frequencies from the frequency spectrum. In this paper we use this latter implementation of the NAFF method.
We refer readers to the above papers and to Section 3.7 of Binney & Tremaine (2008) for a detailed discussion of the main idea behind the recovery of fundamental frequencies. For completeness we provide a brief summary here. The NAFF algorithm for frequency analysis allows one to quickly and accurately compute the fundamental frequencies that characterise the quasi-periodic motion of regular orbits. The entire phase space at a given energy can then be represented by a frequency map which is a plot of ratios of the fundamental frequencies of motion. A frequency map is one of the easiest ways to identify families of orbits that correspond to resonances between the three degrees of freedom.
The structure of phase space in 3-dimensional galactic potentials is quite complex and we summarize some of its properties here to enable the reader to more fully appreciate the results of the analysis that follows. When an integrable potential is perturbed, its phase space structure is altered, resulting in the appearance of resonances (Lichtenberg & Lieberman 1992). Resonances are regions of phase space where the three fundamental frequencies are not linearly independent of each other, but two or more of them are related to each other via integer linear relations. As the perturbation in the potential increases, the potential deviates further and further from integrability, and a larger and larger fraction of the phase space becomes associated with resonances.
In a three dimensional potential, orbits that satisfy one resonance condition such as lωx + mωy + nωz = 0 are referred to as "thin orbit" resonances since they cover the surface of a two dimensional surface in phase space (Merritt & Valluri 1999). If two independent resonance conditions between the fundamental frequencies exist, then the orbit is a closed periodic orbit. Orbits that have frequencies close to the resonant orbit frequencies are said to be resonantly trapped. Such orbits tend to have properties similar to that of the parent resonance, but get "thicker" as their frequencies move away from the resonance. At the boundary of the region of phase space occupied by a resonant family is a region called the "separatrix". The separatrix is the boundary separating orbits with different orbital characteristics. In this case it is the region between orbits that have frequencies that are similar to the resonant orbits and orbits that are not resonant. Chaotic orbits often occur in a "stochastic layer" close to resonances and at the intersections of resonances. In fact one of the primary factors leading to an increase in the fraction of chaotic orbits is the overlap of resonances (Chirikov 1979). Chaotic orbits that are close to a resonant family are referred to as "resonantly trapped" or "sticky orbits" (Habib et al. 1997) and are often only weakly chaotic. Orbits that are "sticky" behave like the resonant parent orbit for extremely long times and therefore do not diffuse freely over their energy surface or undergo significant chaotic mixing.
The frequency analysis method allows one to map the phase space structure of a distribution function and to easily identify the most important resonances by plotting ratios of pairs of frequencies (e.g. ωx/ωz vs. ωy/ωz) for many thousands of orbits in the potential. In such a frequency map, resonances appear as straight lines. Stable resonances appear as filled lines about which many points cluster, and unstable resonances appear as "blank" or depopulated lines. The strength of the resonances can be determined by the number of orbits that are associated with them.

Overcoming microchaos in N -body simulations
In N -body systems like those considered in this paper, the galactic mass distribution is realised as a discrete set of point masses. The discretization of the potential is known to result in exponential deviation of nearby orbits, even in systems where all orbits are expected to be regular (Miller 1964;Goodman et al. 1993;Kandrup & Smith 1991;Valluri & Merritt 2000;Hemsendorf & Merritt 2002). However as the number of particles in a simulation is increased, and when point masses are softened, the majority of orbits begin to appear regular despite the fact that their non-zero Lyapunov exponent implies that they are chaotic. Hemsendorf & Merritt (2002) showed that this Lyapunov exponent saturates at a finite value beyond a few hundred particles and corresponds to an e-folding timescale of 1/20 of a system crossing time (for systems with N ∼ 10 5 particles). Despite having large Lyapunov exponents (i.e. short e-folding times) these orbits behave and look much like regular orbits (Kandrup & Sideris 2001;Kandrup & Siopis 2003). This property of N -body orbits to have non-zero Lyapunov exponents has been referred to as "microchaos" (Kandrup & Sideris 2003) or the "Miller Instability" (Hemsendorf & Merritt 2002;Valluri et al. 2007) and suggests that Lyapunov exponents, while useful in continuous potentials, are not a good measure of chaotic behavior resulting from the global potential when applied to N -body systems. This is a strong motivation for our use of a frequency based method which, as we demonstrate, is extremely effective at distinguishing between regular and chaotic orbits and is apparently largely unaffected by microchaos.
We now discuss how frequency analysis can be used to distinguish between regular and chaotic orbits. In realistic galactic potentials most chaotic orbits are expected to be weakly chaotic and lie close to regular orbits mimicking their behaviour for long times. The rate at which weakly chaotic orbits change their orbital frequencies can be used as a measure of chaos. Laskar (1993) showed that the change in the fundamental frequencies over two consecutive time intervals can be used as a measure of the stochasticity of an orbit. This method has been used to study the phase space structure in galactic potentials (Papaphilippou & Laskar 1996Valluri & Merritt 1998). Examples of frequency spectra for each component of motion, and their resolution into 3 fundamental frequencies are given by Papaphilippou & Laskar (1998) for different types of orbits. For each time series the spectrum is analysed and the three fundamental frequencies are obtained. In Cartesian coordinates the frequencies would be ωx, ωy, ωz.
For each orbit we therefore divide the integration time of 50 Gyr time into two consecutive segments and use NAFF to compute the fundamental frequencies ωx, ωy, ωz (note that all frequencies in this paper are in units of Gyr −1 , therefore units are not explicitly specified everywhere). We compute the three fundamental frequencies ωx(t1), ωy(t1), ωz(t1) and ωx(t2), ωy(t2), ωz(t2) in each of the two intervals t1 and t2 respectively. We compute the "frequency drift" for each frequency component as: We define the frequency drift parameter log(∆f ) (logarithm to base 10) to be the value associated with the largest of the three frequencies fx, fy, fz . The larger the value of the frequency drift parameter, the more chaotic the orbit. Identifying truly chaotic behavior however also requires that we properly account for numerical noise. In previous studies orbits were integrated with high numerical precision for at least 100 orbital periods, resulting in highly accurate frequency determination. For instance, Valluri & Merritt (1998) found that orbital frequencies in a triaxial potential could be recovered with an accuracy of 10 −10 for regular orbits and 10 −4 −10 −6 for stochastic orbits using integration times of at least 50 orbital periods per orbit.
In order to use frequency analysis to characterise orbits as regular or chaotic in N -body systems, it is necessary to assess the numerical accuracy of orbital frequencies obtained by the NAFF code. To quantify the magnitude of frequency drift that arises purely from discretization effects (the microchaos discussed above) we select a system that is spherically symmetric and in dynamical equilibrium. All orbits in a smooth spherically symmetric potential are rosettes confined to a single plane (Binney & Tremaine 2008) and are regular. Hence any drift in orbital frequencies can be attributed entirely to discretization errors (including minute deviations of the N -body potential from perfect sphericity). As a test of our application of the NAFF code to N -body potentials we analyse orbits in spherical NFW halos of two different concentrations (c = 10 and c = 20). The halos are represented by 10 6 particles and have mass ∼ 2 × 10 12 M⊙. Particles in both cases come in two species with softening of 0.1 kpc and 0.5 kpc. We carried out the frequency analysis of 1000 randomly selected orbits which were integrated for 50 Gyr in the frozen N -body realisations of each of the NFW halos. Figure 1 shows the distribution of values of log(∆f ) for both halos. In both cases the distribution has a mean value of log(∆f ) = −2.29, with standard deviations of 0.58 (for the c = 10 halo) and 0.54 (for the c = 20 halo). Both distributions are significantly skewed toward small values of log(∆f ) (skewness = -0.85) and are more peaky than Gaussian (kurtosis = 1.95). Despite the fact that the two NFW halos have different concentrations, the distributions of log(∆f ) are almost identical, indicating that our chaotic measure is largely independent of the central concentration.
To define a threshold value of log(∆f ) at which orbits are classified as chaotic we note that 99.5% of the orbits have values of log(∆f ) < −1.0. Since all orbits in a stationary spherical halo are expected to be regular, we attribute all larger values of log(∆f ) to numerical noise arising from the discretization of the potential. Henceforth, we classify an orbit in our N -body simulations to be regular if it has log(∆f ) < −1.0.
To accurately measure the frequency of an orbit it is necessary to sample a significant part of its phase-space structure (i.e. the surface of a 2-torus in a spherical potential or the surface of the 3-torus in a triaxial potential). Valluri & Merritt (1998) showed that the accuracy of the frequency analysis decreases significantly when orbits were integrated for less than 20 oscillation periods. Inaccurate frequency determination could result in misclassifying orbits as chaotic (since inaccurate frequency measurement can also lead to larger frequency drifts). We test the dependence of log(∆f ) on the number of orbital periods np by plotting the frequency drift parameter against the largest orbital frequency (for orbits with np > 20) in both NFW halos in Figure 2 2 . We use ω instead of np since 2 We use the fractional change in the largest of the three fundamental fre- np ∝ ω but is harder to compute accurately. Particles are binned in equal intervals in ω and the error bars represent the standard deviation in each bin. The straight-lines are best fits to the data-points. The slopes of the correlation for the c = 10 halo (solid line) and for the c = 20 halo (dot-dashed line) are both consistent with zero, indicating that log(∆f ) is largely independent of ω (and hence of np). Henceforth we only use orbits which execute more than 20 orbital periods in the 50 Gyr over which they are integrated. The excluded orbits lie predominantly at large radii and are not significantly influenced by the changes in the inner halo that are investigated here. This rejection criterion affects about 25% of the orbits in the triaxial dark matter halos that we consider later.
The effect of central concentration on the accuracy of frequency estimation is of particular concern during phase b, when the potential is deepened due to the growth of a baryonic component. In this phase, frequencies of those orbits which are strongly influenced by the deepened potential are increased. Consequently some orbits execute many more orbital periods during phase b than they do in phase a or phase c. However we have fixed the orbital sampling time period (not integration timestep) to 1 Myr in all phases. In principle coarse time sampling should not be a concern since the long integration time can still ensure a proper coverage quencies measured over two contiguous time intervals (frequency drift) as a measure of chaos (Laskar 1990). For situations where a large fraction of orbits is resonant, it may be more appropriate to use the smallest of the three frequencies or the component with the largest amplitude. of the phase-space torus. To ensure that the sampling frequency per orbital period does not significantly alter the frequency estimation we re-simulated one model (SA1) in phase b and stored the orbits 5 times more often (i.e. at time intervals of 0.2 Myr). We compared the frequencies of orbits computed for the low (ωL) and high (ωH ) time-sampling runs. Figure 3 (left) shows that there is a strong correspondence between frequencies obtained with the two different samplings. We also found (Fig. 3 middle) that the frequency drift parameter log(∆f ) obtained from the two runs are highly correlated although there is some increase in scatter for orbits with values of log(∆f ) > −2. Since the scattered points lie roughly uniformly above and below the 1:1 correlation line, there is no evidence that the higher sampling rate gives more accurate frequencies. The right panel shows that the overall distribution of log(∆f ) is identical for the two runs. We find that 95% of the particles showed a frequency difference < 0.1% between the two different sampling rates. From these tests we conclude that our choice of sampling rate in the phase b is unlikely to significantly affect the frequency measurements of the majority of orbits. We therefore adopt the lower orbit sampling frequency for all the analysis that follows. Carpintero & Aguilar (1998, hereafter CA98) showed that once a frequency spectrum of an orbit is decomposed into its fundamental frequencies, the relationships between the values of the frequencies (ωx, ωy, ωz) can be used to classify the orbits in a triaxial potential into the major orbit families as boxes, long (x) axis tubes and short (z) axis tubes. (CA98 point out that it is difficult to distinguish between the inner long-axis tubes and outer long-axis tubes from their frequencies alone. Therefore we do not attempt to distinguish between these two families with our automatic classification scheme.)

Orbit classification
In addition to classifying orbits into these three broad categories, they showed that if one or more of the fundamental frequencies is an integer linear combination of the other frequencies, the orbit can be shown to be resonant (either a periodic orbit or an open resonance). We followed the scheme outlined by CA98 to develop our own algorithm to classify orbits as boxes, long-axis tubes (abbreviated as L-tubes) and short-axis-tubes (abbreviated as S-tubes) and to also identify orbits that are associated with low-order resonances. We do not describe the classification scheme here since it is essentially identical to that described by CA98, the main difference lies in that we use NAFF to obtain the fundamental frequencies of orbits in the N -body model, whereas they used a method based on that of Binney & Spergel (1984). We tested our automated classification by visually classifying 60 orbits that were randomly selected from the different models. We then ran our automated orbit classifier on this sample, and compared our visual classification with that resulting from the automated classifier. The two methods agreed for 58/60 orbits (a 96% accuracy rate assuming that the visual classification is perfectly accurate). Hereafter we assume that our automated classification is accurate 96% of the time and therefore any orbit fractions quoted have an error of ±4%.

Quantifying orbital shapes
In any self-consistent potential the distribution of shapes of the majority of the orbits match the overall shape of the density profile. The elongation along the major axis is provided either by box orbits or by inner L-tubes. The ratios of the fundamental frequencies of orbits can be used to characterise their overall shape. Consider a triaxial potential in which the semi-major axis (along the x-axis) has a length ax, the semi-intermediate axis has length ay, and the semiminor axis has a length az. The fact that ax > ay > az implies that the oscillation frequencies along these axes are |ωx| < |ωy| < |ωz| for any (non resonant) orbit with the same over-all shape as the density distribution (we consider only the absolute values of the frequencies since their signs only signify the sense of oscillation). We can use this property to define an average "orbit shape parameter" (χs) for any orbit. For an orbit whose overall shape matches the shape of the potential, The orbit shape parameter χs is positive for orbits with elongation along the figure. The larger the value of χs, the greater the degree of elongation along the major axis. Very close to the centre of the potential it is possible for orbits to have greater extent along the y axis than along the x axis, as is sometimes the case with outer Ltubes. For such orbits χs is slightly negative. An orbit for which all frequencies are almost equal would enclose a volume that is almost spherical. For such an orbit, χs ∼ 0 (which we refer to as "round").
Note that orbits which are close to axisymmetric about the short (z) axis (i.e. the S-tubes) also have χs ∼ 0 because ωx ∼ ωy regardless of the value of ωz. Our definition of shape parameter does not permit us to distinguish between truly round orbits for which ωx ∼ ωy ∼ ωz and S-tubes, but both contribute to a more oblate axisymmetric potential.

RESULTS
For every model in Table 1 the three fundamental frequencies ωx, ωy, ωz of each of the orbits in a selected subsample were computed separately in each of the three phases. For each orbit the largest of the three fundamental frequencies is assumed to represent the dominant frequency of motion. The absolute value of this quantity is referred to as the largest fundamental frequency: we use ωa, ω b , ωc to refer to the largest fundamental frequencies of an orbit in each of the three phases a, b, c. In addition to computing the fundamental frequencies over the entire 50 Gyr interval, we split the interval into two equal halves and computed the frequencies in each to compute the frequency drift parameter log(∆f ) defined in §3. All orbits with log(∆f ) < −1.0 are identified as regular and the rest are identified as chaotic. For each orbit we also compute the total energy (E), the absolute value of the total specific angular momentum (|jtot|), the number of orbital periods (np), and the pericenter and apocenter distance from the centre of the potential (rperi, rapo).
In this section we consider results of five simulations, SA1 (Triax+Disk), P l A3 (Triax+Bulg), P l A4 (Triax+hardpt), P f B2 (Prolt+Ellip), and P l B3 (Prolt+hardpt). We shall show that halos A and B have very different initial and final orbital properties despite the fact that their shapes in the presence of baryonic components are very similar (D08).

Distributions of orbital frequencies
The frequency distribution of randomly selected orbits in a triaxial halo can be used to characterise the orbital structure of phase space. It is useful to begin by discussing our expectations for how orbital frequencies change in response to growth of a central baryonic component. The potential is significantly deeper in phase b compared to phase a, consequently the most tightly bound orbits in the initial potential increase their frequencies. In contrast orbits which largely lie outside the central mass concentration do not experience much deformation or much change in their frequencies.
The higher the initial frequency the greater will be the frequency increase. Hence we expect a faster-than-linear increase in frequency in phase b relative to the frequency in phase a.
When the baryonic component is evaporated, the halo expands once more and the halos regain their triaxiality in models SA1, P l A3 and P f B2 but are irreversibly deformed in runs P l A4 and P l B3. One way to investigate the cause of the difference in these behaviours is to look for correlations between largest fundamental frequencies of each orbit in each of the three phases. When the growth of the baryonic component causes an adiabatic change in orbits, one expects that their frequencies ω b will change in a regular (i.e. monotonic) way so that the particles which are deepest in the potential experience the greatest frequency increase. In Figures 4 and 5 we plot correlations between the frequencies ωa, ω b and ωc.
In Figures 4 we show results for the three models whose baryonic scale length is greater than 1 kpc: the left panels show that ω b increases faster-than-linearly with ωa, as expected with fairly small scatter. The right hand panels show that ωc is quite tightly correlated with ωa in all three models with the tightest correlation for simulation SA1 (the dashed line shows the 1:1 correlation between the two frequencies). The deviation from the dashed line and the scatter is only slightly larger in simulations P l A3 and P f B2. The strong correlation between ωc and ωa in these models supports the argument by D08 that the growth of the baryonic component resulted in regular rather than chaotic evolution. In all three models only a small fraction of points deviate from the dashed line for the highest frequencies.
In contrast, the two models with a hard central point mass, P l B3 and P l A4, (Fig. 5) show a non-monotonic change in ω b in response to the growth of the central point mass as well as a higher degree of scattering in frequency space. In particular we note that orbits with small values of ωa (i.e. those which are most weakly bound and have large apocenters) have the largest values of ω b , which is in striking contrast to the situation in Figure 4. We also see that ω b sometimes decreased instead of increasing -again evidence for a scattering in frequency space rather than an adiabatic change.
There is also greater scatter in the right-hand panels pointing to c 2009 RAS, MNRAS 000, 1-24 What fraction of orbits experience a large fractional change in frequency of an orbit from phase a to phase b, and from phase a to phase c? To investigate this we define: The first quantity is a measure of the change in frequency distribution of orbits induced by the presence of the baryonic component, while the latter quantity measures the irreversibility of the evolution following the "evaporation of the baryonic component".  In Figure 6 we plot kernel density histograms of the distribution of the frequency change ∆ω ab (left panel) and ∆ωac (right panel) for orbits in all five models as indicated by the line-legends. Each curve is normalized so that the area under it is unity.
The distribution of ∆ω ab is much wider for models P l A3, P l A4 and P l B3 than for the other two models. In these three models the scale-length of the baryonic component is 1 kpc and results in a broad distribution of ∆ω ab , indicating that orbits over a wide c 2009 RAS, MNRAS 000, 1-24 range of frequencies experience significant frequency change. For model P l A4 (dotted line) the histogram of values of ∆ω ab appears almost flat on the scale of this figure because it is spread out over a much larger range of abscissa values indicating that many more particles are significantly scattered in phase b. The full distribution for P l A4 (see inset panel) is similar in form to P l A3 and P l B3.
The right panels show that only a small number of orbits in models SA1, P l A3 and P f B2 experience an irreversible frequency change ∆ωac > 20%, with the majority of particles experiencing less than 10%. In contrast in models P l A4 and P l B3, the distribution of ∆ωac is much broader: a significant fraction of particles have experienced a large (20-50%) permanent change in their frequencies, reflecting the fact that the models with a hard-compact point mass are the only ones which experience irreversible chaotic scattering.
Are there specific orbital characteristics that contribute to a large permanent frequency change, ∆ωac, between the two triaxial phases? We address this by determining how this quantity relates to other orbital properties. In Figure 7 we plot ∆ωac versus ωa (left panels) and versus ω b (right panels) for four of our models. In the top two panels (SA1, and P f B2 -models with an extended baryonic component) there is no evidence of a dependence of frequency change on ω b and only a slight increase in ∆ωac at the highest values of ωa (results for P l A3 are not shown but are very similar to those for SA1).
On the other hand, the lower two panels (P l A4 and P l B3models with a compact hard baryonic component) show that there is a strong correlation between ∆ωac and orbital frequency ωa indicating that the orbits with the highest frequencies (ωa) experience the largest frequency change, ∆ωac. This is evidence that scattering by the hard central point mass is greatest for particles that are most tightly bound and therefore closest to the central potential, confirming previous expectations (Gerhard & Binney 1985;Merritt & Valluri 1996). The absence of an appreciable correlation with ω b is the consequence of scattering of orbits in frequency.
In Figure 8 we plot ∆ωac versus rperi (left panels) and versus |jtot| (the total specific angular momentum of an orbit averaged over its entire orbit in phase a) (right panels) for orbits in the two models with compact central point mass (P l A4 and P l B3). (We do not show plots for models SA1, P l A3 and P f B2, because they show no correlation between ∆ωac and either |jtot| or rperi.) The left panels of Figure 8 shows that orbits which pass closest to the central point mass experience the most significant scattering. The absence of a correlation with |jtot| however indicates that scattering is independent of the angular momentum of the orbit. In the next section we show that halo A (the initial triaxial halo for model P l A4) is dominated by box orbits while halo B is initially prolate, and is dominated by L-tubes which circulate about the long axis. Contrary to the prevailing view that centrophilic box orbits are more strongly scattered by a central point mass than centrophobic tube orbits, these figures provide striking evidence that chaotic scattering is equally strong for the centrophobic L-tubes that dominate model P l B3 as it is for the centrophilic box orbits that dominate model P l A4. We return to a fuller discussion of the cause of this scattering in § 5.

Changes in Orbital Classification
As we discussed in § 3.2, relationships between the fundamental frequencies of a regular orbit can be used to classify it as a box orbit, a L-tube or a S-tube orbit. Quantifying the orbital composition of the two different halos A and B and how their compositions change in response to the growth of a baryonic component yields further insight into the factors that lead to halo shape change. Orbits were first classified as regular or chaotic based on their drift parameter log(∆f ) as described in § 3.1. Regular orbits were then classified into each of three orbital families using the classification scheme outlined in § 3.2. The results of this orbit classification for each model, in each of the three phases, are given in Table 2.
The most striking difference between the initial triaxial models is that halo A (phase a of models SA1, P l A3 and P l A4) is dominated by box orbits (84-86%) while halo B (phase a of models P f B2 and P l B3) is dominated by L-tubes (78%). (The small differences between models SA1, P l A3 and P l A4 in phase a is purely a consequence of the selection of different subsets of orbits from halo A.) None of the initial models has a significant fraction of S-tubes or of chaotic orbits.
The very different orbit compositions of halos A and B in phase a results in rather different evolutions of their orbital populations in response to the growth of a central baryonic component.
Although the growth of the disk results in a significant decrease in the box orbit fraction (from 86% to 43%) with boxes being converted to either S-tubes or becoming chaotic in phase b, model SA1 is highly reversible suggesting an adiabatic change in the potential. In model P l A3 and P l A4, the more compact spherical baryonic components decrease the fraction of box orbits even more dramatically (from 84% down to 16-17%), pointing to the vulnerability of box orbits to perturbation by a central component. Despite the similar changes in the orbital populations of the two models, P l A3 is much more reversible than model P l A4, indicating that both the shape of the central potential and its compactness play a role in converting box orbits to other families and that the change in orbit type is not evidence for chaotic scattering. It is striking that the more compact point mass in model P l A4 results in significantly more chaotic orbits (21%) compared to 8% in P l A3.
While halo A is initial dominated by box orbits, halo B is initially dominated by L-tubes, which dominate the orbit population in halo B in all three phases. The growth of the baryonic component in phase b causes the box orbit fraction to decrease (especially in model P l B3) while the fraction of chaotic orbits increases slightly. The more extended point mass in P f B2 causes a larger fraction of L-tubes to transform to orbits of another type than does the harder point mass in P l B3, despite the fact that there is much greater scattering in the latter model. A comparison between the model P l A4  and P l B3 show that their orbit populations in the presence of a baryonic component differ significantly due to the different original orbit populations, while their degree of irreversibility is identical (e.g. Fig. 6) since the point mass in the two models is identical.
A significant fraction (21%) of the orbits in phase b of model SA1 and P l A4 are classified as chaotic (orbits with drift rate log(∆f ) > −1.0), in comparison with 9%, 7% and 8% in models P l A3, P f B2 and P l B3 respectively. While the presence of such a large fraction of chaotic orbits in phase b of model P l A4 may be anticipated from previous work, the high fraction of chaotic orbits in SA1 (Triax+Disk) is puzzling. To address concerns about classification error that could arise from errors in the accuracy of our frequency computation, we showed, in Figure 3, that changing the frequency at which orbits were sampled by a factor of five did not result in any change in the overall distribution of log(∆f ), and hence should not affect our classification of orbits as regular or chaotic. Another puzzling fact is that, although model SA1 in phase b has such a significant fraction of chaotic orbits, the orbit fractions essentially revert almost exactly to their original ratios once the disk is evaporated in phase c. Hence, the large fraction of chaotic orbits in phase b do not appear to cause much chaotic mixing. We will return to a more complete investigation of this issue in § 4.4.
In Figure 9 we investigate how orbits of different types (boxes, L-tubes, S-tubes, chaotic) are distributed with rperi, and how this distribution changes from phase a (black curves) to phase b (red curves). In phase a the initially triaxial halo A models (black curves) are dominated by box orbits. The fraction of box orbits is significantly decreased in phase b. In particular box orbits with large rperi are transformed equally to short axis tubes and chaotic orbits, while some box orbits at small rperi are converted to longaxis tubes. In contrast halo B models are dominated by L-tubes in both phases. Rather striking is how little the fraction of L-tubes in the halo B models changes, despite the fact that the halos are significantly more oblate axisymmetric in phase b than in phase a. We saw in Figure 6 that a significant fraction of orbits experience strong scattering that manifests as a change in their orbital frequencies, and in Figure 8 we noted that the orbits with the smallest pericenter radii experience the largest change in frequency. In both models P l A4 and P l B3 the compact central point mass significantly reduced the box orbits. However model P l B3 only has a small fraction (15%) of box orbits and it seems unlikely that the chaotic scattering of this small fraction of orbits off the cen-tral point mass is entirely responsible for driving the evolution of halo shape. Also P f B2 which has a much more extended baryonic component shows a change in orbit population which closely parallels P l B3 and we saw that P f B2 is quite reversible and shows little evidence for chaotic scattering. It is clear (from Fig. 9) that in the prolate models (halo B) the majority of the orbits are L-tubes with large pericenter radii ( rperi ∼ 3 kpc) and these remain Ltubes in phase b. How then do these prolate models evolve to more spherical models while retaining their dominant orbit populations? To address this question we will now investigate the distribution of orbital shapes in each model at each phase of the evolution.

Changes in orbital shape
A parameter, χs, to quantify the shape of an orbit was defined in Equation 3 of § 3.3. Recall that this quantity is positive when the orbit is elongated along the major axis of the triaxial figure, is negative when elongated along the intermediate axis, and almost zero when the orbit is "round" (ωx ∼ ωy ∼ ωz) or roughly axisymmetric about the minor axis (ωx ∼ ωy). In Figure 10 we show the shape distributions for the orbits in four of our five models. For each model we show kernel density histograms for models in phase a (solid curves), phase b (dot-dashed curve), and phase c (dashed curves). In each plot the curves are normalized such that the integral under each curve is unity. We define orbits to be elongated if χs 0.25, and to be "round" if |χs| 0.1.
Before the growth of the baryonic component (phase a: solid curves) the halo A models (left panels SA1, P l A4) have a distribution of orbital shapes that has a large peak at χs ∼ 0.35, arising from elongated orbits and a very small peak at χs ∼ 0 due to round orbits (model P l A3 is not shown but is similar to P l A4.) In halo B models, (right panels P f B2, P l B3) on the other hand, the distribution of shapes is double peaked with about one third of all orbits contributing to the peak at χs ∼ 0. This implies one third of its orbits in the initially prolate halo B are "round". In both halo A and B however, the larger of the two peaks has a value of χs ∼ 0.35 corresponding to quite elongated orbits. Despite the quite different underlying orbital distributions (halo A models dominated by box orbits while the halo B models are dominated by L-tubes). This illustrates that despite having different orbital compositions, a significant fraction of their orbits are similarly elongated.
The dot-dashed curves in all the panels show the distribution of orbital shapes in phase b. In all four models there is a dramatic increase in the peak at |χs| ∼ 0, pointing to a large increase in the fraction of round (or S-tubes) at the expense of the elongated (L-tube or box) orbits. In the halo B models the elongated orbits are significantly diminished indicating that the elongated Ltubes in phase a are easily deformed to "round" orbits in phase b (most likely squat inner-L-tubes). However, in model SA1 there is a large fraction of orbits with intermediate values of elongation 0.1 χs 0.4.
In phase c (dashed curves) all models show the dominant peak shifting back to quite high elongation values of χs ∼ 0.3 (although this is slightly lower than χs ∼ 0.35 in phase a). The downward shift in the peak is most evident in model P l B3 (Prolt+hardpt), which as we saw before, exhibits the greatest irreversibility in shape. The scattering of a large fraction of the orbits by the hard central potential in model P l B3 seen in Figures 7 and 8 is the major factor limiting reversibility of the potential. The smallest shift is for model SA1 (Triax+Disk), which exhibited the greatest reversibility. We can also investigate how the shapes of orbits vary with pericentric radius. We expect that orbits closer to the central potential should become rounder (χs → 0) than orbits further out. We see that this expectation is borne out in Figure 11 where we plot orbital shape parameter χs versus rperi in both phase a (left hand plots) and phase b (right hand plots). In each plot the dots show values for individual orbits. The solid curves show the mean of the distribution of points in each of 15 bins in rperi. Curves are only plotted if there are more than 30 particles in a particular orbital family (P l A4 is not shown since it is similar to P l A3). For models SA1 (Triax+Disk) and P l A3 (Triax+Bulg) the figure confirms that elongated orbits in the initial halo A were box orbits (black dots and curves) and L-tubes (red dots and curves). The S-tubes (blue dots and curves) are primarily responsible for the "round" population at χs ∼ 0. In phase b (right-hand panels) of both SA1 and P l A3 there is a clear tendency for the elongated orbits (boxes, Ltubes and chaotic) to become rounder at small pericenter distances, but they continue to be somewhat elongated at intermediate to large radii. Chaotic orbits in phase b of model SA1 appear to span the full range of pericentric radii and are not confined to small radii. (Note that the density of dots of a given colour is indicative of the number of orbits of a given type but the relative fractions are better judged from Fig. 9 and Table 2.) For phase a in the models P f B2 (Prolt+Ellip) and P l B3 (Prolt+hardpt) (left panels of each plot), boxes and L-tubes are elongated (χs 0.25), except at large rperi 8 kpc where they become rounder. We see a trend for the average orbital shape (as indicated by the curves) in phase b to become round at small pericenter radii.
Note that in all the plots the curves only show the average shape of orbits of a given type at any radius. The points show that in the case of the L-tubes in particular, the red dots tend to be distributed in two "clouds": one with large elongations χs > 0.3 and one with small elongation χs ∼ 0.1.
Thus in all four models it is clear that orbits that are elongated along the major axis of the triaxial potential in phase a become preferentially rounder at small pericenter radii in phase b. It is this change in orbital shape that plays the most significant role in causing the overall change in the shape of the density in the baryonic phase 3 . Figure 10. Kernel density histograms of the distribution of orbital shape parameter χs for each of the four models: SA1 (top left), P l A4 bottom left , P f B2 (top right ) and P l B3 (bottom right) (P l A3 is not shown since it is very similar to P l A4.) Distributions of χs in phase a are shown by solid curves, in phase b by dot-dashed curves, and in phase c by dashed curves. In all models, a large fraction of orbits in phase b are "round" (χs ≃ 0..)

Frequency maps and chaotic orbits
We saw in Table 2 that phase b of model SA1 (Triax+Disk) and of model P l A4 (Triax+hardpt) have a significant fraction (21%) of chaotic orbits (i.e. orbits with log(∆f ) > −1). While P l A4 shows significant lack of reversibility, which we can attribute to the presence of this high fraction of chaotic orbits, model SA1 does not show evidence for irreversibility. Figure 12 shows kernel density histograms of the chaotic drift parameter log(∆f ) for orbits in each of the three phases in model SA1. It is obvious that in phases a and c there is only a small fraction of chaotic orbits (i.e. orbits with log(∆f ) > −1), whereas a much more significant fraction of orbits lie to the right of this value in phase b. Even the peak of the distribution in phase b is quite significantly shifted to higher drift values.
In this section we investigate the surprising evidence that the chaotic orbits in phase b of model SA1 do not appear to mix. One possible reason for the lack of diffusion of the chaotic orbits is that the timescale for evolution is not long enough. Indeed D08 report that evolving run SA1 with the disk at full mass for an additional 5 Gyr after the growth of the disk is complete, leads to a larger irreversible evolution (see their Figure 3a). Nonetheless, even in that case the irreversible evolution was only marginally larger than when the disk was evaporated right after it grew to full mass. Moreover the growth time was 5 Gyr which means that the halo was exposed to a massive disk for a cosmologically long time.
A second possible reason for the lack of chaotic diffusion is that most of the chaotic orbits in this phase of the simulation are "sticky". The properties of "sticky chaotic orbits" were discussed in § 3. In a series of experiments designed to measure the rate of chaotic mixing, Merritt & Valluri (1996) showed that while ensembles of strongly chaotic orbits diffused and filled an equipotential c 2009 RAS, MNRAS 000, 1-24 Figure 11. For models SA1, P l A3, P f B2 and P l B3, the orbital shape parameter χs for each orbit is plotted against its pericentric radius r peri as a small dot. The orbits of each of the four major orbital families are colour coded as in the figure legends. Left hand panels are for phase a and right hand panels are for phase b. The solid curves show the mean value of χs for all particles of that particular family, in 15 bins in r peri . Curves are not plotted if there are fewer than 30 orbits in a given orbital family. We used a kernel regression algorithm to smooth the curves. surface on timescales between 30-100 dynamical times, similar ensembles of "sticky" or resonantly trapped orbits diffused much less quickly and only filled a small fraction of the allowed surface after very long times. Laskar (1990) showed that frequency maps are a powerful way to identify resonances in dynamical systems. Frequency maps are obtained by plotting ratios of the 3 fundamental frequencies for each individual orbit. If a large and representative orbit population is selected, they can provide a map of the phase space structure of the potential including all the resonances. Resonances appear as straight lines on the frequency map since their fundamental frequencies satisfy a condition like lωx + mωy + nωz = 0. This method of mapping the phase space has the advantage that since it only depends on the ratios of the frequencies and not on the frequencies themselves, it can be used to map phase space for large ensembles of particles without requiring them to be iso-energetic. This is a significant advantage over mapping schemes like Poincaré surfaces-of-section, when applied to an N -body simulation where particles, by design, are initialised to be smoothly distributed in energy. Thus one can use the method to identify global resonances spanning a large range of orbital energies in N -body simulations.
In Figure 13 we present frequency maps for four of the five models in phase a (left panels) and phase b (right panels) (P f B2 is not shown since the frequency maps for this model are indistinguishable from those for P l B3). For each orbit the ratios of the fundamental frequencies ωy/ωz and ωx/ωz are plotted against each other. Particles are colour coded by their energy in phase a. The energy range in phase a was divided into three broad energy bins, with equal numbers of particles per bin. The most tightly bound particles are coloured blue, the least bound particles are coloured red and the intermediate energy range is coloured green.
Resonance lines are seen in the clustering of particles in all the maps. The most striking of the frequency maps is that for phase b of model SA1 (Triax+Disk). This map has significantly more prominent resonance lines, around which many points cluster, than any of the other maps. Three strong resonances and several weak resonances are clearly seen as prominent straight lines. The horizontal line at ωy/ωz = 1 corresponds to the family of orbits associated with the 1:1 closed (planar) orbit that circulates around the x-axis, namely the family of "thin shell" L-tubes. The diagonal line running from the bottom left corner to the top right corner with a slope of unity (ωy/ωz = ωx/ωz) corresponds to the family of orbits that circulates about the z-axis: the family engendered by the "thin shell" S-tubes. Since this latter family shares the symmetry axis of the disk, it is significantly strengthened in phase b by the growth of the disk. In addition to having many more orbits associated with it, this resonance extends over a much wider range in energy as evidenced by the color segregation along the resonance line (blue points to the bottom left and red points at the top right). This segregation is the result of an increase in the gradient of the potential along the z-axis due to the growth of the disk, which results in an increase in ωz. The more tightly bound a particle, the greater the increase in ωz, and the greater the decrease in both its ordinate and abscissa. The most bound particles (blue points) therefore move away from their original positions towards the bottom left hand corner of the plot. The least bound particles (red points) are furthest from the center of the potential and these points experience the least displacement -although these points also shift slightly toward the resonance lines.
A third prominent resonance is the vertical line at ωx/ωz = 0.5 that corresponds to orbits associated with the family of banana (1:2 resonant) orbits. This banana (boxlet) resonance is also enhanced by the growth of the disk since this family of orbits, while not axisymmetric, is characterised by large excursions along the x-axis and smaller excursions in the z direction. Several shorter resonance lines are seen but are too sparsely populated in this plot to properly identify.
The frequency map for model P l A3 in phase b shows that a spherical baryonic component produces a rather different phase space structure than that produced by the disk. In particular it is striking that the most tightly bound (blue) points are now clustered at the intersection of the horizontal and diagonal resonances namely around the closed period orbits 1:1:1. This may be understood as the consequence of the growth of the spherical baryonic point mass around which all orbits are rosettes and since no direction is preferred all orbits are "round". The 1:2 banana resonance is also less prominent in this model (largely because the deep central potential destabilises this boxlet family).
The frequency map for model P l A4 phase b shows the greatest degree of scattering, as evidenced by the thickest resonance lines. We attribute this to the large number of chaotic orbits in this model. Apart from the broad clustering of points around the diagonal (Stube) and horizontal (L-tube) resonances there are no strong resonance lines seen in this map. Unlike the map for P l A3 which shows a clustering of tightly bound (blue) points at the 1:1:1 periodic or- Figure 13. Frequency maps of particles in phase a and phase b for four models. For each particle the ratio of the fundamental frequencies ωy/ωz is plotted versus ωx/ωz is plotted by a single dot. The dots are colour coded by the energy of the particle in phase a. The most tightly bound particles are coloured blue, and the least bound particles are coloured red. Model SA1 has 6000 particles, model P l A4 and P l A3 have 5000 particles, while P l B3 has 1000 particles each. bit resonances, the blue points are widely scattered in the frequency map of P l A4.
The frequency maps for model P l B3 shows that most of the orbits in this model are associated with the (1:1) L-tube family (horizontal line). A smaller number of orbits is associated with the 1:1 S-tube resonance (diagonal line). We saw previously that the growth of the baryonic components in this prolate halo caused little change in the orbit families. This is confirmed by the fact that the frequency maps in both phases are remarkably similar except for an increase in the clustering of points at the intersection of the horizontal (L-tube) and the diagonal (S-tube) resonance, which occurs for the same reason as in P l A3. Since halo B is initially a highly prolate model, it has (as we saw previously) only a small fraction of box (and boxlet) orbits and in particular no banana orbits.
It is quite striking that in phase a the frequency maps show significantly less segregation by energy, and only a few resonances. This is because the initial triaxial models were generated out of mergers of spherical NFW halos which were initially constructed so that orbits were smoothly distributed in phase space. The increase in the number of resonances following the growth of a baryonic component is one of the anticipated consequences of resonant trapping that occurs during the adiabatic change in a potential (e.g. Tremaine & Yu 2000;Binney & Tremaine 2008).
To test the conjecture that the majority of chaotic orbits in model SA1 (phase b) are resonantly trapped, we compute the number of chaotic orbits that lie close to a major resonance line. We define "closeness" to the resonance by identifying those orbits whose frequency ratios lie ±α of the resonant frequency ratio. For example we consider an orbit to be close to the (1:1) L-tube resonance (horizontal line in map), if |ωy/ωz − 1| α. We find that the fraction of chaotic orbits in phase b, that lie close to one of the three major resonances identified above, is 51% when α = 0.01 and 62% when α = 0.03. Weaker resonances lines (which are hard to recognise due to the sparseness of the data points) may also trap some of the chaotic orbits. This supports our conjecture that the main reason that model SA1 does not evolve in phase b, despite the presence of a significant fraction of chaotic orbits, is that the majority of the chaotic orbits are trapped around resonances and therefore behave like regular orbits for very long times.
In Figure 14 we plot four examples of chaotic orbits in phase b of SA1, which illustrate how resonantly trapped or "sticky" chaotic orbits look. Each panel of four sub-plots shows a single orbit plotted in two Cartesian projections (side-by-side). The top pair of subplots show the orbit over the first 10 Gyr long time segment, while the bottom pair shows the same orbit over a second 10 Gyr time segment. The two time segments were separated by 10 Gyr. For illustration we selected orbits with a range of drift parameters. The orbit in the top-left panel is an example of an orbit that conforms to our notion of a chaotic orbit that explores more phase space as time progresses, and has a large drift parameter of log(∆f ) = −0.48.
The top-right panel shows a S-tube orbit that suddenly migrates to a box orbit (this orbit has a log(∆f ) = −0.54) and was probably in the separatrix region between the S-tube and box families.
The bottom-left panel shows an orbit that is originally a S-tube that becomes trapped around a resonant boxlet ("fish") family (with log(∆f ) = −0.66), while the bottom right-hand panel shows a weakly chaotic box orbit (with log(∆f ) = −0.94). Of the 21% of orbits in phase b that are chaotic (log(∆f ) −1), only ∼ 5% have log(∆f ) −0.5. This fraction is small enough that one does not expect it to result in significant chaotic mixing.

SUMMARY AND DISCUSSION
Since it was first proposed, the idea that a central black hole would scatter centrophilic box orbits in triaxial galaxies resulting in more axisymmetric potentials (Gerhard & Binney 1985) has frequently been used to explain the shape change in a variety of systems from the destruction of bars by central black holes (Norman et al. 1996) to the formation of more oblate galaxy clusters in simulations with gas .
Experiments on chaotic mixing indicated that the timescales for such mixing is about 30 -100 dynamical times (Merritt & Valluri 1996), which is much longer than the timescales for evolution of dark matter halos in simulations with gas ). In addition recent detailed studies of N -body simulations with controlled experiments have shown that the role of chaotic mixing may be less dramatic than conjectured by these previous studies. A study of relaxation of collisionless systems following the merger of two spherical galaxies showed that despite the fact that a large fraction of the orbits in a system undergoing violent relaxation are chaotic, the timescales for chaotic diffusion and mixing are too long for this process to play a significant role (Valluri et al. 2007). In fact, even after violent relaxation, orbits retain strong memories of their initial energies and angular momenta. D08 argued that since chaos introduces an irreversible mixing, numerical experiments in which evolution is driven by chaotic orbits should not be reversible. These authors studied the macroscopic shapes of triaxial dark matter halos in response to the growth of a baryonic component. Unless the baryons were too centrally concentrated, or transported angular momentum to the halo, the evolution they saw was reversible, from which they concluded that much of the shape change arises from deformations in the shape of individual orbits rather than significant chaotic scattering. In this paper we investigated this issue in significantly greater detail by applying the Numerical Analysis of Fundamental Frequencies (NAFF) technique that allows us to quantify the degree to which chaotic diffusion drives evolution and to identify the primary physical processes that cause halo shape change. The frequency based method is able to distinguish between regular and chaotic orbits, making it more useful than Lyapunov exponents which are known to be sensitive to discretization effects in N -body systems (Hemsendorf & Merritt 2002). We use the method to quantify the drift in frequencies of large representative samples of orbits, thereby quantifying the degree of chaos in the systems we study. It also allowed us to map the phase space structure of the initial and final halos and to quantify the relationship between the change in the shapes of individual orbits and the shape of the halo as a whole.
Applying various analysis methods to orbits in five systems we demonstrated that the conclusion reached by D08 that chaos is not an important driver of shape evolution when the baryonic component is extended is indeed valid. As did D08, we also found that significant chaotic scattering does occur when the baryonic component is in the form of a hard central point mass (of scale length ∼ 0.1 kpc). It is interesting that regardless of the original orbital composition of the triaxial or prolate halo, and regardless of the shape of radial scale length of the baryonic component, halos become more oblate following the growth of a baryonic component. Thus two quite different processes (chaotic scattering and adiabatic deformation) result in similar final halo shapes even in halos with very different orbital compositions.
We explored two different initial halos, one in which box orbits were the dominant elongated population (halo A) and the other in which L-tubes dominated the initial halo (halo B). Despite the different orbit compositions both models exhibit similar overall evolution with regard to the shapes of orbits. In the halo A models, the box orbits were much more likely to change to either L-axis or S-tubes, whereas in the halo B models, the dominant family of L-tubes largely retained their orbital classification while deforming their shapes.
Below we list the main results of this paper: (i) Correlations between the orbital frequencies in the three different phases ωa, ω b and ωc are a useful way to search for regular versus chaotic evolution of orbits. The orbital frequencies in the three phases are found to be strongly correlated with each other when the baryonic component is extended (Fig. 4), but show significant scattering when the baryonic component is a compact point mass (of scale length about 0.1 kpc) (Fig. 5). In the more extended distributions, only a small fraction of the orbits experience significant change in their original orbital frequencies when the baryons are evaporated, while both the magnitude of scattering in orbital frequency as well as the fraction of orbits experiencing scattering, increases as the baryonic component becomes more compact (Fig. 6).
(ii) In the three models with relatively extended baryonic components, the change in orbital frequency between phase c and phase a (∆ωac) is not correlated with orbital frequency (Fig. 7), pericenter distance or orbital angular momentum. When the baryonic component is a hard point mass, however, the frequency change is greater for orbits that are deeper in the potential and therefore have both a higher initial orbital frequency (Fig. 7) and smaller pericentric radius (Fig. 8). Scattering in frequency affects both the centrophilic box orbits as well as centrophobic L-tubes.
(iii) The growth of a baryonic component in halo A (either disk or softened point mass) causes box orbits with large pericentre radii to be converted to S-tubes, L-tubes or become chaotic ( Fig.9 top panel). While this change is almost completely reversible in the case of the disk or a diffuse point mass, it is less so when the baryonic component is a hard point mass. In halo B, which is dominated by L-tubes, the growth of the baryonic component causes almost no change in the orbital composition of the halo, indicating that the Ltubes are not destroyed but deformed ( Fig.9 bottom panel). Even though P l B3 (Prolt+hardpt) is a model with significant orbit scattering by the hard central point mass, the process appears to mainly convert elongated inner L-tube orbits to somewhat rounder outer L-tubes. In model P l A4 (Triax+hardpt) the box orbits are scattered onto S-tubes or chaotic orbits. The significant amount of scattering seen for even centrophobic L-tube orbits shows that the evolution is not due to direct scattering by a central point mass as sometimes assumed. Two alternative possibilities are more likely to account for the significant scattering in frequency. First the change in the symmetry and depth of the central potential is a perturbation to the potential that gives rise to an increase in the region of phase space occupied by resonances (Kandrup 1998 -private communication). As the resonances overlap there is an increase in the degree of chaotic behavior (Chirikov 1979). The second option is that the point mass attains equipartition with the background mass distribution, resulting in Brownian motion (Merritt 2005). The Brownian motion can cause the centre of the point mass to wander within a region of radius ∼ 0.1 − 1 kpc which can result in a significant change in the maximum gravitational force experienced by an orbit from one pericentre encounter to the next. This change in the maximum gravitational force manifests as scattering of the orbit which is equally effective for both box orbits and long-axis tubes. Indeed a small wandering of the central massive point is seen in the N -body simulations of P l A4; this motion is not included in our orbit calculations since all particles are frozen in place when calculating orbits. While it is beyond the scope of this paper to explore this issue further, we caution that the motion of the point mass in our N -body simulations is likely to over-estimate the magnitude of the Brownian motion, since this depends on the mass resolution of the background particles. This suggests that in a real galaxy the evolution of the shape is much more likely to be driven by smooth adiabatic deformation of orbits than chaotic scattering.
(iv) In triaxial halos, the orbital shapes sharply peaked distribution with the most elongated orbits (χs > 0.25) are either boxes or L-tubes. In the prolate halos the second peak at χs ∼ 0 contains a third of the orbits and is composed of squat outer L-tubes and some box orbits. The growth of a baryonic component of any kind causes orbits of all types to become "rounder", especially at small pericenter radii. This change in orbital shape distribution with radius is the primary cause of the change of halo shapes in response to the growth of a baryonic component. This is consistent with the findings of D08 who also found that the orbits in the models became quite round.
(v) The growth of a disk causes a large fraction of halo orbits to become resonantly trapped around major resonances. The three most important resonances are those associated with the 1:1 tube (thin shell) orbit that circulates about the short axis in the x − yplane, the 1:1 tube (thin shell) orbit that circulates about the long axis in the y − z-plane and the 1:2 banana resonance in the x − zplane. We saw from the frequency maps that the resonant trapping of the halo particles depends both on the form of the baryonic component grown in the halo as well as on the initial orbital population of the halo.
Thus we conclude that the evolution of galaxy and halo shapes following the growth of a central component occurs primarily due to regular adiabatic deformation of orbital shapes in response to the changing central potential. Chaotic scattering of orbits may be important particularly for orbits with small pericentre radii but only when the central point mass is extremely compact. Contrary to previous expectations, chaotic scattering is only slightly more effective for centrophilic box orbits than it is for centrophobic L-tubes. Boxes can be scattered onto both L-and S-tube orbits and a significant fraction become chaotic. When the compact central point mass scatters L-tubes as it does in a prolate halo, they are scattered onto other L-tube orbits rather than onto S-tube orbits. The strong chaotic scattering that we see on centrophobic L-tube orbits has not been previously anticipated.
An important implication of our analysis is that while the shapes of halos (and by extension elliptical galaxies) become more oblate (especially at small radii), following the growth of a baryonic component, the majority of their orbits are not S-tubes as might be predicted from their shapes. Instead our analysis shows that orbits prefer to maintain their orbital characteristics, and the majority of the orbits are those which would be generally found in triaxial galaxies. This is particularly important for studies of the internal dynamics of elliptical galaxies since the fact that their shapes appear nearly axisymmetric need not imply that their orbital structure is as simple as the structure of oblate elliptical galaxies. Modifying the shapes to slightly triaxial could result in significant changes in their orbit populations and consequently could affect both the inferred dynamical structure as well as the estimates of the masses components such as the supermassive black holes in these galaxies (van den Bosch & de Zeeuw 2009).
Finally, our finding that the growth of a stellar disk can result in a large fraction of halo orbits becoming trapped in resonances could have important implications for observational studies of the Milky Way's stellar halo. The computation expense of the orbit calculations forced us to restrict the size of the frequency map for model SA1 to 6000 particles. This is only a tiny fraction of the particles in the original simulation. Despite the smallness of the sample, the frequency maps (Fig. 13) shows a rich resonant structure which implies that the particles (either stars or dark matter) in the stellar and dark matter halos of our Galaxy, particularly those close to the plane of the disk, are likely to be associated with resonances, rather than being smoothly distributed in phase space (this is in addition to structures arising due to tidal destruction of dwarf satellites). Although significantly greater resolution is required to resolve such resonances than is currently available, this could have significant implications for detection of structures in current and upcoming surveys of the Milky Way such as SDSS-III (Segue) and Gaia (Perryman et al. 2001;Wilkinson & et al. 2005) and in ongoing direct detection experiments which search for dark matter candidates.

ACKNOWLEDGMENTS
M.V. is supported by the University of Michigan. V.P.D. thanks the University of Zürich for hospitality during part of this project. Support for one of these visits by Short Visit Grant # 2442 within the framework of the ESF Research Networking Programme entitled 'Computational Astrophysics and Cosmology' is gratefully acknowledged. Support for a visit by M.V. to the University of Central Lancashire at an early stage of this project was made possible by a Livesey Grant held by V.P.D. All simulations in this paper were carried out at the Arctic Region Supercomputing Center. We would like to thank the referee Fred Adams for his very thoughtful and constructive report which helped improve the paper.