Estimating Lyapunov exponents in billiards

Dynamical billiards are paradigmatic examples of chaotic Hamiltonian dynamical systems with widespread applications in physics. We study how well their Lyapunov exponent, characterizing the chaotic dynamics, and its dependence on external parameters can be estimated from phase space volume arguments, with emphasis on billiards with mixed regular and chaotic phase spaces. We show that the leading contribution to the Lyapunov exponent is generically inversely proportional to the chaotic phase space volume. We also present a numerical scheme and software implementation to calculate the Lyapunov exponents for billiards in external magnetic fields.

A billiard consists of a finite (or periodic) domain in which a point particle performs free flight with unit velocity.Upon collision with the boundary of the domain the particle (typically) is specularly reflected.In Fig. 1 we are showing the two example billiards we will be considering in this paper: the mushroom billiard (MB) and the periodic Sinai billiard without (PSB) and with magnetic field (MPSB).
An essential characterization of the chaotic dynamics of a billiard is of course provided by its Lyapunov exponents.For a two dimensional billiard the Lyapunov exponents are four numbers λ 1−4 that measure how "chaotic" the billiard is, in terms of the average exponential rates of expansion (and contraction) of the phase space along certain characteristic directions.Due to the Hamiltonian nature of the dynamics, the Lyapunov exponents fulfil λ 1 = −λ 4 , λ 2 = λ 3 = 0. Therefore in the remainder of the text we will be only considering the largest exponent λ ≡ λ 1 .The fundamental mathematical properties of the Lyapunov exponents in billiards, including rigorous proofs of their existence, have been studied in literature and can e.g.be found in Ref. [19][20][21] and references therein.
Quantitative studies of the Lyapunov exponent in actual physical billiards are surprisingly rare however.A computational framework for calculating λ in billiard systems was formulated by Dellago, Posch, and Hoover in FIG. 1.(a) A regular (blue) and chaotic (red) orbit in the mushroom billiard (MB), whose cap radius is a constant set to 1.The cap and the stem are separated with different background colors (orange, green).(b) chaotic orbits without (orange) and with (red) magnetic field and regular orbits (blue, purple) in the magnetic periodic Sinai billiard (MPSB).
Refs. [22,23] (to which we refer to as DPH framework in the following, and which we will extend to the dynamics in magnetic fields).In literature, Lyapunov exponents have been computed for the PSB on square [24] and hexagonal [22,24] lattices, as well as for for the stadium billiard [22,25], which is related to the mushroom billiard.Furthermore, there are results for the magnetic elliptical billiard [26] and the inverse magnetic stadium [27].
All these quantitative calculations rely on detailed numerical simulations of the complex billiard dynamics.In this paper, however, we want to follow a different approach exploring approximate expressions for the parameter dependence of the Lyapunov exponents in some paradigmatic cases, especially of billiards with mixed phase space, where regions of regular and chaotic dynamics coexist.Our work is motivated by a recent study that has shown that magnetoresistance measurements in graphene and semiconductor nanostructures directly reflect the parameter dependence of the chaotic phase space volume [9].This is due to the fact, that characteristic transport times in the chaotic sea are fundamentally linked to the respective volume fraction of the chaotic phase space in the corresponding billiards.Namely, for the magnetic periodic Sinai billiard (MPSB) in Fig. 1b it was analytically shown that the mean collision time κ(B) between successive collisions with the discs (of radius r) as a function of an applied external magnetic field B is equal to the varying chaotic phase space portion g c (B) times the value of κ(0) at zero magnetic field, i.e.
(The value of κ(0) is obtained from Eq. ( 3) below.) The Lyapunov exponents in billiards are also linked to the mean collision times as the following back-of-theenvelope calculation motivates.Let us study the perturbation growth, i.e. the exponential growth of the phase space distance |δΓ(t)| of two initially infinitesimally close by trajectories.The origin of the exponential perturbation growth and thus of chaos in billiards are collisions with curved boundaries [28].Assuming an average perturbation growth increase of C between collisions with curved boundaries and an average time ∆t = κ between such collisions, one would expect a perturbation growth of |δΓ(t)| ≈ C t/κ |δΓ(0)|.This means for the Lyapunov exponent (see definition in Eq. ( 4) below) we expect In general billiards have non-curved boundaries as well as curved ones, thus κ = τ .Here τ is the mean collision time with the boundary of the billiard that is known analytically for any billiard and is given by where |Q| is the area of the billiard and |∂Q| the total length of its boundary [19].The aim of this paper is to explore how far phase space volume considerations like Eq. (1) allow us to estimate the parameter dependence of the Lyapunov exponent in billiard systems.To this end we will analyze the contributions to the approximate expression (2) and compare it with detailed numerical simulations.In sec.II we provide the basic framework we will use for computing λ, as well as apply the aforementioned back-of-the-envelope calculation to realistic perturbation growth.Following in sec.III we present our results and then conclude by pointing out the generality of our approach.

II. LYAPUNOV EXPONENTS IN BILLIARDS
In this section we first give a brief overview of the DPH framework [22,23] for numerically computing λ, reciting the equations that will be relevant for our study.We will then extend the framework to motion in an external magnetic field.(In the following we will assume that the point particle in the billiard has unit mass and momentum and velocity are the same.) The (maximum) Lyapunov exponent is defined based on the evolution of the four-dimensional perturbation vector δΓ = (δq, δp) T as with δΓ evolving according to the evolution equations in tangent space, δ Γ = J(Γ(t)) • δΓ where J is the Jacobian matrix of the equations of motion.For all Γ(0) inside an ergodic component of phase space the value of λ does not depend on the initial condition.

A. Without magnetic field
The time evolution in tangent space for a particle moving in a straight line is At discrete time points Eq. ( 5) is interrupted by collisions with the boundary and the perturbation vector changes discontinuously.The perturbation vector just after the collision (we'll use to label quantities right after the collision) is derived from the one just before the collision as [23] for a collision with a boundary segment of curvature γ r .
The two types of boundaries we will encounter in this work are straight walls and the boundaries of circular discs.For a straight wall section γ r = 0 and for a disc of radius r we have γ r = ± 1 r , with − for collisions happening from the inside of the disc (as in the MB) and + otherwise (as in the PSB).Here n denotes the normal vector of the boundary segment at the collision point q, and φ is the angle of incidence.The vectors e, e are unit vectors orthogonal to the incident and reflected momenta p and p respectively (for more details see [23]).
Notice that a collision with a straight wall does not change the norm of δΓ because both the coordinate and the velocity components are reflected specularly.

B. With magnetic field
We now extend the DPH formalism for a particle experiencing a magnetic field perpendicular to the billiard plane.In this section we present only the final expressions.The full calculations are presented in appendix A.
The free evolution of the perturbation vector δΓ(t) in the presence of a perpendicular magnetic field is given by with the cyclotron frequency ω = 2B and the cyclotron radius ρ = 1/ω.In the billiard the particle always moves with unit velocity by convention.The expressions that give the discontinuous change of the perturbation vector at a collision with a wall or disc are where γ r again is the curvature of the wall segment (0 for a straight wall, ± 1 r for a disc).

C. The "toy model"
Before finding an approximate expression for the value of λ in our model systems, it is worthwhile to get an impression of how the norm of the perturbation vector evolves with time.In Fig. 2 we show typical plots for the three different billiards.We computed the perturbation growth using the DPH framework, sampling the perturbation vector immediately before and after every collision to resolve the instantaneous jumps.As the DPH evolution is linear, in the actual numerical simulations we renormalized the perturbation vector after sampling to prevent numerical error due to the rapid perturbation growth.
Let us first examine Fig. 2(a, b).We see that the norm of the perturbation vector changes in two ways.Let the j-th collision with a disc happen at time t j = j−1 i=0 ∆t i = t j−1 + ∆t j−1 .There a discontinuous change of the perturbation vector norm happens, so that |δΓ j | = a j |δΓ j | (in general a j is a function of δΓ j ).The collision event is followed by a time-interval ∆t j , in which the perturbation norm changes continuously because there are no collisions with curved boundaries.Just before the next collision with a disc the perturbation norm takes the value |δΓ j+1 |.In the following we will refer to this segment of the growth curve as a unit cell, starting with one collision event with a disc and ending just before the next one.
A crucial simplification we do in deriving an approximate expression for the Lyapunov exponent will be to express the perturbation growth in the time-interval ∆t j as a function z(∆t j ) of the interval length.The actual precise value |δΓ j+1 /δΓ j | is not a simple scalar function of the time interval since for example in the case of the PSB we can obtain from Eq. ( 5) (due to the linearity of the equations of motion of the tangent space we can assume a norm of |δΓ j | = 1 in Eq. ( 9)).Eq. ( 9) depends on the initial orientations of both the momentum and position deviation vectors and thus is not a function of just ∆t.We will show, however, by analysing numerical data, that a reasonable approximation can be obtained by assuming that such a function z(∆t) exists.For each unit cell we thus write We then recursively apply Eq. ( 10) to get For the billiards considered here we have found that furthermore log(z(∆t)) = log(z( ∆t )) = log(z(κ)) is a good approximation that we will use.
In the remainder of the text we will call Eq. ( 12) the "toy model".It is the more detailed version of the back-of-the-envelope calculation given in the introduction, i.e.C = a × z(∆t) .In the following sections we will apply this toy model to specific billiards and see how well it approximates the Lyapunov exponent and its parameter dependence.

D. Software
All numerical computations presented in this paper were performed with an open source software to simulate billiards, DynamicalBilliards.jl[29].The software also has implemented the DPH framework for computing λ, as well as its extension for magnetic fields.

A. Periodic Sinai billiard
We start our analysis with the PSB because it is the simplest case and we are able to give a fully analytical expression for the Lyapunov exponent in the simplified toy model.We note that in the absence of a magnetic field the PSB is ergodic and its phase space is not mixed.Nevertheless, it will serve as a pedagogic example of how the toy model approximates the Lyapunov exponent.
For the PSB τ = κ and the value of τ is known from Eq. ( 3) An approximation for z(∆t) is easily found as well from Eq. ( 9), namely which uses the assumption that after a collision with a disc the momentum contribution to the perturbation vector is much larger than the position contribution, i.e. |δp | |δq |.Numerical calculations show that this approximation is valid for small enough radii (see Fig. 3).The instantaneous change factor log(a) is rather large for all but very large disc radii.Also as seen in panel (c) and its inset (d), Eq. ( 15) almost perfectly approximates the perturbation norm increase during the free flight part.
We still need an approximate expression for a, the instantaneous change factor, which we can derive from Eq. ( 6).Since the norm of the position deviation does not change at the collision, we focus on the momentum deviation δp = δp (r) − 2 r δq•e cos(φ) e with δp (r) = δp−2 (δp • n) n (if not explicitly written otherwise all quantities in this paragraph have a collision-time index j, which we suppress to make the symbols simpler).We define A = We now need to average log(a) .We start our approximation by replacing the inner product (δq • e) by an averaged quantity b(r) (we show below how b depends on r).It is expected that perturbations will grow and orient themselves perpendicular to the particle's direction of motion.Since e is a unit vector perpendicular to the particle's momentum and thus parallel to δq, this means that Which portion of |δΓ| is contained in |δq| is answered based on how perturbations evolve during the free flight part.Starting from the j-th collision the perturbations evolve for time ∆t j .Using Eq. ( 9) (and assuming that the cross-terms δp • δq will drop out in the averaging) we obtain To analytically resolve eq. ( 17) we use the same assumption as in eq. ( 15), |δp | |δq |.We then average, replacing ∆t by κ, which leads to We discard the term 2A(e • δp) in Eq. ( 16) again assuming that the inner product averages to 0. Now the only variable left to average over is φ, the angle with respect to the normal vector.This is distributed in [−π/2, π/2] with probability distribution of cos(φ)/2.Therefore with csch −1 the inverse hyperbolic cosecant.
We then put eqs. ( 14), ( 15), ( 18) and ( 19) into the toy model of Eq. ( 12) and obtain an analytic approximation for the Lyapunov exponent The result is shown in Fig. 3(a), compared with the numerical value of λ using the DPH framework as well as with the result of computing the term log(a) in the toy model numerically from the evolution of the perturbation vector norm.All three curves are in excellent agreement for small and intermediate r, only for large r does Eq.( 20) slightly deviate from the numerical values because the approximation for b(r) in Eq. ( 18) and thus for log(a) becomes less accurate.

B. Magnetic periodic Sinai billiard
We now want to apply the same process to the MPSB, which, however, has a mixed phase space: there exist collision-less orbits like those seen in Fig. 1(b) that constitute the regular part of phase space.We are of course only considering the Lyapunov exponent of the chaotic part of the phase space, which means that we initialise particles only in the chaotic phase space region.The mean collision time between discs (κ MPSB below) is also only defined for the chaotic phase space part (as the regular trajectories do not collide with the discs).
The free flight evolution in the MPSB is fundamentally different from the PSB.Not only are the functional forms different but in addition due to magnetic focusing it is possible (and in fact quite frequent) for the perturbation norm to decrease during the evolution, as can be seen in Fig. 4(e,f).In addition, as visible in Fig. 2(d), it is also possible for the norm to decrease during the instantaneous change as well.
This more complex behaviour is of course hidden in the more complicated formulas of our extension to the DPH framework for magnetic fields.For example, explicitly writing out Eq. (7) gives (where again we assumed |δΓ(0)| = 1).The consequences of Eq. ( 21) can be seen in Fig. 4(d, e).Using a uni-variate scalar function z(∆t) to approximate these distributions appears to be a bold move, but in the end it will turn out to give a good approximation.To obtain z(∆t) we simplify Eq. ( 21) to which is also plotted in Fig. 4(d, e).
In the next step we compute log(a) numerically and use its value in the toy model along with z MPSB (κ(B)).We remind that the value of κ, the mean collision time between discs in MPSB, is not known analytically but it is connected with the chaotic phase space portion through Eq. ( 1).The results of the toy model are presented in Fig. 4.
Besides the fact that our toy model approximates λ very accurately, Fig. 4 nicely shows the impact of phase space restriction on λ.In our toy model the value of λ is composed of five contributions, the first being the denominator κ.The value of log(a) itself has two contributions, one again stemming from κ (as shown in the previous section) and the other from B. The function z MPSB also has two contributions, one from B and one from κ. Therefore three out of five contributions to λ are inherently linked to the restriction of the chaotic phase space by regular orbits.

C. Mushroom billiard
Because the volume fractions of the regular and the chaotic phase space regions are not known analytically in the MPSB we have turned to a billiard that also has a mixed phase space but allows to calculate this fractions, and, as we will show, the relevant average time scales analytically: the mushroom billiard (MB).The regular orbits in the MB are orbits forever staying in the cap, evolving exactly like they would as if they were in a circular billiard [30].The tangential circle to these orbits has a radius ≥ w/2 as shown in Fig. 1(a).The rest of the orbits, which do not satisfy this criterion, eventually enter the stem and are chaotic.The tangential circle argument was used in [5] to obtain an analytic expression for the regular phase space volume V REG of the MB as a function of the billiard parameters where V TOT is the parameter dependence of the total phase space volume (all lengths are scaled to the cap radius r which therefore is fixed to r = 1).We then obtain the chaotic phase space portion as g c = 1 − V REG /V TOT , which we plot in Fig. 5(c,d).Interestingly, g c does not appear to vanish for small h, although it is obvious that there are no chaotic orbits for h = 0.This discontinuity is due to the fact that the volume of chaotic phase space in the cap is independent of stem height for nonzero h, but drops to zero for h = 0.In Fig. 6(c,d) we present a scatter plot of various possible increases of the perturbation norm during the unit cells.We found that there are clearly distinct contributions to the increase, each seemingly approximated as linear function of ∆t.By analysing the dynamics in more detail, it turns out that the different contributions of Fig. 6(c,d) come from the trapping of the chaotic orbits in the regular phase space.In coordinate space this means that the particles get trapped in the cap and mimic the motion of the regular phase space there until eventually escaping after some time.This effect is often called intermittency, and is known to occur in mushroom billiards [31,32].
We therefore have to separate the unit cell into two different "episodes": the chaotic episode c and the laminar episode , where the particle is trapped in the cap.We show these episodes in Fig. 5(a).Numeric calculations shown in Fig. 6(e, f) show that each episode has a different average time, τ c , τ respectively.
During the chaotic episodes the picture is very similar to the PSB.A collision with the cap head gives an instantaneous increase to |δΓ|, followed by an approximately linear increase until the next collision with the cap head.Here the linear increase approximation is valid because for the chaotic episodes ∆t 2h + 2 − w.After colliding with the cap head the particle may return to the stem immediately which initialises another chaotic episode.Occasionally, after ending a chaotic episode, the particle will get trapped in the cap (see Fig. 5(a) orange), starting a laminar episode.Even though there are successive collisions with the cap head in this episode, the perturbations do not increase exponentially.The successive instantaneous increases are very quickly becoming insignificant (see Fig. 2(e,f)) due to the fact that cap collisions have an initially focusing effect which only becomes de-focusing if the consecutive free motion is long enough, which is not the case in the laminar episodes.Therefore the overall perturbation growth inside the cap trapping episodes is linear.
Let n c and n be the counts of chaotic and laminar episodes up to time T .Notice that n is strictly less than n c since a chaotic episode always follows a laminar episode, but the inverse is only occasionally true.In the limit T → ∞ we define f = n /(n + n c ) to be the frequency of the laminar episodes.We then write the function z as with o i the offset and s i the slope of the linear approximation (we obtain these values with least squares fit to Fig. 6(c, d)).For the chaotic episodes s, o are constant versus h, w while for the laminar episodes o depends strongly on w.Also, for the chaotic episodes o has negative value (of around -0.8) which is expected due to the focusing effect.We once again compute log(z(∆t)) simply by replacing ∆t by its average values τ c , τ in Eq. ( 25).
The instantaneous change factor a is the same between the laminar and chaotic episodes so we do not need to separate it.Notice that for the laminar episode we only consider the first jump as the instantaneous increase.Subsequent jumps that decrease rapidly are encoded in the linear growth approximation.After computing log(a) numerically, we still need a value for κ MB , the unit cell average time, to apply our toy model λ = ( log(a) + log(z MB (τ c , τ ))/κ MB .Numerically we can estimate However we can estimate κ MB analytically as well, using Kac's lemma [33][34][35].The key to this is understanding that the mean cell time is equivalent with the mean return time to the stem bottom, since all phases in the end of the day have to go there, since all phases are part of the chaotic phase space.We present the full proof in appendix B. The final expression is given by with g c here being the portion of chaotic phase space of the MB.We compare the analytic formula with the numeric result in Fig. 5(b) and find the expected perfect agreement, since Eq. ( 27) is exact.In appendix C we also present an analytic approximation for τ c .Since we know κ MB and τ c analytically, we also know the product f τ (but we don't have an expression for f or τ individually).
We can now use our toy model to compare with λ, which we do in Fig. 6(a, b).Again we find good agreement between toy model and the numerical simulation using the DPH framework.The model mildly diverges for very small w, probably because the mean laminar time τ diverges as seen in Fig. 6(e).
As was the case in the MPSB, the average unit cell time κ is inversely proportional to λ and directly proportional to the chaotic phase space portion g c .This shows that phase space restrictions have an immediate impact in the value of the Lyapunov exponent even for billiards with intermittent dynamics.Notice that in the MB both g c and λ increase as w increases.This is for two reasons; firstly κ depends on the chaotic phase space volume and not just the relative proportion.The total volume of the MB changes with w (which is not the case in the MPSB: the total volume is independent of B).Secondly, it happens here that besides g c there are also other factors that depend on w in Eq. ( 27), which was not the case in Eq. ( 1) for the MPSB.

IV. CONCLUSIONS
To summarise, we have examined the value of the Lyapunov exponent λ in chaotic billiards.We were able to create a conceptually simple model that approximates λ very well.The model is based on how perturbations evolve in billiards on average and helps to understand how each part of the dynamics contributes to the perturbation increase.The simple model is written as λ = 1 κ ( log (a) + log (z(∆t)) ) (28) where a the instantaneous change of |δΓ| at a collision with a curved boundary and z(t) the continuous change of |δΓ| in between collisions with curved boundaries.κ is the average unit cell time equal to the mean collision time between curved boundaries.We believe Eq. ( 28) is a valid approximation for λ for generic chaotic billiards.Furthermore, for the billiards considered here we found that log (z(∆t)) ≈ log(z(κ)) is a good approximation as well.We used Eq. ( 28) to find an analytic expression for λ in the PSB.We have also shown that this simplified expression can closely replicate the numerically exact values of λ in the MPSB.We could follow the same approach for the MB, even though the process is complicated in this case by intermittent dynamics.In both billiards with mixed phase space we were able to connect the chaotic phase space portion with λ through κ.Due to the generality of Kac's lemma, we believe this connection to be universal.Thus in conclusion we see that λ is in general roughly inversely proportional to the chaotic phase space portion, allowing an estimation of the Lyapunov exponent and its parameter dependence in mixed phase space systems.We have seen, however, that in a more detailed analysis κ and thus the phase space fractions additionally enter the collision term log (a) and that magnetic focusing and intermittency can lead to additional corrections to the Lyapunov exponent.

FIG. 2 .
FIG.2.Typical time evolution of the logarithm of the norm of the perturbation vector |δΓ(t)| for the periodic Sinai billiard without (a) and with (c) magnetic field (B = 1) and for the mushroom billiard (e).Zoom-ins are below each panel.For (a-d) red markers mean collision with the disc while blue markers mean "collision" with the periodic walls (not a true collision, but a way of recording the value of |δΓ|).For (e-f) red is collision with cap head, orange with cap walls, blue with stem sides and green with stem bottom.The coloured background stripes denote the unit cells (random colors are used) with hatched orange color used for the laminar phases.

FIG. 3 .
FIG. 3. (a) Lyapunov exponent of the periodic Sinai billiard for different radii, compared with the toy model.The dashed curve obtains log(a) by numeric average while the red curve uses Eq. (19).The blue curve is using the DPH framework.(b) Average value of log(a) used in panel (a).(c) Perturbation norm increase during the free flight part.In the zoom in (d) we also plot the curve √ 1 + ∆t 2 as a dashed line.

FIG. 4 .
FIG. 4. (a) Lyapunov exponent of the magnetic periodic Sinai billiard (MPSB) for different radii versus the magnetic field, compared with the toy model.The solid curves are the numeric result λ, the dashed curves are the toy model using log(a) .(b) Numeric average of log(a) versus the magnetic field.(c) Chaotic phase space portion (gc = κ/κ(0)) of the MPSB ((a-c) share legend).(d, e) Perturbation norm change during the free flight in the MPSB.Plotted with dashed lines are Eq.(22) (data for r = 0.2).

FIG. 5 .
FIG. 5. (a) A laminar episode (orange, dashed) and two chaotic episodes (blue, purple) in the mushroom billiard.Starts and ends of each episode are denoted with closed and open circles and the blue episode starts directly after the orange.The cap head is plotted in dark color to differentiate.(b) Mean return time to the stem bottom, which is equivalent with the average unit cell time, eqs.(26), (27).(c,d) Portion of chaotic phase space for the mushroom billiard.

FIG. 6 .
FIG. 6. (a, b) Lyapunov exponents in the mushroom billiard (MB) versus the width w or height h of the stem.Solid lines are numeric results using the DPH framework, dashed lines are using the toy model.(c, d) Perturbation norm increase during the chaotic c and laminar episodes.(e, f) Parameters of the toy model versus w or h (for constant h = 1 and w = 1 respectively); legend is shared.