{'fraction_non_alphanumeric': 0.0367686119754117, 'fraction_numerical': 0.01624041891174015, 'mean_word_length': 4.348011363636363, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Orbital Dynamics of Triaxial Black-Hole Nuclei

Orbital motion in triaxial nuclei with central point masses, representing supermassive black holes, is investigated. The stellar density is assumed to follow a power law, rho ~ 1/r^gamma, with gamma=1 or gamma=2. At low energies the motion is essentially regular; the major families of orbits are the tubes and the pyramids. Pyramid orbits are similar to box orbits but have their major elongation parallel to the short axis of the figure. A number of regular orbit families associated with resonances also exist, most prominently the banana orbits, which are also elongated with the short axis. At a radius where the enclosed stellar mass is a few times the black hole mass, the pyramid orbits become stochastic. The energy of transition to this ``zone of chaos'' is computed as a function of gamma and of the shape of the stellar figure; it occurs at lower energies in more elongated potentials. Our results suggest that supermassive black holes may place tight constraints on departures from triaxiality in galactic nuclei, both by limiting the allowed shapes of regular orbits and by inducing chaos.


Introduction
An important change in our thinking about galaxy dynamics took place during the last decade, when it was recognized that the central densities of early-type galaxies and spheroids are generically very high. Evidence for large central masses came from high-resolution kinematical studies of nuclear stars and gas, which revealed the presence in roughly a dozen galaxies of compact dark objects with masses 10 6.5−9.5 M ⊙ , presumably supermassive black holes (Ford et al. 1998). Observations with HST also demonstrated very high stellar densities at the centers of early-type galaxies (Crane et al. 1993;Ferrarese et al. 1994). Low luminosity ellipticals have density profiles that increase as unbroken power laws at small radii, ρ ∼ r −γ , with γ ≈ 2. But Kormendy pointed out already in 1985 that bright galaxies also have central brightness profiles that deviate systematically from that of an isothermal core. The significance of this deviation was not recognized for ten more years due to an optical illusion associated with projection onto the plane of the sky. A luminosity density that varies as r −γ at small radii generates a power-law cusp in projection only if γ > 1. When γ = 1, the central surface brightness is logarithmically divergent (e.g. Dehnen 1993, Fig. 1), and the surface brightness profile differs only subtly from that of a galaxy with an isothermal core. Nonparametric deprojection of the luminosity profiles of bright galaxies (Merritt & Fridman 1995;Gebhardt et al. 1996) revealed that they too harbor power-law cusps but with indices γ < ∼ 1.
While the origin of the power-law cusps is not clear, there are hints that they may be associated with black holes. Steep cusps, γ ≈ 2, form naturally in stellar systems where the black holes grow on time scales long compared to crossing times (Peebles 1972;Quinlan, Hernquist & Sigurdsson 1995). No universally accepted model has yet been proposed for the origin of the weak cusps, but we note here one feature that suggests a link to black holes. Two galaxies, NGC 3379 and M87, have weak cusps with well-determined structural parameters and also have black holes with accurately determined masses. Table 1 gives for each galaxy the "break" radius r b at which the central power law turns over to a steeper outer profile; and also the stellar mass M * contained within r b . In both galaxies, M * (r b ) is identical within the uncertainties with M • , the mass of the black hole; this is in spite of a factor of ∼ 30 difference in M • and ∼ 6 in r b . This rough equality, and the exclusive association of weak cusps with bright galaxies, is consistent with models in which weak cusps are produced following galaxy mergers by the ejection of stars from the nucleus by binary black holes (Ebisuzaki, Makino & Okumura 1991;Makino 1997;Quinlan & Hernquist 1997).
In a fixed stellar potential, the gravitational influence of the black hole is limited to stars with pericenters r p < ∼ r g where M * (r g ) = M • . In an axisymmetric galaxy, r p is bounded from below by the orbital angular momentum L z about the symmetry axis and only a small fraction of the stars, most of which are confined to the nucleus, can be strongly affected by the black hole. But the gravitational influence of a black hole can extend far beyond the nucleus in a non-axisymmetric galaxy since orbital angular momenta are not conserved and stars with arbitrarily large energies can pass close to the center (Gerhard & Binney 1985). In a triaxial potential containing a central point mass, the phase space divides naturally into three regions depending on distance from the center. In the innermost region, r < ∼ r g , the potential is dominated by the black hole and the motion is essentially regular. Sridhar & Touma (1999) and Sambhus & Sridhar (2000) demonstrated the regularity of the motion in black hole nuclei with constant stellar densities, and Merritt & Valluri (1999) found similar results for motion in triaxial models with weak cusps, γ = 0.5. Here we extend those results to cusps with the steeper power laws characteristic of most galaxies. At intermediate radii, the black hole acts as a scattering center rendering almost all of the center-filling orbits stochastic. This "zone of chaos" extends outward from a few times r g to a radius where the enclosed stellar mass is roughly 10 2 times the mass of the black hole (Merritt 1999). In the outermost region, the phase space is a complex mixture of chaotic and regular trajectories, including resonant box orbits that remain stable by avoiding the center (Carpintero & Aguilar 1998;Papaphilippou & Laskar 1998;Valluri & Merritt 1998;Wachlin & Ferraz-Mello 1998).
The focus of the present study is on the inner two regions and specifically on the transition from ordered motion near the black hole to chaotic motion at r > ∼ r g . As noted above, the break radius r b is of order r g in the two bright elliptical galaxies where both radii can be accurately measured. The sudden change in the orbital behavior near r g might therefore imply a change in the three-dimensional shapes of galaxies near r b . In fact there is some evidence for changes in ellipticity and boxiness in bright elliptical galaxies at r ≈ r b (Ryden 1999;Quillen 1999;Bender & Saglia 1999). The work presented here is a prelude to full self-consistency studies which will place more rigorous constraints on the allowed shapes of triaxial black-hole nuclei.
In §2 we present our model for the stellar density and calculate the gravitational potential and forces. The families of orbits are discussed in §3, and the transition from regular to chaotic motion is discussed in §4. §5 sums up.

Mass Model, Potential, Forces
We model the stellar nucleus as a triaxial spheroid with a power-law dependence of density on radius: within a bounding ellipsoid m = m max . The isodensity surfaces are concentric ellipsoids with fixed axis ratios a : b : c; without loss of generality, we assume a > b > c. Triaxiality is measured by the quantity T which is defined as so that prolate galaxies (b = c) have T = 1 and oblate galaxies (a = b) have T = 0. Since the models are scale-free, we are free to assign a mass of 1 to the central point representing the black hole. Most of the discussion below will be restricted to models with γ = 1 ("weak cusp") and γ = 2 ("strong cusp"), typical of the values of γ in bright and faint galaxies respectively.
Our assumption of a power-law density dependence with fixed index γ is reasonable for galaxies with strong cusps, in which ρ ∼ r −γ , γ ≈ 2, even at radii well outside the sphere of influence of the black hole. In weak-cusp galaxies, the shallow inner power law (γ ≈ 1) eventually turns over to a steeper dependence at radii r > ∼ r b . However, as argued above, r b ≈ r g , and we show below that r g is approximately the radius at which a transition to chaos occurs. Thus our γ = 1 models are expected to yield an accurate description of the dynamics of weak-cusp nuclei out to at least the inner edge of the "zone of chaos." The gravitational potential corresponding to the stars can be obtained using Chandrasekhar's theorem (Chandrasekhar 1969 , p. 52, theorem 12) which says that for a density law that is stratified on similar ellipsoids, the gravitational potential can be written as and For the weak-cusp (γ = 1) case, we have while for the strong-cusp (γ = 2) case, The constant terms, depending on m max , will be ignored in what follows since they have no effect on the forces.
For convenience of numerical calculation, the potential in the strong-cusp case may be expressed in terms of a new set of coordinates {r, µ * , ν * } which have the following definitions (de Zeeuw & Pfenniger 1988): where In terms of these variables, the stellar potential in the strong-cusp case can be written as where is the Carlson elliptic integral (Carlson 1988).
Forces may be obtained in analytical form in Cartesian coordinates for both the strong and weak cusp cases. By taking partial derivatives of (6), the weak-cusp force components are found to be where For the strong-cusp forces, we take partial derivatives of (10) in Cartesian coordinates (de Zeeuw & Pfenniger 1988) and the force components are given by where and is the Carlson elliptic integral.
The magnitude of the radial force in the weak cusp case is constant as a function of distance from the center, while in the strong cusp case, the force diverges as r −1 .
The dynamical time T D (E) is defined below as the period of a circular orbit of the same energy in the equivalent spherical potential, which is defined to have a scale length a ave = 3 √ abc. The energy of a circular orbit in the spherical models is and its period is Given the energy E in the triaxial potential, we set E c (r c ) = E and solve for r c and T c (r c ). The latter is equated to T D (E).
Equations (18a) and (18b) are based on the same choices for the zero point of the potential as was made in equations (6) and (7).
The largest Liapunov exponent was computed for all orbits in the standard way, by integrating the equations of motion of an infinitesimal perturbation. Analytical partial derivatives of the forces may be found for both the weak and strong cusp cases. In the strong cusp case, the expressions may be simplified using the identity where The equations of motion were integrated using the routine "RADAU" of Hairer & Wanner (1996). RADAU is a variable time step, implicit Runge-Kutta scheme which automatically switches between orders of 5, 9 and 13. Energy conservation was extremely good; energy was typically conserved to a few parts in 10 9 over 100 dynamical times.
As found also in earlier studies (e.g. Merritt & Valluri 1996), a histogram of Liapunov exponents of orbits at a given energy evolves toward a characteristic form as the orbital integration time is increased. The regular orbits produce a spike at small values of σ, σT D < ∼ 10 −1.1 , whose location moves toward the left roughly as the inverse of the integration time. The chaotic orbits produce a broader peak with a well-defined maximum, typically at σT D ≈ 0.2, and a tail that extends almost to the regular orbits. The tail corresponds to weakly chaotic orbits that are trapped near regular phase space regions for long periods of time. As the integration time is increased, the histogram tends toward two well-separated and narrow peaks as the chaotic orbits become increasingly indistinguishable. In what follows, the identification of chaotic orbits was based both these histograms and on the configuration-space pictures of orbits.

Orbit Families
Although the stellar distribution in our models is scale-free, the presence of the black hole imposes a scale. We expect the orbital population to change systematically with energy, i.e. with distance from the black hole. For each mass model, we defined a grid of energy values as follows. We first adopted a set of values M ⋆ /M • , the ratio of enclosed stellar mass to black hole mass : 1, 0.0, 0.1, 0.2, · · · 1.7, 1.8}.
Each value of M ⋆ /M • defines an ellipsoidal surface, m = m * , such that We then defined the energy E corresponding to this shell as Table 2 (weak cusp) and Table 3 (strong cusp) give x * , E, log (M ⋆ /M • ), and T D for the mass models with T = 0.5 and c/a = 0.5.
We followed the standard practice (Schwarzschild 1993;Merritt & Fridman 1996) of defining two sets of initial-condition spaces. Stationary start space consists of initial conditions lying on an equipotential surface with zero velocity. X − Z start space consists of starting points in the x − z plane with v x = v z = 0. In an integrable triaxial potential, stationary start space generates box orbits while X − Z start space generates mostly tube orbits. These two start spaces probably contain most of the orbits in reflection-symmetric triaxial potentials (Schwarzschild 1993).
As in any non-integrable potential, orbits may be ranked in a hierarchy depending on their phase-space dimensionality (Merritt & Valluri 1999). Stochastic orbits fill a five-dimensional region and in configuration space populate the entire accessible volume within an equipotential surface. Regular, non-resonant orbits occupy 3-tori and densely fill some more restricted volume. Resonant orbits satisfy a single relation of the form between the fundamental frequencies ω i , where the m i are integers, not all of which are zero. Resonant orbits occupy 2-tori and densely fill sheets in configuration space; when stable, they have associated with them families of non-thin orbits which mimic the shape of the parent thin orbit. Orbits satisfying two independent resonance relations are reduced in dimensionality yet again to closed, or periodic, orbits. Periodic orbits are characterized by a single base frequency in terms of which the frequency of motion in any coordinate (e.g. x, y, z) can be expressed as an integer multiple.
We label families of orbits associated with a single resonance by the integer vector (m 1 , m 2 , m 3 ) that defines the resonance. For orbital families associated with doubly resonant, or periodic, orbits we use the notation ν 1 : ν 2 : ν 3 , the ratios of the frequencies in x, y and z. Figures 1 -4 show the major families of orbits and their starting points in our triaxial potentials. Stochastic orbits are present at all energies in both start spaces. They are most prevalent in stationary start space, particularly from starting points near the x− and y− axes. Motion in the vicinity of the z (short) axis tends to be stable. In X − Z start space, stochastic orbits are mostly associated with starting points near the zero-velocity curve, or in a few cases with the transition regions between the different families of tube orbits. As the energy is increased, stationary start space becomes more and more dominated by stochastic orbits; this transition is discussed in more detail in the next section.
Regular motion in stationary start space is dominated by the pyramid orbits. Sridhar & Touma (1999) first described similar orbits, which they called "lenses," in planar, harmonic-oscillator potentials containing central point masses. Merritt & Valluri (1999) demonstrated the existence of the corresponding 3D orbits in triaxial black-hole potentials. Pyramid orbits can be described as Keplerian ellipses with one focus lying near the black hole, and which precess in x and y due to torques from the background stellar potential. They have a roughly rectangular base whose dimensions are fixed by the amplitudes of oscillation in x and y. At low energies, pyramids with a range of shapes exist, having bases elongated parallel to both the x− and y− axes. However their major elongation (more precisely, the elongation of a symmetrical pair of pyramid orbits oriented above and below the x − y plane) is parallel to the z (short) axis. This fact makes pyramid orbits less useful than classical box orbits for reinforcing the shape of the figure.
Close inspection of the pyramid orbits in our numerical integrations reveals many resonant pyramid families. Two of the most important are shown in Figure 1: a (3, 0, −4) resonance between motion in x and z, and a (0, 6, −5) resonance between motion in y and z.
As the energy is increased, a 2 : 1 resonance appears in the narrowest pyramid orbits lying near the short (z) axis. The opening angle of these "banana" orbits increases rapidly with increasing energy, as the stationary point moves along the equipotential surface from the short to the long (x) axis. Many of the orbits from the banana family are found to be associated with a second resonance. Two such, doubly-resonant familes are illustrated in Figure 1: the (2 : 3 : 4) (banana-fish) and (3 : 4 : 6) (banana-pretzel) orbits.
The variation in shape of the bananas with energy is shown in Figure 5. As discussed in the next section, the bananas sometimes persist throughout the "zone of chaos" and sometimes disappear, then reappear at high energies, with their major elongation parallel to the x-axis. Inside of the chaotic zone, their elongation is counter to that of the figure.
One resonant family that is apparently not associated with the pyramids is the (1, −2, 1) family first discussed by Merritt & Valluri (1999). The major elongation of these orbits is parallel to the intermediate (y) axis; they appear both in the weak-and strong-cusp potentials.
The orbit families in X − Z start space are very similar to those in triaxial potentials without black holes (de Zeeuw 1985;Schwarzschild 1993;Merritt & Fridman 1996). Tube orbits avoid the center due to a primary, 1 : 1 resonance in one of the principal planes and are relatively unaffected by the presence of the black hole. The inner long-axis tubes are important only in nearly prolate potentials. The most important resonance among the tube orbits in our potentials is the 2 : 1 resonance in the meridional plane which generates saucer orbits (Lees & Schwarzschild 1992). The saucers appear most prominantly in highly flattened potentials.
We note an interesting feature of the motion in our models. The dynamical roles of the long and short axes at low energies are approximately reversed compared to their roles at large energies, or in triaxial potentials without black holes. The major families of regular, boxlike orbits near the black hole -the pyramids and the bananas -are generated from Keplerian ellipses oriented along the short (z) axis, while in triaxial potentials without central black holes, it is the long (x) axis orbit that generates the boxes and bananas. Similarly, stochastic orbits in our models derive mostly from starting points near the (x − y) plane, while in non-singular potentials the instability strip lies near the (y − z) plane (Goodman & Schwarzschild 1981). This reversal is important because it means that most of the regular orbits near the black hole have the wrong elongation for supporting a triaxial mass distribution.

Transition to Stochasticity
The most dramatic effect of the black hole is to induce a sudden change to stochasticity in stationary start space as the energy is increased. We investigated this transition as a function of γ, c/a and T . Accurate orbital integrations were found to be time-consuming, particularly at low energies and in the strong-cusp model. We therefore used the E10K supercomputer at Rutgers University to distribute the computations over 64 independent processing units. At each energy in each potential, 192 orbits were integrated starting from the equipotential surface for a time of 100T D and the Liapunov exponents were computed. Using all 64 processors, the elapsed time for each set of 192 orbits was ∼ 10 min for γ = 1 and ∼ 30 min for γ = 2.
The results are summarized in Figure 6 and 7. At each energy, the fraction F of the 192 orbits that were found to be stochastic was computed and plotted. Stochastic orbits were identified both by their location in the histogram of Liapunov exponents at a given energy, and by plots of the configuration-space trajectories. While this fraction is not an accurate reflection of the fraction of phase-space associated with chaotic motion, the transition from regularity to chaos is so sudden that there is no need for a more accurate measure.
The basic character of these plots is always the same. At low energies, log (M * /M • ) < ∼ 0, the motion is almost completely regular, consisting mostly of pyramid orbits. Starting at an energy between log (M * /M • ) ≈ 0 and log (M * /M • ) ≈ 0.5, F increases suddenly to ∼ 1 and remains near unity over a range of energies. Finally, at high energies -log (M * /M • ) > ∼ 1 for γ = 1 and log (M * /M • ) > ∼ 1.5 for γ = 2 -F begins to drop and the motion returns to a mixture of regular and chaotic orbits.
The existence of a "zone of chaos" near the centers of triaxial potentials containing black holes was first noted by Merritt & Valluri (1999). Based on our more complete set of numerical experiments, we can make the following statements about how the properties of this zone vary with the parameters of the potential.
1. For a given triaxiality T , chaos sets in at lower energies in more highly elongated models. For instance, for T = 0.5 and γ = 1, F = 0.8 is reached at log (M * /M • ) ≈ 0.3 for c/a = 0.5 and log (M * /M • ) ≈ 0.6 for c/a = 0.8.
2. The transition from F ≈ 0 to F ≈ 1 takes place more rapidly as a function of log (M * /M • ) in the more elongated models.
3. The transition to chaos is interrupted by the appearance of the banana orbits, particularly in the more elongated potentials with γ = 2. For instance, for T = 0.5 and c/a = 0.5, the chaotic fraction first increases to F ≈ 0.7 at log (M * /M • ) = 0.3, then decreases again to ∼ 0.4 at log (M * /M • ) = 0.6 due to the bananas before finally increasing to F ≈ 0.8 at log (M * /M • ) > ∼ 1. The banana orbits in the most highly flattened models (c/a < ∼ 0.6) manage to persist throughout the chaotic zone and keep the chaotic orbit fraction from reaching 100% 4. The transition to chaos depends only weakly on triaxiality T for a given elongation c/a. Figure 7 shows stationary start space as a function of energy for a model (γ = 2, c/a = 0.5) in which the banana orbits persist, with an associated family of regular orbits, throughout the zone of chaos. In mass models where the bananas disappear, the zone of chaos ends with the appearance of a stable resonant family at high energies: typically the 2 : 3 x − z (fish) resonance for γ = 1, and the 2 : 1 x − z (banana) resonance for γ = 2. At still higher energies, the orbital populations are similar to those described by other authors (Carpintero & Aguilar 1998;Papaphilippou & Laskar 1998;Valluri & Merritt 1998;Wachlin & Ferraz-Mello 1998): a complex mix of resonant box orbits, stochastic orbits, and tubes.

Summary
We have investigated the orbital motion in triaxial nuclei with power-law density profiles, ρ ∼ r −γ , γ = (1, 2), and central point masses representing supermassive black holes. The presence of the central point mass divides the phase space into three radial regions. At the lowest energies, the motion is essentially regular. The major orbit families are the tubes, the pyramids, and a number of families associated with resonances, most prominently the 2 : 1 banana resonance. The pyramid orbits are similar in shape to the box orbits of integrable triaxial potentials but have their major elongation parallel to the short axis, making them less useful for reconstructing an elongated figure. At intermediate energies, the tube orbits persist but the pyramid orbits become increasingly chaotic. The transition to a "zone of chaos" occurs rapidly in all of the potentials investigated here, at an energy where the enclosed stellar mass is a few times the mass of the central point.
In the most elongated models with γ = 2, the bananas can persist throughout the zone of chaos; their axis of elongation gradually shifts from the short axis to the long axis of the figure. At higher energies, stable resonant boxlike orbits begin to appear in stationary start space, generated either from closed orbits like the fish or bananas, or from thin orbits corresponding to a 3D resonance.
Our results are limited in their applicability to the central regions of galaxies where the stellar density profile can be approximated as a single power law, ρ ∼ r −γ . In bright elliptical galaxies and bulges, this is the region within r ≈ r b , the break radius, where the central, shallow power law turns over to a steeper dependence at large radii. However we argued (Table 1) that r b is approximately the radius within which the gravitational force from the black hole dominates that from the stars; and the results of §4 show that this is also approximately the radius of transition to the "zone of chaos" induced by the black hole. Thus the onset of chaos in the phase space of real triaxial galaxies should occur at approximately the same radius or energy as calculated here.
Our results highlight two different ways in which central black holes would be expected to limit the degree of triaxiality of real galactic nuclei. First, the regular orbits associated with motion within the "zone of chaos," the pyramids and the tubes, are mostly poorly suited to reinforcing the major elongation of the figure. Second, the black hole induces chaos in the motion of filled-center orbits like the pyramids, causing them to occupy a region that is rounder than that defined by the equidensity contours of the model. We would therefore expect the degree of triaxiality to be limited inside the zone of chaos by the shapes of the regular orbits, and within this region by the lack of regular orbits. These expectations will be tested in a future study where self-consistent triaxial models will be constructed. This work was supported by NSF grant AST 96-17088 and NASA grant NAG 5-6037.    Figure 1. Open circles mark the 1 : 1 closed orbits in the principal planes. Right panels show histograms of Liapunov exponents computed over 100 dynamical times. As the elongation of the model is increased, stochastic and resonant orbits (e.g. the saucers) become more prominent.