††thanks: Currently at AWS Quantum Technologies

# Integrable symplectic maps with a polygon tessellation

T. Zolkin<EMAIL_ADDRESS>Fermilab, PO Box 500, Batavia, IL 60510-5011  Y.
Kharkov  Joint Center for Quantum Information and Computer Science,
NIST/University of Maryland, College Park, MD 20742  Joint Quantum Institute,
NIST/University of Maryland, College Park, MD 20742  S. Nagaitsev
<EMAIL_ADDRESS>Jefferson Lab, Newport News, VA 23606  Old Dominion
University, Norfolk, VA 23529

###### Abstract

The identification of integrable dynamics remains a formidable challenge, and
despite centuries of research, only a handful of examples are known to date.
In this article, we explore a special form of area-preserving (symplectic)
mappings derived from the stroboscopic Poincaré cross-section of a kicked
rotator — an oscillator subjected to an external force periodically switched
on in short pulses. The significance of this class of problems extends to
various applications in physics and mathematics, including particle
accelerators, crystallography, and studies of chaos. Notably, Suris’ theorem
constrains the integrability within this category of mappings, outlining
potential scenarios with analytic invariants of motion. In this paper, we
challenge the assumption of the analyticity of the invariant, by exploring
piecewise linear transformations on a torus ($\mathbb{T}^{2}$) and associated
systems on the plane ($\mathbb{R}^{2}$), incorporating arithmetic
quasiperiodicity and discontinuities. By introducing a new automated
technique, we discovered previously unknown scenarios featuring polygonal
invariants that form perfect tessellations and, moreover, fibrations of the
plane/torus. In this way, this work reveals a novel category of planar tilings
characterized by discrete symmetries that emerge from the invertibility of
transformations and are intrinsically linked to the presence of integrability.
Our algorithm relies on the analysis of the Poincaré rotation number and its
piecewise monotonic nature for integrable cases, contrasting with the noisy
behavior in the case of chaos, thereby allowing for clear separation. Some of
the newly discovered systems exhibit the peculiar behavior of “integrable
diffusion,” marked by infinite and quasi-random hopping between tiles while
being confined to a set of invariant segments. Finally, through the
implementation of a smoothening procedure, all mappings can be generalized to
quasi-integrable scenarios with smooth invariant motion, thereby opening doors
to potential practical applications. Interactive simulations featuring
visualizations of some of the discovered maps are accessible at
http://integrable-maps.streamlit.app.

###### Contents

  1. I Introduction
  2. II  Definitions and notations
    1. II.1 McMillan-Hénon form of the map
    2. II.2 Force function with arithmetical quasiperiodicity
  3. III  Dynamics with arithmetic quasiperiodicity
  4. IV  Chaos detection
    1. IV.1 Search Algorithm
  5. V  Results
    1. V.1 Map $\mathcal{M}_{\mathrm{G}}^{6}$
    2. V.2 Mappings $\mathcal{M}_{\alpha 3}^{4}$ and $\mathcal{M}_{\mathrm{F}1^{\prime}}^{4}$
    3. V.3 Mappings $\mathcal{M}_{\mathrm{E}1^{\prime}}^{3}$, $\mathcal{M}_{\mathrm{a}2}^{6}$ and $\mathcal{M}_{\mathrm{b}2}^{3}$
  6. VI  Integrable polygon maps on tor
    1. VI.1  Unwrapping and integrable diffusion
  7. VII  Applications of mappings with periodic and arithmetically quasiperiodic forces.
    1. VII.1  Conection to McMillan-Suris mappings
    2. VII.2  Smoothening procedure
  8. VIII  Summary
  9. IX Acknowledgments

##  I Introduction

Integrable systems, regarded as fundamental pillars in our understanding of
complex dynamics, manifest across diverse disciplines such as celestial
mechanics, hydrodynamics, nonlinear optics, and condensed matter physics.
Within the domain of classical mechanics, archetypal integrable Hamiltonians
include well-recognized problems, such as the harmonic oscillator, Kepler
problem, Toda lattice, and the Euler top. Considered mathematical “gems” in
the realm of physics, integrable systems are exceedingly rare, with only a
handful of textbook examples known to date.

In the study of Hamiltonian mechanics, symplectic mappings, such as Poincaré
or stroboscopic cross-sections derived from higher-dimensional continuous
systems, serve as essential bridges between discrete and continuous dynamics.
On one hand, they serve as insightful tools for comprehending geometric nature
of motion, and on the other, they can be applied as numerical integrators for
systems governed by Hamilton’s equations of motion. Symplectic mappings are
pervasive in physical applications and naturally emerge in various contexts,
including cyclic particle accelerators Meiss (1992); McMillan (1971); Danilov
(2008), Lagrange descriptions of fluid dynamics and billiards, to name a few
(e.g., see Veselov (1991)).

Among all area-preserving transformations of the plane, our focus is on a
specific representation known as the McMillan-Hénon (MH) form of the map where
dynamics is governed by a force function, $f(q)$ (see McMillan (1971); Turaev
(2002)). This form arises from the stroboscopic Poincaré cross-section (time
discretization) of a kicked rotator, a crucial toy model for exploring
classical and quantum chaos in non-autonomous Hamiltonian systems. It has a
very clear analytical structure, and despite its seemingly simple appearance,
incorporates an exceptionally rich set of scenarios, including (but not
limited) to chaotic Hénon Henon (1969) and Chirikov Chirikov (1969, 1979), and
integrable McMillan mappings McMillan (1971). Of particular significance is
Suris’ theorem Suris (1989), which illuminates potential integrable scenarios
within this category. Suris demonstrated that if an invariant of motion
$\mathcal{K}(p,q)$ takes the form of an analytic function, its possibilities
are confined to being regular, exponential, or trigonometric polynomials of
degree two. This theorem serves as a robust constraint, shaping the landscape
of plausible integrable symplectic maps.

By challenging the assumption of analyticity of the invariant of motion, one
opens the possibility of discovering novel types of integrable symplectic
mappings that go beyond the constraints outlined in Suris’ theorem. Notable
examples of systems exhibiting linear dynamics, yet possessing a nontrivial
phase space structure, include the Brown-Knuth map Brown (1983); Brown and
Slifker (1985); Brown (1993), McMillan beheaded and two-headed ellipses
McMillan (1971), and CNR maps Cairns  _et al._ (2014, 2016). In all instances,
the force is a piecewise linear function, while the integral of motion adopts
a piecewise nature, characterized by segments encompassing hyperbolas,
ellipses, and/or line segments. When all pieces are line segments, the phase
space is fibrated with polygons or polygonal chains, for stable and unstable
cases respectively. In our earlier study Zolkin  _et al._ (2022), we drew
attention to the existence of nonlinear systems with polygonal structure,
which led to the development of an automated search algorithm and the
subsequent discovery of a diverse set of new integrable maps.

Here, we present the next phase of our exploration, building upon the
foundation laid by previous results. In this work, we examine piecewise linear
forces comprised of an infinite number of segments, adhering to a specific
condition of arithmetic quasiperiodicity. In addition, we explore periodic
piecewise linear forces that allow for discontinuities, as well as a more
general class of systems constrained on a torus, which are related to planar
dynamics through a wrapping/unwrapping procedure. These systems can be
considered as a generalization of the famous Bernoulli map Schuster and Just
(2006) (also known as the dyadic transformation or bit-shift map) within
planar symplectic dynamics. Moreover, not only did we uncover new large
families of previously unknown integrable systems but we also succeeded in
once again automating the search process. This time, we employed a novel
algorithm rooted in identifying the rotation number ($\nu$) and utilizing
analytical insights into the piecewise monotonic dependence of $\nu$ on
amplitude. We would like to note that automated searches for integrable
systems are still in their early stages. While some intriguing approaches,
including machine learning, have been recently developed (see Krippendorf  _et
al._ (2021)), they are currently only capable of re-identifying systems that
are already known.

What is intriguing is that, beyond the significance of our results in the
arenas of integrable dynamics and algorithm theory, we unexpectedly uncover a
connection with the problem of periodic tessellations by polygons. Such
challenges, like the famous Euclidean plane tilings through convex regular
polygons, have captivated minds since antiquity. The mathematical journey
began with Kepler, who initiated a systematic treatment in his monumental work
“Harmonices Mundi,” and over time, this pursuit extended to encompass
Archimedean, aperiodic, Coxeter, and various other forms of tessellation.
Different types of tilings are characterized by specific rules governing the
permissible polygons and their arrangements, commonly referred to as the
“tessellation rule.”

For all newly discovered systems, there exists a lattice constant $L$
determining the translational symmetry of the invariant
$\mathcal{K}(q,p)=\mathcal{K}(q+m\,L,p+n\,L)$ for integer values of $m$ and
$n$. This symmetry is attributed to the periodic nature of $f(q)$.
Consequently, the constant level sets of $\mathcal{K}(p,q)$ corresponding to
isolating integrals of motion (separatrices) perfectly tessellate the plane
with polygons, while all other sets contribute to a continuous fibration by
polygons within all tiles necessary for integrability. In both instances, the
rules governing tessellation and fibration are intricately tied to the map’s
invertibility and its decomposition into anti-area-preserving reflections
Lewis Jr (1961), as elucidated in Section II. While there is a comprehensive
understanding of certain types of polygonal tessellations, in general, more
sophisticated or aperiodic patterns continue to be areas of active research
Smith  _et al._ (2023a, b); Zeller  _et al._ (2021); Manolas  _et al._ (2022),
and many key aspects have not achieved full clarity. Perhaps, in this work, by
establishing a connection between polygonal tiling and integrability, we are
adding a new, important contribution regarding all possible systems of this
type.

The present paper is organized as follows. We start in Section II by
introducing the McMillan-Hénon form of the map, key symmetries related to
integrability, and the idea of arithmetic quasiperiodicity. In Section III, we
give clear examples illustrating how integrable tessellations work. In Section
IV, we reveal a new method for spotting chaos with easy-to-follow examples and
describe the automated search algorithm. Sections V and VI then unveil the
main discoveries, explaining both plane and torus tessellations. We analyze
possible dynamical scenarios and report on an intriguing form of motion termed
integrable diffusion. Finally, Section VII discusses some of the possible
applications of the map, including a “smoothening” procedure that allows the
removal of all non-physical features, such as lack of differentiability or
discontinuities, and yields quasi-integrable dynamics with suppressed chaotic
behavior.

##  II  Definitions and notations

In this section, we begin by introducing key concepts such as symplecticity,
orbits, and integrability of the map, setting the stage for a clear
understanding of the system’s dynamics.

  1. 1.

A mapping $\mathcal{M}:X\rightarrow X$ on a $2m$-dimensional manifold $X$ is
symplectic if it satisfies the condition

$\mathrm{J}_{\mathcal{M}}^{\mathrm{T}}\,\Omega\,\mathrm{J}_{\mathcal{M}}=\Omega,\qquad\qquad\qquad\mathrm{J}_{\mathcal{M}}=\begin{pmatrix}\partial\mathbf{q}^{\prime}/\partial\mathbf{q}&\partial\mathbf{q}^{\prime}/\partial\mathbf{p}\\\\[5.69046pt]
\partial\mathbf{p}^{\prime}/\partial\mathbf{q}&\partial\mathbf{p}^{\prime}/\partial\mathbf{p}\end{pmatrix},\qquad\qquad\qquad\Omega=\begin{pmatrix}0&\mathrm{I}_{m}\\\\[5.69046pt]
-\mathrm{I}_{m}&0\end{pmatrix},$ (1)

where a pair of $m$-dimensional variables, representing configuration
$\mathbf{q}$ and momenta $\mathbf{p}$, collectively form a phase space
coordinate denoted as $\zeta=(\mathbf{q},\mathbf{p})$, the application of the
map is denoted as $(^{\prime})$

$\zeta^{\prime}=(\mathbf{q}^{\prime},\mathbf{p}^{\prime})=\mathcal{M}(\mathbf{p},\mathbf{q}),$

$\mathrm{J}_{\mathcal{M}}$ is the Jacobian matrix of $\mathcal{M}$, $\Omega$
is a $2m\times 2m$ nonsingular skew-symmetric matrix, and, $\mathrm{I}_{m}$ is
a $m\times m$ identity matrix. In essence, symplecticity guarantees that the
transformation preserves the sum of areas projected on individual planes
$(q_{i},p_{i})$ within the phase space, emphasizing the geometric structure of
the dynamical system. For mappings of the plane ($m=1$), the symplectic
condition (1) implies that $\det\mathrm{J}_{\mathcal{M}}=1$, indicating that
the map preserves both area and orientation.

  2. 2.

The orbit of a map $\mathcal{M}$ is defined by a sequence of points,
$\zeta_{k}$, obtained through the repeated application of the map to an
initial condition, commonly referred to as the seed, $\zeta_{0}$:

$\left\\{\zeta_{0},\zeta_{1},\zeta_{2},\ldots\right\\},\qquad\qquad\qquad\zeta_{k}=\mathcal{M}^{k}(\zeta_{0}).$

We say that the set of points
$\left\\{\zeta_{0},\zeta_{1},\ldots,\zeta_{n-1}\right\\}$ forms an $n$-cycle
if

$\mathcal{M}^{n}(\zeta_{0})=\zeta_{0}.$

An $n$-cycle is a type of periodic orbit, and the smallest positive integer
$n$ for which the above condition holds is called the period of the cycle. If
$n=1$, the map leaves the point unchanged after a single iteration, and we
refer to this state as a fixed point, indicated by $\zeta_{0}=\zeta_{*}$. If
$n=2$, it is a 2-cycle, and so on.

  3. 3.

A map of the plane is said to be integrable if there exists a non-constant
continuous function $\mathcal{K}(p,q)$, called integral or invariant of
motion, that remains conserved upon the application of the map for all initial
conditions, i.e.,

$\forall\,(q,p)\in\mathbb{R}^{2}:\qquad\mathcal{K}(p^{\prime},q^{\prime})-\mathcal{K}(p,q)=0.$
(2)

This implies that the orbits originating from $\zeta_{0}$ lie on a constant
level of $\mathcal{K}=\mathcal{K}(\zeta_{0})$, whether it manifests as a
singular point, a set of points, or as a single (or collection) of closed/open
curve(s).

###  II.1 McMillan-Hénon form of the map

In this article, we focus on a particular representation of the symplectic
mapping of the plane referred as the McMillan-Hénon (MH) form. A detailed
overview can be found in McMillan (1971); Turaev (2002); Zolkin  _et al._
(2022), while here we briefly list some of its relevant properties. The
transformation and its inverse are explicitly expressed as:

$\begin{array}[]{ll}\mathcal{M}_{\mathrm{MH}}:&\quad
q^{\prime}=p,\\\\[7.11317pt] &\quad
p^{\prime}=-q+f(p),\end{array}\qquad\qquad\qquad\qquad\qquad\begin{array}[]{ll}\mathcal{M_{\mathrm{MH}}}^{-1}:&\quad
q^{\prime}=-p+f(q),\\\\[7.11317pt] &\quad p^{\prime}=q,\end{array}$ (3)

where $f(q)$ is called the force function. This map maintains symplecticity
for any $f(q)$, offering a framework for various well-known examples. These
include the extensively studied Chirikov standard map Chirikov (1969, 1979),
the area-preserving Hénon quadratic map Henon (1969), or a collection of
integrable mappings discovered by E. McMillan McMillan (1971) and Y. Suris
Suris (1989).

The map in MH form can be written as a composition
$\mathcal{M}_{\mathrm{MH}}=\mathcal{F}\circ\mathrm{Rot}(\pi/2)$, where
$\mathrm{Rot}$ is a rotation about the origin followed by a thin lens
transformation $\mathcal{F}$:

$\mathrm{Rot}(\theta):\quad\begin{bmatrix}q^{\prime}\\\
p^{\prime}\end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}\begin{bmatrix}q\\\
p\end{bmatrix}\qquad\qquad\mathrm{and}\qquad\qquad\mathcal{F}:\quad\begin{bmatrix}q^{\prime}\\\
p^{\prime}\end{bmatrix}=\begin{bmatrix}q\\\ p\end{bmatrix}+\begin{bmatrix}0\\\
f(q)\end{bmatrix}.$

This breakdown unveils the physical nature of the system, depicting a linear
oscillator subject to periodic kicks in time, influenced by the force $f(q)$.
In specific applications, a model accelerator featuring one degree of freedom,
comprising a linear optics insert and a thin nonlinear lens (e.g., sextupole
or octupole), can be redefined in the MH form.

Alternatively, $\mathcal{M}_{\mathrm{MH}}$ can be decomposed as
$\mathcal{M}_{\mathrm{MH}}=\mathcal{G}\circ\mathrm{Ref}(\pi/4)$, where
$\mathrm{Ref}(\theta)$ is a reflection about a line passing through the origin
at an angle $\theta$ with the $q$-axis, and $\mathcal{G}$ is a nonlinear
vertical reflection:

$\mathrm{Ref}(\theta):\quad\begin{bmatrix}q^{\prime}\\\
p^{\prime}\end{bmatrix}=\begin{bmatrix}\cos 2\theta&\sin 2\theta\\\ \sin
2\theta&-\cos 2\theta\end{bmatrix}\begin{bmatrix}q\\\
p\end{bmatrix}\qquad\qquad\mathrm{and}\qquad\qquad\mathcal{G}:\quad\begin{bmatrix}q^{\prime}\\\
p^{\prime}\end{bmatrix}=\begin{bmatrix}q\\\
-p\end{bmatrix}+\begin{bmatrix}0\\\ f(q)\end{bmatrix}.$

Both, $\mathrm{Ref}(\theta)$ and $\mathcal{G}$, are anti-area-preserving
involutions, which means that a double application of the map is equivalent to
the identity transformation $\mathrm{Ref}^{2}(\theta)=\mathcal{G}^{2}=I_{2}$
and the determinant of their Jacobians is equal to $-1$. Each reflection
possesses a line of fixed points:

$l_{1}:\,\,p=q\quad\text{for}\quad\mathrm{Ref}(\pi/4),\qquad\qquad\qquad\qquad\qquad
l_{2}:\,\,p=f(q)/2\quad\text{for}\quad\mathcal{G},$

designated as the first and second symmetry lines Lewis Jr (1961); McMillan
(1971).

Any existing invariant level set, even amid chaos and the absence of a global
invariant for all initial conditions, remains unchanged not only under
$\mathcal{M}_{\mathrm{MH}}$, but also when subjected to either one of the
reflections:

$\mathcal{K}(p,q)=\mathcal{K}(q,p)\qquad\qquad\qquad\mathrm{and}\qquad\qquad\qquad\mathcal{K}(p,q)=\mathcal{K}(-p+f(q),q).$
(4)

If set of Eqs. (4) is satisfied for any point in the plane, the map
$\mathcal{M}_{\mathrm{MH}}$ is integrable, providing an easy analytical test
to verify integrability by checking that every invariant set is symmetric with
respect to the main diagonal $p=q$ and consists of parts that are vertically
equidistant from $p=f(q)/2$.

###  II.2 Force function with arithmetical quasiperiodicity

In Zolkin  _et al._ (2022), we introduced numerous nonlinear integrable
mappings with polygon invariants and corresponding piecewise linear force
functions, denoted here as $f_{\mathrm{p.l.}}(q)$. To further expand upon this
result, one can pursue a straightforward approach by exploring forces that
incorporate an increasing number of segments. Alternatively, we can take a
different path to extend our findings by imposing the arithmetical
quasiperiodicity, which implies that for a given $L$

$\forall\,q:\quad
f_{\mathrm{a.q.}}(q+L)=f_{\mathrm{a.q.}}(q)+F\qquad\qquad\mathrm{where}\qquad\qquad
F=f_{\mathrm{a.q.}}(L)-f_{\mathrm{a.q.}}(0)=\mathrm{const}.$ (5)

As evident, a function with arithmetic quasiperiodicity can be expressed as
the sum of a purely periodic function with a period of $L$ and a linear
function with a slope of $k_{0}=F/L$

$\quad\,\,f_{\mathrm{a.q.}}(q)=d+k_{0}\,q+f_{\mathrm{per}}(q)\qquad\qquad\mathrm{where}\qquad\qquad\forall\,q:\,f_{\mathrm{per}}(q+L)=f_{\mathrm{per}}(q).\,\,\,$
(6)

Note, it is crucial to avoid confusing the concept of arithmetic
quasiperiodicity with the idea of quasiperiodic motion in the context of a
symplectic map of the plane. In the latter case, the trajectory of the system
forms a dense curve (or collection of curves) in phase space, and although it
never exactly repeats, it approaches previously visited points arbitrarily
closely. In other words, quasiperiodicity of oscillations implies that the
rotation number is irrational, unlike in periodic motion, where the trajectory
repeats after a fixed number of iterations with a rational rotation number.

Figure 1:  Left plot displays an example of a regular piecewise linear force
function made of 3 segments. Considering the interval $q\in[0,L]$, we can
construct two new forces by imposing arithmetic quasiperiodicity, as shown in
the plot in the middle, or by “unwrapping” it from a torus, the plot on the
right. The solid and dashed curves in cyan represent the periodic and linear
parts of the force with arithmetic quasiperiodicity. For the last plot, the
dashed curve in cyan corresponds to linear flow on the torus, i.e.,
$k_{0}\,(q\,\,\mathrm{mod}\,L)$.

Starting with a regular piecewise linear function that is continuous in the
interval $q\in[0,L]$, after imposing Eq. (5), we obtain a new continuous
function with infinitely many segments and a special “periodic” condition,
refer to the left and middle plots in Fig. 1 for an illustration. Assuming
that one period of $f_{\mathrm{a.q.}}(q)$ has $n$ distinct segments with the
first line segment beginning at $q=0$, for the $j$-th segment we have:

$f_{\mathrm{a.q.}}(q)=k_{j}\,q+\sum_{i=1}^{j-1}k_{i}\,l_{i}+d,\qquad\qquad\qquad
0\leq\sum_{i=1}^{j-1}l_{i}\leq q<\sum_{i=1}^{j}l_{i}\leq
L,\qquad\qquad\qquad(j=1,\ldots,n).$

In this article, we will utilize notations that define the force function as
$\left[\\{k_{1},l_{1}\\},\ldots,\\{k_{n},l_{n}\\};d\right]$, where each of the
$n$ tuples $\\{k_{j},l_{j}\\}$ represents the slope and length of the $j$-th
segment, followed by the value of the vertical shift parameter,
$d=f_{\mathrm{a.q.}}(0)$.

Finally, we will explore two additional generalizations of the
$f_{\mathrm{p.l.}}(q)$ by constraining the dynamics within a torus, where the
equations of motion are subject to a modulo $L$ condition, and then
“unwrapping” it back onto an Euclidean plane as

$f_{\mathrm{tor}}(q)=f_{\mathrm{p.l.}}(q\,\,\mathrm{mod}\,\,L)\quad\mathrm{mod}\,\,L.$
(7)

If $F=0$, we obtain the same outcome as we would for $f_{\mathrm{a.q.}}(q)$,
otherwise, $f_{\mathrm{tor}}(q)$ results in a purely periodic function with
discontinuities, see the plot on the right in Fig. 1. This can be seen as the
inclusion of additional segments with their lengths approaching zero and
slopes tending towards infinity.

The form of the map and the repetitive nature of the force function give rise
to several sources of redundancy in the parameter space. These redundancies
are as follows:

  1. 1.

Each map has a corresponding “twin” map, which can be obtained by using the
transformation $f(q)\rightarrow-f(-q)$. Geometrically, this is equivalent to
rotating the phase space by an angle of $\pi$, i.e. $(q,p)\rightarrow-(q,p)$.

  2. 2.

Cyclic permutation of segments in the force function results in the same map,
provided that the shift parameter $d$ is adjusted to maintain the relative
location of the fixed point, $q_{*}$.

  3. 3.

For any given set of parameters, there exists a set of equidistant values of
the parameter $d$, expressed as $d=d_{0}+m\,(2\,L-F)$, where $m\in\mathbb{Z}$
is an integer. This generates an equivalent dynamical system with the fixed
point $q_{*}\rightarrow q_{*}+m\,L$, owing to the translational symmetry of
the invariant.

  4. 4.

Lastly, due to the linear nature of each segment, the rescaling of dynamical
variables $(q,p)\rightarrow\epsilon\,(q,p)$ leaves the form of the map
invariant. This rescaling allows us to choose one of the segments to have a
length equal to one.

In order to avoid reporting duplicated results, we will eliminate any
redundant solutions and scaling parameters whenever possible during our
analysis. This will ensure that we focus on distinct and meaningful outcomes
in our exploration of the parameter space.

##  III  Dynamics with arithmetic quasiperiodicity

Exploring forces in the form (5), which are obtained from a piecewise linear
functions consisting of two segments, we encountered two integrable maps:

$\mathcal{M}_{\mathrm{a}1}^{6}:\quad
f(q)=\left[\\{0,1\\},\\{2,1\\};2\,m\right]\qquad\qquad\mathrm{and}\qquad\qquad\mathcal{M}_{\mathrm{b}1}^{3}:\quad
f(q)=\left[\\{0,1\\},\\{-2,1\\};4+6\,m\right].$

The corresponding phase space diagrams are presented in plots (a.1) and (b.1)
of Fig. 2. Notably, due to the arithmetic quasiperiodicity of the force, the
constant level sets of the invariant fill the phase space plane with tiles
arranged in a regular pattern and separated by a net representing a
separatrix/isolating integral of motion (shown in red). Adopting terminology
from famous Archimedean tilings, we say that the tiling is regular if the
phase space is tessellated with a single type of polygon separatrices, and
that it is semi-regular when the nets consist of multiple polygon types,
maintaining a coherent pattern across the entire plane.

Figure 2:  Phase space diagrams for maps $\mathcal{M}_{\mathrm{a}1}^{6}$ (a.1)
and $\mathcal{M}_{\mathrm{b}1}^{3}$ (b.1). The black line segment in the
middle of the central cell [cc] consists of fixed points (a.1) or 2-cycles
with a fixed point in the middle (b.2). The blue polygons represent invariant
level sets separated by red separatrices. The solid and dashed green lines in
the diagrams correspond to the second and first symmetry lines, respectively,
$p=f(q)/2$ and $p=q$. Three solid symmetry lines in plot (b.2) illustrate
different configurations of the map. The connected series of line segments
shown in brown depict one period of the force function $f(q)\in[0;L]$ along
with the values for their slopes $k_{i}$ and length $l_{i}$.

While all tiles are geometrically identical, the dynamics on constant level
sets of the invariant varies from one cell of periodicity to another. First
and foremost, there is a unique tile containing a fixed point, which we refer
to as the central cell (or cc for brevity). In both mappings, this cell
contains an invariant line segment in the middle (shown in black), while the
rest of the phase space is filled with invariant polygons (shown in blue)
bounded by the rectangle separatrix (shown in red). Within the central cell,
the dynamics revolve around the fixed point and correspond to the integrable
fragment of the map with a standard (i.e., without arithmetic
quasiperiodicity) piecewise linear force function. For the examples provided,
cc is based on mappings a1 and b1, as indicated by the subscript in the map’s
name (see Zolkin  _et al._ (2022) for further details).

Figure 3:  Schematic illustration of the mode-locking around the central cell
[cc] or central node [cn] for different values of the global rotation number
$\nu_{P}$. Tiles of the same color are grouped in chains of islands, and the
arrows indicate the order in which particular tiles are visited.

Figure 4:  Plot (a.1) displays several iterations under
$\mathcal{M}_{\mathrm{a}1}^{6}$ for two initial conditions (dots with white
circle) separated by one lattice period $L$. Plots (a.1′) and (a.1′′) show the
corresponding invariant level sets. The blue dots/polygon represent the orbit
inside cc, while the cyan dots/polygons indicate a similar orbit starting from
the nearby cell on the right.

Furthermore, all initial conditions on the separatrix isolating cc correspond
to $P$-fold periodic orbits, leading to the non-central tiles grouping into
chains of $P$ islands. Consequently, initial conditions from non-central tiles
exhibit mode locking with a rational rotation number defined by

$\nu_{P}=\frac{1}{2\,\pi}\,\arccos\frac{k_{0}}{2}=\frac{R}{P},\qquad\qquad\qquad\qquad\qquad
k_{0}=\frac{F}{L},$ (8)

where $k_{0}$ is the tangent of the average slope of $f(q)$ determined by Eqs.
(5,6) and $R\in\mathbb{Z}$ is an integer. We refer to $\nu_{P}\in\mathbb{Q}$
and $P\in\mathbb{N}$ as the global rotation number and global period of the
map, respectively; the value of $P$ is provided in the superscript of the
map’s name for convenience. Fig. 3 schematically shows the pattern in which
tiles are grouped for global periods equal to 6, 3, and 4, from left to right,
respectively. A more detailed illustration is given in Fig. 4, providing an
example of two orbits for the map $\mathcal{M}_{\mathrm{a}1}^{6}$ generated
from initial conditions separated by one lattice period in plot (a.1), along
with their corresponding invariant level sets in plots (a.1′) and (a.1′′). The
blue and cyan dots/polygons represent the initial conditions within the
central cell $(q_{0},p_{0})\in\mathrm{cc}$ and $(q_{0}+L,p_{0})$,
respectively.

Additionally, for the map $\mathcal{M}_{\mathrm{b}1}^{3}$, there are two more
ways (configurations) in which the dynamics can be organized:

$f(q)=\left[\\{0,1\\},\\{-2,1\\};6\,m\right]\qquad\qquad\qquad\mathrm{or}\qquad\qquad\qquad
f(q)=\left[\\{0,1\\},\\{-2,1\\};6\,m+2\right].$

In fact, the two systems above are twin maps. Now, instead of being inside the
cc, the fixed point is placed into the node/vertex of the separatrix net. We
will refer to the map with a fixed point in the node as having a central node
(or cn for short). In this case, there is no central cell, and initial
conditions within all tiles exhibit mode-locked behavior. An illustration of
how tiles are grouped into chains for $P=3$ and $4$ is given in the second row
of plot (c.) in Fig. 2.

Finally, to analyze dynamical systems, we will use the rotation number inside
cc as a function of the horizontal distance from the initial condition to the
fixed point, denoted by $s=q_{0}-q_{*}$, along both of the symmetry lines:

$\nu_{1}=\left.\nu(s)\right|_{\zeta_{0}=(q_{0},q_{0})}\qquad\qquad\mathrm{and}\qquad\qquad\nu_{2}=\left.\nu(s)\right|_{\zeta_{0}=(q_{0},f(q_{0})/2)}.$

For example, for the systems mentioned above, we have

$\mathcal{M}_{\mathrm{a}1}^{6}:\quad\nu_{1}=\left\\{\begin{array}[]{ll}\displaystyle
0,&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{2\,(|s|-1/2)}{2+8\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq
1,\end{array}\right.\qquad\qquad\qquad\nu_{2}=\left\\{\begin{array}[]{ll}\displaystyle
0,&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{|s|-1/2}{2+4\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq\frac{3}{2},\end{array}\right.$

and

$\mathcal{M}_{\mathrm{b}1}^{3}:\quad\nu_{1}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{2},&\qquad\displaystyle
s=0,\\\\[9.95863pt]
\displaystyle\frac{1+2\,(|s|-1/2)}{2+8\,(|s|-1/2)},&\qquad\displaystyle
0<|s|\leq\frac{1}{2},\end{array}\right.\qquad\qquad\qquad\nu_{2}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{2},&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{1+(|s|-1/2)}{2+4\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq\frac{3}{2}.\end{array}\right.$

##  IV  Chaos detection

After demonstrating the possibility of constructing nonlinear integrable
mappings with polygon invariants and corresponding $f(q)$ displaying
arithmetic quasiperiodicity, a new question emerged: “How can we find more
systems and automate the search process”? In our previous work Zolkin  _et
al._ (2022), when $f(q)$ was a standard piecewise linear force function, we
relied on the assumption that the phase space should be fibrated with
concentric polygons. By arranging points of an individual orbit by polar angle
and counting the number of vertices, we were able to reject non-integrable
systems where the polygons produced after rearrangement of chaotic
trajectories had an increasing number of vertices that did not converge with
the number of iterations.

Unfortunately, as mentioned in Zolkin  _et al._ (2022), this algorithm fails
when the phase space involves islands, i.e., in the presence of mode-locking.
In such cases, the trajectory sequentially visits each island under iterations
and traces out a collection of invariant tori. Consequently, arranging the
points of the corresponding orbit by polar angle becomes problematic, as
$\mathcal{K}(r,\theta)=\mathrm{const}$ becomes a multiple-valued function of
$\theta$ defined on disjoint intervals, as seen in the example in plot (a.1′′)
of Fig. 4.

To address this issue, we decided to rely on intrinsic dynamical variables
rather than analyzing the structural aspects of the phase space. In
particular, we recall that the canonical form of an integrable system is
Arnold and Avez (1968); Arnol’d (2013); Lichtenberg and Lieberman (2013):

$\begin{array}[]{l}J^{\prime}=J,\\\\[7.11317pt]
\psi^{\prime}=\psi+2\,\pi\,\nu,\end{array}\qquad\qquad\qquad\qquad\begin{array}[]{l}J_{k}=J_{0},\\\\[7.11317pt]
\psi_{k}=\psi_{0}+2\,\pi\,k\,\nu,\end{array}$

where $J_{k}$ and $\psi_{k}$ are the action and angle variables, respectively,
at iteration $k$, while $\nu$ denotes the rotation number. This rotation
number captures the average angular increment of the angle variable during a
single iteration of the map and remains invariant regardless of the mapping’s
representation. For instance, in $(q,p)$-coordinates, it can be analytically
expressed using Danilov’s theorem Zolkin  _et al._ (2017); Nagaitsev and
Zolkin (2020); Mitchell  _et al._ (2021) or as a limit:

$\nu=\frac{1}{2\,\pi}\,\lim_{N\rightarrow\infty}\frac{\mathcal{M}^{N}(\theta_{0})-\theta_{0}}{N}=\frac{1}{2\,\pi}\,\lim_{N\rightarrow\infty}\frac{\theta_{N}-\theta_{0}}{N},\qquad\qquad\qquad\qquad\theta_{k}=\arctan\frac{p_{k}-p_{*}}{q_{k}-q_{*}},\,\,\,k\in[0,N],$
(9)

where $\theta$ is the polar angle measured with respect to the fixed point.
For chaotic trajectories, this limit is non-existent due to the erratic and
unpredictable behavior of the trajectory’s angle, which prevents the
establishment of a well-defined average angular increment.

Furthermore, it is important to recognize that in the case of a stable
integrable map featuring polygonal invariants, the phase space is effectively
segmented into layers defined by the pieces of the force function. These
layers are “demarcated” by separatrices possessing rational rotation numbers
and are composed of polygons with identical numbers of vertices/sides. The
utilization of Danilov’s theorem uncovers that the Poincaré rotation number
must be of the form Zolkin  _et al._ (2022)

$\nu(J)=\frac{a_{0}+a_{1}\,J}{b_{0}+b_{1}\,J}$ (10)

within layers characterized by nonlinear integrable motion, or, it is
simplified to

$\nu(J)=\frac{a_{0}}{b_{0}}$

if the motion within a layer is mode-locked around a fixed point. Here
$a_{0,1}$ and $b_{0,1}$ represent specific constants.

Observing that Eq. (10) constitutes a monotonic function of amplitude enables
us to distinguish integrable layers from those involving chaotic trajectories.
In the latter scenario, the rotation number will not exhibit monotonic
behavior, but rather showcase discontinuous fluctuations during numerical
experiments. To provide a visual representation of the concept underlying this
new approach, we employ three mappings characterized by force functions
derived from the subsequent parameter sets:

$\begin{array}[]{ll}\left[\\{0;1\\},\\{-1;1\\},\\{-2;1\\},\\{-1;1\\};-5\right],&\qquad\qquad\qquad(\text{integrable
map}),\\\\[7.11317pt]
\left[\\{-2;1\\},\\{1;2\\};1\right],&\qquad\qquad\qquad(\text{chaos in
cell}),\\\\[7.11317pt]
\left[\\{1;1\\},\\{0;2\\},\\{-1;1\\},\\{0;1\\};1/2\right],&\qquad\qquad\qquad(\text{chaotic
diffusion}).\end{array}$

The top row of Figure 5 presents the corresponding phase-space portraits. Each
plot encloses a rectangle with sides measuring two lattice periods $2\,L$,
surrounding the central cell (as indicated by dashed cyan lines). The two
bottom rows then illustrate the rotation number within this area and evaluated
along both: the first and second symmetry lines as functions of the distance
from the fixed point.

Figure 5:  Phase space portraits (top row) and spectra along the first (middle
row) and second symmetry lines (bottom). The rotation number is evaluated
within the dashed cyan square around the fixed point, which encloses 2 lattice
periods. $\nu_{1,2}$ are depicted in solid black for trajectories mode-locked
around a fixed point, in blue for corresponding layers with nonlinear
integrable motion, gray for chaotic regions, and dashed black for global mode-
locking of the lattice. Red dots represent the locations of corresponding
separatrices.

First, let’s consider a slightly more intricate example of a nonlinear
integrable system featuring a semi-regular tiling pattern. The central cell of
this system corresponds to a fragment of an integrable map with a standard
piecewise linear function made up of five segments. Surrounding a fixed point,
there are three layers: a line segment with 2-cycles (indicated by a black
line), followed by two layers comprising concentric polygons with nonlinear
dynamics (blue contours separated by red separatrices). The rotation number
along the symmetry lines inside the central cell can be expressed as follows:

$\nu_{1}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{2},&\qquad\displaystyle
s=0,\\\\[9.95863pt]
\displaystyle\frac{1+4\,|s|}{2+12\,|s|},&\qquad\displaystyle
0<|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{3+2\,(|s|-1/2)}{8+8\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq
1,\end{array}\right.\qquad\qquad\mathrm{and}\qquad\qquad\nu_{2}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{2},&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{1+2\,(|s|-1/2)}{2+6\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq\frac{3}{2},\\\\[9.95863pt]
\displaystyle\frac{3+(|s|-3/2)}{8+4\,(|s|-3/2)},&\qquad\displaystyle\frac{3}{2}<|s|\leq\frac{5}{2}.\end{array}\right.$

The central cell is replicated across the phase plane, with gaps between the
bounding octagons being filled with invariant triangles. However, in contrast
to the amplitude-dependent dynamics within the replicated cells, the gaps with
invariant triangles are degenerate in nature: initial conditions on such
triangles visit each island once and return to the exact position where they
started. Both the replicated cells and the invariant triangles are organized
into chains of three islands, and the global rotation number is mode-locked at
the value of $1/3$.

Figure 6:  Diffusion in the phase space for two chaotic orbits originating
from the top-right plot in Fig. 5. Each plot shows the phase space area
occupied after $N$ iterations. The value of $N$ is provided in the frame at
the top-left corner, while the cyan rectangle corresponds to the area occupied
in the adjacent plot to the left.

Moving on to the next example in the middle, we delve into a scenario that
addresses the notion of “moderate perturbation”, primarily affecting the small
and intermediate amplitudes of the system. To illustrate this situation, we
have opted for an extreme case where only a particular set of surviving tori
remains intact, forming a network of separatrices. These isolating invariants
wield a significant influence on the system’s dynamics by effectively
segregating cells of periodicity. All points located on the separatrix possess
a period equal to 4, which gives rise to a global mode locking of $1/4$.
Within the central cell of periodicity, two additional invariant structures of
unit measure emerge: a line segment featuring 2-cycles around the fixed point
(depicted as a black line) and a chain of three islands composed of line
segments (highlighted in magenta). The clusters of dots, varying in shades of
gray, represent three chaotic trajectories that are either “locked” within the
central cell or enclosed by the surrounding cells grouped in fours (as
illustrated in Fig. 3). This particular scenario is termed as the chaos in
cell, denoting the presence of chaotic behavior within a confined region.

When examining the corresponding spectra along the symmetry lines, it becomes
evident that, for nearly all initial conditions within the central cell
(excluding non-isolated 2- and 3-cycles), the rotation numbers exhibit chaotic
signatures akin to noise. Outside the central cell, as anticipated, the
rotation number averages out and becomes mode-locked to $1/4$.

Finally, we consider the scenario where perturbations primarily impact the
larger amplitudes of the system (within the lattice period), causing a
complete disruption of the invariant isolating structures. In this selected
example, the central cell displays two distinct layers surrounding the fixed
point: concentric hexagons (depicted in black) and octagons (highlighted in
blue). The inner layer experiences degeneracy, with every initial condition
being periodic with a period of 6. On the other hand, the second layer
represents nonlinear integrable trajectories with a rotation number described
by the form (10):

$\nu_{1}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{6},&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{1+4\,(|s|-1/2)}{6+16\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq
1,\end{array}\right.\qquad\qquad\mathrm{and}\qquad\qquad\nu_{2}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{6},&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{1+2\,(|s|-1/2)}{6+8\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq\frac{3}{2}.\end{array}\right.$

The tiling follows a semi-regular pattern with another fragment of an
integrable system composed of layers made of triangles, pentagons, and
octagons. Beyond the larger surviving octagons that encompass both types of
tiles (marked in red), a region of chaos unfolds. Notably, some invariant
structures still persist, such as a chain of 9 islands comprising a line
segment with a rotation number of $2/9$ (indicated in magenta). However, for
most initial conditions, the trajectory fills the space between the integrable
fragments in an irregular manner (depicted by black and gray dots representing
two sample orbits).

However, due to the absence of a global isolating invariant, chaotic
trajectories can diffuse between lattice periods, in contrast to the previous
example of chaos in a cell. Figure 6 illustrates the trajectories of the same
two samples as shown in Fig. 5, but on a different scale in both space and
time. Initially, as evident from the first plot, for the first two thousand
iterations, one of the trajectories (gray) remains close to the limiting
octagon, while the other one (black) seems to be attracted to the area around
the 9 islands marked in magenta. However, as the number of iterations
increases, both orbits progressively diffuse out of the central cell and
expand further and further in the phase space, as depicted in the next four
plots of Fig. 6. For this reason, we will refer to this scenario as chaotic
diffusion.

Despite the differences between chaos in a cell and chaotic diffusion, an
important similarity exists: both spectra reveal chaotic trajectories through
the erratic behavior of the rotation number. In the case of chaotic diffusion,
symmetry lines periodically intersect integrable fragments alternating with
chaotic regions. This alternation is reflected in the spectra by the presence
of horizontal lines indicating global mode-locking, and noisy patterns
representing chaotic zones. Refer to the two bottom-right plots in Fig. 5 for
illustration.

###  IV.1 Search Algorithm

To summarize the concept discussed above, we introduce a variable
$\mathrm{jit}$ that reveals rapid oscillations or deviations in the rotation
number from the anticipated piecewise monotonic pattern characterizing an
integrable scenario. In the context of a completely integrable system, the
value of $\mathrm{jit}$ is approximately equivalent to the number of layers
present

$\mathrm{jit}\approx\text{number of layers}.$

In contrast, for scenarios involving chaotic behavior, $\mathrm{jit}$ is
roughly analogous to the number of chaotic samples within the layers

$\mathrm{jit}\approx\text{number of chaotic samples in layers}.$

By densely populating symmetry lines with sample points, we can effectively
distinguish between these two behaviors.

In our pursuit of discovering new mappings, we expand the force function by
incorporating more than two linear segments within a single period $L$, and
then systematically explore the parameter space. Our algorithm bears
resemblance to the one utilized in Zolkin  _et al._ (2022) although the
primary emphasis lies in accurately computing the rotation number rather than
merely enumerating vertices of individual polygonal orbits. The specific
processes and decision points for the algorithm are detailed below:

  1. 0.

Input The search process involves exploring different force functions
$f(q,\mathbf{k},\mathbf{l},d)$ with integer vectors $\mathbf{k}$,
$\mathbf{l}$, and scalar $d$. While the initial scan focused on integer slopes
$\mathbf{k}$, some results can be generalized analytically for more general
values of the shift parameter $d\in\mathbb{R}$ and segment lengths
$\mathbf{l}\in\mathbb{R}_{+}^{n}$.

  2. 1.

Fixed point A crucial requirement for organizing the phase space into cells of
periodicity is the existence of at least one fixed point
$\zeta_{*}=(q_{*},p_{*})$, where $q_{*}=p_{*}=f(q_{*})/2$. Maps lacking a
fixed point are excluded from consideration.

  3. 2.

Orbit tracking For the chosen maps, we trace their orbits
$\vec{\zeta}=(\zeta_{0},\ldots,\zeta_{N})$ for various initial conditions
$\zeta_{0}=(q_{0},p_{0})$ along both symmetry lines:

$p_{0}=q_{0}\qquad\qquad\qquad\mathrm{and}\qquad\qquad\qquad
p_{0}=f(q_{0})/2.$

Initial points are placed within two lattice periods around the fixed point,
$|q_{0}-q_{*}|<L$, ensuring multiple seeds for every $j$-th layer
$q_{0}\in(\sum_{i=1}^{j}l_{i};\sum_{i=1}^{j+1}l_{i})$.

  4. 3.

Stability Check Although mappings possess fixed points, they might still
exhibit unstable trajectories at larger amplitudes. If we observe unbounded
growth of the radial distance from the fixed point for certain initial
conditions, i.e., $\exists\,i:\,\,|\zeta_{i}-\zeta^{*}|>r_{\mathrm{max}}$, the
map is excluded from further analysis. In particular, mappings with chaotic
diffusion might be rejected at this stage. We set $\displaystyle
r_{\mathrm{max}}=4\,L$.

  5. 4.

Rotation Number Calculation Next, we evaluate the rotation number for each
initial condition as:

$\nu\approx\frac{\theta_{N}-\theta_{0}}{2\,\pi\,N}$

where sufficiently large $N$ ensures convergence for invariant level sets.

  6. 5.

Jitter count To quantify the variable $\mathrm{jit}$, we count the number of
times the rotation number changes its monotonicity as a function of the
coordinate:

$\mathrm{jit}=\mathrm{max}\left[\mathrm{jit}_{1},\mathrm{jit}_{2}\right],\qquad\qquad\qquad\qquad\mathrm{jit}_{1,2}=\sum_{|q_{i}-q_{*}|\leq
L}\frac{1-\mathrm{sgn}\,\left\\{[\nu_{1,2}(q_{i})-\nu_{1,2}(q_{i-1})][\nu_{1,2}(q_{i+1})-\nu_{1,2}(q_{i})]\right\\}}{2}.$

  7. 6.

Selection In the final stage, we compare the values of $\mathrm{jit}_{1,2}$
calculated along both symmetry lines with the number of samples in each layer.
If

$\mathrm{jit}\approx n\ll\text{number of samples in layers}$

we select the map and then analytically verify its integrability.

##  V  Results

###  V.1 Map $\mathcal{M}_{\mathrm{G}}^{6}$

During our parameter scan for functions composed of two pieces and integer
coefficients, no additional integrable cases were discovered. However, as we
proceeded to consider $f(q)$ formed from three line segments, we encountered
an intriguing map:

$\mathcal{M}_{\mathrm{G}}^{6}:\qquad
f(q)=[\\{1;2\\},\\{2;1\\},\\{0;1\\};4\,m].$

This specific map doesn’t allow for any choice in selecting the configuration
of the center and doesn’t contain any free parameters except for scaling and
trivial translations. The phase space diagram and rotation number along
symmetry lines are depicted in Figure 7, where $\nu_{1,2}$ are given by:

$\nu_{1}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{5},&\qquad\displaystyle|s|\leq
1,\\\\[9.95863pt] \displaystyle\frac{1}{5+(|s|-1)},&\qquad\displaystyle
1<|s|\leq
2,\end{array}\right.\qquad\mathrm{and}\qquad\nu_{2}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{5},&\qquad\displaystyle|s-1/2|\leq
3/2,\\\\[9.95863pt]
\displaystyle\frac{1}{5+(|s-1/2|-3/2)},&\qquad\displaystyle 3/2<|s-1/2|\leq
5/2.\end{array}\right.$

The central cell has two layers. The inner core consists of nested pentagons
surrounded by a chain of 5 degenerate triangular islands (level sets shown in
black). All initial conditions falling within this layer exhibit periodic
behavior with a rotation number equal to $1/5$. On the other hand, the outer
layer is composed of concentric dodecagons (level sets shown in blue) and
exhibits nonlinear integrable motion, which includes periodic or quasiperiodic
orbits. The surrounding cells are organized into groups of 6, thus
characterizing global mode-locking at a rotation number of $1/6$.

Figure 7:  The left plot shows the phase space portrait for the map
$\mathcal{M}_{\mathrm{G}}^{6}$. Level sets within the central cell are
depicted in solid black for trajectories mode-locked around a fixed point, in
blue for the corresponding layer with nonlinear integrable motion, and in red
for separatrices isolating different layers or tiles. The solid and dashed
green lines in the diagram correspond to the second and first symmetry lines,
respectively, $p=f(q)/2$ and $p=q$. The connected series of line segments
shown in brown depict one period of the force function $f(q)\in[0;L]$ along
with the values for their slopes $k_{i}$ and length $l_{i}$. The two plots on
the right provide corresponding spectra along the first (top) and second
(bottom) symmetry lines. The rotation number is evaluated within 2 lattice
periods around the fixed point (phase space region enclosed by the dashed cyan
square) and depicted with the same colors as corresponding level sets or with
dashed black for initial conditions globally mode-locked around cc.

###  V.2 Mappings $\mathcal{M}_{\alpha 3}^{4}$ and
$\mathcal{M}_{\mathrm{F}1^{\prime}}^{4}$

The next integrable system we discovered for forces composed of three segments
is the mapping:

$\begin{array}[]{rl}\displaystyle\mathcal{M}_{\alpha
3}^{4}:\qquad\mathrm{[cc]}&\quad
f(q)=[\\{0;l_{1}\\},\\{-1;1\\},\\{1;1\\};2\,m\,(2+l_{1})-2\,\,],\\\\[7.11317pt]
\qquad\mathrm{[dc]}&\quad
f(q)=[\\{0;l_{1}\\},\\{-1;1\\},\\{1;1\\};2\,m\,(2+l_{1})+l_{1}],\end{array}\qquad
l_{1}\geq 1,$

with one real parameter $l_{1}$ and two possible configurations of the central
cell. By adding a fourth piece, it can be generalized to another map with two
parameters as follows:

$\begin{array}[]{rl}\displaystyle\mathcal{M}_{\mathrm{F}1^{\prime}}^{4}:\qquad\mathrm{[cc]}&\quad
f(q)=[\\{0;l_{1}\\},\\{-1;1\\},\\{0;l_{3}\\},\\{1;1\\};2\,m\,(2+l_{1}+l_{3})-2-l_{3}],\\\\[7.11317pt]
\qquad\mathrm{[dc]}&\quad
f(q)=[\\{0;l_{1}\\},\\{-1;1\\},\\{0;l_{3}\\},\\{1;1\\};2\,m\,(2+l_{1}+l_{3})+l_{1}],\end{array}\qquad
l_{1}\geq 1,\,\,l_{3}>0.$

Both systems have $F=0$, resulting in a global period and rotation number
equal to $4$ and $1/4$, respectively. The tiling is semi-regular with two
types of cells grouped in chains of four islands around the central cell (see
Fig. 3 for the reference). The configuration labeled as dc stands for
degenerate center and can be seen either as a central cell where all initial
conditions are mode-locked at a global rotation number, or as a “thick”
central node that is part of a net of separatrices. Phase space plots,
bifurcation diagrams for tiles, and examples of the rotation number for
different values of parameters are provided in Fig. 8. Changes in the
parameter $l_{1}$ only affect the degenerate cell, while changes in $l_{3}$
only impact the cell with nonlinear dynamics. For analytical expressions of
rotation numbers, the reader can refer to mappings $\alpha 3$ and
$\mathrm{F}1^{\prime}$ in Zolkin  _et al._ (2022).

Figure 8:  Phase space portraits for mappings $\mathcal{M}_{\alpha 3}^{4}$ and
$\mathcal{M}_{\mathrm{F}1^{\prime}}^{4}$ (top left) and bifurcation diagram
for the central cell and degenerate center as functions of $l_{3}$ and $l_{1}$
respectively (top right). The bottom row shows spectra for the cc
configuration and different values of $l_{3}$ (shown in different colors).

###  V.3 Mappings $\mathcal{M}_{\mathrm{E}1^{\prime}}^{3}$,
$\mathcal{M}_{\mathrm{a}2}^{6}$ and $\mathcal{M}_{\mathrm{b}2}^{3}$

The final group of systems to be presented in this section involves the
extension of the mapping $\mathcal{M}_{\mathrm{a}1}^{6}$ to a force function
composed of four pieces and $\mathcal{M}_{\mathrm{b}1}^{3}$ to forces made of
three or four segments. Respectively

$\mathcal{M}_{\mathrm{a}2}^{6}:\qquad\,\mathrm{[cc]}\quad
f(q)=[\\{0,1\\},\\{1,l\\},\\{2,1\\},\\{1,l\\};2\,m\,(1+l)-l],\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
l>0,$

where the new segments should have equal length, $l_{2}=l_{4}=l$, and

$\begin{array}[]{rlr}\displaystyle\mathcal{M}_{\mathrm{E}1^{\prime}}^{3}:\qquad\mathrm{[cc]}&\quad
f(q)=[\\{0,1\\},\\{-2,1\\},\\{-1,l_{4}\\};3\,m\,(2+l_{4})+4+l_{4}],&\\\\[7.11317pt]
\qquad\mathrm{[cn]}&\quad
f(q)=[\\{0,1\\},\\{-2,1\\},\\{-1,l_{4}\\};3\,m\,(2+l_{4})+2],&\qquad
l_{4}>0,\\\\[7.11317pt] \qquad\mathrm{[dc]}&\quad
f(q)=[\\{0,1\\},\\{-2,1\\},\\{-1,l_{4}\\};3\,m\,(2+l_{4})-l_{4}],&\\\\[7.11317pt]
\displaystyle\mathcal{M}_{\mathrm{b}2}^{3}:\qquad\mathrm{[cc]}&\quad
f(q)=[\\{0,1\\},\\{-1,l_{2}\\},\\{-2,1\\},\\{-1,l_{4}\\};3\,m\,(2+l_{2}+l_{4})+4+2\,l_{2}+l_{4}],&\\\\[7.11317pt]
\qquad\mathrm{[dc]}&\quad
f(q)=[\\{0,1\\},\\{-1,l_{2}\\},\\{-2,1\\},\\{-1,l_{4}\\};3\,m\,(2+l_{2}+l_{4})+2+l_{2}],&\qquad
l_{4}\geq l_{2}>0,\\\\[7.11317pt] \qquad\mathrm{[dc]}&\quad
f(q)=[\\{0,1\\},\\{-1,l_{2}\\},\\{-2,1\\},\\{-1,l_{4}\\};3\,m\,(2+l_{2}+l_{4})-l_{4}],&\end{array}$

where the condition $l_{4}\geq l_{2}$ ensures the exclusion of twin maps.

Figure 9:  Phase space portraits and bifurcation diagram for mappings
$\mathcal{M}_{\mathrm{E}1^{\prime}}^{3}$ and $\mathcal{M}_{\mathrm{b}2}^{3}$.

Figures 9 and 10 provide bifurcation diagrams for cells with nonlinear
dynamics along with phase space plots. In all systems, the tiling is semi-
regular: the plane is tessellated with octagon or hexagon separatrices (red
contours), enclosing cells with nonlinear dynamics for force functions
composed of four or three segments, respectively. The gaps in between are
filled with triangles representing the degenerate center and/or chains of
linear islands (black contours). The map $\mathcal{M}_{\mathrm{a}2}^{6}$ has a
singular configuration, featuring a central cell surrounded by a chain of six
identical islands. On the other hand, $\mathcal{M}_{\mathrm{b}2}^{3}$ allows
for three distinct configurations and encompasses two independent groups of
chains, each composed of three islands. In the case of the latter, the fixed
point can be situated either within the central cell or within one of the two
triangular tiles representing degenerate centers. When $l_{2}=0$ (map
$\mathcal{M}_{\mathrm{E}1^{\prime}}^{3}$), one of the triangular chains
vanishes, leading to the corresponding degenerate center evolving into a
central node. Moreover, the right-hand column of plots in Figure 10 showcases
sample spectra for various parameter values of $l$. The expressions for
$\nu_{1,2}$ are as follows:

$\nu_{j}=\left\\{\begin{array}[]{ll}\displaystyle
0,&\qquad\displaystyle|s|\leq\frac{1}{2},\\\\[9.95863pt]
\displaystyle\frac{|s|-1/2}{2+6\,(|s|-1/2)},&\qquad\displaystyle\frac{1}{2}<|s|\leq\frac{1}{2}+l,\\\\[9.95863pt]
\displaystyle\frac{l\,j+2\,(|s|-1/2-l)}{(2+6\,l)\,j+8\,(|s|-1/2-l)},&\qquad\displaystyle\frac{1}{2}+l<|s|\leq\frac{1}{2}+\frac{j}{2}+l,\end{array}\right.\qquad\qquad\qquad\qquad
j=1,2.$

Figure 10:  Phase space portrait for the map $\mathcal{M}_{\mathrm{a}2}^{6}$.
Bifurcation diagram and spectra for different values of parameter $l$ (shown
with different colors) are provided to the left and right respectively.

##  VI  Integrable polygon maps on tor

Until now, our focus has solely centered on continuous forces with arithmetic
quasiperiodicity, as defined by Eq. (5). This section broadens our
investigation, by studying the dynamics on a torus $(q,p)\in\mathbb{T}^{2}$,
where the equations of motion are considered modulo $L$:

$\begin{array}[]{lll}\mathcal{M}^{\mathrm{tor}}:&\quad
q^{\prime}=p&\quad\mathrm{mod}\,L,\\\\[7.11317pt] &\quad
p^{\prime}=-q+f_{\mathrm{p.l.}}(p)&\quad\mathrm{mod}\,L.\end{array}$ (11)

A systematic scan in the parameter space has revealed a significant number of
new integrable systems. At this juncture, the challenge for mappings
integrable on a torus lies more in the classification of new results and,
ultimately, in comprehending the totality of all possibilities than in merely
expanding the number of new findings. For instance, while 2-piece force
functions with arithmetic quasiperiodicity offered only two systems for any
ratio of $l_{2}/l_{1}$, exploration of dynamics modulo $L$ unveils 6
integrable scenarios, specifically for $l_{1}=l_{2}=L/2$. Fig. 11 presents
corresponding phase space portraits along with a projection of the phase space
onto a toroidal surface, where $p$ and $q$ are treated as angular variables.
We will confine our presentation to these examples, progressing to discuss
another potential generalization of the force function.

Figure 11:  Plots (a.1) – (f.1) show phase space portraits for mappings
integrable on a torus with $l_{1}=l_{2}=L/2$, and corresponding triples of
slopes and shift parameter $[k_{1},k_{2};d]$: $[0,2;L]$, $[-2,0;2\,L]$,
$[1,2;0]$, $[0,-1;L]$, $[0,1;L]$ and $[0,1;L/2]$. Blue and black level sets
represent initial conditions in layers with nonlinear and linear dynamics
respectively. Level set in cyan (e.1) shows a chain of three islands such that
all initial conditions are 9-cycles and set in orange (c.1) are 4-cycles.
Plots (d.2) – (f.2) provide a projection onto a toroidal surface for
corresponding cases.

Figure 12:  Plot (d.1) reproduces the phase space portrait on a torus from
Fig. 11 along with a sample orbit, 7-cycle shown in blue. Plots (d.1′) and
(d.1′′) illustrate the outcomes obtained by constructing a continuous
quasiperiodic function. The last plot (d.1′′′) presents the integrable
tessellation corresponding to $f(q)$ with discontinuity along with a sample
orbit, $7\times 4$-cycle, which differs from the sample in (d.1) by one
lattice period, $q_{0}\rightarrow q_{0}+L$.

###  VI.1  Unwrapping and integrable diffusion

Initially, it is noteworthy that due to the translational symmetry of the
invariant

$\mathcal{K}_{\mathrm{a.q.}}(p,q)\rightarrow\mathcal{K}_{\mathrm{a.q.}}(p+m\,L,q+n\,L),\qquad\qquad\qquad
m,\,n\in\mathbb{Z},$

every integrable map with a continuous force function exhibiting arithmetic
quasiperiodicity can be “wrapped” on a torus ($\mathrm{mod}\,L$), preserving
integrability, even if $F\neq 0$. In simpler terms, if we start with a regular
piecewise linear force $f_{\mathrm{p.l.}}(q)$ that maintains continuity in the
interval $[0;L]$ and we are able to create an integrable map with arithmetic
quasiperiodicity, then the associated system on a torus will also exhibit
integrability. In fact, the two examples mentioned above, namely (a.1) and
(b.1), can be concluded from $\mathcal{M}_{\mathrm{a}1}^{6}$ and
$\mathcal{M}_{\mathrm{b}1}^{3}$.

Conversely, the reverse is not assured. Having an integrable map on a torus,
when considering the dynamics in $\mathbb{R}^{2}$ by “unwrapping” the force
from the interval $[0;L]$ with Eq. (5), the outcome may lack integrability.
This is precisely why mappings (c.1) – (f.1) were absent in our exploration of
2-piece forces with arithmetic quasiperiodicity. Plots (d.1′) and (d.1′′) in
Fig. 12 illustrate two attempts to maintain the continuity of the force
function while unwrapping case (d.1) with respect to both possible
configurations of the center. As evident, for small amplitudes, approximately
on the order of $L$, there exists a continuous fibration of the phase space
with invariant level sets. However, at larger amplitudes, chaotic trajectories
emerge, indicating that the global dynamics lack integrability.

Meanwhile, as mentioned at the beginning, we can forsake continuity by using
equation (7)

$f_{\mathrm{tor}}(q)=f_{\mathrm{p.l.}}(q\,\,\mathrm{mod}\,\,L)\quad\mathrm{mod}\,\,L$

instead of (5). This procedure ensures integrability in $\mathbb{R}^{2}$ for
the map in MH form by construction: by replicating invariant level sets from a
torus across the plane, we automatically satisfy both Eqs. (4) with respect to
$f_{\mathrm{tor}}(q)$. Plot (d.1′′′) in Fig. 12 illustrates a case of
“correct” unwrapping, resulting in discontinuity of force. Phase portraits
(d.1) and (d.1′′′) both contain sample orbits separated by one lattice period,
$L$. In the case of map (d.1), due to degeneracy Zolkin  _et al._ (2022), all
initial conditions form either 7- or 3-cycles. As one would expect, in
situation (d.1′′′), the sample orbit departing from a non-central cell will be
mode-locked to $1/4$ since the average slope $k_{0}$ (and $F$) is equal to
zero, see Eq. (8). It hops between 4 islands, forming $7\times 4$\- or
$3\times 4$-cycles respectively.

While unwrapping from the torus (7) might seem somewhat trivial, as the new
invariant level set is obtained simply by translation, we encounter unexpected
behavior in some cases. Refer to Figure 13, exploring the already familiar
example (a.1) from Figure 11. The first plot reproduces the phase space,
providing a few iterations of the sample orbit. The next plot, (a.1′),
illustrates several iterations of the initial condition within a central cell
for the case of $f_{\mathrm{tor}}(q)$. Unexpectedly, after several iterations,
the orbit escapes the central cell and starts hopping between tiles in a
complicated way, resembling chaotic diffusion. However, in this scenario,
instead of being “random,” the orbit remains on a particular level set, as
shown in Plot (a.1′′). Instead of a finite number of closed rectangles as we
had for the force with arithmetic quasiperiodicity (Figure 4), now we have
infinitely many sets of rectangles split into 2 pieces. Finally, Plot (a.1′′′)
provides $10^{4}$ iterations of the same initial condition as in (a.1′) on a
different scale, showing diffusion similar to what we observed in Figure 6.
Notice how progression in time reveals more and more segments of the
invariant, which matches the pattern in (a.1′′). We will refer to this
behavior as integrable diffusion.

Figure 13:  Plot (a.1) reproduces the phase space portrait on a torus from
Fig. 11 along with a sample orbit, shown with blue dots. Plot (a.1′) presents
the integrable tessellation corresponding to $f_{\mathrm{tor}}(q)$ with
discontinuity and sample orbit originating form the central cell. Next plots
show the corresponding invariant level sets (a.1′′) and $10^{4}$ iterations of
the same orbit on a different scale. Cyan rectangle corresponds to the area
occupied in (a.1′).

Now, we can pose two questions: “Which examples from Fig. 11 will behave
similarly to (a.1) or to (d.1) when unwrapped according to Eq. (7)?” and “Why
is that?” The answer to the first question is that the group (a.1) – (c.1)
behaves like (a.1), and the group (d.1) – (f.1) behaves like (d.1). At this
point, the reader is encouraged to take a break and revisit Fig. 11 in an
effort to uncover the answer to the question “Why?”

The resolution to this puzzle lies in how the network of separatrices (shown
in red) is organized within the square representing torus. In all cases (d.1)
– (f.1), separatrices form closed loops and are confined within
$\mathbb{T}^{2}$, thus providing isolating invariants when the system is
unwrapped. On the other hand, in all cases (a.1) – (c.1), the boundary of the
torus cuts the separatrix, disrupting isolating invariants and allowing
diffusion.

##  VII  Applications of mappings with periodic and arithmetically
quasiperiodic forces.

Applications of integrable maps with polygonal invariants were extensively
discussed in Zolkin  _et al._ (2022). Here, we will briefly highlight the main
points and provide a few illustrations.

###  VII.1  Conection to McMillan-Suris mappings

The simplest 1-piece mappings on a torus with integer coefficients,

$f_{\mathrm{tor}}=k_{0}\,q\,\mod L,\qquad\qquad\qquad k_{0}=-2,-1,0,1,2,$

can be regarded as ultra-discrete or limiting cases of McMillan-Suris mappings
with a trigonometric polynomial,

$\mathcal{K}(p,q)=\Lambda_{1}\left[\cos(p-\Phi)+\cos(q-\Phi)\right]+\Lambda_{2}\cos(p+q-\Psi)+\Lambda_{3}\cos(p-q)$

and force function,

$f(q)=2\,\arctan\left[\frac{\Lambda_{1}\sin\Phi+\Lambda_{2}\sin(\Psi-q)+\Lambda_{3}\sin
q}{\Lambda_{1}\cos\Phi+\Lambda_{2}\cos(\Psi-q)+\Lambda_{3}\cos q}\right].$

For instance, by setting $\delta=\Lambda_{3}/\Lambda_{1}$ and
$\Lambda_{2}=\Phi=0$, we obtain

$\mathcal{K}(p,q)=\cos p+\cos q+\delta\,\cos(p-q),\qquad\qquad
f(q)=2\,\arctan\left[\frac{\delta\,\sin q}{1+\delta\,\cos q}\right].$ (12)

As seen, when $\delta=1$, the force function reduces to
$f(q)=2\,\arctan\left[\tan(q/2)\right]$. Depending on whether we assume the
continuity of $f(q)$ or take the principal value of $\arctan$, we obtain
either $f(q)=q$ or $f(q)=q\,(\mathrm{mod}\,L)$ with $L=2\,\pi$. This
corresponds to unwrapping with arithmetic quasiperiodicity (and zero periodic
part) or to a purely periodic function with discontinuities. The associated
invariant of motion possesses degeneracy, resulting in infinitely many
possible invariants of motion, see Zolkin  _et al._ (2022), including the
smooth form (12) and the polygon tessellation. The middle plot in Fig. 14
shows the polygon invariant for $\delta=1$, surrounded by complimentary cases
with $\delta\lessgtr 1$.

Figure 14:  Phase space portraits for McMillan-Suris mappings with
trigonometric polynomial 12. Different plots correspond to different values of
$\delta$. When $\delta\neq 1$, the dashed brown line shows linear flow on a
torus (case $\delta=1$) for the reference.

###  VII.2  Smoothening procedure

While some of these mappings offer insight into the limiting integrable
dynamics with a smooth $f(q)$, as seen in the previous section, there is a
more important application applicable to all discovered systems. By
“smoothening” $f(q)$ Rychlik and Torgerson (1998); Zolkin  _et al._ (2022), we
can obtain quasi-integrable systems where the phase space volume occupied by
chaotic trajectories decreases with the smoothness parameter, $\epsilon$. In
addition to the ability to suppress chaos, this procedure allows the physical
realization not only for forces with arithmetic quasiperiodicity but even for
periodic functions with discontinuity, where vertical segments can be
approximated by very steep and short parts of the force function. In fact,
Fig. 14 provides two examples of “perfect smoothening” (since resulting
systems are perfectly integrable) with respect to the parameter
$\epsilon=1-\delta$.

By using other rules, as long as the new function is close enough to its
piecewise linear prototype, we can obtain near-integrable dynamics. For
example, for $l_{1,2,3,4}=1$, we can approximate forces for mappings
$\mathcal{M}^{6}_{\mathrm{a2}}$ and $\mathcal{M}^{4}_{\mathrm{F1^{\prime}}}$
as $f_{\mathrm{tr.w.}}(q+1)+q$ and $f_{\mathrm{tr.w.}}(q)$, respectively,
where we introduced smooth differentiable forms (when $\epsilon\neq 0$) of the
trapezoid and the square waves with period $L=4$:

$f_{\mathrm{tr.w.}}(q)=\frac{1}{\pi}\left\\{\arccos\left[(1-\epsilon)\,\cos\frac{\pi\,q}{2}\right]-\arcsin\left[(1-\epsilon)\,\sin\frac{\pi\,q}{2}\right]\right\\}\qquad\mathrm{and}\qquad
f_{\mathrm{sq.w.}}(q)=\frac{4}{\pi}\,\arctan\left[\frac{1}{\epsilon}\,\sin\frac{\pi\,q}{2}\right].$

Figure 15 provides an illustration of phase space portraits for these systems
obtained by tracking for $\epsilon=10^{-1}$. The last example on the right
depicts tracking for the smooth square wave force function, which, for
$\epsilon=0$, represents another limiting case of the McMillan-Suris map with
$k_{0}=0$.

Figure 15:  Illustration of “smoothening” procedure. Examples from left to
right correspond to forces with arithmetic quasiperiodicity
($\mathcal{M}^{6}_{\mathrm{a2}}$), purely periodic
($\mathcal{M}^{4}_{\mathrm{F1^{\prime}}}$), and periodic function with
discontinuity (square wave).

In particular, this method allows the design of nonlinear quasi-integrable
kick functions corresponding to thin RF cavities when considering a simplified
picture of longitudinal motion in particle accelerators. Even for a single
degree of freedom, as highlighted by Chirikov’s early findings Chirikov (1969,
1979), the typical kick from the RF field on a particle ($\propto\sin q$,
taking $q$ as the longitudinal coordinate) results in dynamics that exhibit
chaos and can lead to global stochasticity. The smoothening procedure, indeed,
provides an alternative to the only known integrable case of the McMillan-
Suris map. Cases with $F=0$ represent longitudinal dynamics for the stationary
bucket (indicating no acceleration), while $F\neq 0$ represents regimes with
net acceleration/deceleration.

##  VIII  Summary

In this paper, we have uncovered a novel class of integrable symplectic maps
characterized by a polygonal tessellation and corresponding fibration of the
phase space. While polygonal tilings of the plane have found applications in
diverse fields such as crystallography, the study of topological phases of
matter, and group theory, our work introduces a distinct category of periodic
regular and semi-regular tessellations directly linked to the presence of
integrability. Our investigation focuses on systems in the McMillan-Hénon
form, featuring two specific types of piecewise linear force functions,
namely, (i) functions with arithmetic quasiperiodicity and (ii) periodic
functions with discontinuities. The mentioned subsets collectively form a
broader category of systems defined on a torus.

These nonlinear maps manifest non-trivial dynamical regimes, encompassing
various configurations of the center, mode-locking, and integrable diffusion.
In the latter scenario, which has not been previously documented to the best
of our knowledge, a particle quasi-randomly hops between tiles, gradually
“diffusing” away while remaining constrained to a constant level set of
invariant.

To conduct a systematic search for such maps, we developed an automated
algorithm that evaluates the rotation number along symmetry lines. Integrable
systems of this kind are distinguished by their piecewise monotonic rotation
number behavior concerning the initial amplitude. In contrast, chaotic systems
display irregular and “jittery” behavior, enabling a clear distinction between
the two. Our search focused on force functions with integer slope
coefficients, however, it remains an open question whether integrable maps
with polygonal invariants also exist for rational, and especially irrational,
slope coefficients. While we have presented numerous new mappings throughout
this paper, our algorithm, especially when not restricted to arithmetic
quasiperiodicity, unveils many more systems than we can present in a
methodical and coherent manner. This emphasizes the importance of a more
fundamental challenge: providing all possible integrable configurations from
the first principles.

##  IX Acknowledgments

The authors would like to thank Taylor Nchako (Northwestern University) for
carefully reading this manuscript and for her helpful comments. 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.

## References

  * Meiss (1992) J. D. Meiss, Rev. Mod. Phys. 64, 795 (1992).
  * 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.
  * Danilov (2008) V. Danilov, Phys. Rev. ST Accel. Beams 11, 114001 (2008).
  * Veselov (1991) A. P. Veselov, Russian Mathematical Surveys 46, 1 (1991).
  * Turaev (2002) D. Turaev, Nonlinearity 16, 123 (2002).
  * Henon (1969) M. Henon, Quarterly of applied mathematics , 291 (1969).
  * Chirikov (1969) B. Chirikov, Nuclear Physics, Novosibirsk.(Engl. Trans. CERN Transl. 31, 1971) (1969).
  * Chirikov (1979) B. V. Chirikov, Physics reports 52, 263 (1979).
  * Suris (1989) Y. B. Suris, Functional Analysis and Its Applications 23, 74 (1989).
  * Brown (1983) M. Brown, Am. Math. Mon 90, 569 (1983).
  * Brown and Slifker (1985) M. Brown and J. F. Slifker, Am. Math. Mon 92, 218 (1985).
  * Brown (1993) M. Brown, Lecture notes in pure and applied mathematics , 83 (1993).
  * Cairns  _et al._ (2014) G. Cairns, Y. Nikolayevsky, and G. Rossiter, arXiv preprint arXiv:1407.3364 (2014).
  * Cairns  _et al._ (2016) G. Cairns, Y. Nikolayevsky, and G. Rossiter, The American Mathematical Monthly 123, 363 (2016).
  * Zolkin  _et al._ (2022) T. Zolkin, Y. Kharkov, and S. Nagaitsev, “Machine-assisted discovery of integrable symplectic mappings,” (2022), arXiv:2201.13133 [nlin.SI] .
  * Schuster and Just (2006) H. G. Schuster and W. Just,  _Deterministic chaos: an introduction_ (John Wiley & Sons, 2006).
  * Krippendorf  _et al._ (2021) S. Krippendorf, D. Lüst, and M. Syvaeri, Fortschritte der Physik 69, 2100057 (2021).
  * Lewis Jr (1961) D. C. Lewis Jr, Pacific Journal of Mathematics 11, 1077 (1961).
  * Smith  _et al._ (2023a) D. Smith, J. S. Myers, C. S. Kaplan, and C. Goodman-Strauss, arXiv preprint arXiv:2303.10798 (2023a).
  * Smith  _et al._ (2023b) D. Smith, J. S. Myers, C. S. Kaplan, and C. Goodman-Strauss, “A chiral aperiodic monotile,” (2023b), arXiv:2305.17743 [math.CO] .
  * Zeller  _et al._ (2021) R. Zeller, O. Delgado-Friedrichs, and D. H. Huson, Computer Aided Geometric Design 90, 102027 (2021).
  * Manolas  _et al._ (2022) I. Manolas, F. Laccone, G. Cherchi, L. Malomo, and P. Cignoni, Computers and Graphics 106, 141 (2022).
  * Arnold and Avez (1968) V. I. Arnold and A. Avez,  _Ergodic problems of classical mechanics_ , Vol. 9 (Benjamin, 1968).
  * Arnol’d (2013) V. I. Arnol’d,  _Mathematical methods of classical mechanics_ , Vol. 60 (Springer Science & Business Media, 2013).
  * Lichtenberg and Lieberman (2013) A. J. Lichtenberg and M. A. Lieberman,  _Regular and chaotic dynamics_ , Vol. 38 (Springer Science & Business Media, 2013).
  * 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, Physical Review Accelerators and Beams 23, 054001 (2020).
  * Mitchell  _et al._ (2021) C. E. Mitchell, R. D. Ryne, K. Hwang, S. Nagaitsev, and T. Zolkin, Physical Review E 103, 062216 (2021).
  * Rychlik and Torgerson (1998) M. Rychlik and M. Torgerson, New York J. Math 4, 74 (1998).