Machine-assisted discovery of integrable symplectic mappings

We present a new automated method for finding integrable symplectic maps of the plane. These dynamical systems possess a hidden symmetry associated with an existence of conserved quantities, i.e. integrals of motion. The core idea of the algorithm is based on the knowledge that the evolution of an integrable system in the phase space is restricted to a lower-dimensional submanifold. Limiting ourselves to polygon invariants of motion, we analyze the shape of individual trajectories thus successfully distinguishing integrable motion from chaotic cases. For example, our method rediscovers some of the famous McMillan-Suris integrable mappings and discrete Painlev\'e equations. In total, over 100 new integrable families are presented and analyzed; some of them are isolated in the space of parameters, and some of them are families with one parameter (or the ratio of parameters) being continuous or discrete. At the end of the paper, we suggest how newly discovered maps are related to a general 2D symplectic map via an introduction of discrete perturbation theory and propose a method on how to construct smooth near-integrable dynamical systems based on mappings with polygon invariants.

We present a new automated method for finding integrable symplectic maps of the plane.These dynamical systems possess a hidden symmetry associated with an existence of conserved quantities, i.e. integrals of motion.The core idea of the algorithm is based on the knowledge that the evolution of an integrable system in the phase space is restricted to a lower-dimensional submanifold.Limiting ourselves to polygon invariants of motion, we analyze the shape of individual trajectories thus successfully distinguishing integrable motion from chaotic cases.For example, our method rediscovers some of the famous McMillan-Suris integrable mappings and ultra-discrete Painlevé equations.In total, over 100 new integrable families are presented and analyzed; some of them are isolated in the space of parameters, and some of them are families with one parameter (or the ratio of parameters) being continuous or discrete.At the end of the paper, we suggest how newly discovered maps are related to a general 2D symplectic map via an introduction of discrete perturbation theory and propose a method on how to construct smooth near-integrable dynamical systems based on mappings with polygon invariants.

I. INTRODUCTION
Integrable systems play a fundamental role in mathematics and physics, particularly in dynamical system theory, due to well-understood dynamics.Integrability is generally defined as the existence of a sufficient number of functionally independent conserved quantities (integrals of motion) in involution, which are related to intricate hidden symmetries.Hence, in many cases the construction of integrable systems often requires some special fine-tuning.This makes them very difficult to discover.Although integrable models occupy measure-zero in the space of parameters, they play a crucial role in the theory of dynamical systems.In fact, they shed light on the properties of more generic nearintegrable cases via Kolmogorov-Arnold-Moser (KAM) theorem and various perturbation theories.
While the concept of integrability is most often discussed in the context of continuous-time systems, it also arises in discrete-time maps [1].Mappings are of profound importance in the field of dynamical system theory exhibiting a broad range of behaviours including dynamical chaos, bifurcations, strange attractors, and integrability [2].On the one hand, they can be used as an approximation of a continuous flow (since in many cases the discretized form can be easier to study) or as a reduction to a lower dimensional manifold, e.g., considering Poincaré cross-sections.A particularly valuable class of discrete dynamical systems are symplectic maps, associated with Hamiltonian flows, e.g., a 2D area-preserving map is equivalent to a Hamiltonian system with 1 degree of freedom and a potential energy that is periodic in time.Symplectic maps arise as effective models in a wide range of domains, such as celestial mechanics, plasma hydrodynamics, fluid flows, and particle accelerators [3].
Although area-preserving maps of the plane have been studied for decades, only a handful of integrable cases have been discovered.The existence of integrals of motion restricts phase space dynamics to a set of nested tori.Any discovery of novel integrable maps is highly non-trivial and historically has been accomplished by scientists making inspired guesses and applying relevant analytical tools.In many cases, even the proof of the integrability of a specific mapping is quite challenging.
An automated search for integrable systems remains a nontrivial task: only a few successful examples exist to date where machine-assisted searches resulted in new discoveries in nonlinear dynamics.For example, computer-assisted searches combined with high precision simulations found new families of periodic orbits in a gravitational three-body problem [4][5][6].Recent works utilized machine-learning methods to discover integrals of motion of classical dynamical systems from phase space trajectories, e.g.Refs.[7][8][9][10].These examples employed a computer-assisted search to either discover new families of trajectories or to rediscover previously known invariants of motion in integrable systems.In contrast, our algorithm has found previously unknown integrable symplectic maps, advancing the theory of integrable systems.
In this paper, we propose and describe this new algorithm for discovering area-preserving integrable maps of the plane.While scanning over a large space of potential candidate maps, our algorithm traces out individual orbits and subsequently analyzes shapes of invariant manifolds.We restricted our search to the family of mappings which possess integrals of motion with a polygon geometry, thus constraining our search to piecewise linear transformations of the plane.Piecewise linear mappings could be regarded as the simplest possible dynamical systems with rich nonlinear dynamics, which makes them an attractive model to study.They could be viewed as fundamental building blocks of generic integrable symplectic maps.An arbitrary smooth 2D symplectic map could be approximated as an n-piece piecewise linear map, but only in the uniform sense, see Section VIII C. By applying a machine-assisted search, we found over a hundred new integrable mappings; some of them are isolated in the space of parameters, and some of them are families with one parameter (or the ratio of parameters) being continuous (∈ R (+) ) or discrete (∈ Z (+) ).
The structure of this article is as follows.In Section II we provide a historical overview of integrability of symplectic mappings of the plane.In Sections III and IV we review linear integrable mappings produced by 1-and 2-piecewise linear force functions.At the end of Section IV we introduce nonlinear integrable mappings with polygon invariants, and in Section V we propose an algorithm for the machine-assisted discovery of integrable systems with more complicated force functions.Sections VI and VII presents the main results of our search, and finally, in Section VIII we discuss some practical applications of mappings with polygon invariants.Several appendices contain supplementary materials, tables with dynamical properties of mappings, and a pseudo-code for the search algorithm.

II. 2D SYMPLECTIC MAPPINGS AND INTEGRABILITY: HISTORICAL OVERVIEW
The most general form of an autonomous 2D mapping of the plane, M : R 2 → R 2 , can be written as q ′ = F (q, p), p ′ = G(q, p), where ′ denotes an application of the map, and F and G are functions of both variables.If the Jacobian of a 2D map, J ≡ ∂ q F ∂ p G − ∂ p F ∂ q G, is equal to unity, then the map is area-preserving and thus also symplectic (for higher dimensional mappings, the volume-preservation is necessary but not sufficient for a map to be symplectic).A map is called integrable in a Liouville sense [11] if there exists a real-valued continuous function K(p, q) called an integral of motion (or invariant), for which the level sets K(p, q) = const are curves (or sets of points) and The first nonlinear symplectic integrable map of the plane was discovered by E. McMillan [12] when he studied a very specific form of the map where f (p) will be called the force function.We will refer to this mapping as the McMillan-Hénon form or the MH form for short (please do not confuse with McMillan integrable map [12], which is a special case when f (p) is given by a ratio of two quadratic polynomials).At first, the choice of the MH form (1) might look very restrictive and unmotivated.McMillan's original idea was based on its simplicity and clear symmetries (which will be described in detail below), which still allowed some degree of freedom in modifying the dynamics by varying f (p).As it was demonstrated later by D. Turaev [13], almost every symplectic map of the plane (and even 2nD symplectic map) can be approximated by iterations of MH maps.In this article, we will restrict our consideration to this form; however, our results can be easily generalized to different representations of the map.
Below we summarize several important results regarding the integrability of mappings in a McMillan-Hénon form: 1. First, the family of integrable mappings of the plane in the MH form were discovered by McMillan [12] for a biquadratic invariant in the form (I) : K(p, q) = A p 2 q 2 + B (p 2 q + p q 2 ) + Γ (p 2 + q 2 ) + ∆ p q + E (p + q) corresponding to the force function and integrable for any values of parameters A, B, Γ, ∆ and E (Fig. 1a).
2. Later, Y. Suris [14] studied the integrability of a difference equation in the form x n+1 + x n−1 = f (x n ) referred to as the Suris form of the map (sometimes referred to as the standard type).The Suris form of the mapping has a clear physical interpretation, since in the continuous-time limit it takes a form of a 1D Newton equation ẍ = f (x) − 2 x, where the right hand side plays a role of a mechanical force.Mapping in the Suris form is just another representation of the MH form by setting q n = x n and p n = x n+1 (one can think of an analogy between the Lagrangian approach with a single differential equation of the second order, and the Hamiltonian approach with two differential equations of the first order, e.g.see [11]).Suris proved that for analytic force functions f (p) and K(p, q), the invariant of the integrable mapping can take only one of the forms: (I) biquadratic function of p and q (McMillan map), (II) biquadratic exponential or (III) trigonometric polynomial: (II) : K(p, q) = A e α p e α q + B (e α p + e α q ) + Γ (e α p e −α q + e −α p e α q ) + ∆ (e −α p + e −α q ) + E e −α p e −α q , (III) : and Mappings (I -III) are integrable for any choice of parameters, see some examples with invariant level sets in Fig. 1

In 1983 Morton
Brown proposed an interesting problem in the section "Advanced Problems" of Amer.Math.Monthly [15].
6439.Let {x n } be a sequence of real numbers satisfying the relation One can notice that this map is in the Suris form (so it can be written in the MH form as well) with f (p) = |p|.
In 1985 the solution by J. F. Slifker was published in [16] with a comment below: Also solved by the proposer and sixty-one others.The problem turned out to be rather elementary.
Despite its "simplicity," the map has proven to be interesting for dynamical system theory, combinatorics and topology.10 years later, in 1993, Brown extended his results to the formal publication [17].He reconsidered the equation as the homeomorphism of the plane and provided a geometrical proof of periodicity including a description of the polygon invariant of motion (see Fig. 1d): D. Knuth, who provided his own combinatorical proof of periodicity, wrote a letter to Brown: When I saw advanced problem 6439, I couldn't believe that it was "advanced": a result like that has to be either false or elementary!But I soon found that it wasn't trivial.There is a simple proof, yet I can't figure out how on earth anybody would discover such a remarkable result.Nor have I discovered any similar recurrence relations having the same property.So in a sense I have no idea how to solve the problem properly.Is there an "insightful" proof, or is the result simply true by chance?
While quite often this map is referred to as Knuth map, we would like to call it Brown-Knuth map to acknowledge all contributors.
4. Finally, in his original article [12], Edwin McMillan suggested another intricate integrable system with a piecewise quadratic invariant and a piecewise linear force function where a ∈ R is an arbitrary coefficient.He referred to the corresponding invariant level sets as beheaded and twoheaded ellipses (Fig. 1e).As we will see in Section IV A 2, this map as well as Brown-Knuth map are examples of more general integrable systems with the force function f (p) = a p + b |p| and invariants being collection of ellipses, hyperbolas and straight lines (segments).

5.
At the end, we should mention a special class of finite difference equations, which are integrable (in the sense that they admit a Lax pair representation), known as discrete Painlevé equations.Discrete Painlevé equations are closely related to integrable symplectic maps of the plane.We discuss the connection in Section VIII A.
Several properties of MH mappings (1) which will help us to understand and interpret further results are: • A map in the MH form is symplectic and reversible for any continuous f (p), with the inverse • For any map with f (p) there is a twin map with the mirrored force function −f (−p), so that their dynamics are identical up to a rotation of phase space by an angle of π, (p, q) → −(p, q).Thus, we will omit twin maps unless it is necessary.
• A MH map can be decomposed into two transformations M = F • Rot(π/2), where Rot(π/2) : q ′ = p, p ′ = −q, and F : with Rot(θ) being a rotation (about the origin by an angle θ): and F being a thin lens transformation.The name "thin lens" comes from optics and reflects that the transformation is localized so coordinates remain unchanged, q ′ = q, and momentum (or angle in optics) is changed according to location, ∆p ≡ p ′ − p = f (q).Both transformations are area-preserving and symplectic.This decomposition is very important and explains the physical nature of the map representing a particle kicked periodically in time with the force f (q).A general model of kicked oscillators includes famous examples of Chirikov map [18,19], Hénon map [20] and many others.In particular, this decomposition is often employed in accelerator physics; a model accelerator with one degree of freedom consisting of a linear optics insert and a thin nonlinear lens (such as sextupole or McMillan integrable lens) can be rewritten in the MH form.
• A MH map can be decomposed into superposition of two transformations M = G • Ref(π/4) where Ref(π/4) : q ′ = p, p ′ = q, and G : with Ref(θ) being a reflection about a line passing through the origin at an angle θ with the q-axis Ref(θ) : and G being a special nonlinear vertical reflection.Transformation G leaves the horizontal coordinate invariant, q ′ = q, while vertically a point is reflected with respect to the line p = f (q)/2, i.e. equidistant condition p ′ − f (q)/2 = f (q/2) − p is satisfied.Both reflections are anti-area-preserving (the determinant of Jacobian is equal to −1) and are involutions, i.e. double application of the map is identical to transformation Ref 2 (θ) = G 2 = I 2 .In addition, each reflection has a line of fixed points: p = q for Ref(π/4) (further referred to as the first symmetry line) and p = f (q)/2 for G (further referred to as the second symmetry line), see Fig. 2.
• Fixed points of the map are given by the intersection of the two symmetry lines and always in I or III quadrants: • 2-cycles (if any) are given by the intersection of the second symmetry line p = f (q)/2 with its inverse q = f (p)/2 (so the intersections is always in II and IV quadrants).Additionally, if f (p) is an odd function, 2-cycles are restricted to anti-diagonal p = −q.
• Finally, the last and most important property we would like to list here: any constant invariant level set of an integrable MH map (either a point, set of points, closed and open curves, or even collection of closed curves representing islands) is invariant under both reflections, Ref(π/4) and G.If f (p) is an odd function, invariant level sets possess additional symmetry K(p, q) = K(−p, −q), meaning that the map coincides with its twin.
Invariance under these transformations results in two geometrical consequences, first noticed in relation to integrability by E. McMillan [12].First, invariant level sets are symmetric with respect to transformation Ref(π/4): Second, nontrivial symmetry with respect to transformation G is In order to understand how this symmetry works, let's consider a closed invariant curve K(p, q) = const encompassing the fixed point (blue curve in Fig. 2).Solving implicit equation K(p, q) = const for p results in two roots p = Φ(q) or p = Φ −1 (q), which correspond to the upper or lower branches in Fig. 2. Next, using the mapping equation (MH form of map) the reader can verify that the following important condition holds ( The condition (2) implies that the symmetry line p = f (q)/2 "splits" closed curve K(p, q) = const into two halves: upper Φ(q) and lower Φ −1 (q) branches with the endpoints located at ∂ q Φ(q) = ±∞.For each value of q where K(p, q) = const is defined, the upper and lower parts are equidistant from the second symmetry line McMillan wrote in regards to this condition: "This is a remarkable result, of startling simplicity, which fell out almost without effort on my part.It leads to not just one, but to great families of functions f (p) giving regions of stability.It also has a simple geometrical interpretation." Here we would like to make two comments.(i) The origin of these two symmetries is due to reversibility of symplectic maps [21]; any symplectic map with an inverse can be decomposed into two involutions similar to a map in the MH form with the integral of motion being invariant under both transformations.(ii) Both symmetries also hold for survived invariant tori in chaotic maps.Fig. 2 illustrates symmetry properties for two classes of survived tori in Hénon map: a closed curve and a set of islands.
2. Two constant level sets K(p, q) = const for the chaotic Hénon map, f (p) = a p + b p 2 .The black curve shows a closed trajectory encompassing the origin and the gray curves show a set of islands (i.e.under iterations, points hop from island to island covering the closed gray curves labeled i1 -i5).While the first (reflection) symmetry is quite obvious from the figure, we will focus on the second one.The upper branch of the blue curve, Φ(q), is vertically equidistant from the lower branch, Φ −1 (q).Island i3 satisfies the second symmetry in the same manner as does the blue curve.Islands i1 and i2 satisfy the symmetry in a sense that they are vertical reflections of islands i5 and i4, respectively (e.g. the lower part of i1 is equidistant from the upper part of i5 and the upper part of i1 is equidistant from the lower part of i5).The dashed and solid green lines are the first and second symmetry lines.

III. 1-PIECE MAPS
In this section, we will consider mappings in the MH form with f being a simple linear function, f (q) = k q.We will refer to such mappings as 1-piece maps or mappings with linear force function, purposely avoiding the term "linear map."We will reserve the term linear map (in contrast to nonlinear map) in order to refer to independence of dynamics on amplitude (i.e.rotation number is independent on amplitude for all stable trajectories); mappings in the MH form with nonlinear force function can possess linear dynamics.In order to avoid confusion, we will follow the terminology established above.As we will see, all 1-piece maps are integrable and linear.

A. Polygon maps with integer coefficients
According to crystallographic restriction theorem [22,23], if A is an integer 2×2 matrix and A n = I 2 for some natural n ∈ N, where I 2 = diag (1, 1) is an identity matrix, then the only possible solutions have periods n = 1, 2, 3, 4, 6, which corresponds to 1-,2-,3-,4-and 6-fold rotational symmetries.Transformations with a period n = 1, 2 are simply ±I 2 , which can be considered as special cases of rotation of the plane, Rot(θ), through the angles equal to θ = 0 or π respectively.Three other cases, n = 3, 4, 6, are given by mappings in the MH form further referred to as α, β and γ respectively, and, with invariant polygons being concentric triangles, squares and hexagons (see Fig. 3).3. Integer mappings with linear force function and polygon invariants α, β and γ (we use the same naming convention as Ref. [24,25]).Note that the invariant triangle was modified in order to satisfy both symmetries.Each figure shows one invariant polygon, with each colored region i being mapped to region i + 1.Here n is the period and ν is the rotation number of the map, which represents the average increase in the angle per map iteration.The dashed and solid dark green lines show the first and second symmetry lines.

B. Fundamental polygons of the first kind
In order to understand the nature of invariant polygons from crystallographic restriction theorem, we first will consider a more general 1-piece map in the MH form where parameter k ∈ R is not restricted to integers.This map is known to be integrable with quadratic invariant of motion [12] K(p, q) = p 2 − k p q + q 2 .Mapping (3) is unstable for |k| > 2, that results in level sets of the invariant being hyperbolas.On the other hand, the mapping is stable for |k| < 2 with level sets corresponding to a family of concentric ellipses, see Fig. 4. Within the parameter range corresponding to stable maps, the rotation number ν is independent of the amplitude and given by ν = arccos(k/2)/(2π).
Constant level sets of invariant for 1-piece map in the MH form, K(p, q) = p 2 − k p q + q 2 = const .For |k| < 2 the level sets are ellipses, for |k| > 2 the level sets are hyperbolas and for |k| = 2 they are parallel lines.The dashed and solid green lines are the first and second symmetry lines.
Thus, one can see that if arccos(k/2) is incommensurate with π, the motion is quasi-periodic, and according to the Poincaré recurrence theorem, iterations starting from the initial condition (q 0 , p 0 ) will result in points densely covering the ellipse K(p, q) = K(p 0 , q 0 ).However, if arccos(k/2) is commensurate with π, motion is strictly periodic and ν becomes a rational number.As a result, point (q 0 , p 0 ) will return to its initial position after a finite number of iterations and will not depict an ellipse.In fact, linear mappings with rational ν are more than integrable, they posses superintegrability, meaning that they have more than one functionally independent invariant of motion K(p, q).
For example, for each map with rational ν ∈ Q, one can construct two sets of polygon invariants by iterating an initial point starting on the first, p = q, or second, p = k q/2, symmetry lines and using iterations as vertices (Fig. 5).We will refer to such polygons as first and second sets of fundamental polygons of the first kind; first or second sets refer to the position of the initial point on first or second symmetry lines respectively, and "the first kind" refer to the fact that the force function consists of one piece only.In the case of even period, two sets of polygons are actually different while for odd period sets, polygons are identical up to rotation by an angle π, (q, p) → (−q, −p).In the former case, the first set of polygons has no vertices on the second symmetry line (and two vertices on the first symmetry line) while the second set has no vertices on the first symmetry line (and two vertices on the second).In the latter case, both sets of polygons have a single vertex on each symmetry line.
In addition, for |k| = 2 invariant level sets are straight lines.We will denote these mappings as "a" for k = 2 and "b" for k = −2.While being unstable, both maps have polygonal chains (connected series of line segments, in this case just line itself) as its invariant.In the following sections, we will see that unstable and stable polygon maps can be combined together to form more complicated systems.Each plot shows the invariant ellipse K(p, q) = p 2 − k p q + q 2 = const (black) and two invariant fundamental polygons inscribed into it (blue and red).The dashed and solid green lines correspond to the first and second symmetry lines.

Polygon maps with integer coefficients
Another remarkably interesting result is given by CNR Theorem (Cairns, Nikolayevsky and Rossiter [24,25]).Suppose that M is a periodic continuous mapping of the plane that is a linear transformation with integer coefficients in each half plane q ≥ 0 and q < 0. Then M has a period n = 1, 2, 3, 4, 5, 6, 7, 8, 9 or 12.  [24]).Each figure shows one invariant polygon, with each colored region i being mapped to the region i + 1.Here n and ν are the period and rotation number of the map.The solid and dashed dark green lines show the second and first symmetry lines.
All maps discovered by Cairns and others are in the MH form with force function consisting of two linear segments In addition to previously described 1-piece maps (trivial cases k 1 = k 2 ) with n = 1, 2, 3, 4, 6, now we have 5 more mappings with n = 5, 7, 8, 9, Fig. 6.These maps are stable, linear and integrable, with polygons being their invariant sets, and denoted with capital Latin letters D -H.Note how Cairns, Nikolayevsky and Rossiter rediscovered Knuth map, H, from the first principles.

Fundamental polygons of the second kind
Similar to the case of 1-piece mappings, in order to understand the origin of polygon invariants we should consider a less restricted map, corresponding to force function (4): where parameters k 1,2 ∈ R.This is a very nontrivial dynamical system considered in great details in [26][27][28][29].Here, we will provide some results contributing to our understanding of polygon maps.
• Map ( 5) is linear, i.e., there is no dependence of the rotation number on the amplitude.Since we are free to choose units of p and q, using scaling transformation (ϵ > 0) one can reconstruct dynamics on any ray of initial conditions (starting at the origin) from a single trajectory.
• The invariant of motion is a combination of segments of ellipses, straight lines or hyperbolas "glued" together and defined by C 1 p 2 + C 2 p q + C 3 q 2 = const.In the plane of parameters (k 1 , k 2 ) there are lines of constant topology.Along these lines, the number of segments remains constant.
• Dynamics of map ( 5) can be either stable or unstable depending on the values of k 1,2 , that creates a nontrivial fractal area on the stability diagram (see Ref. [26] for more details).While we are not considering the fine structure of the fractal area, we would like to note that all pairs (k 1 , k 2 ) resulting in stable mappings are within the area defined by (Fig. 7) p, q < 2 ∧ p q < 4 for p, q < 0.
Stable maps are either periodic when the rotation number ν is rational or quasi-periodic for irrational ν.
- 7. The plane of parameters (k1, k2) for 2-piece maps.The dashed lines show the area of global stability, p, q < 2 ∧ p q < 4: the map can be stable (but not necessarily) only within the area of stability, shown in white.The solid curves show all lines of "constant topology" passing through integer nodes (k1, k2) ∈ Z 2 within the stable zone: the red line corresponds to a single ellipse k1 = k2, blue lines corresponds to two ellipses (McMillan beheaded and two-headed ellipse) k1,2 = 0, green lines refer to three ellipses defined by k1,2 = 1 and purple and brown hyperbolas are k1k2 = 2 and k1k2 = 3 respectively.Integer nodes corresponding to stable mappings with polygon invariants are shown in black, unstable mappings with an invariant being a polygonal chain in purple, and the remaining integer nodes within the global region of stability are shown with crosses (unstable maps with invariant consisting of hyperbolas).
• For example, the simplest case corresponds to the family of 1-piece maps, k 1 = k 2 , (red line in Fig. 7).As we learned above, in this case the phase space consists of only one ellipse (or two sets of hyperbolas conjugate to each other) defined by K(p, q) = p 2 − k 1 p q + q 2 .As we move along the line k 1 = k 2 , ellipses bifurcate into hyperbolas at k 1 = |2| and trajectories lose stability.
Next line of constant topology is given by k 2 = 0 and line k 1 = 0 corresponding to family of twin maps, see blue solid lines in Fig. 7.This is the aforementioned beheaded and two-headed ellipses discovered by McMillan with invariant defined by two ellipses glued together: Here is his explanation on how this map works: In an earlier paragraph a promise was made that cases could be constructed with boundaries that more than double valued.I shall now fulfill the promise by describing a procedure by which this can be done.Take any known case with a center of symmetry at the origin, and erase the part of the boundary lying in the +, + quadrant.Replace f (q) by f (q) = 0 to the right of the origin.Fill in the erased part of the boundary by reflecting the part in the +, − quadrant about the q-axis.The resulting completed boundary is an invariant under the transformation with the original f (q) to the left of the boundary, and f (q) = 0 to the right.Starting with ellipses, we can get beheaded ellipses (first row of Fig. 8, case 0 or double headed ellipses (first row of Fig. 8, case −2 < k 1 < 0), and here we see a four-valued function acting as an invariant boundary.Note that in order to have only two ellipses being invariants of motion, the second ellipse have to be added to the +, + quadrant (or to −, − quadrant for twin maps with k 1 = 0).Otherwise, an addition of an ellipse into the +, − quadrant will result in its reflection in the quadrant −, + (due to the first symmetry line p = q) making number of ellipses equal to three.
In the case when the second segment of the force function has zero slope, k 2 = 0, the map in Eq. ( 5) is stable for |k 1 | < 2 (see first row of Fig. 8).For the slopes k 1 = ±2 the invariant level sets become polygonal chains with 1-and 3-fold symmetries respectively.When stable, the rotation number is given by As in the case k 1 = k 2 , when arccos(k/2) is commensurate with π, we have a rational rotation number ν ∈ Q resulting in a periodic map.The only integer values of k 1 satisfying this condition are k 1 ∈ {−1, 0, 1} (mappings F, β and G respectively).
FIG. 8. Invariant level sets along lines of constant topology for a 2-piece map in the MH form.The first and second rows illustrate the transformation of the invariant for cases k2 = 0 and k2 = 1 (blue and green lines in Fig. 7).The three left plots in the third row show bifurcation along hyperbola k1k2 = 2 and the three right plots are for k1k2 = 3 (purple and brown hyperbolas in Fig. 7).The solid and dashed green lines are the second and first symmetry lines.
In fact, this line of constant topology is just one of many defined by k 2 = 2 cos(π/n) for n ≥ 2. On each line the map is stable for |k 1 | < 2, it produces a polygonal chain as invariant for k 1 = ±2 and when stable has a rotation number Along each line of constant topology the invariant consists of n ellipses (hyperbolas) with n − 1 of them being in the first quarter.
Among these lines there is another line that we should focus on: case n = 3 with k 2 = 1 (and twin k 1 = 1).As we can see, in addition to line k 2 = 0, this is the only line with integer k 1 .The invariant of motion is second row of Fig. 8. From this line we have two more mappings with polygonal chain (k 1 = ±2) and three periodic maps with integer coefficients (k 1 = −1, 0, 1, i.e maps α F and H).
• As one can notice, lines k 1 = k 2 and k 1,2 = 0, 1 are covering most of integer nodes (k 1 , k 2 ) ∈ Z 2 within the area of global stability (see Fig. 7).The remaining nodes belongs to two lines of constant topology from the family given by k 1 k 2 = 4 cos 2 (π/(2n)) with k 1,2 < 0 and n ≥ 2. Those lines are for n = 2 and 3. Along those lines map is always stable and periodic with the rotation number ν = (2 n − 1)/(4 n).The invariant of motion is glued of (2 n − 1) ellipse: the central ellipse, (n − 1) ellipses in II quadrant, and another (n − 1) in IV quadrant (placed symmetrically with respect to p = q).For example, for n = 2 the invariant is Third row of Fig. 8 illustrates this case.For k 1,2 = −1 from lines n = 2, 3 we have two remaining periodic maps with integer coefficients, E and D.
Thus, similarly to the case of the 1-piece map (degenerate line k 1 = k 2 ), we can define fundamental polygons of the second kind: for each line of constant topology, its stable maps with a rational rotation number are degenerate and have more than one invariant of motion including polygons.In contrast to the case k 1 = k 2 , when we had two sets of polygons, now we only have only a single set of fundamental polygons.These polygons have two vertices on the second symmetry line for maps with an even period or a single vertex on both symmetry lines for odd periods.Fig. 9 shows a few examples of this degeneracy illustrating fundamental polygons of the second kind inscribed in alternative invariant set of ellipses.Note that the first example is for Brown-Knuth map showing an alternative three-headed ellipse which is also an invariant of motion.

Ingredients of mappings with polygon invariant
Hereby, one can see that there is one condition and two mechanisms determining polygon shapes (including unstable mappings with invariant level set being connected series of line segments -non-intersecting polygonal chains).The condition is that the force function f should be in a piecewise linear form that follows from the second symmetry: the sum of vertically opposed polygon sides is proportional to the value of the force function.Two mechanisms are: • Superintegrability in linear maps.Rational rotation number causes degeneracy when more than one invariant of motion exists, including polygons (e.g.fundamental polygons of the first and second kind).
• Loss of stability.Ellipses are transitioned to hyperbolas via the straight lines in linear 2-piece maps, e.g.see cases |k 1 | = 2 for first two rows in Fig. 8.
These mechanisms are respectively responsible for stable and unstable polygons.As we will see further, more complicated nonlinear mappings with polygon invariants can be constructed, but at lower ((q, p) → 0) or higher amplitudes ((q, p) → ∞), they reduce to one of the already described transformations; we will observe that both mechanisms can be combined in one map resulting in very nontrivial dynamics.Below, in Fig. 10 we illustrate all integer linear mappings with 2-piece force function and polygon invariants (stable and unstable). .Integer linear mappings with 2-piece force function and polygon invariants (stable and unstable maps).The black polygons (or polygonal chains for unstable mappings) show constant level sets of the invariant, red level sets are separatrix for unstable mappings with stable periodic dynamics, and the green line is the second symmetry line.

B. Nonlinear mappings
Here we begin our exploration of nonlinear symplectic mappings.So far we considered only mappings with force function in the form (4).The simplest modification we can make is to add a constant shift parameter, d: The map becomes nonlinear and we can not reconstruct dynamics along the ray of initial conditions using scaling transformation (q, p) → ϵ(q, p), ϵ > 0. Instead, now we have a choice of natural units ϵ = d, and one can show that without loss of generality for 2-piece maps, d can be restricted to ±1 or 0. Introduction of d not only adds the dependence on amplitude, but might result in a loss of integrability causing chaotic behavior.A famous example is the Gingerbreadman [30] map -a chaotic two-dimensional map in the MH form with f (q) = |q| + 1 (see first plot in Fig. 11).The remarkable property of this map is that some invariant tori survived perturbation and after deformation remained polygons.The chaotic dynamics is sectioned by concentric polygons making this map the simplest model for area-preserving twist mappings with zones of instability.It shares many properties of the smooth mappings and has an advantage of having an exact arithmetic on rational numbers.Another mapping we would like to mention here is the Rabbit map for f (q) = |q| − 1.To the best of our knowledge this map has not been mentioned before, but we are certain that Robert Devaney, the discoverer of Gingerbreadman, knew about it due to the close connection between the maps (second plot in Fig. 11).

Zoo mappings
After reading the article of Cairns, Nikolayevsky and Rossiter [24], our first idea, which predetermined discoveries made in this article, was to produce objects similar to Devaney map using newly discovered mappings D -G.First we tried mapping D, and amazingly it worked!We introduce two new dynamical systems, the Octopus and the Crab mappings for f (q) = |q|−2q ±1 respectively [31] (last two plots in Fig. 11).Informally, we grouped them together with Rabbit and Gingerbreadman, and, called them zoo maps.Mappings are again chaotic and "sectioned" by infinitely many concentric invariant polygons.Basic dynamics of these maps is described in the caption to Figure 11  Within each zone of instability there are three different scenarios.First, orbits that result from initial conditions (q0, p0) within any gray zone such that both q0, p0 ∈ Q.Those are stable orbits with an exact periodic motion and rational rotation number resulting in an invariant set of points.The second scenario includes all other initial conditions within the gray zones with any or both q0, p0 being irrational, q, p0 ∈ R \ Q.Those are chaotic orbits never returning to initial conditions and densely covering the zone of instability.The third and last scenario is for any initial condition within the colored regions of phase space.Those are chains of linear islands; the initial condition from these regions of phase space hops from island to island of the same group (each group is shown with its own color related to the relative orbit period).Such initial condition returns to itself and the motion is strictly periodic with the rotation number being mode-locked within the island.To view an "animal" you must rotate each image 135 degrees clockwise.
Dynamical properties of integer nonlinear integrable mappings with polygon invariants and 2-piece force function.
Here ν0 is the rotation number in the inner linear layer, ν1, J1 and J ′ 1 are respectively rotation number, action and partial action in the outer nonlinear layer, x is the linear amplitude in a nonlinear layer defined as horizontal or vertical distance from separatrix to the trajectory under consideration.

Nonlinear integrable maps with polygon invariants
When d = ±1, not all mappings experience chaotic behavior like Gingerbreadman or zoo maps.Unexpectedly, transformations E, F and G produced nonlinear integrable systems with polygon invariants if the vertical shift parameter d is added [31] (see Fig. 12).Since d ̸ = 0, fixed point moves from the vertex to one of the linear pieces (depending on a sign of d) of the force function.As a result, the area of mode-locked motion around fixed point is formed with invariant level sets being triangles, squares, hexagons or line segment p = −q with period 2 for map b0.For larger amplitudes, the motion is still integrable, but rotation number becomes amplitude dependent (trajectories shown in blue).The inner linear and outer nonlinear layers are separated by separatrix matching dynamics (shown in red).Dynamical properties for all mappings are listed in Table I as functions of amplitude along q-axis x, and, methods for calculation of rotation number as well as action-angle variables are provided in Appendix A.  10); note that the invariant level set of the outer nonlinear layer at the large amplitude is asymptotic to the level set of the corresponding map.The black (blue) level sets are linear (nonlinear) trajectories.The red level sets represent the separatrix isolating the inner linear from the outer nonlinear layers.The green line shows the second symmetry line, p = f (q/2).

V. MACHINE-ASSISTED DISCOVERY OF INTEGRABLE MAPS WITH POLYGON INVARIANTS
After discovering 6 nonlinear polygon maps shown in Fig. 12 and performing few more numerical experiments, we realized that the force function can be further generalized to have n > 2 linear segments, potentially producing new families of integrable systems.In order to automate the search, we designed a relatively simple algorithm that finds integrable mappings based on the analysis of individual trajectories; first, we verified some of its principles in [32], and then we significantly modified it and performed an extensive search.The detailed pseudo-code of the search algorithm is presented in Appendix B while its Python implementation is available at github.com/yourball/polygon-maps.The algorithm is illustrated via flowchart in Fig. 13 and its blocks with processes/decisions are described below.0. Input We decided to focus on piecewise linear forces with integer coefficients, f (q, k, l, d) with integer vectors k, l and d.While we performed a scan only for integer slopes k, we were able to analytically generalize some of the results for more general values of the shift parameter d ∈ R and lengths of segments l ∈ R + .
1. Fixed point In order for a map to be stable, i.e. all phase space orbits are bounded, it should have at least one fixed point ζ * = (q * , p * ).Mappings that do not have a fixed point, ∀ q ∈ R : q ̸ = f (q)/2, are out of scope.

Stability
Mappings with fixed point still can have unstable trajectories at larger amplitudes.If for some initial condition we detect an unbounded growth of the radial distance from the fixed point, ∃ i : |ζ i − ζ * | > r max , the map is excluded from further analysis.

4.
Ordering At the next stage the points of each individual trajectory are ordered according to their polar angle ϕ = arctan [(p − p * )/(q − q * )], such that the angle ϕ is monotonically increasing.After arranging trajectory points according to ϕ we join consecutive points resulting in an orbit polygon.

Vertex count
In order to count the number of vertices, V , for each orbit polygon we use scalar products of two vectors build on three consequent points (v where θ(x) = 1 x>0 is a Heaviside step function.In numeric experiments, when v 1 ∥ v 2 , the normalized scalar product might be slightly different from unity, so a small cutoff parameter should be introduced, ϵ < 6. Vertex cutoff At the final stage we are left with stable integrable or chaotic maps which we would like to separate.The key insight is that we are looking for transformations with all invariant tori being polygons and assume that the phase space consists of concentric layers; in each layer the orbit polygons have the same finite number of vertexes which remains bounded when increasing the number of iterations N → ∞.Conversely, for a chaotic trajectory the number of vertices linearly grows with the number of iterations, V = O(N ).Such drastic difference in the growth of V with the iteration number allows us to automatically separate between the two cases.In numeric experiment we check that for all initial conditions the number of vertices is bounded by some cutoff parameter, V < V cutoff ≤ N .
In addition, for each mapping found by the search algorithm, we manually verify the McMillan integrability condition, thus rigorously proving their integrability.Although the search algorithm described here is quite powerful, there are some caveats to bear in mind: • Maps with strongly non-convex layers of polygon invariants can be missed by the algorithm.In cases when the orbit polygon is not convex with respect to the fixed point (q * , p * ), the points of the trajectory can not be ordered according to the polar angle.While we do not have a particular example, we should keep it in mind.
• Inside integrable nonlinear layers fibrated with polygons there are two types of orbits: those with rational and those with irrational rotation numbers.The latter are quasi-periodic with entire polygon being densely covered.On the other hand, the former are strictly periodic with orbit visiting only finite number of points, thus potentially showing a smaller number of vertexes.In Appendix A we show that all rotation numbers within a layer are in the form ν(q 0 ) = (α+β q 0 )/(γ +δ q 0 ), where α, β, γ and δ are integer parameters.In order to avoid periodic orbits which complicate analysis, we can add a small irrational perturbation η to all initial conditions, q 0 → q 0 + η.In fact, from this formula one can see that in numeric experiments all observed orbits are periodic but this is not a problem as long as period is sufficiently large compare to the number of polygon sides, so points of orbit visit all sides of polygon covering them densely enough.13.Flowchart representing the machine-assisted discovery of integrable maps.The algorithm verifies if the map defined by the input (0.) is integrable with the invariants being concentric polygons.If the map has a fixed point (1.), the algorithm performs tracking for different initial conditions (2.), and, if all orbits are bounded (3.) it proceeds to the construction of orbit polygons via ordering (4.).Finally, after counting the number of vertices V (5.), if for each polygon V < V cutoff , the map is sent for analytic verification.The search process terminates when any conditional operation is false (NO).

Yes
• Finally, while we've been only looking for maps with invariant layers being concentric polygons, while phase space can involve polygon islands.In the case of completely mode-locked islands (linear islands) with all orbits visiting each island only once, our algorithm will successfully handle this case.But in the case of nonlinear islands (i.e. a trajectory visits sequentially each island under iterations and traces out invariant tori around centers of islands) our algorithm indeed will fail; as in the case of non-convex polygons, points of orbit can not be arranged by polar angle since now they are a multiple valued function.
We performed the search for 3-and 4-piece integer force functions using our algorithm across a wide range of the parameter values.The integer valued slopes of the segments of the force function are selected from the range −10 ≤ k i ≤ 10 and the shift parameter −50 ≤ d ≤ 50.In the case of four-piece maps, the additional degree of freedom in the space of parameters corresponds to the ratio of the lengths of segments, r ∈ [1/10, 1/9, . . ., 1, 2, . . ., 10].However, the search space can be further reduced when relying on our prior knowledge about properties of piecewise linear mappings with two segments.A generic n-piece mapping at small and large amplitudes reduces to a two-piece mapping.As illustrated in Fig. 10 two-piece mappings with integer slope coefficients and polygon invariants are confined to the parameter range −3 ≤ k 1,2 ≤ 2. Therefore the values of the slopes of the outer segments in mappings with polygon invariants are also restricted to −3 ≤ k 1,n ≤ 2. Similarly, in the vicinity of the fixed point (q * , p * ) at small amplitudes |q − q * | → 0 and |p − p * | → 0, the mapping behaves as a single-piece or two-piece map, depending on whether the fixed point belongs to the segment of the force function or coincides with one of the vertexes separating adjacent segments.Therefore, in order to guarantee that the mapping has polygon tori in the vicinity of the fixed point, the range for the slope k i for the corresponding segment is also confined to the same range.The values of the shift parameter d and the length ratio r are less restrictive, because they are only responsible for the position of the fixed point.While we performed a scan only for integer d and r we were able to analytically generalize some of these results for d, r ∈ R.

VI. 3-PIECE MAPS
In order to construct more complex piecewise linear symplectic maps, one can add more segments to the force function f .When the number of segments is equal to three, we have four independent parameters: three slopes k 1,2,3 and a shift parameter d.By performing the rescaling of coordinates, one can see that without loss of generality, the length of the middle segment can always be assumed to be equal to one.In addition, unlike the case of 2-piece force function when d = 0, ±1, now the shift parameter belongs to reals.Below, we will classify search results as isolated with respect to d (finite set of integer values), discrete (all integer values, d ∈ Z), and continuous (map is integrable for any d ∈ R).Results of the machine-assisted search of new integrable maps are briefly summarized in Fig. 14 (the symmetry with respect to k 1 = k 3 is due to "twin" mappings).

A. Isolated d
Integer integrable polygon mappings with 3-piece force function and isolated d are shown in Fig. 15 and rotation numbers are provided in Table II.As one can see, the addition of new pieces can stabilize previously unstable 2 pieces maps a (a1,2 * ), b, c and d (respectively mappings a1-3, b1-3, c1-2 and d1-2).Mappings α3, β4-5 and H1 can be seen as previously considered linear transformations (α, β and H) with outer layer of nonlinear stable trajectories determined by the new piece in the force function.
Something interesting can be observed for mappings b0.1, β3 and G1.For the map b0.1 we see that the inner and outer nonlinear layers are separated by a chain of islands.Unlike the "usual" chaotic chains observed in systems like Hènon [20] or Chirikov [18,19] mapping (where under iterations phase-space point jumps from island to island and trace out elliptic orbits around centers of islands), these are linear islands.They show true mode-locking similar to one observed in Arnold tongue [33] or even more similar to Gingerbread man map [30]; when iterated, any initial condition within the chain returns exactly to the same position after visiting each island only once.Thus, the rotation number of the inner layer and the one of the outer layer where x and y represent the linear amplitudes, are matched on the outer and inner separatrices, respectively x = 1 2 and y = 0.They are equal to the rotation number inside the linear chain consisting of 5 islands Black and blue lines correspond to linear and nonlinear invariant phase space trajectories, red lines are separatrices, and the thick green line is a second symmetry line, p = f (q)/2.Mappings b2 and β5 are also integrable for any d ∈ R and mappings a2 and β4 are also integrable for any d ∈ Z. Mappings α4 * and a2.1 * are unstable and were missed by the search algorithm.
Rotation numbers of the inner linear, ν0, and outer nonlinear layers, ν1, for integer integrable polygon mappings with 3-piece force function and isolated d.Last two columns show action, J, and partial action (Hamiltonian), J ′ , in nonlinear layer.Parameter x ∈ [0, ∞) is horizontal (vertical) distance from the separatrix to considered trajectory.
Note that polygons inside the chain of islands satisfy both symmetries as a level set of the entire group and are invariant under transformation; for each island, there is another one which is its vertical reflection with respect to the second symmetry line, and one reflected with respect to the main diagonal (unless the island is its own reflection with respect to the symmetry line).Due to the fact that dynamics inside the island chain is degenerate (all rotation numbers are the same), invariant contours can not be simply traced out in the numerical experiment by looking at the mapping iterations, but can be reconstructed from the symmetries in order to match the separatrices.
Mappings β3 and G1 don't have an inner nonlinear layer, and the chain of islands is attached to a linear layer, so the entire combined phase space is mode-locked.As expected, the number of islands is equal to the period inside linear zone (4 for β3 and 5 for G1) and rotation numbers are 1/4 and 1/5 respectively.
Finally, two more mappings, α4 * and a2.1 * , are shown at the end of the last row.These cases were not found by an algorithm, since the outer layer is unstable, and were constructed by authors using symmetries.
TABLE III.Rotation number, νi, action, Ji, and partial action, J ′ i , for the i-th layer of the integer integrable polygon mappings with 3-piece force function and discrete d. ν isl and #C are the rotation number in mode-locked area and number of chains in group of linear islands, respectively.Mappings a2 and β4 are integrable for any d ∈ Z, forming two families of transformations with discrete parameter.Invariant level sets are shown in Fig. 16.The rotation number is provided in Table III and Fig. 17 shows it as a function of amplitude.Inner and outer layers become separated with chains of linear islands, with increment (or decrement) of d from cases a2 and β4, and the number of chains increases linearly with d.On these groups of linear islands, the rotation number is again mode-locked.
Rotation number ν as a function of amplitude x (horizontal or vertical distance measured from fixed point to invariant polygon) for families of mappings a2-k and β4-k with k = 1, . . ., 4.

C. Continuous d, (C)
In this subsection we present the last class of integer mappings with 3-piece force function which is integrable for d ∈ R. Previously mentioned mappings b2 and β5 as well as two more, γ2 and δ1, are integrable for any real value of the shift parameter.As d is continuously changed, the fixed point is moving along the second symmetry line passing each segment and vertex of the force function.Bifurcation diagrams for mappings' phase space is shown in Fig. 18 and rotation numbers in different layers are given in Table IV.
Something interesting can be observed for the mapping δ1.When d > 3 map has three layers: inner linear area with period 4 and two outer nonlinear layers.As d approaches 3 from above, linear area shrinks to zero and first nonlinear layer "mode-locks" to linear area with rotation number being 2/7; it reflects the bifurcation in the system since fixed point moved to the second vertex of the force function.Decreasing value of d further form 3 to 2, one see how separatrix split forming three layers again with inner layer being period 3.At d = 2 inner nonlinear layer disappears again in another bifurcation.Finally, in region of d in between 2 and 3/2 inner linear area transforms from triangle to hexagon in a very nontrivial way involving appearance of the chain of linear islands (see Fig. 18).The surprising fact is that rotation number in nonlinear layer for d ∈ [3/2, 3] does not depend on the value of d.
x ∈ [0, 1] ν2 1+2y 4+6y d+2y 3d+1+6y 2+2y 7+6y x ∈ [0, 1] ν2 1+2y 6+8y d+y 4d+1+4y 1+y 5+4y Rotation number, νi, for the i-th layer of the integer integrable polygon mappings with 3-piece force function and continuous d.Different columns corresponds to different values (or intervals) of the shift parameter d.For each nonlinear layer, the rotation number is provided along with the interval of amplitude x, y.Amplitude in each layer is measured as a horizontal or vertical distance from the inner separatrix to the invariant curve under consideration.

VII. 4-PIECE MAPS
In this section we will consider force function with four segments.Adding an extra segment results in total of 6 independent parameters defining the force function: 4 slopes k 1,2,3,4 , shift parameter d and ratio r = l 2 /l 3 .The scale transformation allows us to normalize one of two finite pieces (l 2 or l 3 ) to unity so only one of them is independent; alternatively we will use the ratio r when it appears to be more convenient.In next three subsections we will present mappings isolated with respect to d (finite set of integer values) and isolated with respect to r as well, discrete r (all positive integer values, r ∈ N + ) or continuous r (map is integrable for any r ∈ R + ).Results of a search are summarized in Tables V and VI.In addition, more cases which were harder to classify are presented in the last Subsection VII D.
Integrable polygon mappings with 4-piece force function and isolated integer d and r.Such mappings were found only for r = 1, 2, 3 and corresponding values of parameter d (in units of l2) are presented for each quadruplet of integer slopes.4-piece mappings with isolated integer values of d and r are shown in Fig. 19 and action angle variables are provided in Tables VIII and IX; these and all further tables with dynamical properties of 4-piece maps are moved to Appendix C for the reader's convenience.
Naively we expected that layers of phase space change along with the force function and separatrices of motion are associated with vertices of f (q), for example look at mappings b1.1, α3.1, α3.3 or β4.1.But results of our search are showing that dynamics can be more complicated.As one can see for mappings a1.1, a3.2, c1.1 and few others, an additional layer can appear on outer pieces of force function.This is due to the fact that polygons in one of them have sides decreasing with the amplitude; as a result, at some critical amplitude length of the side becomes equal to zero and invariant level sets transition to new layer.We would like to name this phenomena as long-range catastrophe.In this subsection we present systems integrable with any positive integer ratio of r = l 1 /l 2 ∈ N + but isolated integer values of the shift parameter d.Five different families of mappings with 2 prototypes from each are shown in Fig. 20.These maps are similar to 3-piece force function mappings with discrete d (D) (Subsection VI B): they all have chains of linear islands separating outer layers with nonlinear motion, and, the number of chains within the linear island is linearly growing with r.Rotation number ν i for different layers are listed in Tables X and XI.
H1-1 Here we report systems which are integrable with any positive real ratio of r = l 1 /l 2 ∈ R + but still isolated integer values of the shift parameter d.Some of them experience bifurcations of phase space with the change of parameter r (mappings with ′ index) while the rest of systems have same set of layers for any value of r ∈ (0; ∞).Corresponding level sets of invariant and its bifurcation diagrams are shown below in Fig. 21.All dynamical properties are listed in Tables XII and XIII.Mapping α3.1 ′ presents an interesting example: when r = 1 outer otherwise nonlinear layer becomes mode-locked

VIII. APPLICATIONS OF MAPPINGS WITH POLYGON INVARIANTS
In the next three subsections we will discuss the importance and some practical applications of mappings with polygon invariants as well as its connection to other dynamical systems.

A. Connection to Painlevé equations and Suris mappings
Some integrable mappings with polygon invariants have a direct connection to discrete Painlevé equations (dP).Painlevé equations are nonlinear difference equations with canonical form being [34,35]: We will be interested in the autonomous (i.e.mappings do not have explicit dependence on n) limit of dP I−III corresponding to λ → 1.It is worth mentioning that autonomous limit of Eqs. ( 6) could be written in the Suris form by introducing a new set of variables X n = e xn/ϵ , α = e A/ϵ , β = e B/ϵ : dP S I : )(e xn/ϵ + e −A/ϵ ) (e xn/ϵ + e B/ϵ )(e xn/ϵ + e −B/ϵ ) .
While the force function for dP S II is simply a linear function, two other cases are specific examples of integrable exponential Suris mapping.
The ultradiscretized form of Painlevé equations (udP) could be obtained by taking the limit ϵ → 0 in Eqs.(7).At this point we should clarify the terminology: by saying "discrete" Painlevé equation we refer to the difference nature of the equation in contrast to nonlinear second-order ordinary differential Painlevé equations.The addition "ultra" refers to another level of discretization when the force function is transitioning from a smooth continuous to a piecewise linear function.Using identity lim ϵ→0 ϵ log [exp (x/ϵ) + exp (y/ϵ)] = max(x, y) we arrive to udP I : x The udP equations ( 8) have a symplectic Suris form x n+1 + x n−1 = f (x n ) with following ultradiscrete piecewise linear force functions The force function f I (corresponding to the first equation udP I ) is the simplest two-piece linear force without shift parameter, d = 0, and with slopes defined as (k 1 , k 2 ) = (−σ, 1 − σ).Since σ can only take values 0,1 and 2, we obtain previously considered mappings G, F and E respectively (see Fig. 10).The second ultradiscretized Painlevé equation has a trivial linear force function f II and is integrable for any values of parameters A and σ.In the case of integer slopes this results in 3 mappings with stable trajectories [24] α, β, γ and 2 unstable mappings a, b (see Fig. 10).Finally, for the third ultradiscretized Painlevé equation udP III (assuming B > A > 0), we obtain force function containing five pieces.Taking the limit A = const and B → ∞, from Eq. ( 9) we have the linear map α with slope of force function k 1 = −1.Considering another limiting case A → 0, one can see that f (x) could be reduced to one of our 4-piece polygon maps α3.1 ′ (slope coefficients {k 1 , k 2 , k 3 , k 4 } = {0, −1, 1, 0} and d = 2B) with ratio parameter restricted to r = 1 (note that our map α3.1 ′ is integrable for arbitrary r, see Fig. 21). Figure 22 shows invariant level sets for this last case for different values of ϵ.As one can see, approaching the ultradiscrete level ϵ → 0, the force function transitions to piecewise linear and level sets transform to polygons.When ϵ = ∞ map becomes linear with ν = 1/4; note that invariant level sets are chosen as circles and not polygons.As we previously discussed linear mappings with rational rotation numbers are degenerate and allow infinitely many invariants of motion, but in this case the shape is chosen as a limiting case ϵ → ∞ for the exponential Suris invariant.
We would like to emphasize that our polygon mappings represent a much richer family of integrable systems beyond ultradiscretezations of known Painlevé equations udP I-III ; ultradiscretized equations of Painlevé type are limited to only 5-piece force functions and the set of parameters is very restricted (e.g.slopes of outer pieces of 5-piece force functions can only be equal to zero).

B. Near-integrable systems
Since force function is piecewise linear and isn't smooth (not differentiable at discrete set of points), it might seems that these polygon mappings are somewhat impractical.However, in a similar manner to the example in Fig. 22 where the polygon mapping can be smoothed using parameter ϵ, below we will demonstrate how one can construct not integrable, but quasi-integrable systems.One of the possible ways to perform this was first done by Cohen [36,37] in application to Brown-Knuth map (map H) by introducing a new force function as As one can see it tends to |q| as ϵ → 0 or at large amplitudes as q → ∞.Fig. 23 shows a comparison between Cohen and Brown-Knuth maps.Looking at the tracking results one can see that new f (q) produces a system with surprisingly regular behavior.However, detailed numerical experiments reveal multiple island structures at small scales (see Fig. 23 case (c.)) indicating that analytic invariant doesn't exist.Thus, in a similar manner, one can first rewrite a general n-piece linear function using one linear and (n − 1) modulus functions in a form k i |q − q i | where q i is a location of i − th vertex, then substitute all modulus with square roots k Fig. 24 below shows examples of how the "smoothening" procedure can be used to produce new nearly integrable systems.More importantly, in addition to being near integrable, these mappings are stable since limiting polygon at larger amplitudes provides an isolating integral.Note that from a practical standpoint, the fact that produced systems are quasi-and not exactly integrable is not a problem as long as ϵ is sufficiently small and chaotic formations are within needed tolerances.Even if a system is originally designed as integrable, small errors (e.g.rounding errors for simulations, or tolerance errors for physical experimental setups) will destroy integrability on a small scale.Also, the square root function is only one of infinitely many ways of smoothing force function (e.g. one can use arctan, tanh or others) and is used here as an example.FIG.24.Examples of quasi-integrable systems produced by "smoothening" 3-piece integrable polygon maps using ϵ = 0.05.For each plot the label shows parameters {k1, k2, k3; d}.All level sets are obtained by tracking.The green curve is the second symmetry line f (q)/2 and the black dashed curve is the force function f (q).

C. Discrete perturbation theory
Finally, knowledge of all possible integrable mappings with polygon invariants will not only help to construct new near integrable systems, but to understand dynamics and phase space of more general non-integrable mappings with smooth f (p).This can be done by discretizing f (p) by a piecewise linear function; so we have a perturbation theory where order is determined by number of pieces and smallness parameter is defined by a smallest piece of force function.In order to demonstrate the idea, we provide a comparison of a phase space for chaotic quadratic f (p) = a q + q 2 and cubic Hénon maps f (p) = a p + q 3 around specific resonances (ν = 1/4 and stop band ν = 1/2) and corresponding polygon maps.As one can see in Fig. 25 below, even low order approximations can reveal important qualitative dynamical features, including shapes of trajectories in phase space and the topological structure of the phase portrait.
An interesting observation can be made looking at the smaller alignment index (SALI) color maps.SALI is an indicator (alternative to traditional maximal Lyapunov exponent) distinguishing between ordered and chaotic motion [38].As one can see, there is a significant area around a fixed point with almost the same color (white).It means that in this area the motion is quasi-integrable and almost linear (very small variation of rotation number).One can note that this area resembles a separatrix of inner linear layer from corresponding polygon map; thus, we learn that linear dynamics in the inner layer of integrable polygon map is not just a consequence of oversimplified discretization of force function, but a reflection of actual dynamics -negligible spread of rotation numbers.In contrast to smooth perturbation theory when in zeroth order we have an ellipse K(p, q) = p 2 − a p q + q 2 , here we have a polygon with degenerate motion.While in both approaches the rotation number is independent of amplitude, discrete perturbation theory gives better boundary of linear area when system is close to some resonant condition, while smooth perturbation theory will work better far away from any major resonances or very close to the origin where ellipses are not distorted.
FIG. 25.The middle row of plots (a.1 -c.1) shows tracking for Hénon map with quadratic f (q) = a q + q 2 (cases a. and c.), and cubic force functions f (q) = a q + q 3 (case b.).The value of parameter is (a.) a = a 1/4 − 0.005, (b.) a = a 1/4 + 0.005 and (c.) a = a 1/2 + 0.05, where a 1/4 = 0 is the main octupole resonance with ν = 1/4 and a 1/2 = −2 is the linear-half integer stopband ν = 1/2.The green curve is the second symmetry line f (q)/2 and the black dashed curve is the force function f (q).All plots correspond to a range p, q ∈ [−1, 1].The top row shows corresponding low order approximation via mappings with polygon invariants.The bottom row of plots (a.2 -c.2) shows log of SALI values for the same cases, by courtesy of Ivan Morozov (BINP).

IX. SUMMARY
In this article, we presented an algorithm for automated discovery of new integrable symplectic maps of the plane with polygon invariants, and as a result, we reported over 100 new integrable families.The algorithm successfully rediscovered some famous discrete Painlevé equations as well as some of McMillan-Suris integrable mappings.Fig. 26 presents most, but not all, systems discovered by the search algorithm.They are grouped by similarity of dynamics around the fixed point while arrows are pointed toward systems with more segments in f , thus showing hierarchy.FIG.27.Calculation of action-angle variables for map β1.Plot (a.) shows three level sets of invariant: level set from linear layer (black), separatrix (red) and level set from nonlinear layer (blue).Solid and dotted green lines are second (p = f (q)/2) and first (p = q) symmetry lines.Segments of nonlinear trajectory are labeled with Roman numerals.Plot (b.) shows partition of the area under nonlinear trajectory into fundamental rectangles and triangles.Area of green triangles is proportional to initial condition squared, area of blue rectangles is linearly proportional and red square is independent of q0.
Appendix A: Calculation of action-angle variables for polygon mappings

Action
Since invariant curves are polygons, the calculation of the area enclosed under the trajectory is a trivial task.Each polygon can be split into collection of simple rectangles and (or) triangles.Thus, the action integral (as well as partial action integral) has the form J = 1 2 π p dq = C 1 q 2 0 + C 2 q 0 + C 3 , where 1,2,3 are some constants and (q 0 , 0) is an initial condition on the trajectory.The first term represents rectangles (or triangles) with both sides being proportional to initial conditions, the next term with only one side of rectangle being proportional to q 0 and the other side determined by map parameters, and the final term represents contribution from fundamental rectangles independent of initial condition.
As an example, below we will calculate action and partial action for the outer layer of nonlinear trajectories of the map β1, see Fig. (A2)
FIG.3.Integer mappings with linear force function and polygon invariants α, β and γ (we use the same naming convention as Ref.[24,25]).Note that the invariant triangle was modified in order to satisfy both symmetries.Each figure shows one invariant polygon, with each colored region i being mapped to region i + 1.Here n is the period and ν is the rotation number of the map, which represents the average increase in the angle per map iteration.The dashed and solid dark green lines show the first and second symmetry lines.

7 FIG. 5 .
FIG.5.Fundamental polygons of the first kind for different values of rotation number ν.Each plot shows the invariant ellipse K(p, q) = p 2 − k p q + q 2 = const (black) and two invariant fundamental polygons inscribed into it (blue and red).The dashed and solid green lines correspond to the first and second symmetry lines.

FIG. 6 .
FIG.6.Integer mappings with 2-piece force function and polygon invariants D -H (named after[24]).Each figure shows one invariant polygon, with each colored region i being mapped to the region i + 1.Here n and ν are the period and rotation number of the map.The solid and dashed dark green lines show the second and first symmetry lines.

12 FIG. 9 .
FIG. 9. Shapes of fundamental polygons of the second kind corresponding to different families of maps represented in Fig.7.Each plot shows invariant ellipses (black) and invariant polygons inscribed into it (red).From left to right: line k2 = 1 for k1 = −1 (Brown-Knuth map, H), line k1k2 = 2 for k1 = −1, line k1k2 = 3 for k2 = −3 (map E).The dashed and solid green lines correspond to the first and second symmetry lines.

FIG. 10
FIG.10.Integer linear mappings with 2-piece force function and polygon invariants (stable and unstable maps).The black polygons (or polygonal chains for unstable mappings) show constant level sets of the invariant, red level sets are separatrix for unstable mappings with stable periodic dynamics, and the green line is the second symmetry line.

FIG. 11 .
FIG.11.Zoo maps: Gingerbreadman (a.), Rabbit (b.), Octopus (c.) and Crab (d.) mappings.The light and dark gray areas represent alternating zones of instability which are sectioned by two families of concentric invariant polygons resembling animals.Within each zone of instability there are three different scenarios.First, orbits that result from initial conditions (q0, p0) within any gray zone such that both q0, p0 ∈ Q.Those are stable orbits with an exact periodic motion and rational rotation number resulting in an invariant set of points.The second scenario includes all other initial conditions within the gray zones with any or both q0, p0 being irrational, q, p0 ∈ R \ Q.Those are chaotic orbits never returning to initial conditions and densely covering the zone of instability.The third and last scenario is for any initial condition within the colored regions of phase space.Those are chains of linear islands; the initial condition from these regions of phase space hops from island to island of the same group (each group is shown with its own color related to the relative orbit period).Such initial condition returns to itself and the motion is strictly periodic with the rotation number being mode-locked within the island.To view an "animal" you must rotate each image 135 degrees clockwise.

1 FIG. 12 .
FIG.12.Integer nonlinear integrable maps with polygon invariants and 2-piece force function.The top and bottom rows of the figures show level sets of the invariant for d = ∓1.The middle row shows cases for d = 0, and, in each column the transformations have the same slopes k1,2 as for d = 0 (see Fig.10); note that the invariant level set of the outer nonlinear layer at the large amplitude is asymptotic to the level set of the corresponding map.The black (blue) level sets are linear (nonlinear) trajectories.The red level sets represent the separatrix isolating the inner linear from the outer nonlinear layers.The green line shows the second symmetry line, p = f (q/2).
FIG.13.Flowchart representing the machine-assisted discovery of integrable maps.The algorithm verifies if the map defined by the input (0.) is integrable with the invariants being concentric polygons.If the map has a fixed point (1.), the algorithm performs tracking for different initial conditions (2.), and, if all orbits are bounded (3.) it proceeds to the construction of orbit polygons via ordering (4.).Finally, after counting the number of vertices V (5.), if for each polygon V < V cutoff , the map is sent for analytic verification.The search process terminates when any conditional operation is false (NO).
FIG.14.Integrable polygon mappings corresponding to a piecewise linear force function with three segments.Each table corresponds to its own value of k2 and shows the value of the shift parameter d.Black cells represent degenerate cases when 3-piece force function becomes 2-piece (k1 = k2 or k2 = k3), white cells represent cases with isolated integer d showing its values, orange cells represent mappings integrable for any integer d ∈ Z ("D"), cyan cells represent mappings integrable for continuous d ∈ R ("C") and red cells represent unstable nonlinear polygon mappings where the numbers in the cell shows the values of d corresponding to integrable maps.

FIG. 15 .
FIG.15.Integrable polygon mappings with 3-piece force function and isolated d.Black and blue lines correspond to linear and nonlinear invariant phase space trajectories, red lines are separatrices, and the thick green line is a second symmetry line, p = f (q)/2.Mappings b2 and β5 are also integrable for any d ∈ R and mappings a2 and β4 are also integrable for any d ∈ Z. Mappings α4 * and a2.1 * are unstable and were missed by the search algorithm.

3 FIG. 16 .
FIG. 16.Integer integrable polygon mappings with 3-piece force function and discrete d ∈ Z.Each row shows first 4 members of a family.B. Discrete d, (D)

2 FIG. 18 .
FIG.18.Bifurcation diagram for integer integrable polygon mappings with 3-piece force function and continuous d ∈ R. Each row shows changes of phase space with shift parameter d.Phase space plots at exact values of d correspond to bifurcations in a system and plots in between them show prototype of map for some value of d in the corresponding interval of parameter.

FIG. 19 .
FIG.19.Integer integrable polygon mappings with 4-piece force function and isolated values of shift parameter, d, and ratio of finite pieces, r.

1 FIG. 21 .
FIG. 21.Integer integrable polygon mappings with 4-piece force function and isolated values of shift parameter, d, and continuous ratio of finite pieces, r ∈ R + .Mappings with ′ index experience bifurcations of phase space with the change of parameter r.

FIG. 22 .
FIG. 22. Constant level sets of invariant on (xn, xn+1) plane for discrete Painlevé map dP S III in Suris form with A = 0 and x measured in units of B. For ϵ = ∞ map degenerate to linear with f (x) = 0 and level sets degenerate to concentric circles.When ϵ = 0 we observe ultradiscrete limit udPIII corresponding to polygon map α3.1 ′ with ratio of finite pieces restricted to unity, r = 1.Black and blue level sets show trajectories in linear and nonlinear areas of phase space, red level sets are separatrices of motion between different layers.Green curve is the second symmetry line f (x)/2 and black dashed curve is the force function f (x).

FIG. 23 .
FIG.23.The left plot (a.) illustrates invariant level sets for the Brown-Knuth map, force function f (q) = |q|.The middle plot (b.) displays invariant level sets for the Cohen map, f (q) = q 2 + 1.The right plot (c.) again provides invariant level sets for the Cohen map, but on a different scale showing one of the island structures.The level sets for Cohen map are obtained by tracking.The green curve is the second symmetry line p = f (q)/2.

2 FIG. 26 .
FIG.26.Most (but not all) mappings discovered by an algorithm: 1-through 4-piece force functions.Maps are organized by dynamics around fixed point; systems with the same rotation number around fixed points are presented in the same solid box.Dashed boxes indicates that map belongs to discrete or continuous family.

27 .
Introducing a new variable x = |q 0 − 1/2| which represents a horizontal distance from separatrix to nonlinear trajectory, we have )+9y , y ∈ [0; ∞] TABLE X. Rotation number, νi, for the i-th layer of the integer polygon mappings with 4-piece force and integrable for isolated integer value of shift parameter d and discrete ratio of finite pieces, r ∈ N + .ν isl and #C are the rotation number in mode-locked area and number of chains in group of linear islands, respectively.Part 1.

TABLE VI .
Integrable polygon mappings with 4-piece force function and isolated integer d and discrete of continuous r.Values of parameters d (in units of l2) and r are presented for each quadruplet of integer slopes.
A. Isolated d and r

TABLE IX .
Rotation number, νi, action, Ji, and partial action, J ′ i , for the i-th layer of the integer polygon mappings with 4-piece force and integrable for isolated integer values of shift parameter, d, and ratio of finite pieces, r.Part 2.

TABLE XI .
Rotation number, νi, for the i-th layer of the integer polygon mappings with 4-piece force and integrable for isolated integer value of shift parameter d and discrete ratio of finite pieces, r ∈ 2N + and n = r/2.ν isl and #C are the rotation number in mode-locked area and number of chains in group of linear islands, respectively.Part 2.