Transforming to actions, angles, and frequencies#

Introduction#

Regular orbits permit a (local) transformation to a set of canonical coordinates such that the momenta are independent, isolating integrals of motion (the actions, \(\boldsymbol{J}\)) and the conjugate coordinate variables (the angles, \(\boldsymbol{\theta}\)) linearly increase with time. Action-angle coordinates are useful for a number of applications because the equations of motion are very simple:

\[\begin{split}H &= H(\boldsymbol{J})\\ \dot{\boldsymbol{J}} &= -\frac{\partial H}{\partial \boldsymbol{\theta}} = 0\\ \dot{\boldsymbol{\theta}} &= \frac{\partial H}{\partial \boldsymbol{J}} = \boldsymbol{\Omega}(\boldsymbol{J}) = {\rm constant}\end{split}\]

Analytic transformations from phase-space to action-angle coordinates are only known for a few simple cases where the gravitational potential is separable or has many symmetries. However, astronomical systems can often be approximately axisymmetric or triaxial, or have complex radial profiles that are not captured by these simple gravitational potentials where the transformations are known.

Several numerical methods have been developed over recent years to enable approximate transformations between ordinary position and velocity to action-angle coordinates – see [sanders16] for a summary of these methods. In Gala, we have implemented the method described in [sanders14] – later in [sanders16] named the “O2GF” method – for computing actions and angles from numerically integrated orbits. Gala also implements the “Staeckel Fudge” method [binney12], which is much faster but only useful for axisymmetric or spherical potentials. In general, the O2GF method is more accurate, but slower, but we typically recommend using the Staeckel Fudge method (see below) for most use cases.

The O2GF action solver#

As mentioned above, this method was first introduced in [sanders14] and later described in [sanders16]. This method is very general in that it works with any numerically-integrated orbital time series. However, it is slower than other approximate methods: If your system is spherical or axisymmetric, other methods will perform much better. If your system is triaxial, this method is your best option. We demonstrate this method below with two qualitatively different orbits:

A tube orbit in an axisymmetric potential
A tube orbit in a triaxial potential

(see also [binneytremaine] and [mcgill90] for more context). For the examples below, we will use the galactic unit system and assume the following imports have been executed:

>>> import astropy.coordinates as coord
>>> import astropy.units as u
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> import gala.dynamics as gd
>>> import gala.integrate as gi
>>> import gala.potential as gp
>>> from gala.units import galactic

For many more options for action calculation, see tact.

A tube orbit in an axisymmetric potential#

For an example of an axisymmetric potential, we use a flattened logarithmic potential:

\[\Phi(x,y,z) = \frac{1}{2}v_{\rm c}^2\ln (x^2 + y^2 + (z/q)^2 + r_h^2)\]

with parameters

\[\begin{split}v_{\rm c} &= 150~{\rm km}~{\rm s}^{-1}\\ q &= 0.9\\ r_h &= 0\end{split}\]

For the orbit, we use initial conditions

\[\begin{split}\boldsymbol{r} &= (8, 0, 0)~{\rm kpc}\\ \boldsymbol{v} &= (75, 150, 50)~{\rm km}~{\rm s}^{-1}\end{split}\]

We first create a potential and set up our initial conditions:

>>> pot = gp.LogarithmicPotential(
...     v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
...     units=galactic)
>>> w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
...                            vel=[75, 150, 50.]*u.km/u.s)

We will now integrate the orbit and plot it in the meridional plane:

>>> w = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000)
>>> cyl = w.represent_as('cylindrical')
>>> fig = cyl.plot(['rho', 'z'], linestyle='-')

(Source code, png, pdf)

To solve for the actions in the true potential, we first compute the actions in a “toy” potential – a potential in which we can compute the actions and angles analytically. The two simplest potentials for which this is possible are the IsochronePotential and HarmonicOscillatorPotential. We will use the Isochrone potential as our toy potential for tube orbits and the harmonic oscillator for box orbits.

We start by finding the parameters of the toy potential (Isochrone in this case) by minimizing the dispersion in energy for the orbit:

>>> toy_potential = gd.fit_isochrone(w)
>>> toy_potential
<IsochronePotential: m=1.24e+11 solMass, b=4.02 kpc (kpc,Myr,solMass,rad)>

Alternatively, method="force" chooses the parameters so that the radial force of the toy potential matches that of the true potential at the minimum and maximum radius of the orbit, as recommended by Sanders & Binney (2016). This is much faster than the default energy-variance fit, but we have not found it to give more accurate actions:

>>> toy_potential_force = gd.fit_isochrone(w, method="force")

The actions and angles in this potential are not the true actions, but will only serve as an approximation. This can be seen in the angles: the orbit in the true angles would be perfectly straight lines with slope equal to the frequencies. Instead, the orbit is wobbly in the toy potential angles:

>>> toy_actions,toy_angles,toy_freqs = toy_potential.action_angle(w)
>>> fig,ax = plt.subplots(1,1,figsize=(5,5))
>>> ax.plot(toy_angles[0], toy_angles[2], linestyle='none', marker=',')
>>> ax.set_xlim(0,2*np.pi)
>>> ax.set_ylim(0,2*np.pi)
>>> ax.set_xlabel(r"$\theta_1$ [rad]")
>>> ax.set_ylabel(r"$\theta_3$ [rad]")

(Source code, png, pdf)

This can also be seen in the value of the action variables, which are not time-independent in the toy potential:

>>> fig,ax = plt.subplots(1,1)
>>> ax.plot(w.t, toy_actions[0], marker='')
>>> ax.set_xlabel(r"$t$ [Myr]")
>>> ax.set_ylabel(r"$J_1$ [rad]")

(Source code, png, pdf)

We can now find approximations to the actions in the true potential. We have to choose the maximum integer vector norm, N_max, which here we arbitrarily set to 8. This will change depending on the convergence of the action correction (the properties of the orbit and potential) and the accuracy desired:

>>> result = gd.find_actions_o2gf(w, N_max=8, toy_potential=toy_potential)
>>> result.keys()
dict_keys(['Sn', 'nvecs', 'freqs', 'dSn_dJ', 'angles', 'actions'])

The value of the actions, frequencies, and the angles at t=0 are returned in the result dictionary:

>>> result['actions']
<Quantity [ 0.12472277, 1.22725461, 0.05847431] kpc2 solMass / Myr>

To visualize how the actions are computed, we again plot the actions in the toy potential and then plot the “corrected” actions – the approximation to the actions computed using this machinery:

>>> nvecs = gd.generate_n_vectors(8, dx=1, dy=2, dz=2)
>>> act_correction = nvecs.T[...,None] * result['Sn'][None,:,None] * np.cos(nvecs.dot(toy_angles))[None]
>>> action_approx = toy_actions - 2*np.sum(act_correction, axis=1)*u.kpc**2/u.Myr
>>>
>>> fig,ax = plt.subplots(1,1)
>>> ax.plot(w.t, toy_actions[0].to(u.km/u.s*u.kpc), marker='', label='$J_1$')
>>> ax.plot(w.t, action_approx[0].to(u.km/u.s*u.kpc), marker='', label="$J_1'$")
>>> ax.set_xlabel(r"$t$ [Myr]")
>>> ax.set_ylabel(r"[kpc ${\rm M}_\odot$ km/s]")
>>> ax.legend()

(Source code, png, pdf)

Above the blue line represents the approximation of the actions in the true potential.

A tube orbit in a triaxial potential#

The same procedure works for regular orbits in more complex potentials. We demonstrate this below by repeating the above in a triaxial potential. We again use a logarithmic potential, but with flattening along two dimensions:

\[\Phi(x,y,z) = \frac{1}{2}v_{\rm c}^2\ln ((x/q_1)^2 + (y/q_2)^2 + (z/q_3)^2)\]

with parameter values:

\[\begin{split}v_{\rm c} &= 150~{\rm km}~{\rm s}^{-1}\\ q_1 &= 1\\ q_2 &= 0.9\\ q_3 &= 0.8\\ r_h &= 0\end{split}\]

and the same initial conditions as above:

\[\begin{split}\boldsymbol{r} &= (8, 0, 0)~{\rm kpc}\\ \boldsymbol{v} &= (75, 150, 50)~{\rm km}~{\rm s}^{-1}\end{split}\]

import astropy.coordinates as coord
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import gala.potential as gp
import gala.dynamics as gd
from gala.units import galactic

# define potential
pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=0.9, q3=0.8, r_h=0,
                              units=galactic)

# define initial conditions
w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
                           vel=[75, 150, 50.]*u.km/u.s)

# integrate orbit
w = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000)

# solve for toy potential parameters
toy_potential = gd.fit_isochrone(w)

# compute the actions,angles in the toy potential
toy_actions,toy_angles,toy_freqs = toy_potential.action_angle(w)

# find approximations to the actions in the true potential
import warnings
with warnings.catch_warnings(record=True):
    warnings.simplefilter("ignore")
    result = gd.find_actions_o2gf(w, N_max=8, toy_potential=toy_potential)

# for visualization, compute the action correction used to transform the
#   toy potential actions to the approximate true potential actions
nvecs = gd.generate_n_vectors(8, dx=1, dy=2, dz=2)
act_correction = nvecs.T[...,None] * result['Sn'][0][None,:,None] * np.cos(nvecs.dot(toy_angles))[None]
action_approx = toy_actions - 2*np.sum(act_correction, axis=1)*u.kpc**2/u.Myr

fig,axes = plt.subplots(3,1,figsize=(6,14))

for i,ax in enumerate(axes):
    ax.plot(w.t, toy_actions[i].to(u.km/u.s*u.kpc), marker='', label='$J_{}$'.format(i+1))
    ax.plot(w.t, action_approx[i].to(u.km/u.s*u.kpc), marker='', label="$J_{}'$".format(i+1))
    ax.set_ylabel(r"[kpc ${\rm M}_\odot$ km/s]")
    ax.legend(loc='upper left')

ax.set_xlabel(r"$t$ [Myr]")
fig.tight_layout()

(Source code, png, pdf)

The Staeckel Fudge#

find_actions_staeckel_fudge() computes approximate actions, angles, and frequencies in axisymmetric potentials with the “Staeckel Fudge” [binney12], using the expressions for angles and frequencies from [sanders12]. See also the Staeckel Fudge implementation in Agama [vasiliev19]. This method is only applicable for axisymmetric systems, but is much faster than the O2GF method, and works directly on phase-space positions (i.e. no orbit integration is needed). The focal length parameter, Δ, of the prolate spheroidal coordinate system can be passed in with the delta argument; otherwise it is estimated with get_staeckel_fudge_delta():

>>> from gala.dynamics.actionangle import find_actions_staeckel_fudge
>>> mw = gp.MilkyWayPotential()
>>> w = gd.PhaseSpacePosition(pos=[8.0, 0, 0.5] * u.kpc,
...                           vel=[30, 220, 40] * u.km / u.s)
>>> aaf = find_actions_staeckel_fudge(mw, w)
>>> aaf["actions"]

As a more complete example, below we will compute vertical actions for a set of phase-space positions that have a range of initial vertical velocities and are orbiting in a two-component model for a galactic potential (a disk + halo model).

We will start by defining the potential model (this is just a simple example):

>>> halo = gp.NFWPotential.from_M200_c(
...     M200=1e12*u.Msun, c=15,
...     units=galactic
... )
>>> disk = gp.MN3ExponentialDiskPotential(
...     m=8e10*u.Msun, h_R=3.5*u.kpc, h_z=0.4*u.kpc,
...     units=galactic
... )
>>> pot = halo + disk

We next define a grid of orbital initial conditions with close to the circular velocity but varying vertical velocities:

>>> vcirc = pot.circular_velocity([8, 0, 0])
>>> vz_grid = np.linspace(0.5, 200, 64) * u.km/u.s
>>> xyz = np.repeat([[8., 0, 0]], len(vz_grid), axis=0).T * u.kpc
>>> vxyz = np.repeat([[0, 1.1, 0]], len(vz_grid), axis=0).T * vcirc
>>> vxyz[2] = vz_grid
>>> w0 = gd.PhaseSpacePosition(xyz, vxyz)

To compute the actions, angles, and frequencies for these phase-space positions, we can use the Staeckel Fudge:

>>> aaf = find_actions_staeckel_fudge(pot, w0)
>>> J = aaf["actions"]

Since we only varied the vertical velocity, we expect there to be a smooth and simple relationship between the vertical action and the initial vertical velocity. We can visualize this below:

>>> fig, ax = plt.subplots(figsize=(6, 6), constrained_layout=True)
>>> ax.plot(w0.v_z, J[:, 2].to(u.kpc * u.km / u.s))
>>> ax.set_xlabel(f"$v_z$ [{w0.v_z.unit:latex_inline}]")
>>> ax.set_ylabel(r"$J_z$ [kpc km/s]")

(Source code)

One limitation (or feature, if you’re Binney) of the Staeckel Fudge is that it effectively interpolates over resonances in the true potential. This means that the Staeckel Fudge will return global actions, angles, and frequencies even for orbits that are trapped in resonances. In these cases, the returned values may be misleading. To illustrate this, we’ll compute orbits for the initial conditions above and compute the approximate actions, angles, and frequencies along each orbit:

>>> orbits = pot.integrate_orbit(
...     w0, dt=1, t1=0, t2=6*u.Gyr,
...     Integrator='dopri853'
... )
>>> orbits.cylindrical.plot(['rho', 'z'], alpha=0.5, marker=',')

(Source code, png, pdf)

Now the actions - here, we’ll compute the mean actions over the orbit and make the same plot of vertical velocity vs. vertical action but now with the mean values:

>>> aaf_orbits = find_actions_staeckel_fudge(pot, orbits, mean=False)
>>> mean_J = np.mean(aaf_orbits["actions"], axis=0)
>>> plt.plot(w0.v_z, mean_J[:, 2].to(u.kpc * u.km / u.s))

(Source code)

The overall trend looks right, but what is that weird break that occurs around \(v_z\) ~ 120 km/s?! Let’s visualize orbits with initial conditions on either side of this value:

(Source code)

Aha – this region is special. It is a resonance in the potential! Orbits in this region of phase-space have qualitatively different behavior than those outside of this region because they are trapped by the resonance. However, Staeckel potentials do not have resonances (are globally integrable) and so, when we use the Staeckel Fudge on phase-space positions only, we gloss over the resonance structure. This is something to be aware of when using the Staeckel Fudge!