A UNIFIED FRAMEWORK FOR THE ORBITAL STRUCTURE OF BARS AND TRIAXIAL ELLIPSOIDS

The N-body disk models studied in this paper are almost identical to those presented in Shen & Sellwood (2004, hereafter SS04) and were previously analyzed by Brown et al. (2013), who examined the effects of the growth of SMBHs on the measurement of the observable velocity dispersion within the effective radius (σ_e) and the resultant consequences for the scatter and slope of the well-known scaling relationship between SMBH mass (${M}_{\mathrm{BH}}$) and stellar velocity dispersion (σ): the "${M}_{\mathrm{BH}}\;\mbox{--}\quad \sigma$ relationship." Below we give a brief description of the simulations and refer the reader to the above papers for a more detailed discussion of the setup and simulation process.

In this paper we restrict our orbital analysis to a single set of initial conditions composed of a Kuzmin disk embedded in a static logarithmic halo (for details see SS04). Following standard practice, the units used in the simulations are $G={M}_{d}={R}_{d}=1$, where G is Newton's gravitational constant, M_d is the mass of the disk, and R_d is the initial disk scale radius. By dimensional analysis the unit of time is ${t}_{\mathrm{dyn}}={({R}_{d}^{3}/{{GM}}_{d})}^{1/2}$. In this paper all figures are presented in these units. Physically relevant scalings can be obtained by choosing observationally motivated values for M_d and R_d. For example, ${M}_{d}=5\times {10}^{10}\;{M}_{\odot }$ and R_d = 1 kpc would yield a unit of time ${t}_{\mathrm{dyn}}\;\sim \;$ 2.1 Myr. In these units the semimajor axis length of the bar is about 3–4 kpc in length and the disk has a total radius of about 25 kpc, making it similar to the Milky Way.

The initial disk distribution function consisted of 2.8 × 10⁶ disk particles set up with a Toomre Q parameter ≃1.5, a condition that ensures that it is unstable to bar formation (Athanassoula & Sellwood 1986). The initial conditions were evolved using a three-dimensional, cylindrical, polar-grid-based N-body code described in Sellwood & Valluri (1997) and Sellwood (2014a).

The disk formed a bar that subsequently experienced vertical thickening via the buckling instability (Toomre 1966, p. 111; Raha et al. 1991; Sellwood & Wilkinson 1993) and developed the peanut-shaped bulge that is characteristic of this instability. After the disk reached a nearly steady state (simulation time t₁ = 700), a central point mass⁵ representing an SMBH of mass 0.2% ${M}_{d}$ was grown adiabatically (for details see Brown et al. 2013). At t₂ = 1200 the transients due to the changing central point mass (SMBH) potential had dissipated and the simulation was frozen, and orbits were examined and compared with orbits from the frozen potential at t₁.

Hereafter we refer to the frozen snapshots at t₁ and t₂ as Model A and Model B, respectively. The bar pattern speed at t₁ and t₂ was computed from several successive time steps before and after each snapshot. The galaxy potential for each snapshot is derived from the full N-body distribution and is not described by an analytic form. Since our objective is to understand the actual orbits that populate the N-body distribution function, we randomly select 10,000 particles (from the total of ∼2.8 million) uniformly sampling the entire distribution function. The instantaneous positions and velocities of the selected particles (transformed to the rotating frame of the bar) were used to advance each of the 10,000 particles individually while keeping the rest of the particles fixed, by using a modification of the N-body code used to run the simulations. Orbit integration was carried out as described in footnote 10. The orbits were integrated in Cartesian coordinates⁶ in the rotating frame using the accelerations derived from the frozen potential of that snapshot and adding the appropriate Coriolis and centrifugal pseudoforces determined by the bar pattern speed (see Equations (2)–(4)). Each orbit was integrated for 1000 time units (corresponding to $\sim 2.1\times {10}^{9}$ yr for the units adopted above) and sampled at 20,000 equally spaced time intervals. Since orbits at different radii have different orbital periods, this corresponds to hundreds of orbital periods for the innermost particles but to just tens of orbital periods for orbits in the outskirts of the disk.

Since the 10,000 particles for which orbits were integrated were selected at random from the entire distribution function, disk particles were included. The bar length in N-body simulations is generally estimated to be the radius at which the strength of the bar mode A₂ drops to some fraction (e.g., 1/2) of its maximum value (Athanassoula & Misiriotis 2002). By this criterion the bar length in this simulation is estimated to be about 3 (in program length units). However, visual classification of all 10,000 orbits in Model A showed that a significant fraction of bar-like orbits continue to exist all the way out to ∼4 units in length (see the Appendix).

We compare orbits from the bar snapshots with orbits in stationary and rotating triaxial potentials that are generalizations of the spherical "Dehnen models" (Dehnen 1993; Tremaine et al. 1994). These models have density profiles that provide good fits to luminosity profiles of elliptical galaxies and the bulges of spiral galaxies. The density distribution in these models is stratified on concentric ellipsoids with principle axes aligned with Cartesian coordinates x, y, z and with the semimajor, semi-intermediate, and semiminor axis lengths a, b, c, respectively. The parameter γ determines the logarithmic slope of the central density cusp and ranges observationally from γ = 0.1−1 in bright galaxies with shallow cusps to γ = 2 in fainter galaxies with steep cusps (Gebhardt et al. 1996; Lauer et al. 2007). We restrict ourselves to models with γ = 0.1 since more cuspy models produce many resonant boxlet orbit families (Miralda-Escude & Schwarzschild 1989; Merritt & Valluri 1999) that are not found in our N-body bars. Equations describing the potential and the accelerations for this model are given in Merritt & Fridman (1996) and Deibel et al. (2011).

In a rotating frame, the Jacobi integral (E_J) is a conserved quantity (equivalent to energy in a stationary frame):

$$\begin{eqnarray}&&{E}_{{\rm{J}}}=\displaystyle \frac{1}{2}| \dot{{\boldsymbol{x}}}{| }^{2}+{\rm{\Phi }}-\displaystyle \frac{1}{2}| {{\rm{\Omega }}}_{p}\times {\boldsymbol{x}}{| }^{2},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{x}}$ and $\dot{{\boldsymbol{x}}}$ are three-dimensional spatial and velocity vectors, respectively. When the rotation is about the short axis (z), the equations of motion in Cartesian coordinates are given by

$$\begin{eqnarray}&&\ddot{x}=-\displaystyle \frac{\partial {\rm{\Phi }}}{\partial x}+2{{\rm{\Omega }}}_{p}\dot{y}+{{\rm{\Omega }}}_{p}^{2}x,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\ddot{y}=-\displaystyle \frac{\partial {\rm{\Phi }}}{\partial y}-2{{\rm{\Omega }}}_{p}\dot{x}+{{\rm{\Omega }}}_{p}^{2}y,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\ddot{z}=-\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z},\end{eqnarray} \tag{ 4 }$$

where $2{{\rm{\Omega }}}_{p}\dot{y}$ and $-2{{\rm{\Omega }}}_{p}\dot{x}$ are components of the Coriolis force and ${{\rm{\Omega }}}_{p}^{2}x$ and ${{\rm{\Omega }}}_{p}^{2}y$ are components of the centrifugal force. Following standard practice, we adopt the right-handed rule where motion anticlockwise about the z-axis (as viewed from positive z) has positive angular momentum, while clockwise motion about the z-axis has negative angular momentum. In this terminology the direction of the pattern rotation of the simulated bars and the direction in which the triaxial Dehnen model is rotating when viewed from an inertial frame are both positive (anticlockwise). Similarly, when we discuss short-axis (long-axis) tubes, anticlockwise motion about the z-axis (x-axis) corresponds to positive angular momentum J_z (J_x).

Deibel et al. (2011) studied orbits in triaxial Dehnen models with a variety of shapes subjected to a range of pattern speeds. Here we restrict our comparison to orbits in a prolate triaxial potential with axis ratios $c/a=0.4$, $b/a=0.48$, and hence triaxiality parameter $T=({a}^{2}-{b}^{2})/({a}^{2}-{c}^{2})=0.916$ and a weak cusp (γ = 0.1). This shape is quite close to that of the bars in our N-body simulations. As in the case of the N-body simulation, we adopt a set of units where the total mass of the model M, the semimajor axis scale length a, and the gravitational constant G are set to unity.

Following standard practice, the patten speed used to describe the tumbling of the triaxial figure is defined in terms of the "co-rotation radius," hereafter R_CR. In a nearly axisymmetric potential the co-rotation radius is the radius at which the azimuthal frequency of a closed (almost circular) orbit in the equatorial (x–y) plane of the potential is the same as the pattern frequency (generally called "pattern speed") Ω_p. The pattern speeds of bars have been measured by applying the Tremaine–Weinberg method (Tremaine & Weinberg 1984; Meidt et al. 2008) to stellar kinematical velocity fields (as well as Hα and H i). These measurements find that the ratio of the co-rotation radius to the length of the bar is R_CR/R_bar = [1, 1.4] for a bar of length ${R}_{{\rm{bar}}}$ (e.g., Merrifield & Kuijken 1995; Debattista et al. 2002; Debattista & Williams 2004; Aguerri et al. 2003; Corsini 2011; Aguerri et al. 2015). For the bar in the simulations described in the previous section, ${R}_{\mathrm{CR}}/{R}_{{\rm{bar}}}\sim 1$. For the orbits in the triaxial Dehnen model the pattern speed was chosen so that the co-rotation radius was at roughly the same radial distance from the center as it is in the case of the bar simulations (i.e., at ∼3 units or one bar length). Note that the triaxial ellipsoidal surface with major-axis length ∼3 in the Dehnen model encloses only 1/2 of the total mass of the model. To ensure a fair comparison, the orbits from the stationary and rotating triaxial Dehnen model were selected to have similar radial extents to the bar orbits that they are compared with.

Orbits from the N-body bars were classified both visually and using our new automated orbit classification algorithm. All 10,000 orbits in Models A and B were visually classified by C.A. In addition, the innermost ∼4000 particles of Model A were visually classified by M.V. Visual classification was based on x–y, x–z, y–z projections and plots of angular momenta as a function of time (${J}_{x}(t),{J}_{z}(t)$) (examples of such plots are given in Figures 1, 4, 5, 6). We adopted classifications based on the full orbit (integration time t = 1000), which is ∼200 orbital periods for the innermost bar orbits but only ∼30 orbital periods for the outermost disk particles. Although all orbits in an N-body system are mildly chaotic (e.g., Miller 1964; Hemsendorf & Merritt 2002), for the purpose of testing our automated classification scheme, we also attempt to visually classify chaotic orbits, although this is not possible to do in a robust manner. An orbit was visually classified as chaotic only if it could not be easily identified with a major orbit family (box, short-axis tube, long-axis tube), or if it showed signs of changing from one family to another during the integration time (indicating that it is sticky chaotic). For the ∼4000 orbits that were visually classified by two of us (C.A. and M.V.) we found that fewer than 3% of orbits were classified differently. Where classifications disagreed, the orbits were generally transitional orbits—probably lying close to the separatrix between two families.

Zoom In Zoom Out Reset image size

In Model A visual classification resulted in fewer chaotic orbits than were found by the automatic classification method (which relies on orbital frequency drift; see Appendix A); however, both methods yielded similar numbers of chaotic orbits in Model B. It is difficult to visually identify weakly chaotic (sticky chaotic) orbits since these typically lie at the edge of a resonant island or between two resonances and can appear regular for long times.

In Section 4.2 we show frequency maps with several minor resonant families. In frequency maps resonant orbits appear clustered along thin lines that satisfy a resonance condition $l{{\rm{\Omega }}}_{x}+m{{\rm{\Omega }}}_{y}+n{{\rm{\Omega }}}_{z}=0$ (where l, m, n are small integers). The signs of the integers indicate the phase relationship between the frequency components. The automatic classification method easily identifies resonant orbits, but these are visual identified only if they are quite close to the resonant parent orbit. Weakly resonant orbits are visually classified as boxes.

Voglis et al. (2007) used the orbital fundamental frequencies in cylindrical polar coordinates (${{\rm{\Omega }}}_{R},{{\rm{\Omega }}}_{\phi },{{\rm{\Omega }}}_{z}$) to classify orbits in an N-body bar. In Appendix B we show that a much clearer separation of orbit families is obtained with fundamental frequencies computed in Cartesian coordinates. Our new automatic bar orbit classification scheme is described in detail in Appendix B, and a detailed comparison between the automatic and visual classification schemes is deferred to Section 4.1.

The family of orbits parented by the periodic x1 orbit is considered to be the main bar-supporting orbit family (see Binney & Tremaine 2008, hereafter BT08). As mentioned previously, a particle on such an orbit makes two excursions in radius during each complete circuit in the azimuthal angle ϕ (e.g., Contopoulos & Papayannopoulos 1980), and therefore the radial oscillation frequency and tangential frequency are resonant ${{\rm{\Omega }}}_{R}$:${{\rm{\Omega }}}_{\phi }$ = 2:1. Since its apocenter radius increases with increasing Jacobi integral E_J, at large E_J, x1 orbits develop loops at their extremities. A particle on an x1 orbit travels in a prograde sense about the center of the galaxy (i.e., in the same sense that the bar pattern is rotating) except in the loops where the motion is retrograde relative to the figure.

In triaxial ellipsoids the main orbit family responsible for providing the high density along the long axis is the three-dimensional box orbit family (BT08). In a stationary potential boxes have no net angular momentum about any axis. However, in a rotating frame Coriolis forces result in "envelope doubling," which imparts a small net angular momentum to the parent x-axis orbit, as well as box orbits (Schwarzschild 1982; de Zeeuw & Merritt 1983). The "x1 orbit" is therefore also the parent of the quasi-periodic box orbit family in a rotating triaxial potential (Martinet & de Zeeuw 1988). In this section we select orbits from our self-consistent N-body bar simulations to illustrate that the vast majority of orbits in these simulations are boxes with little net rotation in the bar frame, unlike their prograde parent x1 orbit.

Figure 1 shows six different orbits from the pure N-body bar (Model A) computed in the frame of reference co-rotating with the bar. On each row we show two projections (x–y and x–z) and the angular momentum as a function of time ${J}_{z}(t)$, for a single orbit. ${J}_{z}(t)$ is normalized relative to its maximum absolute value over the entire orbit. The top two rows show orbits near the closed periodic x1 parent (it is difficult to find strictly periodic x1 orbits in an N-body simulation). For these two orbits we see that the angular momentum J_z oscillates between two positive values throughout its orbit. The next three rows show orbits that travel increasingly farther from the parent x1 orbit. For these orbits ${J}_{z}(t)$ does not remain positive but becomes negative for those portions of the orbit when the motion is retrograde. The examples in Figure 1 demonstrate that as orbits deviate farther from the periodic x1 parent (i.e., as their extent in the y direction increases) they spend more and more time moving retrograde. In fact, such orbits are found to be the primary building blocks of the central parts of orbit-based bar models studied by Patsis & Katsanikas (2014b).

One way to quantify the deviation from the parent x1 orbit is to measure the extent of an orbit in the y direction relative to its extent in the x direction. We do this by computing the ratio of the maximum absolute y value attained over the orbit, $| y{| }_{{\rm{max}}}$, relative to the maximum absolute x value attained over the orbit, $| x{| }_{{\rm{max}}}$. Based on the visual classification of orbits, we find that classical x1 orbits have $| y{| }_{{\rm{max}}}/| x{| }_{{\rm{max}}}\lt 0.35$ (they are significantly more elongated along the x-axis than along the y-axis).

Figure 2 shows $\langle {J}_{z}\rangle$ (the time average of the angular momentum J_z normalized to its maximum value) as a function of $| y{| }_{{\rm{max}}}/| x{| }_{{\rm{max}}}$. Each dot represents an orbit in Model A that was classified as x1 or box (we exclude the resonant boxlet orbits that are discussed in the next section). The solid line is the mean value of the distribution represented by dots. Orbits with the smallest values of $| y{| }_{{\rm{max}}}/| x{| }_{{\rm{max}}}$ have the highest net angular momentum and are classical x1 orbits. As orbits get thicker in y, this average angular momentum decreases. For reference, an orbit with constant angular momentum would have $\langle {J}_{z}\rangle =\pm 1;$ the orbits in the top two rows of Figure 1 have $\langle {J}_{z}\rangle =0.6$, while the orbit in the bottom row of Figure 1 has $\langle {J}_{z}\rangle =-0.09$. Only a small fraction of orbits have both the high $\langle {J}_{z}\rangle$ and the small $| y{| }_{{\rm{max}}}/| x{| }_{{\rm{max}}}$ that are characteristic of classical x1 orbits. In fact, from Figure 2 we see that most orbits have $\langle {J}_{z}\rangle \lt 0.25$ and a significant fraction have negative $\langle {J}_{z}\rangle$, implying that they spend more time traveling retrograde than prograde (for example, the orbit in the last row of Figure 1).

Zoom In Zoom Out Reset image size

We illustrate this point further by launching orbits from box orbit initial conditions in a prolate triaxial ellipsoid with the Dehnen density profile (Section 2.2), with shape parameters close to that of the N-body bar in Model A ($c/a=0.4,b/a=0.48,T=0.916,\gamma =0.1$). In Figure 3 each row shows two orbits launched from the same initial conditions in a stationary Dehnen model (three left columns) and in the Dehnen model with a fast pattern speed (${R}_{\mathrm{CR}}/{r}_{1/2}=1.1$, where ${r}_{1/2}$ is the half-mass radius of the model; three right columns). For each orbit we show projections of the orbit (x–y, x–z) and its angular momentum with time (${J}_{z}(t)$). The orbits in the top row were launched very close to the long (x) axis of the model. On the left is a long-axis orbit (parent of the box family) with no net angular momentum J_z or J_x. The three right-hand columns show what happens to this orbit in a rotating frame: the Coriolis force now causes the orbit to loop around the center in an anticlockwise sense, and hence ${J}_{z}(t)$ becomes strictly positive; i.e., in the rotating frame this orbit has become an "x1 orbit." In the second row the three left columns show a standard box orbit in a stationary potential while the three right columns show the same orbit in the rotating frame. This quasi-periodic orbit behaves in a manner identical to the "thick x1" orbits from Model A (e.g., in rows 4 and 5 of Figure 1). (The box orbits in the third and fourth rows of Figure 3 will be discussed later.)

Zoom In Zoom Out Reset image size

The differences seen in the shapes and angular momentum distributions of the orbits in the stationary potential (three left-hand columns) and in the rapidly rotating potential (three right-hand columns) are entirely a consequence of the pseudoforces in the rotating frame. A comparison of orbits in the first two rows of this figure with the x1 family orbits drawn from the N-body bar in Figure 1 shows that both box orbits and x1 orbits are parented by the linear long-axis orbit. Figure 10 of SS04 gives very clear examples of quasi-periodic orbits parented by an x1 orbit, as well as orbits that they refer to as "fat x1" (but that we call boxes). Despite their visual difference (classical x1 versus box-like orbits), SS04 show that these orbits form a continuous sequence in an SoS (e.g., upper panel of Figure 9 in SS04). We discuss this further in Section 4.2.1.

To summarize, it is well known that the x1 orbit in bars is the same as the parent of the box family in rotating triaxial potentials (Schwarzschild 1982; Martinet & de Zeeuw 1988). What is less well appreciated is that in bars, as in triaxial ellipsoids, the dominant bar-supporting family is box-like with little angular momentum about the short axis, unlike their periodic x1 parent, which rotates prograde.

3.1.1. Boxlets: Fish, Pretzels, and Banana Orbits

Miralda-Escude & Schwarzschild (1989) and Merritt & Valluri (1999) showed that when an integrable triaxial potential is perturbed, e.g., by the introduction of a cusp or a central point mass, almost all orbits that originate from "box-like" initial conditions (i.e., launched with zero initial velocity from an equipotential surface) become either resonant or chaotic. Resonant orbits are easily identified by frequency mapping. In stationary triaxial potentials they include the well-known 3:0:−2 "fish" resonance and the 4:−3:0 "pretzel" resonance—both of which are absent in our rapidly rotating Dehnen model and N-body bars. Commonly found resonant orbits in the rotating Dehnen model are shown in the third row of the three right columns of Figure 3 (${{\rm{\Omega }}}_{x}$:${{\rm{\Omega }}}_{y}$:${{\rm{\Omega }}}_{z}$ = 3:−2:0, fish/pretzel) and in the first and second rows of the three right columns of Figure 3 (the "banana" resonance). A few other resonances are seen in the frequency maps in Figure 12.

Our examination of orbits from the N-body bar showed that they can also be associated with resonances. Figure 4 shows examples of resonant orbits found in the N-body bar. Orbits in the first three rows are all characterized by ${{\rm{\Omega }}}_{x}$:${{\rm{\Omega }}}_{y}$:${{\rm{\Omega }}}_{z}$ = 3:−2:0 and look like fish/pretzels (and correspond to the same resonance as in the third row of Figure 3). The fourth row shows an orbit associated with "x1-banana" resonance (${{\rm{\Omega }}}_{x}$:${{\rm{\Omega }}}_{y}$:${{\rm{\Omega }}}_{z}$ = 2:−2:−1),⁷ and the fifth row is an "x1-anti-banana" (Miralda-Escude & Schwarzschild 1989; Pfenniger & Friedli 1991), which satisfies the same resonance conditions as a banana orbit, but passes through $x=0,z=0$ (the automatic classification code is unable to distinguish between banana and anti-banana, which are treated as one group). As pointed out previously by Pfenniger & Friedli (1991), the banana and anti-banana orbits are vertical bifurcations of the x1 family (as can be seen from their x–y projections in the left-hand column). These orbits are referred to as the x1v1 family by Patsis et al. (2002). Orbits like those in the second row of Figure 3 are boxy-bananas: they have a banana shape in x–z associated with 2:0:−1 resonance. The third and fourth columns of Figure 4 show that like box orbits in triaxial potentials, the sign of J_z changes at each extremum, implying that the orbits have little or no net angular momentum about the z-axis, despite being parented by x1 orbits. The frequency maps in Figure 13 show several other resonances. The presence of resonant boxlet orbits was noted by SS04 (see their Figure 13), especially at high E_J values.

Zoom In Zoom Out Reset image size

Vertical bifurcations of the x1 orbit associated with the 2:−2:−1 resonance may be responsible for the "X-shaped" structures seen in buckled edge-on bars (Patsis et al. 2002; Athanassoula 2005). Patsis & Katsanikas (2014a, 2014b) also present a dynamical mechanism for building X-shaped peanuts with families of periodic orbits that are not bifurcations of x1 orbits. It was suggested by Portail et al. (2015) that a resonant family of ${{\rm{\Omega }}}_{x}$:${{\rm{\Omega }}}_{y}$:${{\rm{\Omega }}}_{z}$ = 3:0:−5 "brezel" orbits (bottom row of Figure 4) are the backbone of the "X-shaped" structures in their made-to-measure N-body bar models, but we found only 88 "brezels" (1.5% of the bar orbits) in both Models A and B. A more detailed analysis of the orbits contributing to the X-shape will be discussed in a future work (C. Abbott et al. 2016, in preparation).

De Zeeuw (1985a) showed that for integrable triaxial potentials two types of long-axis tubes exist, the inner long-axis tubes and the outer long-axis tubes. This family is "parented" by closed 1:1 periodic orbits that lie in the y–z plane with angular momentum along the x-axis. Heisler et al. (1982) studied the stability of these periodic orbits and found that they are stable to figure rotation, but the Coriolis force tips them about the y-axis in a direction that depends on the sign of J_x. Two such stable periodic orbits exist rotating clockwise and anticlockwise about the x-axis, the orbits with positive J_x are tipped clockwise about the y-axis, while orbits with negative J_x are tipped anticlockwise about the y-axis. These orbits were termed "anomalous" by van Albada et al. (1982) since they both acquire retrograde motion about the z-axis in the rotating frame as a result of being tipped. Since the long-axis tube family is "parented" by anomalous orbits, they too are expected to be stable and "tipped" about the y-axis in a rotating frame, and may acquire some retrograde angular momentum about the z-axis.

Figure 5 shows examples of inner long-axis tubes (top two rows) and outer long-axis tubes (third and fourth rows) selected from bar Model A. The second column (y–z projection) shows that these orbits circulate about the long (x) axis. The rightmost column shows the angular momentum J_x about the x-axis, which, as expected, is strictly positive or strictly negative. An inspection of the first column shows that the orbits are tipped about the y-axis exactly as predicted by Heisler et al. (1982), further supporting our claim that bars and rotating triaxial ellipsoids have fundamentally similar orbital building blocks. As far as we are aware, previous studies of self-consistent orbits in three-dimensional N-body bars (Voglis et al. 2007) have not found long-axis tubes, suggesting that their existence could be sensitive to bar shape and the details of its formation mechanism.

Zoom In Zoom Out Reset image size

Short-axis (z) tubes constitute an important family in oblate-triaxial ellipsoids but are somewhat less important in more prolate systems. Short-axis tubes are parented by closed periodic 1:1 orbits that lie in the x–y plane and circulate about the z-axis. These orbits are known to be stable to figure rotation (de Zeeuw & Merritt 1983). However, Deibel et al. (2011) found that as the pattern speed of a triaxial figure was increased, prograde short-axis tubes were replaced by retrograde ones. This was also found by Martinet & de Zeeuw (1988) for loop orbits in the x–y plane of a rotating triaxial ellipsoid and occurs because tube/loop orbits launched prograde "fall behind" the figure, becoming retrograde as the pattern speed increases.

Periodic x2 and x4 orbits (like the classical x1 orbit) satisfy the condition ${{\rm{\Omega }}}_{R}$:${{\rm{\Omega }}}_{\phi }$ = 2:1. In the characteristic diagrams (e.g., Sellwood & Wilkinson 1993) the prograde x2 family is found at low E_J and is elongated along the y-axis of the bar, while the x4 family is retrograde and is found over a wide entire range of E_J values. Both the x2 and x4 orbits are elongated along the y-axis at small E_J. The x4 family is elongated along the y-axis at small E_J but becomes rounder at large E_J. Since x2 and x4 orbits are planar periodic orbits, they are not expected to be found in significant numbers in an N-body model, but they do parent families of prograde and retrograde short-axis tubes, which can be easily identified if they exist. Henceforth, we refer to all nonplanar orbits with a fixed sign of angular momentum about the short axis as "z-tubes."

The visual classification of 20,000 orbits in Model A and Model B identified ∼1.5% of bar orbits in Model A as retrograde z-tubes, but none of the orbits in this model were found to be prograde z-tubes. The top two rows of Figure 6 show examples of retrograde z-tubes from Model A. The fact that we were unable to find a single example of a prograde z-tube orbit is consistent with the findings of Sparke & Sellwood (1987) and Voglis et al. (2007), who only found retrograde x4 orbits in their N-body bar models. x2 orbits (and orbits parented by them) are expected from the locations of the inner Lindblad resonance and the co-rotation resonance in Model A; however, this family appears to be severely underpopulated. A significant decline in prograde z-tubes at high pattern speeds is, however, predicted by Martinet & de Zeeuw (1988) and Deibel et al. (2011).

Zoom In Zoom Out Reset image size

We examined whether it was possible to generate x2 and x4 orbits from perturbations to a linear y-axis orbit in the same way as it was possible to generate the x1 orbit from the linear x-axis orbit. The y-axis linear orbit is known to be unstable (Adams et al. 2007); hence, launching orbits from the y-axis of the Dehnen model yielded unstable chaotic orbits. An example of an orbit launched from the "stationary start space" (i.e., with no initial velocity) near (but not on) the y-axis is shown in the fourth row of Figure 3. In the stationary frame (three left columns) the orbit is a stable box elongated along the y-axis. In the rotating frame (three right columns) this orbit becomes an even thicker box but now with $\langle {J}_{z}\rangle \lt 0$, similar to the orbit in the last row of Figure 1. If $| {J}_{z}|$ was smaller than a minimum value, box-like orbits elongated along the y-axis always resulted. It was only possible to generate retrograde x4 orbits by launching an orbit from near the y-axis with a substantial $| {J}_{z}|$. All the retrograde z-tubes in Model A are found at small values of E_J and are more elongated along the y-axis than along the x-axis.

In Model B, however, we do find prograde and retrograde z-tube orbits. Some are more elongated along the y-axis, while others are more elongated along the x-axis than along the y-axis. We defer a discussion of this to Section 4.1, where we argue, based on the work of Brown et al. (2013), that this new population of prograde z-tubes is a result of the growth of the SMBH in the bar, which induces angular momentum redistribution.

Figure 6 shows three prograde z-tubes from Model B. In the third row the orbit is elongated along the y-axis (i.e., it is parented by the periodic x2 orbit); in the fifth row it is elongated along the x-axis, and the fourth row shows a transitional orbit. We assert that all these orbits belong to the same family since they must have significant initial net angular momentum to avoid becoming boxes, in contrast with x1 orbits, which are related to boxes and derive their prograde motion from the Coriolis forces.

Previous studies of orbits in triaxial and axisymmetric potentials have shown that as one approaches the massive central point mass in the region of the potential where the mass of stars is comparable to the mass of the SMBH, orbits begin to resemble Keplerian ellipses, which are perturbed by the large-scale triaxial or axisymmetric stellar potential and therefore precess slowly (Sridhar & Touma 1999; Poon & Merritt 2001; Merritt & Vasiliev 2011; Li et al. 2015). Orbits associated with black holes include "saucers," "pyramids," and a variety of resonant families (Merritt & Valluri 1999; Merritt & Vasiliev 2011). In this work we will refer to these collectively as PKOs.

The spatial resolution of our simulations was high enough⁸ that random sampling of the distribution function of Model B yielded a small number of PKOs in the vicinity of the black hole. Figure 7 shows examples of PKOs that were found in Model B. Each row shows a single orbit: the three panels on the left show the orbit integrated for about 20 orbital periods,⁹ while the three panels on the right show the projected density distribution of the orbit integrated over hundreds of orbital periods (for t = 1000 units). Notice that in the short integrations (left three columns) the orbits in the top three rows look very similar to each other: they are all slowly precessing ellipses. However, the projected density distributions of the same orbits show that they are morphologically quite different from each other. The orbit in the top row is a box (although ${J}_{z}(t)$ and ${J}_{x}(t)$ are not shown, it also has no net angular momentum about either axis), the second row shows a short-axis (z) tube for which ${J}_{z}(t)\gt 0$, and the third row shows a resonant short-axis tube. The bottom row shows a long-axis tube (${J}_{x}(t)\lt 0$) that is tipped anticlockwise about the y-axis. (The orbits in the last two rows appear to pass through the origin, but zoomed-in plots show that they do not.) In total, 113 orbits (1.8% of bar orbits) in Model B were found to be PKOs. At the very small radii at which these orbits are found, the effects of the centrifugal forces from the rotating pattern should be quite small, but the Coriolis forces are large enough (because the velocities are high) to produce effects similar to those seen on orbits at larger radii. A detailed study of the effects of figure rotation on these nearly Keplerian orbits is beyond the scope of this paper and will be investigated in a future paper.

Zoom In Zoom Out Reset image size

Table 1 compares the overall orbit fractions obtained by both classification methods (visual and automatic) for both models. The integer in each column represents the total number of orbits in that family (from a total of 10,000 orbits in each model), while the percentage of bar orbits (i.e., excluding disk orbits) contributed by this family is given in parentheses. Note that for Model B we exclude the 113 visually classified "PKO" orbits since the automated method is unable to distinguish these from boxes, long-axis tubes, and short-axis tubes at large radii.

In Figures 8, 10, and 11, and Table 1 we present x1 orbits, 3:−2:0 orbits, and boxes separately but remind readers that they are all members of the box orbit family.

Zoom In Zoom Out Reset image size

The biggest difference between the automatic and visual classification is the fraction of orbits that are classified as "chaotic" or "x1/x1+banana." Many of the orbits that are visually classified as "x1/x1+banana" are automatically classified as chaotic (and, as we see in Figure 8, tend to lie toward the outer part of the bar). In addition, while 90 orbits were visually classified as retrograde z-tube orbits, only 31 were automatically classified as such. Examination of the remaining 60 showed that they were also identified as chaotic by the automatic classifier.

In Figure 8 we compare how the orbit populations in Models A and B change with R_apo (the apocenter radius of orbits in the x–y plane) and how the orbit population at each radius depends on the classification scheme.

We distribute bar orbits into 6 bins in ${R}_{\mathrm{apo}}\lt 4\;{\rm{kpc}}$. Figure 8 shows the percentage of orbits (on a logarithmic scale) in each bin as a function of R_apo in each of the main families (as indicated by the legends). The left-hand column shows percentage of orbits in each radial bin as determined by the automated classification scheme, while the right-hand column shows the same plot using visual classification. At all radii the orbits classified as boxes (blue triangles, connected by blue dot-dashed lines) dominate the population in both models. The fractions of boxes obtained by both methods are also nearly identical. The x-tube orbits (pink circles connected by pink lines) decline with radius in both models, and there are slightly more of them in Model A than in Model B (both classification methods yield similar fractions of these orbits). The short-axis tube family shown by green stars connected by solid lines is slightly more important in Model B than in Model A. The trend of the short-axis tube population with radius is similar in both classification schemes. There are somewhat more chaotic orbits (black circles connected by black lines) in the inner regions of Model B than in Model A (with both methods). Visual classification yields fewer chaotic orbits in the outer regions of both models than does the automatic classification. The fractions of x1+banana orbits (red triangles connected with solid red lines) and 3:−2:0 resonant orbits (yellow squares connected with solid yellow lines) differ between automatic (left) and visual classification methods (right), but both populations are small (less than 10%) at most radii. The biggest differences are in the outer part of the bar, where the visual classification identifies orbits as x1 or 3:−2:0, while the automated classification classifies these orbits as chaotic.

Overall Figure 8 shows that the populations in Models A and B resemble each other and that by both methods of classification box orbits and chaotic orbits make the most significant contribution to the overall population within the bar in both models. All other families (x-tubes, z-tubes, x1+banana, and resonant 3:−2:0) each contribute less than 10% at any radius. The overall trends with radius are similar, but there are differences of a few (∼2%–3%) between the two models and between the numbers obtained by visual and automated classification schemes. These differences are comparable to the differences in classification obtained by two different visual classifiers. Since these differences are small and unavoidable, in the rest of this paper we use automated classification as the basis for our analysis.

4.2.1. Poincáre SoSs

The distribution of E_J values for all the orbits in each of the two models is shown in Figure 9. Since the particles were selected at random from the simulations, each distribution is representative of the distributions of E_J values for the entire model. Five bins containing 100 orbits each were defined to lie roughly at equal intervals in E_J between the minimum value of E_J and the value at the co-rotation radius of the bar. The bin limits are shown as gray bands overlying the histograms. These bins were used to select orbits to plot the SoSs that follow.

Zoom In Zoom Out Reset image size

A standard way to represent the phase-space distribution of orbits in a bar is to plot a Poincáre SoS for a number of orbits at a single value of the Jacobi Integral E_J. For two-dimensional bars it is customary to construct an SoS by plotting the velocity component V_y as a function of y every time an orbit intersects the x-axis with a negative value of V_x (e.g., Sellwood & Wilkinson 1993). In three-dimensional bars one might choose to restrict orbits plotted on an SoS to those that are confined to the equatorial (i.e., x–y) plane of the model (e.g., Shen & Sellwood 2004). When orbits of a single value of E_J are plotted on an SoS, regular orbits follow thin closed or broken curves, resonant orbits form groups of islands, while chaotic orbits tend to fill larger areas.

Two problems arise in plotting an SoS for orbits from an N-body simulation: first, all orbits occupy three spatial dimensions and are not well represented on a two-dimensional SoS; second, since orbits have a range of E_J values (see Figure 9), curves in the SoS appear fuzzy. Figure 10 show SoSs centered on five different E_J values in Model A. The left-hand column shows y versus V_y when orbits cross the y–z plane with positive V_x; the SoS is appropriate for x1 orbits (red dots), box orbits (blue dots), 3:−2:0 orbits (orange dots), and x4 orbits (green dots). The SoSs in the right-hand column show V_z versus z when orbits cross the x–z plane with negative V_y; these coordinates are more meaningful for x-axis tubes (pink dots). Chaotic orbits (black dots) were plotted in the right-hand plots simply to avoid crowding.

Zoom In Zoom Out Reset image size

In Figure 10 the spread in E_J (as indicated in the legends) results in few clear curves in the SoSs. Nonetheless, one can see that each orbit family occupies a slightly different region of the SoS. x1 orbits (red dots) that are prograde about the z-axis cluster around V_y = 0 at positive y values (see dense red core at the bull's-eye of the SoS in Figure 10). Box orbits (blue regions) lie on oval curves surrounding the x1 orbits, indicating that they belong to the same sequence. (This is much more clearly seen in Figure 9 of SS04, which shows that the x1 orbits form the smallest ovals forming a "bull's-eye," and this sequence of ovals extends to larger radii.) Since the bigger circles (and the blue regions in our plots) can have negative y values, it implies that their J_z becomes negative as they acquire a box-like appearance. In an SoS resonant orbits appear as multiple islands. Seen from the point of view of dynamical tools like the SoS, box orbits and x1 orbits do belong to the same family, although they can have slightly different appearances. SS04 show that a large number of resonances of various orders can (in principle) exist in an N-body bar potential. However, our SoSs (and frequency maps in the next section) show that only one prominent resonant family is populated in the distribution function. The 3:−2:0 family appears in the SoS in Figure 10 as orange points that appear in multiple discontinuous islands in the region occupied by box orbits. In Figure 10 retrograde z-tubes (green dots) lie to the left of y = 0, and there are no prograde (x2) orbits in this model.

Figure 11 shows SoSs after the growth of the SMBH (Model B). The most striking difference between Model B (after SMBH growth) and Model A is the appearance of prograde z-tube orbits (parented by the x2 orbit). It is worthwhile reviewing briefly the results of Brown et al. (2013), which will shed light on the appearance of these orbits. Brown et al. (2013) carried out a detailed analysis of the intrinsic velocity distributions (velocity dispersion profiles and velocity anisotropy) in the models analyzed here, as well as an axisymmetric disk in which an SMBH was grown in an identical manner. In both a barred and axisymmetric galaxy, the growth of the SMBH increases the depth of the central potential ("adiabatic contraction"), but the density increase was significantly higher in the barred galaxy than in the axisymmetric one. Furthermore, as matter is dragged inward in the axisymmetric galaxy, the distribution function becomes tangentially biased and inflow is limited by angular momentum conservation. However, in a time-dependent barred galaxy the inward flow of mass is enhanced by outward transport of angular momentum. This results in a higher central stellar density in the bar than in the axisymmetric model. The differences between Model A and Model B arise as a result of two effects: (a) elongated orbits, e.g., x1 orbits and x-tubes are unable to survive the growth of the central point mass, which makes the potential more axisymmetric (Deibel et al. 2011); and (b) inflow of matter from large radii with positive angular momentum results in the formation of prograde z-tube orbits seen in the SoS, which did not exist in the original bar.

Zoom In Zoom Out Reset image size

4.2.2. Frequency Maps

An alternative way to represent the phase-space distribution function is via a "frequency map." Frequency mapping exploits the fact that most orbits in galaxies are approximately quasi-periodic (BT08); hence, a Fourier transform of their space and/or velocity coordinates can be used to obtain the fundamental orbital frequencies that characterize each regular orbit (for a detailed description of the orbital spectral analysis method see Appendix A). A frequency map is useful for representing the phase-space structure of a large sample of orbits even if they have a wide range of energies (or E_J). On such a map one can plot a large number of orbits that are representative of an entire distribution function (instead of just a subset at discrete values of E_J as in the case of SoSs). In an appropriate coordinate system, orbits belonging to different families will appear in different regions of the map, providing an easy way to visually assess the importance of various orbit families to the distribution function and the range of E_J values for which each family is important. Frequency maps are also useful for identifying different resonances and their importance to the distribution function (Valluri et al. 2010, 2012).

A frequency map is obtained by plotting pairs of fundamental orbital frequencies or ratios of such frequencies against each other. We integrated a large number of orbits in the Dehnen model at a single energy level launched from "stationary start space" initial conditions (i.e., zero initial velocity) to generate box orbits. The same initial conditions were integrated with figure rotation. Figure 12 (top) shows the three-dimensional distribution of these 3888 initial conditions, with each point given a unique color determined by its coordinates (x, y, z). The RGB color indices are determined such that points that lie on the x, y, z axes are red, green, and blue, respectively, and the colors of other points reflect their distances from the principal axes.¹⁰ The middle panel shows the frequency map in the stationary Dehnen model, while the bottom panel shows the frequency map of the same orbits subjected to figure rotation. The colors of points in the frequency maps enable the reader to visually map points on the frequency maps to their initial launch positions in the top panel.

Zoom In Zoom Out Reset image size

In the stationary Dehnen model (middle panel) most of the points lie in a fairly regular grid above the diagonal (1, −1, 0) resonance line and below the horizontal (0, 1, −1) resonance line. In the rotating Dehnen model (bottom panel) most of the points move to the left, but a few points appear along the (1, −1, 0) diagonal near the label "x1." Their red color indicates that they were launched from near the x-axis, and as expected they are converted to "x1" orbits in the rotating frame. There are box-like red points associated with the (2, 0, −1) "banana resonance," which is bifurcation of the x1 family (these are the box-like orbits that also have the banana shape in x–z projection, e.g., second row of Figure 3, but are launched from some height above the x–y plane).

As predicted by Martinet & de Zeeuw (1988), z-axis orbits (and orbits launched close to the z-axis, shown in blue) with adequate angular momentum generate both retrograde z-tubes (parented by the x4 family) and the long-axis tubes (parented by the stable anomalous orbits). These orbits appear as blue dots along the horizontal (0, 1, −1) resonance and at the top of the diagonal (1, −1, 0) resonance. The label "I" shows the location of the inner long-axis tubes, and the label "O" shows the location of the outer long-axis tubes. Finally, orbits launched close to the intermediate y-axis in the rotating Dehnen model with sufficient angular momentum generate retrograde z-tube orbits (which appear as green dots along the (1, −1, 0) diagonal).

We now compare these frequency maps with those generated for the orbits drawn from N-body bar distribution functions (Figure 13) for Model A (left) and Model B (right). Accurate determination of fundamental orbital frequencies requires that orbits are integrated for at least 20 orbital periods. Consequently, we exclude all orbits that were integrated for less than 20 orbital periods from the frequency maps. After excluding orbits with short integration times, the frequency maps in the top row of this figure show 5704 orbits for Model A and 6168 orbits in Model B.¹¹ The two lower panels zoom into the region of each map occupied by boxes. In both models the excluded particles were visually classified as belonging to the disk, so the frequency maps are still excellent representations of the distribution functions of the bars.

In Figure 13 the dots indicate regular orbits, while plus signs indicate chaotic orbits. The colors signify the value of an orbit's Jacobi Integral E_J (and not its initial position as in the previous maps). Particles in each map were divided into three equal groups in E_J: most negative values are colored blue, the least negative E_J values are colored red, and green indicates the intermediate range. In Model A most of the chaotic orbits are red and lie in the region labeled SPO/LPO/Disk—which marks the transition region between the disk and the bar (this is expected, e.g., from the work of Voglis et al. 2007, Contopoulos & Harsoula 2013, and Patsis & Katsanikas 2014b). These are the short-period orbits (SPOs) and long-period orbits (LPOs) that are temporarily trapped around the L4 and L5 Lagrange points. These orbits are not present in the Dehnen model (a pure ellipsoid) since they arise at the transition between the bar and disk. In Model B chaotic orbits appear in all energy ranges as a consequence of scattering by the central point mass.

Zoom In Zoom Out Reset image size

In both maps points cluster along "resonances" marked by dashed lines labeled with frequency ratios. The x1 orbits and z-tube orbits parented by x2 and x4 orbits lie in the same part of the frequency map as they do in the rotating Dehnen model (Figure 12, lower panel).

A "cloud" of box orbits appears to the left of the (1, −1, 0) line and below the (0, 1, −1) line in the same part of the frequency map as box orbits in the rotating Dehnen model in Figure 12 (bottom panel). The lower two panels in Figure 13 zoom into this region to more clearly show the various resonant boxlets. The vertical line labeled "2:0:−1" marks the boxy banana orbits. The x1 bananas lie at the intersection between "2:0:−1" and "1:−1:0" and are periodic orbits that satisfy the condition "${{\rm{\Omega }}}_{x}:{{\rm{\Omega }}}_{y}:{{\rm{\Omega }}}_{z}\quad =$ 2:−2:−1" (these are also referred to as "x1v1" orbits by Skokos et al. 2002a). Also seen is the 3:−2:0 resonance. Although the "double pretzel"-like orbits (Figure 4, top row) look slightly different from the "double fish"-like orbits (Figure 4, second and third rows), they both belong to this resonance family.

The frequency map for Model B shows essentially the same features as Model A except that a number of orbits with low values of E_J (blue) are scattered around the map. These orbits were scatted by the growing central point mass and associated with inner x-tubes (near the label "I"), outer x-tubes (near label "O"), and prograde z-tubes.