13.4. Orbits in triaxial systems¶
Unfortunately, in static triaxial potentials we cannot easily use surfaces of section, because static triaxial potentials only have a single classical integral of the motion (the specific energy) and this does not suffice to reduce the dimensionality of phase-space to three dimensions at a given energy. The study of orbits in triaxial systems is therefore far more difficult. However, there is also greater urgency in understanding the orbital structure of triaxial systems, because it is a priori not at all clear that static, long-term-stable triaxial mass configurations can be self-consistently constructed out of their orbits. That is, it is not clear that the density distribution obtained by letting a large number of bodies move along their orbits in a given triaxial gravitational potential can be the density that satisfies the Poisson equation for this potential for any setting of the relative number of bodies on different orbits.
The issue of whether or not a given mass distribution can be self-consistently constructed out of their orbits is not one that we have given much importance to until now. We discussed equilibrium configurations of spherical systems in Chapter 5, where we derived the Eddington formula for determining the ergodic (isotropic) distribution function that self-consistently generates the density in Section 5.6.1. We saw there that it is possible for this ergodic distribution function to be negative, meaning that there isn’t a self-consistent ergodic distribution function for all spherical systems. But once we relax the requirement of isotropy, it is clear that there is always an anisotropic distribution function that can self-consistently generate a given spherical density: we can simply build the density with shells made up of bodies on circular orbits, with the number of bodies proportional to the surface density of each shell. Similarly, we did not consider self-consistency much when discussing distributions of orbits in disk in Chapter 10 for two reasons. First of all, galactic disks are embedded in massive dark-matter halos and the disk therefore only ever contributes part of the gravitating mass density (albeit a not-insubstantial part). Second of all, as for spherical systems, we can build razor-thin disk distributions out of circular orbits with the number of bodies on each orbit set to match the disk’s surface density (e.g., Chapter 10.2.1) and it is quite plausible that we can generate any realistic warm, thick disk by heating such a cold distribution function in the radial and vertical direction (e.g., Chapter 10.2.2 and 10.4). Thus, it is clear that spherical and axisymmetric disky systems can be relatively-easily self-consistently generated (although we will side-step the question of whether such configurations are stable until Chapter 20.1, which especially for disks can be an issue).
Whether or not self-consistent triaxial configurations are possible is not immediately clear. The obvious generalization of the existence “proof” of self-consistent solutions for spherical and disk systems from the previous paragraph is to build a self-consistent solution out of the generalizations of circular orbits. As we saw in Section 13.2 above, in two-dimensional non-axisymmetric systems (an easy toy stand-in for full three-dimensional triaxial systems that we will rely on further below), generalized circular orbits are the closed loop orbits. But let’s see what happens when we look at a closed loop orbit in, for example, the non-axisymmetric logarithmic potential that we employed in Section 13.2 above. An example is shown in the left panel of Figure 13.18.
[61]:
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
lp= LogarithmicHaloPotential(normalize=True,b=0.8,core=0.2)
ts= numpy.linspace(0.,10.,601)
figure(figsize=(8,4))
subplot(1,2,1)
lp.plotSurfaceDensity(xmin=-0.43,xmax=0.43,ymin=-0.43,ymax=0.43,cntrls='--',
justcontours=True,log=True,ncontours=5,gcf=True)
o= Orbit([0.147,0.,1.07630163,0.])
o.integrate(ts,lp)
o.plot(gcf=True,lw=2.)
o= Orbit([0.147,0.15,1.07630163,0.])
o.integrate(ts,lp)
o.plot(gcf=True,zorder=0)
xlim(-0.43,0.43)
ylim(-0.43,0.43)
gca().set_aspect('equal')
subplot(1,2,2)
lp.plotSurfaceDensity(xmin=-0.43,xmax=0.43,ymin=-0.43,ymax=0.43,cntrls='--',
justcontours=True,log=True,ncontours=5,gcf=True)
o= Orbit([0.147,1.07630163,0.,0.])
o.integrate(ts,lp)
o.plot(gcf=True,lw=2.)
o= Orbit([0.147,1.07630163,0.2,0.])
o.integrate(ts,lp)
o.plot(gcf=True,zorder=0)
xlim(-0.43,0.43)
ylim(-0.43,0.43)
gca().set_aspect('equal')
tight_layout();
Figure 13.18: Loop (left) and box (right) orbits vs. equidensity contours in a non-axisymmetric potential.
The closed loop orbit is shown as the solid curve elongated in the \(y\) direction and we have also drawn three dashed contours of constant (surface) density. From this we see that the closed loop orbit is elongated in the direction perpendicular to the direction along which the density is elongated! This remains the case if we go to non-closed loop orbits (an example orb it is also shown in Figure 13.18) and it is generally true for loop orbits in non-axisymmetric potentials. Because they are elongated in the wrong direction, loop orbits cannot be a significant contributor to the density of self-consistent non-axisymmetric systems. Instead, box orbits need to contain a substantial fraction of the mass, because they are elongated in the correct direction. For example, a closed and non-closed box orbit with a similar energy as the loop orbits in the left panel are displayed in the right panel of Figure 13.18. Thus, box orbits are crucial for building triaxial mass distributions and whether or not a long-term stable, self-consistent triaxial density is possible therefore depends on whether or not box orbits exist and are stable.
In this section, we therefore discuss the properties of orbits in three-dimensional triaxial systems in detail with special attention paid to the box orbits. We start by discussing the major orbit families in triaxial systems, which will turn out to be box orbits and different generalizations of the loop orbits. Then we will move on to a more in-depth discussion of when box orbits exist, how they may be destroyed, and whether they can exist in realistic models for elliptical galaxies and dark matter halos.
13.4.1. The major orbit families in triaxial gravitational potentials: orbits in the perfect ellipsoid¶
Visualizing and studying orbits in six-dimensional phase-space is difficult and in the previous parts, we made use of the fact that spherical and axisymmetric systems are separable in various ways to simplify the investigation of their orbits. Unfortunately, the lack of symmetries in static triaxial mass distributions beyond the time invariance that leads to energy conservation, means that orbits in triaxial systems cannot in general be simplified by separating components. The exception to this rule are the Staeckel potentials. These are the most general class of potentials for which the Hamilton-Jacobi equation from Chapter 3.4.3 can be solved using the separation-of-variables technique. As discussed in Chapter 3.4.3, if the Hamilton-Jacobi equation can be separated, then orbits are regular (i.e., not chaotic) and are fully described as the combination of three independent oscillations in a two-dimensional phase-space each. We already made use of this in Chapter 9.4.3 to approximate disky potentials as oblate Staeckel potentials and we employed this approximation to compute action-angle coordinates. The most general coordinate system in which the Hamilton-Jacobi equation can be solved using the separation-of-variables technique is that of ellipsoidal coordinates, a generalization of the oblate and prolate spheroidal coordinates that we discussed in Chapter 12.2.1. Ellipsoidal coordinates \((\lambda,\mu,\nu)\) for a given Cartesian point \((x,y,z)\) are defined as the three roots of the equation \begin{equation}\label{eq-triaxialorbs-ellipsoidal-coords-definition} {x^2 \over \tau + \alpha} + {y^2 \over \tau + \beta} +{z^2 \over \tau + \gamma} = 1\,, \end{equation} where \((\alpha,\beta,\gamma)\) are constants that define the coordinate system that satisfy \(\alpha < \beta < \gamma\). We further order the roots such that \(\nu < \mu < \lambda\); the roots then satisfy \begin{equation} -\gamma \leq \nu \leq -\beta \leq \mu \leq -\alpha \leq \lambda\,. \end{equation} Because \(\lambda\), \(\mu\), and \(\nu\) cover different ranges, in many situations we can also represent them using a single variable \(\tau\) that either denotes \(\lambda\), \(\mu\), or \(\nu\) depending on which range it falls in. For completeness, we note that the inverse transformation from ellipsoidal to Cartesian coordinates is given by \begin{align} x^2 & = {(\lambda + \alpha)(\mu + \alpha)(\nu + \alpha) \over (\alpha-\beta)(\alpha-\gamma)}\,;\quad y^2 = {(\lambda + \beta)(\mu + \beta)(\nu + \beta) \over (\beta-\alpha)(\beta-\gamma)}\,;\quad z^2 = {(\lambda + \gamma)(\mu + \gamma)(\nu + \gamma) \over (\gamma-\alpha)(\gamma-\beta)}\,. \end{align} Each of \((x,y,z)\) is therefore only specified up to a sign. The transformation equations for the momenta can be derived using the coordinate-transformation generation function from Chapter 3.4.2. Staeckel potentials are then potentials that in ellipsoidal coordinates can be written as \begin{equation}\label{eq-triaxialorbs-staeckel} \Phi(\lambda,\mu,\nu) = {f(\lambda) \over (\lambda-\mu)(\nu-\lambda)} + {f(\mu) \over (\mu-\nu)(\lambda-\mu)} + {f(\nu) \over (\nu-\lambda)(\mu-\nu)} \,. \end{equation} For this form of the potential, the Hamilton-Jacobi equation in ellipsoidal coordinates can be solved by separating the variables.
Few galactic potentials can be written in the form of Equation (13.8) in any ellipsoidal coordinate system, but one example is the perfect ellipsoid that we introduced in Chapter 12.2.3. In fact, the perfect ellipsoid is the only gravitational potential of Staeckel form for which the density is stratified on similar ellipsoids (\(\rho\equiv \rho(m)\) for \(m\) given by Equation 12.32; de Zeeuw & Lynden-Bell 1985). Thus, it is the only realistic triaxial galactic potential for which the Hamilton-Jacobi equation can be solved using the separation-of-variables technique and where orbits, therefore, can be easily studied by separating them into lower-dimensional components. As such, the perfect ellipsoid is the perfect potential to learn about the major orbit families that exist in triaxial galactic mass distributions!
More specifically, the perfect ellipsoid potential for a given set of axis lengths \(a > b > c\) can be written in Staeckel form when we set \(\alpha = -a^2\), \(\beta = -b^2\), and \(\gamma = -c^2\). Because the orbits separate into decoupled motions in \(\lambda\), \(\mu\), and \(\nu\), we can get a first sense of the orbital structure by considering what surfaces with constant values of one of these coordinates look like: orbits oscillate between such surfaces. For this, we first implement the transformation between Cartesian and ellipsoidal coordinates:
[62]:
def cartesian_to_ellipsoidal(x,y,z,alpha,beta,gamma):
x= numpy.atleast_1d(x)
y= numpy.atleast_1d(y)
z= numpy.atleast_1d(z)
N= len(x)
out= numpy.empty((N,3))
for ii,(tx,ty,tz) in enumerate(zip(x,y,z)):
these_coords= numpy.polynomial.polynomial.Polynomial(\
(beta*gamma*tx**2.+alpha*gamma*ty**2.+alpha*beta*tz**2.
-alpha*beta*gamma,
(beta+gamma)*tx**2.+(alpha+gamma)*ty**2.+(alpha+beta)*tz**2.
-alpha*beta-alpha*gamma-beta*gamma,
tx**2.+ty**2.+tz**2.-alpha-beta-gamma,
-1.)).roots()
out[ii]= sorted(these_coords)[::-1]
return out
def ellipsoidal_to_cartesian(l,m,n,alpha,beta,gamma):
x= numpy.sqrt((l+alpha)*(m+alpha)*(n+alpha)/(alpha-beta)/(alpha-gamma))
y= numpy.sqrt((l+beta)*(m+beta)*(n+beta)/(beta-alpha)/(beta-gamma))
z= numpy.sqrt((l+gamma)*(m+gamma)*(n+gamma)/(gamma-beta)/(gamma-alpha))
return numpy.array([x,y,z]).T
Then we can look at surfaces of constant \(\lambda\), \(\mu\), and \(\nu\) for a potential with \(a = 1\), \(b = 0.7\), and \(c = 0.5\) in Figure 13.19.
[63]:
from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.colors import Normalize
def draw_coordinate_axes(ax,size=2.5):
# Draw centered axes (https://stackoverflow.com/a/61927456)
val = [size,0,0]
labels = [r'$x$', r'$y$', r'$z$']
for v in range(3):
x = [val[v-0], -val[v-0]]
y = [val[v-1], -val[v-1]]
z = [val[v-2], -val[v-2]]
ax.plot(x,y,z,'k-', linewidth=1)
ax.text(val[v-0]-0.1,val[v-1],val[v-2]-0.25,
labels[v],color='k',fontsize=20)
return None
a, b, c= 1., 0.7, 0.5
alpha, beta, gamma= -a**2., -b**2., -c**2.
fig= figure(figsize=(16.,4.8))
ax= fig.add_subplot(1,3,1,projection='3d')
# Surface of constant lambda
la,_,_= cartesian_to_ellipsoidal(1.,0.,0.,alpha,beta,gamma).T
mus= numpy.linspace(-beta,-alpha,11)
nus= numpy.linspace(-gamma,-beta,11)
munus= numpy.meshgrid(mus,nus,indexing='ij')
xyz= ellipsoidal_to_cartesian(la[0]+numpy.zeros_like(munus[0]),
munus[0],munus[1],alpha,beta,gamma)
norm= Normalize(-numpy.nanmax(xyz[:,:,2]),numpy.nanmax(xyz[:,:,2]))
for sgnx in [-1,1]:
for sgny in [-1,1]:
for sgnz in [-1,1]:
ax.plot_trisurf(sgnx*xyz[:,:,0].flatten(),
sgny*xyz[:,:,1].flatten(),
sgnz*xyz[:,:,2].flatten(),
cmap=cm.plasma,norm=norm,
linewidth=10.2,antialiased=True)
draw_coordinate_axes(ax,size=2.5)
ax.view_init(40,50)
ax.set_xlim(-1.2,1.2)
ax.set_ylim(-1.2,1.2)
ax.set_zlim(-1.2,1.2)
ax._axis3don= False
annotate(r'$\mathrm{constant}\ \lambda$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
ax= fig.add_subplot(1,3,2,projection='3d')
# Surface of constant mu
_,ma,_= cartesian_to_ellipsoidal(0.3,0.4,0.3,alpha,beta,gamma).T
las= numpy.linspace(-alpha,-alpha+3.,11)
nus= numpy.linspace(-gamma,-beta,11)
lanus= numpy.meshgrid(las,nus,indexing='ij')
xyz= ellipsoidal_to_cartesian(lanus[0],ma[0]+numpy.zeros_like(lanus[0]),
lanus[1],alpha,beta,gamma)
norm= Normalize(-numpy.nanmax(xyz[:,:,2]),numpy.nanmax(xyz[:,:,2]))
for sgnx in [-1,1]:
for sgny in [-1,1]:
for sgnz in [-1,1]:
ax.plot_trisurf(sgnx*xyz[:,:,0].flatten(),
sgny*xyz[:,:,1].flatten(),
sgnz*xyz[:,:,2].flatten(),
cmap=cm.plasma,norm=norm,
linewidth=10.2,antialiased=True)
draw_coordinate_axes(ax,size=2.5)
ax.view_init(40,50)
ax.set_xlim(-1.2,1.2)
ax.set_ylim(-1.2,1.2)
ax.set_zlim(-1.2,1.2)
ax._axis3don= False
annotate(r'$\mathrm{constant}\ \mu$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
ax= fig.add_subplot(1,3,3,projection='3d')
# Surface of constant nu
_,_,na= cartesian_to_ellipsoidal(0.5,0.6,0.7,alpha,beta,gamma).T
las= numpy.linspace(-alpha,-alpha+2.,11)
mus= numpy.linspace(-beta,-alpha,11)
lamus= numpy.meshgrid(las,mus,indexing='ij')
xyz= ellipsoidal_to_cartesian(lamus[0],lamus[1],
na[0]+numpy.zeros_like(lamus[0]),alpha,beta,gamma)
norm= Normalize(-numpy.nanmax(xyz[:,:,2]),numpy.nanmax(xyz[:,:,2]))
for sgnx in [-1,1]:
for sgny in [-1,1]:
for sgnz in [-1,1]:
ax.plot_trisurf(sgnx*xyz[:,:,0].flatten(),
sgny*xyz[:,:,1].flatten(),
sgnz*xyz[:,:,2].flatten(),
cmap=cm.plasma,norm=norm,
linewidth=10.2,antialiased=True)
draw_coordinate_axes(ax,size=2.5)
ax.view_init(40,50)
ax.set_xlim(-1.2,1.2)
ax.set_ylim(-1.2,1.2)
ax.set_zlim(-1.2,1.2)
ax._axis3don= False
annotate(r'$\mathrm{constant}\ \nu$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.);
Figure 13.19: Surfaces of constant \(\lambda\), \(\mu\), and \(\nu\) in the ellipsoidal coordinate system.
Surfaces of constant \(\lambda\) are ellipsoids with the long axis in the \(z\) direction and the short axis in the \(x\) direction. For \(\lambda \gg -\alpha\), these ellipsoids are approximately spherical, while for \(\lambda = -\alpha\), the ellipsoid collapses to a 2D ellipse at \(x=0\). Surface of constant \(\mu\) are hyperboloids of one sheet and surfaces of constant \(\nu\) are hyperbolae of two sheets. The one-sheet hyperboloids collapse to a 2D hyperbola at \(y=0\) for \(\mu = -\beta\), while the two-sheet hyperboloid collapses to a two-dimensional area at \(z=0\). At small distances from the origin, the ellipsoidal coordinate system is approximately the Cartesian, while at large distances it is approximately spherical; the ellipsoidal coordinate system therefore smoothly interpolates between these two extremes. This is similar to how the oblate and prolate coordinate systems that we used in Chapter 12.2.1 go from being approximately the cylindrical coordinate frame at small distances from the origin to being approximately spherical at large distances.
Because of the separability of the Hamilton-Jacobi equation for the perfect ellipsoid, one can compute two integrals of the motion in addition to the specific energy that are generalizations of the \(z\)-component of the angular momentum and the third integral for axisymmetric systems; we denote these as \(I_2\) and \(I_3\). Similar to what is possible for spherical and axisymmetric systems, the motion in the perfect ellipsoid can be described as that in a one-dimensional effective potential \(\Phi_{\mathrm{eff}}(\tau)\) for a given value of \(I_2\) and \(I_3\) (where we write this potential as a function of \(\tau\) to cover all three coordinates as explained above). It turns out that motion is then possible in regions where \(\Phi_{\mathrm{eff}}(\lambda) \leq E\), \(\Phi_{\mathrm{eff}}(\mu) \leq E\), and \(\Phi_{\mathrm{eff}}(\nu) \geq E\); note that the inequality in the last equation is opposite that of the others (and opposite to the equivalent condition for spherical or axisymmetric systems). Turning points of the orbit are then those \(\tau\) where \(\Phi_{\mathrm{eff}}(\tau) = E\) and this equation is satisfied for all orbits at three points \(\tau_i\) (up from two in the spherical and axisymmetric case). Depending on where the \(\tau_i\) lie relative to \(-\gamma\), \(-\beta\), and \(-\alpha\), we get different types of orbits.
From Equation (13.5) it is clear that the origin \((x,y,z) = (0,0,0)\) is mapped onto \((\lambda,\mu,\nu) = (-\alpha,-\beta,-\gamma)\). Orbits with the three turning points \(\tau_i\) arranged as \(-\gamma < \tau_1 < -\beta < \tau_2 < -\alpha < \tau_3\) are therefore able to pass through the origin. In each ellipsoidal coordinate direction, they oscillate between the two-dimensional area at \(x,y,z=0\) at \(\lambda,\mu,\nu = -\alpha,-\beta,-\gamma\) and the ellipsoid or hyperboloid of their turning point. The three-dimensional volume is therefore a box and these are the three dimensional box orbits. It turns out that there are three other possible arrangements of the turning points and they all have two turning points either between \(-\beta\) and \(-\alpha\) or above \(-\alpha\). Such orbits do not include \((\lambda,\mu,\nu) = (-\alpha,-\beta,-\gamma)\) and therefore do not pass through the origin. These orbits have a definite sense of rotation around one of the axes and are the generalizations of the loop orbits in axisymmetric potentials; in the context of triaxial systems they are known as tube orbits. Orbits with \(-\gamma < \tau_1 < -\beta < -\alpha < \tau_2 < \tau_3\) oscillate between two ellipsoids in \(\lambda\) and between the two branches of the turning hyperboloid of two sheets corresponding to \(\tau_1\). From the constant \(\nu\) and constant \(\lambda\) shapes in the figure above, it is clear that such orbits loop around the minor, \(z\) axis. These orbits are known as the short-axis tube orbits. The remaining two options both have one or two turning hyperboloids of one sheet, which avoid the \(x\) axis and these orbits therefore loop around the major, \(x\) axis; they are known as long-axis tube orbits. They further separate into two classes depending on whether there are one or two turning hyperboloids of one sheet: \(-\gamma < -\beta < \tau_1 < \tau_2 < -\alpha < \tau_3\) is an inner long-axis tube orbit, while \(-\gamma < -\beta < \tau_1 < -\alpha < \tau_2 < \tau_3\) gives an outer long-axis tube orbit. Figure 13.20 displays a representative example of these four orbit families.
[80]:
from galpy.potential import PerfectEllipsoidPotential
from galpy.orbit import Orbit
from mpl_toolkits.mplot3d import Axes3D
def draw_coordinate_axes(ax,size=2.5,shrinkx=1.,shrinkz=1.):
# Draw centered axes (https://stackoverflow.com/a/61927456)
val = [size,0,0]
labels = [r'$x$', r'$y$', r'$z$']
for v in range(3):
x = [val[v-0]/shrinkx,-val[v-0]/shrinkx]
y = [val[v-1],-val[v-1]]
z = [val[v-2]/shrinkz, -val[v-2]/shrinkz]
ax.plot(x,y,z,'k-',linewidth=1)
ax.text(val[v-0]/shrinkx-0.1/shrinkx,
val[v-1],val[v-2]/shrinkz-0.25/shrinkz,
labels[v],color='k',fontsize=20)
return None
a, b, c= 1., 0.7, 0.5
pep= PerfectEllipsoidPotential(normalize=1.,a=a,b=b,c=c)
fig= figure(figsize=(10.,10),dpi=150)
ax= fig.add_subplot(2,2,1,projection='3d')
# box orbit
o= Orbit([0.5,-0.01,0.991,0.15,0.,numpy.pi/2.])
ts= numpy.linspace(0.,10000.,100001)
o.integrate(ts,pep)
ax.scatter(o.x(ts),o.y(ts),o.z(ts),marker='o',
s=6.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(40,50.)
ax.set_xlim(-1.,1.)
ax.set_ylim(-1.,1.)
ax.set_zlim(-0.25,0.25)
ax._axis3don= False
draw_coordinate_axes(ax,size=1.5,shrinkz=2.50)
annotate(r'$\mathrm{box}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
# Short axis tube
ax= fig.add_subplot(2,2,2,projection='3d')
o= Orbit([1.2,-0.21,0.61,0.4,-0.2,numpy.pi/2.])
ts= numpy.linspace(0.,10000.,100001)
o.integrate(ts,pep)
ax.scatter(o.x(ts),o.y(ts),o.z(ts),marker='o',
s=6.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(40,50.)
ax.set_xlim(-1.,1.)
ax.set_ylim(-1.,1.)
ax.set_zlim(-.6,.6)
ax._axis3don= False
draw_coordinate_axes(ax,size=1.85,shrinkz=1.25)
annotate(r'$\mathrm{short-axis\ tube}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
# Inner long-axis tube orbit
ax= fig.add_subplot(2,2,3,projection='3d')
o= Orbit([0.4,0.51,0.51,.6,-0.195,0.])
ts= numpy.linspace(0.,2000.,100001)
o.integrate(ts,pep)
ax.scatter(o.x(ts)*1.25,o.y(ts)*1.25,o.z(ts)*1.25,
marker='o',s=9.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(20,40)
ax.set_xlim(-.75,.75)
ax.set_ylim(-1.,1.)
ax.set_zlim(-1.,1.)
ax._axis3don = False
draw_coordinate_axes(ax,size=1.65)
annotate(r'$\mathrm{inner\ long-axis\ tube}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
# Outer long-axis tube orbit
ax= fig.add_subplot(2,2,4,projection='3d')
o= Orbit([2.9,0.31,0.41,2.4,-0.675,numpy.pi/2.])
ts= numpy.linspace(0.,10000.,100001)
o.integrate(ts,pep)
ax.scatter(o.x(ts)/4.,o.y(ts)/4.,o.z(ts)/4.,
marker='o',s=6.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(20,30.)
ax.set_xlim(-1.,1.)
ax.set_ylim(-1.,1.)
ax.set_zlim(-1.,1.)
ax._axis3don = False
draw_coordinate_axes(ax,size=1.6,shrinkx=0.8)
annotate(r'$\mathrm{outer\ long-axis\ tube}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
tight_layout();
Figure 13.20: The four main orbit families in the perfect ellipsoid.
These four orbit families constitute the main orbit families in all triaxial mass distributions. The box orbits and the short-axis tube orbits are the generalizations of the box and loop orbits in two-dimensional, non-axisymmetric potentials that we saw in Section 13.2 above and those loop orbits were themselves the generalizations of the loop orbits in axisymmetric mass distributions. However, there is no analog to the long-axis tube orbits in axisymmetric systems and these orbits only exist as stable orbit families in triaxial mass distributions. Observationally, these manifest themselves as apparent rotation along the minor axis and this minor-axis rotation is therefore a tell-tale observational signature of triaxiality (see the discussion in Chapter 12.1). That there is no intermediate-axis tube orbit family means that there are no stable orbits that loop around the intermediate axis (Heiligman & Schwarzschild 1979).
As an example of these four main orbit families in other triaxial potentials, we can look at similar orbits as those shown above, but now computed in a triaxial NFW potential; this is shown in Figure 13.21.
[81]:
from galpy.potential import TriaxialNFWPotential
from galpy.orbit import Orbit
from mpl_toolkits.mplot3d import Axes3D
def draw_coordinate_axes(ax,size=2.5,shrinkx=1.,shrinkz=1.):
# Draw centered axes (https://stackoverflow.com/a/61927456)
val = [size,0,0]
labels = [r'$x$', r'$y$', r'$z$']
for v in range(3):
x = [val[v-0]/shrinkx,-val[v-0]/shrinkx]
y = [val[v-1],-val[v-1]]
z = [val[v-2]/shrinkz, -val[v-2]/shrinkz]
ax.plot(x,y,z,'k-',linewidth=1)
ax.text(val[v-0]/shrinkx-0.1/shrinkx,
val[v-1],val[v-2]/shrinkz-0.25/shrinkz,
labels[v],color='k',fontsize=20)
return None
a, b, c= 1., 0.7, 0.5
nep= TriaxialNFWPotential(normalize=1.,a=a,b=b,c=c)
fig= figure(figsize=(10.,10),dpi=150)
ax= fig.add_subplot(2,2,1,projection='3d')
# box orbit
o= Orbit([0.1,-0.7,0.9,0.15,0.,numpy.pi/2.])
ts= numpy.linspace(0.,10000.,100001)
o.integrate(ts,nep,method='dop853_c')
ax.scatter(o.x(ts)*2,o.y(ts)*2,o.z(ts)*2,marker='o',
s=6.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(40,50.)
ax.set_xlim(-1.,1.)
ax.set_ylim(-1.,1.)
ax.set_zlim(-1.,1.)
ax._axis3don= False
draw_coordinate_axes(ax,size=1.5,shrinkz=0.8)
annotate(r'$\mathrm{box}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
# Short axis tube
ax= fig.add_subplot(2,2,2,projection='3d')
o= Orbit([1.2,-0.21,0.61,0.4,-0.2,numpy.pi/2.])
ts= numpy.linspace(0.,10000.,100001)
o.integrate(ts,nep,method='dop853_c')
ax.scatter(o.x(ts),o.y(ts),o.z(ts),marker='o',
s=6.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(40,50.)
ax.set_xlim(-1.,1.)
ax.set_ylim(-1.,1.)
ax.set_zlim(-.6,.6)
ax._axis3don= False
draw_coordinate_axes(ax,size=1.85,shrinkz=1.25)
annotate(r'$\mathrm{short-axis\ tube}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
# Inner long-axis tube orbit
ax= fig.add_subplot(2,2,3,projection='3d')
o= Orbit([0.4,0.51,0.51,.6,-0.195,0.])
ts= numpy.linspace(0.,2000.,100001)
o.integrate(ts,nep,method='dop853_c')
ax.scatter(o.x(ts)*1.25,o.y(ts)*1.25,o.z(ts)*1.25,
marker='o',s=9.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(20,40)
ax.set_xlim(-.75,.75)
ax.set_ylim(-1.,1.)
ax.set_zlim(-1.,1.)
ax._axis3don = False
draw_coordinate_axes(ax,size=1.65)
annotate(r'$\mathrm{inner\ long-axis\ tube}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
# Outer long-axis tube orbit
ax= fig.add_subplot(2,2,4,projection='3d')
o= Orbit([2.9,0.31,0.41,2.4,-0.675,numpy.pi/2.])
ts= numpy.linspace(0.,10000.,100001)
o.integrate(ts,nep,method='dop853_c')
ax.scatter(o.x(ts)/4.,o.y(ts)/4.,o.z(ts)/4.,marker='o',
s=6.,edgecolors=None,alpha=0.01,rasterized=True)
ax.view_init(20,30.)
ax.set_xlim(-1.,1.)
ax.set_ylim(-1.,1.)
ax.set_zlim(-1.,1.)
ax._axis3don = False
draw_coordinate_axes(ax,size=1.6,shrinkx=0.8)
annotate(r'$\mathrm{outer\ long-axis\ tube}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
tight_layout();
Figure 13.21: The four main orbit families in a triaxial NFW profile.
We see that these orbits are qualitatively similar to those in the perfect ellipsoid and the main orbit families are therefore the same in a triaxial NFW potential as in the perfect ellipsoid. But the two orbit families that approach the center the closest—the box orbits and the inner long-axis tubes—look somewhat different. This is because the inner density of the NFW potential is very different from that of the perfect ellipsoid: the perfect ellipsoid has a constant-density core at \(m \ll 1\), while the NFW density diverges as \(\rho(m) \propto m^{-1}\). We see in particular that the box orbit is more concentrated towards the center and more fanned out at large distances than a similar box orbit in the perfect ellipsoid. Because box orbits are so important for building self-consistent models of triaxial mass distributions, this spells trouble.
13.4.2. Orbits in non-axisymmetric potentials with central density cusps¶
13.4.2.1. General considerations¶
As we discussed in Chapter 12.1, elliptical galaxies have inner surface-density profiles that are power-laws: \(I(R) \propto R^{-\gamma}\) with \(0 < \gamma \lesssim 1\). In particular, the surface density keeps increasing at the smallest scale that can be resolved in all elliptical galaxies, albeit very shallowly in luminous elliptical galaxies (\(0 < \gamma \lesssim 0.3\)). Thus, unlike the perfect ellipsoid from the previous section or the two-dimensional cored logarithmic potential from Section 13.2 above, elliptical galaxies have cuspy rather than cored profiles. The same holds for dark-matter halos, where \(\rho(r) \propto r^{-1}\) (Section 2.4.6). However, we derived the existence of box orbits in Section 13.2 from the approximately-harmonic motion that occurs in a constant-density mass distribution. If the gravitational potential is not harmonic however close we get to the center, it is therefore not clear that box orbits exist or are long-term stable.
Another reason that elliptical galaxies have a central non-harmonic potential is the existence of supermassive black holes at their centers (we will discuss the evidence for this further in Chapter 16.1). The presence of a supermassive black hole has various dynamical consequences. The most obvious is that the black hole constitutes a point-mass potential and that the motion near the black hole is therefore Keplerian rather than harmonic. Stars that orbit within the region where the black hole dominates the mass distribution—the black hole’s so called sphere of influence, which is about \(5\) to \(50\,\mathrm{pc}\) depending on its mass—are therefore Keplerian ellipses rather than box orbits. But the larger-scale orbital distribution is also affected by the supermassive black hole, because any box orbit eventually enters the black hole’s sphere of influence on a hyperbolic trajectory and exits it on a different orbit due to the deflection from the interaction with the black hole. This tends to disrupt the skinny box orbits that are crucial for building triaxial systems. A secondary effect occurs when the black hole grows smoothly: then it is predicted to adiabatically grow a steep cusp of stars with \(\rho(r) \propto r^{-7/4}\) (Bahcall & Wolf 1976) that further pushes the potential away from harmonicity (however, no such cusps have ever been observed). Numerical studies of orbits in two-dimensional (e.g., Gerhard & Binney 1985) and three-dimensional non-axisymmetric mass distributions (e.g., Merritt & Valluri 1996; Merritt & Quinlan 1998; Valluri & Merritt 1998) demonstrate that supermassive black holes with masses exceeding \(\approx 0.5\,\%\) of the galaxy’s mass significantly disrupt the box orbits even in an otherwise cored gravitational potential. However, enough box orbits survive such that combined with weakly-chaotic orbits they can support the global triaxiality of such models even when the black hole is as massive as \(1\,\%\) of the galaxy’s mass (Holley-Bockelmann et al. 2002). Because supermassive black holes typically have masses of \(\approx 0.2\) to \(0.3\,\%\) of the galaxy’s mass (e.g., Magorrian et al. 1998; Häring & Rix 2004; Gültekin et al. 2009), orbital evolution induced by them is therefore likely not that important, at least on scales beyond the inner tens of parsecs where the black hole dominates the gravitational potential. We leave a more detailed exploration of this topic therefore as an exercise and focus here on the effect of density cusps.
The density cusps we are interested in behave as \(\rho(r) \propto r^{-\gamma}\) with \(\gamma \approx 1\) for dark-matter halos (see Chapter 2.4.6), and \(0 \lesssim \gamma \lesssim 1\) or \(1 \lesssim \gamma \lesssim 2\) for the shallow or steep cusps at the centers of luminous and less luminous elliptical galaxies (see Chapter 12.1). From the behavior of spherical power-law density profiles that we discussed in Chapter 2.4.5, it is clear that the force at the origin is finite for \(\gamma \leq 1\), while it diverges for \(\gamma > 1\). While the force diverges, the circular velocity remains finite until \(\gamma > 2\); this indicates that while the central force may diverge, the orbital effects it gives rise to are limited. At \(\gamma > 2\), both the central force and the circular velocity diverge and orbits that pass near the center are strongly affected by the central cusp. Therefore, there are three important regimes (Gerhard & Binney 1985):
\(\gamma \leq 1\) (\(\rho[r]\) shallower than \(1/r\)): a density cusps exists, but the central force and circular velocity remain finite. In a triaxial system, the dynamics is therefore not qualitatively different from that in a constant-density core and we expect box orbits to exist and be stable. Dark matter halos and luminous elliptical galaxies fall in this regime and we therefore expect few problems in building steady-state, stable triaxial configurations for these systems.
\(\gamma > 2\) (\(\rho[r]\) steeper than \(1/r\)): the central force and circular velocity increase without bound as \(r \rightarrow 0\). Passages near the origin strongly affect a body’s trajectory, which gets deflected every time it passes near the origin. We therefore expect box orbits to be destroyed. Only some of the less luminous elliptical galaxies fall within this range.
\(1 < \gamma \leq 2\) (\(\rho[r]\) between than \(1/r\) and \(1/r^2\)): the central force diverges, but the circular velocity remains finite. Orbits that pass near the origin are affected by the central cusp, but the extent to which they are affected is unclear without further investigation. Less luminous elliptical galaxies fall within this range.
What happens when \(\gamma\) is increased significantly above one is that box orbits get replaced or destroyed by other types of orbits. These orbits fall in two categories: (i) regular orbits that avoid the center (centrophobic orbits, as opposed to the centrophilic box orbits)—we will refer to these as boxlets—and (ii) chaotic orbits. As \(\gamma\) increases, the chaotic orbits start appearing by originating from box orbits that pass through the center and that are deflected by the central density cusp. For \(\gamma \approx 1\), the deflection of a single passage through the origin is weak and the resulting orbit is only weakly chaotic, appearing to be regular over long times (that is, equal to or longer than the age of the Universe); as \(\gamma\) increases, the deflections become larger and orbits become more strongly chaotic on time scales similar to the dynamical time. Both of these types of orbits are bad news for building a steady-state, self-consistent triaxial model. As we saw above, box orbits, and especially very narrow ones, are crucial for supporting the triaxial figure of the mass distribution. Because they pass close to the center more often than other box orbits, the very narrow ones are the first to go once a central density cusp is introduced. Chaotic orbits are particularly problematic for supporting the triaxial figure: strongly chaotic orbits in a static triaxial potential will tend to quickly fill the entire volume of phase-space allowed by their single integral of the motion, the specific energy \(E\). The spatial part of this volume is the equipotential surface \(\Phi(\vec{x}) = E\). However, in general the gravitational potential is much rounder than the density (see, e.g., the discussion around Equation 7.20) and strongly chaotic orbits can therefore not significantly contribute to a steady, self-consistent model for a triaxial system. Weakly-chaotic orbits do not fill their equipotential surface as quickly and can therefore be employed to build self-consistent models that are stable for the age of the Universe (e.g., Schwarzschild 1993; Merritt & Fridman 1996).
The dynamical conclusion of these arguments is therefore that we only expect to see strongly triaxial mass distributions in systems that have only a shallow cusp, \(\gamma \lesssim 1\), while systems with steeper cusps should be axisymmetric. This is indeed what is observed for elliptical galaxies. As we discussed in Chapter 12.1, luminous elliptical galaxies have shallow central-density cusps and triaxial shapes, while less luminous elliptical galaxies have steep central densities, but they are approximately axisymmetric. We therefore have a dynamically-coherent picture of the observed structure of elliptical galaxies. Why it is the case that luminous ellipticals are triaxial, while less luminous ones are not is a question for galaxy formation rather than galaxy dynamics.
Because the predicted central-density cusps of dark matter halos are shallow \(1/r\) cusps, dynamically, dark matter halos could be triaxial. Whether or not they are again depends on the details of their formation in the cosmological paradigm. As we discussed in Chapter 12.1, dark-matter-only simulations predict that dark matter halos are strongly triaxial, while simulations that include the growth of a baryonic galaxy component find axisymmetric halos. This happens because the triaxiality-supporting box orbits are deformed and destroyed by the growth of the baryonic component, as we will discuss in further detail below.
13.4.2.2. Boxlets in the logarithmic potential¶
To gain further understanding of what happens to box orbits as the central-density cusp becomes stronger, we look at what happens in the two-dimensional logarithmic potential from Section 13.2 when the core radius shrinks to zero and the central potential transitions from a constant-density core to a \(\rho(r) \propto r^{-2}\) cusp. To start, let’s look at the surface of section for the same potential as we used in Section 13.2 —a cored logarithmic potential with a core radius of \(0.2\) and a flattening of \(0.7\) in the \(y\) direction—but at the higher energy \(E=0\). We consider the \((y,v_y)\) surface of section (when orbits pass from negative to positive \(x\) through the \(x=0\) plane) in Figure 13.22.
[68]:
from galpy.potential import LogarithmicHaloPotential, evaluatePotentials
from galpy.orbit import Orbit
from scipy import optimize
figure(figsize=(9,6),dpi=300)
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,y,vy) with given E"""
return Orbit([y,vy,numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)))
lp= LogarithmicHaloPotential(amp=1.,b=0.7,core=0.2)
E= 0.
# Location of closed orbits in [y,vy]
# Note that the closed 'pretzel' isn't really
# a pretzel, and the closed 'box' is what's typically
# referred to as a pretzel...
# Don't worry too much about the names...
closed= {'banana': [0.0565,0.],
'pretzel': [0.3235,0.],
'pretzel2': [1e-4,1.52651722],
'box': [0.4005,0.],
'chaos': [],
'harbinger': [0.427,0.],
'loop': [0.5715,0.]}
nonclosed= {'banana': [[0.0001,0.]],
'pretzel': [[0.34,0.],[0.04,1.55],[0.04,1.525],[0.33,0.],[0.343,0.]],
'pretzel2': [],
'box': [[0.3,0.],[0.1,0.],[0.2,0.],[0.15,0.],[0.25,0.],[0.25,0.],
[0.34775,0.],[0.42,0.],[0.41,0.],[0.39,0.]],
'harbinger': [[0.429,0.]],
'chaos':[[0.44,0.]],
'loop': [[0.48,0.],[0.455,0.],[0.52,0.],[0.4475,0.]]}
colors= {'banana': 'tab:orange',
'box': 'tab:blue',
'pretzel': 'tab:purple',
'pretzel2': 'tab:purple',
'loop': 'tab:cyan',
'chaos': 'tab:pink',
'harbinger': 'tab:brown'}
for orb in closed:
for sgn in [-1,1]:
if not orb == 'chaos':
o= orbit_yvyE(sgn*closed[orb][0],closed[orb][1],E,pot=lp)
o.plotSOS(lp,ncross=1_000,method='dop853_c',surface='x',
gcf=True,overplot=(orb != 'banana') + (sgn != -1),
marker='.',mec='none',zorder=2,
mfc=colors[orb],rasterized=True,
xrange=[-0.75,0.75],yrange=[-1.9,1.9])
for nonc in nonclosed[orb]:
o= orbit_yvyE(sgn*nonc[0],nonc[1],E,pot=lp)
o.plotSOS(lp,ncross=1_000+(orb == 'chaos')*3_000,
method='dop853_c',surface='x',
overplot=True,rasterized=True,
marker='.',ms=3.,mec='none',zorder=2,mfc=colors[orb])
# Also plot the zero-velocity curve, first find min x at which it exists
miny= optimize.brentq(lambda y: E-evaluatePotentials(lp,y,0.,phi=numpy.pi/2.),
0.01,0.9)-1e-10
ys= numpy.linspace(-miny,miny,201)
plot(ys,zvc(ys,E,pot=lp),color='k',zorder=0)
plot(ys,-zvc(ys,E,pot=lp),color='k',zorder=1)
galpy_plot.text(r'$b=0.7$'+'\n'+r'$R_c=0.2$',top_left=True,size=18.);
Figure 13.22: Surface of section of the cored, non-axisymmetric logarithmic potential at \(E=0\).
Comparing this surface of section to that at \(E=-0.87\) from Figure 13.9 shows that this surface of section looks quite different. The \(E=-0.87\) surface of section was dominated by box orbits, with only a small fraction of loop orbits. At \(E=0\), however, loop orbits make up a substantial fraction of the section at large distances from the center (the cyan curves above). Box orbits dominate the central part of the section (the dark blue curves above), but very close to the center a new family of orbits appears that avoids the center (orbits parented by the orange dots). Two more small families of regular orbits appear as well (the purple and brown curves), and there is a closed orbit inside of the box orbit region (the blue dots) and a chaotic orbit between the loop and box orbits (the pink dots). All of the closed regular orbits shown here aside from that in the box-orbit region are stable: they are surrounded by qualitatively similar, non-closed orbits.
To further understand the new orbit families, the four new types of closed orbits that we find in this surface of section are shown in Figure 13.23.
[69]:
figure(figsize=(10,3.5))
ts= numpy.linspace(0.,100.,100001)
for ii,orb in enumerate(['banana','pretzel','box','harbinger']):
o= orbit_yvyE(closed[orb][0],closed[orb][1],E,pot=lp)
o.integrate(ts,lp,method='dop853_c')
subplot(1,4,ii+1)
o.plot(gcf=True,
xrange=[-1.,1.] if ii == 0 else [-0.675,0.675],
yrange=[-0.75,0.75])
plot([0.],[0.],'ko',ms=7.)
if ii > 0:
gca().get_yaxis().get_label().set_visible(False)
yticks([])
tight_layout();
Figure 13.23: Closed boxlets in the surface of section of Figure 13.22.
The orbit on the left corresponds to the orange dot in the surface of section and it is what the closed long-axis box orbit becomes at this energy. That is, the straight \(y=0\) orbit at \(E=-0.87\) becomes a bent orbit at \(E=0\). Because of its shape, this type of orbit is known as a banana orbit. The closed banana orbit avoids the center and is therefore centrophobic. The other closed orbits are more complicated. We see that two of them are centrophilic and pass through the origin, while the orbit that is surrounded by box orbits (the third orbit here, which corresponds to the blue points in the surface of section) is centrophobic. We will refer to all orbits like this as boxlets. Some examples of boxlets in the 1/r` cusps of dark-matter halos are shown in Section 13.4.3 below.
To get a better sense of how these orbits develop, let’s look at animations for some of these closed orbits. Animating the banana orbit, we get
[13]:
orb= 'banana'
ts= numpy.linspace(0.,100.,3001)
o= orbit_yvyE(closed[orb][0],closed[orb][1],E,pot=lp)
o.integrate(ts,lp,method='dop853_c')
o.animate(staticPlot=True); # remove the ; to display the animation
Similarly, for the pretzel-like orbit, we get
[14]:
orb= 'box'
ts= numpy.linspace(0.,100.,3001)
o= orbit_yvyE(closed[orb][0],closed[orb][1],E,pot=lp)
o.integrate(ts,lp,method='dop853_c')
o.animate(staticPlot=True); # remove the ; to display the animation
We can now increase the central density by reducing the size of the core. First, we will keep the core radius at a finite value, such that there is still a small harmonic core to the potential at \(R \ll R_c\), but from the vantage of the \(E=0\) orbits that we are considering, this region is almost negligible. For \(R_c = 0.05\), the surface of section is shown in Figure 13.24.
[70]:
from galpy.potential import LogarithmicHaloPotential, evaluatePotentials
from galpy.orbit import Orbit
from scipy import optimize
figure(figsize=(9,6),dpi=300)
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,y,vy) with given E"""
return Orbit([y,vy,numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)))
lp= LogarithmicHaloPotential(amp=1.,b=0.7,core=0.05)
E= 0.
# Location of closed orbits in [y,vy]
closed= {'banana': [0.1252,0.],
'box': [], # not closed!
'fish': [0.1741,1.099995],
'pretzel': [0.001,1.795],
'loop': [0.56,0.]}
nonclosed= {'banana': [[0.02,0.],[0.05,0.]],
'box': [[0.3,0.],[0.25,0.],[0.33,0.],[0.06,1.25],[0.35,0.],[0.39,0.]],
'fish': [[0.18,1.],[0.09,1.5]],
'pretzel': [[0.04,1.55]],
'loop': [[0.48,0.],[0.44,0.],[0.52,0.],[0.39,0.],[0.37,0.]]}
colors= {'banana': 'tab:orange',
'fish': 'tab:green',
'pretzel': 'tab:purple',
'box': 'tab:blue',
'loop': 'tab:cyan'}
for orb in closed:
for sgn in [-1,1]:
if not orb == 'box':
o= orbit_yvyE(sgn*closed[orb][0],closed[orb][1],E,pot=lp)
o.plotSOS(lp,ncross=1_000,method='dop853_c',surface='x',
gcf=True,overplot=(orb != 'banana') + (sgn != -1),
marker='.',mec='none',zorder=2,ms=3.,
mfc=colors[orb],rasterized=True,
xrange=[-0.75,0.75],yrange=[-2.6,2.6])
for nonc in nonclosed[orb]:
o= orbit_yvyE(sgn*nonc[0],nonc[1],E,pot=lp)
o.plotSOS(lp,ncross=2_000,method='dop853_c',surface='x',
overplot=True,rasterized=True,
marker='.',mec='none',zorder=2,ms=3.,
mfc=colors[orb])
# Also plot the zero-velocity curve, first find min x at which it exists
miny= optimize.brentq(lambda y: E-evaluatePotentials(lp,y,0.,phi=numpy.pi/2.),
0.01,0.9)-1e-10
ys= numpy.linspace(-miny,miny,201)
plot(ys,zvc(ys,E,pot=lp),color='k',zorder=0)
plot(ys,-zvc(ys,E,pot=lp),color='k',zorder=1)
galpy_plot.text(r'$b=0.7$'+'\n'+r'$R_c=0.05$',top_left=True,size=18.);
Figure 13.24: Like Figure 13.22, but for a smaller core radius \(R_c\).
The surface of section has now dramatically changed. Loop orbits and banana orbits take up a very large fraction of the section, with box orbits relegated to the small area between the loop and banana orbits. But in this area, two new minor orbit families appear that take up a significant fraction and squeeze the box orbits further. There are thus very few box orbits at this energy in this potential and the few that exist are quite extended in \(y\) and therefore not very suitable for sustaining the non-axisymmetric structure of the potential. To further understand the new minor orbit families, we look at the closed orbits in the surface of section that are the parents of non-loop orbits in Figure 13.25.
[71]:
figure(figsize=(8,3))
ts= numpy.linspace(0.,100.,100001)
labels= [r'$\mathrm{banana}$',
r'$\mathrm{fish}$',
r'$\mathrm{pretzel}$']
for ii,orb in enumerate(['banana','fish','pretzel']):
o= orbit_yvyE(closed[orb][0],closed[orb][1],E,pot=lp)
o.integrate(ts,lp,method='dop853_c')
subplot(1,3,ii+1)
o.plot(gcf=True,
xrange=[-1.,1.],
yrange=[-0.75,0.75])
plot([0.],[0.],'ko',ms=7.)
if ii > 0:
gca().get_yaxis().get_label().set_visible(False)
yticks([])
galpy_plot.text(labels[ii],top_left=True,size=24.)
tight_layout();
Figure 13.25: Closed boxlets in the surface of section of Figure 13.24.
We see that we again have a banana orbit that parents the large family of banana orbits near the center of the surface of section (the orange curves). The parent of the green curves in the surface of section is the second closed orbit, which looks like a fish and is therefore known as a fish orbit. The third closed orbit, the parent of the purple curves, twists around a few times and we’ll refer to it as a pretzel orbit. The remarkable property of these orbits is that they all avoid the center and are thus centrophobic; there are no closed orbits that pass through the center.
Finally, we can set the core radius to zero and look at orbits in the cuspy logarithmic potential. The surface of section is displayed in Figure 13.26 (see Miralda-Escude & Schwarzschild 1989).
[72]:
from galpy.potential import LogarithmicHaloPotential, evaluatePotentials
from galpy.orbit import Orbit
from scipy import optimize
figure(figsize=(9,6),dpi=300)
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,y,vy) with given E"""
return Orbit([y,vy,numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)))
lp= LogarithmicHaloPotential(amp=1.,b=0.7,core=0.)
E= 0.
# Location of closed orbits in [y,vy]
closed= {'banana': [0.12775,0.],
'fish': [0.16,1.15],
'pretzel': [0.34775,0.],
'53': [0.2905,0.],
'loop': [0.5575,0.]}
nonclosed= {'banana': [[0.03,0.],[0.07,0.],[0.0125,0.]],
'fish': [[0.05,2.],[0.12,1.2]],
'pretzel': [[0.343,0.],[0.105,1.8]],
'53': [[0.31,0.],],
'loop': [[0.5,0.],[0.4,0.],[0.36,0.]]}
colors= {'banana': 'tab:orange',
'fish': 'tab:green',
'pretzel': 'tab:purple',
'53': 'tab:blue',
'loop': 'tab:cyan'}
for orb in closed:
for sgn in [-1,1]:
o= orbit_yvyE(sgn*closed[orb][0],closed[orb][1],E,pot=lp)
o.plotSOS(lp,ncross=1_000,method='dop853_c',surface='x',
gcf=True,overplot=(orb != 'banana') + (sgn != -1),
marker='.',mec='none',zorder=2,ms=4.,
mfc=colors[orb],rasterized=True,
xrange=[-0.75,0.75],yrange=[-3.2,3.2])
for nonc in nonclosed[orb]:
o= orbit_yvyE(sgn*nonc[0],nonc[1],E,pot=lp)
o.plotSOS(lp,ncross=2_000,method='dop853_c',surface='x',
overplot=True,rasterized=True,
marker='.',mec='none',zorder=2,ms=4.,
mfc=colors[orb])
# Also plot the zero-velocity curve, first find min x at which it exists
miny= optimize.brentq(lambda y: E-evaluatePotentials(lp,y,0.,phi=numpy.pi/2.),
0.01,0.9)-1e-10
ys= numpy.linspace(-miny,miny,201)
plot(ys,zvc(ys,E,pot=lp),color='k',zorder=0)
plot(ys,-zvc(ys,E,pot=lp),color='k',zorder=1)
galpy_plot.text(r'$b=0.7$'+'\n'+r'$R_c=0.0$',top_left=True,size=18.);
Figure 13.26: Like Figure 13.22 and 13.24, but for a singular potential.
We see that box orbits have now completely disappeared and are instead replaced by different families of boxlets. All of these are regular orbits, although there is some mild chaos in this surface of section that is not shown here and has only a minor effect on the orbital structure. The closed orbits at the center of the boxlet orbit families are again centrophobic. They give rise to non-closed orbit families with qualitatively similar behavior to the closed orbits. All of these orbits are quite thick in the \(y\) direction and are therefore not very suitable for supporting the non-axisymmetric structure of the potential. Thus, we see that as the cored logarithmic potential becomes the singular, cuspy logarithmic potential, all box orbits cease to exist and that they are replaced by centrophobic, thick boxlets.
13.4.3. The shape evolution of dark matter halos¶
As discussed in the previous section, the \(1/r\) inner-density cusps of dark-matter halos are weak enough that supporting long-term-stable triaxial configurations is relatively straightforward. Thus, it is not surprising that dark-matter halos formed in simulations that only include dark matter and its collisionless evolution are able to remain triaxial for the age of the Universe (e.g., Frenk et al. 1988; Dubinski & Carlberg 1991; Jing & Suto 2002). Simulations that include the dissipationless physics of the gaseous component, as well as star formation and its associated feedback processes, form a central galaxy at the bottom of the dark-matter potential well. In such more realistic simulations, dark matter halos are generally rounder, that is, much closer to axisymmetric and gently oblate rather than strongly prolate (e.g., Dubinski 1994; Gustafsson et al. 2006; Abadi et al. 2010; Kazantzidis et al. 2010; Chua et al. 2019). Detailed analyses of the changes to halo orbits induced by the growth of the central galaxy demonstrate that the reason for these shape changes are undramatic changes to the shapes of the orbit, at least for realistic gaseous and stellar density profiles (e.g., Debattista et al. 2008; Valluri et al. 2010). Only if a very compact central galaxy forms does strong scattering significantly change the nature of the dark matter orbits in an irreversible manner (that is, chaotically).
Box orbits form the backbone of stable triaxial mass configurations and what happens to the box orbits when a central galaxy forms therefore determines how the triaxial figure evolves. To investigate this, we can see what happens when we grow a stellar disk inside of a triaxial dark-matter halo with typical values of the axis ratios \(b/a = 0.7\) and \(c/a = 0.5\) (triaxiality parameter \(T=0.68\); see Equation 12.7). We use values such that the dark matter halo and the disk are similar to the Milky Way’s, in units where the Sun is at \(R=1\) from the center; we use a Hernquist profile rather than a NFW profile such that the model has finite mass, but this does not affect the central structure. In Figure 13.27, we look at what happens to a typical narrow box orbit at \(z=0\) that extends out to the radius of the Sun.
[73]:
from galpy.potential import (TriaxialHernquistPotential,
MiyamotoNagaiPotential,
DehnenSmoothWrapperPotential)
from galpy.orbit import Orbit
hp= TriaxialHernquistPotential(amp=28.,a=2.,b=0.7,c=0.5)
mp= MiyamotoNagaiPotential(amp=0.6,a=1./3.,b=1./30.) # approx 1./10 mass within r=20.
o= Orbit([1.,0.78,0.1,0.])
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,hp)
figure(figsize=(9,3.5))
o.plot(label=r'$\mathrm{Before}$',gcf=True);
growts= numpy.linspace(0.,300.,1001) # approx. adiabatic growth
o.integrate(growts,hp+DehnenSmoothWrapperPotential(pot=mp,tform=0.,tsteady=growts[-1]))
o= o(growts[-1])
o.integrate(ts,hp+mp)
o.plot(label=r'$\mathrm{After\ adiabatic\ disk\ growth}$',overplot=True)
ylim(-0.6,0.6)
gca().set_aspect('equal')
legend(frameon=False,fontsize=15);
Figure 13.27: The evolution of a narrow box orbit in a triaxial Hernquist potential caused by the adiabatic growth of a disk.
We show the orbit before and after the adiabatic growth of the disk—adiabatic meaning that we grow the disk on a time scale that is much larger than the typical orbital time scales (see Chapter 3.4.3). There are two noticeable changes to the orbit. The first is that the orbit after the disk is grown does not extend to as large radii as before the disk growth. This is called adiabatic contraction (Blumenthal et al. 1986) and it occurs whenever a centrally-concentrated distribution of mass grows within a dark-matter halo (whether it is a disk, bulge, or elliptical galaxy). Adiabatic contraction causes the density of the dark matter to increase in the inner part of the halo, but it does not directly affect the shape, unless it also steepens the inner density cusp significantly away from \(1/r\), such that box orbits would be significantly deformed or destroyed. However, realistic central galaxies have a shallow enough density profile themselves that they do not significantly change the cusp, which remains \(\approx 1/r\) (Gnedin et al. 2004).
The second noticeable change between the box orbit before and after the adiabatic disk growth is that it fans out more after disk growth. All of the mass associated with this orbit would therefore be more spread out along the direction \(y\) along which the potential is flattened. This happens to all box orbits and the overall effect of the growth of the disk is therefore to make the triaxiality-supporting box orbits, and therefore the entire mass distribution, rounder. To see this more clearly, we can look at an illustrative surface of section. Considering orbits within the \(z=0\) plane with the same energy as the orbit above, the surface of section \((y,v_y)\) is shown in Figure 13.28.
[74]:
from galpy.potential import TriaxialHernquistPotential, evaluatePotentials
from galpy.orbit import Orbit
from scipy import optimize
figure(figsize=(9,6),dpi=300)
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,y,vy) with given E"""
return Orbit([y,vy,numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)))
hp= TriaxialHernquistPotential(amp=28.,a=2.,b=0.7,c=0.5)
E= -1.8
# Location of closed orbits in [y,vy]
closed= {'banana': [0.00655,0.],
'fish': [0.1127,0.49],
'pretzel': [0.33338,0.],
'pretzel2': [0.44355,0.],
'box': None,
'loop': [0.905,0.]}
nonclosed= {'banana': [[0.01,0.]],
'fish': [[0.14,0.5],[0.12,0.2]],
'pretzel': [[0.345,0.],[0.365,0.]],
'box': [[0.03,0.],[0.07,0.],[0.1,0.],[0.15,0.],[0.2,0.],
[0.25,0.],[0.275,0.],[0.29,0.],[0.4,0.],[0.42,0.]],
'pretzel2': [[0.44,0.]],
'loop': [[0.7,0.],[0.6,0.],[0.5,0.],[0.47,0.],[0.46,0.]]}
colors= {'banana': 'tab:orange',
'fish': 'tab:green',
'pretzel': 'tab:purple',
'pretzel2': 'tab:brown',
'box': 'tab:blue',
'loop': 'tab:cyan'}
for orb in closed:
for sgn in [-1,1]:
if not closed[orb] is None:
o= orbit_yvyE(sgn*closed[orb][0],closed[orb][1],E,pot=hp)
o.plotSOS(hp,ncross=100,method='dop853_c',surface='x',
gcf=True,overplot=(orb != 'banana') + (sgn != -1),
marker='.',mec='none',zorder=2,ms=4.,
mfc=colors[orb],rasterized=True,
xrange=[-1.35,1.35],yrange=[-2.3,2.3])
for nonc in nonclosed[orb]:
o= orbit_yvyE(sgn*nonc[0],nonc[1],E,pot=hp)
o.plotSOS(hp,ncross=1_200,method='dop853_c',surface='x',
overplot=True,rasterized=True,
marker='.',mec='none',zorder=2,ms=4.,
mfc=colors[orb])
# Also plot the zero-velocity curve, first find min x at which it exists
miny= optimize.brentq(lambda y: E-evaluatePotentials(hp,y,0.,phi=numpy.pi/2.),
0.1,1.9)-1e-10
ys= numpy.linspace(-miny,miny,201)
plot(ys,[zvc(y,E,pot=hp) for y in ys],color='k',zorder=0)
plot(ys,[-zvc(y,E,pot=hp) for y in ys],color='k',zorder=1);
Figure 13.28: A typical surface of section of the triaxial Hernquist potential from Figure 13.27.
This surface of section has a similar structure to the cored and singular logarithmic potentials that we considered in the previous section. There is a small family of banana orbits at the very center and other minor orbit families consisting of ‘fish’- and ‘pretzel’-like boxlets, but a substantial fraction of the surface of section is occupied by the regular box orbits (the blue curves above) that provide the main support for triaxiality (there are also many loop orbits). The closed orbits in the minor orbit families are displayed in Figure 13.29.
[75]:
figure(figsize=(12,3.5))
ts= numpy.linspace(0.,100.,100001)
for ii,orb in enumerate(['banana','fish','pretzel','pretzel2']):
o= orbit_yvyE(closed[orb][0],closed[orb][1],E,pot=hp)
o.integrate(ts,hp,method='dop853_c')
subplot(1,4,ii+1)
o.plot(gcf=True,
xrange=[-1.8,1.8],
yrange=[-1.3,1.3])
plot([0.],[0.],'ko',ms=7.)
if ii > 0:
gca().get_yaxis().get_label().set_visible(False)
yticks([])
tight_layout();
Figure 13.29: Closed boxlets in the surface of section of Figure 13.28.
Comparing them to those in the singular logarithmic potential, we find that they are similar, but substantially flatter in \(y\) and therefore better able to provide support for the non-axisymmetric figure.
To see what these orbits evolve to, we create a similar surface of section of the orbits after we have grown the disk. Note that the result is not a constant-energy surface of section any longer, because the energy is not conserved during the adiabatic growth. But it is useful to look at this to see what a population of orbits at a given energy evolves into. The surface of section is displayed in Figure 13.30.
[76]:
from galpy.potential import TriaxialHernquistPotential
from galpy.orbit import Orbit
from scipy import optimize
figure(figsize=(9,6),dpi=300)
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,y,vy) with given E"""
return Orbit([y,vy,numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,y,0.,phi=numpy.pi/2.)))
hp= TriaxialHernquistPotential(amp=28.,a=2.,b=0.7,c=0.5)
mp= MiyamotoNagaiPotential(amp=0.6,a=1./3.,b=1./30.) # approx 1./10 mass within r=20.
E= -1.8
growts= numpy.linspace(0.,300.,1001)
# Location of closed orbits in [y,vy]
closed= {'banana': [0.00655,0.],
'fish': [0.1127,0.49],
'pretzel': [0.33338,0.],
'pretzel2': [0.44355,0.],
'box': None,
'loop': [0.905,0.]}
nonclosed= {'banana': [[0.01,0.]],
'fish': [[0.14,0.5],[0.12,0.2]],
'pretzel': [[0.345,0.],[0.365,0.]],
'box': [[0.03,0.],[0.07,0.],[0.1,0.],[0.15,0.],[0.2,0.],
[0.25,0.],[0.275,0.],[0.29,0.],[0.4,0.],[0.42,0.]],
'pretzel2': [[0.44,0.]],
'loop': [[0.7,0.],[0.6,0.],[0.5,0.],[0.47,0.],[0.46,0.]]}
colors= {'banana': 'tab:orange',
'fish': 'tab:green',
'pretzel': 'tab:purple',
'pretzel2': 'tab:brown',
'box': 'tab:blue',
'loop': 'tab:cyan'}
for orb in closed:
for sgn in [-1,1]:
if not closed[orb] is None:
o= orbit_yvyE(sgn*closed[orb][0],closed[orb][1],E,pot=hp)
o.integrate(growts,hp+DehnenSmoothWrapperPotential(pot=mp,tform=0.,tsteady=growts[-1]),
method='dop853_c')
o= o(growts[-1])
o.plotSOS(hp+mp,ncross=100,method='dop853_c',surface='x',
gcf=True,overplot=(orb != 'banana') + (sgn != -1),
marker='.',mec='none',zorder=2,ms=4.,
mfc=colors[orb],rasterized=True,
xrange=[-1.35,1.35],yrange=[-2.3,2.3])
for nonc in nonclosed[orb]:
o= orbit_yvyE(sgn*nonc[0],nonc[1],E,pot=hp)
o.integrate(growts,hp+DehnenSmoothWrapperPotential(pot=mp,tform=0.,tsteady=growts[-1]),
method='dop853_c')
o= o(growts[-1])
o.plotSOS(hp+mp,ncross=1_200,method='dop853_c',surface='x',
overplot=True,rasterized=True,
marker='.',mec='none',zorder=2,ms=3.,
mfc=colors[orb]);
Figure 13.30: Evolution of the surface of section of Figure 13.28 caused by the adiabatic growth of a disk.
We have plotted the surface of section for the same axis limits to make it easily comparable to the original one. This surface of section confirms that the orbits by and large contract radially and expand in the \(y\) direction along which the mass distribution is flattened, because they pass \(x=0\) at smaller \(y\), but much higher velocity. The overall effect is thus to significantly round out the mass distribution. Indeed, it is clear both from the single orbit that we looked at above and the entire surface of section that the box orbits expand by a factor of \(\approx 1.5\) to \(2\) in the \(y\) direction. This is plenty to turn the \(b/a = 0.7\) flattening into a \(b/a \approx 1\) axisymmetric distribution. A similar effect in the \(z\) direction turns the \(c/a = 0.5\) into \(c/a \approx 0.75\) and the resulting shape is therefore an mildly oblate, axisymmetric dark matter halo, as seen in full \(N\)-body simulations.
In Figure 13.30, we have kept the color-coding of the different orbit families from the original surface of section. This makes it clear that some orbits change family due to the adiabatic growth of the disk. A few examples of this are shown in Figure 13.31.
[77]:
def plot_orbit_change(y,vy):
o= orbit_yvyE(y,vy,E,pot=hp)
o.integrate(ts,hp,method='dop853_c')
o.plot(gcf=True,xrange=[-1.49,1.49],yrange=[-1.49,1.49]);
o= orbit_yvyE(y,vy,E,pot=hp)
o.integrate(growts,hp+DehnenSmoothWrapperPotential(pot=mp,tform=0.,tsteady=growts[-1]),
method='dop853_c')
o= o(growts[-1])
o.integrate(ts,hp+mp,method='dop853_c')
o.plot(overplot=True)
gca().set_aspect('equal')
figure(figsize=(12,3.5))
subplot(1,4,1)
plot_orbit_change(nonclosed['banana'][0][0],nonclosed['banana'][0][1])
subplot(1,4,2)
plot_orbit_change(nonclosed['fish'][-1][0],nonclosed['fish'][-1][1])
gca().get_yaxis().get_label().set_visible(False)
yticks([])
subplot(1,4,3)
plot_orbit_change(nonclosed['box'][-1][0],nonclosed['box'][-1][1])
gca().get_yaxis().get_label().set_visible(False)
yticks([])
subplot(1,4,4)
plot_orbit_change(nonclosed['pretzel2'][-1][0],nonclosed['pretzel2'][-1][1])
gca().get_yaxis().get_label().set_visible(False)
yticks([])
tight_layout();
Figure 13.31: Evolution of selected boxlet orbits caused by the adiabatic growth of a disk.
This figure demonstrates a few important aspects of the family changes. Boxlet families are largely turned into regular box orbits: we see how a banana orbit and a fish orbit turn into regular box orbits. However, this does not significantly change the overall shape of the spatial extent of the original orbit. Similarly, we see how a box orbit and a pretzel orbit become loop orbits. But again, while the orbits contract, the overall shape of their spatial extent is similar. Therefore, the effect of orbit-family changes on the shape of the dark matter halo is limited. Shape changes due to the growth of baryonic disks are therefore largely caused by the fact that they make narrow box orbits rounder rather than by significant changes to the types of orbits that are populated.
The discussion here has focused on the effect of the adiabatic growth of a central galactic disk. But in addition to deepening the gravitational potential well, gas cooling and star formation is also accompanied by energetic feedback processes from stellar winds and supernova explosions for galaxies like the Milky Way and smaller and by feedback from accretion of gas onto the central supermassive black hole in larger galaxies (see Chapter 19.3.3) . These processes can transfer energy and momentum to the dark matter halo and this happens in a largely non-adiabatic manner (that is, fast compared to the orbital time scale). This transfer of, e.g., energy will further effect the orbits of halo particles. The adiabatic contraction that we witnessed in this section will be particularly affected by these processes, because increasing the energy of the dark matter particles will tend to enlarge their orbits and reduce the central density (we discussed this in the context of dwarf galaxies in Chapter 6.8). It is therefore likely that there are no lasting effects from adiabatic contraction in galaxies that we observe today. Feedback processes, however, do not strongly affect the shape further and certainly cannot reverse the trend towards more axisymmetric and more spherical shapes, because energy transfer from the baryons to the dark matter does not flatten or further populate the triaxiality-supporting box orbits. It is therefore a robust prediction that the shapes of the inner density profiles of dark matter halos with massive galaxies at their centers should be axisymmetric, oblate ellipsoids with a flattening that in detail depends on the relative mass and radial profile of the disk (e.g., Kazantzidis et al. 2010).
.
.
.