Integrable symplectic maps with a polygon tessellation

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 Poincare cross-section of a kicked rotator. 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 and associated systems on the plane, 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 Poincare 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.


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 [1][2][3].
In the study of Hamiltonian mechanics, symplectic mappings, such as Poincaré or stroboscopic crosssections 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 the 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 [4][5][6], Lagrange descriptions of fluid dynamics and billiards, to name a few (e.g., see [7]).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 [5,8,9]).This form arises from the stroboscopic Poincaré crosssection (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 [10] and Chirikov [11,12], and integrable McMillan mappings [5].Of particular significance is Suris' theorem [13], which illuminates potential integrable scenarios within this category.Suris demonstrated that if an invariant of motion K[q, p] 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 systems that go beyond the constraints outlined in Suris' theorem.Notable examples of mappings exhibiting linear dynamics, yet possessing a nontrivial phase space structure, include the Brown-Knuth map [14][15][16], McMillan beheaded and two-headed ellipses [5], and CNR maps [17,18].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 [9], 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 [19] (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 (ν) and utilizing analytical insights into the piecewise monotonic dependence of ν 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 [20]), they are currently only capable of re-identifying systems that are already known.
What adds intrigue to the discoveries is that, beyond contributing to the arenas of integrable dynamics and algorithm theory, they establish 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 K[q, p] = K[q + m L, p + n L] for m, n ∈ Z.This symmetry is attributed to the periodic nature of f (q).Consequently, the constant level sets of K[q, p] 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 [21] (see 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 [22][23][24][25], 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 paper is structured as follows.In Section II, we introduce the MH form of the map, provide definitions of force functions on both a torus and a plane, and discuss various symmetries of the invariant.Section III offers clear examples illustrating how integrable tessellations operate through the exploration of dynamic regimes, including global mode-locking on the lattice and previously unknown, to the best of our knowledge, regime of integrable diffusion.In Section IV, we introduce a new method for identifying chaos with easy-to-follow examples, while the detailed description of the automated search algorithm is left for Appendix A. Section V reveals the main discoveries, and a detailed analysis, including spectra and analytical expressions for the rotation number, is presented at the end of the article in Appendix B. Section VI discusses some possible applications of the map, including a "smoothening" procedure that allows the removal of all non-physical features, such as a lack of differentiability or discontinuities, resulting in quasiintegrable 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.
• A mapping M : X → X on a 2m-dimensional manifold X is symplectic if it satisfies the condition with Here J M is the Jacobian matrix of M, Ω is a 2m × 2m nonsingular skew-symmetric matrix, I m is a m × m identity matrix, a pair of m-dimensional variables, representing configuration q and momenta p, collectively form a phase space coordinate ζ = (q, p), and, ( ′ ) denotes the application of the map 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 J M = 1, indicating that the map preserves both area and orientation.
• The orbit of a map M is defined by a sequence of points, ζ k , obtained through the repeated application of the map to an initial condition, commonly referred to as the seed, ζ 0 : We say that the set of points {ζ 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 , it is a 2-cycle, and so on.
• A map of the plane is said to be integrable if there exists a non-constant continuous function K[q, p], called integral or invariant of motion, that remains conserved upon the application of the map for all initial conditions, i.e., This implies that the orbits originating from ζ 0 lie on a constant level of K = K(ζ 0 ), whether it manifests as a singular point, a set of points, or as a single (or collection) of closed/open curve(s).

A. McMillan-Hénon form of the map
This article focuses 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 [5,8,9], while here we briefly list some of its relevant properties.The transformation and its inverse are explicitly expressed as: 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 [11,12], the areapreserving Hénon quadratic map [10], or a collection of integrable mappings discovered by E. McMillan [5] and Y. Suris [13].
The map in MH form can be written as a composition M = F • Rot(π/2), where Rot(θ) is a rotation about the origin Rot(θ) : followed by a thin lens transformation 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 and comprising a linear optics insert in tandem with a thin nonlinear lens (e.g., sextupole, octupole or RF-station), can be redefined in the MH form.Alternatively, M can be decomposed as G • Ref(π/4), where Ref(θ) is a reflection about a line passing through the origin at an angle θ with the q-axis Ref(θ) : and G is a nonlinear vertical reflection Both, Ref(θ) and G, are anti-area-preserving involutions, which means that a double application of the map is equivalent to identity transformation Ref 2 (θ) = G 2 = I 2 and the determinant of their Jacobians is equal to −1.
Each reflection possesses a line of fixed points: designated as the first and second symmetry lines [5,21].

B. Force functions with patterns of regularity
In [9], we introduced numerous nonlinear integrable mappings with polygon invariants and corresponding piecewise linear force functions, denoted here as f 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 and examine dynamics with piecewise linear forces on a torus ζ ∈ T 2 , i.e., when the equations of motion (3) are considered modulo L, where L is referred as the lattice constant.
Furthermore, we return our study to the planar dynamics ζ ∈ R 2 , through two possible "unwrappings".The first one, results in force that consists of an infinite number of segments, strictly periodic f per (q) = f per (q + L), and might include discontinuities, that can be seen as the inclusion of additional segments with their lengths approaching zero and slopes tending towards infinity.Otherwise, we can impose the arithmetical quasiperiodicity, which implies that ∀ q : f a.q.(q + L) = f a.q.(q) + F where As evident, a function with arithmetic quasiperiodicity can be expressed as the sum of a purely periodic and linear functions where k 0 = F/L provides the average slope of f a.q.(q).In this article, we will utilize notations that define the force function as [{k 1 , l 1 }, . . ., {k n , l n }; d], where each of the n tuples {k j , l j } represents the slope and length of the j-th segment for q ∈ [0, L], followed by the value of the vertical shift parameter, d = f p.l. (0).The top row of Fig. 1 provides illustrative examples of piecewise linear forces on a plane: regular, periodic with discontinuities, and a function with arithmetic quasiperiodicity.Their representation on a torus (mod L) is shown at the bottom, offering visual cues for the introduced notations.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.
Top row illustrates a regular piecewise linear force function composed of three segments (left), along with its periodic (middle) and arithmetically quasiperiodic (right) unwrapping from the torus (mod L), depicted with a cyan square.The bottom plot displays the corresponding force on a torus, providing visual references for the introduced notations.
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: • Each map has a corresponding "twin" map, which can be obtained by using the transformation f (q) → −f (−q).Geometrically, this is equivalent to rotating the phase space by an angle of π, i.e. (q, p) → −(q, p).
• 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 * .
• For any given set of parameters, there exists a set of equidistant values of the parameter d, expressed as where m ∈ Z is an integer.This generates an equivalent dynamical system with the fixed point q * → q * +m L, owing to the translational symmetry of the invariant.
• Lastly, due to the linear nature of each segment, the rescaling of dynamical variables (q, p) → ϵ (q, p) leaves the form of the map invariant.This rescaling allows us to choose one of the segments (or, e.g., L or d) 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.

C. Symmetries of the invariant
Primarily, for the mapping in the MH form (3), any invariant level set (even amid chaos and the absence of a global invariant for all initial conditions) remains unchanged not only under M, but also when subjected to either one of the reflections, Ref(π/4) and G, respectively If the set of Eqs. ( 7) is satisfied for all points in the plane/torus, the map M is integrable.This offers a straightforward analytical test to verify integrability by confirming 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.
At the same time, when unwrapping the integrable system from the torus to the Euclidean plane, we introduce the lattice translation symmetry of the invariant which arises from the patterns of regularity in both f per (q) and f a.q.(q).Consequently, the constant level sets of K[q, p], corresponding to isolating integrals of motion/separatrices on a torus, now perfectly tessellate the plane with polygons.Using terminology from famous Archimedean tilings, we designate the tiling as regular if the phase space is tessellated with a single type of polygon separatrices, and as semi-regular when the net of separatrices consists of multiple polygon types while maintaining a coherent pattern across the entire plane.Simultaneously, all other invariant level sets contribute to a continuous fibration by polygons within all tiles necessary for integrability.In both instances, the governing tessellation/fibration rules are provided by Eq. ( 7) and are satisfied Finally, some of the newly discovered systems possess an additional discrete translational symmetry of the invariant: it is possible that, for a given tessellation by polygons, Eqs. ( 7) can be satisfied for some values of d that differ from d 0 + m (2 L − F ), where m is an integer.This symmetry will be called the configuration symmetry because it is associated with the relative position of the fixed point, which for the mapping in MH form is provided by the intersection of symmetry lines, l 1 and l 2 : For instance, in a regular tiling, the fixed point can be located either inside one of the tiles, which we refer to as the central cell (or cc for brevity), or it can be placed into the node/vertex of the separatrix net.In the latter case, we say that the map has a central node configuration (or cn for short).
It is worth noting that owing to the lattice translation symmetry of the invariant, for every integrable map with f per (q) or f a.q.(q), we can define a corresponding integrable system on a torus.However, the converse guarantees integrability only for purely periodic functions.In Throughout this paper, we employ the following color scheme to represent invariant level sets in phase space: red corresponds to isolating invariants/separatrices, while blue and black indicate phase space orbits within layers with nonlinear (amplitudedependent) and linear (mode-locked) dynamics, respectively.The solid and dashed green lines in all diagrams correspond to the second, p = f (q)/2, and first symmetry lines, p = q.other words, the "unwrapping" from the torus while attempting to enforce arithmetical quasiperiodicity can result in the second symmetry line f a.q.(q)/2 that is incompatible with the lattice.
To effectively illustrate the concepts discussed, we use a simple integrable map on a torus with a 2-piece force function defined as: Figure 2 displays the phase space portrait on the torus (left) and the corresponding "unwrapping" of the invariant to the Euclidean plane (right).Organizing the phase space portrait (mod L) on a regular square grid in R 2 reveals a regular tessellation.Each tile contains an invariant line segment with 2-cycles (depicted in black), while the remaining phase space is filled with concentric polygons (depicted in blue), bounded by the rectangular separatrix (depicted in red).This tessellation not only aligns with periodic unwrapping, as intended in its construction, but also accommodates unwrapping with arithmetic quasiperiodicity.Additionally, the system allows for different fixed point locations, either inside one of the tiles (cc) or at the separatrix nodes (cn), owing to the presence of configuration symmetry.

III. EXPLORING DYNAMICAL REGIMES
In order to describe dynamics, we consider two additional systems that are integrable on a torus, each associated with corresponding forces derived from piecewise linear functions consisting of two segments: M tor a1 : f tor (q) = [{0, L/2}, {2, L/2}; L].
Figure 3 presents phase space portraits on a torus (top row) along with their possible unwrappings to the plane.
For each integrable scenario, we provide sample orbits represented by series of connected points.In the case of mappings involving chaotic behavior, scattered dots obtained through numerical iterations are used instead.
FIG. 3. Phase space portraits for the mappings on a torus, M tor d1 and M tor a1 (top row), along with their corresponding unwrapping with arithmetical quasiperiodicity (middle row) and periodic unwrapping with discontinuities (bottom row).The left plot in the middle, for the case (d.1), is non-integrable; chaotic trajectories are depicted with points in varying shades of gray.All other systems maintain integrability, showcasing sample orbits represented by series of connected dots.

A. Dynamics (mod L)
We begin our exploration by addressing dynamics modulo L. While, in general, the main aspects closely resemble maps with regular piecewise linear forces f p.l. [9], subtle nuances exist.Primarily, for regular mappings on the plane, the orbits can be either stable (polygons) or unstable (simple polygonal chains).However, due to the finite and compact nature of the torus surface in contrast to the infinite and unbounded plane, we refrain from explicitly considering stability.Instead, we discern between oscillations, where the trajectory rounds the fixed point, tracing out closed-loop paths (as can be seen for the map (d.1)), and, drift or flow, where the orbit "weaves" and "loops" around the torus upon reaching its border.One can draw an analogy to the two regimes of motion observed in a classical pendulum.Nonlinear flow implies trajectories that are more intricate than simple circular or closed-loop paths, as they can traverse different regions of the torus, as seen in example (a.1) in Fig. 3.
Moreover, concerning planar dynamics, we categorize stable orbits as either purely periodic (n-cycles) or quasiperiodic, contingent on whether the rotation number is rational or irrational.This concept extends to both types of motion on a torus, where the term winding number is sometimes employed to distinguish drift from oscillations with rotation number.Mapping (d.1), being degenerate, yields all orbits as 7-or 3-cycles, while the dynamics of map (a.1) is nonlinear, leading to a dense coverage of most orbits across their invariant level sets.

B. Global mode-locking
Next, we turn our attention to planar dynamics with continuous forces exhibiting arithmetical quasiperiodicity.As previously noted, unwrapping the system in this way while maintaining integrability may be impossible, as illustrated by the left plot in the middle of Fig. 3, where the presence of chaotic orbits is evident.On the contrary, mapping (a.1) facilitates integrability for f a.q.(q) and generates a regular tessellation with a central cell, similar to the previously examined map (b.1).
While all tiles of the same type are geometrically identical, the dynamics on similar level sets of the invariant vary from one cell of periodicity to another.Foremost, there exists a unique tile (or a unique node in the case of the cn configuration) containing a fixed point.Within the cc, the dynamics revolve around the fixed point, akin to the scenario observed in regular polygon maps.
Furthermore, all initial conditions on the separatrix isolating cc correspond to P -fold periodic orbits (Pcycles), causing the non-central tiles to form chains of P islands.As a result, initial conditions from non-central tiles exhibit global mode locking on the lattice with a rational rotation number defined by: Here, we refer to P ∈ N and ν P ∈ Q as the global period and global rotation number of the map.These values are in accordance with Crystallographic theorem [9,17,18] and naturally emerge when dealing with integer coefficients.It remains an open question whether it is possible to organize segments of the force function to produce a tessellation with P values different from the aforementioned ones, requiring k 0 ∈ Q.A more detailed illustration is presented in the right plot in the middle of Fig. 3, providing an example of two orbits generated from initial conditions separated by one lattice period.The blue and cyan dots/polygons correspond to initial conditions within the central cell (q 0 , p 0 ) ∈ cc and (q 0 + L, p 0 ), respectively.In the latter case, the trajectory sequentially visits six islands under iterations and traces out a collection of invariant tori, see additional plot in Fig. 5.

C. Integrable diffusion
Here, we proceed to the last scenario of periodic functions resulting from simple unwrapping from the torus, Eq. ( 4).While the continuity of f a.q.(q) from the previous subsection guarantees the continuity of the invariant of motion, allowing global mode-locking, the inclusion of discontinuities can lead to significant differences.
The left plot at the bottom row of Fig. 3 provides a corresponding unwrapping for the map (d.1), showing another sample orbit, this time departing from a noncentral cell located one lattice period to the right.As one would expect, the sample orbit is mode-locked to 1/4, since the average slope k 0 is equal to zero, and it hops between four islands, forming 7 × 4-cycles (or 3 × 4cycles if departing from a non-central triangle).
While unwrapping from the torus might appear somewhat trivial, involving the creation of new invariant level sets through simple translations, unexpected behavior arises in the case of mapping (a.1), illustrated at the bottom right corner of Fig. 3. Surprisingly, after several iterations, the orbit originating inside the central cell escapes and begins to hop between tiles in a rather intricate manner.However, instead of exhibiting a "random" trajectory, the orbit remains on a specific level set, as further demonstrated in Fig. 5.
Comparing this case to the previously considered unwrapping with arithmetic quasiperiodicity, as shown in Fig. 5, reveals that instead of a finite number of closed rectangles, we now have infinitely many sets of rectangles split into two pieces.The two plots at the bottom provide 10 4 and 10 5 iterations of the same sample orbit as in Fig. 3 but on a different scale, clearly depicting diffusion.Notice how the progression in time unveils more and more segments of the invariant, aligning with the observed pattern.We will refer to this behavior as "integrable diffusion." Now, we can raise the question: "Why do cases (a.1) and (d.1) exhibit such different behavior when considered on the plane?"The resolution to this puzzle lies in how the network of separatrices (shown in red) is organized within the square representing the torus: in case (d.1), separatrices form closed loops and are confined within T 2 , thus providing isolating invariants when the system is unwrapped.On the other hand, for case (a.1), the boundary of the torus cuts the separatrix, creating gaps in isolating invariants and allowing for diffusion.In fact, one can observe that both behaviors fit into the larger picture where we had stable and unstable orbits for regular piecewise mappings of the plane, then the corresponding pair of rotations around a fixed point (oscillations) and flow (drift or winding) on a torus, and finally, global mode locking and integrable diffusion when force exhibits patterns of regularity.

FIG. 5. Top row displays the replication of invariant level
sets traced out by sample trajectories for the unwrappings of (a.1), Fig. 3.The bottom row illustrates integrable diffusion for the top-right plot, depicting 10 4 and 10 5 iterations obtained by tracking, with a brown square having a side length of 2 L for reference.

IV. CHAOS DETECTION
After demonstrating the possibility of constructing nonlinear integrable mappings allowing tessellation, a new question emerged: "How can we find more systems and automate the search process"?In [9], where f p.l. (q) represented a standard piecewise linear force function, our approach relied on the assumption that the phase space should be fibrated with concentric polygons.By organizing points of an individual orbit based on polar angles and counting the number of vertices, we could identify non-integrable systems where the polygons produced after rearranging chaotic trajectories had an increasing number of vertices that did not converge with the number of iterations.
Unfortunately, as mentioned in [9], this algorithm fails when the phase space involves islands with nonlinear dynamics, i.e., in the presence of global mode-locking: arranging the points of the corresponding orbit by polar angle becomes problematic, as K(r, θ) = const transforms into a multiple-valued function of θ defined on disjoint intervals, as can be seen in top left plot in Fig. 5.
To address this issue, we decided to rely on intrinsic dynamical variables rather than analyzing the structural aspects of the phase space.We first recall the canonical form of an integrable system [26][27][28]: where J k and ψ k are the action and angle variables at the k-th iteration, and ν denotes the rotation number which captures the average angular increment of the angle variable during a single iteration of the map.ν remains invariant regardless of the mapping's representation, and for instance, in (q, p)-coordinates, it can be analytically expressed using Danilov's theorem [29][30][31] or as a limit: where θ k = arctan[(p k − p * )/(q k − q * )] 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 the phase space is effectively segmented into layers defined by the pieces of the force function and "demarcated" by separatrices.The utilization of Danilov's theorem uncovers that within each layer the Poincaré rotation number must be of the form [9] where 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 f a.q.(q) derived from the subsequent parameter sets: To analyze these dynamical systems, we will utilize the rotation number as a function of the horizontal distance from the initial condition to the fixed point, denoted by s = q 0 − q * , along both symmetry lines: Both spectra are evaluated within two lattice periods, as indicated on the phase portraits by a dashed cyan square.

A. Integrable map
First, let's explore a more intricate example of a nonlinear integrable system featuring a semi-regular tiling pattern, illustrated in Fig. 6.Around the fixed point, there are three layers: a line segment with 2-cycles (black line), succeeded by two layers consisting of concentric polygons with nonlinear dynamics (depicted by blue contours separated by red separatrices).The cc is replicated throughout the phase plane, and the gaps between the bounding octagons are filled with invariant triangles.The average slope k 0 of the force f a.q (q) is equal to −1, leading to the organization of non-central tiles into chains of three islands, with a corresponding global rotation number of ν P = 1/3.As anticipated, the spectra display piecewise monotonic behavior inside the cc, indicated by corresponding colors, with global mode-locking outside matching ν P , represented by dashed black lines.FIG. 6. Phase space portrait and spectra along symmetry lines for the case with integrable dynamics.FIG. 7. Phase space portrait and spectra along symmetry lines for the case of "chaos in cell."

B. Chaos in cell
Moving on to the next example shown in Fig. 7, 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 chosen an extreme case where the network of separatrices is the only set of invariant tori that remains intact.These isolating invariants wield a significant influence on the system's dynamics by effectively segregating cells of periodicity.The average slope is k 0 = 0, resulting in a global mode locking of 1/4.Consequently, all points located on the separatrix possess a period equal to P = 4.
Within the central cell of periodicity, there are two additional invariant structures of unit measure that 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. 4).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 nonisolated 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.Red dots represent the locations of corresponding separatrices.

C. Chaotic diffusion
Finally, we consider a scenario where perturbation primarily impacts the larger amplitudes of the system (within the lattice period), causing a complete disruption of the invariant isolating structures, as illustrated in Fig. 8.In this selected example, the phase space consists of two integrable fragments following a semi-regular pattern and separated by a "chaotic sea."The central cell (and corresponding tiles) exhibits two distinct layers surrounding the fixed point: concentric hexagons (depicted in black) and octagons (highlighted in blue).The second integrable fragment is composed of three layers made of triangles, pentagons, and octagons.All layers shown with contours in black experience degeneracy, with every initial condition being periodic, while contours in blue represent layers with amplitude-dependent dynamics, ν(s) ̸ = const.
Beyond the larger surviving octagons that encompass both types of tiles (marked in red), a region of chaos unfolds, where, due to the absence of a global isolating invariant, chaotic trajectories can diffuse between lattice periods.For this reason, we will refer to this scenario as chaotic diffusion.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) and does not "trace out" invariant level sets, as observed in the case of integrable diffusion.Figure 9 illustrates the trajectories of the same two samples as shown in Fig. 8, but on a different scale in both space and time.As evident from the sequence of plots, as the number of iterations increases, both orbits progressively diffuse out of the central cell and expand further in the phase space.
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 at 1/4 and noisy patterns representing chaotic zones.Refer to the two plots on the right in Fig. 8 for illustration.
Flowchart and step-by-step decision list for the algorithm designed to distinguish chaotic and integrable spectra based on the assumed piecewise-monotonicity of ν(s) are presented in Appendix A.

V. NEW MAPPINGS
A systematic exploration within the same parameter space as in [9] has unveiled a notably larger number of new integrable scenarios.Almost all of over a hundred previously discovered systems with regular piecewise linear forces can be in some way wrapped on a tor.At this point, the challenge for mappings integrable on a torus shifts towards the classification of these new findings and ultimately comprehending the entirety of all possibilities, rather than merely expanding the count of discoveries.Fig. 10 illustrates several examples of phase space portraits for systems integrable on a torus.All forces are derived from the simplest piecewise linear functions comprised of two segments of an equal length, L/2.As continuous forces with arithmetic quasiperiodicity occupy a significantly smaller set in the parameter space, we were able to classify most new findings and present them in a compact manner.During our parameter scan for functions composed of two pieces and integer coefficients, no additional integrable cases were discovered, besides the aforementioned tessellations that can be produced by unwrapping mappings (a.1) and (b.1) from the torus, as shown in Fig. 11.However, as we proceeded to consider f a.q.(q) formed from three line segments, we discovered additional mappings.Some of these systems allowed for further extension to four pieces and/or generalization of the length of line segments to the field of rational numbers.Phase space portraits, along with bifurcation diagrams for the tiles, are presented in Fig. 12, while a detailed description of the dynamics, along with analytical expressions for rotation numbers and spectra, is provided in Appendix B. Each diagram is accompanied by a sketch of one lattice period of the force function (shown in brown), along with correponding values of slopes and length of segments.

VI. APPLICATIONS OF MAPPINGS WITH PERIODIC AND ARITHMETICALLY
QUASIPERIODIC FORCES.
Applications of integrable maps with polygonal invariants were extensively discussed in [9].In this section, we specifically highlight some points related to tessellations.
12. Phase space portraits and bifurcation diagrams for integrable maps with continuous arithmetically quasiperiodic fa.q(q).Each case is provided with one period of the force function (brown line) along with values of slopes and length of segments.Throughout this paper, we use red lines to depict separatrices, and blue and black to show invariant level sets within the areas with amplitude-dependent or linear (mode-locked) motion, respectively.Dashed and solid green lines correspond to the first and second symmetry lines.For instance, by setting δ = Λ 3 /Λ 1 and Λ 2 = Φ = 0, we obtain and f (q) = 2 arctan δ sin q 1 + δ cos q .
As seen, when δ = 1, the force function reduces to f (q) = 2 arctan [tan(q/2)].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 (mod L) with L = 2 π.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 [9], including the smooth form (11) and the polygon tessellation.The middle plot in Fig. 13 shows the polygon invariant for δ = 1, surrounded by complimentary cases with δ ≶ 1.

B. 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 that can be used with all discovered systems.By "smoothening" f (q) [9,32], we can obtain quasiintegrable systems where the phase space volume occupied by chaotic trajectories decreases with the smoothness parameter, ϵ.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. 13 provides two examples of "perfect smoothening" (since resulting systems are perfectly integrable) with respect to the parameter ϵ = 1 − δ.
By using other rules, as long as the new function is close enough to its piecewise linear prototype, we can obtain near-integrable dynamics.Using 4-piece mappings M 6 a2 and M 4 F1 ′ for illustration, if all segments have equal unity length, we can approximate forces as where we introduced smooth differentiable form (when ϵ ̸ = 0) of the trapezoid 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 [11,12], the typical kick from the RF field on a particle (∝ 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 ̸ = 0 represents regimes with net acceleration/deceleration.

VII. 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 with patterns of regularity, 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 piecewise monotonic behavior of rotation number concerning 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.Another interesting future direction would be to study a quantum version of the discovered integrable maps.Proceeding to force functions formed from three line segments, we encountered an intriguing map This system 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 16, where ν 1,2 are given by: and 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.Both systems have k 0 = 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. 4 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. 17.Note that 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 α3 and F1 ′ in [9].The final group of systems to be presented in this section involves the extension of the mapping M 6 a1 to a force function composed of four pieces and M 3 b1 to forces made of three or four segments.The first extension is M The l 4 ≥ l 2 condition ensures the exclusion of twin maps.
Figures 18 and 19 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.The gaps in between are filled with triangles representing the degenerate center and/or chains of linear islands (black contours).The map M 6 a2 has a singular configuration, featuring a central cell surrounded by a chain of six identical islands.On the other hand, M 3 b2 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 M 3 E1 ′ ), one of the triangular chains vanishes, leading to the corresponding degenerate center evolving into a central node.Moreover, the bottom plots in Figure 19 showcases sample spectra for various parameter values of l.The expressions for ν 1,2 are as follows:

FIG. 2 .
FIG.2.Phase space portrait for the map M tor b1 along with the corresponding unwrapping of the invariant to R 2 .Various second symmetry lines illustrate possible configurations of the map: periodic force with discontinuities (per), or, arithmetically quasiperiodic force allowing for the central cell (cc) and central nodes (cn).Throughout this paper, we employ the following color scheme to represent invariant level sets in phase space: red corresponds to isolating invariants/separatrices, while blue and black indicate phase space orbits within layers with nonlinear (amplitudedependent) and linear (mode-locked) dynamics, respectively.The solid and dashed green lines in all diagrams correspond to the second, p = f (q)/2, and first symmetry lines, p = q.

4 FIG. 4 .
FIG. 4. Schematic illustration of the global mode-locking around the central cell [cc] or central node [cn] for different values of ν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
Figure 4 schematically depicts the arrangement in which tiles are grouped into islands, with P values of 6, 3, and 4. The occurrence of these specific periods is not arbitrary and is determined by the integer values of k 0 = 1, −1 and 0.

FIG. 8 .
FIG.8.Phase space portrait and spectra along symmetry lines for the case of "chaotic diffusion."

FIG. 13 .
FIG. 13.Phase space portraits for McMillan-Suris mappings with trigonometric polynomial 11.Different plots correspond to different values of δ.When δ ̸ = 1, the dashed brown line shows linear flow on a torus (case δ = 1) for the reference.

Figure 14
Figure 14 provides an illustration of phase space portraits for these systems obtained by tracking for ϵ = 10 −1 .The last example at the bottom depicts tracking for the smooth square wave force Π[q] = 4 π arctan 1 ϵ sin π q 2 which, for ϵ = 0, represents another limiting case of the McMillan-Suris map with k 0 = 0.

FIG. 16 .
FIG.16.Phase space portrait (top) and the corresponding spectra along the first (ν1) and second (ν2) symmetry lines for the map M 6 G .

FIG. 17 .
FIG. 17. Bifurcation diagram (top) and phase space portraits (middle) for mappings M 4 α3 and M 4 F1 ′ .The bottom row shows spectra for the cc configuration and different values of l3 (shown in different colors).

6 FIG. 19 .
FIG. 19.Phase space portrait for the map M 6 a2 .Bifurcation diagram and spectra for different values of parameter l (shown with different colors) are provided to the left and right respectively.