# Action-Angle Variables for Axisymmetric Potentials via Birkhoff
Normalization

Sam Hadden Canadian Institute for Theoretical Astrophysics, 60 St George St
Toronto, ON M5S 3H8, Canada

(April 2024)

###### Abstract

We describe a method for calculating action-angle variables in axisymmetric
galactic potentials using Birkhoff normalization, a technique from Hamiltonian
perturbation theory. An advantageous feature of this method is that it yields
explicit series expressions for both the forward and inverse transformations
between the action-angle variables and position-velocity data. It also
provides explicit expressions for the Hamiltonian and dynamical frequencies as
functions of the action variables. We test this method by examining orbits in
a Miyamoto-Nagai model potential and compare it to the popular Stäckel
approximation method. When vertical actions are not too large, the Birkhoff
normalization method achieves fractional errors smaller than a part in
$10^{3}$ and outperforms the Stäckel approximation. We also show that the
range over which Birkhoff normalization provides accurate results can be
extended by constructing Padé approximants from the perturbative series
expressions developed with the method. Numerical routines in Python for
carrying out the Birkhoff normalization procedure are made available.

##  1 Introduction

While a comprehensive dynamical description of a galaxy must take into account
mutual gravitational interactions among stars, dark matter, and gas,
dramatically simplified models that consider the orbits of stars in static,
smooth potentials provide surprisingly powerful tools for understanding the
dynamics of these complex systems. Numerical experiments show that typical
stellar orbits in model galactic potentials execute quasi-periodic motions
(e.g., Binney & Spergel, 1982, 1984). Such motions are restricted to
3-dimensional tori embedded in the full 6-dimensional phase space spanned by
stars’ spatial positions and velocities. For completely integrable Hamiltonian
systems, action-angle (AA) variables provide a particularly convenient set of
canonical coordinates for describing motion on such tori (e.g., Goldstein et
al., 2002). Fixing the values of three action variables, $\boldsymbol{J}$,
specifies an individual torus as a 3-dimensional submanifold of the
6-dimensional phase space, while the canonically conjugate angle variables,
$\boldsymbol{\theta}$, serve as coordinates on the torus. Furthermore, the
time evolution of trajectories on a given torus are described simply by the
linear advance of the angle variables with time at rates
$\frac{d}{dt}\boldsymbol{\theta}=\nabla_{\boldsymbol{J}}H$, where $H$ is the
Hamiltonian function.

While AA variables are a powerful theoretical tool, it is seldom possible to
compute them explicitly from stellar position and velocity data for potentials
commonly adopted when modeling our Milky Way and other galaxies. This is
because these potentials generally do not admit a set of three globally-
defined isolating integral which are required in order to have a well-defined
transformation to AA variables everywhere in phase space. Nonetheless, the
quasi-periodic nature of stellar orbits over much of the phase space of
interest in galactic potentials makes its possible to construct local
transformations to AA variables via numerical means.

Numerous numerical methods have been developed for the purposes of
constructing AA variables in galactic potentials (see review by Sanders &
Binney, 2016). One of the most widely-used methods is the Stäckel
approximation or “Stäckel fudge”, originally proposed by Binney (2012). In
axisymmetric potentials, this method allows radial and vertical actions to be
estimated via numerical quadratures by implicitly treating the underlying
potential as though it were in Stäckel form, meaning separable in prolate
spheroidal coordinates. While the Stäckel approximation usually provides quite
accurate determinations of AA variables, it does have some shortcomings.
First, the transformation cannot be readily inverted to provide positions and
velocities as functions of actions and angle nor express the original
Hamiltonian in terms of AA variables. Furthermore, because the method requires
evaluation of integrals by numerical quadratures, it can become
computationally expensive when evaluating transformations for large amounts of
data.

In this paper, we present a perturbative method to determine action-angle
variables in axisymmetric galactic potentials. This method is based on a
method originally proposed by Birkhoff (1927) for developing successively
refined approximations to the motion of Hamiltonian systems in the vicinity of
elliptic equilibria. It is variously referred to as to as “Birkhoff
normalization” (e.g., Deprit, 1969) or the “Birkhoff-Gustavson method” (e.g.,
Lowenstein, 2012). One of the principal advantages of this method is that it
provides explicit expressions for the transformation from position and
velocity data to AA variables. Furthermore, this method simultaneously yields
explicit expressions for the inverse transformation, meaning stellar positions
and velocities as well as the Hamiltonian can be represented explicitly in
terms of the AA variables. Expressions are constructed as multi-variate
polynomials in complex canonical coordinates, introduced in Section 2, and
therefore can be evaluated at minimal computational cost. Together, these
features can make Birkhoff normalization an attractive alternative to the
Stäckel approximation, depending on the application.

The use of perturbation theory in galactic dynamics to derive approximate
integrals of motion is not new (e.g., Contopoulos, 1963; Saaf, 1968; de Zeeuw
& Merritt, 1983). In fact, Gustavson (1966) implements and applies an
algorithm very similar to the one presented here to derive a formal second
integral of the Hénon–Heiles potential (Henon & Heiles, 1964). The current
paper builds on this past work in at least a couple respects. First, we follow
Deprit (1969) and develop perturbative series using the Lie transform method.
This method allows any functions defined on phase space in terms of the
original canonical variables to be readily expressed in transformed variables
and vice versa. This fact is exploited to construct an explicit expression for
the canonical angle variable, $\theta_{\phi}$, associated with orbits’
azimuthal degree of freedom, thereby providing a complete transformation from
6D position and velocity data to AA variables. Second, we show how Padé
approximants can be deployed to extend the range of validity of the series
expressions furnished by the method. Finally, a modern computational
implementation of the method in Python is made publicly
available.111github.com/shadden/AA_variables_via_Birkhoff_Normalization The
availability of computer algebra and automatic differentiation packages in
modern languages like Python make it feasible to apply Birkhoff normalization
to more realistic potential models than the simple Hénon–Heiles model
originally considered by Gustavson (1966).

The plan of this paper is as follows: we introduce the algorithm in Section 2.
In Section 3, we demonstrate the method with an example application to stellar
orbits in Miyamoto-Nagai potential (Miyamoto & Nagai, 1975), a classic model
for the potential of a galactic disk. Finally, we summarize and discuss our
results in Section 4.

##  2 The Algorithm

###  2.1 Setting

We consider orbits in an axisymmetric potential, $\Phi$. The Hamiltonian,
expressed in terms of polar coordinates $(R,z,\phi)$ along with their
canonically conjugate momenta, $(p_{R},p_{z},L)$, is

$\mathcal{H}(R,z,\phi,p_{R},p_{z},L)=\frac{1}{2}\left(p_{R}^{2}+p_{z}^{2}+\frac{L^{2}}{R^{2}}\right)+\Phi(R,z)~{}.$
(1)

Since the Hamiltonian in (1) has no explicit $\phi$ dependence, the angular
momentum, $L$, is conserved and we can consider the dynamics in the reduced
phase space comprised of $(z,R,p_{z},p_{R})$, treating $L$ as a parameter. The
Hamiltonian governing the dynamics in this reduced phase space is then given
by

${H}(R,z,\phi,p_{R},p_{z};L)=\frac{1}{2}\left(p_{R}^{2}+p_{z}^{2}\right)+{\Phi}_{\mathrm{eff}}(R,z;L)$
(2)

where we have defined
${\Phi}_{\mathrm{eff}}={\Phi}({R,z})+\frac{L^{2}}{2R^{2}}$ as the effective
potential. For a given value of angular momentum, $L$, there is a (typically
unique) circular orbit of radius of $R_{C}(L)$ which satisfies

$\frac{\partial}{\partial
R}{\Phi}_{\mathrm{eff}}({R,z})\bigg{|}_{(R,z)=(R_{C},0)}=0~{}.$ (3)

For nearly circular and planar orbits, the Hamiltonian (2) is given
approximately by

$H\approx
H_{2}(z,R,p_{R},p_{z};L):=\frac{1}{2}\left(p_{R}^{2}+\kappa^{2}\delta
R^{2}\right)+\frac{1}{2}\left(p_{z}^{2}+\nu^{2}z^{2}\right)~{}.$ (4)

where $\kappa^{2}=\frac{\partial^{2}}{\partial R^{2}}{\Phi}_{\mathrm{eff}}$
and $\nu^{2}=\frac{\partial^{2}}{\partial z^{2}}{\Phi}_{\mathrm{eff}}$, with
derivatives evaluated at $(R,z)=(R_{C},0)$, are the epicyclic and vertical
frequencies, respectively and $\delta R=R-R_{C}$. Equation (4) is, of course,
the Hamiltonian of two un-coupled harmonic oscillators of frequency $\kappa$
and $\nu$ and the equations of motion derived from it provide the so-called
“epicyclic approximation” for stars’ orbits (e.g., Binney & Tremaine, 2008).

###  2.2 Complex canonical variables

The equations of motion derived from Equation (4) are completely integrable
and we can perform a canonical transformation,
$T:(R,z,p_{R},p_{z})\mapsto(\phi_{R},\phi_{z},I_{R},I_{z})$, to new variables
such that $H_{2}\circ T^{-1}=\kappa I_{R}+\nu I_{z}$. Throughout this paper,
we will frequently work instead with the associated complex canonical
variables

$\displaystyle x_{R}$ $\displaystyle=\sqrt{\frac{\kappa}{2}}\left(\delta
R+\mathrm{i}\frac{p_{R}}{\kappa}\right)=\sqrt{I_{R}}e^{-\mathrm{i}\phi_{R}}$
(5) $\displaystyle x_{z}$
$\displaystyle=\sqrt{\frac{\nu}{2}}\left(z+\mathrm{i}\frac{p_{z}}{\nu}\right)=\sqrt{I_{z}}e^{-\mathrm{i}\phi_{z}}~{}$
(6)

and their complex conjugates, $\bar{x}_{R}$ and $\bar{x}_{z}$. The Poisson
bracket of two functions, $f$ and $g$, of phase space coordinates is expressed
in terms of these variables by

$[f,g]=-\mathrm{i}\sum_{j}\left(\frac{\partial
f}{\partial{x}_{j}}\frac{\partial g}{\partial\bar{x}_{j}}-\frac{\partial
f}{\partial\bar{x}_{j}}\frac{\partial g}{\partial{x}_{j}}\right)$ (7)

where the subscript $j$ denotes either $R$ or $z$. In terms of the complex
canonical variables,222 Here and elsewhere we will make a slight abuse of
notation and use the same symbol to denote the Hamiltonian function expressed
in terms of complex variables as we do for the Hamiltonian expressed in terms
of polar coordinates. Equation (4) becomes $H_{2}=\kappa x_{R}\bar{x}_{R}+\nu
x_{z}\bar{x}_{z}$ and Hamilton’s equations give
$\frac{d}{dt}x_{R}=[x_{R},H_{2}]=-\mathrm{i}\kappa x_{R}$ and
$\frac{d}{dt}x_{z}=[x_{R},H_{2}]=-\mathrm{i}\nu x_{z}$. Thus, under the
epicyclic approximation, the complex canonical variables evolve at constant
magnitude with a uniformly rotating phase. This will generally no longer be
true when higher order terms are included in the Hamiltonian. Our goal will be
to develop an algorithm that constructs a canonical transformation to new
complex canonical variables, $x_{R}^{\prime}$ and $x^{\prime}_{z}$, for which
this is again true up to some specified order in the new variables and their
complex conjugates.

###  2.3 An illustrative example

Before giving the general algorithm, we first demonstrate its basic principles
by way of an example application to the pendulum. This example provides a
pedagogical introduction to Birkhoff normalization and the construction of
canonical transformations via the Lie series method with complex canonical
variables. The example illustrates the basic ideas behind the method while
avoiding some of the cumbersome mathematical notation and manipulations
required to formulate the procedure in a more general context. Readers
uninterested in the technical details involved in an algorithmic
implementation of the Birkhoff normalization procedure can therefore read this
section while skipping or skimming Section 2.4.

We start with the Hamiltonian governing the dynamics of a simple pendulum,

$h(\theta,p)=\frac{p^{2}}{2}-\omega_{0}^{2}\cos{\theta}.$ (8)

For small $\theta$, the truncated Hamiltonian,
$h_{2}(\theta,p)=(p^{2}+\omega_{0}^{2}\theta^{2})/2$, gives a harmonic
oscillator approximation of the dynamics. In terms of the complex canonical
variable $x=\sqrt{\frac{\omega_{0}}{2}}(\theta+\mathrm{i}p/\omega_{0})$, the
Hamiltonian (8) becomes

$h(x,\bar{x})=\omega_{0}x\bar{x}+\sum_{n=2}^{\infty}h_{2n}(x,\bar{x})$ (9)

where

$h_{2n}=\frac{(-1)^{n+1}\omega_{0}^{2}}{(2n)!}\left(\frac{1}{2\omega_{0}}\right)^{n}\left(x+\bar{x}\right)^{2n}~{}.$
(10)

We will seek a canonical transformation to a new complex canonical variable,
$x^{\prime}$, so that the transformed Hamiltonian, $h^{\prime}$, depends only
on $|x^{\prime}|^{2}$ up to a specified power, $N$, in the transformed
variable, $x^{\prime}$, and its complex conjugate, $\bar{x}^{\prime}$.

The Lie series method (e.g., Ferraz-Mello, 2007; Lichtenberg & Lieberman,
1992) provides a particularly convenient means of constructing such a
transformation. This method involves the construction of a generating
function, $\chi$, such that the transformation from the new variables to the
old variables is given by the time-1 flow of the Hamiltonian vector field
generated by $\chi$. The time derivative of a function, $f$, under this flow
is given by the Poisson bracket $[f,\chi]$ and we define the differential
operator $\mathcal{L}_{\chi}:=[\cdot,\chi]$, referred to as the Lie
derivative. This allows us to express a function, $f$, given in terms of the
original phase space coordinates $(x,\bar{x})$, in terms of the new phase
space variables, $(x^{\prime},\bar{x}^{\prime})$, formally as

$f(x,\bar{x})=(\exp[\mathcal{L}_{\chi}]f)(x^{\prime},\bar{x}^{\prime})=\left(\sum_{n=0}^{\infty}\frac{1}{n!}\mathcal{L}_{\chi}^{n}f\right)(x^{\prime},\bar{x}^{\prime})~{}.$
(11)

Returning to our pendulum example, we seek a generating function $\chi$ such
that
$h^{\prime}(x^{\prime},\bar{x}^{\prime})=(\exp[\mathcal{L}_{\chi}]h)(x^{\prime},\bar{x}^{\prime})=\sum_{n=1}^{N/2}h^{\prime}_{2n}(|x^{\prime}|^{2})+R_{N}(x^{\prime},\bar{x}^{\prime})$
where the remainder terms
$R_{N}(x^{\prime},\bar{x^{\prime}})\sim\mathcal{O}(|x^{\prime}|^{N+1})$. Our
approach will be to write our generating function as
$\chi=\sum_{n=2}\chi_{2n}$ where
$\chi_{2n}\sim{\mathcal{O}(|x^{\prime}|^{2n})}$. This will allow us to develop
an iterative procedure for determining $\chi_{2n}$ and $h^{\prime}_{2n}$ at
each stage in terms of the functions $h_{2},h_{4}...,h_{2n}$ defined in
Equation (10) and the functions $\chi_{4},\chi_{6},...,\chi_{2n-2}$ determined
in previous stages. Writing out the transformation Equation (11) explicitly in
terms of Poisson brackets, we obtain

$\displaystyle\exp[\mathcal{L}_{\chi}]h$
$\displaystyle=h+[h,\chi]+\frac{1}{2}[[h,\chi],\chi]+...$ (12)
$\displaystyle=\underbrace{[h_{2},\chi_{4}]+h_{4}}_{h^{\prime}_{4}}+\underbrace{[h_{2},\chi_{6}]+[h_{4},\chi_{4}]+\frac{1}{2}[[h_{2},\chi_{4}],\chi_{4}]+h_{6}}_{h^{\prime}_{6}}+\mathcal{O}(|x^{\prime}|^{8})$
(13)

where, in the second line, we have grouped terms by common powers of
$|x^{\prime}|$. Consider the equation

$h^{\prime}_{4}=[h_{2},\chi_{4}]-\left(\frac{x^{\prime 4}}{96}+\frac{x^{\prime
3}\bar{x}^{\prime}}{24}+\frac{x^{\prime 2}\bar{x}^{\prime
2}}{16}+\frac{x^{\prime}\bar{x}^{\prime 3}}{24}+\frac{\bar{x}^{\prime
4}}{96}\right)$ (14)

defining the fourth order term in the transformed Hamiltonian, where we have
used Equation (10) to write out $h_{4}$ explicitly. We wish to define
$\chi_{4}$ so that $h^{\prime}_{4}$ depends only on $|x^{\prime}|^{2}$. Notice
that the Poisson bracket between $h_{2}$ and any monomial, $x^{\prime
k}\bar{x}^{\prime\bar{k}}$, gives $[h_{2},x^{\prime
k}\bar{x^{\prime}}^{\bar{k}}]=\mathrm{i}\omega_{0}(k-\bar{k})x^{\prime
k}\bar{x}^{\prime\bar{k}}$. Therefore, if we choose $\chi_{4}$ as the sum of
monomial terms

$\displaystyle\chi_{4}=\frac{1}{\mathrm{i}\omega_{0}}\left(\frac{x^{\prime
3}\bar{x}^{\prime}}{48}-\frac{x^{\prime}\bar{x}^{\prime
3}}{48}-\frac{\bar{x}^{\prime 4}}{384}+\frac{x^{\prime 4}}{384}\right)$ (15)

then each term in the polynomial $[h_{2},\chi_{4}]$ cancels a corresponding
term appearing in Equation (14) except the term $\propto x^{\prime
2}\bar{x}^{\prime 2}$. We therefore have

$h^{\prime}_{4}=-\frac{1}{16}x^{\prime 2}\bar{x}^{\prime 2}~{},$ (16)

which depends only on $|x^{\prime}|^{2}$ as desired.

With $\chi_{4}$ determine, we can write $h^{\prime}_{6}$ according to Equation
(13) as

$h^{\prime}_{6}=[h_{2},\chi_{6}]+\frac{x^{\prime
6}}{2560\omega_{0}}-\frac{x^{\prime
5}\bar{x}^{\prime}}{3840\omega_{0}}-\frac{7x^{\prime 4}\bar{x}^{\prime
2}}{1536\omega_{0}}-\frac{x^{\prime 3}\bar{x}^{\prime
3}}{256\omega_{0}}-\frac{7x^{\prime 2}\bar{x}^{\prime
4}}{1536\omega_{0}}-\frac{x^{\prime}\bar{x}^{\prime
5}}{3840\omega_{0}}+\frac{\bar{x}^{\prime 6}}{2560\omega_{0}}~{}.$ (17)

We will again choose $\chi_{6}$ so that the Poisson bracket $[h_{2},\chi_{6}]$
cancels all of the monomial terms appearing on the right-hand side of Equation
(17) except the term $\propto x^{3}\bar{x}^{3}$. Proceeding, we find

$\displaystyle\chi_{6}$
$\displaystyle=\frac{1}{\mathrm{i}\omega_{0}^{2}}\left(-\frac{x^{\prime
6}}{15360}+\frac{x^{\prime 5}\bar{x}^{\prime}}{15360}+\frac{7x^{\prime
4}\bar{x}^{\prime 2}}{3072}-\frac{7x^{\prime 2}\bar{x}^{\prime
4}}{3072}-\frac{x^{\prime}\bar{x}^{\prime 5}}{15360}+\frac{\bar{x}^{\prime
6}}{15360}\right)~{}.$ (18) $\displaystyle h^{\prime}_{6}$
$\displaystyle=-\frac{x^{\prime 3}\bar{x}^{\prime 3}}{256\omega_{0}}~{}.$ (19)

If we carried on with the expansion of Equation (13) to higher order, we would
find that at each step, we need to solve an equation of the form

$h^{\prime}_{2n}=[h_{2},\chi_{2n}]+\Psi_{2n}$ (20)

where $\Psi_{2n}$ is a collection of terms of order $2n$ involving Poisson
brackets of the known functions $h_{4},h_{6},...,h_{2n}$ and
$\chi_{4},...,\chi_{2n-2}$. In the next section we will detail an algorithm
for iteratively constructing and solving _homological equations_ like Equation
(20).

Let us conclude our present example by approximating the pendulum’s
oscillation frequency as a function of maximum angular displacement,
$\theta_{\mathrm{max}}$, in order to illustrate how functions of the new,
transformed variables, $x^{\prime}$ and $\bar{x}^{\prime}$, may be calculated
explicitly from the original coordinates and vice versa. Our transformed
Hamiltonian is given by

$h^{\prime}=\omega_{0}|x^{\prime}|^{2}-\frac{1}{16}|x^{\prime}|^{4}-\frac{1}{256\omega_{0}}|x^{\prime}|^{6}+\mathcal{O}(|x^{\prime}|^{8})$
(21)

so the oscillation frequency will be given, in terms of the action
$J=|x^{\prime}|^{2}$, by
$\omega(J)=\omega_{0}-\frac{1}{8}J-\frac{3}{256\omega_{0}}J^{2}$. We therefore
need the value of $J$ in terms of $\theta_{\mathrm{max}}$.

Given any function, $f$, on the phase space expressed in terms of the
transformed variables, $x^{\prime}$ and $\bar{x}^{\prime}$, its expression in
terms of the original, un-primed variables $x$ and $\bar{x}$ is given by

$f\left(x^{\prime}(x,\bar{x}),\bar{x}^{\prime}(x,\bar{x})\right)=(\exp[-\mathcal{L}_{\chi}]f)(x,\bar{x})$
(22)

since $\exp[-\mathcal{L}_{\chi}]f$ gives the function $f$ evaluated at the
time-($-1$) flow of the Hamiltonian vector field generated by $\chi$. Applying
Equation (22) to the function
$\omega:(x^{\prime},\bar{x}^{\prime})\mapsto\omega_{0}-\frac{1}{8}|x^{\prime}|^{2}-\frac{3}{256\omega_{0}}|x^{\prime}|^{4}$,
we obtain the oscillation frequency in terms of the original, un-transformed
complex variable
$x=\sqrt{\frac{\omega_{0}}{2}}(\theta+\mathrm{i}p/\omega_{0})$ as

$\displaystyle(\exp[-\mathcal{L}_{\chi}]\omega)(x,\bar{x})$
$\displaystyle=\omega_{0}-\frac{|x|^{2}}{8}-\frac{3}{256}|x|^{4}+\frac{1}{8}\left[|x|^{2},\chi_{4}\right]+\mathcal{O}(|x|^{6})$
(23)
$\displaystyle=\omega_{0}-\frac{x\bar{x}}{8}+\frac{1}{\omega_{0}}\left(\frac{\bar{x}^{4}}{768}+\frac{\bar{x}^{3}x}{192}-\frac{3\bar{x}^{2}x^{2}}{256}+\frac{\bar{x}x^{3}}{192}+\frac{x^{4}}{768}\right)+\mathcal{O}(|x|^{6})~{}.$
(24)

Substituting $x=\sqrt{\frac{\omega_{0}}{2}}\theta_{\mathrm{max}}$ into
Equation (23) gives

$\omega\approx\omega_{0}\left(1-\frac{\theta_{\max}^{2}}{16}+\frac{\theta_{\max}^{4}}{3072}+...\right),$
(25)

which matches the first three terms of the Taylor expansion of the exact
result,
$\omega=\frac{\pi\omega_{0}}{2\mathbb{K}\left(\sin^{2}\left(\frac{\theta_{\max}}{2}\right)\right)}$,
where $\mathbb{K}$ denotes the complete elliptic integral of the first kind
(e.g., Lichtenberg & Lieberman, 1992).

###  2.4 The general algorithm

The Birkhoff normalization procedure, illustrated above in the case of the
pendulum, readily generalizes to systems with $d>1$ degrees of freedom.
Specifically, the procedure can be applied in the vicinity of an elliptic
equilibrium to construct a formal series for a generating function, $\chi$,
that transforms the Hamiltonian so that it only depends on $d$ action-like
variables. While these formal series are, in general, divergent, their finite-
order truncations are nonetheless often very useful for approximating a given
system’s dynamics (e.g., Contopoulos, 2003; Efthymiopoulos et al., 2004). The
numerical results presented in later in Section 3 show that the approximate
integrals computed for our chosen order of truncation are very nearly
constants over a large region of the phase space.

We apply the Birkhoff normalization procedure to orbits in the vicinity of
planar, circular orbits of axisymmetric galactic potentials, which are
elliptic equilibria of the two degree-of-freedom Hamiltonian given in Equation
(2). In terms of the complex canonical variables $x_{R}$ and $x_{z}$
introduced in Section 2.2, the Hamiltonian can be written as

$\displaystyle H=\kappa x_{R}\bar{x}_{R}+\nu\kappa
x_{z}\bar{x}_{z}+\sum_{n=3}^{\infty}H_{n}(x_{R},x_{z},\bar{x}_{R},\bar{x}_{z})$
(26)

where

$\displaystyle
H_{n}(x_{R},x_{z},\bar{x}_{R},\bar{x}_{z})=\sum_{m=0}^{n}\frac{\Phi_{\mathrm{eff}}^{(m,n-m)}}{(n-m)!m!}\left(\sqrt{\frac{1}{2\kappa}}\right)^{n-m}\left(\sqrt{\frac{1}{2\nu}}\right)^{m}\sum_{p=0}^{m}\sum_{q=0}^{n-m}\binom{m}{p}\binom{n-m}{q}x_{R}^{p}x_{z}^{q}\bar{x}_{R}^{m-p}\bar{x}_{z}^{n-m-q}$
(27)

and $\Phi_{\mathrm{eff}}^{(m,n-m)}=\frac{\partial^{n}}{\partial
R^{m}z^{n-m}}\Phi_{\mathrm{eff}}\big{|}_{(R,z)=(R_{C},0)}$. Values of the
partial derivatives of the effective potential,
$\Phi_{\mathrm{eff}}^{(m,n-m)}$, up to the desired maximum order of the
perturbative expansion are required to carry out the Birkhoff normalization
procedure. Calculation of these derivatives is most conveniently done with the
aid of a computer algebra system or automatic differentiation packages.

Equation (26) represents the Hamiltonian as a series of terms grouped by their
order in the complex canonical variables. We will construct similar series for
our generating function, $\chi=\sum_{k=k_{\mathrm{min}}}^{\infty}\chi_{k}$,
and the transformed Hamiltonian,
$H^{\prime}=\sum_{m=1}^{\infty}H^{\prime}_{2m}$ through an iterative process.
Note that because $H^{\prime}$ depends only on the action variables, it only
contains even-order terms.

####  2.4.1 Grouping Lie Series Terms by Order

To construct our series solution, it will be necessary to group terms of the
same order in the expansion of the Lie series
$H^{\prime}=\exp[\mathcal{L}_{\chi}]H$. This is most readily accomplished
using recursion formulae originally derived by Deprit (1969). Below, we give a
derivation of these formulae that closely follows the presentation of Ferraz-
Mello (2007, Chapter 5).

If $f_{k}$ and $g_{l}$ are homogeneous polynomials of degree $k$ and $l$,
respectively, in complex canonical variables then $[f_{k},g_{l}]$ is a
homogeneous polynomial of degree $k+l-2$ in those variables. Define
$D_{k}:=[\cdot,\chi_{k}]$ so that, given series expansions
$f=\sum_{k^{\prime}=k^{\prime}_{\mathrm{min}}}^{\infty}f_{k^{\prime}}$, and
$\chi=\sum_{k=k_{\mathrm{min}}}^{\infty}\chi_{k}$, we can write

$\displaystyle\mathcal{L}_{\chi}f$
$\displaystyle=\sum_{k=k_{\mathrm{min}}}^{\infty}\sum_{k^{\prime}=k^{\prime}_{\mathrm{min}}}^{\infty}D_{k}f_{k^{\prime}}=\sum_{l=k_{\mathrm{min}}+k^{\prime}_{\mathrm{min}}-2}^{\infty}\sum_{k=k_{\mathrm{min}}}^{l+2-k^{\prime}_{\mathrm{min}}}D_{k}f_{l+2-k}:=\sum_{l=l_{\mathrm{min}}^{(1)}}^{\infty}\Upsilon^{1}_{l}(f,\chi)$

where $l_{\mathrm{min}}^{(1)}=k_{\mathrm{min}}+k^{\prime}_{\mathrm{min}}-2$
and $\Upsilon^{1}_{l}(f,\chi)$ represents all terms of order $l$ appearing in
an expansion of $\mathcal{L}_{\chi}f$. If we now substitute
$\mathcal{L}_{\chi}f$ for $f$ in Equation (LABEL:eq:Df_defn), we obtain

$\mathcal{L}^{2}_{\chi}f=\sum_{l=k_{\mathrm{min}}+l_{\mathrm{min}}^{(1)}-2}^{\infty}\sum_{k=k_{\mathrm{min}}}^{l+2-l_{\mathrm{min}}^{(1)}}D_{k}\Upsilon^{1}_{l+2-k}(f,\chi):=\sum_{l=l_{\mathrm{min}}^{(2)}}^{\infty}\Upsilon^{2}_{l}(f,\chi)$
(29)

where $l_{\mathrm{min}}^{(2)}=k_{\mathrm{min}}+l^{(1)}_{\mathrm{min}}-2$. More
generally, we have

$\displaystyle l_{\mathrm{min}}^{(n)}$
$\displaystyle=k_{\mathrm{min}}+l^{(n-1)}_{\mathrm{min}}-2$ (30)
$\displaystyle\mathcal{L}^{n}_{\chi}f$
$\displaystyle=\sum_{l=l_{\mathrm{min}}^{(n)}}^{\infty}\sum_{k=k_{\mathrm{min}}}^{l+2-l_{\mathrm{min}}^{\mathrm{(n-1)}}}D_{k}\Upsilon^{(n-1)}_{l+2-k}(f,\chi):=\sum_{l=l_{\mathrm{min}}^{(n)}}^{\infty}\Upsilon^{n}_{l}(f,\chi)$
(31)

In other words, $\Upsilon^{n}_{l}$ is the collection of order $l$ terms
arising from the $n$th Lie derivative in our expansion. Values of
$\Upsilon^{n}_{l}$ for a given $n$ and $l$ can be computed using the recursion
formulae

$\displaystyle\Upsilon^{0}_{l}(f,\chi)$ $\displaystyle=f_{l}$ (32)
$\displaystyle l_{\mathrm{min}}^{(0)}$
$\displaystyle=k^{\prime}_{\mathrm{min}}$ (33)
$\displaystyle\Upsilon^{n+1}_{l}(f,\chi)$
$\displaystyle=\sum_{k=k_{\mathrm{min}}}^{l+2-l^{(n)}_{\mathrm{min}}}D_{k}\Upsilon^{n}_{l+2-k}(f,\chi)$
(34) $\displaystyle l_{\mathrm{min}}^{(n+1)}$
$\displaystyle=k_{\mathrm{min}}+l^{(n)}_{\mathrm{min}}-2~{}.$ (35)

####  2.4.2 Expansion of the Transformed Hamiltonian

We can now write the transformed Hamiltonian,
$\exp[\mathcal{L}_{\chi}]H=\sum_{n=0}^{\infty}\frac{1}{n!}{\mathcal{L}^{n}_{\chi}}H$,
in terms of $\Upsilon_{l}^{n}(H,\chi)$. Since the lowest order terms appearing
in $H$ are order $k^{\prime}_{\mathrm{min}}=2$, Equation (35) becomes
$l_{\mathrm{min}}^{(n)}=2+(k_{\mathrm{min}}-2)n$. We have

$\exp[\mathcal{L}_{\chi}]H=\sum_{n=0}^{\infty}\sum_{l=l_{\mathrm{min}}^{(n)}}^{\infty}\frac{1}{n!}\Upsilon^{n}_{l}(H,\chi)=\sum_{l=2}^{\infty}\sum_{n=0}^{N^{(l)}}\frac{1}{n!}\Upsilon^{n}_{l}(H,\chi):=\sum_{l=2}^{\infty}(\Psi_{l}+[H_{2},\chi_{l}])$
(36)

where we define $N^{(l)}=\lfloor\frac{l-2}{k_{\mathrm{min}}-2}\rfloor$ and

$\Psi_{l}=H_{l}+\sum_{k=k_{\mathrm{min}}}^{l-1}D_{k}H_{l+2-k}+\sum_{n=2}^{N^{(l)}}\frac{1}{n!}\Upsilon^{n}_{l}(H,\chi)~{}.$
(37)

The expression in Equation (37) defining $\Psi_{l}$ only contains terms,
$\chi_{k}$, in the expansion of the generating function $\chi$ with $k<l$. The
terms of order $l$ gathered together in the sum in Equation (36) are, by
definition, the terms of order $l$ in the transformed Hamiltonian,
$H^{\prime}$. In other words, we have

$\Psi_{l}+[H_{2},\chi_{l}]=H^{\prime}_{l}~{}.$ (38)

Since $\Psi_{l}$ contains only terms $\chi_{k}$ with $k<l$, we can iterate
over $l$, solving for $H^{\prime}_{l}$ and $\chi_{l}$ at each stage. Equation
(38) is referred to as the homological equation.

####  2.4.3 Solution of Homological Equation

Recall that $\Psi_{l}$ is a homogeneous polynomial of degree $l$ in the
complex variables
$(x^{\prime}_{R},x^{\prime}_{z},\bar{x}^{\prime}_{R},\bar{x}^{\prime}_{z})$,
i.e., a sum of terms $\propto x_{R}^{\prime k_{R}}x_{z}^{\prime
k_{z}}\bar{x}_{R}^{\prime\bar{k}_{R}}\bar{x}_{z}^{\prime\bar{k}_{z}}$ with
$k_{R}+k_{z}+\bar{k}_{R}+\bar{k}_{z}=l$. Let $\langle\Psi_{l}\rangle$ denote
the collection of terms in $\Psi_{l}$ with $k_{R}=\bar{k}_{R}$ and
$k_{z}=\bar{k}_{z}$ and let

$\\{\Psi_{l}\\}:=\Psi_{l}-\langle\Psi_{l}\rangle=\sum_{i}C_{i}x_{R}^{\prime
k_{R,i}}x_{z}^{\prime
k_{z,i}}\bar{x}_{R}^{\prime\bar{k}_{R,i}}\bar{x}_{z}^{\prime\bar{k}_{z,i}}~{}.$
(39)

Also note that

$[H_{2},\chi_{k}]=-\mathrm{i}\kappa\left(\bar{x}^{\prime}_{R}\frac{\partial\chi_{k}}{\partial\bar{x}^{\prime}_{R}}-{x}^{\prime}_{R}\frac{\partial\chi_{k}}{\partial{x}^{\prime}_{R}}\right)-\mathrm{i}\nu\left(\bar{x}^{\prime}_{z}\frac{\partial\chi_{k}}{\partial\bar{x}^{\prime}_{z}}-{x}^{\prime}_{z}\frac{\partial\chi_{k}}{\partial{x}^{\prime}_{z}}\right)~{}.$
(40)

At each iteration, we solve Equation (38) by taking

$\displaystyle H^{\prime}_{l}$ $\displaystyle=\langle\Psi_{l}\rangle$ (41)
$\displaystyle\chi_{l}$
$\displaystyle=\sum_{i}\frac{\mathrm{i}C_{i}}{(k_{R,i}-\bar{k}_{R,i})\kappa+(k_{z,i}-\bar{k}_{z,i})\nu}x_{R}^{\prime
k_{R,i}}x_{z}^{\prime
k_{z,i}}\bar{x}_{R}^{\prime\bar{k}_{R,i}}\bar{x}_{z}^{\prime\bar{k}_{z,i}}~{}.$
(42)

Equation (42) assumes a non-resonance condition on $\kappa$ and $\nu$: we
cannot have terms333Such resonant terms can instead be incorporated into
$H^{\prime}_{l}$ by simply omitting the corresponding terms from the sum in
Equation (42) defining $\chi_{l}$. The transformed Hamiltonian would, however,
no longer be a function of actions alone. with
$(k_{R,i}-\bar{k}_{R,i})\kappa+(k_{z,i}-\bar{k}_{z,i})\nu=0$.

###  2.5 Azimuthal degree of freedom

Thus far we have only considered the orbits of stars in the reduced phase
space, ignoring their azimuthal degree of freedom. Quasi-periodic orbits in
the full 6 dimensional phase space will be characterized by a third dynamical
frequency, $\Omega_{\phi}$, in addition to the two frequencies,
$\Omega_{R}=\frac{1}{x_{R}}\frac{\partial}{\partial\bar{x}_{R}}H^{\prime}$ and
$\Omega_{z}=\frac{1}{x_{z}}\frac{\partial}{\partial\bar{x}_{z}}H^{\prime}$,
readily obtained from the Birkhoff normalization procedure. Here we detail how
this third frequency can be determined and used to define an additional angle
variable, $\theta_{\phi}$, once the Birkhoff normalization procedure has been
applied to solve for orbits’ dynamical evolution in the reduced phase space.

The time derivative of the azimuthal angle, $\phi$ is given in terms of
complex canonical variables by

$\displaystyle\frac{d}{dt}{\phi}=\frac{L}{R_{C}^{2}}\left(1+\frac{x_{R}+\bar{x}_{R}}{R_{C}\sqrt{2\kappa}}\right)^{-2}~{}.$
(43)

For non-circular orbits, the right-hand side of Equation (43) is not constant
in time but will oscillate about a mean value that defines the azimuthal
frequency, $\Omega_{\phi}$. Introducing the canonical angle variable
$\theta_{\phi}=\Omega_{\phi}t+\phi_{0}$ we can write the $\phi$ coordinate
explicitly as

$\displaystyle\phi=\theta_{\phi}+\rho_{\phi}(x^{\prime}_{R},x^{\prime}_{z},\bar{x}^{\prime}_{R},\bar{x}^{\prime}_{z},L)$
(44)

where $\rho$ is a zero-mean $2\pi$-periodic function of
$\theta_{R}=-\arg{x^{\prime}_{R}}$ and $\theta_{z}=-\arg{x^{\prime}_{z}}$. ,
We can combine Equations (43) and (44) to express $\Omega_{\phi}$ and
$\rho_{\phi}$ in terms of the variables, $x^{\prime}_{R}$ and
$x^{\prime}_{z}$, and their complex conjugates. Recalling the notation
$\langle\cdot\rangle$ and $\\{\cdot\\}$ introduced above in Section 2.4.3 to
denote mean and oscillating terms of a Poisson series, the equations

$\displaystyle\Omega_{\phi}$ $\displaystyle=$
$\displaystyle\frac{L}{R_{C}^{2}}\left\langle\left(\exp[\mathcal{L}_{\chi}]\sum_{m=0}^{\infty}\binom{-2}{m}\left(\frac{x_{R}+\bar{x}_{R}}{R_{C}\sqrt{2\kappa}}\right)^{m}\right)\right\rangle$
(45) $\displaystyle\frac{d}{dt}\rho_{\phi}$ $\displaystyle=$
$\displaystyle\frac{L}{R_{C}^{2}}\left\\{\left(\exp[\mathcal{L}_{\chi}]\sum_{m=0}^{\infty}\binom{-2}{m}\left(\frac{x_{R}+\bar{x}_{R}}{R_{C}\sqrt{2\kappa}}\right)^{m}\right)\right\\}~{}$
(46)

define $\Omega_{\phi}$ and the time derivative of $\rho_{\phi}$ as series in
the variables $x^{\prime}_{R}$ and $x^{\prime}_{z}$ and their complex
conjugates that can be expanded up to some maximum desired order. In order to
derive an expression for $\rho_{\phi}$, we express the series on the right
hand side of Equation (46) as a sum of monomial terms

$\displaystyle\left\\{\left(\exp[\mathcal{L}_{\chi}]\sum_{m=0}^{\infty}\binom{-2}{m}\left(\frac{x_{R}+\bar{x}_{R}}{R_{C}\sqrt{2\kappa}}\right)^{m}\right)\right\\}=\sum_{i}A_{i}x_{R}^{\prime
k_{R,i}}x_{z}^{\prime
k_{z,i}}\bar{x}_{R}^{\prime\bar{k}_{R,i}}\bar{x}_{z}^{\prime\bar{k}_{z,i}}$
(47)

an then integrate with respect to time to obtain

$\displaystyle\rho(x^{\prime}_{R},x^{\prime}_{z},\bar{x}^{\prime}_{R},\bar{x}^{\prime}_{z},L)=\mathrm{i}\frac{L}{R_{C}^{2}}\sum_{i}\frac{A_{i}x_{R}^{\prime
k_{R,i}}x_{z}^{\prime
k_{z,i}}\bar{x}_{R}^{\prime\bar{k}_{R,i}}\bar{x}_{z}^{\prime\bar{k}_{z,i}}}{(k_{R,i}-\bar{k}_{R,i})\kappa+(k_{z,i}-\bar{k}_{z,i})\nu}~{}.$
(48)

It follows immediately from Hamilton’s equations that $J_{\phi}=L$ serves as
canonical action coordinate conjugate to the angle $\theta_{\phi}$. We
therefore have established a complete canonical transformation from the
original variables $(R,z,\phi,p_{R},p_{z},L)$ of the 6 dimnensional phase
space into to a complete set AA variables,
$(\theta_{R},\theta_{z},\theta_{\phi},J_{R},J_{z},J_{\phi})$ where
$J_{i}=|x^{\prime}_{i}|^{2}$ and $\theta_{i}=-\arg x^{\prime}_{i}$ for the
subscript $i\in\\{z,R\\}$.

###  2.6 Numerical Implementation

The Birkhoff normalization algorithm described above is implemented in the
poisson_series module of the open-source Python package celmech (Hadden &
Tamayo, 2022). This module provides routines for representing and manipulating
so-called “Poisson series” (Brumberg, 1995). These are series in $N$ complex
canonical variables $x_{i}$ and their complex conjugates, $\bar{x}_{i}$ along
with $M$ real action-angle pairs, $(P_{i},Q_{i})$, that are given as sums of
individual monomial terms of the form

${\cal
M}(k,\bar{k},p,q)=\prod_{i=1}^{N}x_{i}^{k_{i}}\bar{x}_{i}^{\bar{k}_{i}}\prod_{j=1}^{M}P_{j}^{p_{j}}\exp[iq_{j}Q_{j}]~{}.$
(49)

The Poisson bracket of any pair of Poisson series is another Poisson series,
as can be readily seen by considering the Poisson bracket of two monomial
terms. Defining the length-$N$ vectors $o^{N}_{i}$ such that
$[o^{N}_{i}]_{j}=\delta_{ij}$ where $\delta_{ij}$ is the Kroenecker delta, the
Poisson bracket of two monomial terms is then a Poisson series consisting of
(at most) $N+M$ new monomial terms and given by

$[{\cal M}(k,\bar{k},p,q),{\cal
M}(l,\bar{l},r,s)]=\mathrm{i}\sum_{i=1}^{N}(\bar{k}_{i}{l}_{i}-k_{i}\bar{l}_{i}){\cal
M}(k+l-o^{N}_{i},\bar{k}+\bar{l}-o^{N}_{i},p+r,q+s)+\\\
\mathrm{i}\sum_{j=1}^{M}(q_{j}r_{j}-s_{j}p_{j}){\cal
M}(k+l,\bar{k}+\bar{l},p+r-o^{M}_{j},q+s)~{}.$ (50)

The poisson_series module defines the PoissonSeries class to store sums of
monomial terms of the form given in Equation (49) multiplied by numerical
coefficients. It also defines a number of routines for manipulating these
series, including performing scalar multiplication, multiplying or adding
pairs of series, and calculating Poisson brackets. Note that the Poisson
series used in this paper involve only complex canonical variables (i.e.,
$M=0$).

The poisson_series module additionally defines routines for carrying out the
Birkhoff normalization procedure and evaluating finite-order approximations of
exponential operators, $\exp[\mathcal{L}_{\chi}]$, applied to series. These
routines require grouping terms in the expansions of various functions, such
as the Hamiltonian, $H$, and the generating function, $\chi$, by their order
in powers of the complex canonical variables. In the framework of the
poisson_series module, this grouping is done by using Python dictionary
objects to associate PoissonSeries objects representing terms in an expansion
of a given order to key values that indicate their order. Example Jupyter
notebooks illustrating the use of celmech’s poisson_series module to carry out
the Birkhoff normalization procedure are available at
github.com/shadden/AA_variables_via_Birkhoff_Normalization.

###  2.7 Padé Approximants

Since the Birkhoff normalization procedure is based upon a Taylor expansion of
the effective potential, the method can only be expected to provide accurate
results for orbits that remain within regions of phase space where this
expansion is convergent. In practice, this can make the procedure inapplicable
to orbits that experience significant excursions in the vertical direction.
For example, the Miyamoto-Nagai potential,

$\Phi(R,z)=-\frac{1}{\sqrt{R^{2}+(a+\sqrt{z^{2}+b^{2}})^{2}}}~{}$ (51)

with parameters $a$ and $b$, is a commonly adopted as a model for the galactic
disk’s contribution to the gravitational potential (Miyamoto & Nagai, 1975).
The parameter $b$ effectively sets the thickness of the galactic disk and the
branch points of Equation (51) occurring in the complex plane at $z=\pm ib$
restrict the convergence Taylor series expansions to $|z|<b$.

Accurate local approximations of functions can sometimes be extended beyond
the domain of convergence of their Taylor series with the use of Padé
approximants. The $(n,m)$ Padé approximant of a function is the rational
function with numerator of degree $m$ and denominator of degree $n$ that
agrees with the Taylor series expansion of the function up to degree $x^{m+n}$
(e.g., Teukolsky et al., 1992). If the first $n+m$ terms of a function’s
Taylor series are $\sum_{k=0}^{n+m}c_{k}x^{k}$ then the coefficients of an
$(n,m)$ Padé approximant are readily obtained by multiplying both left and
right hand sides of

$\displaystyle\sum_{k=0}^{n+m}c_{k}x^{k}=\frac{a_{0}+a_{1}x+...a_{m}x^{m}}{1+b_{1}x+...+b_{n}x^{n}}\mod{x^{n+m+1}}$
(52)

by $1+b_{1}x+...+b_{n}x^{n}$. The coefficients, $b_{i}$, satisfy the linear
equation

$\displaystyle\begin{pmatrix}c_{n+m-1}&\ldots&c_{m}\\\ \vdots&\ddots&\vdots\\\
c_{m}&\ldots&c_{m-n+1}\end{pmatrix}\cdot\begin{pmatrix}b_{1}\\\ \vdots\\\
b_{n}\end{pmatrix}=-\begin{pmatrix}c_{n+m}\\\ \vdots\\\ c_{m+1}\end{pmatrix}.$
(53)

and the coefficients, $a_{i}$, are then given by
$a_{i}=\sum_{k=0}^{n}c_{i-k}b_{k}$ where $b_{0}:=1$ and $c_{k}$ is taken to be
$0$ for $k<0$.

Below in Section 3, we show that constructing Padé approximants in the
variable $I_{z}=|x_{z}|^{2}$, often improves the accuracy of the Birkhoff
normalization procedure. The Padé approximants involve coefficients, $b_{i}$
and $a_{i}$, that are functions of $x_{R},~{}\bar{x}_{R}$ and
$\phi_{z}=-\arg(x_{z})$. These coefficients are determined by evaluating
finite-order truncations of the series
$x^{\prime}_{R}=\exp[-\mathcal{L}_{\chi}]x_{R}$ and
$x^{\prime}_{z}=\exp[-\mathcal{L}_{\chi}]x_{z}$ and then grouping terms by
powers of $I_{z}$. Note that the generating function, $\chi$, is developed as
a series in powers of $x_{z}$ and $\bar{x}_{z}$, which are both proportional
to $\sqrt{I_{z}}$. However, it is straightforward to show that for any
$\Phi(R,z)$ like Equation (51) that is an even function of $z$, only integer
powers of $I_{z}$ will appear in $\chi$ as well as finite-order truncations of
$\exp[-\mathcal{L}_{\chi}]x_{R}$. Thus, these finite-order truncations are
polynomials in $I_{z}$ (rather than $\sqrt{I_{z}}$) and can be used to obtain
Padé approximants in $I_{z}$. Additionally, finite-order truncations of
$\exp[-\mathcal{L}_{\chi}]x_{z}$ yield expressions that can be written as
$x_{z}P(I_{z})$ where $P$ is a polynomial in $I_{z}$ with coefficients that
depend on $x_{R},\bar{x}_{R},$ and $\phi_{z}$. Padé approximants are derived
from the polynomial $P$ and used to calculate $\exp[-\mathcal{L}_{\chi}]x_{z}$
in Section 3.

##  3 Tests

In this section we apply the Birkhoff normalization procedure to orbits in the
Miyamoto-Nagai potential in Equation (51). The parameters $a=3$, $b=0.3$, and
$L=3$ are chosen, with a corresponding circular orbit radius of
$R_{C}(L)\approx 10.4$. With distance units taken to be kpc, these parameters
approximately match the parameters adopted for the Miyamoto-Nagai component of
the Milky Way potential model, MWPotential2014, implemented in the galpy
Python package (Bovy, 2015). The Birkhoff normalization procedure is carried
out by developing the generating function, $\chi$, up to 10th order in complex
canonical variables. The sympy symbolic mathematics package (Meurer et al.,
2017) is used to calculate the requisite partial derivatives of the effective
potential.

We also compute actions and frequencies using the Stäckel approximation method
introduced by Binney (2012). These calulations are done using the
actionAngle.actionAngleStäckel function implemented in galpy (Bovy & Rix,
2013; Bovy, 2015). This function requires as input an estimated focal length
parameter, $\Delta$, in addition to a potential and orbital initial
conditions. For all orbits considered, The galpy function
actionAngle.estimateDeltaStäckel is used to estimate the focal length based on
partial derivatives of the potential evaluated at $(R,z)=(R_{C},0)$ following
Sanders (2012).

Figure 1: Color maps showing r.m.s. fractional variations in estimated values
of actions on a grid of orbits in a Miyamoto-Nagai potential (Equation 51)
with parameters $a=3$ and $b=0.3$. Actions are computed along orbits starting
with $(L,R,z)=(3,R_{C},0)$ and initial velocity values indicated by the axes,
which are plotted in units of the circular velocity, $v_{C}=L/R_{C}$. The top
row shows variations in the radial action, $J_{R}$, while the bottom row shows
variations in the vertical action, $J_{z}$. Columns are labeled according to
the method used to compute action values. See text for additional details.
Dash-dot, dashed, and solid contours mark fractional error levels of
$0.1\%,1\%$ and $10\%$, respectively.

Figure 1 illustrates the accuracy of actions determined by Birkhoff
normalization and Stäckel approximation across a range of initial conditions.
Orbits are initialized at $(R,z)=(R_{C},0)$ over the plotted range of
$(p_{R},p_{z})$ values and integrated for 10 times the radial epicyclic period
(i.e., $T=10\times(2\pi/\kappa)$) with 512 outputs generated at equally-spaced
intervals along the orbits. Transformed complex canonical variable,
$x^{\prime}_{i}$, are computed as functions of the un-transformed variables
along a given orbit using various approximations to the series

$\displaystyle
x^{\prime}_{i}(x_{R},x_{z},\bar{x}_{R},\bar{x}_{z})=(\exp[-\mathcal{L}_{\chi}]x_{i})(x_{R},x_{z},\bar{x}_{R},\bar{x}_{z})$
(54)

where the subscript $i$ denotes $R$ or $z$. Action values are then computed
from the transformed complex canonical variable as
$J_{i}=|x^{\prime}_{i}|^{2}$. The root-mean-square (r.m.s.) fractional
variations in these action values computed along each orbit are recorded in
the color scale.

The leftmost column of Figure 1 show the r.m.s. variations in actions computed
using the epicyclic approximation, where Equation (54) is approximated simply
as $x_{i}^{\prime}\approx x_{i}$. The variations in actions computed via this
approximation provide a measure of how significantly orbits are affected by
the anharmonicity of the potential. The second column of Figure 1 shows the
r.m.s. variations for actions computed by expanding the right hand side of
Equation (54) as a Taylor series. With $\chi$ determined up to 10th order in
the canonical variables, the terms of these series are determined up to 9th
order in the complex canonical variables. The third and fourth columns show
r.m.s. variations for actions computed by constructing $(3,1)$ and $(2,2)$
Padé approximants in $I_{z}$ as described in Section 2.7. The final columns
the r.m.s. variations for actions computed using the Stäckel approximation.

The Taylor series approximations of actions in Figure 1 provide good agreement
with the numerical results for mid-plane vertical velocity values
$p_{z}\lesssim 0.1v_{C}$, where $v_{C}=L/R_{C}$ is the circular orbital
velocity. The poor agreement for greater vertical velocities is attributable
to the fact that orbits with greater vertical velocities reach heights $z>b$,
beyond the radius of convergence of the Taylor series expansion of the
potential in Equation (51). An estimate of the maximum vertical height,
$z_{\mathrm{max}}$, achieved by an orbit with a given vertical velocity,
$p_{z,0}$, in the midplane can be obtained by neglecting any motion in the
radial direction and using conservation of energy to equate
$\frac{1}{2}p_{z,0}\approx\Phi_{\mathrm{eff}}(R_{C},z_{\mathrm{max}})-\Phi_{\mathrm{eff}}(R_{C},0)$.
Taking $z_{\mathrm{max}}=b$, the corresponding midplane velocity is
$p_{z,0}\approx 0.09v_{C}$, closely matching the value of mid-plane velocity
above which the Taylor series approximations fail in Figure 1. The Padé
approximants clearly extend the range over which the Birkhoff normalization
procedure accurately predicts actions, achieving r.m.s. variations in radial
actions of $\lesssim 10\%$ for initial vertical velocities $p_{z}\lesssim
0.2v_{C}$. The (2,2) Padé approximant is slightly more accurate than the (3,1)
approximant. Performance of the algorithm for vertical action determination is
slightly worse, with $\lesssim 10\%$ accuracy restricted to initial vertical
velocities $p_{z}\lesssim 0.15v_{C}$.

Figure 2 compares the accuracy of action determinations made using Birkhoff
normalization with (2,2) Padé approximants against the Stäckel approximation.
Birkhoff normalization provides a more accurate determination of radial
actions over the majority of the plotted range. Vertical actions are also more
accurately determined for $p_{z}\lesssim 0.1v_{C}$.

Figure 2: Comparison of the accuracy of actions determined using Birkhoff
normalization and (2,2) Padé approximants versus the Stäckel approximation.
The color scale indicates the r.m.s. variation in actions computed using the
Stäckel approximation divided by r.m.s. variations computed via Birkhoff
normalization so that red colors indicate orbits where Birkhoff normalization
yields better accuracy while blue colors indicate orbits for which the Stäckel
approximation is more accurate. The grid of orbits is the same as shown in
Figure 1.

Figure 3 illustrates the accuracy of the determination of the dynamical
frequencies $\Omega_{R}$, $\Omega_{z}$, and $\Omega_{\phi}$ using various
methods. Exact values of frequencies are determined numerically using the
frequency modified Fourier transform algorithm of Šidlichovský & Nesvorný
(1996), implemented in the celmech package as the function
miscellaneous.frequency_modified_fourier_transform. The method was applied to
the time series of complex canonical variables, $x_{R}$ and $x_{z}$, for each
orbit integration and the frequencies of the largest Fourier components
determined for each variable were identified as $-\Omega_{R}$ and
$-\Omega_{z}$, respectively. Values of $\Omega_{\phi}$ were determined by
applying the same method to time series of $e^{\mathrm{i}\phi}$. The first
column of Figure 3 records the fractional errors $\kappa/\Omega_{R}-1$,
$\nu/\Omega_{z}-1$, and $R_{C}^{-2}L/\Omega_{\phi}-1$ incurred by adopting the
epicyclic approximation for the dynamical frequencies. These panels illustrate
that the frequencies $\Omega_{R}$ and $\Omega_{\phi}$ vary gradually with both
initial radial and vertical velocities while the frequency $\Omega_{z}$ is a
strong function of the initial vertical velocity.

The second and third columns of Figure 3 record fractional errors incurred by
using Birkhoff normalization and Padé approximants to determine frequencies.
Frequencies were determined as follows. First, values of complex canonincal
variables were assigned based on orbits’ initial conditions as
$x_{R}=\mathrm{i}\sqrt{\frac{1}{2\kappa}}p_{R}$ and
$x_{z}=\mathrm{i}\sqrt{\frac{1}{2\nu}}p_{z}$ (see Equations 5 and 6). Next,
values of transformed complex canonical variables, $x^{\prime}_{R}$ and
$x^{\prime}_{z}$, were computed using the Padé approximation method described
in Section 2.7. Then, series expressions for frequencies,
$\Omega_{R}(J_{R},J_{z})=\partial_{J_{R}}H^{\prime}$ and
$\Omega_{z}(J_{R},J_{z})=\partial_{J_{z}}H^{\prime}$, derived from the
transformed Hamiltonian, $H^{\prime}$, and the series expression for
$\Omega_{\phi}$ given by Equation (45) were computed. Finally, these series
expressions were used to construct a second set of (3,1) and (2,2) Padé
approximants, this time in $J_{z}=|x^{\prime}_{z}|^{2}$, and evaluated.
Comparing Figures 1 and 3, we see that dynamical frequencies are accurately
recovered in essentially the same regions of phase space that Birkhoff
normalization and Padé approximants provide accurate action variables. The
fourth column of Figure 3 shows frequency errors incurred using the Stäckel
approximation for comparison.

Figure 3:  Color maps showing fractional errors in the estimated values of
radial, vertical, and azimuthal frequencies on a grid of orbits in the
Miyamoto-Nagai potential (Equation (51)) with parameters $a=3$ and $b=0.3$.
Frequencies are computed for orbits starting with $(L,R,z)=(3,R_{C},0)$ and
initial velocity values indicated by the axes, which are plotted in units of
the circular velocity, $v_{C}=L/R_{C}$. The top, middle, and bottom rows shows
fractional errors in the radial ($\Omega_{R}$), vertical ($\Omega_{z}$), and
azimuthal ($\Omega_{\phi}$) frequencies, respectively. Columns are labeled
according to the method used to compute the frequency values. See text for
details. Dash-dot, dashed, and solid contours mark fractional error levels of
$0.1\%,1\%$ and $10\%$, respectively.

##  4 Summary & Discussion

We have shown how Birkhoff normalization, a technique from Hamiltonian
perturbation theory, can be algorithmically implemented and applied to derive
AA variables for stellar orbits in axisymmetric galactic potentials. One
significant advantage of this method is that it provides explicit expressions
for actions and angle variables as functions of positions and velocities and
vice versa. It similarly yields an explicit expression for the Hamiltonian as
a function of the action variables. The only input required by the method is
the partial derivatives of the effective potential up to the maximum desired
order of the series expansions.

Since Birkhoff normalization relies on expansion of the effective potential
about its value for circular, planar orbits, its range of applicability is
nominally limited by the radius of convergence of this expansion. However, as
illustrated by the tests in Section 3, the series expressions generated by the
procedure can be used to construct Padé approximants and extend the range over
which it provides accurate results. The tests presented in Section 3 showed
that Birkhoff normalization provides comparable or better accuracy than the
popular Stäckel approximation, at least for quasi-planar orbits.

The expressions provided by Birkhoff normalization take the form of
polynomials involving the complex canonical variables introduced in Section
2.2. This makes evaluating both forward and inverse transformations between
position-velocity data in AA variables exceptionally computationally
efficient. This computational efficiency represents another potential
advantage of the method over the Stäckel approximation, which requires
evaluation of integrals via numerical quadrature. The most computationally
expensive parts of the Birkhoff normalization algorithm will generally be the
initial calculation of potential derivatives and construction of series.
However, these calculations represent a fixed cost and need only be performed
once, after which the Birkhoff normalization procedure can be applied to as
many orbits as desired. For the tests presented in Section 3, computing
partial derivatives of the Miyamoto-Nagai potential up to 10th order via
symbolic means required $\sim 20~{}\mathrm{s}$ on a 2.5 GHz Intel i7
processor. The initial construction of the series expression for $\chi$ up to
10th order in complex canonical variables took $\sim 3~{}\mathrm{s}$ and
calculations of additional series expressions approximating the exponential
operators $\exp[\pm\mathcal{L}_{\chi}]$ applied to various functions, such as
Equation (54), took approximately $\sim 1~{}\mathrm{s}$ each. A potential
downside to the Birkhoff normalization method, in terms of computational cost,
is that the effective potential partial derivatives and series must nominally
be recalculated for each value of angular momentum considered. A simple
strategy for applying the method to orbits with a range of angular momenta
would be to carry out series constructions on a 1 dimensional grid of angular
momenta and simply interpolate results. Alternatively, the construction of
various series could be carried out symbolically with the aid of a computer
algebra system. We leave the development of such methods to future work as the
choice of optimal strategy is likely to be application-specific.

We conclude by mentioning a few potential future applications of the Birkhoff
normalization procedure presented in this paper. First, analytic expressions
furnished by the Birkhoff normalization procedure could make it well-suited
for problems requiring the representation of equilibrium distribution
functions, which by Jean’s theorem must depend only on integrals of motion, in
terms of spatial positions and velocities. Additionally, the explicit
expression for the Hamiltonian in terms of action variables provided by the
Birkhoff normalization procedure can be used as a starting point for using
perturbation theory to understand the non-equilibrium signatures of dynamical
perturbations like those induced by a stellar bar (e.g., Binney, 2018) or
satellite galaxy (e.g., Banik et al., 2022). Finally, even though the Birkhoff
normalization procedure may not supply accurate AA variables for orbits with
large vertical actions, it could provide useful initial ‘toy’ AA variables
that serve as the starting point of the torus mapping method of McGill &
Binney (1990).

##  5 Acknowledgments

I thank Rimpei Chiba and Neige Frankel for helpful discussions. I thank Scott
Tremaine for suggesting the use of Padé approximants. Simulations were
performed on the Sunnyvale cluster at the Canadian Institute for Theoretical
Astrophysics. I acknowledge support by the Natural Sciences and Engineering
Research Council of Canada (NSERC), funding references CITA 490888-16 and
RGPIN-2020-03885. This project was developed in part at the Gaia Hike, a
workshop hosted by the University of British Columbia and the Canadian
Institute for Theoretical Astrophysics in 2022 June.

## References

  * Banik et al. (2022) Banik, U., Weinberg, M. D., & van den Bosch, F. C. 2022, ApJ, 935, 135, doi: 10.3847/1538-4357/ac7ff9
  * Binney (2012) Binney, J. 2012, MNRAS, 426, 1324, doi: 10.1111/j.1365-2966.2012.21757.x
  * Binney (2018) —. 2018, MNRAS, 474, 2706, doi: 10.1093/mnras/stx2835
  * Binney & Spergel (1982) Binney, J., & Spergel, D. 1982, ApJ, 252, 308, doi: 10.1086/159559
  * Binney & Spergel (1984) —. 1984, MNRAS, 206, 159, doi: 10.1093/mnras/206.1.159
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton University Press)
  * Birkhoff (1927) Birkhoff, G. D. 1927, Dynamical systems, Vol. 9 (American Mathematical Soc.)
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29, doi: 10.1088/0067-0049/216/2/29
  * Bovy & Rix (2013) Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115, doi: 10.1088/0004-637X/779/2/115
  * Brumberg (1995) Brumberg, V. A. 1995, Analytical Techniques of Celestial Mechanics (Berlin: Springer-Verlag)
  * Contopoulos (1963) Contopoulos, G. 1963, AJ, 68, 1, doi: 10.1086/108903
  * Contopoulos (2003) —. 2003, Journal of Physics A: Mathematical and General, 36, 8639.  http://resolver.scholarsportal.info/resolve/03054470/v36i0032/8639_nofiom.xml
  * de Zeeuw & Merritt (1983) de Zeeuw, T., & Merritt, D. 1983, ApJ, 267, 571, doi: 10.1086/160894
  * Deprit (1969) Deprit, A. 1969, Celestial Mechanics, 1, 12, doi: 10.1007/BF01230629
  * Efthymiopoulos et al. (2004) Efthymiopoulos, C., Giorgilli, A., & Contopoulos, G. 2004, Journal of Physics A Mathematical General, 37, 10831, doi: 10.1088/0305-4470/37/45/008
  * Ferraz-Mello (2007) Ferraz-Mello, S. 2007, Canonical perturbation theories: degenerate systems and resonance, Vol. 345 (Springer Science & Business Media)
  * Goldstein et al. (2002) Goldstein, H., Poole, C., & Safko, J. 2002, Classical mechanics (Addison-Wesley)
  * Gustavson (1966) Gustavson, F. G. 1966, AJ, 71, 670, doi: 10.1086/110172
  * Hadden & Tamayo (2022) Hadden, S., & Tamayo, D. 2022, AJ, 164, 179, doi: 10.3847/1538-3881/ac8d01
  * Henon & Heiles (1964) Henon, M., & Heiles, C. 1964, AJ, 69, 73, doi: 10.1086/109234
  * Lichtenberg & Lieberman (1992) Lichtenberg, A., & Lieberman, M. 1992, Regular and Chaotic Dynamics (New York: Springer-Verlag)
  * Lowenstein (2012) Lowenstein, J. H. 2012, Essentials of Hamiltonian dynamics (Cambridge, UK: Cambridge University Press)
  * McGill & Binney (1990) McGill, C., & Binney, J. 1990, MNRAS, 244, 634.  https://ui.adsabs.harvard.edu/abs/1990MNRAS.244..634M
  * Meurer et al. (2017) Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Computer Science, 3, e103, doi: 10.7717/peerj-cs.103
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
  * Saaf (1968) Saaf, A. F. 1968, ApJ, 154, 483, doi: 10.1086/149776
  * Sanders (2012) Sanders, J. 2012, MNRAS, 426, 128, doi: 10.1111/j.1365-2966.2012.21698.x
  * Sanders & Binney (2016) Sanders, J. L., & Binney, J. 2016, MNRAS, 457, 2107, doi: 10.1093/mnras/stw106
  * Teukolsky et al. (1992) Teukolsky, S. A., Flannery, B. P., Press, W., & Vetterling, W. 1992, SMR, 693, 59
  * Šidlichovský & Nesvorný (1996) Šidlichovský, M., & Nesvorný, D. 1996, Celestial Mechanics and Dynamical Astronomy, 65, 137, doi: 10.1007/BF00048443