4.4. Actions and angles in spherical potentials¶
We saw in Chapter 3.4.3 that many orbits in galactic potentials can be described using a set of coordinates consisting of actions \(\vec{J}\) and angles \(\boldsymbol{\theta}\) in which the dynamics is given by the simple set of equations (3.69)-(3.70). These coordinates can be obtained by solving the Hamilton-Jacobi equation in its form of Equation (3.63) for static potentials such as we will assume here. As already discussed in the previous section, this equation can be solved for any static spherical potential and we present the general solution here.
For a static spherical potential \(\Phi(r)\), Equation (3.63) can be written as \begin{equation}\label{eq-spherorbit-HJ-Econserved-spherical} \left(\frac{\partial W}{\partial r}\right)^2 + \frac{1}{r^2\,\sin^2 \theta}\,\left(\frac{\partial W}{\partial \phi}\right)^2+ \frac{1}{r^2}\,\left(\frac{\partial W}{\partial \theta}\right)^2+ 2m^2\,\Phi(r) = 2m^2\,E\,, \end{equation} where \(E\) here denotes the specific energy. Note that while it is somewhat confusing, we continue to use \(\theta\) for the polar coordinate in the spherical coordinate system \((r,\phi,\theta)\) and \(\boldsymbol{\theta}\) for the angles in the action-angle coordinate system. We will only ever use non-bolded angle coordinates when subscripting them with their associated spatial coordinates (e.g., \(\theta_\phi\); see below), so there should be no ambiguity in what \(\theta\) means in any context. Because we are dealing with a fixed potential, a body’s trajectory and thus its action-angle coordinates, do not (non-trivially) depend on its mass and we therefore set \(m=1\) in this section.
Because the gravitational potential is only a function of one of the variables, we can solve Equation (4.45) using the separation-of-variables technique. That is, we posit that the solution has the form \(W(r,\theta,\phi;\vec{C}) = W_r(r;\vec{C})+W_\phi(\phi;\vec{C})+W_\theta(\theta;\vec{C})\) and we can then write Equation (4.45) as \begin{equation} r^2\,\left(\frac{\partial W_r}{\partial r}\right)^2 + 2r^2\,\left[\Phi(r)-E\right] = - \frac{1}{\sin^2 \theta}\,\left(\frac{\partial W_\phi}{\partial \phi}\right)^2- \left(\frac{\partial W_\theta}{\partial \theta}\right)^2\,. \end{equation} Because the left-hand side of this equation only depends on \(r\) and the right-hand side only depends on \((\phi,\theta)\), each side has to be equal to a constant for this equation to hold; we write this constant as \(-L^2\) (see below). Working with the right-hand side of the equation, we can then similarly derive \begin{equation} \left(\frac{\partial W_\phi}{\partial \phi}\right)^2 = \sin^2\theta\,L^2- \sin^2\theta\left(\frac{\partial W_\theta}{\partial \theta}\right)^2\,. \end{equation} This equation’s left-hand side now only depends on \(\phi\) and the right-hand side only depends on \(\theta\), so each side must again equal a constant and we write this constant as \(L_z^2\) (again, see below). Thus, we end up with a system of three ordinary differential equations for the functions \(\left[W_r(r;\vec{C}),W_\phi(\phi;\vec{C}),W_\theta(\theta;\vec{C})\right]\): \begin{align} \left(\frac{\partial W_r}{\partial r}\right)^2 & = 2E-2\Phi(r) -{L^2\over r^2}\,,\\ \left(\frac{\partial W_\phi}{\partial \phi}\right)^2 & = L_z^2\,,\\ \left(\frac{\partial W_\theta}{\partial \theta}\right)^2 & = L^2 - { L_z^2 \over \sin^2\theta}\,. \end{align} Because the partial derivatives of \(W\) are the canonical momenta, using the Lagrangian formalism it is straightforward to demonstrate that \(L\) is the magnitude of the specific angular momentum and \(L_z\) is its \(z\) component, which explains why we wrote the constants as we did (note that we have the freedom to choose the sign of the \(L_z\) separation constant such that this is the case). These differential equations are easy to solve through direct integration and we then obtain the solution of the Hamilton-Jacobi equation \begin{align}\label{eq-spherorbit-HJ-Econserved-solution} W&(r,\theta,\phi;E,L,L_z) = W_\phi(\phi;L_z)+W_\theta(\theta;L,L_z)+W_r(r;E,L) \\ & = \phi\,L_z \pm \int_{\pi/2}^\theta\mathrm{d}\theta'\,\sqrt{L^2-{L_z^2 \over \sin^2\theta'}} \pm \int_{r_p}^r\mathrm{d}r'\,\sqrt{2\,\left[E-\Phi(r')\right]-\frac{L^2}{r'^2}}\,,\nonumber \end{align} where we must choose the signs of the second and third terms on the right such that the integrals increase monotonically along the orbit and \(r_p\) is the pericenter.
We can then evaluate the actions using Equation (3.67). The azimuthal action associated with \(\phi\) is \(J_\phi = (1/2\pi)\,\int_0^{2\pi}\mathrm{d}\phi\,L_z = L_z\). We leave it as an exercise to demonstrate that the latitudinal action \(J_\theta = L-|L_z|\). The radial action also immediately follows: \begin{equation}\label{eq-spherorbit-Jr} J_r = {1\over \pi}\,\int_{r_p}^{r_a}\mathrm{d}r\,\sqrt{2\,\left[E-\Phi(r)\right]-\frac{L^2}{r^2}}\,. \end{equation}
To obtain the angle coordinates \(\boldsymbol{\theta}\), we need the partial derivatives of \(W(r,\theta,\phi;E,L,L_z)\) with respect to the actions (Equation 3.68). Because the expression for \(W(r,\theta,\phi;E,L,L_z)\) from Equation (4.51) is not an explicit function of \(J_\theta\), this is more difficult than it has to be. But we can get around this by using a canonical transformation to a new set of action-angle coordinates \((\vec{J}^L,\boldsymbol{\theta}^L)\) for which the actions are \(\vec{J}^L = (J_r,L,L_z)\) such that \(W\) is an explicit function of the actions \(L\) and \(L_z\). This can be accomplished using the following generating function of the second kind: \begin{equation}\label{eq-spherorbit-HJ-Econserved-generatingfunc-to-L} F_2(\boldsymbol{\theta},\vec{J}^L) = \theta_r\,J^L_r+\theta_\phi\,J^L_\phi +\theta_\theta\,(J^L_\theta-|J^L_\phi|)\,. \end{equation} This generating function preserves the \(r\) and \(\phi\) actions: \(J_r = \partial F_2 / \partial \theta_r = J^L_r\) and similarly \(J_\phi = J^L_\phi = L_z\). The original \(\theta\) action is \(J_\theta = \partial F_2 / \partial \theta_\phi = J^L_\theta-|J^L_\phi|= J^L_\theta-|L_z|\) and we then have \(J^L_\theta = L\). The new angles are \(\theta^L_r = \partial F_2 / \partial J^L_r = \theta_r\), \(\theta^L_\theta = \theta_\theta\), and \begin{equation}\label{eq-spherorbit-HJ-Econserved-solution-thetaphi-as-thetaphiL} \theta^L_\phi = \theta_\phi-\mathrm{sign}(L_z)\,\theta_\theta \end{equation} and transforming between the two sets of angle coordinates is, thus, straightforward. Equation (4.52) demonstrates that \(E\) is only a function of \((J_r,L)\) and we can then write Equation (4.51) in a way that explicitly shows all dependence on the actions \((J_r,L,L_z)\) \begin{align}\label{eq-spherorbit-HJ-Econserved-solution-explicitactions} W&(r,\theta,\phi;J_r,L,L_z) = \\ & \phi\,L_z \pm \int_{\pi/2}^\theta\mathrm{d}\theta'\,\sqrt{L^2-{L_z^2 \over \sin^2\theta'}} \pm \int_{r_p}^r\mathrm{d}r'\,\sqrt{2\,\left[E(J_r,L)-\Phi(r')\right]-{L^2\over r'^2}}\,.\nonumber \end{align} The \(\theta^L_\phi\) angle is then \begin{align}\label{eq-spherorbit-HJ-Econserved-solution-thetaLphi-1} \theta^L_\phi = {\partial W \over \partial L_z} & = \phi \mp \int_{\pi/2}^\theta\mathrm{d}\theta'\,{L_z\over \sin^2\theta'\,\sqrt{L^2-L_z^2/ \sin^2\theta'}} \\& = \phi \mp \mathrm{sign}(L_z)\int_{\pi/2}^\theta\mathrm{d}\theta'\,{1\over \sin\theta'\,\sqrt{\sin^2\theta'\,\sec^2 i-1}} \,, \end{align} using Equation (4.16), we can write this in terms of the inclination angle. Figure 4.1 demonstrates that for an orbit with \(L_z > 0\), \(\dot{\theta} < 0\) as the orbit passes upwards through the line of nodes and, thus, we need to choose the first sign in Equation (4.55) to be negative and, consequently, the sign in Equation (4.56) to be positive. It turns out that the integral can be solved analytically and we find \begin{equation}\label{eq-spherorbit-HJ-Econserved-solution-thetaLphi} \theta^L_\phi = \phi -\mathrm{arcsin}\left(\cot i\,\cot \theta\right) = \Omega\,, \end{equation} where the final result follows, because we can use Figure 4.1 to determine that \(\mathrm{arcsin}\left(\cot i\,\cot \theta\right) = \phi-\Omega\) (because the length of the dash-dotted line is \(\cot i\,\cos \theta\) and this line is perpendicular to the line of nodes). Thus, the \(\theta^L_\phi\) angle coordinate is actually a constant equal to the longitude of the ascending node! This result makes sense, because a static spherical potential has four independent integrals of the motion (the three components of the angular momentum and the energy) and, thus, cannot have three independent angle coordinates. Physically, the angle associated with the oscillation of the orbital plane around an average orbital plane, as discussed near the end of Chapter 3.4.3, is constant because the orbital plane is fixed in a spherical potential.
Next, we obtain the radial angle as \begin{equation}\label{eq-spherorbit-HJ-Econserved-solution-thetaLr} \theta_r = \theta^L_r =\mathrm{sign}(v_r)\,{2\pi \over T_r}\,\int_{r_p}^r\mathrm{d}r'\,{1 \over \sqrt{2\,\left[E(J_r,L)-\Phi(r')\right]-L^2/r'^2}}\,, \end{equation} where we used that \(\partial E / \partial J_r = \Omega_r\) (Equation 3.72) and, as discussed further below, \(\Omega_r = 2\pi/T_r\) with \(T_r\) the radial period from Equation (4.13). From Equation (4.13), it is clear that this angle increases by \(2\pi\) as a body moves from pericenter to apocenter and back. The sign is chosen such that the angle increases monotonically along the orbit. Thus, the radial angle is essentially given by the indefinite version of the integral involved in calculating the radial period.
Finally, the \(\theta_\theta\) angle is \begin{equation}\label{eq-spherorbit-HJ-Econserved-solution-thetaLtheta-1} \theta_\theta = \theta^L_\theta = {\partial W \over \partial L} = \pm \int_{\pi/2}^\theta\mathrm{d}\theta'\,{L \over \sqrt{L^2-L_z^2/ \sin^2\theta'}} \pm \int_{r_p}^r\mathrm{d}r'\,{\partial E / \partial L - L /r'^2\over\sqrt{2\,\left[E(J_r,L)-\Phi(r')\right]-L^2/r'^2}}\,. \end{equation} Again employing Equation (4.16), we can show that the first integral equals \(\mathrm{arcsin}(\cos \theta/ \sin i)\) and Figure 4.1 shows that this angle is in fact the angle \(\psi\) in the orbital plane. Using the fact that \(\partial E / \partial L = \Omega_\theta\) (which we determine below) and using Equation (4.59), we can derive \begin{equation}\label{eq-spherorbit-HJ-Econserved-solution-thetaLtheta} \theta_\theta = \theta^L_\theta = \psi+ {\Omega_\theta \over \Omega_r}\,\theta_r- \mathrm{sign}(v_r)\,\int_{r_p}^r\mathrm{d}r'\,{L\over r'^2\,\sqrt{2\,\left[E(J_r,L)-\Phi(r')\right]-L^2/r'^2}}\,, \end{equation} where we have again chosen the sign of the final term to make sure this angle changes monotonically around the orbit. Thus, the \(\theta\) angle is a combination of the angle in the orbital plane, the radial angle weighted by the frequency ratio, and the indefinite version of the integral appearing in the calculation of \(\Delta \psi\) from Equation (4.14). Using these angles, the original \(\theta_\phi\) angle can be obtained using Equation (4.54).
The frequencies can be obtained using Equation (3.72). Because the Hamiltonian is only a function of \(J_r\) and \(L\) (see Equation 4.52, where \(H[J_r,L]=E\)), we have that \(\Omega_\phi^L = 0\) in accordance with Equation (4.58). Taking the partial derivative of each side of Equation (4.52) with respect to \(J_r\), it is easy to show that \(\Omega^L_r = 2\pi/T_r\), with \(T_r\) the radial period from Equation (4.13). Similarly, taking the partial derivative with respect to \(L\), we can show that \(\Omega^L_\theta = 2\pi/T_\psi\), with \(T_\psi\) the azimuthal period from Equation (4.15). Using that \(\theta_r^L = \theta_r\), \(\theta_\theta^L = \theta_\theta\), and the relation between \(\theta_\phi\) and \(\theta_\phi^L=\Omega\) from Equation (4.54), we have that \(\Omega_r = \Omega^L_r\), \(\Omega_\theta = \Omega_\theta^L\), and \(\Omega_\phi = \mathrm{sign}(L_z)\Omega_\theta\).
Thus, for any spherical potential, we can compute all actions, angles, and frequencies for any orbit using a few one-dimensional integrals. The isochrone potential from Chapter 2.4.4 and of Section 4.2.3 is the most general galactic potential for which all of these integrals can be performed analytically (see Binney & Tremaine 2008 for full expressions). For the isochrone potential, we can in fact invert the action-angle transformation and obtain position–velocity \((\vec{x},\vec{v})\) given \((\vec{J},\boldsymbol{\theta})\) requiring only a single transcendental equation to be solved numerically—a generalization of Kepler’s equation (4.39) (McGill & Binney 1990). Because the time evolution of an orbit in action-angle coordinate is such that actions are conserved and angles increase linearly in time, this means that orbits in the isochrone potential can be obtained without solving a differential equation and we, therefore, refer to them as analytic. We cannot easily invert the action-angle transformation for general spherical potentials.
To illustrate action-angle coordinates in spherical potentials, we compute the angles for the orbit in the isochrone potential shown in Figure 4.8. The following animation displays the orbit in \((x,y)\) (left panel) and \((\theta_r,\theta_\theta)\) (left panel). It is clear that the complicated time behavior from the left panel translates to the simple linear behavior in angle space in the right panel:
[24]:
from galpy import potential
from galpy.orbit import Orbit
from galpy.actionAngle import actionAngleIsochrone
ip= potential.IsochronePotential(normalize=1.,b=1.)
aAI= actionAngleIsochrone(ip=ip)
oi= Orbit([1.,0.25,1.4,0.])
ts= numpy.linspace(0.,40.,3001)
oi.integrate(ts,ip)
oi.animate(
d1=['x',lambda t: aAI.actionsFreqsAngles(oi(t))[6]],
d2=['y',lambda t: aAI.actionsFreqsAngles(oi(t))[8]],
width=700.,
height=350.,
xlabel=['x','θ<sub>r</sub>'],
ylabel=['y','θ<sub>θ</sub>']); # remove ; to show the animation (note that this animation will show jumps when angles go from 2pi --> 0 that are removed for the website animation
A static version of this animation is:
[25]:
from galpy import potential
from galpy.orbit import Orbit
from galpy.actionAngle import actionAngleIsochrone
ip= potential.IsochronePotential(normalize=1.,b=1.)
aAI= actionAngleIsochrone(ip=ip)
oi= Orbit([1.,0.25,1.4,0.])
ts= numpy.linspace(0.,40.,301)
oi.integrate(ts,ip)
figure(figsize=(8,4))
subplot(1,2,1)
oi.plot(d1='x',d2='y',gcf=True,ls='-',lw=0.5,marker='.',ms=3.)
subplot(1,2,2)
oi.plot(d1=lambda t: aAI.actionsFreqsAngles(oi(t))[6],
d2=lambda t: aAI.actionsFreqsAngles(oi(t))[8],
xlabel=r'$\theta_r$',
ylabel=r'$\theta_\theta$',
gcf=True,
ls='none',marker='.',ms=3.)
gca().set_aspect('equal')
tight_layout();
Figure 4.11: An orbit in a spherical potential in configuration space \((x,y)\) and in angle space \((\theta_r,\theta_\theta)\). Points show the orbit at equally-spaced time intervals, illustrating that angles increase linearly in time.
.
.
.