{'fraction_non_alphanumeric': 0.0333655346659731, 'fraction_numerical': 0.01252136434569369, 'mean_word_length': 4.125690344696248, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
About ergodicity in the family of limacon billiards

By continuation from the hyperbolic limit of the cardioid billiard we show that there is an abundance of bifurcations in the family of limacon billiards. The statistics of these bifurcation shows that the size of the stable intervals decreases with approximately the same rate as their number increases with the period. In particular, we give numerical evidence that arbitrarily close to the cardioid there are elliptic islands due to orbits created in saddle node bifurcations. This shows explicitly that if in this one parameter family of maps ergodicity occurs for more than one parameter the set of these parameter values has a complicated structure.


Introduction
In Hamiltonian systems ergodicity is a rather rare property [1]. Most systems are "mixed", in the sense that there are regions of regular motion of positive measure and regions with stochastic motion, which are expected to be of positive measure, too; see [2] for a review on the coexistence problem.
In this paper we investigate the possibility of ergodicity in the family of billiards introduced by Robnik [3], which is given by the simplest non-trivial conformal mapping of the unit circle; in polar coordinates the boundary of the family of limaçon billiards is parametrised by where ε ∈ [0, 1] is the family parameter. This family of billiards has received a lot of attention, especially in the context of quantum chaos, see e.g. [4,5,6], and the references therein. Corresponding to ε = 0 is the integrable billiard in the circle, while ε = 1 gives the ergodic billiard inside the cardioid [7,8,9]. For a small perturbation of the integrable case the KAM theorem implies that most invariant tori persist, whereas for the neighborhood of the hyperbolic case there is no general result about persistence of non-uniformly hyperbolicity. From this point of view one of the most interesting questions for this billiard family is whether it is hyperbolic or even ergodic for some ε < 1.
For ε ∈ [0, 1 2 [ the boundary of the billiards is convex, C ∞ -smooth and the curvature is strictly positive. As a consequence of a general theorem by Lazutkin [10] these billiards have a collection of caustics near the boundary, which correspond to rotational invariant circles in the Poincaré map from bounce to bounce. Since rotational invariant circles divide the Poincaré section into invariant regions these billiards are not ergodic. At ε = 1 2 the boundary has zero curvature at φ = ±π and therefore the invariant circles corresponding to the caustics near the boundary no longer persist [11]. For ε > 1 2 the billiards are no longer convex, and the billiard map becomes discontinuous. From these observations and from numerical experiments Robnik suggested that for ε > 1 2 the billiards might be ergodic [3], see also [12,13]. It turned out, however, that at least close to ε ≥ 1 2 this is not true because still a family of Liapunov stable periodic orbits exists [14]. This orbit is involved in a cascade of period doubling bifurcations, where the extrapolated limit is ε = 0.5582 . . . [14]. A first example of a stable periodic orbit above this value was given in [15] where a bifurcation is discussed for ε ≈ 0.66948.
One way to prove that the system is not ergodic is to establish the existence of an invariant curve. This is possible by using KAM theory in the neighborhood of an elliptic fixed point of the map. A sufficient condition for the existence of an invariant curve in the neighborhood of an elliptic fixed point is that the eigenvalue of the linearised map at the fixed point is not a third or forth root of unity and that the twist does not vanish, see [16]. In [14] this was rigorously established for a certain orbit of the limaçon family. In the following we will rely on numerical results in which we cannot easily check the twist condition. However, as in one parameter families the twist generically vanishes at points [17], the existence of an elliptic orbit for one given parameter value gives strong numerical evidence that the map is not ergodic at this parameter.
As for ε larger than the extrapolated limit of the period-doubling cascade there are no obvious stable islands in the Poincaré section, it might be conjectured that at this point all the orbits that are part of the bifurcation tree that is connected to the integrable circle billiard have turned unstable. But even if this was true (for which we have not found any counterexamples), it does not imply ergodicity of the system. The reason is that there exist stable orbits which are not connected to the integrable limit. They appear in saddle-node bifurcations, in which a stable and unstable orbit are created "out of nowhere", in principle at any parameter value. At each saddle-node bifurcation the newly created elliptic orbit exists in a certain parameter interval until it turns unstable in a period doubling bifurcation. We give numerical evidence that there are saddle-node bifurcations arbitrarily close to the ergodic limit at ε = 1. Moreover, our numerical computations show that the sum of the widths of the stable intervals of elliptic orbits born in saddle-node bifurcations has considerable measure in parameter space. For a related study of the phenomenon of the appearance of elliptic islands in systems close to ergodic billiards see [18] or in the standard map with large parameter see [19].
However, we cannot rule out that ergodicity holds for some values ε < 1, because the stable intervals of these orbits are so small that they do not seem to cover the whole interval [ 1 2 , 1[. The method we use is continuation of periodic orbits away from the hyperbolic limit. At this limit we know all periodic orbits. For the regular orbits, i.e. those orbits which never hit the cusp of the cardioid, a binary symbolic dynamics [5,20] was proven in [21,22], and for the cusp orbits a similar coding with an additional letter can be used. The idea to continue periodic orbits from hyperbolic a hyperbolic limit is similar to continuation from an antiintegrable limit [23], which has recently be applied to the Hénon map [24]. However for an unbounded system like the Hénon map, the map can actually be proven to be hyperbolic near its limit [25,26]. For a system with bounded phase space this is rather unlikely due to the conservative Newhouse phenomenon [27]. Moreover, the cardioid billiard is not uniformly hyperbolic. There is no positive lower bound for the Liapunov exponent of certain orbits. Nevertheless, for the cardioid every finite periodic orbit is hyperbolic and therefore can be continued. But due to the nonuniform hyperbolicity, the system can become non-hyperbolic when perturbed.

Finding and continuing orbits
Our approach is to continue periodic orbits of the cardioid billiard (ε = 1) to smaller values of the parameter ε. For the continuation we use the discrete variational principle, which is numerically much more stable than iterating the map, and also allows for the continuation of singular orbits. The discrete analogue of Hamilton's principle states that the action is stationary on orbits, where L is the discrete Lagrangian of the system. Performing the variation δW = 0 gives a second order difference equations. Assuming the twist condition holds the analogue of the Legendre transformation can be used to rewrite it as an area preserving map, see e.g. [28] for a review. For billiards the discrete Lagrangian L(x, x ′ ) is just the geometric length between the two points on the boundary parameterised by x and x ′ , so that W is the total length of the orbit. The Poincaré map from bounce to bounce of a smooth billiard is a map of a finite cylinder to itself. The global winding number of a periodic orbit is given by m n , where n is the number of reflections and m counts how often the orbit goes around the cylinder.
We use polar coordinates (ρ, φ) to describe the billiard boundary, so that the polar angle φ replaces x. With (1) the discrete Lagrangian becomes For a doubly periodic Lagrangian periodic orbits of period n are critical points of the periodic action where X n = (x 0 , x 1 , . . . , x n−1 ). In general the critical points of W p n are found using Newton's method applied to the system of equations ∇W p n = 0. In [20,22] it was shown that all regular orbits, i.e. those which never hit the cusp, of the cardioid are maxima of W p n , such that much more stable numerical methods can be used to find them.
Let us now briefly recall the definition of the symbolic dynamics for the cardioid billiard [5,20] which gives a complete description of all possible orbits in the cardioid, see [21,22] for proofs. Consider two consecutive points on the boundary with polar The method described in [22] requires to find initial conditions consistent with a given code word, for which we use the following procedure. Rotate the string until it starts with BA. In step 1 choose (φ 0 , φ 1 ) = (z, δ) with z > δ > 0 as the initial pair corresponding to the letter A. Here δ is a small positive number and z will be determined at the end. In step i take φ i−1 from the previous step, and construct the pair In the last step n we obtain by construction the letter B, which is the first letter of the word. We have to fulfill the conditions φ n−1 < z and z > δ, e.g. by the choice z = (max(φ n−1 , δ) + π)/2. Note that all φ i are in [−π, π] for arbitrary sequences of A and B. This shows that we can choose an admissible initial condition for each code word, i.e. there is no intrinsic pruning [22].
If the maximization of W for a given symbol sequence leads to a solution for which one or more segments run outside of the billiard, this periodic orbit is not realised by the system; the orbit is called "f-pruned" [20]. It may also occur that the solution possesses a point in the singularity, a situation in which we say the orbit is "s-pruned", or in the more general context of [22] extrinsically pruned. Also combined types of pruning occur, see [20] for examples of the different types of pruning.
In order to find all orbits for ε < 1 we also have to find finite orbits in the cardioid billiard, i.e. cusp orbits, which begin and end with an orbit segment in the cusp, since most of them turn into regular orbits for ε < 1. This might seem surprising at first; see below for a proof. The main difference between these singular orbits and the regular ones it that there is no reflection condition to be fulfilled for orbit segments in-and outgoing to the singularity of the cardioid at φ i = ±π because In order to find cusp orbits we can therefore look for critical points of an action with fixed initial and final point Note that W f n describes an orbit with n reflections not counting reflections at α and ω. It is sufficient to restrict to cusp orbits hitting the cusp just once, as these can be combined to multiple cusp orbits by concatenation. Similar to the case of regular orbits in the cardioid we can find cusp orbits by looking for maxima of W f n . The proof is as for W p n , the only difference is that L(φ n−1 , φ 0 ) is replaced by L(φ n−1 , ω) + L(α, φ 0 ). For the cardioid with α = ω = π the Hessian of this function is diag(− cos φ 0 , − cos φ n−1 ). Thus it is negative definite for orbits which hit the cusp from the right (both angles, before and after hitting the cusp, have modulus smaller than π/2) and the same proof as in [22] works. Unfortunately we have not been able to show that orbits with modulus of φ 0 or φ n−1 larger then π/2 are also minima. Even if such an orbit is initially pruned, e.g. A 2n CC for n ≥ 2, it can become physical under continuation. The shortest code for which we cannot find a corresponding minimum in the action is A 6 B 3 CC. However, using the symmetry line method described in [6] to find this orbit shows that it is in fact s-pruned.
The initial conditions for the numerical determination of the maximum of W f n for a given cusp orbit are obtained similarly to the case of regular orbits: first the code word is rotated so that it starts with CC. The same procedure as above is applied for intermediate φ i , except in the end we just put φ 0 = z = π. To follow orbits under parameter variation we use a predictor-corrector method, see e.g. [29]. For a slightly changed parameter ε + ∆ε an initial guess for the orbit is computed in the predictor step, which is used in the corrector step to find the actual orbit at the new parameter. We continue the orbit Φ = (φ 0 , . . . , φ n−1 ) under parameter variation in ε by expanding Because Φ is an orbit at parameter ε we have ∇W p n = 0. Hence we can solve for ∆Φ provided that D 2 W p n is nonsingular. This gives the predictor Φ + ∆Φ for the continued orbit at parameter ε + ∆ε. The corrector step consists in Newton's method taking the predictor as the initial guess for the orbit at the new parameter ε + ∆ε. Note that for continuing orbits away from the cardioid it is necessary to use Newton's method, because for ε < 1 not all the critical points of W p n are maxima. If the orbit cannot be found for the new parameter ε + ∆ε the step size is decreased. If under parameter variation det D 2 W p n = 0 is reached, we have found the point (in most cases, see Sec. 3) at which the orbit is created. In order to be able to continue orbits through this singular point in the parametrization by family parameter ε we consider the augmented problem, see e.g. [29, chpt. 4.4], in which the parameter is itself a variable, and the continuation uses the arc length along the curve in (Φ, ε) space as a new parameter.
A cusp orbit of length n has n − 1 points Ψ n−1 = (ψ 0 , . . . , ψ n−2 ) and one point at the cusp. Ψ n−1 satisfies the equation ∇W f n−1 (Ψ n−1 ; π, π) = 0 at ε = 1. To continue a cusp orbit away from the cardioid we make the fixed angle π in W f n−1 (Ψ n−1 ; π, π) an additional variable and extremise W p n (Φ n ) for ε < 1, where initially Φ n = (Ψ n−1 , π). Note that because of (5) we have ∇W p n (Φ n ) = 0 at ε = 1 as before. However, the Hessian of W p n is not definite at this point. For the continuation the important fact is that it is nondegenerate, which we will now show. First note that and that outside the cusp the second derivatives of W f n and W p n coincide. Therefore we find where ϑ = cos ψ 0 + cos ψ n−2 .
We have already shown that D 2 W f n−1 is negative definite, and we also noted that for physical orbits the cosines in ϑ are both positive. Therefore D 2 W p n is nondegenerate at cusp orbits. The difference to orbits that do not hit the cusp is that here the Hessian is not negative definite, instead it has one positive eigenvalue. In any case these orbits can be continued away from the cardioid. By obvious extension of the above the case where a finite orbit hits the cusp more than once can be included as well.
It is surprising that the continuation of cusp orbits away from the cardioid is not singular, because these orbits are "infinitely unstable" in the cardioid. This can be seen using a formula for the residue R of an orbit [30] with cyclic boundary conditions. Observe that for ε → 1 we have that B 0 → 0, B n−1 → 0 due to (8) and thus R diverges. The limiting behaviour can be found by expansion around ε = 1. This gives B i = (1 − ε) cos φ i for i = 0 and i = n − 1, so that we find near ε = 1. Note that for multiple cusp orbits a similar result holds with the power of 1 − ε decreased by one for each point in the cusp.

Bifurcations of orbits
To study the bifurcations of a periodic orbit we start in the cardioid and continue it from ε = 1 towards the circle ε = 0 as long as it exists. We include all cusp orbits and also all pruned orbits, because they typically become physical for smaller ε when the billiard becomes closer to being convex. An orbit can cease to exist if it collides with another one at residue R = 0, corresponding to an eigenvalue 1 of the linearised Poincaré map. For all other finite R orbits cannot disappear, which can be seen by applying the implicit function theorem to (7). This only holds for smooth maps, so the singular bifurcations which occur in a non-convex billiard when there is a tangency of an orbit segment and the boundary can occur with any residue. The variational approach is blind to these singular bifurcations, except if there is a reflection at the point of tangency. In this case the stability of the orbit diverges at the tangency. Numerically the continuation is done for decreasing ε. However, when describing bifurcations we always consider an increasing parameter. The following possibilities arise: • The orbit exists down to ε = 0. In the integrable circle billiard periodic orbits can be labeled by their rational rotation number ω = p/q. Under small perturbations the corresponding rational torus breaks into a Poincaré-Birkhoff chain consisting of one unstable and one stable orbit of period q. The existence of rotational p/q orbits for larger ε is guaranteed by Aubry-Mather theory [31,32] (see [28] for a review)    as long as the billiard is strictly convex, i.e. for ε < 1/2. The unstable p/q orbits maximise the action W , the (initially, i.e. for ε > 0, ε small) stable ones are called minimax, they are saddles of W with one up direction. The minimax orbits up to period 7 together with the ε value at which they turn unstable are listed in Table 1.
• The orbit only exists down to a nonzero parameter value ε c at which it has R = 0 and ∂R/∂ε = 0. In systems with a reversing symmetry generically one of the following bifurcations occurs; the orbit collides with 1. another orbit with opposite sign of R in a saddle node bifurcation (sn); there is no parent orbit (Fig. 1).
2. another orbit with the same sign of R in a symmetry breaking pitchfork bifurcation (pf); the parent orbit has R = 0 (Fig. 3). The subcritical pitchfork bifurcation in which the parent turns stable is denoted by upf (Fig. 4).
3. itself in a period doubling bifurcation (pd); the parent orbit has R = 1. The subcritical period doubling in which the parent turns stable is denoted by upd (Fig. 2).
4. another orbit with opposite sign of R in a p/q-tupling bifurcation, q ≥ 3; the parent orbit has R = sin 2 (πp/q). The "touch and go" bifurcation with ∂R/∂ε = 0 for q = 3 and possibly 4 does not terminate the orbit because it also exists for smaller ε.
The above cases exhaust the list of generic bifurcations with one symmetry [33]. The family of limaçon billiards has two reversing symmetries, X and T , see [20]. Therefore there are more generic bifurcations of symmetric orbits in this family [34,35,36]. E.g., the period tripling of an orbit with full symmetry creates four orbits with less symmetry whose phase portrait looks like a six-tupling bifurcation. We denote this bifurcation by s1/3 (Fig. 2). All bifurcations for periodic orbits up to period 7 are listed in Table 2. The parent orbit (which does not exist for sn) exists before and after the bifurcation. With the exception of the q-tupling bifurcations it changes its stability at the bifurcation, i.e. it turns from elliptic to hyperbolic (for pf, pd) or vice versa (upf, upd). If there are two partners in a bifurcation the stable one is listed first. In the last column of Table 2 Table 1: Minimax orbits with winding number p/q, the corresponding code in the cardioid billiard and the ε value for which they have R = 1 and turn unstable. The last column lists the ε value where the orbit has a tangency with the boundary and becomes unphysical for larger ε. parent Table 2: Table of bifurcations for orbits up to period 7 in the family of limaçon billiards. Since maximizing orbits do not bifurcate they are not included. Note the sequences of symmetric saddle nodes followed by a pitch fork bifurcation, see Fig. 3,4. At the end of the table there are the first two pairs of orbits born in a saddle node bifurcation without symmetry, see Fig. 1. For higher periods most orbits are created in this type of bifurcation.
stable interval is listed, where a negative number means that it extends to the left from the listed value of ε. For symmetric orbits usually one of the partners has a second interval of stability. These orbits are listed as a single partner without a bifurcation type together with the value ε for which R = 1. The same families will also be parents of pitch-fork bifurcations when R = 0. This is listed next, so that for most symmetric orbits there are three consecutive entries that describe the essentials of the corresponding bifurcation tree (Fig. 3,4). For orbits with full symmetry the two intervals of stability join at a point where the graph of R(ε) has a tangency with 1 (Fig. 2). In this way all stable intervals of orbits up to period 7 and their bifurcations are listed in Table 2.
The relation between the winding number p/q and the symbolic code for the maximizing and minimax orbits is simple. The code of the maximizing p/q orbit is given by the word with symbols where [x] denotes the integer part of x and i = 0, . . . , q − 1. The corresponding minimax orbit of the same rotation number is obtained by replacing the first and the last letter by C, see Table 1 for some examples. With these rules we can predict at ε = 1 which orbit is going to have a particular rotation number at the integrable limit. A typical simple saddle node bifurcation is shown in Fig. 1. Depending on the symmetry we observe more complicated sequences of bifurcations, see Fig. 2, Fig. 3, and Fig. 4 for examples. These are the four basic types of (sequences of) bifurcations that occur up to period 7, compare to the corresponding entries in Table 2.

Statistics of Bifurcations
We numerically computed all stable intervals of all orbits up to period 13. This is a complicated task that involves continuing 3697 periodic orbits. In particular it came as a surprise that the width of the stability intervals for some of the orbits could not be calculated using double precision, i.e. approximately 15 significant digits. E.g. already at period 10 the smallest interval of stability occurs for the saddle node pair AAACCBCCCC, AAAABBCCCC, and is only ∆ε ≈ 4.1 · 10 −18 . In order to be able to calculate such small intervals we used the doubledouble package [37], which offers 30 significant digits. This enabled us to determine all stable intervals of continued orbits up to period 13. Studying higher period would require even more precision. We now report on the results of our numerical study.
Up to period 7 the results are listed in Table 2. Fig. 5 shows the width ∆ε of intervals of stability as a function of ε (at the left boundary of the interval) for orbits of period 11 and for all orbits up to period 13. The number of points grows with the period like the number of orbits (including cusp orbits) in the cardioid, i.e. exponentially with topological entropy ln(1 + √ 2) [20]. Smaller intervals occur for larger ε. The trivial upper bound is ∆ε(ε) = 1 − ε. For fixed period a lower bound is approximately given by a straight line. There are conspicuous patterns that can be related to the structure of the code word. The detailed structure is quite complicated, and it is difficult to decide whether for each parameter there exist elliptic orbits. Even though the number of stable intervals grows exponentially with the period this does not mean that they are easily visible in phase space. In [12,13] it is referred to an unpublished numerical study by Li and Robnik   (1994), giving that if there is a stability island for ε > 0.555 its relative area has to be smaller than 5 · 10 −7 . This is in agreement with examples considered in [6].
To obtain a quantitative estimate of the behaviour of width of stability intervals with the period we now consider the average ∆ε of ∆ε for each period, see Fig. 6. Here it is important to treat orbits with and without symmetry separately because orbits with symmetry tend to have larger intervals of stability. A reason for this is that many symmetric orbits undergo symmetry breaking pitchfork bifurcations, see Fig. 3. Therefore they have to intersect TrM = 2 with finite slope, so that the second interval of stability is larger than the first. There is a difference between even and odd periods. On average nonsymmetric orbits with odd period have smaller stable width than those with even period. This is mainly due to the fact that there are no period doubled orbits among them, which tend to have a larger stable width because of their finite slope in the bifurcation diagram at creation (see e.g. Fig. 2). The average width of the parameter interval on which the orbits are stable follows a power law with the period. For the non-symmetric orbits the exponent is close to the topological entropy 0.881. The actual growth of the number of intervals of stability is shown in Fig. 7. Therefore the width of the intervals decreases at the same rate as the number of orbits increases, such that it is not possible from these estimates to determine the asymptotic behaviour of the sum of all intervals. For the symmetric orbits a approximately half the exponent is expected and measured.
We can conclude that the behaviour of this family of billiards is much more complicated than one would expect on first sight. Even arbitrarily close to the ergodic limit of the cardioid there exist tiny elliptic islands. This can be seen by considering the family of orbits A n BB that limits onto the whispering gallery orbit. This family illustrates the nonuniform hyperbolicity of the billiard map, its Liapunov exponent tends to 0 as n tends to infinity, see [6]. In Fig. 9 a plot of the parameter ε n where A n BB is born in a saddle node bifurcation is shown. The data is well described by the fitted curve ε = 1 − an −3/2 where a = 1.65. This shows that arbitrarily close to the cardioid there are stable elliptic orbits created in saddle node bifurcations. In Fig. 5 this family of orbits appears approximately as an upper bound for ∆ǫ, i.e. they tend to have the largest stability intervals for given ǫ.
Whether there are parameter values for which the corresponding limaçon billiard is ergodic cannot be decided by a numerical experiment like the one we performed. The results indicate, however, that the set of parameter values for which the map is ergodic (if larger than one point) has a very complicated structure.
To our knowledge this is the first time that the scaling of stability intervals with the period has been extensively studied. It is well known that in hyperbolic systems there is a scaling of the instability of periodic orbits with the period. This expresses the fact the the phase space volume is a constant so that growth of instability has to be exactly compensated by growth of the number of orbits [38]. Our results indicate that there is a similar relation for stability intervals in parameter space for families limiting to hyperbolic system.