9.4. Action-angle coordinates in and around disks¶
9.4.1. General considerations¶
We introduced action-angle coordinates as a particularly-useful set of canonical coordinates in general in Chapter 3.4.3 and demonstrated in Chapter 4.4 that we can compute them for all spherical potentials. Action-angle coordinates are in many ways even more useful for orbits in and around disks, because of the general difficulty of determining the third integral of motion for disk orbits. Because there are three conserved action coordinates, if we can compute the three actions, we have computed a third integral. Actions are for this reason very useful for labeling disk orbits and for building equilibrium models of galactic disks, because realistic disk models must depend on the third integral as we will discuss at the start of Chapter 10. The combination of advances in computing action-angle coordinates for disky potentials discussed below and the explosion of data on Milky-Way disk orbits from the Gaia survey, has meant that action-angle coordinates have become a very popular tool in Galactic dynamics and archaeology in recent years (e.g., Trick et al. 2019; Vasiliev 2019; Naidu et al. 2020).
Computing action-angle coordinates in disky potentials—potentials that are overall disk-shaped due to the presence of a massive disk, but can include a bulge and halo, such as the MWPotential2014 model used above—is complicated, because as soon as we leave spherical symmetry behind, there is no general way to solve the Hamilton-Jacobi equation and obtain \((\vec{J},\boldsymbol{\theta})\) for a given
single phase-space position \((\vec{x},\vec{v})\). The only practical method for solving a partial differential equation such as the Hamilton-Jacobi equation (3.63) is using the separation-of-variables technique, but for general axisymmetric potentials the Hamilton-Jacobi equation cannot be separated. We are, therefore, limited to specific gravitational potentials for which the Hamilton-Jacobi equation separates, or we need approximate the
Hamiltonian using a separable version and we are then restricted to using approximate action-angle coordinates. As we will discuss below, practical approximations allow actions to be computed that are precise enough to be useful. A somewhat different method for obtaining actions and angles is to avoid solving the Hamilton-Jacobi equation altogether, and instead directly solve for the generating function that maps \((\vec{x},\vec{v}) \rightarrow (\vec{J},\boldsymbol{\theta})\). As we will
discuss below, the generating function can be determined using an orbit integration. Because numerically integrating an orbit is a somewhat slow process, this is therefore a laborious way of obtaining action-angle coordinates, but it has the advantage that the actions and angles are highly precise.
We start by discussing techniques that solve the Hamilton-Jacobi equation by approximating the Hamiltonian as being separable. As in Chapter 10, we set the mass equal to one, because a body’s mass does not (non-trivially) affect its action-angle coordinates in a fixed potential. Writing out the Hamilton-Jacobi equation (3.63) in cylindrical coordinates, we find \begin{equation}\label{eq-diskorbit-HJ-Econserved-cylindrical} \left(\frac{\partial W}{\partial R}\right)^2 +{1 \over R^2}\,\left(\frac{\partial W}{\partial \phi}\right)^2 +\left(\frac{\partial W}{\partial z}\right)^2 +2\Phi(R,\phi,z) = 2E\,. \end{equation} For any axisymmetric potential, we can write this as \begin{equation}\label{eq-diskorbit-HJ-Econserved-cylindrical-sepphi} R^2\,\left(\frac{\partial W}{\partial R}\right)^2 +R^2\,\left(\frac{\partial W}{\partial z}\right)^2 -2R^2\,\left[E-\Phi(R,z)\right]=-\left(\frac{\partial W}{\partial \phi}\right)^2\,. \end{equation} If we posit that \(W(R,\phi,z;\vec{C}) = W_{Rz}(R,z;\vec{C})+W_\phi(\phi;\vec{C})\), then the left-hand side of this equation only depends on \((R,z)\) while the right-hand side only depends on \(\phi\) and each must, therefore, equal a constant that we can identify as \(-L_z^2\) (see Section 9.1). We then immediately find that \(W_\phi(\phi;\vec{C}) = \phi\,L_z\). The separability of the \(\phi\) part of Hamilton’s characteristic function \(W(\cdot)\) should come as no surprise given the discussion earlier in this chapter. The equation for \(W_{Rz}(R,z;\vec{C})\) becomes \begin{equation}\label{eq-diskorbit-HJ-Econserved-cylindrical-rzcomponents} \left(\frac{\partial W_{Rz}}{\partial R}\right)^2 +\left(\frac{\partial W_{Rz}}{\partial z}\right)^2 =2\left[E-\Phi(R,z)\right]-{L_z^2 \over R^2}\,. \end{equation} This partial differential equation cannot, in general, be solved using a separation-of-variables approach.
9.4.2. The adiabatic approximation¶
To make progress towards solving Equation (9.36), we can make use of the fact that orbits approximately separate into planar and vertical components as discussed in Section 9.2 above; in the context of action-angle coordinates this is known as the adiabatic approximation (Binney 2010). Positing that the potential separates as in Equation (9.8), Equation (9.9) shows that the Hamiltonian separates into separate \(R\) and \(z\) components. Applying the same potential separation to Equation (9.36) for some to-be-determined \(\tilde{R}\), we find that we can separate the solution as \(W_{Rz}(R,z;\vec{C}) = W_R(R;\vec{C})+W_z(z;\vec{C})\) satisfying \begin{equation}\label{eq-diskorbit-HJ-Econserved-cylindrical-rzcomponents-separate} \left(\frac{\partial W_R}{\partial R}\right)^2 -2\left[E-\Phi_R(R)\right] +{L_z^2 \over R^2} =-\left(\frac{\partial W_z}{\partial z}\right)^2-2\Phi_z(z;\tilde{R})\,. \end{equation} The left-hand side of this equation only depends on \(R\) while the right-hand side only depends on \(z\) and, thus, each needs to equal a constant. This constant is equal to \(-2E_z\), where \(E_z = v_z^2/2 + \Phi_z(z;\tilde{R})\) is the vertical (specific) energy from Section 9.2. The full solution of the Hamilton-Jacobi equation then follows \begin{align}\label{eq-diskorbit-HJ-Econserved-solution} W&(R,\phi,z;E,L_z,E_z) = W_\phi(\phi;L_z)+W_z(z;E_z)+W_R(R;E_R,L_z) \\ & = \phi\,L_z \pm \int_{-z_\mathrm{max}}^z\mathrm{d}z'\,\sqrt{2\left[E_z-\Phi_z(z';\tilde{R})\right]} \pm \int_{R_p}^R\mathrm{d}R'\,\sqrt{2\,\left[E_R-\Phi_R(R')\right]-{L_z^2\over R'^2}}\,,\nonumber \end{align} where \(z_\mathrm{max}\) is the maximum height an orbit with vertical energy \(E_z\) reaches in \(\Phi_z[z';\tilde{R}]\) (the positive value of \(z\) where the integrand of the second integral equals zero), \(R_p\) is the pericenter of the orbit in the planar potential \(\Phi_R(R)\) (the smaller of the two \(R\) where the integrand of the third integral equals zero; the bigger value is the planar apocenter \(R_a\)), and \(E_R = E-E_z\) is the planar contribution to the energy. The signs of the integrals must be chosen such that \(W\) increases monotonically around the orbit. The actions then follow directly from Equation (3.67) and we find that \begin{align}\label{eq-diskorbit-HJ-Econserved-solution-actions-Jphi} J_\phi & = L_z\,,\\ J_R & = {1\over \pi}\,\int_{R_p}^{R_a}\mathrm{d}R\,\sqrt{2\,\left[E_R-\Phi_R(R)\right]-\frac{L_z^2}{R^2}}\,,\label{eq-diskorbit-HJ-Econserved-solution-actions-JR}\\ J_z & = {2\over \pi}\,\int_{0}^{z_\mathrm{max}}\mathrm{d}z\,\sqrt{2\left[E_z-\Phi_z(z;\tilde{R})\right]}\,.\label{eq-diskorbit-HJ-Econserved-solution-actions-Jz} \end{align}
When computing the actions for a given phase-space point at position \((R,\phi,z)\), we can simply set \(\tilde{R} = R\). Alternatively, we can use \(\tilde{R} = R_g\) with \(R_g\) the guiding-center radius from Equation (9.7). In practice, there is little difference in the actions one obtains. Frequencies and angles can be determined using a similar approach as we used for spherical potentials (see Chapter 4.4).
As an example of actions computed using the adiabatic approximation, we compute the radial and vertical actions \(J_R\) and \(J_z\) for the orbit displayed in Figure 9.2, which is quite eccentric compared to most disk orbits, and compare them to the actions of the more typical disk orbit from Figure 9.8 in Figure 9.9.
[24]:
from galpy import potential
from galpy.orbit import Orbit
from astropy import visualization
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.)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o_highjz= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
R, E, Lz= 0.65, -1.385, 0.6
vR= 0.02
o_lowjz= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,101)
o_highjz.integrate(ts,potential.MWPotential2014)
o_lowjz.integrate(ts,potential.MWPotential2014)
with visualization.quantity_support():
figure(figsize=(8,4))
subplot(1,2,1)
line_highjz = o_highjz.plot(
d1='t',
d2=lambda t: o_highjz(t).jr(pot=potential.MWPotential2014,type='adiabatic'),
ylabel=r'$J_R\,(\mathrm{kpc\,km\,s}^{-1})$',
yrange=[0.,65.],
gcf=True)
line_lowjz = o_lowjz.plot(
d1='t',
d2=lambda t: o_lowjz(t).jr(pot=potential.MWPotential2014,type='adiabatic'),
overplot=True)
hijz_jR= o_highjz(ts).jr(pot=potential.MWPotential2014,type='adiabatic')
lojz_jR= o_lowjz(ts).jr(pot=potential.MWPotential2014,type='adiabatic')
galpy_plot.text(2,50.,rf'$\sigma_{{J_R}} = {100.*numpy.std(hijz_jR)/numpy.mean(hijz_jR):.0f}\,\%$',
size=17.,color=line_highjz[0].get_color())
galpy_plot.text(2,10.,rf'$\sigma_{{J_R}} = {100.*numpy.std(lojz_jR)/numpy.mean(lojz_jR):.0f}\,\%$',
size=17.,color=line_lowjz[0].get_color())
# Also overplot the true action computed using the convergent method from the next section
axhline(
o_highjz().jr(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_highjz[0].get_color(),
ls='--',
zorder=-1
)
axhline(
o_lowjz().jr(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_lowjz[0].get_color(),
ls='--',
zorder=-1
)
subplot(1,2,2)
line_highjz = o_highjz.plot(
d1='t',
d2=lambda t: o_highjz(t).jz(pot=potential.MWPotential2014,type='adiabatic'),
ylabel=r'$J_z\,(\mathrm{kpc\,km\,s}^{-1})$',
yrange=[0.,65.],
gcf=True)
line_lowjz = o_lowjz.plot(
d1='t',
d2=lambda t: o_lowjz(t).jz(pot=potential.MWPotential2014,type='adiabatic'),
overplot=True)
hijz_jz= o_highjz(ts).jz(pot=potential.MWPotential2014,type='adiabatic')
lojz_jz= o_lowjz(ts).jz(pot=potential.MWPotential2014,type='adiabatic')
galpy_plot.text(2,45.,rf'$\sigma_{{J_z}} = {100.*numpy.std(hijz_jz)/numpy.mean(hijz_jz):.0f}\,\%$',
size=17.,color=line_highjz[0].get_color())
galpy_plot.text(2,6.,rf'$\sigma_{{J_z}} = {100.*numpy.std(lojz_jz)/numpy.mean(lojz_jz):.0f}\,\%$',
size=17.,color=line_lowjz[0].get_color())
# Also overplot the true action computed using the convergent method from the next section
axhline(
o_highjz().jz(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_highjz[0].get_color(),
ls='--',
zorder=-1
)
axhline(
o_lowjz().jz(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_lowjz[0].get_color(),
ls='--',
zorder=-1
)
tight_layout()
Figure 9.9: Actions computed using the adiabatic approximation.
The actions for the orbit from Figure 9.2 are those with the larger value. We see that the actions fluctuate by \(\approx 10\%\) around this orbit, reflecting the fact that the separability assumption is not great for this orbit. The assumption is much better for the more typical orbit, which has actions that only fluctuate by \(\approx 2\%\) and are, thus, quite close to being conserved as the actions should. The dashed lines show the correct action computed using the method from Section 9.4.4 below and we see that the adiabatic-approximation actions fluctuate around the true actions. Thus, for typical disk orbits, the adiabatic approximation returns relatively precise and accurate actions.
9.4.3. The Staeckel approximation¶
Because disky galactic potentials in general do not separate in \(R\) and \(z\), it may appear that the adiabatic separability approximation is the best we can do, but this is not the case. The most general coordinate system in which the Hamilton-Jacobi equation can be solved using the separation-of-variables approach is the ellipsoidal coordinate system, which we will discuss in more detail in Chapter 13.4.1 in the context of orbits in triaxial potentials. For axisymmetric systems, this coordinate system reduces to the oblate and prolate spheroidal coordinate systems discussed in Chapter 12.2.1 (see Figure 12.8 for an illustration of these coordinate systems). In order for the Hamilton-Jacobi equation to fully separate, the potential must also be separable in these ellipsoidal or spheroidal coordinates. For oblate potentials, which includes essentially all disky potentials unless the dark-matter halo is extremely prolate, the relevant coordinate system is that of prolate spheroidal coordinates (de Zeeuw 1985), which generalize the cylindrical \((R,z)\) coordinates to \((u,v)\) coordinates using the transformation \begin{align}\label{eq-prolate-coords} R & = \Delta\, \sinh u\,\sin v\,;\quad z = \Delta\, \cosh u\,\cos v\,, \end{align} where \(\Delta\) is a constant focal length that is a free parameter of the transformation. The transformation of the momenta in the meridional plane follows from the generating function of the third kind for this transformation: \(F_3(p_R,p_z,u,v) =\) \(p_R\,R(u,v)+p_z\,z(u,v)\). That disk orbits approximately separate in a well-chosen prolate spheroidal coordinate system can be understood qualitatively by overlaying the orbits shown in Figures 9.2 and 9.3 on the grid lines of a prolate spheroidal coordinate system with \(\Delta = 3.2\,\mathrm{kpc}\) (see below for how this value was chosen); this is shown in Figure 9.10.
[25]:
from galpy import potential
from galpy.orbit import Orbit
from galpy.util import coords
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.)
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.,1001)
o.integrate(ts,potential.MWPotential2014)
figure(figsize=(10,6))
o.plot(lw=0.75,xrange=[2.5,8.6],yrange=[-2.4,2.4],zorder=10,gcf=True)
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
o.integrate(ts,potential.MWPotential2014)
o.plot(lw=0.5,overplot=True,zorder=9)
# Now plot the prolate spheroidal coordinate grid
delta= 3.2 # Focal length
# First plot curves of constant u
us= numpy.linspace(0.1,1.475,9)*1.3
vs= numpy.linspace(0.,numpy.pi,101)
for tu in us:
Rs,zs= coords.uv_to_Rz(tu*numpy.ones_like(vs),vs,
delta=delta)
plot(Rs,zs,'k-')
# Then plot curves of constant v
us= numpy.linspace(0.,2.,101)
vs= numpy.linspace(numpy.pi/3.,2.*numpy.pi/3.,8)
for tv in vs:
Rs,zs= coords.uv_to_Rz(us,tv*numpy.ones_like(us),
delta=delta)
plot(Rs,zs,'k--')
Figure 9.10: Disk orbits approximately align with grid lines of a well-chosen prolate spheroidal coordinate system.
In prolate spheroidal coordinates, the effective Hamiltonian becomes \begin{equation}\label{eq-diskorbits-Hamiltonian-effective-prolatespheroidal} H_\mathrm{eff}(u,v,p_u,p_v;L_z) = {p^2_u+p^2_v \over 2\,\Delta^2\,(\sinh^2 u + \sin^2v)}+\frac{L_z^2}{2\Delta^2\,\sinh^2u\,\sin^2v} + \Phi(u,v)\,. \end{equation} We leave the derivation of this result as an exercise. It is then easy to derive that when the potential is of Staeckel form and said to be a Staeckel potential, \begin{equation}\label{eq-diskorbits-staeckelpot} \Phi(u,v) = {U(u)-V(v) \over \sinh^2 u + \sin^2 v}\,, \end{equation} the Hamilton-Jacobi equation can be solved using the separation-of-variables technique and all orbits are decoupled oscillations in the \(u\) and \(v\) coordinates (de Zeeuw 1985). The separation constant is the third integral \begin{equation}\label{eq-diskorbits-staeckelpot-I3} I_3 = E\,\sinh^2 u - U(u) -{p_u^2 \over 2\Delta^2} - {L_z^2 \over 2\Delta^2\,\sinh^2 u} = -\left[E\,\sin^2 v + V(v) -{p_v^2 \over 2\Delta^2} - {L_z^2 \over 2\Delta^2\,\sin^2 v}\right]\,. \end{equation}
All actions, angles, and frequencies can be calculated using one-dimensional numerical integrals similar in form and computational complexity to those involved in the spherical action-angle transformation or the adiabatic approximation discussed above. In particular, in addition to the usual \(J_\phi = L_z\), the actions are \begin{align}\label{eq-diskorbits-staeckelpot-actions-Ju} J_u & = {1 \over \pi}\int_{u_\mathrm{min}}^{u_\mathrm{max}}\mathrm{d} u\,\sqrt{2\Delta^2 \left[E\,\sinh^2 u - I_3-U(u)\right]-{L_z^2 \over \sinh^2 u}}\,,\\ J_v & = {2 \over \pi}\int_{v_\mathrm{min}}^{\pi/2}\ \mathrm{d} v\,\sqrt{2\Delta^2\,\left[\,E\,\sin^2 v\ +\ I_3+V(v)\kern0.1em\right]-{L_z^2 \over \sin^2 v}}\,,\label{eq-diskorbits-staeckelpot-actions-Jv} \end{align} where \(u_\mathrm{min}\), \(u_\mathrm{max}\), and \(v_\mathrm{min}\) are the positive values at which the relevant integrands vanish (they are, respectively, essentially the pericenter, apocenter, and maximum height of the orbit in prolate spheroidal coordinates; see Mackereth & Bovy 2018). Because \(u\) roughly corresponds to radius and \(v\) to height in the prolate spheroidal coordinate system, we usually refer to \(J_u\) as the radial action and we generally write it as \(J_R\), and similarly we call \(J_v\) the vertical action and write it as \(J_z\), in analogy with the actions computed in the adiabatic approximation. We do not list the equations for the frequencies and the angles here, but relevant expressions can be found in, e.g., de Zeeuw (1985) and Binney (2012a).
Unfortunately, aside from the Kuzmin model, none of the popular disk potentials from Chapter 7 are of Staeckel form, so we cannot directly apply the formalism from the last paragraph to them. But if we can approximate the gravitational potential as following Equation (9.44), we can compute action-angle coordinates using the Staeckel formalism. One option is to fit the functions \(U(u)\) and \(V(v)\) over the volume occupied by an orbit to minimize the difference between the true potential and its Staeckel approximation (Sanders 2012). However, this is a laborious procedure, as integrating an orbit to determine its volume and performing a minimization of the difference between the potential and its Staeckel approximation is computationally expensive and difficult to parallelize. But we can also directly perform a split of an oblate potential \(\Phi(R,z)\) into \(U(u)\) and \(V(v)\) in a manner that is similar in spirit to the split from Equation (9.8) in the adiabatic approximation.
Specifically, choosing a reference value \(u_0\), we can define \begin{align}\label{eq-diskorbits-staeckelapprox-Uu} U(u) & = \cosh^2 u \,\Phi\left(u,{\pi\over 2}\right)\,,\\ V(v) & = \cosh^2 u_0 \,\Phi\left(u_0,{\pi\over 2}\right)-\left(\sinh^2 u_0+\sin^2 v\right)\,\Phi\left(u_0,v\right)\,.\label{eq-diskorbits-staeckelapprox-Vv} \end{align} Compared to the split from Equation (9.8), \(U(u)\) represents the planar potential \(\Phi(R,0)\) like \(\Phi_R(R)\) does in the adiabatic approximation, because \(v = \pi/2\) corresponds to \(z=0\). The function \(V(v)\) is similar to \(\Phi_z(z;\tilde{R})\), but instead of representing the vertical potential at \(\tilde{R}\), \(V(v)\) corresponds to the potential along the curve \(u=u_0\). It is easy to show that for a potential that is actually of Staeckel form \(\Phi(u,v) = \left[\tilde{U}(u)-\tilde{V}(v)\right]/\left[\sinh^2 u + \sin^2 v\right]\), Equation (9.48) reduces to \(U(u) = \tilde{U}(u)-V(\pi/2)\) and Equation (9.49) reduces to \(V(v) = \tilde{V}(v)-V(\pi/2)\); the constant \(V(\pi/2)\) simply shifts \(I_3\) by \(V(\pi/2)\) and, thus, has no effect on the actions, angles, and frequencies. Approximating the gravitational potential using Equations (9.48) and (9.49) is known as the Staeckel approximation (or the Staeckel fudge; Binney 2012a).
The final ingredient in approximating an oblate potential as a Staeckel potential is choosing a good value of the focal length parameter \(\Delta\) that defines the prolate coordinate system. This is an important parameter, because how close a potential is to having Staeckel form strongly depends on \(\Delta\) (e.g., how well an orbit aligns with the prolate spheroidal coordinate system as in Figure 9.10 hinges on \(\Delta\)). One way to determine \(\Delta\) is to use the fact that for a Staeckel potential \(\partial^2 /\partial u \partial v \left[\left(\sinh^2 u + \sin^2 v\right)\,\Phi\right] = 0\) (see Equation 9.44). This can be worked out to give a relation between \(\Delta\) and the first and second derivatives of the potential in cylindrical coordinates (Bienayme 2009; Sanders 2012) \begin{equation}\label{eq-diskorbits-staeckelapprox-delta} \Delta^2 = z^2 - R^2 +\left[3\left(z\,{\partial \Phi \over \partial R} - R\,{\partial \Phi \over \partial z}\right)+Rz\,\left({\partial^2 \Phi \over \partial R^2}-{\partial^2\Phi \over \partial z^2}\right)\right]\Big/ {\partial^2 \Phi \over \partial R \partial z}\,. \end{equation} For example, for the eccentric orbit shown in Figure 9.10, the typical value of \(\Delta\) calculated along the orbit in this way is \(\approx 3.2\,\mathrm{kpc}\). Thus, for a given phase-space position \((\vec{x},\vec{v})\) in a disky potential, we can compute actions, angles, and frequencies efficiently by (a) computing \(\Delta\) by evaluating Equation (9.50) at the location \(\vec{x}\), (b) choosing a value of \(u_0\), and (c) performing the integrals from Equations (9.46) and (9.47) (and similar ones for the frequencies and angles) using the Staeckel approximation of the potential from Equations (9.48) and (9.49). Similar to the adiabatic approximation, we can simply set \(u_0\) equal to the value implied by the position \((R,z)\) or we can, e.g., set it based on the guiding-center radius (see also Binney 2012a).
To demonstrate the precision of the actions calculated using the Staeckel approximation, we compute the radial and vertical actions along the orbit for the eccentric orbit from Figure 9.2 and the more typical disk orbit from Figure 9.8 as in Figure 9.9 using the Staeckel approximation. We see in Figure 9.11 that both actions are much closer to being constant along the orbit than using the adiabatic approximation, while still accurately fluctuating around the true action:
[26]:
from galpy import potential
from galpy.orbit import Orbit
from astropy import visualization
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.)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o_highjz= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
R, E, Lz= 0.65, -1.385, 0.6
vR= 0.02
o_lowjz= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,101)
o_highjz.integrate(ts,potential.MWPotential2014)
o_lowjz.integrate(ts,potential.MWPotential2014)
delta_highjz= 3.2/8. # Internal galpy units
delta_lowjz= 4.1/8. # Internal galpy units
with visualization.quantity_support():
figure(figsize=(8,4))
subplot(1,2,1)
line_highjz = o_highjz.plot(
d1='t',
d2=lambda t: o_highjz(t).jr(pot=potential.MWPotential2014,type='staeckel',delta=delta_highjz),
ylabel=r'$J_u\,(\mathrm{kpc\,km\,s}^{-1})$',
yrange=[0.,65.],
gcf=True)
line_lowjz = o_lowjz.plot(
d1='t',
d2=lambda t: o_lowjz(t).jr(pot=potential.MWPotential2014,type='staeckel',delta=delta_lowjz),
overplot=True)
hijz_jR= o_highjz(ts).jr(pot=potential.MWPotential2014,type='staeckel',delta=delta_highjz)
lojz_jR= o_lowjz(ts).jr(pot=potential.MWPotential2014,type='staeckel',delta=delta_lowjz)
galpy_plot.text(2,42.5,rf'$\sigma_{{J_u}} = {100.*numpy.std(hijz_jR)/numpy.mean(hijz_jR):.1f}\,\%$',
size=17.,color=line_highjz[0].get_color())
galpy_plot.text(2,10.,rf'$\sigma_{{J_u}} = {100.*numpy.std(lojz_jR)/numpy.mean(lojz_jR):.1f}\,\%$',
size=17.,color=line_lowjz[0].get_color())
# Also overplot the true action computed using the convergent method from the next section
axhline(
o_highjz().jr(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_highjz[0].get_color(),
ls='--',
zorder=-1
)
axhline(
o_lowjz().jr(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_lowjz[0].get_color(),
ls='--',
zorder=-1
)
subplot(1,2,2)
line_highjz = o_highjz.plot(
d1='t',
d2=lambda t: o_highjz(t).jz(pot=potential.MWPotential2014,type='staeckel',delta=delta_highjz),
ylabel=r'$J_v\,(\mathrm{kpc\,km\,s}^{-1})$',
yrange=[0.,65.],
gcf=True)
line_lowjz = o_lowjz.plot(
d1='t',
d2=lambda t: o_lowjz(t).jz(pot=potential.MWPotential2014,type='staeckel',delta=delta_lowjz),
overplot=True)
hijz_jz= o_highjz(ts).jz(pot=potential.MWPotential2014,type='staeckel',delta=delta_highjz)
lojz_jz= o_lowjz(ts).jz(pot=potential.MWPotential2014,type='staeckel',delta=delta_lowjz)
galpy_plot.text(2,49.5,rf'$\sigma_{{J_v}} = {100.*numpy.std(hijz_jz)/numpy.mean(hijz_jz):.2f}\,\%$',
size=17.,color=line_highjz[0].get_color())
galpy_plot.text(2,6.,rf'$\sigma_{{J_v}} = {100.*numpy.std(lojz_jz)/numpy.mean(lojz_jz):.2f}\,\%$',
size=17.,color=line_lowjz[0].get_color())
# Also overplot the true action computed using the convergent method from the next section
axhline(
o_highjz().jz(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_highjz[0].get_color(),
ls='--',
zorder=-1
)
axhline(
o_lowjz().jz(
pot=potential.MWPotential2014,
type='isochroneapprox',
b=0.3,
tintJ=1000.,
ntintJ=100_000,
integrate_method="dop853_c"),
color=line_lowjz[0].get_color(),
ls='--',
zorder=-1
)
tight_layout()
Figure 9.11: Actions computed using the Staeckel approximation.
For typical disk orbits, the Staeckel approximation is usually precise to \(\lesssim 1\%\). Because most disk populations have velocity dispersions that are \(\gtrsim 5\%\) of the circular velocity, this precision is good enough to characterize disk populations. Because the prolate spheroidal coordinate system becomes the spherical one as \(\Delta \rightarrow 0\) (see Equations 9.42, where \(\Delta \cosh u \rightarrow \Delta \sinh u \rightarrow r\)), the Staeckel approximation also works well for stars on halo-like orbits as long we set \(\Delta\) appropriately. The reason that the Staeckel approximation is so much better than the adiabatic approximation is that orbits separate much better into decoupled oscillations in the prolate spheroidal coordinate system. For example, comparing the track of the eccentric orbit from Figure 9.2 in the \((z,v_z)\) plane of the adiabatic approximation to the track in the \((v,p_v)\) plane from the Staeckel approximation in Figure 9.12, we see the track is far more one-dimensional in the latter case.
[27]:
from galpy.util import coords
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o_highjz= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,200.,10001)
o_highjz.integrate(ts,potential.MWPotential2014)
figure(figsize=(8,4))
subplot(1,2,1)
o_highjz.plot(d1='z',d2='vz',lw=0.5,zorder=0,gcf=True,
xrange=[-1.74,1.74],yrange=[-145,145.])
subplot(1,2,2)
delta= 3.2 # Focal length
u,v= coords.Rz_to_uv(o_highjz.R(ts,quantity=False),o_highjz.z(ts,quantity=False),delta=delta)
pu,pv= coords.vRvz_to_pupv(o_highjz.vR(ts,quantity=False),o_highjz.vz(ts,quantity=False),u,v,delta=delta,uv=True)
plot(v/numpy.pi,pv,lw=0.5)
xlabel(r'$v/\pi$')
ylabel(r'$p_v\,(\mathrm{km\,s}^{-1})$')
tight_layout();
Figure 9.12: Quality of the decoupled “vertical” orbit approximation in the adiabatic approximation (left panel) and the Staeckel approximation (right panel).
The same happens for the track in the \((u,p_u)\) plane versus that in the \((R,v_R)\) plane.
9.4.4. Convergent methods¶
While the adiabatic and Staeckel approximations are highly useful for quickly determining action-angle coordinates for large numbers of bodies, the main issue with them is that they are limited in precision in a way that cannot be improved within the context of their method. This is the reason that the Staeckel approximation is often referred to as a fudge instead. The most general axisymmetric coordinate system in which the Hamilton-Jacobi equation for a disky potential can be solved is the prolate spheroidal coordinate system from the Staeckel approximation, and there is, therefore, no hope of exactly solving the Hamilton-Jacobi equation using separation-of-variables for disky potentials. We thus need a different approach to compute actions, angles, and frequencies to arbitrary precision for disky potentials.
One alternative way of obtaining action-angle coordinates is to directly constrain a canonical transformation \((\vec{x},\vec{v}) \rightarrow (\vec{J},\boldsymbol{\theta})\) by demanding that the new coordinates are actions and angles of the potential we are working in (often referred to as the target potential in this context). This may seem impossible, but it can be done. Current implementations of this technique all follow the principle of the torus-mapping approach of McGill & Binney (1990) and Binney & Kumar (1993), an approach originally designed to deliver the inverse transformation \((\vec{J},\boldsymbol{\theta}) \rightarrow (\vec{x},\vec{v})\). In the torus-mapping approach, we split the full canonical transformation between \((\vec{x},\vec{v})\) and \((\vec{J},\boldsymbol{\theta})\) into (i) a transformation of \((\vec{x},\vec{v})\) to the actions and angles \((\vec{J}^A,\boldsymbol{\theta}^A)\) of a fixed auxiliary potential and (ii) a canonical transformation \(F_2(\boldsymbol{\theta}^A,\vec{J})\) to the correct action-angle coordinates for the target potential. The first part of this transformation is canonical, because the action-angle coordinates of any potential are a canonical coordinate system (that is, canonical coordinate transformations are not directly tied to the potential). Usually an isochrone potential is used as the auxiliary potential, because its actions and angles can be obtained analytically and the action-angle transformation can be inverted. For the second part of the transformation, the periodicity of the angle coordinates suggests the general form \begin{equation}\label{eq-diskorbits-actions-convergent-generatingfunction} F_2(\boldsymbol{\theta}^A,\vec{J}) = \boldsymbol{\theta}^A \cdot \vec{J} +2\sum_{\vec{n}>0}{S_n(\vec{J})\,\sin\left(\vec{n}\cdot \boldsymbol{\theta}^A\right)}\,, \end{equation} where the \(\vec{n} > 0\) sum runs over the integer three-dimensional half-space that does not include the origin; this form includes all possible mappings for a static potential (McGill & Binney 1990). This generating function gives the following canonical transformation \begin{equation}\label{eq-diskorbits-actions-convergent-generatingfunction-actiontransform} \vec{J}^A = {\partial F_2(\boldsymbol{\theta}^A,\vec{J}) \over \partial \boldsymbol{\theta}^A} = \vec{J} + 2\sum_{\vec{n}>0}{\vec{n}\,S_n(\vec{J})\,\cos\left(\vec{n}\cdot \boldsymbol{\theta}^A\right)}\,, \end{equation} and \begin{equation}\label{eq-diskorbits-actions-convergent-generatingfunction-angletransform} \boldsymbol{\theta} = {\partial F_2(\boldsymbol{\theta}^A,\vec{J}) \over \partial \vec{J}} = \boldsymbol{\theta}^A +2\sum_{\vec{n} > 0} {{\partial S_{\vec{n}}(\vec{J})\over \partial \vec{J}}\,\sin\left(\vec{n}\cdot\boldsymbol{\theta}^A\right)}= \boldsymbol{\theta}_0 + \boldsymbol{\Omega}(\vec{J})\,t \,, \end{equation} where we have specified that the angles \(\boldsymbol{\theta}\) in the target potential should linearly increase in time according to the frequency \(\boldsymbol{\Omega}(\vec{J})\).
For any set of \(S_\vec{n}(\vec{J})\) functions, Equation (9.51) defines a canonical mapping of the tori of the auxiliary potential onto a different set of tori, which are the tori of some potential. It can then be demonstrated that if the target Hamiltonian is constant along the tori defined by Equation (9.52), then the actions \(\vec{J}\) are the correct actions for the target potential (McGill & Binney 1990). Thus, in torus mapping, we can determine the \(S_\vec{n}(\vec{J})\) functions by demanding that \(H(\vec{J}) = \mathrm{constant}\) along a torus defined by \(\vec{J}\) and evaluated on a grid of \(\boldsymbol{\theta}^A\). Angles and frequencies can be obtained by numerically integrating small sections of the orbit in the target potential at a grid of \(\boldsymbol{\theta}^A\) and performing a linear fit for \([\boldsymbol{\theta}_0,\boldsymbol{\Omega}(\vec{J}),\partial S_{\vec{n}} / \partial \vec{J}]\) using Equation (9.53) (Kaasalainen & Binney 1994). In practice we need to, of course, limit the number of terms in the Fourier expansion of the generating function to a large enough number to approximate the torus well and a small enough number to efficiently perform the optimization of the parameters. Aside from potential numerical issues, tori can be mapped arbitrarily precisely in this way by including enough expansion terms, and we thus obtain the mapping \((\vec{J},\boldsymbol{\theta}) \rightarrow (\vec{x},\vec{v})\). However, this mapping is not easily inverted (McMillan & Binney 2008).
To instead map \((\vec{x},\vec{v}) \rightarrow (\vec{J},\boldsymbol{\theta})\), we can simply determine the expansion coefficients \([S_{\vec{n}}(\vec{J}),\partial S_{\vec{n}} / \partial \vec{J}]\) and the actions, angles, and frequencies \((\vec{J},\boldsymbol{\theta},\boldsymbol{\Omega})\) using direct fits of Equations (9.52) and (9.53) to an integrated orbit \((\vec{x},\vec{v})[t]\) (Bovy 2014; Sanders & Binney 2014). That is, starting from \((\vec{x},\vec{v})\), we numerically integrate the orbit—in practice, both forwards and backwards for increased numerical stability—convert each point \((\vec{x}_i,\vec{v}_i)]\) along the orbit to \((\vec{J}^A_i,\boldsymbol{\theta}^A_i)\), and we create a system of equations by plugging each of these into Equations (9.52) and (9.53). These equations are linear in the unknown parameters and can, thus, be efficiently solved using simple linear algebra. Moreover, Bovy (2014) demonstrated that for non-resonant orbits that fill the entire \(0\) to \(2\pi\) range in each auxiliary angle, each action can be determined as \begin{equation} J_i = {\int_0^{2\pi}\mathrm{d}\theta^A_i\,J^A_i \over \int_0^{2\pi}\mathrm{d}\theta^A_i}\,, \end{equation} that is, the true action is the average of the auxiliary actions over the auxiliary angle. By including more and more Fourier expansion terms for the generating function and integrating the orbit longer and more finely, the actions and angles can be determined arbitrarily well. Note that because \(\boldsymbol{\Omega}\) is an unknown that needs to be determined in the fit, we also obtain the frequencies. Similar to the Staeckel approximation, we need to first determine a good auxiliary isochrone potential to use, which requires some trial and error. A good starting point for the search is to match the circular velocity and the slope of the rotation curve of the target potential at the point where one is computing the actions.
The method described in this section applies to any static potential, including non-axisymmetric, triaxial potentials such as those discussed in Chapter 12. For axisymmetric potentials \(J^A_\phi = J_\phi = L_z\) and, therefore, all of the \(S_\vec{n}\) with \(n_\phi \neq 0\) are equal to zero, thus reducing the complexity of the fit. As an example of the application of this method to an axisymmetric potential, the true actions computed using this method are shown as dashed lines in Figures 9.9 and 9.11.
Because the convergent method from this section requires a laborious orbit integration, it is much slower to compute than the adiabatic and staeckel approximations and is, in particular, not practical to apply to large data sets. However, it is useful in contexts where very precise actions, frequencies, and angles are required, for example, when building models of cold tidal streams (see Chapter 19.4.2.5; Bovy 2014).
9.4.5. The action diamond¶
As discussed in Chapter 3.4.3, the fact that actions are adiabatically invariant makes them the best quantities to use for labeling orbits in galactic potentials. All orbits in axisymmetric potentials and, as we will see in Chapter 13 and beyond, almost all orbits in static and rotating triaxial potentials, are orbits that loop around the center with a well-defined sense of rotation, while oscillating in the two directions perpendicular to the loop. The actions directly quantify these three components of the motion: \(J_\phi\) (\(L_z\) in an axisymmetric potential) represents the location of the main loop, while the radial and vertical actions (\(J_R\) and \(J_z\) in the adiabatic approximation, or \(J_u\) and \(J_v\) in the Staeckel approximation; we generally refer to this pair as \(J_R\) and \(J_z\)) characterize the amplitude of the oscillations perpendicular to the loop.
Often we are interested only in an orbit’s shape rather than its absolute extent within the galaxy. For example, when observing stars near the Sun, the overall locations of their orbits is largely set by the volume of the observational survey rather than the intrinsic properties of the Milky Way, so the distribution of absolute extents is not reflective of the Milky Way itself, but rather of how we observe it. It is, therefore, preferable to work with dimensionless representations of the orbital properties, which are less affected by observational biases. Eccentricity is one such dimensionless property, but on its own it reflects only the amplitude of the radial oscillations. To characterize both the radial and vertical extent of the orbit, Vasiliev (2019) introduced a scaled version of action space \((J_R,J_\phi,J_z)\) consisting of \((J_z-J_R)\) and \(J_\phi\) each normalized by the angular momentum \(L_c(E)\) of a circular orbit with the same energy. In subsequent work, different normalizations have been used more often, with the most popular one being \(J_\mathrm{tot} = |J_\phi|+J_R+J_z\) (note that \(J_R\) and \(J_z\) are always positive), but some authors also use \(J_\mathrm{tot} = \sqrt{J_\phi^2+J_R^2+J_z^2}\). When using \(J_\mathrm{tot} = |J_\phi|+J_R+J_z\), orbits are restricted to lying in a square rotated by \(45^\circ\) and for this reason, this scaled action space is known as the action diamond (Lane et al. 2022). As an example, we can plot the action diamond for the Milky Way’s globular clusters, satellite galaxies, and some disk stars in Figure 9.13.
[28]:
from galpy.potential import MWPotential2014
from galpy.orbit import Orbit
def action_diamond_outline(color='0.7',ls='-',labels=True):
plot([0,1],[-1,0],color=color,ls=ls)
plot([0,-1],[-1,0],color=color,ls=ls)
plot([0,1],[1,0],color=color,ls=ls)
plot([0,-1],[1,0],color=color,ls=ls)
if labels:
galpy_plot.text(1.075,0,r'$\mathrm{prograde}$',fontsize=17.,va='center',ha='center',rotation=270.)
galpy_plot.text(-1.075,0,r'$\mathrm{retrograde}$',fontsize=17.,va='center',ha='center',rotation=90.)
galpy_plot.text(0,1.075,r'$\mathrm{polar}$',fontsize=17.,va='center',ha='center',rotation=0.)
galpy_plot.text(0,-1.075,r'$\mathrm{radial}$',fontsize=17.,va='center',ha='center',rotation=0.)
galpy_plot.text(0.575,-.55,r'$\mathrm{in\!-\!plane}$',fontsize=17.,va='center',ha='center',rotation=45.)
galpy_plot.text(-0.575,-.55,r'$\mathrm{in\!-\!plane}$',fontsize=17.,va='center',ha='center',rotation=-45.)
xlim(-1.175,1.175)
ylim(-1.175,1.175)
# Globular clusters
o= Orbit.from_name("MW globular clusters")
JR,Jphi,Jz = o.jr(pot=MWPotential2014,type='staeckel'), o.Lz(), o.jz(pot=MWPotential2014,type='staeckel')
Jtot= numpy.fabs(Jphi)+JR+Jz
figure(figsize=(6.5,6.5))
line_gcs= plot(Jphi/Jtot,(Jz-JR)/Jtot,'o')
annotate_indx= numpy.argmin((Jphi/Jtot+1.)**2.+((Jz-JR)/Jtot+0.)**2.)
annotate(r"$\mathrm{Globular\ clusters}$",((Jphi/Jtot)[annotate_indx],((Jz-JR)/Jtot)[annotate_indx]),
(-0.3625,-1.05),ha='right',
arrowprops=dict(color=line_gcs[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color=line_gcs[0].get_color(),zorder=1)
annotate_indx= numpy.argmin((Jphi/Jtot+0.1)**2.+((Jz-JR)/Jtot+0.8)**2.)
annotate(r"$\mathrm{Globular\ clusters}$",((Jphi/Jtot)[annotate_indx],((Jz-JR)/Jtot)[annotate_indx]),
(-0.63,-1.05),ha='right',
arrowprops=dict(color=line_gcs[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color='w',zorder=-1)
# Satellite galaxies
o= Orbit.from_name("MW satellite galaxies")
JR,Jphi,Jz = o.jr(pot=MWPotential2014,type='staeckel'), o.Lz(), o.jz(pot=MWPotential2014,type='staeckel')
Jtot= numpy.fabs(Jphi)+JR+Jz
line_gcs= plot(Jphi/Jtot,(Jz-JR)/Jtot,'o')
annotate_indx= numpy.argmin((Jphi/Jtot+0.1)**2.+((Jz-JR)/Jtot-0.8)**2.)
annotate(r"$\mathrm{Satellite\ galaxies}$",((Jphi/Jtot)[annotate_indx],((Jz-JR)/Jtot)[annotate_indx]),
(-0.75,1.),ha='center',
arrowprops=dict(color=line_gcs[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color=line_gcs[0].get_color(),zorder=1)
annotate_indx= numpy.argmin((Jphi/Jtot+0.9)**2.+((Jz-JR)/Jtot-0.2)**2.)
annotate(r"$\mathrm{Satellite\ galaxies}$",((Jphi/Jtot)[annotate_indx],((Jz-JR)/Jtot)[annotate_indx]),
(-0.63,1.0),ha='right',
arrowprops=dict(color=line_gcs[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color='w',zorder=-1)
# Also the Sun
o= Orbit()
JR,Jphi,Jz = o.jr(pot=MWPotential2014,type='staeckel'), o.Lz(), o.jz(pot=MWPotential2014,type='staeckel')
Jtot= numpy.fabs(Jphi)+JR+Jz
line_sun= plot(Jphi/Jtot,(Jz-JR)/Jtot,'s',ms=7.)
annotate(r"$\mathrm{Sun}$",(Jphi/Jtot,(Jz-JR)/Jtot),
(0.95,0.65),ha='left',
arrowprops=dict(color=line_sun[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color=line_sun[0].get_color())
# Also Arcturus, a metal-poor standard star
o= Orbit.from_name("Arcturus")
JR,Jphi,Jz = o.jr(pot=MWPotential2014,type='staeckel'), o.Lz(), o.jz(pot=MWPotential2014,type='staeckel')
Jtot= numpy.fabs(Jphi)+JR+Jz
line_arcturus= plot(Jphi/Jtot,(Jz-JR)/Jtot,'s',ms=7.)
annotate(r"$\mathrm{Arcturus}$",(Jphi/Jtot,(Jz-JR)/Jtot),(0.45,-1.),ha='center',
arrowprops=dict(color=line_arcturus[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color=line_arcturus[0].get_color())
# Also Arcturus, a metal-poor standard star
o= Orbit.from_name("LHS 1815")
JR,Jphi,Jz = o.jr(pot=MWPotential2014,type='staeckel'), o.Lz(), o.jz(pot=MWPotential2014,type='staeckel')
Jtot= numpy.fabs(Jphi)+JR+Jz
line_arcturus= plot(Jphi/Jtot,(Jz-JR)/Jtot,'s',ms=7.)
annotate(r"$\mathrm{LHS\ 1815}$",(Jphi/Jtot,(Jz-JR)/Jtot),(0.925,-0.85),ha='center',
arrowprops=dict(color=line_arcturus[0].get_color(),width=.5,headwidth=0.,shrink=0.005),
fontsize=17.,color=line_arcturus[0].get_color())
# Plot the diamond itself
action_diamond_outline()
xlabel(r'$J_\phi / J_\mathrm{tot}$')
ylabel(r'$\left(J_z-J_R\right) / J_\mathrm{tot}$')
gca().set_aspect('equal');
#gca().axis('off'); # uncomment this line to show a pure diamond
Figure 9.13: Action diamond for the Milky Way’s globular clusters, satellite galaxies, the Sun, and a few metal-poor stars.
We have labeled the characteristics of the orbits at the extremes of the diamond. Objects with positive angular momentum have prograde orbits and lie on the right side of the diamond, while retrograde stars lie on the left side. Stars with \(J_z = 0\) are confined to the potential’s mid-plane and they lie along the bottom edges of the diamond. Because disk stars have orbits that remain relatively close to the mid-plane and that are prograde, disk stars occupy the right corner of the diagram, hugging the in-plane lower edge. We have include the location of the Sun, the metal-poor giant Arcturus, and an exoplanet-hosting M dwarf with a particularly large vertical excursion (LHS 1815; Gan et al. 2020) to illustrate where disk stars lie in the action diamond. Outside the small corner occupied by disk stars, halo orbits populate the diagram, with highly radial orbits near the bottom corner and more circular halo orbits that loop over the Galactic poles near the top of the diagram. We see that globular clusters populate the diagram almost evenly, while satellite galaxies avoid planar orbits. This makes sense, as satellite galaxies on planar orbits are quickly destroyed by tidal forces, while the denser globular clusters are more resilient to tidal disruption.
.
.
.