# Dynamics of McMillan mappings I.
McMillan multipoles

T. Zolkin<EMAIL_ADDRESS>Fermilab, PO Box 500, Batavia, IL 60510-5011  S.
Nagaitsev<EMAIL_ADDRESS>Brookhaven National Laboratory, Upton, NY 11973
Old Dominion University, Norfolk, VA 23529  I. Morozov<EMAIL_ADDRESS>Novosibirsk State Technical University, Novosibirsk 630073, Russia

###### Abstract

In this article, we consider two dynamical systems: the McMillan sextupole and
octupole integrable mappings, originally proposed by Edwin McMillan. Both
represent the simplest symmetric McMillan maps, characterized by a single
intrinsic parameter. While these systems find numerous applications across
various domains of mathematics and physics, some of their dynamical properties
remain unexplored. We aim to bridge this gap by providing a comprehensive
description of all stable trajectories, including the parametrization of
invariant curves, Poincaré rotation numbers, and canonical action-angle
variables.

In the second part, we establish connections between these maps and general
chaotic maps in standard form. Our investigation reveals that the McMillan
sextupole and octupole serve as first-order approximations of the dynamics
around the fixed point, akin to the linear map and quadratic invariant (known
as the Courant-Snyder invariant in accelerator physics), which represents
zeroth-order approximations (referred to as linearization). Furthermore, we
propose a novel formalism for nonlinear Twiss parameters, which accounts for
the dependence of rotation number on amplitude. This stands in contrast to
conventional betatron phase advance used in accelerator physics, which remains
independent of amplitude. Notably, in the context of accelerator physics, this
new formalism demonstrates its capability in predicting dynamical aperture
around low-order resonances for flat beams, a critical aspect in beam
injection/extraction scenarios.

##  I Introduction

A generic dynamical system typically endures the presence of chaos. The motion
is chaotic when there is an insufficient amount of integrals of motion, also
known as invariants or conserved quantities. Such systems are challenging to
analyze and often considered intractable due to their high sensitivity to
initial conditions, rounding errors, and external perturbations, rendering
practical solutions computationally demanding.

Conversely, when a system possesses an adequate number of first integrals, it
is termed integrable, featuring relatively simple and predictable dynamics.
Historically, completely integrable problems in physics have played a crucial
role in understanding more complex scenarios. For instance: (i) A multi-
dimensional harmonic oscillator serves as a good approximation for small
amplitude oscillations in general nonlinear systems (linearization). (ii) The
Kepler problem provides insight into celestial mechanics, such as Earth’s
orbit, despite the chaotic nature of the entire solar system. (iii) Planetary
motion around two fixed centers has applications in satellite trajectory
calculations Alexeev (1965); Marchal (1966, 1986), electron acceleration in
atomic collisions Jakas (1995, 1996), and determining energy levels in
${\mathrm{H}_{2}}^{+}$ Strand and Reinhardta) (1979).

In this article, we focus on another renowned example known as the McMillan
integrable map McMillan (1971). Originally proposed by Edwin McMillan as a
simplified accelerator lattice with one degree of freedom, consisting of
linear optics (corresponding to simple linear transformation in phase space)
and a special thin nonlinear lens (nonlinear vertical shear transformation).
As later demonstrated by Suris Suris (1989), the McMillan mapping is among the
few possible integrable symplectic transformations of the plane with an
analytic integral of motion. Specifically, for maps in the form (1), the
invariant must be a regular, exponential, or trigonometric polynomial of
degree two in coordinate and momentum.

We demonstrate that the McMillan map is not merely a model lattice but is also
related to general symplectic mappings of the same form, featuring chaotic
dynamics as its first and second-order approximations (Section VI). This
insight aids in “integrating-out” additional nonlinear features of small
amplitude oscillations, such as detuning or the frequency’s dependence on the
action variable. As a consequence, we propose the concept of nonlinear Twiss
parameters for accelerator physics, serving as a natural extension to the
existing Courant-Snyder formalism. This new framework enables a deeper
understanding of amplitude-dependent shifts in the betatron tune and
facilitates predicting the dynamical aperture around isolated low order
resonances, thereby representing a nonlinear integrable model for slow
resonant beam extraction for a flat beam.

Finally, in first part of the article, we fill the gap in analytical results
by providing a canonical transformation to action-angle variables (Section
III) and complimentary bifurcation diagrams for McMillan sextupole and
octupole mappings (Section B); these are the simplest symmetric McMillan maps
with only one intrinsic parameter and quadratic (sextupole) or cubic
(octupole) nonlinearities. Octupole, and even more general asymmetric McMillan
maps, have been studied in great detail, especially by Iatrou and Roberts
Iatrou and Roberts (2002) who provided parametrization of individual curves.
However, explicit expressions for rotation number as a function of canonical
action integral have been missing until now. This knowledge is critical in
understanding underlying dynamics and design of integrable accelerator
lattices with desirable properties. Additionally, these analytical results can
be applied to more realistic integrable accelerator lattices based on McMillan
map, e.g. they describe a particle with zero angular momentum in the
integrable (in 4D phase space) accelerator ring with axially symmetric
McMillan electron lens Cathey  _et al._ (2021).

The article is structured as follows: We commence by establishing definitions
for McMillan multipoles and drawing parallels with auxiliary classical
Hamiltonians in Section II. Subsequently, Section III outlines the general
method for dynamic analysis based on the Danilov theorem, using the McMillan
sextupole as an example. Additionally, it establishes a connection between the
parametrization of invariant curves Iatrou and Roberts (2002) and the set of
canonical action-angle variables for the McMillan octupole. Sections IV and V
delve into a detailed analysis of dynamics by examining amplitude-frequency
dependencies and stability diagrams for different motion regimes. Next,
Section VI introduces low orders of perturbation theory that relates McMillan
multipoles to general chaotic systems of the same map form. We investigate
small and large amplitudes using the Hénon quadratic and cubic mappings as
examples for study. Additionally, we present a general result for a model
accelerator lattice with a thin nonlinear lens and compare our findings with
other techniques such as Lie algebra. Finally, Section VII discusses possible
generalizations and extensions of our results, including the utilization of a
more general form of the map. Here, we apply perturbation theory to the
horizontal dynamics within the Fermilab delivery ring, which is employed for a
third-integer resonant extraction in the Mu2e experiment.

At the conclusion of the article, we provide supplementary materials.
Appendices A and B offer a detailed description of the symmetries of motion
invariants and the stability of fixed points and $n$-cycles for both, McMillan
sextupole and octupole mappings. The last Appendix, C, contains a list of
analytical expressions for action integrals and their power series.

##  II McMillan multipoles

The most general form of the symmetric McMillan map is given by McMillan
(1971); Quispel  _et al._ (1988, 1989); Iatrou and Roberts (2002)

$\begin{array}[]{l}q^{\prime}=p,\\\\[7.11317pt]
p^{\prime}=-q+f(p),\end{array}$ (1)

where the prime (′) denotes the application of the map, and $f(q)$ is the
force function defined as

$f(p)=-\,\frac{\mathrm{B}\,p^{2}+2\,\epsilon\,p+\Xi}{\mathrm{A}\,p^{2}+\mathrm{B}\,p+\Gamma}.$

This map is integrable for any set of parameters, with the invariant of motion
$\mathcal{K}[p,q]$ being a biquadratic function of the coordinate $q$ and
momenta $p$, that can be written in matrix form as:

$\mathcal{K}[p,q]=\begin{bmatrix}p^{2}\\\ p\\\
1\end{bmatrix}^{\mathrm{T}}\left(\mathrm{M}\cdot\begin{bmatrix}q^{2}\\\ q\\\
1\end{bmatrix}\right),\quad\mathrm{M}=\begin{bmatrix}\mathrm{A}&\mathrm{B}&\Gamma\\\
\mathrm{B}&2\,\epsilon&\Xi\\\ \Gamma&\Xi&\Lambda\end{bmatrix}.$ (2)

In order to gain a deeper understanding of the physical meaning behind the
mapping’s parameters, we introduce the auxiliary Hamiltonian of a classical
particle:

$\mathcal{H}[p,q;t]=\frac{p^{2}}{2\,m}+\lambda+\xi\,q+\gamma\,\frac{q^{2}}{2}+\beta\,\frac{q^{3}}{3}+\alpha\,\frac{q^{4}}{4}.$
(3)

Rewriting Eq. (2) explicitly

$\begin{array}[]{l}\mathcal{K}[p,q]=\underbrace{\mathrm{A}\,p^{2}q^{2}}_{\text{cubic
nonlin.}}+\,\underbrace{\mathrm{B}\,(p^{2}q+p\,q^{2})}_{\text{quadratic
nonlin.}}\,+\\\\[19.91684pt]
\qquad\qquad\qquad+\,\underbrace{\Gamma\,(p^{2}+q^{2})+2\,\epsilon\,p\,q}_{\text{harmonic
oscillator}}\,\,+\\!\underbrace{\Xi\,(p+q)}_{\text{defines
origin}}+\,\,\Lambda\end{array}$

allows us to establish a meaningful analogy between the terms. Foremost,
adding a constant term to the Hamiltonian ($\lambda$) or to the invariant of
the map ($\Lambda$) only results in a shift of the energy scale. This shift
doesn’t affect the relative positions or momenta of the particles in the
system, and therefore, it doesn’t alter the physical trajectories or dynamics
of the system. Next, assuming an oscillatory regime (i.e., where $\mathcal{H}$
and $\mathcal{K}$ have a local extremum), we can set $\xi$ and $\Xi$ to zero
by positioning the origin at the stationary point of the Hamiltonian or the
fixed point of the map. Further, when examining quadratic components:

$\frac{p^{2}}{2\,m}+\gamma\,\frac{q^{2}}{2}\qquad\text{and}\qquad\Gamma\,(p^{2}+q^{2})+2\,\epsilon\,p\,q$

it is apparent that both contribute to inducing harmonic linear oscillations
in the system. However, while we can normalize units of time and coordinate to
set $m=\gamma=1$ in the Hamiltonian system, the mapping introduces an
intrinsic irreducible parameter $a=-2\,\epsilon/\Gamma$ associated with the
discrete nature of time (see Appendix A). Finally, the pair of parameters
$\mathrm{A}$ and $\mathrm{B}$ or $\alpha$ and $\beta$ contributes to higher-
order nonlinear effects, playing the roles of cubic and quadratic
nonlinearities respectively.

If only one of nonlinear terms is present, after performing non-
dimensionalization of Hamiltonian (3), we obtain 4 possible situations with
stable trajectories:

$\begin{array}[]{lc}\displaystyle\mathcal{H}_{\mathrm{sxt\,}}[p,q;t]=\frac{p^{2}}{2}+\frac{q^{2}}{2}+\frac{q^{3}}{3},&\mathrm{(SX)}\\\\[7.11317pt]
\displaystyle\mathcal{H}^{-}_{\mathrm{oct\,}}[p,q;t]=\frac{p^{2}}{2}+\frac{q^{2}}{2}-\frac{q^{4}}{4},&\mathrm{(DO)}\\\\[7.11317pt]
\displaystyle\mathcal{H}^{+}_{\mathrm{oct\,}}[p,q;t]=\frac{p^{2}}{2}+\frac{q^{2}}{2}+\frac{q^{4}}{4},&\mathrm{(FO)}\\\\[7.11317pt]
\displaystyle\mathcal{H}_{\mathrm{Duf}}[p,q;t]=\frac{p^{2}}{2}-\frac{q^{2}}{2}+\frac{q^{4}}{4}.&\mathrm{(Df)}\\\\[7.11317pt]
\end{array}\vspace{-0.25cm}$ (4)

These Hamiltonians serve as universal models for nonlinear oscillators,
finding utility across various applications. For instance, the first system
serves as a prototype for the 1D Hénon-Heiles potential Hénon and Heiles
(1964) or represents horizontal dynamics inside long sextupole magnet combined
with linear focusing (SX). The subsequent two models describe planar
defocusing (DO) and focusing octupoles (FO), while the final one is associated
with the well-known unforced undamped Duffing oscillator (Df) Duffing (1918).
Throughout this article, we will use lowercase subscripts “sxt” and “oct” to
describe the leading nonlinearity in the Hamiltonian or invariant, while
uppercase labels (SX), (DO), (FO), and (Df) will be used to specify the
underlying regime of motion.

We further expand our analogy by introducing two dynamical systems: the
McMillan sextupole and octupole integrable mappings, defined by corresponding
forces

$f_{\text{sxt}}(p)=-p\,\frac{p+2\,\epsilon}{p+\Gamma},\qquad\text{and}\qquad
f_{\text{oct}}(p)=-\frac{2\,\epsilon\,p}{p^{2}+\Gamma},$

along with matrices of coefficients in the form (2)

$\mathrm{M}_{\text{sxt}}=\begin{bmatrix}0&1&\Gamma\\\ 1&2\,\epsilon&0\\\
\Gamma&0&0\end{bmatrix},\qquad\qquad\mathrm{M}_{\text{oct}}=\begin{bmatrix}1&0&\Gamma\\\
0&2\,\epsilon&0\\\ \Gamma&0&0\end{bmatrix}.$

After eliminating all dependent parameters (Appendix A), we establish the
following correspondence:

$\begin{array}[]{llc}\displaystyle\mathcal{K}_{\mathrm{sxt}}[p,q]=\mathcal{K}_{0}+p^{2}q+q\,p^{2},&\,\,|a|\leq
2,&\mathrm{(SX)}\\\\[7.11317pt]
\displaystyle\mathcal{K}^{-}_{\mathrm{oct}}[p,q]=\mathcal{K}_{0}-p^{2}q^{2},&\,\,|a|\leq
2,\,\Gamma<0,&\mathrm{(DO)}\\\\[7.11317pt]
\displaystyle\mathcal{K}^{+}_{\mathrm{oct}}[p,q]=\mathcal{K}_{0}+p^{2}q^{2},&\,\,|a|\leq
2,\,\Gamma>0,&\mathrm{(FO)}\\\\[7.11317pt]
\displaystyle\mathcal{K}^{+}_{\mathrm{oct}}[p,q]=\mathcal{K}_{0}+p^{2}q^{2},&\,\,|a|>2,\,\Gamma>0,&\mathrm{(Df)}\end{array}$
(5)

where

$\mathcal{K}_{0}[p,q]=p^{2}-a\,p\,q+q^{2}.$

While the systems above are not precisely identical to the Hamiltonians (4),
we will observe substantial similarities between them, aiding in constructing
a comprehensive understanding of the symmetric McMillan map dynamics. For
instance, (I) both Hamiltonians and mappings exhibit similar amplitude-
frequency dependencies below the intrinsic resonance. (II) The solutions for
$q$ involve Jacobi elliptic functions: $\mathrm{sn}$, $\mathrm{cn}$, and
$\mathrm{dn}$ for (DO), (FO), and (Df), respectively Iatrou and Roberts
(2002), and rational function of $\mathrm{sn^{2}}$ for the (SX) case Zolkin
_et al._ (2022). Most importantly, we demonstrate that, akin to Hamiltonian
(3) being an approximation for a more general potentials $U(q)$, McMillan
multipoles serve as first and second-order approximations for the map in a
McMillan-Hénon form (1) with a smooth $f(q)$. This is particularly useful for
“integrating out” leading nonlinear effects around main resonances for chaotic
systems, akin to Hénon Henon (1969) or Chirikov standard mappings Chirikov
(1969, 1979).

##  III Solving the Equations of Motion

###  III.1 Sextupole map

Here, we will illustrate the general method of obtaining a parametrization of
an individual curve and a set of action-angle variables based on Danilov’s
theorem Zolkin  _et al._ (2017); Nagaitsev and Zolkin (2020); Mitchell  _et
al._ (2021), using McMillan sextupole as an example. We limit ourselves to
cases with stable trajectories around fixed point at the origin, considering
only $-2<a<2$.

$\bullet$ We start by introducing a formal Hamiltonian

$\mathrm{H}[p,q;t]\equiv\mathcal{K}_{\mathrm{sxt}}[p,q]=p^{2}q+p\,q^{2}+p^{2}-a\,p\,q+q^{2}$

which satisfies the system of Hamilton’s equations

$\begin{array}[]{r}\displaystyle\qquad\frac{\mathrm{d}q}{\mathrm{d}t}=\frac{\partial\mathrm{H}}{\partial
p}=q^{2}-a\,q+2\,p+2\,p\,q,\\\\[9.95863pt]
\displaystyle-\frac{\mathrm{d}p}{\mathrm{d}t}=\frac{\partial\mathrm{H}}{\partial
q}=p^{2}-a\,p+2\,q+2\,p\,q.\end{array}$ (6)

$\bullet$ Since Hamiltonian function does not change along the trajectory, we
can solve for $p$ from $\mathrm{H}[p,q]=\mathrm{const}$:

$p=\frac{1}{2}\,\left[f_{\text{sxt}}(q)\pm\frac{\sqrt{\mathcal{P}(q)}}{q+1}\right],$
(7)

where

$\mathcal{P}(q)=q^{4}-2\,(2+a)\,q^{3}+(a^{2}-4)\,q^{2}+4\,\mathrm{H}\,q+4\,\mathrm{H}.$

When the fixed point at the origin is stable ($|a|<2$), for all closed
invariant curves we have $\mathrm{H}>0$, and, polynomial $\mathcal{P}(q)$ has
four distinct real roots $q_{i}$ such that

$\sum_{i=1}^{4}\frac{1}{q_{i}}=-1,\qquad\qquad q_{1}<q_{2}\leq q\leq
q_{3}<q_{4}.$

Using Vieta’s formulas for the quartic $\mathcal{P}(q)$, one can express map
parameter and energy level using $q_{i}$ as:

$a=\sum_{i=1}^{4}\frac{q_{i}}{2}-2,\qquad\qquad\qquad\mathrm{H}=\frac{1}{4}\prod_{i=1}^{4}q_{i}.$

$\bullet$ Next, we should relate continuous flow given by Hamiltonian with the
mapping equations. Substitution of (7) into the first Hamilton’s equation (6)
provides

$\mathrm{d}\,t=\pm\mathrm{d}\,q/\sqrt{\mathcal{P}(q)}.$ (8)

Taking the integral of Eq. (8) along the invariant curve, we can find the
period of motion

$\mathrm{T}\equiv\oint\mathrm{d}\,t=2\,\int_{q_{2}}^{q_{3}}\frac{\mathrm{d}\,q}{\sqrt{\mathcal{P}(q)}}=\frac{4\,\mathrm{K}[\kappa]}{\sqrt{(q_{4}-q_{2})(q_{3}-q_{1})}}$

and the time of one-step of the map

$\displaystyle\mathrm{T^{\prime}}\equiv\int_{0}^{\mathrm{T^{\prime}}}\mathrm{d}\,t=\int_{q_{0}}^{q_{0}^{\prime}}\frac{\mathrm{d}\,q}{\sqrt{\mathcal{P}(q)}}=\int_{q_{2,3}}^{q_{2,3}^{\prime}}\frac{\mathrm{d}\,q}{\sqrt{\mathcal{P}(q)}}.$

Figure 1:  Phase space diagram for the McMillan sextupole map with $a=-3/2$
(left), and parametrization (9) of invariant curve with initial conditions
$(\\{q_{0}\\},\\{p_{0}\\})=(1/3,0)$ (right).

The last integral is independent of choice of lower bound, so we can choose
one of the turning points as $q_{0}$, where with the help of mapping equations
we have

$q_{2,3}^{\prime}=p(q_{2,3})=f_{\text{sxt}}(q_{2,3})/2=-q_{2,3}\,\frac{q_{2,3}-a}{q_{2,3}+1}.$

$\text{T}^{\prime}$ provides the discretization time interval and defines the
rotation number of the map as the ratio

$\nu=\frac{\mathrm{T}^{\prime}}{\mathrm{T}}=\frac{\mathrm{F}\left[\arcsin\sqrt{\frac{q_{3}-q_{1}}{q_{3}-q_{2}}\frac{q_{2}(3\,q_{2}+2-a)}{2\,q_{1}(q_{2}+1)+q_{2}\,(q_{2}-a)}},\kappa\right]}{2\,\mathrm{K}[\kappa]}.$

Here $\mathrm{F}[\phi,\kappa]$ and $\mathrm{K}[\kappa]$ are complete and
incomplete elliptic integrals of the first kind with elliptic modulus

$\kappa=\sqrt{\frac{(q_{3}-q_{2})(q_{4}-q_{1})}{(q_{3}-q_{1})(q_{4}-q_{2})}}.$

$\bullet$ Taking integral (8) from $q(0)=\\{q_{0}\\}$ to $q(t)$, we can obtain
parametrization of the map
$\\{q_{n}\\},\\{p_{n-1}\\}=q(n\,\mathrm{T}^{\prime})$:

$q(t)=\frac{q_{2}-q_{1}\,\frac{q_{3}-q_{2}}{q_{3}-q_{1}}\,\mathrm{sn}_{1}^{2}}{1-\frac{q_{3}-q_{2}}{q_{3}-q_{1}}\,\mathrm{sn}_{1}^{2}}=\frac{q_{3}-q_{4}\,\frac{q_{3}-q_{2}}{q_{4}-q_{2}}\,\mathrm{sn}_{2}^{2}}{1-\frac{q_{3}-q_{2}}{q_{4}-q_{2}}\,\mathrm{sn}_{2}^{2}},$
(9)

where $\mathrm{sn}_{1,2}$ stands for one of two possible Jacobi elliptic sine
functions (see Fig. 1)

$\mathrm{sn}_{i}=\mathrm{sn}\left[\sqrt{(q_{4}-q_{2})(q_{3}-q_{1})}\,\frac{t-t^{*}_{i}}{2},\kappa\right]$

with initial phase shifts

$\displaystyle t^{*}_{1,2}=\mp
2\,\frac{\mathrm{F}\left[\arcsin\sqrt{\frac{(q_{3,4}-q_{1,2})(\\{q_{0}\\}-q_{2,3})}{(q_{3}-q_{2})(\\{q_{0}\\}-q_{1,4})}},\kappa\right]}{\sqrt{(q_{4}-q_{2})(q_{3}-q_{1})}}.$

$\bullet$ Finally, we can rewrite the map in its canonical form

$\begin{array}[]{ll}J^{\prime}=J,&\qquad\\{J_{n}\\}=\\{J_{0}\\},\\\\[7.11317pt]
\psi^{\prime}=\psi+2\,\pi\,\nu(J),&\qquad\\{\psi_{n}\\}=\\{\psi_{0}\\}+2\,\pi\,n\,\nu(\\{J_{0}\\}),\end{array}$

where $J$ and $\psi$ are action-angle variables. The action integral is given
in Appendix C, while initial phase $\\{\psi_{0}\\}$ can be chosen e.g., as
$\arctan(\\{p_{0}\\}/\\{q_{0}\\})$ or $-2\,\pi\,t^{*}_{1,2}/\mathrm{T}$.

|  |  DO ($\Gamma<0$, $|a|<2$) |  FO ($\Gamma>0$) |  Df ($\Gamma>0$, $|a|>2$)
---|---|---|---|---
Elliptic function |  ef |  sn |  cn |  dn
Domain of invariant |  $\mathcal{K}\in$ |  $\left[0;(|\frac{a}{2}|-1)^{2}\right]\geq 0$ |  $[0;\infty)\geq 0$ |  $\left[-(|\frac{a}{2}|-1)^{2};0\right]\leq 0$
Amplitude |  $h(\kappa)$ |  $\displaystyle\sqrt{k}$ |  $\displaystyle\sqrt{k/k^{\prime}}$ |  $\displaystyle\sqrt{1/k^{\prime}}$
Phase advance ($a>0$) |  $\eta_{+}$ |  $\displaystyle\mathrm{arcns}\,\frac{\sqrt{k}}{\sqrt[4]{\mathcal{K}}}$ |  $\displaystyle\mathrm{arcds}\,\frac{\sqrt{k\,k^{\prime}}}{\sqrt[4]{\mathcal{K}}}$ |  $\displaystyle\mathrm{arccs}\,\frac{\sqrt{k^{\prime}}}{\sqrt[4]{|\mathcal{K}|}}$
Initial phase |  $\\{u_{0}\\}$ |  $\displaystyle\mathrm{arcsn}\,\frac{\\{q_{0}\\}}{h(\kappa)\,\sqrt[4]{\mathcal{K}}}$ |  $\displaystyle\mathrm{arccn}\,\frac{\\{q_{0}\\}}{h(\kappa)\,\sqrt[4]{\mathcal{K}}}$ |  $\displaystyle\mathrm{arcdn}\,\frac{\\{q_{0}\\}}{h(\kappa)\,\sqrt[4]{|\mathcal{K}|}}$
Elliptic modulus |  $k(B)$ |  $\displaystyle\frac{B-\sqrt{B^{2}-4}}{2}$ |  $\displaystyle\frac{1}{\sqrt{2}}\sqrt{1+\frac{B}{\sqrt{B^{2}+4}}}$ |  $\displaystyle\frac{\sqrt{B+2}-\sqrt{B-2}}{2\,(B^{2}-4)^{-1/4}}$
Complimentary modulus |  $k^{\prime}(B)$ |  $\displaystyle\frac{\sqrt{B+2}-\sqrt{B-2}}{2\,(B^{2}-4)^{-1/4}}$ |  $\displaystyle\frac{1}{\sqrt{2}}\sqrt{1-\frac{B}{\sqrt{B^{2}+4}}}$ |  $\displaystyle\frac{B-\sqrt{B^{2}-4}}{2}$
$B$-function |  $B(a,\mathcal{K})$ |  $\displaystyle-\frac{(a/2)^{2}-1-\mathcal{K}}{\sqrt{\mathcal{K}}}$ |  $\displaystyle\frac{(a/2)^{2}-1+\mathcal{K}}{\sqrt{\mathcal{K}}}$ |  $\displaystyle\frac{(a/2)^{2}-1+\mathcal{K}}{\sqrt{|\mathcal{K}|}}$
|  |  $\displaystyle=\frac{1}{k}+k\geq 2$ |  $\displaystyle=\frac{k}{k^{\prime}}-\frac{k^{\prime}}{k}\in\mathbb{R}$ |  $\displaystyle=\frac{1}{k^{\prime}}+k^{\prime}\geq 2$
Table 1: Elliptic parametrization of stable trajectories for McMillan octupole
map (after Iatrou and Roberts (2002)); for $a<0$ use
$\eta_{-}=2\,\mathrm{K}[\kappa]-\eta_{+}$.

###  III.2 Octupole map

The parametrization of individual curves for the octupole map was obtained by
Iatrou and Roberts in Iatrou and Roberts (2002). Additionally, they provided
two different methods demonstrating how an arbitrary constant level set of an
asymmetric McMillan map can be transformed first into a symmetric, and then
into the octupole (canonical) form. Consequently, the authors described the
dynamics on each symmetric biquadratic
$\mathcal{K}^{\mp}_{\mathrm{oct}}[p,q]=\mathrm{const}$ using the elliptic
modulus $\kappa$ and the argument $\\{u_{n}\\}=\\{u_{0}\\}+n\,\eta$ of the
elliptic function as a pair of “action-angle” variables. However, despite
capturing the essence of the dynamics, the change of variables from $(q,p)$ to
$(\kappa,u)$ does not conserve the phase space area. Below, we summarize their
results using our notations and complement it with a form-invariant set of
canonical action-angle variables Kolmogorov (1954); Moser (1962); Arnol’d
(1963), thereby extending the description from an individual curve to a
continuum of stable trajectories.

Due to additional symmetries of the invariant, the functional dependence is
much simpler compared to the sextupole case and all trajectories can be
parametrized as follows:

$\left\\{\begin{array}[]{l}\displaystyle\\{q_{n}\\}=\sqrt[4]{|\mathcal{K}|}\,h(\kappa)\,\mathrm{ef}\left[\\{u_{0}\\}+\eta\,n,\kappa\right],\\\\[7.11317pt]
\\{p_{n}\\}=\\{q_{n+1}\\},\end{array}\right.$ (10)

where $\mathrm{ef}$ is an appropriate elliptic function with elliptic modulus
$\kappa$, $h(\kappa)$ is an amplitude function, $\eta$ is the phase advance,
and, $\\{u_{0}\\}$ is the initial phase shift defined by initial condition
$\\{q_{0}\\}$. For defocusing octupole ($\Gamma<0$), the proper Jacobi
function for stable trajectories is elliptic sine $\mathrm{sn}$, while for
focusing system ($\Gamma>0$) one should use elliptic cosine $\mathrm{cn}$ for
trajectories rounding the origin, or delta amplitude $\mathrm{dn}$ for closed
orbits inside the figure-eight/lemniscate shape separatrix (case $|a|>2$),
akin regimes with finite motion in classical Hamiltonians

$\mathcal{H}[p,q;t]=\frac{p^{2}}{2}+k\,\frac{q^{2}}{2}\mp\frac{q^{4}}{4}.$

Expressions for all parameters are listed in Table 1 with
$\eta_{-}=2\mathrm{K}[\kappa]-\eta_{+}$.

Using parametrization (10), we can evaluate the integral for canonical action
as in zol (2017):

$\displaystyle J_{\text{ef}}=\frac{1}{2\,\pi}\,\oint
p\,\mathrm{d}q=\frac{\sqrt{|\mathcal{K}|}}{2\,\pi}\times\left\\{\begin{array}[]{l}\displaystyle\qquad\kappa\,\,S_{\mathrm{sn}},\\\\[5.69046pt]
\displaystyle(\kappa/\kappa^{\prime})\,S_{\mathrm{cn}},\\\\[5.69046pt]
\displaystyle(1/\kappa^{\prime})\,S_{\mathrm{dn}},\end{array}\right.$ (11)

where $S_{\mathrm{ef}}$ is the area of corresponding elliptic Lissajous curve
with matching frequencies and phase difference $\eta$.

Finally, by rescaling the argument
$u\rightarrow\psi=2\,\pi\,u/\mathrm{T}_{\text{ef}}$, we define the canonical
angle variable

$\\{\psi_{n}\\}=\\{\psi_{0}\\}+2\,\pi\,\nu\,n,$

where corresponding period of elliptic functions $\mathrm{T}_{\text{ef}}$ is
equal to $4\,\mathrm{K}[\kappa]$ for oscillations rounding the origin
($\mathrm{sn}$ and $\mathrm{cn}$), or, $2\,\mathrm{K}[\kappa]$ for
trajectories inside the figure-eight curve (dn). The Poincaré rotation number
again given by $\mathrm{T^{\prime}}/\mathrm{T_{\text{ef}}}$ where one-step
time interval $\mathrm{T^{\prime}}\equiv\eta$ such that

$\eta=\left\\{\begin{array}[]{ll}\eta_{+}&\quad\text{for }a>0,\\\\[7.11317pt]
2\,\mathrm{K}[\kappa]-\eta_{+}&\quad\text{for }a<0.\end{array}\right.$

Explicit expressions for $J_{\text{ef}}$, $S_{\text{ef}}$ and power series
$\nu(J_{\text{ef}})$ are given in Appendix C.

Figure 2:  Rotation number as a function of the action variable $\nu(J)$ for
McMillan sextupole (SX), defocusing (DO) and focusing (FO) octupoles, along
with a focusing octupole in the Duffing regime (Df) with $a\geq 2$. The blue
diagrams schematically depict phase space portraits corresponding to different
curves, with the values of the map parameter $a$ indicated in red. The action
variable for trajectories inside the figure-eight separatrix (case Df) is
multiplied by a factor of 2 to ensure continuity of the graph with large
amplitude trajectories rounding the origin.

##  IV Amplitude-frequency dependence

In this section, we explore intrinsic dynamical properties that remain
invariant regardless of the mapping’s representation and examine limiting
cases of both small and large amplitudes. Figure 2 illustrates the dependence
of the rotation number on the action variable $\nu(J)$ for all the cases under
consideration. To deepen our understanding, we compare the qualitative
behavior of mappings with our reference classical Hamiltonians (4) by setting
Fig. 2 side by side with corresponding dependencies of frequency on action
$\omega(J)$, as depicted in Figure 3. For the readers’ convenience, all exact
expressions along with their power series are listed at the end of the article
in Appendix C.

$\bullet$ Small amplitude oscillations can be analyzed using a linearization
approach, where the system’s behavior near an equilibrium or fixed point is
approximated. In all cases with stable orbits around the origin (SX, DO, and
FO), such motion is determined by quadratic terms:

$\mathcal{H}_{0}[p,q;t]=\frac{p^{2}}{2}+\frac{q^{2}}{2},\qquad\quad\mathcal{K}_{0}[p,q]=p^{2}-a\,p\,q+q^{2},$

for Hamiltonians or mappings, respectively. This yields the corresponding
frequency and rotation number evaluated at $J=0$:

$\omega_{0}=1,\qquad\qquad\qquad\qquad\nu_{0}=\frac{1}{2\,\pi}\,\arccos\frac{a}{2}.$

For the Duffing oscillator, linearization around stable symmetric
stationary/fixed points provides:

$\omega_{0}=\sqrt{2},\qquad\qquad\nu_{0}=\frac{1}{2\,\pi}\,\arccos\left(\frac{4}{a}-1\right).$

Although we can always normalize $\omega_{0}$ to unity by adjusting units of
time for Hamiltonians, mappings correspond to systems with $1\frac{1}{2}$
degrees of freedom, resulting in an irreducible intrinsic parameter even for
the case of linear oscillations.

Figure 3:  Frequency as a function of the action variable, $\omega(J)$, for
all Hamiltonians (4). The corresponding phase space diagrams are displayed to
the right, color-coded consistently. The action variable for trajectories
inside the figure-eight separatrix (case Df) is multiplied by a factor of 2 to
ensure continuity of the graph with large amplitude trajectories rounding the
origin.

$\bullet$ When all trajectories are stable (cases FO and Df), the frequency of
oscillations for large amplitudes in Hamiltonian systems experiences unbounded
growth, as seen in the blue and cyan curves in Fig. 3. Conversely, with the
rotation number constrained to the range $\nu\in[0,1/2]$, the behavior of
$\nu(J)$ remains monotonic but may increase or decrease towards a value of
intrinsic resonance $\nu_{r}=1/4$, contingent upon whether the linear tune is
below $\nu_{0}<\nu_{r}$ or above $\nu_{0}>\nu_{r}$ with

$\lim_{J\rightarrow\infty}\omega(J)=\infty,\qquad\qquad\qquad\qquad\lim_{J\rightarrow\infty}\nu(J)=\nu_{r},$

as depicted in the two plots on the right in Fig. 2.

In all other instances (cases SX, DO, and small oscillations for Df), the
limiting amplitude is determined by the separatrix. In Hamiltonian dynamics,
this results in a vanishing frequency as the period on the separatrix tends to
infinity, while for mappings, the limiting rotation number tends to
$\frac{1}{2}$ for systems above the intrinsic resonance ($\nu_{0}>\nu_{r}$)
due to the presence of 2-cycles:

$\omega(J_{\text{sep}})=0,\qquad\nu(J_{\text{sep}})=\frac{1+\mathrm{sign}\,[\nu_{0}-\nu_{r}]}{4},\,\,(\nu_{0}\neq\nu_{r})$

where $\nu_{r}=1/3$ for the (SX) and $1/2$ for the (DO) cases. Notice that for
the respective resonant values of the parameter $a_{r}=2\cos(2\pi\nu_{r})$,
the mapping degenerates to linear with force function $f=a_{r}\,q$ and no
dependence of rotation number on amplitude:

$\nu(J)=\nu_{r},\qquad(\nu_{0}=\nu_{r}).$

Refer to the two left plots in Fig. 2.

$\bullet$ Within Hamiltonian systems (4), focusing and defocusing octupoles
display constant detuning values, assessed at the equilibrium point at the
origin, albeit with opposing signs:

$\left.\frac{\mathrm{d}\omega_{\mathrm{FO}}}{\mathrm{d}J_{\mathrm{FO}}}\right|_{J=0}=-\left.\frac{\mathrm{d}\omega_{\mathrm{DO}}}{\mathrm{d}J_{\mathrm{DO}}}\right|_{J=0}=\frac{3}{4}.$

The terms “focusing” and “defocusing” in the systems’ names correspond to
their optical properties. However, for McMillan octupoles, the detuning is a
function of the parameter $a$, with a reversal of optical properties upon
crossing the intrinsic resonance:

$\left.\frac{\mathrm{d}\nu_{\mathrm{FO,DO}}}{\mathrm{d}J_{\mathrm{FO,DO}}}\right|_{J=0}=\pm\frac{3}{2\,\pi}\frac{a}{4-a^{2}}=\pm\frac{3}{4\,\pi}\frac{\cot(2\,\pi\,\nu_{0})}{\sin(2\,\pi\,\nu_{0})}.$
(12)

Similarly, the detuning for McMillan sextupole becomes positive for
$\nu_{0}>1/3$

$\begin{array}[]{l}\displaystyle\left.\frac{\mathrm{d}\nu_{\mathrm{SX}}}{\mathrm{d}J_{\mathrm{SX}}}\right|_{J=0}=\displaystyle-\frac{1}{2\,\pi}\frac{(1+a)(8+a)}{(2+a)(2-a)^{2}}\\\\[11.38092pt]
\displaystyle\qquad=-\frac{1}{32\,\pi}\frac{8\,\sin(2\,\pi\,\nu_{0})+\sin(4\,\pi\,\nu_{0})}{\sin^{3}(\pi\,\nu_{0})\,\sin^{3}(2\,\pi\,\nu_{0})}\,\sin(3\,\pi\,\nu_{0})\end{array}$
(13)

in contrast to the always negative value in its reference Hamiltonian (SX):

$\left.\frac{\mathrm{d}\omega_{\mathrm{SX}}}{\mathrm{d}J_{\mathrm{SX}}}\right|_{J=0}=-\frac{5}{6}.$

Additionally, above the resonance, the invariant
$\mathcal{K}_{\text{sxt}}[p,q]$ has 3 critical points, making its level sets
topologically different from $\mathcal{H}_{\text{sxt}}[p,q]$ (see Appendix B
for details). Finally, considering the limit for the case (Df):

$\displaystyle\left.\frac{\mathrm{d}\nu_{\mathrm{Df}}}{\mathrm{d}J_{\mathrm{Df}}}\right|_{J=0}=\frac{1}{2\,\pi}\frac{4+a}{a\,(2-a)}=-\frac{1}{4\,\pi}\frac{2+\cos(2\,\pi\,\nu_{0})}{\tan^{2}(\pi\,\nu_{0})}.$
(14)

Figure 4 provides graphs for the detunings (12 – 14) and the second
derivatives of the rotation number with respect to the action variable at the
fixed point as functions of $a$ and $\nu_{0}$.

Figure 4:  Detuning $\mathrm{d}\nu/\mathrm{d}J$ and the second derivative
$\mathrm{d}^{2}\nu/\mathrm{d}J^{2}$ evaluated at the fixed point ($J=0$) for
all McMillan multipoles. The left and right columns show these quantities as
functions of the map parameter $a$ and $\nu_{0}$ respectively.

##  V Stability diagrams

To summarize the dynamical properties and establish their connection to
physical variables, we present stability diagrams, Figs. 5 and 6, utilizing
color to depict $\nu(q_{0})$ for finite trajectories at different parameter
values $a$. Each case includes two plots, corresponding to initial conditions
set up on the first (left column) and second (right column) symmetry lines,
respectively:

$l_{1}:\,p_{0}=q_{0},\qquad\quad\text{and}\quad\qquad
l_{2}:\,p_{0}=f(q_{0})/2.$

Figure 5:  Stability diagrams for McMillan sextupole. The left and right plots
correspond to initial conditions along the first or second symmetry lines. The
color represents the rotation number for stable trajectories, ranging from 0
(red) to 1/2 (violet). Tick marks at the top provide resonant values of
$\nu_{0}$.

Figure 6:  Same as in Fig. 5 but for defocusing (top) and focusing (bottom)
McMillan octupoles.

These lines traverse through key features, and in the context of accelerator
physics, the second diagram defines the dynamical aperture (the largest
possible coordinate), delineating the phase space available for the beam of
particles. Both symmetry lines intersect all fixed points, whereas the second
symmetry line and its inverse additionally pass through the 2-cycle. For more
details, please refer to McMillan (1971); DeVogelaere  _et al._ (1958); Lewis
Jr (1961); Zolkin  _et al._ (2023a).

##  VI McMillan map and perturbation theory

In Hamiltonian dynamics, a common technique involves expanding the potential
energy function around its minimum, especially when examining small
oscillations or perturbations. This expansion enables us to approximate the
potential and gain insights into the particle’s behavior in proximity to its
equilibrium position. For a smooth potential, denoted as $U(q)$ with a minimum
at the origin and satisfying $U(0)=0$, its Taylor expansion is given by:

$U(q)=\frac{U^{\prime\prime}(0)}{2!}\,q^{2}+\frac{U^{\prime\prime\prime}(0)}{3!}\,q^{3}+\frac{U^{\prime\prime\prime\prime}(0)}{4!}\,q^{4}+\mathcal{O}(q^{5}).$

As a result, the Hamiltonian $\mathcal{H}_{\mathrm{sxt}}[p,q]$ serves as the
first-order approximation for a generic asymmetric potential, while
$\mathcal{H}^{\pm}_{\mathrm{oct}}[p,q]$ are the first-order approximations for
a symmetric $U(q)$, or more generally, for an expansion with
$U^{\prime\prime\prime}(0)=0$. These model Hamiltonians play a crucial role in
capturing nonlinear effects, especially when higher-order terms in the
expansion come into play.

Beyond its significance as one of the very few known exact integrable models,
the symmetric McMillan map holds deeper meaning by serving as the first and
second-order approximations for a more general and even chaotic mappings in
McMillan-Hénon form (1), with a differentiable but otherwise arbitrary force
function $f(q)$, resembling situation in continuous dynamics. To demonstrate
this, we introduce a small positive parameter $\varepsilon$ representing the
amplitude of oscillations. This is achieved through a change of variables
$(q,p)\rightarrow\varepsilon\,(q,p)$:

$\begin{array}[]{l}q^{\prime}=p\\\\[7.11317pt] \displaystyle
p^{\prime}=-q+\frac{f(\varepsilon\,p)}{\varepsilon}=-q+a\,p+\varepsilon\,b\,p^{2}+\varepsilon^{2}c\,p^{3}+\ldots\end{array}$
(15)

where the force function is expanded in a power series of $(\varepsilon\,p)$

$a\equiv\partial_{p}f(0)/1!,\quad b\equiv\partial_{p}^{2}f(0)/2!,\quad
c\equiv\partial_{p}^{3}f(0)/3!,\quad\ldots,$

and we assume the fixed point to be at the origin, necessitating $f(0)=0$.
Subsequently, we seek an approximated invariant of motion that is conserved
with an accuracy of order $\mathcal{O}(\varepsilon^{n+1})$:

$\mathcal{K}^{(n)}[p^{\prime},q^{\prime}]-\mathcal{K}^{(n)}[p,q]=\mathcal{O}(\varepsilon^{n+1}).$
(16)

The invariant is sought in the form of a polynomial:

$\mathcal{K}^{(n)}=\mathcal{K}_{0}+\varepsilon\,\mathcal{K}_{1}+\varepsilon^{2}\,\mathcal{K}_{2}+\ldots+\varepsilon^{n}\,\mathcal{K}_{n},$

where $\mathcal{K}_{m}$ consists of homogeneous polynomials in $p$ and $q$ of
$(m+2)$ degree

$\begin{array}[]{l}\mathcal{K}_{0}=C_{2,0}\,p^{2}+C_{1,1}\,p\,q\,+C_{0,2}\,q^{2},\\\\[7.11317pt]
\mathcal{K}_{1}=C_{3,0}\,p^{3}+C_{2,1}\,p^{2}q+C_{1,2}\,p\,q^{2}+C_{0,3}\,q^{3},\\\\[7.11317pt]
\cdots,\end{array}$

and $C_{i,j}$ are coefficients to be determined to satisfy Eq. (16). The
reader can check that, in the first two orders of this perturbation theory, a
general result is provided

$\begin{array}[]{l}\displaystyle\mathcal{K}^{(2)}[p,q]=\mathcal{K}_{0}[p,q]-\varepsilon\,\frac{b}{a+1}\,(p^{2}q+p\,q^{2})\,+\\\\[9.95863pt]
\qquad\displaystyle+\,\varepsilon^{2}\left(\left[\frac{b^{2}}{a\,(a+1)}-\frac{c}{a}\right]p^{2}q^{2}+C\,\mathcal{K}_{0}^{2}[p,q]\right)\end{array}$
(17)

where $C$ is a coefficient such that Eq. (16) is satisfied for any value it
takes.

Here, if we set $C=0$, one recognizes the invariant for the symmetric McMillan
map. Thus, we learn that mappings (5) are not just integrable models. Similar
to how the quadratic invariant $\mathcal{K}_{0}[p,q]$ with the corresponding
force $f(p)=a\,p$ allows us to “integrate out” linear dynamics around the
fixed point (zeroth order), McMillan sextupoles and octupoles provide the next
order of approximation, defining detuning and even approximating the dynamical
aperture around low-order resonances. Below, we separately analyze cases of
small ($\varepsilon\rightarrow 0$) and large amplitudes, providing several
examples and discussing applicability.

###  VI.1 Small amplitudes

####  VI.1.1 McMillan mapping in Duffing regime

As our first example, let’s consider the familiar case of the McMillan
octupole (5) in the Duffing regime (Df). Since the fixed point
$\zeta^{(1-1)}=(0,0)$ is unstable, we can shift the origin to one of the
symmetric fixed points (25):

$(q,p)\rightarrow(\widetilde{q},\widetilde{p})=(q,p)-\zeta^{(1-2)},\quad\zeta^{(1-2)}=\sqrt{\frac{a-2}{2}}\,(1,1).$

This preserves the form of the map (1) and results in the transformation of
the force function:

$f(q)=\frac{a\,q}{1+q^{2}}\,\,\rightarrow\,\,\widetilde{f}(\,\\!\widetilde{q}\,)=f\left(\widetilde{q}+q^{(1-2)}\right)-2\,q^{(1-2)}.$

Expanding the new force in a series

$\begin{array}[]{l}\displaystyle\widetilde{f}(\,\\!\widetilde{q}\,)=\alpha\,\widetilde{q}-\frac{(1+\alpha)\,\sqrt{4-\alpha^{2}}}{2}\,\widetilde{q}^{\,2}+\\\\[9.95863pt]
\displaystyle\qquad\qquad\,\,\,\,\,+\,\frac{(2+\alpha)(2-\alpha^{2})}{4}\,\widetilde{q}^{\,3}+\mathcal{O}(\widetilde{q}^{\,4}),\end{array}$

where we introduced the parameter

$\alpha=\frac{8}{a}-2\,:\qquad\nu_{0}=\frac{\arccos(\alpha/2)}{2\,\pi},$ (18)

we can write an invariant in the form (17)

$\mathcal{K}[\widetilde{p},\widetilde{q}\,]=\mathcal{K}_{0}[\widetilde{p},\widetilde{q}\,]+\frac{\sqrt{4-\alpha^{2}}}{2}\,(\widetilde{p}^{\,2}\widetilde{q}+\widetilde{p}\,\widetilde{q}^{\,2})+\frac{2+\alpha}{4}\,\widetilde{p}^{\,2}\widetilde{q}^{\,2}.$

In this example, the second-order integral of motion from above is an exact
invariant (i.e., not approximated), which, when written in new coordinates,
provides coefficients in front of two basic nonlinearities. We see that
oscillations inside the figure-eight separatrix are equivalent to the
symmetric McMillan map carrying both sextupole (SX) and octupole (FO)
components. Thus, we can find detuning by adding Eqs. (13) and (12) that are
multiplied by appropriate scaling coefficients from the invariant:

$\begin{array}[]{l}\displaystyle\left.\frac{\mathrm{d}\nu_{\mathrm{Df}}}{\mathrm{d}J_{\mathrm{Df}}}\right|_{J=0}=-\frac{1}{2\,\pi}\frac{(1+\alpha)(8+\alpha)}{(2+\alpha)(2-\alpha)^{2}}\times\left(\frac{\sqrt{4-\alpha^{2}}}{2}\right)^{2}+\\\\[14.22636pt]
\displaystyle\qquad+\,\frac{3}{2\,\pi}\frac{\alpha}{4-\alpha^{2}}\times\left(\frac{2+\alpha}{4}\right)=-\frac{1}{8\,\pi}\frac{(2+\alpha)(4+\alpha)}{2-\alpha}.\end{array}$

Using Eq. (18), we can convert this to the form independent of the mapping’s
representation:

$\displaystyle\left.\frac{\mathrm{d}\nu_{\mathrm{Df}}}{\mathrm{d}J_{\mathrm{Df}}}\right|_{J=0}=-\frac{1}{4\,\pi}\frac{2+\cos(2\,\pi\,\nu_{0})}{\tan^{2}(\pi\,\nu_{0})}$

which matches the previously obtained result (14). Notice that this time we
didn’t solve for action-angle variables, then take a limit of corresponding
elliptic functions, but rather used our knowledge of underlying
nonlinearities.

####  VI.1.2 Hénon quadratic and qubic mappings

It is important to recognize that for an arbitrary $f(p)$, the global dynamics
of equation (1) is known to be chaotic Henon (1969); Chirikov (1969, 1979),
except for a few very special integrable cases Veselov (1991); McMillan
(1971); Quispel  _et al._ (1988, 1989); Suris (1989); Brown (1993); Cairns
_et al._ (2016); Zolkin  _et al._ (2023a, b). For a brief review on the
integrability of mappings in the form (1), please refer to Zolkin  _et al._
(2023a).

As a case study for chaotic dynamics, we use the aforementioned area-
preserving Hénon mappings (H) with quadratic Henon (1969) and cubic force
functions:

$f^{\text{(H)}}_{\text{sxt}}(q)=a\,q+q^{2},\qquad\qquad\qquad
f^{\text{(H)}}_{\text{oct}}(q)=a\,q\pm q^{3}.$ (19)

$\bullet$ Starting with the octupole forces, by comparing
$f^{\text{(H)}}_{\text{oct}}$ with Eq. (15), we have $b=0$ and $c=\pm 1$,
providing an approximated McMillan-Hénon (MH) invariant of the second order:

$\mathcal{K}_{\text{DO,FO}}^{\text{(MH)}}[p,q]=\mathcal{K}_{0}[p,q]\mp\frac{1}{a}\,p^{2}q^{2}$

corresponding to symmetric McMillan mappings with forces

$f_{\text{DO,FO}}^{\text{(MH)}}(q)=\frac{a^{2}q}{a\mp q^{2}}=a\,q\pm
q^{3}+\mathcal{O}(q^{5}).$

Therefore, with the use of the scaling, from Eq. (12) we have

$\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(MH)}}_{\text{DO,FO}}}{\mathrm{d}J^{\text{(MH)}}_{\text{DO,FO}}}\right|_{J=0}=\frac{1}{a}\times\displaystyle\left.\frac{\mathrm{d}\nu_{\text{DO,FO}}}{\mathrm{d}J_{\text{DO,FO}}}\right|_{J=0}=\mp\frac{3}{8\,\pi}\frac{1}{\sin^{2}(2\,\pi\,\nu_{0})}.$
(20)

$\bullet$ For the quadratic Hénon map, we can also go up to the second order
of perturbation theory by first matching the quadratic term in the force
function ($b=1$)

$f^{\text{(MH)}}_{\text{SX-1}}(q)=\frac{a\,(a+1)+q}{(a+1)-q}\,q=a\,q+q^{2}+\frac{q^{3}}{a+1}+\mathcal{O}(q^{4}),$

and then removing the cubic term from the expansion ($c=0$) by using

$f^{\text{(MH)}}_{\text{SX-2}}(q)=\frac{a\,(a+1)+q}{(a+1)-q+\frac{1}{a}\,q^{2}}\,q=a\,q+q^{2}+\mathcal{O}(q^{4}).$

While the first order is simply a rescaled McMillan sextupole (SX), the second
order represents a specific mixture of both McMillan sextupole and focusing
octupole (FO), as can be seen from the invariants:

$\begin{array}[]{l}\displaystyle\mathcal{K}_{\text{SX-1}}^{\text{(MH)}}[p,q]=\mathcal{K}_{0}[p,q]-\frac{p^{2}q+p\,q^{2}}{a+1},\\\\[11.38092pt]
\displaystyle\mathcal{K}_{\text{SX-2}}^{\text{(MH)}}[p,q]=\mathcal{K}_{0}[p,q]-\frac{p^{2}q+p\,q^{2}}{a+1}+\frac{p^{2}\,q^{2}}{a\,(a+1)}.\end{array}$

Applying the same logic, we evaluate the detuning for the quadratic Hénon map
as

$\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(MH)}}_{\text{SX-1}}}{\mathrm{d}J^{\text{(MH)}}_{\text{SX-1}}}\right|_{J=0}=\left(\frac{-1}{1+a}\right)^{2}\times\displaystyle\left.\frac{\mathrm{d}\nu_{\text{SX}}}{\mathrm{d}J_{\text{SX}}}\right|_{J=0}$

and

$\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(MH)}}_{\text{SX-2}}}{\mathrm{d}J^{\text{(MH)}}_{\text{SX-2}}}\right|_{J=0}=\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(MH)}}_{\text{SX-1}}}{\mathrm{d}J^{\text{(MH)}}_{\text{SX-1}}}\right|_{J=0}+\frac{1}{a\,(1+a)}\times\displaystyle\left.\frac{\mathrm{d}\nu_{\text{FO}}}{\mathrm{d}J_{\text{FO}}}\right|_{J=0}$

resulting in

$\begin{array}[]{l}\displaystyle\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(MH)}}_{\text{SX-1}}}{\mathrm{d}J^{\text{(MH)}}_{\text{SX-1}}}\right|_{J=0}=-\frac{1}{16\,\pi}\frac{9\,\cos(\pi\,\nu_{0})+\cos(3\,\pi\,\nu_{0})}{\sin^{3}(2\,\pi\,\nu_{0})\,\sin(3\,\pi\,\nu_{0})},\\\\[15.6491pt]
\displaystyle\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(MH)}}_{\text{SX-2}}}{\mathrm{d}J^{\text{(MH)}}_{\text{SX-2}}}\right|_{J=0}=-\frac{1}{16\,\pi}\frac{3\,\cot(\pi\,\nu_{0})+\cot(3\,\pi\,\nu_{0})}{\sin^{3}(2\,\pi\,\nu_{0})}.\end{array}$
(21)

To check the accuracy of Eqs. (21) and (20) in describing actual dynamics, we
compare them with the dependencies obtained numerically from tracking. Despite
the chaotic nature of the global dynamics with forces (19), the KAM theorem
assures the existence of action-angle variables in the vicinity of stable
fixed points Kolmogorov (1954); Moser (1962); Arnol’d (1963). We introduce a
test particle that is initially placed, for example, on one of the symmetry
lines with a small deviation from the origin, $q_{0}=\delta q$. First, we
perform tracking for $N$ iterations:

$\\{q_{0},p_{0}\\}\rightarrow\\{q_{1},p_{1}\\}\rightarrow\ldots\rightarrow\\{q_{N},p_{N}\\},$

and then, we rewrite it in polar phase space coordinates

$\\{\rho_{0},\phi_{0}\\}\rightarrow\\{\rho_{1},\phi_{1}\\}\rightarrow\ldots\rightarrow\\{\rho_{N},\phi_{N}\\}$

where

$q=\rho\,\cos\phi,\qquad\qquad\qquad p=\rho\,\sin\phi.$

From the sub-array of polar angles $\\{\phi_{0},\phi_{q},\ldots,\phi_{N}\\}$,
we evaluate the rotation number as:

$\nu^{*}(\delta
q)=\frac{1}{2\,\pi}\,\sum_{j=0}^{N-1}\frac{\phi_{j+1}-\phi_{j}}{N}.$

and by sorting tuples $\\{\rho_{j},\phi_{j}\\}$ such that $\forall\,j\in[1,N]$
we have $\phi_{j}>\phi_{j-1}$, we estimate the action integral according to

$J^{*}(\delta
q)=\frac{1}{2\,\pi}\,\sum_{j=0}^{N-1}\,\frac{\rho_{j}^{2}}{2}\,(\phi_{j+1}-\phi_{j}).$

Then, detuning can be numerically evaluated as:

$\left.\frac{\mathrm{d}\nu}{\mathrm{d}J}\right|_{J=0}=\lim_{\begin{subarray}{c}\delta
q\rightarrow 0\,\,\,\\\ N\rightarrow\infty\end{subarray}}\frac{\nu^{*}(\delta
q)-\nu_{0}}{J^{*}(\delta q)}.$

By choosing a sufficiently large $N$ and sufficiently small $\delta q$, we
ensure convergence.

Figure 7:  Detuning for the Hénon quadratic (left) and cubic (right) mappings
(19), assessed from tracking data (black curves), juxtaposed with analytical
predictions from perturbation theory (illustrated with colors).

Figure 7 presents the comparison of tracking (black curve) along with
analytical results (shown in colors). As observed, the first order of
perturbation theory (SX-1) for the Hénon quadratic map is not sufficient and
provides an accurate approximation of detuning only in the range
$0\leq\nu_{0}\lesssim 0.2$. However, the second order for both quadratic and
cubic mappings (SX-2, DO and FO respectively) provides coincidence between
both approaches up to a machine-order accuracy.

####  VI.1.3 Thin nonlinear lens

Let’s now consider a more practical example of a simple accelerator lattice
with one degree of freedom consisting of linear optics elements (drift spaces,
dipoles, and quadrupoles) and a single thin nonlinear lens Lee (2018):

$\mathrm{F}:\quad\,\,\begin{bmatrix}x\\\
\dot{x}\end{bmatrix}^{\prime}=\begin{bmatrix}x\\\
\dot{x}\end{bmatrix}+\begin{bmatrix}0\\\ F(x)\end{bmatrix}.$

The effect on a test particle from all linear elements can be represented
using a matrix with Courant-Snyder parametrization Courant and Snyder (1958):

$\mathrm{M}:\quad\begin{bmatrix}x\\\
\dot{x}\end{bmatrix}^{\prime}=\begin{bmatrix}\cos\Phi+\alpha\,\sin\Phi&\beta\,\sin\Phi\\\
-\gamma\,\sin\Phi&\cos\Phi-\alpha\,\sin\Phi\end{bmatrix}\begin{bmatrix}x\\\
\dot{x}\end{bmatrix},$

where $\alpha$, $\beta$ and $\gamma$ are Twiss parameters (also known as
Courant-Snyder parameters) at the thin lens location, and $\Phi$ is the
betatron phase advance over the linear optics insert

$\Phi=\int\frac{\mathrm{d}s}{\beta(s)}.$

Without the nonlinear lens, the Twiss parameters are functions of the
longitudinal coordinate $s$ with $\beta(s)$ referred to as the
$\beta$-function, $\alpha(s)\equiv-\dot{\beta}(s)/2$, and
$\gamma(s)\equiv[1+\alpha^{2}(s)]/\beta(s)$. At any location, the Courant-
Snyder invariant is defined as:

$\gamma(s)\,x^{2}(s)+2\,\alpha(s)\,x(s)\,\dot{x}(s)+\beta(s)\,\dot{x}^{2}(s)=\mathrm{const}.$

The rotation number (or betatron tune in accelerator physics) is independent
of amplitude and given by:

$\nu_{0}=\frac{1}{2\,\pi}\,\oint\frac{\mathrm{d}s}{\beta(s)}.$

When the nonlinear lens is introduced, the combined one-turn map
$\mathrm{M}\circ\mathrm{F}$ can be rewritten in the form (1) using a change of
variables

$\left\\{\begin{array}[]{l}q=x,\\\\[5.69046pt]
p=x\,(\cos\Phi+\alpha\,\sin\Phi)+\dot{x}\,\beta\,\sin\Phi,\end{array}\right.$
(22)

with the force function given by

$f(q)=2\,q\,\cos\Phi+\beta\,\sin\Phi\,F(q).$

Using results from the previous subsection, we know that we can define an
approximated integral of motion (17) with

$\begin{array}[]{l}\displaystyle
a=2\,\cos\Phi+\beta\,\sin\Phi\,\partial_{q}\,F(0),\\\\[7.11317pt]
\displaystyle b=\beta\,\sin\Phi\,\partial_{qq}\,F(0),\\\\[7.11317pt]
\displaystyle c=\beta\,\sin\Phi\,\partial_{qqq}\,F(0).\end{array}$ (23)

This can be seen as a nonlinear analog of the Courant-Snyder invariant that
includes higher-order terms and can be easily propagated through the linear
part of the lattice, thus defined for any azimuth, $\mathcal{K}[p,q;s]$.

Figure 8:  Detuning for the thin sextupole lens as a function of the betatron
tune $\nu_{0}$. The green curve corresponds to the value (SX-2) that matches
numerical simulation, while the orange curve represents the first-order
approximation (SX-1).

If necessary, $(q,p)$ can be inverted back to $(x,\dot{x})$, providing results
in terms of physical variables. Using thin (Th) sextupole and thin octupole
lenses as an example Lee (2018):

$F_{\mathrm{sxt}}(x)=k_{\text{sxt}}\,x^{2},\qquad F_{\mathrm{oct}}(x)=\pm
k_{\text{oct}}\,x^{3},$

we obtain first and second order detuning values:

$\begin{array}[]{l}\displaystyle\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(Th)}}_{\text{SX-1}}}{\mathrm{d}J^{\text{(Th)}}_{\text{SX-1}}}\right|_{J=0}=-\frac{1}{16\,\pi}\,\,\frac{9\,\cos(\pi\,\nu_{0})+\cos(3\,\pi\,\nu_{0})}{\sin(3\,\pi\,\nu_{0})}\,\,\,\beta^{3}\,k_{\text{sxt}}^{2},\\\\[15.6491pt]
\displaystyle\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(Th)}}_{\text{SX-2}}}{\mathrm{d}J^{\text{(Th)}}_{\text{SX-2}}}\right|_{J=0}=-\frac{3}{16\,\pi}\left[3\,\cot(\pi\,\nu_{0})+\cot(3\,\pi\,\nu_{0})\right]\,\beta^{3}\,k_{\text{sxt}}^{2},\\\\[15.6491pt]
\displaystyle\displaystyle\left.\frac{\mathrm{d}\nu^{\text{(Th)}}_{\text{DO,FO}}}{\mathrm{d}J^{\text{(Th)}}_{\text{DO,FO}}}\right|_{J=0}=\mp\frac{3}{8\,\pi}\,\beta^{2}\,k_{\text{oct}},\end{array}$

see Fig. 8. These results are consistent with other derivations using various
perturbation theories including the Deprit perturbation theory Michelotti
(1995) and the Lie algebra treatment Bengtsson (1997); Morozov and Levichev
(2017). Notice that in addition to the scaling provided by Eqs. (23), an
additional factor equal to the Jacobian of the transformation (22),
$\mathbf{J}=\beta\,\sin\Phi$, must be taken into account to obtain the
equations above from the McMillan-Hénon detunings (20) and (21).

###  VI.2 Large amplitudes

While the dependencies related to small amplitudes, such as the rotation
number or detuning around the fixed point, are smooth functions of the map
parameter, the distance to the largest stable trajectory (dynamical aperture)
in the case of general chaotic dynamics forms a fractal curve. Mathematically,
the boundary of the area of stability in the parameter space is defined by
lines associated with $n$-cycles of the system, and is known to be a non-
trivial task since the early works attributed to P. Fatou, G. Julia and B.
Mandelbrot Fatou (1917a, b); Julia (1918); Brooks and Matelski (1981);
Mandelbrot (1980).

Another nontrivial dependence in the case of large amplitudes is the behavior
of the rotation number. Considering $\nu(q_{0})$ for a fixed value of the map
parameters is an example of a singular (continuous, but not absolutely
continuous) function similar to a Cantor function. When plotted against the
initial condition $q_{0}$, the graph reveals a “devil’s staircase” resulting
from mode-locking (chaotic islands) that occurs for every rational $\nu$.

The applicability of any perturbation theory in a particular order greatly
depends on the “leading” nonlinearity. As previously noted, the area of
stability exhibits a complex shape, yet McMillan multipoles offer a reliable
estimate for the separatrix $q_{\mathrm{sep}}$ and $\nu(q_{0})$ near low-order
(integer, half-integer, and third-integer) resonances. While we do note a
substantial alignment between our perturbation theory and other methods in
addressing dynamics around the fixed point, it’s important to highlight that
the dependence of $\nu(q_{0})$ on large amplitudes through elliptic functions
inherently differs from the typical power series of $q_{0}$ often obtained in
methods such as Lie algebra. In particular, analytical expressions for the
rotation number of McMillan multipoles experience very rapid change around the
limiting $n$-cycle, providing a more realistic description of behavior near
the bounding separatrix.

To illustrate the concept, we first employ the Hénon cubic map
$f^{\text{(H)}}_{\text{oct}}(q)=a\,q+q^{3}$ in the defocusing octupole regime
($a=1.8$) above the integer resonance

$\delta\nu=\nu_{0}-0\approx 0.072,$

and in the Duffing regime ($a=-2.1$) with an unstable orbit at the origin and
a stable 2-cycle slightly below the half-integer resonance

$\delta\nu=\nu_{0}-\frac{1}{2}\approx-0.07.$

Both case studies are depicted in Fig. 9, showcasing tracking for stable
trajectories in the top row, the corresponding approximated invariant of the
second order in the middle row, and a comparison of rotation numbers along the
second symmetry line assessed from tracking versus perturbation theory in the
bottom. In both situations, we observe a good estimate for the location of the
limiting separatrix (or figure-eight trajectory), as well as qualitative and
quantitative agreement of rotation numbers from numerical experiments and
analytical predictions. When considering dynamics outside of the lemniscate
(plots b.), perturbation theory fails to predict the next isolating invariant
caused by a high-order resonance, while providing an accurate estimate for
$\nu(q_{0})$ in the appropriate range.

Figure 9:  The top row illustrates phase space diagrams for the Hénon cubic
map $f(q)=a\,q+q^{3}$ obtained through tracking. The middle row displays level
sets for the corresponding approximated McMillan-Hénon invariant of the second
order. Dashed and solid green curves are the first ($p=q$) and second
($p=f(q)/2$) symmetry lines, respectively. The bottom row presents a
comparison of the rotation number as a function of the initial coordinate
along the second symmetry line $\nu(q_{0})$, evaluated from tracking (black
curve) and the analytical approximation (shown in color).  Figure 10:  Same as
Fig. 9, but for the Hénon quadratic map $f(q)=a\,q+q^{2}$. An additional row
(SX-1) and an orange curve in the bottom plots represent the first order of
perturbation theory.

Next, we examine the Hénon quadratic map (19), once again above the integer
resonance ($a=1.6$)

$\delta\nu=\nu_{0}-0\approx 0.1,$

and then above the third-integer resonance ($a=-1.2$)

$\delta\nu=\nu_{0}-\frac{1}{3}\approx 0.02,$

as shown in Fig. 10. This time, the two rows in the middle show the first and
second order approximated invariants, including corresponding dependencies
$\nu(q_{0})$ in the plot at the bottom. For the third-order resonance (plots
d.), the approximation SX-1 fails to predict the proper sign of detuning, as
expected, but in the next order SX-2 provides a quite accurate estimate of
$\nu$ within the range of stable trajectories obtained by tracking. Despite
the discrepancy, both orders provide useful information regarding the general
shape and orientation of the phase space trajectories, which is valuable in
practical applications such as resonant beam extraction, as demonstrated in
the example in the next section.

##  VII Extensions and Generalizations

Numerous avenues exist for further generalizing the obtained results and their
application. While each aspect merits its own dedicated publication, here we
provide a brief overview of the possibilities.

$\bullet$ Higher dimensions. The symmetric McMillan map can be extended to an
axially symmetric map with two degrees of freedom (4D phase space), which is
integrable and allows for the separation of variables in polar coordinates.
This map holds potential for implementation in accelerator physics, as it can
be realized by incorporating a short electron lens into existing linear
lattice structures.

This problem is addressed in the second part of this manuscript, while here,
we provide a summary of two key points. Firstly, the McMillan octupole serves
as an exact solution to radial motion with zero angular momentum $p_{\theta}$,
offering insights into general dynamics. Secondly, of particular significance
is the observation that, similar to perturbation theory for one degree of
freedom, the axially symmetric McMillan electron lens can be regarded as a
second-order approximation for a generic axially symmetric thin lens.

$\bullet$ Basis functions. In constructing the perturbation theory, the
integral of motion was sought in the form of a polynomial. As previously
mentioned, only in the first two orders (in addition to the zeroth order) do
we obtain an approximated invariant that is an exact invariant for another
system (symmetric McMillan map) sharing the same map form. This observation
stems from Suris’ theorem Suris (1989), which establishes that for mappings in
the form (1), integrable systems with smooth invariants are limited to
regular, exponential, and trigonometric polynomials of degree two. This
limitation not only excludes polynomials of degree higher than biquadratic in
$q$ and $p$, but also severely restricts the class of functions where
integrable approximations can be sought.

On the other hand, it should be noted that if trigonometric polynomials are
employed instead of regular ones, integrable approximations via Suris mappings
are still achieved in the first two orders. This alternative approach is more
suitable for systems akin to the Chirikov map Chirikov (1969, 1979), which
feature periodic force functions $f(q)$ with specific regularity properties
Zolkin  _et al._ (2023b).

$\bullet$ Higher orders. Another question that remains beyond the scope of
this article is how to proceed with higher orders of perturbation theory. In
an upcoming manuscript, we will describe the application of Danilov’s theorem
to extract action-angle variables from an approximated invariant, and, an
averaging procedure that facilitates the minimization of (16) in higher orders
by selecting $C\neq 0$ in (17).

$\bullet$ Form of the map. Throughout this article, we have utilized a model
form of the map given by Eq. (1). While this may initially seem restrictive,
we have demonstrated that it encompasses numerous well-known dynamical
systems, ranging from the integrable symmetric McMillan map to chaotic systems
like Hénon and Chirikov, as well as an accelerator lattice with a single thin
nonlinear lens. Furthermore, according to Turaev theorem Turaev (2002), almost
every symplectic map of the plane (and even higher-dimensional) can be
approximated by iterations of mappings in this form.

However, what approach should be adopted for a general symplectic
transformation that represents more realistic systems, such as an accelerator
with multiple nonlinear lenses (not necessarily thin) located at different
positions along the lattice?

In such cases, one approach is to decompose the mapping into the involution of
two consecutive mappings (1), and, then seek an approximated invariant of
motion corresponding to the asymmetric McMillan map. Alternatively, one can
employ physical variables and search for a general polynomial. It is worth
noting that “thick” transformations or lenses can always be incorporated into
the analysis. If an exact transformation is known for each element of the
lattice, the full one-turn map can be expanded into a power series up to the
third order. If some transformations are unknown, a symplectic integrator
Yoshida (1990) consisting of drifts and thin kicks, or the Dirac interaction
picture along with the Magnus expansion Morozov and Levichev (2017), can be
used.

In the first scenario, the dynamics of the asymmetric map can be described
with the help of parametrizations by Iatrou and Roberts and Danilov theorem to
provide canonical variables, similar to what was demonstrated for the McMillan
sextupole. To illustrate the second approach, let’s consider the case of
third-integer resonant extraction for the Mu2e experiment at Fermilab’s muon
campus. The delivery ring consists of six sextupole magnets separated by
linear lattice inserts with a combined tune set up to $\nu_{0}=0.35$, slightly
above the third-integer resonance

$\delta\nu=\nu_{0}-\frac{1}{3}\approx 0.017$

for this specific example. After collecting a one-turn map (that includes all
elements) and expanding it into a power series, we apply perturbation theory
to dynamics in the horizontal plane $(x,\dot{x})$, since the beam is flat.
Seeking approximated invariants from zeroth to the second order, we obtain:

$\begin{array}[]{l}\mathcal{K}^{(0)}(x,\dot{x})=7.2998\,\dot{x}^{2}-2.1672\,\dot{x}\,x+0.2504x^{2},\\\\[7.11317pt]
\mathcal{K}^{(1)}(x,\dot{x})=\mathcal{K}^{(0)}(x,\dot{x})-7.2429\,\dot{x}^{3}+5.8544\,\dot{x}^{2}x\,-\\\\[7.11317pt]
\qquad\qquad\qquad-\,1.0091\,\dot{x}\,x^{2}+0.0339\,x^{3},\\\\[7.11317pt]
\mathcal{K}^{(2)}(x,\dot{x})=\mathcal{K}^{(1)}(x,\dot{x})-0.034\,\dot{x}^{3}x+0.037\,\dot{x}^{2}x^{2}\,-\\\\[7.11317pt]
\qquad\qquad\qquad-\,0.028\,\dot{x}x^{3}+1.5\times 10^{-5}\,x^{4},\end{array}$

where $K^{(0)}(x,\dot{x})$ is proportional to the conventional Courant-Snyder
invariant $\mathrm{C.S.}$. Figure 11 provides corresponding level sets
alongside a comparison with tracking (dots in color) for the second order
SX-2.

Figure 11:  Third-integer resonant extraction for the Fermilab Mu2e
experiment. The top plots present level sets for the zeroth (linear) and
first-order approximated invariants (SX-1), while the plot at the bottom
displays the second-order approximation (SX-2) along with tracking (colored
dots). All plots are presented in physical coordinates $(x,\dot{x})$.

While the conventional Courant-Snyder formalism has been successfully utilized
by the accelerator community for years, it’s important to acknowledge its
limitations as it solely focuses on linearization, disregarding nonlinear
effects and leaving room for improvement (as evidenced by the disparity
between large amplitude tracking with zeroth order $n=0$ invariant). Of
particular interest are the linearly independent amplitude linear tune
$\nu_{0}$ and the estimation of dynamical aperture. Upon comparing tracking
results with higher orders, we observe that key features such as the shape and
orientation of stable phase space trajectories align quite well, mirroring our
observations with the standard form of the map.

The suggested extension is highly natural, given that the original Courant-
Snyder and approximated invariants both follow polynomial forms, allowing for
seamless integration with virtually any existing tracking software. This
extension could prove invaluable for rapid (non-tracking) phase space analysis
and facilitating the incorporation of chromatic effects from sextupoles and
tune shifts due to octupoles. Moreover, if asymmetric/symmetric McMillan
mappings serve as the basis, the newly extended invariant aligns once more
with an exactly integrable system familiar to the accelerator community.

##  VIII Summary

In this article, we revisited an integrable McMillan map and demonstrated its
central role in the general symplectic dynamics of the plane. Through
perturbation theory, we showed that McMillan sextupole and octupole mappings
serve as the first and second-order approximations for planar transformations
in standard form (1). This framework allows for a natural extension of the
optical function formalism employed in accelerator physics, encompassing
chromatic effects and tune shifts induced by sextupole and octupole magnets,
along with dynamical aperture around low-order resonances. As an illustration,
we applied these concepts to a real accelerator lattice utilized for the
third-order resonant extraction at Mu2e experiment (Fermilab), revealing a
strong correspondence between tracking results and analytical predictions of
perturbation theory.

This manuscript marks the first in a series of publications dedicated to the
McMillan map and related integrable systems in higher dimensions, with
applications to accelerator physics. Here, we offer a systematic exposition of
fundamental concepts, providing a comprehensive description of stable
trajectories, including bifurcation diagrams, the parametrization of invariant
curves (see Iatrou and Roberts (2002)), sets of canonical action-angle
variables, global expression for Poincaré rotation numbers (nonlinear betatron
tune) and detuning evaluated at the origin.

##  IX Acknowledgments

The authors would like to thank Eric Stern (FNAL) and Taylor Nchako
(Northwestern University) for carefully reading this manuscript and for their
helpful comments. Moreover, we would like to extend our gratitude to Vladimir
Nagaslaev (FNAL) for multiple discussions and his generous contributions in
preparation of Fig. 11.

This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics. This work was also supported by
Brookhaven Science Associates, LLC under Contract No. DESC0012704 with the
U.S. Department of Energy.

##  Appendix A  Symmetries of invariants of motion

In this section, we discuss the symmetries of the invariants of motion with
respect to system parameters, beginning with the McMillan sextupole

$\mathcal{K}_{\mathrm{sxt}}[p,q]=p^{2}q+p\,q^{2}+\Gamma\,(p^{2}+q^{2})+2\,\epsilon\,p\,q.$

The following propositions are straightforward to verify.

###### Proposition A.0.1

Simultaneous change of map parameters ($\varepsilon\neq 0$)

$\Gamma\rightarrow\varepsilon\,\Gamma,\qquad\qquad\qquad\qquad\quad\epsilon\rightarrow\varepsilon\,\epsilon,$

along with the scaling transformations

$(p,q)\rightarrow\varepsilon\,(p,q)\qquad\text{and}\qquad\mathcal{K}\rightarrow\varepsilon^{3}\,\mathcal{K}$

leaves the form of the map and biquadratic invariant.

###### Proposition A.0.2

Simultaneous change of map parameters

$\Gamma\rightarrow\frac{\Gamma-2\,\epsilon}{3},\qquad\qquad\epsilon\rightarrow-\frac{4\,\Gamma+\epsilon}{3},$

along with the translation transformations

$(p,q)\rightarrow(p,q)+\frac{2}{3}\,(\Gamma+\epsilon)\quad\mathrm{and}\quad\mathcal{K}\rightarrow\mathcal{K}-\left[\frac{2}{3}\,(\Gamma+\epsilon)\right]^{3}$

leaves the form of the map and biquadratic invariant.

There are two consequences stemming from the first proposition. Initially, in
the parameter space of the map $(\Gamma,\epsilon)$, all dynamical systems on a
ray $\epsilon/\Gamma=\mathrm{const}$ that starts at the origin are similar up
to a scaling transformation. This is the only intrinsic parameter
corresponding to the trace of the Jacobian evaluated at the fixed point at the
origin, denoted as $\zeta^{(1-1)}=(0,0)$:

$\mathrm{Tr}\,\mathrm{J}_{\mathrm{sxt}}(\zeta^{(1-1)})=-\frac{2\,\epsilon}{\Gamma}=a.$

Further, Proposition A.0.1 tells us that if we know the dynamics on a ray
$\epsilon/\Gamma=\mathrm{const}$, the dynamics for the systems on the opposite
ray, with $\Gamma\rightarrow-\Gamma$ and $\epsilon\rightarrow-\epsilon$, is
given by the rotation of the phase space by an angle of $\pi$ and inversion of
$\mathcal{K}\rightarrow-\mathcal{K}$, (case $\varepsilon=-1$).

The second proposition arises from two potential choices of the origin: the
translation

$(p,q)\rightarrow(p,q)-\zeta^{(1-2)}$

moves the origin to the second fixed point, denoted as

$\zeta^{(1-2)}=-\frac{2}{3}\,(1,1)\,(\Gamma+\epsilon).$

The parameter transformation,

$\begin{bmatrix}\Gamma\\\\[5.69046pt]
\epsilon\end{bmatrix}\rightarrow\begin{bmatrix}\frac{1}{3}&-\frac{2}{3}\\\\[5.69046pt]
-\frac{4}{3}&-\frac{1}{3}\end{bmatrix}\begin{bmatrix}\Gamma\\\\[5.69046pt]
\epsilon\end{bmatrix},$

is an inversion ($\det=-1$) with two eigenvectors: the eigenvector along the
line $\epsilon=-\Gamma$ with a unit eigenvalue

$v_{1}=(-1,1)/\sqrt{2},\qquad\lambda_{1}=1,$

and another eigenvector with a minus unit egenvalue along the line
$\epsilon=2\,\Gamma$,

$v_{2}=(1,2)/\sqrt{5},\qquad\lambda_{2}=-1.$

For all dynamical systems along the line $\epsilon=-\Gamma$, fixed points
$\zeta^{(1-1,2)}$ undergo transcritical bifurcation and degenerate to a single
point. When $\epsilon=2\,\Gamma$, the phase space portrait displays mirror
symmetry with respect to the line $p=-q-1/2$, with points $\zeta^{(1-1,2)}$
being mirror images of each other. Without loss of generality, and based on
Propositions A.0.1–A.0.2, we will focus our analysis solely on a specific
region of the parameter space (see dashed sector A in the left plot of the
Fig. 12)

$(\epsilon>-\Gamma)\,\cap\,(\epsilon<2\,\Gamma)\qquad\mathrm{or}\qquad
a\in[-4;2],\quad\Gamma\geq 0.$

Figure 12:  Left plot illustrates the $(\Gamma,\epsilon)$-plane of parameters
for the McMillan sextupole map. The dynamics in sectors A and A′ (or B and B′)
are interrelated through Proposition A.0.1, while the dynamics in sectors A
and B (or A′ and B′) are connected through Proposition A.0.2. For the McMillan
octupole (plot on the right), we observe focusing ($\Gamma>0$) and defocusing
($\Gamma<0$) regimes, where the dynamics in sectors A and A’ or B and B’ are
linked by Proposition A.0.5. The filled regions denote the nonredundant
parameter space.

Next, for the invariant of the McMillan octupole,

$\mathcal{K}_{\text{oct}}[p,q]=p^{2}q^{2}+\Gamma\,(p^{2}+q^{2})+2\,\epsilon\,p\,q,$

we have the following propositions:

###### Proposition A.0.3

Simultaneous change of map parameters ($\varepsilon>0$)

$\Gamma\rightarrow\varepsilon\,\Gamma,\qquad\qquad\qquad\qquad\quad\epsilon\rightarrow\varepsilon\,\epsilon,$

along with the scaling transformations

$(p,q)\rightarrow\sqrt{\varepsilon}\,(p,q)\qquad\mathrm{and}\qquad\mathcal{K}\rightarrow\varepsilon^{2}\,\mathcal{K}$

leaves the form of the map and biquadratic invariant.

###### Proposition A.0.4

Change of the map parameter

$\epsilon\rightarrow-\epsilon,$

along with the reflection with respect to $p=0$ (or $q=0$),

$\begin{array}[]{l}q^{\prime}=q\\\\[7.11317pt]
p^{\prime}=-p\end{array}\qquad\left(\text{or}\qquad\begin{array}[]{l}q^{\prime}=-q\\\\[7.11317pt]
p^{\prime}=p\end{array}\right)$ (24)

leaves the biquadratic invariant.

As a consequence of Proposition A.0.3, similar to the sextupole case, the
system still has one intrinsic parameter ($a=-2\,\epsilon/\Gamma$), and the
dynamics on the ray $\epsilon/\Gamma=\mathrm{const}$ with the initial point at
the origin is identical up to the scaling transformation. However, since
$\varepsilon$ is strictly positive, in contrast to the previous situation, we
have two possible forms of the invariant corresponding to regimes with
focusing and defocusing motion, for $\Gamma\gtrless 0$. The next Proposition
A.0.4 is less powerful, as it only guarantees the conservation of the form of
the invariant but not the map. In fact, it can be accompanied by another
proposition:

###### Proposition A.0.5

Change of the map parameter

$\epsilon\rightarrow-\epsilon,$

influences the dynamics of stable trajectory on a given energy level
$\mathcal{K}(p,q)\rightarrow\mathcal{K}(p,-q)=\mathcal{K}(-p,q)$, as follows:

  * •

Symmetric fixed points transform to 2-cycle, and vice versa 2-cycle transforms
to a pair of symmetric fixed points.

  * •

Dynamics around fixed point at the origin transforms according to

$\begin{array}[]{l}J\rightarrow J,\\\\[7.11317pt]
\nu\rightarrow\frac{1}{2}-\nu,\end{array}$

where $J$ is the action and $\nu$ is the rotation number.

$\bigtriangleup$ Proof. The transformation (24) only alters the orientation
without affecting the area under the closed curve, implying $J\rightarrow J$.
Due to the symmetries of the invariant, $\mathcal{K}(p,q)=\mathcal{K}(q,p)$
and $\mathcal{K}(p,q)=\mathcal{K}(-q,-p)$, we can represent the action as
$J=2\,(J_{\mathrm{I}}+J_{\mathrm{II}})$, where $J_{\mathrm{I}}$ denotes the
area over $2\pi$ in the I quadrant, and $J_{\mathrm{II}}$ represents the area
in the II quadrant of the $(q,p)$-plane. Utilizing Danilov’s theorem Zolkin
_et al._ (2017); Nagaitsev and Zolkin (2020), we observe that

$\nu=\frac{\mathrm{d}J_{\mathrm{II}}}{\mathrm{d}J}\qquad\rightarrow\qquad\frac{\mathrm{d}J_{\mathrm{I}}}{\mathrm{d}J}=\frac{\mathrm{d}\left(\frac{J}{2}-J_{\mathrm{II}}\right)}{\mathrm{d}J}=\frac{1}{2}-\nu$

as $J_{\mathrm{I}}$ interchanges $J_{\mathrm{II}}$ under the action of (24).
$\bigtriangledown$

Consequently, we can focus solely on cases where $a>0$ for both the focusing
($\Gamma>0$) and defocusing ($\Gamma<0$) regimes. Subsequently, with the aid
of the propositions, we can deduce the dynamics for $a<0$.

##  Appendix B Stability of fixed points and 2-cycles.
Bifurcation diagrams.

Both fixed points ($n=1$) and $n$-cycles are fundamental concepts that play
significant roles in understanding the behavior of a dynamical system. For a
mapping $\mathcal{M}$, they are respectively defined by

$\zeta^{(n)}\equiv(q,p):\qquad\qquad\mathcal{M}^{n}(q,p)=(q,p)$

and are considered stable if the trace of the Jacobian is smaller than two by
its absolute value

$\left|\mathrm{Tr}\,\mathrm{J}_{\mathcal{M}^{n}}\left[\zeta^{(n)}\right]\right|<2.$

Iatrou and Roberts Iatrou and Roberts (2002) proved that for an integrable map
with a smooth integral, isolated critical points of the integral belong to
(isolated) cycles of the map, and the points of isolated cycles of the map are
(isolated) critical points of the integral. Applying this to a general
asymmetric McMillan map, one can verify that the invariant of motion can have
up to 5 isolated critical points. Without loss of generality, we assume that
one of them corresponds to a fixed point at the origin. The four other
critical points have to be fixed points as well, or two of them can form an
isolated 2-cycle. In an exceptional case, these fixed points and 2-cycles can
appear on the same level of the invariant and then form 3- or 4-cycles; in
this scenario, the map degenerates to a linear and, in fact, all orbits
(except the fixed point at the origin) become periodic with corresponding
periods of 3 and 4.

Furthermore, for the mapping in the form (1), fixed points belong to the main
diagonal and are given by the intersection of two symmetry lines

$p=q\qquad\mathrm{and}\qquad p=f(q)/2.$

2-cycles are given by the intersection of the second symmetry line with its
inverse

$p=f(q)/2\qquad\mathrm{and}\qquad q=f(p)/2.$

If $f(q)$ is an odd function, then fixed points have to appear in symmetric
pairs with respect to the origin, and 2-cycle belong to the intersection of
the anti-diagonal $p=-q$ with $p=f(q)/2$.

Below, we summarize the stability of fixed points and $n$-cycles for McMillan
multipoles and provide associated bifurcation diagrams of the invariant.

###  B.1 McMillan sextupole

The sextupole map exhibits two fixed points defined for any values of
parameters on the $(\Gamma,\epsilon)$-plane: the point at the origin,

$\zeta^{(1-1)}=(0,0),$

and one on the main diagonal,

$\zeta^{(1-2)}=-\frac{2}{3}\,(\Gamma+\epsilon)\,(1,1)=\frac{a-a_{0}}{3}\,(1,1)\,\Gamma.$

Additionally, there exists a 2-cycle defined by

$\begin{array}[]{l}\zeta^{(2)}=(\Gamma-\epsilon)\,(1,1)\pm\sqrt{(\Gamma-\epsilon)(5\,\Gamma-\epsilon)}\,(1,-1)\\\\[7.11317pt]
\displaystyle=\frac{1}{2}\,\left[(a-a_{\frac{1}{2}})\,(1,1)\pm\sqrt{(a-a_{\frac{1}{2}})(a-a^{*}_{\frac{1}{2}})}\,(1,-1)\right]\,\Gamma,\end{array}$

which is restricted to a sector where
$(\Gamma-\epsilon)(5\,\Gamma-\epsilon)>0$, as shown in the bottom plot of the
right column in Fig. 13.

Figure 13:  The left plot illustrates the domain of stability for fixed points
and the 2-cycle in the $(\Gamma,\epsilon)$-plane of parameters for the
McMillan sextupole. The corresponding real domain $\zeta\in\mathbb{R}$ for
each case is depicted in the right column. If the fixed point/2-cycle is
stable or real, then the corresponding area is filled with color.

Stability in the space of parameters $(\Gamma,\epsilon)$ and bifurcation
diagram for fixed points/2-cycle as a function of $a$ are presented in the
main plot of Fig. 13 and left plot of Fig. 15 respectively. Analysis reveals
that the fixed point at the origin is stable for $a_{1/2}=-2<a<a_{0}=2$, where
the lower index for a particular value of $a_{\nu}$ represents the
corresponding rotation number $\nu=\arccos(a_{\nu}/2)/(2\pi)$. When $a=a_{0}$,
$\zeta^{(1-1)}$ undergoes a transcritical bifurcation [T], exchanging
stability with $\zeta^{(1-2)}$. Both points go through a sub-critical period
doubling bifurcation at $a=a_{1/2}$ [SBPD1] or $a^{*}_{1/2}=-10$ [SBPD2] for
$\zeta^{(1-1)}$ and $\zeta^{(1-2)}$ respectively. Finally, for $a=a_{1/3}=-1$,
the second fixed point and 2-cycle approach the same energy level, causing the
map to degenerate [D] to linear with all orbits around the origin being period
3. The 2-cycle is always unstable when defined on a real domain but becomes
stable if $\zeta^{(2)}\in\mathbb{C}^{2}$.

Towards the conclusion of this section, Fig. 16 offers an additional diagram
providing insight into all potential dynamical regimes. It showcases phase
space plots featuring different level sets of the invariant
$\mathcal{K}_{\mathrm{sxt}}[p,q]$. These plots are organized along rays (or
between them) representing constant values of $a=-2\,\epsilon/\Gamma$,
indicative of the system’s degeneracy or bifurcations. The specific values of
$a$ are highlighted on a scale delineated with a dash-dotted circle.

###  B.2 McMillan octupole

Similar to the sextupole scenario, the octupole map also possesses a fixed
point at the origin:

$\zeta^{(1-1)}=(0,0),$

which is defined for all $(\Gamma,\epsilon)\in\mathbb{R}^{2}$ and remains
stable for $|\epsilon|<|\Gamma|$ (or $|a|<2$). Additional symmetric fixed
points and a 2-cycle are provided by:

$\zeta^{(1-2,3)}=\pm\sqrt{-(\Gamma+\epsilon)}\,(1,1)=\pm\sqrt{\pm\frac{a-a_{0}}{2}}\,(1,1)\sqrt{|\Gamma|},$
(25)
$\zeta^{(2)}=\pm\sqrt{-(\Gamma-\epsilon)}\,(1,-1)=\pm\sqrt{\mp\frac{a-a_{\frac{1}{2}}}{2}}\,(1,-1)\sqrt{|\Gamma|},$
(26)

where the $\pm$-sign under the square roots varies based on whether $\Gamma$
is greater or less than 0; note that, in contrast to the sextupole map, the
“natural” unit of distance shifts to $\sqrt{|\Gamma|}$ rather than $\Gamma$.
These points are real only for $\epsilon<-\Gamma$ (or $a\gtrless a_{0}$ for
$\Gamma\gtrless 0$) and $\epsilon>\Gamma$ (or $a\lessgtr a_{1/2}$ for
$\Gamma\gtrless 0$), as visualized in the right column of Fig. 14.

Figure 14:  Same as Fig. 13 but for McMillan octupole.

The primary plot in Fig. 14 alongside the subsequent plots in Fig. 15
illustrates the stability and bifurcation diagrams across all conceivable
scenarios. Now, let’s break down the results: (I) When $\Gamma<0$, fixed
points $\zeta^{(1-2,3)}$ and the 2-cycle are always unstable, and the only
region with stable trajectories is $a_{1/2}<a<a_{0}$. At $a_{1/2,0}$, we
observe that the fixed point at the origin undergoes sub-critical period
doubling [PD] or pitchfork [PF] bifurcations, respectively. (II) For
$\Gamma>0$, at $a_{1/2,0}$ the fixed point $\zeta^{(1-1)}$ goes through super-
critical [PD] or [PF] bifurcations. When it is stable ($|a|<2$), it is the
only isolated cycle with all trajectories being finite. (III) In the Duffing
regime, the point at the origin is locally unstable, becoming a center of
lemniscate separatrix. If $a>a_{0}$, trajectories inside the figure-eight
curve round the centers of nearby “eyes,” while when $a<a_{1/2}$, the orbit
jumps from “eye” to “eye” since their centers form a 2-cycle. (IV) If $a=0$,
it degenerates into a linear system with a map period of 4. Fig. 17
consolidates the aforementioned results, depicting phase space diagrams on the
plane of map parameters $(\Gamma,\epsilon)$, similar to Fig. 16.

Figure 15:  Bifurcation diagrams illustrating fixed points and 2-cycles of the
McMillan sextupole (SX), octupole in defocusing (DO), focusing (FO), and
Duffing (Df, $|a|>2$) regimes. Plots depict normalized coordinates as a
function of the map parameter $a=-2\,\epsilon/\Gamma$. Stable fixed
points/2-cycles are represented by solid lines, while unstable ones are
dashed. Transcritical, pitchfork, and period doubling bifurcations are denoted
by [T], [PF], and [PD] respectively; SB indicates sub-critical, and SP
represents super-critical. Label [D] corresponds to system degeneracy.  Figure
16:  Schematic phase space diagrams for the McMillan sextupole map arranged in
the plane of map parameters $(\Gamma,\epsilon)$. Isolated fixed points,
2-cycles and degenerate 3-cycles, along with their corresponding level sets,
are color-coded according to the legend; other level sets are represented in
black. The first ($p=q$), the second ($p=f(q)/2$) and the inverse ($q=f(p)/2$)
symmetry lines are displayed in thick black, green, and cyan respectively.
Figure 17:  Schematic phase space diagrams for the McMillan octupole map
arranged in the plane of map parameters $(\Gamma,\epsilon)$. Isolated fixed
points, 2-cycles and degenerate 4-cycles, along with their corresponding level
sets, are color-coded according to the legend; other level sets are
represented in black. The first ($p=q$), the second ($p=f(q)/2$) and the
inverse ($q=f(p)/2$) symmetry lines are displayed in thick black, green, and
cyan respectively.

##  Appendix C Action integrals and power series

To derive action integrals for the McMillan octupole, we observe that the
normalized phase space trajectories (10) correspond to elliptic Lissajous
figures

$\begin{array}[]{l}x=\mathrm{ef}\,(\omega_{1}\,t,\kappa),\\\\[5.69046pt]
y=\mathrm{ef}\,(\omega_{2}\,t+\eta,\kappa).\end{array}$

where the Jacobi elliptic functions $\mathrm{ef}$ are employed instead of
regular trigonometric functions for conventionally defined Lissajous curves
and a unit ratio of $\omega_{1}:\omega_{2}$. The area enclosed by such a curve
can be determined using the integral $S=\oint y\,\mathrm{d}x$, which for
“major” elliptic

functions yields Zolkin  _et al._ (2022):

$\begin{array}[]{l}\displaystyle
S_{\text{sn}}=\frac{4}{\kappa^{2}}\frac{\mathrm{sd}^{2}\eta\,\,\mathrm{E}[\kappa]-\mathrm{K}[\kappa]+\mathrm{cn}^{2}\eta\,\,\Pi[k^{2}\mathrm{sn}^{2}\eta,\kappa]}{\mathrm{sn}\,\eta\,\,\mathrm{sd}^{2}\eta},\\\\[8.5359pt]
\displaystyle
S_{\text{cn}}=\frac{4}{\kappa^{2}}\frac{\mathrm{K}[\kappa]-\mathrm{sn}^{2}\eta\,\,\mathrm{E}[\kappa]-\mathrm{cn}^{2}\eta\,\,\Pi[k^{2}\mathrm{sn}^{2}\eta,\kappa]}{\mathrm{sd}\,\eta\,\,\mathrm{sn}^{2}\eta},\\\\[7.11317pt]
\displaystyle
S_{\text{dn}}=\frac{2}{\mathrm{sc}\,\eta\,\,\mathrm{sn}^{2}\eta}\,\times\\\\[9.95863pt]
\left\\{(\mathrm{sn}^{2}\eta+\mathrm{dn}^{2}\eta)\,\mathrm{K}[\kappa]-\mathrm{sn}^{2}\eta\,\,\mathrm{E}[\kappa]-\mathrm{dn}^{2}\eta\,\Pi[k^{2}\mathrm{sn}^{2}\eta,\kappa]\right\\}.\end{array}$

Figure 18 illustrates the behavior of the area as a function of $\eta$ and
$\kappa$ along with several examples of Lissajous-Jacobi orbits. Utilizing the
scaling outlined in (11):

$\bullet$ For $\mathrm{sn}$ and $\mathrm{cn}$ cases, the action integrals and
their power series expressed in terms of elliptic modulus $\kappa$ are:

$\begin{array}[]{rl}J_{\text{sn}}[\kappa]=&\displaystyle\\!\frac{2}{\pi}\,\frac{\mathrm{sd}^{2}\eta_{+}\,\mathrm{E}[\kappa]-\mathrm{K}[\kappa]+\mathrm{cn}^{2}\eta_{+}\,\Pi[k^{2}\,\mathrm{sn}^{2}\eta_{+},\kappa]}{\mathrm{sd}\,\eta_{+}\,\,\mathrm{nd}\,\eta_{+}}\\\\[9.95863pt]
=&\displaystyle\\!\frac{(4-a^{2})^{\textstyle\frac{3}{2}}}{2^{4}}\left[\kappa^{2}+\frac{2-7\,a^{2}}{{2^{4}}}\,\kappa^{4}+\frac{6-68\,a^{2}+33\,a^{4}}{2^{7}}\,\kappa^{6}+\frac{100-2398\,a^{2}+2838\,a^{4}-715\,a^{6}}{2^{12}}\,\kappa^{8}+\mathcal{O}(\kappa^{10})\right],\\\\[11.38092pt]
J_{\text{cn}}[\kappa]=&\displaystyle\\!\frac{2}{\pi}\,\frac{\mathrm{K}[\kappa]-\mathrm{sn}^{2}\eta_{+}\,\mathrm{E}[\kappa]-\mathrm{cn}^{2}\eta_{+}\,\Pi[k^{2}\,\mathrm{sn}^{2}\eta_{+},\kappa]}{\mathrm{sn}\,\eta_{+}\,\,\mathrm{dn}\,\eta_{+}}\\\\[9.95863pt]
=&\displaystyle\\!\frac{(4-a^{2})^{\textstyle\frac{3}{2}}}{2^{4}}\left[\kappa^{2}+7\,\frac{2+a^{2}}{2^{4}}\,\kappa^{4}+\frac{102+44\,a^{2}+33\,a^{4}}{2^{7}}\,\kappa^{6}+\frac{3036+1246\,a^{2}+330\,a^{4}+715\,a^{6}}{2^{12}}\,\kappa^{8}+\mathcal{O}(\kappa^{10})\right],\end{array}$

or as a series of the invariant

$J_{\text{sn,cn}}[\mathcal{K}]=\frac{1}{\sqrt{4-a^{2}}}\,\left[\mathcal{K}\pm\frac{2+a^{2}}{(4-a^{2})^{2}}\,\mathcal{K}^{2}+2\,\frac{6+12\,a^{2}+a^{4}}{(4-a^{2})^{4}}\,\mathcal{K}^{3}\pm
5\,\frac{20+90\,a^{2}+30\,a^{4}+a^{6}}{(4-a^{2})^{6}}\,\mathcal{K}^{4}+\mathcal{O}(\mathcal{K}^{5})\right].$

Furthermore, by expanding the rotation number in the same series up to
$\mathcal{O}(\kappa^{10})$ or $\mathcal{O}(\mathcal{K}^{5})$

$\begin{array}[]{rl}\nu_{\text{sn}}\approx&\displaystyle\\!\nu_{0}-\frac{a\,\sqrt{4-a^{2}}}{2^{4}\,2\,\pi}\left[3\,\kappa^{2}+\frac{86-35\,a^{2}}{2^{5}}\,\kappa^{4}+7\,\frac{46-48\,a^{2}+11\,a^{4}}{2^{7}}\,\kappa^{6}+\frac{39340-74926\,a^{2}+40722\,a^{4}-6435\,a^{6}}{2^{14}}\,\kappa^{8}\right],\\\\[9.95863pt]
\nu_{\text{cn}}\approx&\displaystyle\\!\nu_{0}+\frac{a\,\sqrt{4-a^{2}}}{2^{4}\,2\,\pi}\left[3\,\kappa^{2}+\frac{10+35\,a^{2}}{2^{5}}\,\kappa^{4}+\frac{18-56\,a^{2}+77\,a^{4}}{2^{7}}\,\kappa^{6}+\frac{1364-13\,(26+858\,a^{2}-495\,a^{4})\,a^{2}}{2^{14}}\,\kappa^{8}\right],\end{array}$

$\begin{array}[]{ll}\displaystyle\text{or
}\quad\nu_{\text{sn,cn}}[\mathcal{K}]=&\displaystyle\\!\nu_{0}+\frac{a}{2\,\pi\,(4-a^{2})^{3/2}}\,\left[\mp\,3\,\mathcal{K}-\frac{1}{2}\,\frac{86+13\,a^{2}}{(4-a^{2})^{2}}\,\mathcal{K}^{2}\mp
2\,\frac{322+200\,a^{2}+9\,a^{4}}{(4-a^{2})^{4}}\,\mathcal{K}^{3}\,-\right.\\\\[15.6491pt]
&\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\,\,\,\,-\left.\frac{39340+55506\,a^{2}+10386\,a^{4}+221\,a^{6}}{(4-a^{2})^{6}}\,\mathcal{K}^{4}+\mathcal{O}(\mathcal{K}^{5})\right],\end{array}$

we derive series $\nu(J)$:

$\displaystyle\nu_{\text{sn,cn}}[J]=\nu_{0}+\frac{a}{2\,\pi}\,\left[\mp\frac{3}{4-a^{2}}\,J-\frac{74+7\,a^{2}}{2\,(4-a^{2})^{5/2}}\,J^{2}\mp
5\,\frac{(2+a^{2})(46+a^{2})}{(4-a^{2})^{4}}\,J^{3}+\mathcal{O}(J^{4})\right].$

$\bullet$ For $\mathrm{dn}$ trajectories that orbit symmetric fixed points in
the Duffing regime, case (Df), we have:

$\begin{array}[]{rl}J_{\text{dn}}[\kappa]=&\displaystyle\\!\frac{1}{\pi}\,\frac{(\mathrm{sn}^{2}\eta_{+}+\mathrm{dn}^{2}\eta_{+})\,\mathrm{K}[\kappa]-\mathrm{sn}^{2}\eta_{+}\,\,\mathrm{E}[\kappa]-\mathrm{dn}^{2}\eta_{+}\,\,\Pi[k^{2}\,\mathrm{sn}^{2}\eta_{+},\kappa]}{\mathrm{sn}\,\eta_{+}\,\,\mathrm{cn}\,\eta_{+}}\\\\[9.95863pt]
=&\displaystyle\\!\frac{(a-2)^{3/2}}{2^{4}\sqrt{2}\,a}\left[\kappa^{4}+\kappa^{6}+\left(\frac{119}{8}-\frac{4+3\,a}{a^{2}}\right)\,\frac{\kappa^{8}}{2^{4}}+\left(\frac{55}{8}-\frac{4+3\,a}{a^{2}}\right)\,\frac{\kappa^{10}}{2^{3}}+\mathcal{O}(\kappa^{12})\right],\\\\[12.80365pt]
J_{\text{dn}}[\widetilde{\mathcal{K}}]=&\displaystyle\\!\frac{2^{-3/2}}{\sqrt{a-2}}\,\left[\widetilde{\mathcal{K}}+\frac{8-a}{2^{4}\,(a-2)^{2}}\,\widetilde{\mathcal{K}}^{2}+\frac{192-32\,a+3\,a^{2}}{2^{8}\,(a-2)^{4}}\,\widetilde{\mathcal{K}}^{3}+\frac{12800-5\,(480-72\,a+5\,a^{2})\,a}{2^{13}\,(a-2)^{6}}\,\widetilde{\mathcal{K}}^{4}+\mathcal{O}(\widetilde{\mathcal{K}}^{5})\right].\end{array}$

and

$\begin{array}[]{rl}\nu_{\text{dn}}\approx&\displaystyle\\!\nu_{0}-\frac{1}{2\,\pi}\frac{\sqrt{a-2}}{2^{4}\sqrt{2}\,a^{2}}\left[(4+a)\,(\kappa^{4}+\kappa^{6})-\frac{384+96\,a-964\,a^{2}-235\,a^{3}}{2^{8}\,a^{2}}\,\kappa^{8}-\frac{384+96\,a-452\,a^{2}-107\,a^{3}}{2^{7}\,a^{2}}\,\kappa^{10}\right]\\\\[9.95863pt]
\approx&\displaystyle\\!\nu_{0}-\frac{1}{2\,\pi}\frac{2^{-3/2}}{a\,(a-2)^{3/2}}\left[(4+a)\,\widetilde{\mathcal{K}}-\frac{128-224\,a-68\,a^{2}+5\,a^{3}}{2^{5}\,a\,(a-2)^{2}}\,\widetilde{\mathcal{K}}^{2}\,+\right.\\\\[9.95863pt]
&\displaystyle\\!\left.\qquad\qquad\qquad\qquad\qquad\qquad\quad\,\,\,\,+\,\frac{4096-11264\,a+12672\,a^{2}+4384\,a^{3}-476\,a^{4}+33\,a^{5}}{3\times
2^{8}\,a^{2}\,(a-2)^{4}}\,\widetilde{\mathcal{K}}^{3}\right]\\\\[9.95863pt]
\approx&\displaystyle\\!\nu_{0}-\frac{1}{2\,\pi}\left[\frac{1}{a}\,\frac{4+a}{a-2}\,J-\frac{1}{a^{2}}\frac{128-160\,a-60\,a^{2}+3\,a^{3}}{2^{7/2}(a-2)^{5/2}}\,J^{2}\right],\qquad\quad\,\text{where
}\widetilde{\mathcal{K}}=\mathcal{K}+\left(\frac{a}{2}-1\right)^{2}.\end{array}$

Figure 18:  Sample elliptic Lissajous curves with a unit ratio of
$\omega_{1}:\omega_{2}=1$ (top plots) and the corresponding area
$S_{\text{ef}}$ enclosed under it (plots in the bottom) for major elliptic
functions as a function of modulus $\kappa$ and phase difference $\eta$. Note
that $\eta$ is measured in units of a quarter period $\mathrm{K}[\kappa]$.

$\bullet$ While the roots $q_{1,2,3,4}$ of the characteristic polynomial
$\mathcal{P}(q)$ for the sextupole map represent turning points of the
coordinate, they are not particularly useful for obtaining expansions of $J$
or $\nu$. (Indeed, this can be achieved using modern software.) Therefore, we
solely present an exact expression for the action integral here:

$\begin{array}[]{l}\displaystyle
J_{\mathrm{sxt}}=\frac{1}{2\pi}\,\int_{q_{2}}^{q_{3}}\frac{\sqrt{\mathcal{P}(q)}}{q+1}\,\mathrm{d}q=\frac{q_{2}-q_{1}}{\pi}\,\sqrt{\frac{q_{3}-q_{1}}{q_{4}-q_{2}}}\times\displaystyle\left\\{\frac{a+4+q_{4}-q_{1}}{4}\,\mathrm{K}[\kappa]+\frac{a+4}{4}\frac{q_{4}-q_{2}}{q_{2}-q_{1}}\,\mathrm{E}[\kappa]\,-\right.\\\\[12.80365pt]
\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.-\,\frac{a+1}{q_{3}-q_{1}}\,\Pi\left[\frac{q_{4}-q_{2}}{q_{4}-q_{1}}\,\kappa^{2},\kappa\right]-\frac{(1+q_{3})(1+q_{4})}{q_{3}-q_{1}}\,\Pi\left[\frac{1+q_{1}}{1+q_{2}}\frac{q_{4}-q_{2}}{q_{4}-q_{1}}\,\kappa^{2},\kappa\right]\right\\}.\end{array}$

Figure 19:  Phase space portraits for the McMillan sextupole map with
$a=-3/2$, illustrating level sets of the invariant in the original $(p,q)$ and
new rotated coordinates $(P,Q)$. Turning points are identified by the
intersection of the level set (depicted in black) with either the second
(green curve) or first (dashed green line) symmetry lines, indicated by the
red and blue dots respectively.

A more practical approach involves performing a rotation of the phase space by
an angle $\alpha=\pi/4$ clockwise:

$q\rightarrow Q\,\cos\alpha-P\,\sin\alpha,\qquad\qquad\qquad p\rightarrow
Q\,\sin\alpha+P\,\cos\alpha,$

which allows us to “decouple” one of the roots, as illustrated in Fig. 19.
After implementing this change of variables, the Hamiltonian becomes:

$H[P,Q;t]\quad\rightarrow\quad
K[P,Q;t]=\frac{P^{2}\,(a+2-\sqrt{2}\,Q)-Q^{2}\,(a-2-\sqrt{2}\,P)}{2}.$

Following the same logic as in Section III.1, we again solve for momentum

$P=\sqrt{\frac{\mathcal{G}(Q)}{Q_{4}-Q}},$

which is now expressed through a third-order polynomial $\mathcal{G}$ and a
known root $Q_{4}$ that is independent of $K$:

$\mathcal{G}(Q)=\sqrt{2}\,\mathcal{K}-\frac{2-a}{\sqrt{2}}\,Q^{2}-Q^{3},\qquad\qquad\qquad\qquad\qquad\qquad
Q_{4}=\frac{2+a}{\sqrt{2}}.$

For stable trajectories, we have $Q_{1}<Q_{2}<Q<Q_{3}<Q_{4}$, such that:

$\begin{array}[]{l}\displaystyle
Q_{1}=-\frac{2-a}{3\sqrt{2}}\left(1+2\,\cos\left[\frac{1}{3}\,\arccos\left[1-\frac{54\,K}{(2-a)^{3}}\right]\right]\right),\\\\[12.80365pt]
\displaystyle
Q_{2}=-\frac{2-a}{3\sqrt{2}}\left(1-2\,\sin\left[\frac{1}{3}\,\arcsin\left[1-\frac{54\,K}{(2-a)^{3}}\right]\right]\right),\\\\[12.80365pt]
\displaystyle
Q_{3}=-\frac{2-a}{3\sqrt{2}}\left(1-2\,\cos\left[\frac{1}{3}\,\arccos\left[\frac{54\,K}{(2-a)^{3}}-1\right]\right]\right),\end{array}$

allowing expansion in terms of $K$. Starting with the action integral:

$\begin{array}[]{ll}\displaystyle
J_{\mathrm{sxt}}&\displaystyle\\!\\!=\frac{1}{2\,\pi}\,\oint
P\,\mathrm{d}Q=\frac{1}{\pi}\,\int_{Q_{2}}^{Q_{3}}P\,\mathrm{d}Q=\frac{2}{\pi}\frac{1}{\sqrt{(Q_{4}-Q_{2})(Q_{3}-Q_{1})}}\,\left\\{(1+a)(Q_{4}-Q_{3})\,\Pi\left[\frac{Q_{3}-Q_{2}}{Q_{4}-Q_{2}},\kappa\right]\right.+\\\\[14.22636pt]
&\displaystyle\\!\\!\qquad\qquad+\left.\frac{(2+a)^{2}-K}{4}\,\left[\left(\frac{4+a}{Q_{4}-Q_{1}}-3\sqrt{2}\right)\,\mathrm{K}[\kappa]+\frac{(4+a)(Q_{3}-Q_{1})}{(Q_{4}-Q_{1})(Q_{4}-Q_{3})}\,\mathrm{E}[\kappa]\right]\right\\}\\\\[14.22636pt]
&\displaystyle\\!\\!=\frac{(2-a)(4-a^{2})^{3/2}}{2^{8}}\,\left[\kappa^{4}+\kappa^{6}+\frac{97-27\,a-15\,a-a^{3}}{2^{7}}\,\kappa^{8}+\frac{33-27\,a-15\,a-a^{3}}{2^{6}}\,\kappa^{10}+\mathcal{O}(\kappa^{12})\right]\\\\[14.22636pt]
&\displaystyle\\!\\!=\frac{1}{\sqrt{4-a^{2}}}\,\left[K-3\,\frac{2+2\,a+a^{2}}{(a-2)(a^{2}-4)^{2}}\,K^{2}+10\,\frac{22+48\,a+42\,a^{2}+16\,a^{3}+3\,a^{4}}{(a-2)^{2}(a^{2}-4)^{4}}\,K^{3}+\mathcal{O}(K^{4})\right],\end{array}$

where the new elliptic modulus is given by:

$\kappa=\sqrt{\frac{(Q_{3}-Q_{2})(Q_{4}-Q_{1})}{(Q_{3}-Q_{1})(Q_{4}-Q_{2})}}.$

Next, substitution of $P$ into the equation of motion

$\mathrm{d}t=\left(\frac{\partial K}{\partial
P}\right)^{-1}\,\mathrm{d}q=\mathrm{d}q/\sqrt{2\,(Q_{4}-Q)(Q_{3}-Q)(Q-Q_{2})(Q-Q_{1})}$

provides the period of motion and time of one step of the map:

$\mathrm{T}=\oint\mathrm{d}t=\frac{2\,\sqrt{2}\,\mathrm{K}[\kappa]}{\sqrt{(Q_{4}-Q_{2})(Q_{3}-Q_{1})}},\qquad\qquad\mathrm{T}^{\prime}=\int_{Q_{0}}^{Q_{0}^{\prime}}\mathrm{d}t=\frac{\sqrt{2}\,\mathrm{F}[\arcsin\phi,\kappa]}{\sqrt{(Q_{4}-Q_{2})(Q_{3}-Q_{1})}}.$

Using new turning points as the initial point for the integral,
$Q_{0}=Q_{2,3}$:

$\left\\{Q_{0},P_{0}\right\\}=\left\\{\sqrt{2}\,q_{0},0\right\\}\qquad\rightarrow\qquad\left\\{Q_{0}^{\prime},P_{0}^{\prime}\right\\}=\left\\{\frac{a-q_{0}}{1+q_{0}}\frac{q_{0}}{\sqrt{2}},\frac{a-2-3\,q_{0}}{1+q_{0}}\frac{q_{0}}{\sqrt{2}}\right\\}$

defines

$\phi=\sqrt{Q_{2}\,\frac{Q_{3}-Q_{1}}{Q_{3}-Q_{2}}\frac{2\,\sqrt{2}-\sqrt{2}\,a+3\,Q_{2}}{2\,Q_{1}\,(\sqrt{2}+Q_{2})-Q_{2}\,(\sqrt{2}\,a-Q_{2})}}=\sqrt{Q_{3}\,\frac{(2\,\sqrt{2}+\sqrt{2}\,a-2\,Q_{2})(2\,\sqrt{2}-\sqrt{2}\,a+3\,Q_{3})}{2\,(Q_{3}-Q_{2})(4+2\,a+2\,\sqrt{2}\,Q_{3}+Q_{3}^{2})}},$

and finally, the rotation number:

$\begin{array}[]{rl}\nu_{\mathrm{sxt}}=&\displaystyle\\!\\!\frac{\mathrm{T}^{\prime}}{\mathrm{T}}=\frac{\mathrm{F}[\arcsin\phi,\kappa]}{2\,\mathrm{K}[\kappa]}\\\\[14.22636pt]
=&\displaystyle\\!\\!\nu_{0}-\frac{1+a}{2\,\pi}\,\frac{\sqrt{4-a^{2}}}{2^{8}}\left[(8+a)\,(\kappa^{4}+\kappa^{6}-\kappa^{10})+5\,\frac{768+30\,a-78\,a^{2}-17\,a^{3}-a^{4}}{2^{9}}\,(\kappa^{8}+2\,\kappa^{10})+\mathcal{O}(\kappa^{12})\right]\\\\[14.22636pt]
=&\displaystyle\\!\\!\nu_{0}-\frac{1}{2\,\pi}\frac{1+a}{(2-a)(4-a^{2})^{3/2}}\left[(8+a)\,K+\frac{832+734\,a+258\,a^{2}+13\,a^{3}-a^{4}}{2\,(2-a)(4-a^{2})^{2}}\,K^{2}+\mathcal{O}(K^{3})\right]\\\\[14.22636pt]
=&\displaystyle\\!\\!\nu_{0}-\frac{1}{2\,\pi}\frac{1+a}{(2-a)(4-a^{2})}\left[(8+a)\,J+\frac{736+626\,a+198\,a^{2}+7\,a^{3}-a^{4}}{2\,(2-a)(4-a^{2})^{3/2}}\,J^{2}+\mathcal{O}(J^{3})\right].\end{array}$

## References

  * Alexeev (1965) V. Alexeev, Bull. Inst. Theoret. Astron 10, 241 (1965).
  * Marchal (1966) C. Marchal, Bulletin Astronomique 1, 189 (1966).
  * Marchal (1986) C. Marchal, Celestial mechanics 38, 377 (1986).
  * Jakas (1995) M. M. Jakas, Physical Review A 52, 866 (1995).
  * Jakas (1996) M. M. Jakas, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 115, 255 (1996).
  * Strand and Reinhardta) (1979) M. P. Strand and W. P. Reinhardta), The Journal of Chemical Physics 70, 3812 (1979).
  * McMillan (1971) E. M. McMillan, in  _Topics in modern physics, a tribute to E.V. Condon_ , edited by E. Brittin and H. Odabasi (Colorado Associated University Press, Boulder, CO, 1971) pp. 219–244.
  * Suris (1989) Y. B. Suris, Functional Analysis and Its Applications 23, 74 (1989).
  * Iatrou and Roberts (2002) A. Iatrou and J. A. Roberts, Nonlinearity 15, 459 (2002).
  * Cathey  _et al._ (2021) B. Cathey, G. Stancari, A. Valishev, and T. Zolkin, Journal of Instrumentation 16, P03041 (2021).
  * Quispel  _et al._ (1988) G. Quispel, J. Roberts, and C. Thompson, Physics Letters A 126, 419 (1988).
  * Quispel  _et al._ (1989) G. Quispel, J. Roberts, and C. Thompson, Physica D: Nonlinear Phenomena 34, 183 (1989).
  * Hénon and Heiles (1964) M. Hénon and C. Heiles, The Astronomical Journal 69, 73 (1964).
  * Duffing (1918) G. Duffing,  _Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedeutung_ , 41-42 (Vieweg, 1918).
  * Zolkin  _et al._ (2022) T. Zolkin, S. Nagaitsev, and I. Morozov, “Mcmillan map and nonlinear twiss parameters,” (2022), arXiv:2204.12691 [nlin.SI] .
  * Henon (1969) M. Henon, Quarterly of applied mathematics 27, 291 (1969).
  * Chirikov (1969) B. V. Chirikov, Nuclear Physics, Novosibirsk.(Engl. Trans. CERN-Trans-71-40, 1971) (1969), translated at CERN by A.T. Sanders from the Russian (IYAF-267-TRANS-E).
  * Chirikov (1979) B. V. Chirikov, Physics Reports 52, 263 (1979).
  * Zolkin  _et al._ (2017) T. Zolkin, S. Nagaitsev, and V. Danilov, arXiv preprint arXiv:1704.03077 (2017).
  * Nagaitsev and Zolkin (2020) S. Nagaitsev and T. Zolkin, Phys. Rev. Accel. Beams 23, 054001 (2020).
  * Mitchell  _et al._ (2021) C. E. Mitchell, R. D. Ryne, K. Hwang, S. Nagaitsev, and T. Zolkin, Phys. Rev. E 103, 062216 (2021).
  * Kolmogorov (1954) A. N. Kolmogorov (1954) pp. 527–530.
  * Moser (1962) J. Moser, Matematika 6, 51 (1962).
  * Arnol’d (1963) V. Arnol’d, Russian Mathematical Surveys 18, 9 (1963).
  * zol (2017) _Proc. of North American Particle Accelerator Conference (NAPAC’16), Chicago, IL, USA, October 9-14, 2016_ , North American Particle Accelerator Conference No. 3 (JACoW, Geneva, Switzerland, 2017) presented at NAPAC2016 in Chicago, IL, USA, unpublished.
  * DeVogelaere  _et al._ (1958) R. DeVogelaere, H. ANTOSIEWICZ, W. T. KYNER, R. BASS, S. LEFSCHETZ, D. BUSHAW, L. MARKUS, R. DeVOGELAERE, P. MENDELSON, S. KAKUTANI, G. SEIFERT, and D. SLOTNICK, “On the structure of symmetric periodic solutions of conservative systems, with applications,” in  _Contributions to the Theory of Nonlinear Oscillations (AM-41), Volume IV_ (Princeton University Press, 1958) pp. 53–84.
  * Lewis Jr (1961) D. C. Lewis Jr, Pacific Journal of Mathematics 11, 1077 (1961).
  * Zolkin  _et al._ (2023a) T. Zolkin, Y. Kharkov, and S. Nagaitsev, Phys. Rev. Res. 5, 043241 (2023a).
  * Veselov (1991) A. P. Veselov, Russian Mathematical Surveys 46, 1 (1991).
  * Brown (1993) M. Brown, Lecture notes in pure and applied mathematics , 83 (1993).
  * Cairns  _et al._ (2016) G. Cairns, Y. Nikolayevsky, and G. Rossiter, The American Mathematical Monthly 123, 363 (2016).
  * Zolkin  _et al._ (2023b) T. Zolkin, Y. Kharkov, and S. Nagaitsev, pre-print (2023b), arXiv:2311.17616 [nlin.SI] .
  * Lee (2018) S.-y. Lee,  _Accelerator Physics (Fourth Edition)_ (World Scientific Publishing Company, 2018).
  * Courant and Snyder (1958) E. D. Courant and H. S. Snyder, Annals of physics 3, 1 (1958).
  * Michelotti (1995) L. Michelotti,  _Intermediate Classical Dynamics with Applications to Beam Physics_ (Wiley, 1995).
  * Bengtsson (1997) J. Bengtsson,  _The sextupole scheme for the Swiss Light Source (SLS): an analytic approach_ , Tech. Rep. SLS-TME-TA-1997-0009 (Paul Scherrer Institut, 1997).
  * Morozov and Levichev (2017) I. Morozov and E. Levichev, CERN Proceedings 1, 195 (2017).
  * Fatou (1917a) P. Fatou, Comptes rendus de l’Académie des sciences, Paris 164, 806 (1917a).
  * Fatou (1917b) P. Fatou, Comptes rendus de l’Académie des sciences, Paris 165, 992 (1917b).
  * Julia (1918) G. Julia, Journal de Mathématiques Pures et Appliquées 8, 47 (1918).
  * Brooks and Matelski (1981) R. Brooks and J. P. Matelski, “The dynamics of 2-generator subgroups of psl(2,c),” in  _Riemann Surfaces And Related Topics: Proceedings of the 1978 Stony Brook Conference. (AM-97), Volume 97_ , edited by I. Kra and B. Maskit (Princeton University Press, Princeton, 1981) pp. 65–72.
  * Mandelbrot (1980) B. B. Mandelbrot, Annals of the New York Academy of Sciences 357, 249 (1980).
  * Turaev (2002) D. Turaev, Nonlinearity 16, 123 (2002).
  * Yoshida (1990) H. Yoshida, Physics letters A 150, 262 (1990).