Stable and unstable regimes in higher-dimensional convex billiards with cylindrical shape

We introduce a class of convex, higher-dimensional billiard models which generalise stadium billiards. These models correspond to the free motion of a point-particle in a region bounded by cylinders cut by planes. They are motivated by models of particles interacting via a string-type mechanism, and confined by hard walls. The combination of these elements may give rise to a defocusing mechanism, similar to that in two dimensions, which allows large chaotic regions in phase space. The remaining part of phase space is associated with marginally stable behaviour. In fact periodic orbits in these systems generically come in continuous parametric families, sociated with a pair of parabolic eigen-directions: the periodic orbits are unstable in the presence of a defocusing mechanism, but marginally stable otherwise. By performing the stability analysis of families of periodic orbits at a nonlinear level, we establish the conditions under which families are nonlinearly stable or unstable. As a result, we identify regions in the parameter space of the models which admit non-linearly stable oscillations in the form of whispering gallery modes. Where no families of periodic orbits are stable, the billiards are completely chaotic, i.e.\ the Lyapunov exponents of the billiard map are non-zero.


Introduction
Billiard models, in which a point particle moves uniformly until it undergoes abrupt elastic collisions with a fixed boundary, are the playground of statistical physicists and mathematicians alike, whose work focuses on the interplay between their dynamical and statistical properties [1].
There are two main categories of chaotic billiards. The better-known ones are (semi-)dispersing billiards, of which the hard-sphere gas is the prototypical example. The motion of hard spheres in a bounded region, undergoing elastic collisions with each other, is equivalent to that of a billiard model -a point particle which moves uniformly in the exterior of a collection of spherical cylinders in a high-dimensional phase space, with specular collisions at the boundary [2]. The Sinai billiard [3] and the periodic Lorentz gas [4] are two-dimensional examples of dispersing billiards, in which a particle moves outside a disk on the 2-torus or a periodic configuration of them on the plane. The latter is a useful model of tagged particle diffusion in a binary mixture [5]; it enjoys fast decay of correlations and thus converges to a Brownian motion [6]. The mechanism giving rise to chaos in these billiards, as well as in hard sphere gases, is that of dispersion, where nearby trajectories separate exponentially fast as a function of the number of collisions with surfaces of positive curvature, implying that the system has a positive Lyapunov exponent [7] and is thus chaotic. Indeed, hardsphere gases are essentially the only systems of interacting particles for which rigorous results, ranging from hyperbolicity to exponential decay of correlations, have been firmly established [8,9,10].
The second category is that of defocusing billiards, the paradigm of which is the Bunimovich stadium [11,12]. Here, chaos is due to a mechanism different from dispersion, namely the defocusing mechanism. As opposed to the Sinai billiard, the boundary of the stadium curves inwards with respect to the particle motion, i.e. it is convex. Although nearby trajectories initially focus after a collision with this boundary, if the distance to the next collision is longer than the distance to the focal point then they eventually defocus even more. This mechanism thus leads to an overall expansion in phase space, again measured by a positive Lyapunov exponent [13].
Since its discovery, the mechanism of defocusing has attracted much attention in the physics community, particularly in connection to quantum chaos [14], acoustic experiments in closed chaotic cavities [15], optical microcavity laser experiments [16,17], or quantum conductance experiments [18] to name but a few.
In spite of its potential appeal to a broad range of physical applications, it has, however, remained a difficult problem to establish the conditions under which the defocusing mechanism can be extended to higher dimensions [19,20]. There are still only a few models of higher-dimensional stadia based on three-and higher-dimensional cavities which have lent themselves to a systematic study and are known to be fully chaotic [21,22,23,24,25].
There are arguably two difficulties which must be dealt with in order to obtain chaotic behaviour: (i) two curved dimensions, as in a spherical cap, often generate stable oscillations; and (ii) three-dimensional billiard domains with a single curved dimension, such as a cylindrical surface, have flat components which complicate the stability analysis.
The first of these two difficulties is related to the fact that one cannot simply construct a three-dimensional surface of revolution by rotating, for example, a two-dimensional stadium to obtain a chaotic three-dimensional cavity, since angular momentum is then conserved. Furthermore, as shown by Wojtkowski [20], intersecting three-dimensional spherical caps and flat components may produce stable periodic orbits, giving rise to billiards with mixed phase space. The models studied by Bunimovich and Rehacek [21,22,23] provide an exception to this principle and, in this sense, are rather special. They are furthermore non-convex.
A natural alternative to spherical caps is to use boundary components made out of cylindrical caps and flat planes. It is clear that "extruding" or extending the two-dimensional stadium to a three-dimensional cylindrical-type stadium (a cylindrical shape whose cross section is a stadium) by a perpendicular translation does not suffice to produce a fully chaotic billiard in three dimensions, since the motion in this perpendicular direction is trivial. Thus Papenbrock [26] proposed to consider a three-dimensional stadium with half-cylindrical caps along perpendicular axes at both ends of a cuboidal shape. This construction gives rise to a billiard which is fully chaotic, as was recently shown rigorously by Bunimovich and Del Magno [24,25], and is a particular case of more general constructions discussed by Wojtkowski [27]. In this respect, it is the first known example of a three-dimensional billiard in a convex domain with this property.
The specificity of the Papenbrock three-dimensional stadium is that its cylindrical caps are placed opposite each other in such a way that periodic orbits visiting the caps are isolated. That is, even though each cylindrical component is curved along a single direction, , there is no degree of freedom with which to move the periodic orbits around, since the axes of the two cylinders are mutually transverse.
The Papenbrock stadium thus avoids the second of the two difficulties we alluded to above, which are typical of three-dimensional cavities constructed out of cylindrical components and flat planes: namely, that their periodic orbits generically belong to continuous parametric families. Take, for example, the three-dimensional cylindrical-type stadium discussed above and break its translational symmetry by cutting it with an oblique plane. A periodic orbit of this system can be shifted along the direction of the cylinder axis, remaining basically unchanged, thus giving rise to a one-parameter family. This has the consequence that its stability analysis yields, out of the four phase-space dimensions of the billiard map, two parabolic directions (with unit eigenvalue), corresponding to motion along the family.
Another example of a three-dimensional cavity whose periodic orbits belong to continuous parametric families is the three-dimensional billiard whose domain consists of the intersection of a sphere with a cuboid [28]. Though this billiard was numerically found to have four non-zero Lyapunov exponents in a large interval of its parameter, this property is yet unproven as it does not lend itself to the kind of analysis performed in refs [21,22,23,24,25].
It is the goal of this paper to consider billiards whose phase space combines neutral directions with curved regions, associated with stable and unstable behaviors, and explore the conditions under which they can display hyperbolic regimes, i.e. such that all pairs of opposite Lyapunov exponents are non-zero. To this end, we introduce a class of higher-dimensional convex billiards, in which the particle motion occurs inside cylindrical domains bounded by oblique planes; we call these cylindrical stadium billiards. The reasons for using the term stadium are manifold. As we shall show, these billiards share some of the essential properties of the two-dimensional stadium billiard. In particular, and foremost, the mechanism that drives instability is of the defocusing type; other common features of specific interest to us, and which will be discussed below, are the existence of classes of periodic orbits similar to the bouncing-ball and whispering-gallery modes of the stadium billiard [29].
The motivation for such a general class of models arises from models of particles which interact via virtual strings, as introduced by Papenbrock to model self-bound nuclei [30,31] and which were subsequently modified by Gilbert and Lefevere to model heat conduction by gas particles trapped in a nano-porous solid matrix [32]. The interaction between pairs of neighbouring particles is a sort of hard-core interaction at a distance, and works as follows: particles move freely until they reach a critical relative distance, at which point they interchange the longitudinal components of their momenta. This can be thought of as an interaction mediated by a string, which affects the particles only once the string becomes fully extended, at which point it applies an instantaneous impulse to each particle so they flip directions. A chain of particles with such interactions were considered in [32], but placed inside a rectangular domain divided up into square cells, with each particle confined to its own cell. As their study shows, the combination of specular reflections along flat walls and circular arcs gives rise to chaotic motion in this system. Thus, given a system of N ≥ 2 particles, they find 2N − 1 pairs of non-zero Lyapunov exponents opposite each other (the single pair of vanishing exponents corresponds to the energy conservation and associated time translation symmetry).
As we shall shortly describe further, models of particles with such a string-type interaction, and confined by hard walls, are equivalent to billiards in higher-dimensional configuration spaces inside convex regions bounded by cylinders and planes. This motivates the introduction of a general class of convex billiards bounded by cylinders and planes, which includes such models. Note that this is the opposite of hard-sphere gases, which are equivalent to billiards outside cylinders.
The presence of a flat direction along the axes of the cylinders in such systems implies that periodic orbits come in continuous parametric families, with each periodic orbit displaying marginal stability along the direction of the family (each dimension associated with a flat direction actually gives rise to a pair of parabolic eigenvalues along the corresponding direction and conjugate momentum). Although the existence of such a family was noted in ref. [33], its detailed analysis is notably absent in the models previously considered. Indeed, the stability of such families of periodic orbits must be studied at a nonlinear level in order to determine whether or not part of the family can or cannot be stabilised -a daunting task.
As a first step towards a comprehensive understanding of the dynamics of such systems, we present in this paper a detailed analysis of the dynamical regimes occurring in one of the simplest such systems, which consists of a cylinder in three dimensions with a flat bottom and cut by an oblique plane at its top. This analysis will be performed in terms of the stability of periodic orbits, and comes in two parts. At a linear level, stability is understood in terms of the lengths of the segments of the orbit, and, specifically, that of paths corresponding to two collisions with the cylindrical surface, separated by a single collision with the oblique plane. As we shall show, when this length is large enough, defocusing takes place and the periodic orbit is unstable. If on the contrary, this length remains small, then the orbit may remain stable. However, because such periodic orbits have a pair of parabolic eigenvalues, the effect of nonlinear perturbations along the neutral directions must be taken into account, which is the second part of the stability analysis.
There are essentially two classes of periodic orbits of cylindrical stadium billiards whose segments' lengths can remain small enough that they are potentially stable. We identify these two classes as planar and helical periodic orbits. The former are similar to bouncing-ball orbits in the stadium billiard, and the latter to whispering-gallery modes. By analysing in detail the existence conditions and linear and nonlinear stability properties of these two classes of orbits, we show, at the linear level, the existence of a bifurcation from marginally-stable to unstable (hyperbolic) regimes, and, at a nonlinear level, the existence of a "restoring force" which, in specific regions of the parameter space, is able to stabilise the motion by restricting it to a certain part of the parameter space where the family of periodic orbits remains marginally stable, thus preventing the periodic orbit from escaping to a hyperbolic region through the bifurcation point. The parameter regions where none of the periodic orbit families can be stabilised correspond to chaotic regimes of the system. The paper is organised as follows. In section 2, we provide a detailed description of particles interacting through a string-type mechanism in terms of cylindrical billiard models. In section 3 we introduce the particular three-dimensional cylindrical stadium billiard described above, and we offer a general characterisation of its periodic orbits in section 4. The two classes of planar and helical periodic orbits are studied in further details in sections 5 and 6, respectively, where we identify parametric regions of non-linear stability of these periodic orbits. A summary and discussion of our results are presented in section 7. Two appendices, Appendix A and Appendix B, provide exact forms obtained for the linear and non-linear analysis of some of the periodic orbits studied in sections 5 and 6.

Interacting-particle models and cylindrical stadium billiards
We start by relating interacting-particle models with string-type interactions and billiard models, and then introduce the general class of cylindrical stadium billiards.

Particles with string-type interactions and flat walls
Consider N ≥ 2 particles interacting via a string-type interaction [30,31,32]. The particles are joined by strings in some configuration, and have the following dynamics. Two particles which are joined by a string move freely until they reach a mutual distance r i j = , at which point the string becomes taut (fully extended), and exerts an instantaneous impulse on each particle, which interchanges the components of their momenta along the line between them. The interaction potential between the particles due to the string, as a function of the distance r, is thus v(r) = 0 for r < and v(r) = ∞ for r > .
The dynamics of the N-particle system corresponds to the motion of a point Γ ≡ (q, p) ≡ (q 1 , . . . , q N , p 1 , . . . , p N ) in a 2dN-dimensional phase space, where d is the spatial dimension of the system, so that q i ∈ R d and p i ∈ R d .
If particles i and j are joined by a string of length , then they are constrained by the inequality which specifies a hyper-spherical region S ∈ R 2d in the dimensions spanned by the coordinates q i and q j . The other coordinates are free, so that this corresponds to an allowed cylindrical region S × R d(N−2) in configuration space. This is analogous to the case of hard-sphere interactions [2], for which the constraint has the opposite inequality, q i − q j 2 ≥ 2 . Figure 1. Two interacting particles joined by a string and trapped in a square box. Their relative positions are constrained to a disk of radius , while each particle is reflected elastically by the horizontal and vertical walls. The two continuous lines show the trajectories of each particle as they interact and collide with the walls. The motion of the center of mass (dashed line) is uniform in between collisions of either particle with the walls. The dotted arrows indicate the directions of the velocity vectors at the initial and final positions.
The particles may also be confined by hard walls such as shown in figure 1, where two particles joined by a string are trapped to a square box. Each wall imposes a further constraint on the position q i of each particle, namely (q i − c) · n ≤ 0, where c is a point on the wall and n the outgoing normal vector. Each wall thus imposes N constraints on the total configurationspace position vector q, which are of the form (q −c) ·ñ ≤ 0, whereñ ≡ (0, . . . , n, . . . , 0), with n appearing in the ith position, and similarly forc; these correspond to half-spaces bounded by hyper-planes.
The configuration space vector q is thus confined to the intersection of the cylinders coming from the string-type interactions, and the half-spaces corresponding to the bounding hard walls; this region is the intersection of convex sets, and is thus itself convex. Due to the fact that the interactions are hard-core, the motion in phase space is in fact equivalent to a semi-focusing billiard inside this convex set [24]; this should be compared to the semidispersing billiards which arise from hard-sphere gases [2].

Simplest interacting-particle model
In order to get a grasp on the types of phenomena which can occur in such models, it is instructive to look at the simplest model in this class, which consists of a two-dimensional string-bound diatomic molecule, with string length (maximum separation) , confined in a channel formed by two parallel, infinitely long, hard walls at y = ± d 2 , i.e. similar to the example shown in figure 1, but with the vertical walls removed and horizontal walls infinitely long. We denote the particle positions and velocities by q i = (x i , y i ) and p i = (u i , v i ), respectively, for i = 1, 2. The string constrains the particles by q 1 − q 2 2 ≤ 2 , and the above argument shows that the hard walls give 4 allowed half-spaces, given by extending the inequalities − d 2 ≤ y i ≤ d 2 to the configuration space vector q. Although the phase space of this system is a priori 8-dimensional, translational symmetry along the channel implies that the horizontal component of linear momentum of the centre of mass is conserved, so that we may choose a reference frame in which it vanishes, u 1 + u 2 = 0, and in which we can also fix its position, x 1 + x 2 = 0. It is convenient to rotate the coordinate system by angles π/4 and normalise them according to , setting X ≡ (x 1 − x 2 )/ , Y ≡ (y 1 − y 2 )/ , and Z ≡ (y 1 + y 2 )/ . This gives an equivalent billiard system, consisting of a point particle which bounces elastically inside the region bounded by the cylinder X 2 +Y 2 ≤ 1, and chopped by the four planes Y ± Z = ±h, with h ≡ d/ . This, in turn, is equivalent to a billiard in a domain bounded below by the plane Z = 0 and above by the two planes Z ±Y = h, shown in figure 2.
Similar constructions give rise to cylindrical stadia in higher dimensions. For example, in the case of the diatomic molecule inside a 2-dimensional square box with hard walls shown in figure 1, there is no longer conservation of momentum, so that we obtain a billiard inside a cylinder with two curved and two flat directions in 4 dimensions, bounded by eight hyperplanes.

Cylindrical stadium billiards
The examples of particles interacting via string-type interactions motivate our introduction of a large class of cylindrical stadium billiards. Informally, these consist of one or more cylinders which are cut by (hyper-)planes. More formally, they are the intersection of one or more infinite cylindrical region(s) with two or more half spaces, such that the resulting billiard is bounded and such that at least part of the boundary of the resulting region is curved (that is, comes from a cylinder). This last condition is required since otherwise the region is polyhedral, and the dynamics is known to be non-chaotic [13].
Note that these billiards do not share the property of the Bunimovich and Papenbrock stadia that the cylindrical and flat portions of the boundary meet tangentially (C 1 boundaries). The motion of two interacting particles joined by a string in a two-dimensional infinite channel is equivalent to the motion of a point particle in a three-dimensional cylindrical domain chopped by four planes, or, equivalently, a cylindrical domain with a flat bottom and topped by two interesting planes at z ± y = h, as shown above. This billiard has a single parameter, the height of the oblique planes, given by the ratio of the width of the channel to the length of the string of the interacting di-atomic model. Nonetheless, they are higher-dimensional convex billiards with large chaotic regions, and in this sense are a possible generalisation of the stadium which are of intrinsic interest. Note also that the stadium is only one of a large class of 2D billiards proved to be chaotic, many of which do not share its property of having C 1 boundaries [13].

A three-dimensional cylindrical stadium billiard
In order to investigate in more detail the types of mechanisms which are present and which give rise to chaotic and/or stable motions in cylindrical stadium billiards, we henceforth restrict attention to one among the simplest possible such systems, which is a simplified version of the 3D cylinder found in the previous section.
We consider a cylindrical stadium billiard in R 3 formed by the cylinder x 2 + y 2 ≤ 1 with unit radius and vertical axis, and bounded below by the plane z = 0, which is perpendicular to the axis, and an oblique plane that is inclined away from perpendicular. Here we focus on a one-parameter family of such billiards, where the inclined plane is placed at an angle π/4 with respect to the cylinder axis, and is parameterised by its height h, with equation y + z = 1 + h. The parameter h ≥ −2 thus defined measures the minimal height of the cylinder at y = 1 ("above" its base); h < 0 corresponds to the case where the inclined plane cuts the flat plane inside the cylinder. The convex billiard domain then consists of the intersection of the cylinder with the two half-spaces z ≥ 0 and y + z ≤ 1 + h -see figure 3.
We view this convex billiard as the elementary cell of an expanded square-shaped cylindrical billiard, or more simply squared cylindrical stadium, obtained by "unfolding" the elementary cell, that is, repeatedly reflecting it in its two planes. This cavity is thus made out of the union of four sections of cylinders, which intersect pairwise at right angles -see figure 4. The only relevant parameter is the height of each of these cylinders; their half-height is the parameter h introduced above, which we will refer to as the geometrical parameter. It is bounded below by h = −2, but is not bounded above.   Considering the unfolded billiard, we will refer to the plane generated by the axes of the four cylinders as the plane of the billiard. A transverse plane refers to any plane perpendicular to that plane. These include the planes perpendicular to the cylinders' axes, as well as the oblique planes in which they intersect.
Choosing a Cartesian coordinate system whose center is at the center of symmetry of the square, we take the plane x = 0 to be the plane of the billiard. The cylinder axes are then located at x = 0 and y = ±(1 + h), parallel to the z-axis, and at x = 0 and z = ±(1 + h), parallel to the y-axis.
The unfolded billiard undergoes significant topological changes at h = −1 and h = 0. When h = −1, each pair of parallel axes coalesces, so that the unfolded billiard is a convex cavity for −2 ≤ h ≤ −1. When h = 0, two parallel cylinders intersect along a single line on their surfaces -they intersect for h < 0 and separate for h > 1. The unfolded billiard thus has genus 0 for h < 0 and genus 1 for h > 0.

Periodic orbit families
Much in the same way that one distinguishes bouncing-ball and whispering-gallery periodic orbits in the two-dimensional stadium billiard, we will distinguish two specific types of periodic orbits in the three-dimensional square cylindrical stadium. The first type consists of planar periodic orbits, which are analogous to the bouncing-ball periodic orbits: they bounce from one side of the cavity to the opposite one, within the plane of the billiard. The second type consists of helical orbits, which whirl around the surfaces of the cylinders, going smoothly from one cylinder to the next across their intersection. As we shall show, the analysis of these two classes of periodic orbits is sufficient to determine the parametric regions where stable oscillations in cylindrical stadium billiards occur.
The specificity of cylindrical billiards is that all their surfaces have a flat component. This has the consequence that every periodic orbit has at least two parabolic eigen-directions, i.e. whose corresponding eigenvalues are unity. These are associated to motion parallel to the axes of the cylinders.
Since the squared cylindrical stadium gives rise to a four-dimensional symplectic billiard map, a periodic orbit also has another pair of eigenvalues, in addition to these two parabolic ones, whose product is also unity. This pair is associated to the motion in planes transverse to the cylinders' axes. It can be either hyperbolic, when the two eigenvalues are real and are inverses of each other, or elliptic, when they are complex conjugates with unit modulus. The case with four parabolic eigenvalues is irrelevant.
The consequences of the existence of a pair of parabolic eigenvalues for each periodic orbit are twofold. First, a periodic orbit is unstable when the second pair of eigenvalues is hyperbolic, but only marginally stable when this pair is elliptic. In other words, a given periodic orbit can never be linearly stable by itself. The second consequence is that periodic orbits exist in continuous one-parameter families, which can be parameterised by a displacement along the cylinders' axes. We refer to this displacement as the dynamical parameter, which we will denote in general by ε, and which is defined suitably to parameterise each family of orbits. It is not to be confused with the height of the cylinders, which is the model's geometrical parameter h.
If one were to limit the stability analysis to linear terms, as has previously been the case, then the naive conclusion would be that a family of orbits is stable on condition that every single orbit of this family admits a pair of elliptic eigenvalues, and provided that there are no geometric constraints which prevent this family from being complete because of a discontinuity. The reasoning is that the parabolic eigen-directions behave neutrally, so that given a periodic orbit and a small perturbation along these directions, the trajectory will move along them as it oscillates in the plane spanned by the elliptic eigenvectors, sweeping through the whole family of associated periodic orbits, parameterised by different values of ε, but left unscathed in the absence of a destabilizing mechanism. If, on the other hand, the family of periodic orbits undergoes a transition from elliptic to hyperbolic regimes as ε variesas happens in most cases -then the initially stable oscillations in the elliptic regime would become unstable at the bifurcation, when a value of ε is reached at which the family crosses over from the elliptic to the hyperbolic regions, and thus the motion would destabilise.
As it turns out, when elliptic regions occur for a given family of periodic orbits, transitions from elliptic to hyperbolic regimes are always found as soon as the value of the maximal height 2 + h of the cylinders in the elementary cell is larger than the cylinders' radii, i.e. when h > −1. The reason is essentially that beyond a bifurcation value of the dynamical parameter associated to a given family of periodic orbits, the distance between successive collisions is, at some point on the orbit, large enough to induce defocusing. According to the above arguments, one would then be led to believe that no stable oscillations could persist for square cylindrical stadia with heights larger than this bound.
However, interestingly, the story does not stop there. Extending the stability analysis to the next order beyond linear order, we find that, under some assumptions on the initial conditions, nonlinearities are sometimes able to stabilise a (partial) family of periodic orbits in a limited region of the dynamical parameter ε, leaving the orbit confined to the elliptic region, oscillating back and forth between two values of this parameter. Remarkably, the nonlinear perturbations are quadratic in the coordinates associated to the elliptic eigenvalues. As we show below, the key to understanding whether there exists a nonlinear stabilizing mechanism lies in a qualitative analysis of the eigenvalue spectrum of these quadratic forms: the signs of their eigenvalues tell us whether or not the oscillations in the plane spanned by the elliptic eigenvectors are able to stabilise the oscillations along the parabolic directions.
As our analysis will demonstrate, the conditions for stable oscillations are met only for a restricted set of values of the geometrical parameter h, consisting of a disjoint union of intervals and isolated values. For values of the parameter in that set, the phase space is typically mixed -different stable regimes may coexist, together with unstable ones. On the other hand, away from these values, the square cylindrical stadia are fully chaotic.

Planar orbits
By planar orbits, we mean those that remain in the plane of the billiard (x = 0), corresponding to all the orbits whose initial conditions are in that plane, with velocity components along the yand z-axes only ‡ For h ≤ 0, the periodic orbits in that plane are readily identified as the periodic orbits of a square billiard with sides of length 2(2 + h). Given the relative prime integers m and n, the point q = {0, y, z} is periodic provided the velocity is given by p = {0, cos ω m,n , sin ω m,n }, where ω m,n = arctan(n/m). Over a period, the corresponding orbit makes exactly m collisions on each of the horizontal walls, and n on each of the vertical walls, in a time τ m,n = 4(2 + h) √ m 2 + n 2 max(m/n, n/m). Furthermore, for any pair of relative primes {m, n}, there is a continuous family of orbits, indexed by these integers, obtained by varying the position of the collision impacts on the square sides. Still assuming h < 0, this family extends over all the possible values of y and z, For h > 0, the periodic planar orbits are those of a two-dimensional square billiard of side length 2(2 + h) with a square of side length h removed in the centre. The periodic orbits of this billiard include some that do not collide with the central square, which are common to the h ≤ 0 case, as well as others that do collide with the central square and which are in general more difficult to classify.
We will call (m, n) 0 planar periodic orbits (PPO) those periodic orbits of the plane that make collisions only with the outer surfaces of the cylinders. Examples are shown in figure 5. Except for n = 0 and m = n = 1, we will assume m < n.

. Existence
The simplest class of planar orbits are period-2 orbits which bounce back and forth between two opposite cylinders with velocity perpendicular to the cylinder axes. Assuming they propagate along the y axis, we identify them by the phase point with Cartesian coordinates These orbits exist everywhere along the axis of the cylinder so long as Cartesian coordinates. It is useful to go through the stability analysis explicitly. To do so, we compute the Jacobian of the periodic orbit by following the transformation rules of the tangent vectors δ q i and δ p i along the orbit, where i = x, y, z and each of these vectors has six dimensions. These transformation rules are given according to the two phases of propagation and collision by [34] wheren is the unit vector normal to the surface, and δ n is a triplet of six component vectors which depends on the curvature of the surface, Given a cylinder with axis in the yor z-directions and a trajectory in the plane of the billiard (the x components of the positions and velocities vanish), this reduces to δ n y = δ n z = 0 .
The orbit starting at the coordinates (2) bounces back and forth from q 0 , p 0 to Correspondingly, letting τ = 2(2 + h) denote the time between successive collisions, the tangent vectors are mapped according to equations (3)-(4) as follows, following the propagation (prop.) and collision (coll.) phases: This matrix has six eigenvalues. Two pairs are trivial, the first corresponding to the conservation of energy and the associated time-translation symmetry, and the second being the pair of parabolic eigenvalues which are due to the flat components of the cylinders' surfaces. The only pair of non-trivial eigenvalues of this matrix are then They are elliptic for −2 < h < −1 and hyperbolic otherwise, as expected. We note that a similar conclusion (with simpler eigenvalues) can be reached by considering only half of the orbit. This is because the corresponding periodic orbit in the elementary cell involves only a single collision with the cylindrical surface. However, because the parity of the number of reflections on the transverse planes can be odd in the elementary cell, considering the periodic orbits in the expanded billiard is better suited to the non-linear stability analysis discussed below.
Birkhoff coordinates. The symplectic structure of the collision map is best recovered by considering, instead of Cartesian coordinates, the Birkhoff coordinates of the billiard map, i.e. the coordinates on the constant energy surface in which the map is volume-preserving. They are given by: (i) θ , the position angle, measured in the plane transverse to the cylinders; (ii) ξ , the sine of the associated velocity, measured with respect to the normal to the cylinder's surface; and (iii)-(iv) the pair z, w of the position and velocity measured along the cylinder axis: In terms of the Birkhoff coordinates, the linear map (9) for the perturbations δ θ , δ ξ , δ w and δ z reduces to  It turns out that the structure of this linear map is common to all planar periodic orbits; the only non-trivial coefficients are the four matrix elements corresponding to the motion in the θ -ξ plane and the element of axial velocity contribution to the displacement along the cylinder axis. In the case of the (1, 0) 0 PPO, the matrix of elements m αβ is given by The four eigenvalues of this matrix form the pair of parabolic eigenvalues and the pair λ 1,2 obtained in equation (10). This result is summarised in figure 7, where the parameter space of the (1, 0) 0 PPOs is colored, when the orbit exists, in white where the eigenvalues (10) are elliptic and gray where they are hyperbolic, and in light red where the orbits do not exist. Similar conventions will be used throughout the paper. The conclusion we infer from this calculation is that the parameter region −2 < h < −1 corresponds to a regime with mixed phase space, where the period-2 orbits in the plane of the billiard give rise to an elliptic island. Though each of these orbits is only marginally stable, it belongs to a continuous family of similar orbits, which is moreover complete, meaning that it extends from one corner of the billiard to the opposite one. The family of orbits therefore corresponds to a plane of elliptic stability.
We note that the (1, 0) 0 PPOs are peculiar in the sense that their nonlinear stability analysis yields non-trivial corrections only to cubic order in the perturbations δ θ , δ ξ and δ w. Turning to the class of (1, n) 0 PPOs, we will see that, not surprisingly, the (1, 0) 0 PPOs are not the only planar periodic orbits which give rise to a stable regime. In fact, a perturbation of the (1, 0) 0 PPOs in the plane of the billiard will typically be closer to a family of (1, n) 0 PPOs with n large. Interestingly, such periodic orbits may exhibit a mechanism for non-linear stability quadratic in the perturbations δ θ , δ ξ , δ w, absent for the (1, 0) 0 family.

5.2.
(1, 1) 0 -planar periodic orbits 5.2.1. Existence The orbit with unit speed at angle π/4 with respect to the cylinder axes is the period-4 orbit which cycles through the phase points where

Linear stability
Going through the analysis described in section 5.1.2, we obtain a linear mapping describing the evolution of the perturbation vector along the periodic orbit, which, in Birkhoff coordinates, is identical to equation (12) with matrix elements The two non-trivial eigenvalues are The regions of elliptic and hyperbolic eigenvalues are displayed in figure 8. There are two separate elliptic regions, bounded by h = −1 ± ε and h = −1 ± √ 1 + ε 2 . For h = −1 and 2/3 ≤ h < 2, there are no elliptic regions along the ε axis and, where they exist, the (1, 1) 0 PPOs are everywhere hyperbolic for these values of the geometric parameter.
The transitions from elliptic to hyperbolic regimes are understood as follows. The segment of orbit between, say, q 1 and q 2 lies in the transverse plane z − y = 2 + h − ε, which cuts the cylinders along the yand z-axes in the form of two partial ellipses, glued one against the other. The semi-minor axis of these ellipses is unity (the radius of the cylinders) and the semi-major axis √ 2. The segment of orbit that joins q 1 to q 2 goes along the major axes of these two ellipses. The orbit is thus stable only so long as the distance traveled along the major axes is less than the radius √ 2/4 of the disk inscribed into the ellipses at q 1 and q 2 . This condition yields the critical values ε = ±(1 + h), which separates the hyperbolic and elliptic regions in figure 8. Our analysis shows that the family of (1, 1) 0 PPOs under consideration is linearly stable and complete only for −2 < h ≤ −3/2. For larger values of h, though some orbits display elliptic behavior, they do not correspond to linearly stable oscillations in the vicinity of the periodic orbit. The reason, according to the linear stability analysis, is that a trajectory initially close to an elliptic PPO will oscillate in the planes transverse to the cylinder's axes, and simultaneously move along the z-axis due to perturbations along the w-axis, which remain unchanged under iterations of the map (12), whose coefficient m zw , given by equation (15e), depends solely on h. In other words, the orbit moves in the parameter space of the (1, 1) 0 PPOs along the ε-axis due to the action of the perturbation along the vertical velocity axis. In the absence of nonlinear effects acting on δ w, the trajectory would eventually have to cross over from the elliptic to the hyperbolic region and therefore lose stability as the parameter ε translates away from the interval of elliptic eigenvalues.

Nonlinear stability
Linear stability would have us conclude that periodic orbit families whose non-trivial eigenvalues display bifurcations from elliptic to hyperbolic regimes must correspond to unstable oscillations, even when the reference periodic orbit we initially perturb away from is marginally stable. This, however, is in general incorrect. Going beyond linear stability, we find that the perturbations along the w and z coordinates are nonlinear functions of the oscillations that take place in the θ -ξ plane. At second order in the stability analysis, the θ -ξ oscillations give rise to quadratic terms which act on the perturbations along the parabolic directions according to where Q w (δ θ , δ ξ ) and Q z (δ θ , δ ξ ) are quadratic forms in their variables, which we write as The coefficients a w , b w , c w , a z , b z and c z are functions of the two parameters h and ε. Now, if δ w is large enough to start with, then the quadratic terms in (17a)-(17b) will do little to affect the motion along the w and z axes. Thus the problem we need to address is the following: assuming § that δ w is small with respect to quadratic terms in δ θ and δ ξ , we must find out whether the quadratic forms Q w and Q z are or are not able to induce oscillations of δ w and δ z, i.e. ε. Such oscillations would prevent the contribution m zw δ w in (12) from inducing a translation in the parameter ε (along the z axis) and would thus keep the trajectory oscillating within the elliptic region. In this case the corresponding family of (1, 1) 0 PPOs would remain stable as a result of the nonlinear oscillations.
Fortunately, the solution to this problem is simpler than it may seem, for, since the coefficient m zw in (15e) is positive, a mechanism inducing oscillations of w around ε = 0 (or, for that matter, around any fixed value ε = ε 0 ) exists whenever Q w has the opposite sign from ε (or ε −ε 0 if ε 0 = 0), i.e. whenever Q w is positive where ε is negative and Q w is negative where ε is positive. These oscillations in ε will be further amplified by Q z , constructively so provided Q z also takes the opposite sign from ε. If, on the other hand, Q w is positive where ε is positive and is negative where ε is negative, then, whether or not Q z keeps the opposite sign from ε, δ w will keep increasing or decreasing, steadily if Q w has the sign of ε, or on average if it has the opposite sign from ε, until the trajectory reaches the limits of the elliptic region and thus destabilises, either because of a transition to a hyperbolic regime, or because the limits of the range of allowed ε values is reached.
The quadratic forms (18a)-(18b) are positive-or negative-definite whenever the eigenvalues of their defining matrices have the same signs. (Since these matrices are symmetric and have real entries, their eigenvalues are always real.) Thus, irrespective of the oscillatory motion of the θ -ξ variables, Q w (resp. Q z ) is positive when its two eigenvalues are positive and negative when they are negative. When the two eigenvalues are of opposite signs, on the other hand, we must bear in mind that the oscillations of the θ -ξ variables keep going at their own pace, and they will probe the two eigendirections of Q w as they do so. Thus, as Q w (resp. Q z ) oscillates, it will be expected to take, on average, the sign of the eigenvalue that has the largest magnitude. Figure 9 shows the eigenvalue spectrum of the matrices of the quadratic forms (18a)-(18b), superimposed on the elliptic regions found from the linear stability analysis in figure 8. The coefficients of these quadratic forms are given explicitly in Appendix A.1, equations (1.1a)-(1.1c) and (1.2a)-(1.2c). For convenience, we analyze them by using a representation in four different colors, in terms of the sign and magnitude of their eigenvalues, which are § Identifying the initial position along the z axis with the dynamical parameter ε, we can always set δ z = 0. where both eigenvalues are positive (λ 1 λ 2 > 0, λ 1 + λ 2 > 0), red where both are negative (λ 1 λ 2 > 0, λ 1 + λ 2 < 0), green where one is more positive than the other is negative (λ 1 λ 2 < 0, λ 1 + λ 2 > 0), and cyan where one is more negative than the other is positive (λ 1 λ 2 < 0, The separation between nonlinearly stable and unstable elliptic regions at h = −1 is supported by numerical results such as those shown in figures 10 and 11.

5.3.
(1, 1) k -planar periodic orbits 5.3.1. Existence Orbits of the plane with unit ratio between the y and z velocity components are peculiar in that they remain periodic even when, for h > 0, they hit the inner walls of the cylinders, i.e. at y, z = ±h. This is generally not so of (1, n) PPOs, which typically cease to exist after a collision with the inner walls of the cylinders.
We define the (1, 1) k planar periodic orbits to be those among such orbits which make exactly k collisions with the inner walls of each cylinder over a periodic cycle. Examples are shown in figure 12. We identify these orbits by the initial coordinates    with the corresponding velocity  The values of the geometrical parameter, h, for which these orbits can be observed are restricted to the interval h min < h ≤ h max , whose bounds correspond respectively to the value of h for which the orbit hits the corners at the intersection of outer cylinders walls (in which case the segments of the orbits have constant lengths), and to the value of h for which the longest segment grazes the corner at the intersection of the inner cylinder walls: Given a value of h in this interval, the dynamical parameter, ε, which characterises the orbits of this family, is itself restricted to the interval ε min ≤ ε < ε max , where Interestingly, the (1, 1) k+1 orbit at h = h min is identical to the (1, 1) k−1 orbit at h = h max , except for eight segments which connect the inner to the outer corners back to back.  The results of this analysis are shown in figure 13. The patterns are very similar to those displayed in figure 9 for the (1, 1) 0 PPOs, with the conclusion that elliptic regions recur for every value of k at geometrical parameter values of h near h = 2k. To be more precise, the elliptic regions come in the shape of two tongues, elongated along the lines h = 2k − 1/(2k + 1) ± ε, which are symmetric reflections of each other along the line h = 2k − 1/(2k + 1). Moreover, the lower tongues are nonlinearly stable around ε = 0 in the interval 2k − 2/(2k + 1) ≤ h < 2k − 1/(2k + 1) and the upper tongues unstable.  (14), and velocity

Stability
These orbits exist for geometric parameter values between the bounds Given a value of h in this interval, the dynamical parameter, ε, is restricted, for n odd, to the interval ε min ≤ ε < ε max , where and, for n even and h < 0, to the interval ε min ≤ ε < ε max , where or, for h ≥ 0, to the two intervals ε min ≤ |ε| < ε max , where We note that, as n increases, the elliptic areas fill out the whole area −2 < h < −1, −2 − h < ε < 2 + h, resembling the pattern observed for the (1, 0) 0 PPOs shown in figure 7. This is not surprising, since, as n increases, the (1, n) 0 PPOs can be viewed as perturbations in the plane of the billiard of the the (1, 0) 0 PPOs.

Stability
The interesting feature of the (1, n) 0 PPOs is that, as opposed to the (1, 0) 0 PPOs, their nonlinear stability analysis is not trivial. As with the linear stability analysis, we proceed in a way similar to section 5.2.3 to obtain quadratic corrections to the linear stability map (13), similar to equations (17a)-(17b) and (18a)-(18b). We provide in section Appendix A.2 the explicit coefficients of the quadratic forms associated to the (1, 2) 0 PPOs, see equations (1.4a)-(1.4c) and (1.5a)-(1.5c). The expressions coefficients of the quadratic forms associated to the (1, n) 0 PPOs, n ≥ 3, will not be given explicitly.
The results of the nonlinear stability analysis are displayed in figure 15 and suggest that every elliptic strip may give rise to nonlinearly stable oscillations, at least over some range of values of the geometric parameter h. It is also interesting to note that the coefficients of the quadratic forms are all identically zero at the n values of h which correspond to the straightline lower bounds of the tongues of elliptic stability of these orbits.
Examples of such stable oscillations are shown from actual trajectories in figure 16, including the example of nonlinear oscillations around a (2, 5) 0 PPO. This orbit, as well as other (m, n) 0 PPOs, display features similar to those described above for the (1, n) 0 PPOs, and will we not go into further details.

Helical orbits
We now proceed to analyse orbits of helical shape, which wind around each cylinder in the squared billiard in turn. By symmetry, helical orbits of any pitch may occur along the surfaces of infinitely long, straight cylinders. When these orbits cross the intersection between two cylinders, however, they typically lose focus and go astray as they hit a corner between the two cylindrical surfaces. Nonetheless, there are isolated points where the intersection of the surfaces of the two cylinders is smooth, and a helical orbit will thus be transmitted unscathed if it passes through these specific points. One may further expect these orbits to generate stable oscillations in their immediate vicinity. Indeed, so long as the orbits remain close enough to the surface of the billiard, with a short separation between successive collisions, no defocusing can take place. We emphasise that this phenomenon is specific to higher-dimensional billiards with curved surfaces. The problem of assessing the stability of such orbits is again a fully nonlinear one which requires a suitable treatment.
In the case under consideration, where the cylinders' axes are in the x = 0 plane and intersect at right angles, the points of smooth intersections are the points x = ±1, where the curvature with respect to the other y or z coordinate vanishes. Thus a helical orbit which crosses over from one cylinder to the other at x = ±1, y = ±(1 + h) and z = y, or z = −y, is a smooth periodic orbit of the squared cylindrical stadium. (The conditions must be verified at every intersection.) These are rather stringent conditions and impose restrictive symmetries on the orbits.
We characterise these helical periodic orbits in terms of the number of times that they whirl around specific axes. Let n 1 and n 2 denote the number of half periods the helices make about the cylinders along the y and z axes, respectively. By symmetry, n 1 and n 2 must be odd integers, and it is easily found that the (n 1 , n 2 )-helix exists only at h = √ n 1 n 2 π/2 − 1, with initial condition in the z = 0 plane at x = 0, y = h if (n 1 − 1)/2 is even, y = h + 2 if odd, and with initial velocity components ± n 1 /(n 1 + n 2 ) along the x-axis and ± n 2 /(n 1 + n 2 ) along the z-axis. Technically, these orbits are of infinite period and are therefore impossible to realise. However, they give rise to an infinite sequence of periodic orbits of finite periods which approximate them. We will denote these approximate helices by (n 1 , n 2 ) k -helical periodic orbits (HPO), meaning that each half period (of which there are 2(n 1 + n 2 ) along the orbit) is made out of k segments. Examples of such long periodic orbits are shown in figure 17; see also related media . An important point is that, irrespective of their stability, finite-period approximations to helical periodic orbits can be moved around in some neighborhood, so long as they pass through the intersection between two cylinders away from their boundaries and do not hit a sharp angle. Finite-period approximations to a helical periodic orbit thus belong to oneparameter continuous families of such orbits, indexed by the dynamical parameter ε which we define below.
We now turn to the properties of these orbits and analyze their stability for specific classes.
6.1. (1, 1) k helical periodic orbits 6.1.1. Existence Finite-length (1, 1) k (k ∈ N, k ≥ 2) helical periodic orbits which approximate the (1, 1) helix cross the plane intersecting pairs of cylinders perpendicularly ¶. Their velocity at the intersections thus has the form (0, 1 √ 2 , 1 √ 2 ). The complete orbits can be identified as going through the point with Cartesian coordinates + :  (28) with the corresponding velocity The values of the geometrical parameter, h, for which these orbits can be observed are restricted to the interval h min ≤ h < h max , whose bounds correspond, respectively, to the value of h for which the lengths of the segments of the orbit going through the intersections of the cylinders shrink exactly to zero, and to the value of h for which these segments have their maximal length, grazing the surface of the cylinders at the intersection, Note that, for k large, the width of this interval shrinks to zero like 1/k, with lim k h max −h min 2π/k. The limiting value lim k h max = lim k h min = π/2 − 1 is the only value of the geometrical parameter for which the (1, 1) helix exists. Given a value of h in this interval, the dynamical parameter, ε, which characterises the orbits of this family, is itself restricted to the interval ε min ≤ ε < ε max , where As with h, these bounds correspond to the values of ε for which the periodic orbit either hits a corner at the intersection between two cylinders, or, the other way around, collides with the surface of the cylinders in the middle of its longest segment connecting the orbit between two transverse cylinders. Examples are displayed in figure 18. ¶ This property is shared with other symmetrical (n, n) k HPOs. + Although we list only one, there are in total four such similar orbits, which are easily obtained from each other by symmetry. This property is shared by all helical periodic orbits of any length. Figure 18. Families of (1, 1) k HPOs found at h = 0.5. As opposed to the exact helix, which is unique, its finite-length approximations are parameterised by a continuous parameter in the interval (31a)-(31b). We show them for different values of the parameter taken across the interval, for (a) k = 2, As opposed to the parameter h which fixes the geometry of the billiard, the parameter ε is a dynamical one: it changes as perturbations of the periodic orbit are introduced. Just as was the case with the planar periodic orbits, displacements of the helical periodic orbit along the z axes as the trajectory returns periodically close to its initial condition (28) can be interpreted as changes in the value of the parameter ε; the trajectory comes back closer to another periodic orbit than the one it was perturbed away from. In particular the trajectory moves away from this family of orbits as soon as the value of ε exceeds either of its two bounds, ε min or ε max , since it ceases to exist beyond them.

Stability
Birkhoff coordinates. To analyze the stability of the (1, 1) k orbits, we change coordinates to the Birkhoff coordinates (11a)-(11d), here given by Linear stability analysis. Introducing perturbations δ θ , δ ξ , δ z and δ w, the linear stability analysis yields the symplectic map  This form is similar to that of the planar periodic orbits, (12), but with extra non-zero coefficients m zθ and m zξ . Its eigenvalues remain similar: there are two unit eigenvalues associated to the motion in the z-w plane, and two nontrivial eigenvalues, η 1 and η 2 , associated to the motion in the θ -ξ plane. These two eigenvalues can be either hyperbolic, when they are real, η 1 = η −1 2 , or elliptic, when they have a non-zero imaginary part, η 1 = η * 2 . The matrix elements m αβ are generally complicated functions of the two parameters h and ε.
The cases k = 2 and k = 3 are the simplest and yield results shown in figure 19, exhibiting features rather similar to the (1, 1) k PPOs, which are common to all (1, 1) k HPOs. (Explicit values of the coefficients m αβ in the linear form (33) for these two cases are provided in Appendix B.) Namely, there are regions in the ε-h parameter plane where η 1 and η 2 are elliptic, which are in the shape of two symmetric tongues. For k = 2 the symmetry line corresponds to h = 3 √ 2/4 − 1, and to h = 1/3 for k = 3. This value keeps growing for larger values of k towards h = π/2 − 1.
Here, as with the (1, 1) k PPOs, we see that the coefficient m zw in the linear form (33) depends only on h -see equations (2.1g)-(2.3g) in section Appendix B. Repeating the arguments given in section 5.2.2, we would conclude that all these periodic orbits are linearly unstable because of the existence of a bifurcation from elliptic to hyperbolic regimes in the range of dynamical parameter values the orbit would visit as it translates further away from its initial position under the action of δ w. Nonlinear stability analysis. There are, however, again nonlinear effects which may affect this scenario, due to terms which are quadratic in the small perturbations. The picture here is actually slightly simpler than with the planar periodic orbits. Indeed, since the θ , ξ and z variables pick up nontrivial perturbations at the order of linear contributions, it is enough to perform the quadratic expansion for the w axis alone. This yields the quadratic map where Q w (δ θ , δ ξ ) is a quadratic form, given by similar to (18a). The coefficients a w , b w and c w are here again functions of the two parameters h and ε.
Recapitulating the arguments given in section 5.2.3, the sign of the quadratic form (35) is determined by that of the eigenvalues of its matrix. When the eigenvalues are of the same sign, Q w will be positive when they are positive and negative when they are negative. If, on the other hand, the eigenvalues are of opposite signs, Q w will oscillate, but, on average, will take the sign of the eigenvalue that has the largest magnitude. Stable nonlinear oscillations can take place in the parameter regions where the eigenvalue spectrum of Q w changes sign so as to confine the small oscillations to an interval of the dynamical parameter where the eigenvalues of the linear form (33) are elliptic. Figure 20 shows, for the k = 2 and k = 3 HPO families, the eigenvalue spectrum of the matrix of the quadratic form (35), superimposed on the elliptic regions found from the linear stability analysis in figure 19. The coefficients of these quadratic forms are given in Appendix B. The color convention is the same as the one used in section 5.2.3, in terms of the sign and magnitude of their eigenvalues, which are denoted by λ 1 and λ 2 . We observe that the elliptic regions below and above the lines at h = 3 √ 2/4 − 1 for k = 2 and h = 1/3 for k = 3 display clearly different patterns. Whereas the upper elliptic regions are everywhere unstable under the quadratic perturbations, the lower elliptic regions contain stable regions around ε = 0 and h above h = −1+1/ √ 2 for k = 2 and h = 1/6 for k = 3 (these bounds correspond to the lowest point of intersection between the elliptic and hyperbolic regions at ε = 0). Though the lower bound we read from this analysis is sharp, the upper bound appears perhaps less so. More precise results might be obtained from a more detailed analysis of the eigenspectrum of Q w , but this will not concern us here. Numerical results, such as shown below, support the claim that the upper bound is actually sharp too. The main result here is the existence of a region of parameter space where the helical periodic orbits under consideration are nonlinearly stable.
The stability can be checked numerically by perturbing actual trajectories about these periodic orbits, as shown in figure 21.
Similar results, obtained for larger values of k, are shown in figure 22, there resorting to a numerical computation of the quadratic form Q w on a grid of points of the parameter space. That is to say, we fix numerical values of the two parameters h and ε, identify the corresponding periodic orbit, and perform a perturbation analysis around that orbit, which yields numerical values of the coefficients m αβ in equation (33) and a w , b w and c w in equation (35). As k gets larger, the regions of stability get thinner, but they persist and can be observed. For example, we were able to check the existence of nonlinearly stable periodic orbits for the family of (1, 1) 101 HPO at h = 0.5707, which belong to class of orbits displayed in figure 17(a). The same orbit is, however, found to be unstable at h = 0.5708.  6.2. (n, n) k Helical Periodic Orbits 6.2.1. Existence We can identify, for arbitrary odd integer n, any finite length (n, n) k HPO approximating the (n, n) helices, with the same initial conditions (28)- (29). These are basically the same orbits, only the number of collisions along the periodic orbit changes from 2k to 2nk and the values of the parameters (30a)-(30b) and (31a)-(31b) are accordingly changed to

Stability
The problem of determining the stability of (n, n) k HPOs transposes verbatim from section 6.1.2. The computation gets more challenging as the lengths of the orbits increase, but similar results can be obtained. This is supported by the results shown in figure 25, obtained for the (3, 3) 2 HPO. 6.3. (n 1 , n 2 ) k Helical Periodic Orbits 6.3.1. Existence. Given odd integers n 1 and n 2 , n 1 = n 2 , the (n 1 , n 2 ) k HPOs which approximate the (n 1 , n 2 ) helical periodic orbit can be identified starting with initial positions (28) and velocity components adapted according to the ratio between the numbers of half periods to be completed along the two axes y and z: Because of the asymmetry between the y and z motions, there is here, and for any given k, a unique value of the geometric parameter at which these orbits are realised, The corresponding interval of values of the dynamical parameter, ε, has the bounds ε max = 2 n 2 /n 1 sin[π/(2k)], k even, ε min = −2 sin[π/(2k)], k even, −(1 + n 2 /n 1 ) sin[π/(2k)], k odd. (40b) We will focus our attention on three groups of such orbits, namely the (1, 3) k , the (1, 5) k , and the (3, 5) k HPOs. Specific examples of families of these classes of (n 1 , n 2 ) k orbits are shown in figures 26, 27 and 28, respectively.

Stability
The nonlinear stability analysis proceeds along the lines described in section 6.1.2. However, since each (n 1 , n 2 ) k HPO exists only for a singular value of the geometric parameter, we can consider the problem of identifying the elliptic regime of the map (33) and the quadratic form (34)-(35) as a function of the dynamical parameter ε alone. This computation can be carried out numerically for moderate values of the overall period. The results are displayed in figures 29, 30 and 31, respectively, where we plot the determinant (=λ 1 λ 2 ) and trace (=λ 1 + λ 2 ) of the quadratic form throughout the range of dynamical parameter ε values where the eigenvalues of the linear stability analysis are found to be elliptic. In all of these cases, we observe that the elliptic regime of dynamical parameter values is asymmetric with respect to ε = 0 for k even and symmetric for k odd. In these plots, candidate stable HPOs would correspond to values of ε for which both the determinant and  figure 29, we observe the opposite: namely, the points where both the determinant and the trace of Q w vanish correspond to positive values of the derivative of the trace. Therefore they are unstable, since oscillations around the given HPO will tend to grow away from it. Similar conclusions are drawn for the (3,5) k HPOs whose stability analysis is shown in figure 30.
The case of the stability analysis of the (1, 5) k HPOs, shown in figure 30, is different. There are now two points where both the determinant and trace of Q w vanish. The leftmost one is of the same kind as the one observed for the (1, 3) k HPOs. However, the rightmost one is certainly a stable point, since it appears to correspond to a minimum of the determinant of Q w . It must, however, be noted that it arises near the border of the bifurcation point of ε which separates the elliptic region from the hyperbolic one. Stable oscillations around this point must therefore be confined to a very small region of phase space. We have not been able to obtain conclusive evidence of their existence.

Summary and perspectives
The flat surface component of three-dimensional cylindrical billiards complexifies the stability analysis of their periodic orbits. In the presence of two elliptic eigenvalues which generate linear oscillations in the plane transverse to the cylinder axis, nonlinear effects along the neutral directions may act as a restoring force, allowing for stable oscillations in all phase space directions. Summarizing our results, the three-dimensional cylindrical stadium billiard with a single oblique plane at angle π/4 with respect to its axis displays a large range of different dynamical regimes, stable and unstable, which can be analysed in terms of the spectral properties of families of planar and helical periodic orbits of this billiard. We classify these regimes in terms of the geometry of the expanded square stadium billiard, which is parameterised by its height: −2 < h < −1 Stable oscillations can be observed in large regions of phase space when the height is small enough that opposite circular arcs are closer than their diameter. These oscillations take place near a number of families of stable planar periodic orbits. Nonlinear effects act on these oscillations as a restoring force which for some of these families may constrain the oscillations around the periodic orbit located at the center of the family. Though most of the periodic orbits are unstable in some regions of the parameter space in this range, others remain stable so that there are no parameter values in this region which correspond to a fully chaotic regime.
−1 < h < −1 + 1/ √ 2 All but the (1, 1) 0 planar periodic orbit are unstable in this parameter range. Though the (1, 1) 0 PPOs still display regions of elliptic regimes, they are destabilized by nonlinear effects. Though the (1, 1) 2 HPOs exist in this range, they are unstable. This range therefore corresponds to a fully chaotic regime: the system is ergodic with all four Lyapunov exponents non-zero.
The first familiy of helical periodic orbits, the (1, 1) 2 HPOs, are nonlinearly stable in this range, around the central orbit at ε = 0. This gives rise to a small island of stability in an otherwise chaotic phase space.
The (1, 1) 2 HPOs are nonlinearly unstable in this range and the (1, 1) 3 HPOs are unstable. This is another range of parameter values for which the system is ergodic with all four Lyapunov exponents non-zero.
π/2 − 1 < h < 4/3 The system is ergodic with all four Lyapunov exponents non-zero in this range, which separates the last of the (1, 1) k HPOs from the region of nonlinear stability of the (1, 1) 1 PPOs.
A similar succession of fully chaotic and ergodic regimes, and mixed, non-ergodic, regimes where the phase space is separated into a single chaotic region and small elliptic regions, continues as we keep increasing the height parameter h. Thus the next regions of nonlinear stability correspond to the (3, 3) k HPOs which start for k = 2 near h = 3 and culminate at h = 3π/2 − 1, followed, in the interval 18/5 < h < 19/5, by the (1, 1) 2 PPOs, and so on. Let us also mention the possibility that stable regimes of oscillations exist for isolated values of h, such as with the (1, 5) k HPOs. However, we have not been able to confirm their existence numerically.
We note that the stable oscillation regimes described above are not typical of all billiards in the larger class of convex, cylindrical stadium billiards, as they can be easily destroyed by a change in the angle of the oblique plane or by the insertion of another oblique plane. In particular, these stable oscillations do not exist in the billiards obtained from interactingparticle models confined by hard walls. The three-dimensional cylindrical stadium billiard with two perpendicular planes at angle π/4 with respect its axis is thus ergodic and fully chaotic in its whole range of parameter values. It is indeed easy to see that there are no stable planar periodic orbits and that helical periodic orbits do not exist. A periodic orbit of this billiard always has at least one segment long enough to induce defocusing.

Appendix B. Stability analysis of the Helical Periodic Orbits
We provide below the expressions of the non-trivial coefficients of the linear stability matrix (33)