13.1. Surfaces of section¶
Orbits in three-dimensional potentials are trajectories in six-dimensional phase-space. We have so far only discussed orbits in spherical potentials and axisymmetric potentials, where the symmetry of the potential allows one to get a reasonable projection of the orbit when viewed in just two of the six phase-space coordinates. However, when we consider orbits in triaxial potentials or when we want to more closely inspect the orbital structure of spherical and axisymmetric potentials, a better visualization of the orbits is in order. Here we discuss a standard method for visualizing the orbital structure of gravitational potentials and use it first to reconsider the orbits in disk potentials that we first discussed in Chapter 9. As in that chapter, we consider axisymmetric potentials that are symmetric with respect to rotation around the \(z\) axis and mirror-symmetric around \(z=0\).
Orbits in static axisymmetric potentials have two integrals of the motion: the energy \(E\) and the \(z\) component \(L_z\) of the angular momentum. As discussed in Chapter 9.1 (see also Chapter 4.3), the motion in \((\phi,v_\phi)\) is directly determined from \((L_z,R[t])\) through the conservation of \(L_z\) and we can thus fully investigate the orbital structure by considering only the four-dimensional phase-space \((R,z,v_R,v_z)\) in the meridional plane. In this four-dimensional space, the energy integral confines the orbit to a three-dimensional volume. We can choose to visualize this volume in \((R,z,v_R)\), because the energy fixes \(v_z\) given \((R,z,v_R)\) (the energy only fixes \(v_z\) up to a sign and we thus pick one of the possible signs, say the positive one). We can easily visualize this three-dimensional volume, by taking a slice through it at \(z=0\); this is called a surface of section. It is well-defined for an axisymmetric potential that is mirror-symmetric around \(z=0\), because every orbit passes through \(z=0\). A two-dimensional slice through a three-dimensional volume generically fills a two-dimensional area and this is what we expect to see.
As an example, we consider the same orbit that we displayed in Figure 9.2 in Chapter 9.1. To start, we take the simple function to compute the surface of section defined in Bovy (2015), which simply searches for \(z=0\) crossings with \(v_z>0\) along an integrated orbit. The surface of section of this orbit is shown in the left panel of Figure 13.1.
[6]:
from galpy import potential
from galpy.orbit import Orbit
def orbit_RvRELz(R,vR,E,Lz,pot=None):
"""Returns Orbit at (R,vR,phi=0,z=0) with given (E,Lz)"""
return Orbit([R,vR,Lz/R,0.,
numpy.sqrt(2.*(E-potential.evaluatePotentials(pot,R,0.)
-(Lz/R)**2./2.-vR**2./2)),0.],ro=8.,vo=220.)
def surface_section(Rs,zs,vRs):
# Find points where the orbit crosses z from - to +
shiftzs= numpy.roll(zs,-1)
indx= (zs[:-1] < 0.)*(shiftzs[:-1] > 0.)
return (Rs[:-1][indx],vRs[:-1][indx])
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,10001)
o.integrate(ts,potential.MWPotential2014)
sect1Rs,sect1vRs=surface_section(o.R(ts),o.z(ts),o.vR(ts))
figure(figsize=(10,3.5))
subplot(1,2,1)
plot(sect1Rs,sect1vRs,'o',mec='none')
xlabel(r'$R\,(\mathrm{kpc})$')
ylabel(r'$v_R\,(\mathrm{km\,s}^{-1})$')
subplot(1,2,2)
ts= numpy.linspace(0.,3000.,300001)
o.integrate(ts,potential.MWPotential2014)
sect1Rs,sect1vRs=surface_section(o.R(ts),o.z(ts),o.vR(ts))
plot(sect1Rs,sect1vRs,'o',mec='none')
xlabel(r'$R\,(\mathrm{kpc})$')
ylabel(r'$v_R\,(\mathrm{km\,s}^{-1})$')
tight_layout(w_pad=0.4);
Figure 13.1: A coarse (left) and well-sampled (right) surface of section.
We see that the orbit does not appear to fill a two-dimensional area, but is instead confined to a one-dimensional curve. Even if we integrate the orbit for many more dynamical times this observation does not change, as the well-sampled surface of section in the right panel of Figure 13.1 demonstrates. Thus, the locus of the orbit in the surface of section is one-dimensional, which means that the locus of the orbit in the three-dimensional \((R,z,v_R)\) space is only two-dimensional, a shell. The motion must therefore be confined to this two-dimensional subspace by an integral of the motion in addition to the energy and the \(z\) component of the angular momentum. The surface of section offers more direct support for this than simply looking at the orbit in the \((R,z)\) plane as we did in Chapter 9.1. This third integral is non-classical, because it does not follow from a symmetry of the potential, like for the energy (from time invariance) and \(L_z\) (from axisymmetry).
While the method we have used so far to compute the surface of section is conceptually simple, it does not precisely determine the location of crossings and it can be slow when we investigate more complicated orbits. From now on, we will instead use a more sophisticated way to compute surfaces of section that consists of re-writing the orbit integration in terms of a new independent variable \(\psi\) defined through (Hunter et al. 1998) \begin{align}\label{eq-triaxialorbits-sos-psi} z & = A\,\sin \psi\,;\quad v_z = A\,\cos \psi\,, \end{align} and similar when we use \(x=0\) or \(y=0\) surfaces of section for two-dimensional potentials below. Crossings then occur at integer multiples of \(\psi = 2\pi\) and we can exactly determine these by integrating with \(\psi\) as the independent variable (time then becomes a dependent variable). We leave the derivation of this method and of the conditions under which it works as an exercise, but suffice it to say that it works well for all of the examples in this chapter.
To visualize all of the possible orbits with \(E\) and \(L_z\) of the orbit above, we generate surfaces of section for six orbits ranging from having all available initial kinetic energy in the \(v_R\) component to having this all in the \(v_z\) component. This generates the sequence of curves displayed in Figure 13.2.
[7]:
def maxvR_RELz(R,E,Lz,pot=None):
"""Returns maximum v_R for orbit at (R,phi=0,z=0) with given (E,Lz)"""
return numpy.sqrt(2.*(E-potential.evaluatePotentials(pot,R,0.))-(Lz/R)**2.)
R, E, Lz= 0.8, -1.25, 0.6
# Avoid maxvR orbit that is exactly in the z=0 plane
vRs= numpy.sqrt(numpy.linspace(0.,(maxvR_RELz(R,E,Lz,pot=potential.MWPotential2014)-1e-6)**2.,7))
o= Orbit([orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
for vR in vRs])
figure(figsize=(6,3.5))
o.plotSOS(potential.MWPotential2014,ms=3.,ncross=1_000,gcf=True);
Figure 13.2: Surfaces of section of constant \(E\) and \(L_z\).
The innermost contour, which has all initial available kinetic energy in \(v_z\), is the same as that shown above. The other contours have more and more initial kinetic energy in \(v_R\). All of the orbits at this \((E,L_z)\) are confined between the innermost and outermost curves shown here and we see that all orbits form a closed, one-dimensional curve in the surface of section.
Another way to look at the orbital structure is to consider orbits at constant specific energy \(E\), but that vary \(L_z\) between a value close to zero and the maximum value for this \(E\) (which happens at the circular orbit with this \(E\)). To display only a single orbit at each \(L_z\) among the continuum of orbits at \((E,L_z)\), we always choose the orbit that has half of its initially-available kinetic energy in the \(v_z\) component and half in \(v_R\) (the middle orbit in Figure 13.2). An example of this is shown in the left panel of Figure 13.3.
[8]:
from scipy import optimize
def maxLz_E(E,pot=None):
"""Returns (R,Lz) for the circular orbit with energy E"""
rE= optimize.brentq(lambda R: potential.evaluatePotentials(pot,R,0.)
+potential.vcirc(pot,R)**2./2.-E,0.01,2.)
return (rE,rE*potential.vcirc(pot,rE))
E= -1.25
Rc, maxLz= maxLz_E(E,pot=potential.MWPotential2014)
Lzs= numpy.linspace(0.05,maxLz-0.000001,7) #v. close to circular
figure(figsize=(10,3.5))
subplot(1,2,1)
for Lz in Lzs:
o= orbit_RvRELz(Rc,
maxvR_RELz(Rc,E,Lz,pot=potential.MWPotential2014)\
/numpy.sqrt(2.),
E,Lz,pot=potential.MWPotential2014)
if Lz == Lzs[-1]:
marker= 'o'
else:
marker= '.'
o.plotSOS(potential.MWPotential2014,gcf=True,
marker=marker,overplot=Lz != Lzs[0])
subplot(1,2,2)
E= -1.25
Rc, maxLz= maxLz_E(E,pot=potential.MWPotential2014)
# Previous Lzs, we want to zoom in around the third from the end
Lzs= numpy.linspace(0.05,maxLz-0.000001,7) #v. close to circular
Lzs= numpy.linspace(Lzs[-3]-0.03,Lzs[-3]+0.03,5)
o= Orbit([orbit_RvRELz(Rc,maxvR_RELz(Rc,E,Lz,pot=potential.MWPotential2014)\
/numpy.sqrt(2.),
E,Lz,pot=potential.MWPotential2014)
for Lz in Lzs])
o.plotSOS(potential.MWPotential2014,ncross=1_000,ms=3.,gcf=True)
tight_layout(w_pad=0.4);
Figure 13.3: Surfaces of section of constant \(E\) with \(L_z\) ranging from zero to that of a circular orbit. The right panel shows a closer sampling near the resonant orbit from the left panel.
We can learn a lot about the orbits in the axisymmetric Milky-Way potential MWPotential2014 used to create this diagram. The innermost orbit has the maximum \(L_z\) for this energy and this is a circular orbit. The surface of section therefore reduces to a point. The other orbits generate one-dimensional curves that typically loop around this circular orbit; we say that the circular orbit is the parent orbit
of this family of orbits. The loops get bigger and bigger as \(L_z\) decreases, with the lowest \(L_z\) orbit (the outermost curve) getting quite close to the center of the potential. The bulge component in MWPotential2014 ranges out to \(r = 1.9\,\mathrm{kpc}\); the lowest \(L_z\) orbit enters the bulge region where the potential is quite different from that in the disk and halo region (where the other orbits in this diagram orbit) which leads to the kink toward smaller
\(|v_R|\) at \(R \lesssim 1.9\,\mathrm{kpc}\).
All but one of the orbits in the left panel of Figure 13.3 loop around the parent orbit in a regular fashion. The exception is the third orbit starting from the center. This orbit only covers part of a loop. Such partial loops in a surface of section are indicative of a resonant orbit: an orbit whose fundamental orbital frequencies are commensurate, such that these orbits cannot cover their entire allowed phase-space volume (even after taking into account the third integral). We take a closer look at this orbit in Figure 13.4.
[9]:
ts= numpy.linspace(0.,100.,10001)
Lz= Lzs[-3]
o= orbit_RvRELz(Rc,maxvR_RELz(Rc,E,Lz,pot=potential.MWPotential2014)\
/numpy.sqrt(2.),
E,Lz,pot=potential.MWPotential2014)
o.integrate(ts,potential.MWPotential2014)
figure(figsize=(8,4))
subplot(1,2,1)
o.plot(gcf=True)
subplot(1,2,2)
o.plot(d1='vR',d2='vz',gcf=True)
tight_layout();
Figure 13.4: A resonant orbit.
We see that this orbit is indeed resonant. Unlike the orbits that we displayed in Chapter 9.1, there is a clear correlation between the radial and vertical components of the orbit, due to the commensurability of the radial and vertical frequencies. Another way to display this is in the projection of the orbit in the \((v_R,v_z)\) plane, shown on the right.
Such resonant orbits are a common feature of galactic potentials and they occur more frequently and with much more variety in triaxial potentials. Resonances can make phase-space complicated in their surroundings, but this particular resonant orbit does not “affect” too much of phase space in this way. To demonstrate this, we zoom in around this orbit the right panel of Figure 13.3, which displays two orbits with slightly higher \(L_z\) and two orbits with slightly lower \(L_z\). As is clear, these are regular, non-resonant orbits that loop around the parent circular orbit.
Thus, orbital surfaces of section are highly useful diagnostics of the orbital structure of galactic potentials. For the axisymmetric potentials, they clearly demonstrate that most orbits must have a non-classical third integral of motion and they provide a simple manner to find resonant orbits.
.
.
.