9.2. The separability of disk orbits¶
The orbits that we looked at above occupy different regions of phase-space, even though they have the same \((E,L_z)\), with the difference being due to their different initial ratio of \(v_R/v_z\). If you play around with changing the \((E,L_z)\) in the examples above, you will see that this is a general property in this Milky-Way-like axisymmetric potential and numerical integrations in a wide range of axisymmetric galactic potentials demonstrate that this is generally the case (Ollongren 1962).
In a completely general, time-dependent potential with no symmetries, orbits would cover a fully six-dimensional region of phase-space if integrated for a long enough time. As we discussed in Chapter 4.3, if a potential has symmetries, Noether’s theorem implies that orbits in this potential satisfy integrals of the motion. For example, as we have already discussed a few times, if the potential does not explicitly depend on time, then the energy of an orbit is conserved (Chapter 3.1). Similarly, when the potential is spherical, the full angular-momentum vector is conserved (as we demonstrated explicitly in Chapter 4.1) and when it is axisymmetric, the component of the angular momentum parallel to the symmetry axis is conserved (see above). As we discussed in Chapter 4.3, one aspect of these integrals that is especially important is that they are all isolating integrals: they restrict the volume of phase-space occupied by an orbit to a lower-dimensional subspace.
If we have identified all independent integrals of the motion in a given potential, then we expect to see orbits fill the entire region allowed by these integrals. For example, for general spherical potentials the conservation of energy and the magnitude of the angular momentum imply that orbits are restricted to a radial range \(r_p \leq r \leq r_a\), with the zero-velocity curve in this case given by the two parts: \(r = r_p\) and \(r = r_a\) (“zero velocity” in this case means zero radial velocity; the rotational velocity is always non-zero if the angular momentum is non-zero). Unless the potential is that of a point-mass or a homogeneous sphere (which have an additional isolating integral of motion), the orbit will occupy the entire region \(r_p \leq r \leq r_a\). For example, returning to the example orbit in an isochrone potential with \(b=1\) from Figure 4.8 and integrating it for thousands of dynamical times, we see in Figure 9.5 that it eventually covers the entire region between \(r_p \leq r \leq r_a\).
[12]:
from galpy import potential
from galpy.orbit import Orbit
ip= potential.IsochronePotential(normalize=1.,b=1.)
oi= Orbit([1.,0.25,1.4,0.])
ts= numpy.linspace(0.,10000.,1000001)
oi.integrate(ts,ip)
figure(figsize=(6,4))
oi.plot(xrange=[-2.55,2.55],yrange=[-2.55,2.55],
ls='None',marker='.',gcf=True,
ms=.25,alpha=0.05,rasterized=True)
gca().set_aspect('equal');
Figure 9.5: A long-term integration of a disk orbit demonstrating where it spends most time.
where we have used transparency to make it clear where the orbit spends most of its time.
The fact that orbits in axisymmetric potentials do not fill the entire region of space allowed by their \((E,L_z)\) indicates that there may be another, third integral beyond \((E,L_z)\) that restricts the motion. It is unclear what this additional independent integral would be, because the potential does not have any symmetries beyond time independence and axisymmetry. Therefore, we cannot use Noether’s theorem to identify another integral of the motion.
We can gain some insight into what is going on by looking at how well orbits are represented by the product of their motion in the \(z=0\) plane and that perpendicular to this plane. For the axisymmetric orbit with \((E,L_z,v_R) = (0.8,-1.25,0.35)\) that we looked at in Figure 9.3, we can plot the projection of the orbit in the \(z=0\) plane and compare it to the orbit that we would get if we integrated the orbit with the same initial \((R,v_R)\) in the \(z=0\) gravitational potential (that is, we completely ignore the vertical motion). We look at the orbit in the \((x,y)\) plane and as \(R\) versus \(t\) in Figure 9.6.
[13]:
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.)
figure(figsize=(10,4.3))
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,potential.MWPotential2014)
subplot(1,2,1)
line_full= o.plot(d1='x',d2='y',zorder=1,
label=r'$\mathrm{Full\ 3D}$',
xrange=[-10,12.5],yrange=[-10,12.5],
gcf=True)
subplot(1,2,2)
line_full= o.plotR(yrange=[0.,8.],zorder=1,
label=r'$\mathrm{Full\ 3D}$',
gcf=True)
# Only consider (x,y) part of the orbit
op= o.toPlanar()
op.integrate(ts,potential.MWPotential2014)
subplot(1,2,1)
line_planar= op.plot(d1='x',d2='y',overplot=True,zorder=0,
label=r'$\mathrm{Planar\ (2D)}$')
gca().set_aspect('equal')
legend(handles=[line_full[0],line_planar[0]],
fontsize=18.,loc='upper right',frameon=False)
subplot(1,2,2)
line_planar= op.plotR(overplot=True,zorder=0,
label=r'$\mathrm{Planar\ (2D)}$')
legend(handles=[line_full[0],line_planar[0]],
fontsize=18.,loc='lower right',frameon=False);
Figure 9.6: Approximating a disk orbit’s planar components in the mid-plane potential.
We see that we get an excellent approximation to the full orbit’s path in \((x,y)\) by ignoring the vertical part of the orbit. Comparing the time dependence of \(R(t)\) shows that the planar approximation provides a great match to the radial range that is covered by the full orbit. The oscillations in \(R(t)\) are similar in the full orbit and the 2D approximation, but the oscillations start to go out of phase after about a few dozen oscillations. This means that the frequency of the oscillation is not exactly matched in the approximation. But overall we would derive very similar properties about the orbit’s projection on the \(z=0\) plane from approximating it as never leaving the plane. This is remarkable, given that this particular orbit has a rather large eccentricity \(e \approx 0.35\). Most stars, for example, in the Milky-Way disk have smaller eccentricities than this.
We can similarly look at the projection of the orbit on the \((z,v_z)\) plane in phase space and compare it to integrating the same \((z,v_z)\) initial conditions in the fixed one dimensional potential defined by evaluating \(\Phi(R,z)\) at constant \(R\). To get a better sense of how the orbit covers the \((z,v_z)\) plane, we integrate the full orbit for ten times longer than before and display it in the left panel of Figure 9.7.
[14]:
from galpy import potential
from galpy.orbit import Orbit
figure(figsize=(12,6))
subplot(1,2,1)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,10001)
shortts= numpy.linspace(0.,20.,1001)
o.integrate(ts,potential.MWPotential2014)
o.plot(d1='z',d2='vz',zorder=0,lw=0.6,
xrange=[-0.4,0.4],yrange=[-45,45],
gcf=True)
# Only consider orbit in z
ol= o.toLinear()
ol.integrate(shortts,potential.RZToverticalPotential(\
potential.MWPotential2014,R))
line_olmax= ol.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{max}\, R(t)$')
olmin= o(ts[numpy.argmin(o.R(ts))]).toLinear()
olmin.integrate(shortts,potential.RZToverticalPotential(\
potential.MWPotential2014,
numpy.amin(o.R(ts))))
line_olmin= olmin.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{min}\, R(t)$')
olmed= o(ts[o.R(ts) == numpy.median(o.R(ts))]).toLinear()
olmed.integrate(shortts,potential.RZToverticalPotential(\
potential.MWPotential2014,
numpy.median(o.R(ts))))
line_olmed= olmed.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{med}\, R(t)$')
legend(handles=[line_olmin[0],line_olmed[0],
line_olmax[0]],
fontsize=18.,loc='center',frameon=False)
subplot(1,2,2)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
# 1D little faster if we initialize off of dt=20
ts= numpy.linspace(0.,20.,1001)
o.integrate(ts,potential.MWPotential2014)
ol= o.toLinear()
ol.integrate(shortts,potential.RZToverticalPotential(\
potential.MWPotential2014,R))
olmin= o(ts[numpy.argmin(o.R(ts))]).toLinear()
olmin.integrate(shortts,potential.RZToverticalPotential(\
potential.MWPotential2014,
numpy.amin(o.R(ts))))
olmed= o(ts[o.R(ts) == numpy.median(o.R(ts))]).toLinear()
olmed.integrate(shortts,potential.RZToverticalPotential(\
potential.MWPotential2014,
numpy.median(o.R(ts))))
ts= numpy.linspace(0.,100.,10001)
o.integrate(ts,potential.MWPotential2014)
o.plot(d1='z',d2='vz',zorder=0,lw=0.8,
xrange=[-1.74,1.74],yrange=[-130,130.],
gcf=True)
line_olmin= olmin.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{min}\, R(t)$')
line_olmax= ol.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{max}\, R(t)$')
line_olmed= olmed.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{med}\, R(t)$')
legend(handles=[line_olmin[0],line_olmed[0],
line_olmax[0]],
fontsize=18.,loc='center',frameon=False)
gca().set_ylabel('')
tight_layout(w_pad=0.);
Figure 9.7: Approximating a disk orbit’s vertical components in a one-dimensional, vertical potential at \(R\).
The full orbit’s projection is shown in blue. We show three one-dimensional approximations: one initialized at a point along the full orbit that is at the minimum \(R\) that it reaches, one at a point at the median value of \(R\), and one at the maximum \(R\). We see that the 1D approximation to the orbit started at the median \(R\) does a good job of delineating the typical path of the full orbit and the approximations at the minimum and maximum \(R\) capture the extremes of the region well. Nevertheless, the one-dimensional approximation for the \(z(t)\) motion is much less successful than the two-dimensional approximation of the planar part of the motion. This is partly because the orbit has such high eccentricity: the orbit goes through a wide range of \(R\) where the vertical potential is quite different.
This behavior holds more generally. If we look at the first axisymmetric orbit that we considered (the one with \([E,L_z,v_R] = [-1.250,0.6,0.0]\) from Figure 9.2), we get a similar behavior of the one-dimensional approximations, shown in the right panel of Figure 9.7. It is quite remarkable that we can still delineate the region occupied by the full orbit in the \((z,v_z)\) plane even though this orbit travels up to 1.5 kpc from the plane!
To get a better sense for how this works for more typical disk orbits, which have smaller eccentricities than the orbits that we have looked at so far, we look at the orbit with \((E,L_z,R,v_R) = (-1.385,0.6,0.65,0.02)\). This is orbit is much closer to a circular orbit and has an eccentricity of \(\approx 0.1\). The projection of this orbit onto the \(z=0\) plane and the \(R(t)\) dependence are shown in the left and middle panels of Figure 9.8.
[17]:
from matplotlib import gridspec
from galpy import potential
from galpy.orbit import Orbit
figure(figsize=(13,4.3))
gs = gridspec.GridSpec(1,3,width_ratios=[1.25,1.,1.15])
R, E, Lz= 0.65, -1.385, 0.6
vR= 0.02
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,5001)
o.integrate(ts,potential.MWPotential2014)
subplot(gs[0])
line_full= o.plot(d1='x',d2='y',zorder=1,lw=0.7,
label=r'$\mathrm{Full\ 3D}$',
xrange=[-7.,9],yrange=[-7.,9.],
gcf=True)
subplot(gs[1])
line_full= o.plotR(yrange=[0.,7.],zorder=1,
label=r'$\mathrm{Full\ 3D}$',
gcf=True)
op= o.toPlanar()
op.integrate(ts,potential.MWPotential2014)
subplot(gs[0])
line_planar= op.plot(d1='x',d2='y',overplot=True,zorder=0,
lw=0.7,label=r'$\mathrm{Planar\ (2D)}$')
gca().set_aspect('equal')
legend(handles=[line_full[0],line_planar[0]],
fontsize=18.,loc='upper right',frameon=False)
subplot(gs[1])
line_planar= op.plotR(overplot=True,zorder=0,
label=r'$\mathrm{Planar\ (2D)}$')
legend(handles=[line_full[0],line_planar[0]],
fontsize=18.,loc='lower right',frameon=False)
subplot(gs[2])
line_full= o.plot(d1='z',d2='vz',zorder=0,lw=0.7,
xrange=[-0.3,0.3],yrange=[-35,35.],
gcf=True,
label=r'$\mathrm{Full\ 3D}$')
ts= numpy.linspace(0.,2.,1001)
olmed= o(ts[o.R(ts) == numpy.median(o.R(ts))]).toLinear()
olmed.integrate(ts,potential.RZToverticalPotential(\
potential.MWPotential2014,
numpy.median(o.R(ts))))
line_med= olmed.plot(overplot=True,zorder=1,lw=2.,
label=r'$R = \mathrm{med}\, R(t)$')
legend(handles=[line_full[0],line_med[0]],
fontsize=18.,loc='center',frameon=False)
tight_layout(w_pad=0.5);
Figure 9.8: Like Figures 9.6 and 9.7, but for an orbit with typical eccentricity and vertical excursions.
We again compare the full orbit to the approximation when ignoring the vertical motion. This approximation works similar as before: The radial range covered by the orbit is again well approximated by the 2D approximation, while the oscillation frequency is slightly different, causing the oscillations to drift out of phase as time goes on. Because the orbit covers a smaller range of radii \(R\), the projection of the motion onto the \((z,v_z)\) plane is much better approximated by the one-dimensional approximation at the median \(R\) as shown in the right panel of Figure 9.8.
All of these numerical investigations lead us to conclude that a good description of a disk orbit that is close to circular in an axisymmetric potential—a good approximation to the actual dynamical situation of stars in galactic disks—can be obtained by treating the in-plane motion \((R,\phi,v_R,v_\phi)\) and the perpendicular-to-plane motion in \((z,v_z)\) separately. That is, a good approximation of the orbits can be obtained by splitting the overall potential \(\Phi(R,z)\) into a radial \(\Phi_R(R) = \Phi(R,z=0)\) and a vertical \(\Phi_z(z;\tilde{R}) = \Phi(\tilde{R},z)\) part (where the vertical part is defined at some \(\tilde{R}\) along the orbit): \begin{equation}\label{eq-axipot-sep} \Phi(R,z) = \Phi_R(R) + \Phi_z(z;\tilde{R})\,. \end{equation} When this is the case, the effective Hamiltonian in Equation (9.3) further simplifies to \begin{align} H_\mathrm{eff}(R,z,p_R,p_z;L_z) & = \frac{1}{2}\,\left(p_R^2 +p_z^2\right)+\Phi(R,z)+\frac{L_z^2}{2R^2}\nonumber\\ & = \underbrace{\frac{p_R^2}{2}+\Phi_R(R)+\frac{L_z^2}{2R^2}}_{H_{R,\mathrm{eff}}(R,p_R;L_z)} + \underbrace{\frac{p_z^2}{2}+\Phi_z(z;\tilde{R})}_{H_z(z,p_z)}\label{eq-hamiltonian-effective-sep-cyl}\,. \end{align} Therefore, in this case, the Hamiltonian splits into two independent parts as \(H_\mathrm{eff}(R,z,p_R,p_z;L_z) = H_{R,\mathrm{eff}}(R,p_R;L_z) + H_z(z,p_z)\), which will therefore have two independent motions.
What this means in terms of integrals of the motion is the following: when the Hamiltonian splits in the manner of Equation (9.9), then the specific energy associated with the radial motion: \(E_R = v_R^2/2 + v_\phi^2/2+\Phi_R(R)\) and that associated with the vertical motion \(E_z = v_z^2 + \Phi(z;\tilde{R})\) are each conserved on their own. Thus, in this case we have three integrals of the motion in this region of phase-space near circular orbits of the axisymmetric potential: \(L_z\), \(E_R\), and \(E_z\). These three integrals of the motion act to restrict the orbit to a three-dimensional subspace of phase-space.
We will return to the separability of disk orbits in Section 9.4.3 below, where we will demonstrate that disk orbits often separate even better when considering them in prolate spheroidal coordinates. We will use the separability of disk orbits into planar and vertical components to build equilibrium models of galactic disks in Chapter 10. Finally, we will come back to the question of whether all orbits in an axisymmetric potential have a third integral in Chapter 13.3.
.
.
.