Estimating Lyapunov exponents in billiards

TL;DR

This paper investigates how to estimate Lyapunov exponents in billiards using phase space volume arguments, revealing an inverse proportionality in mixed phase spaces and extending existing methods to include magnetic fields.

Contribution

It introduces a phase space volume-based estimation method for Lyapunov exponents in billiards and extends existing formalism to account for magnetic fields.

Findings

Lyapunov exponent inversely proportional to chaotic phase space volume

Extended formalism to include magnetic fields

Provided software implementation for calculations

Abstract

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 in the very diverse billiards considered here the leading contribution to the Lyapunov exponent is inversely proportional to the chaotic phase space volume, and subsequently discuss the generality of this relationship. We also extend the well established formalism by Dellago, Posch, and Hoover to calculate the Lyapunov exponents of billiards to include external magnetic fields and provide a software implementation of it.

Figures

Equations

$$\kappa (B)=g_c(B)\times \kappa (0)=g_c(B)\times \frac{1-\pi r^2}{2r}.$$

$$\lambda \approx \frac{\log (C)}{\kappa }.$$

$$\tau =\frac{\pi |Q|}{|\partial Q|}$$

$${\lambda }_{\mathbf{\Gamma }(0),\delta \mathbf{\Gamma }(0)}=\underset{t\to \infty }{\lim }\frac{1}{t}\log \frac{|\delta \mathbf{\Gamma }(t)|}{|\delta \mathbf{\Gamma }(0)|}$$

$$\delta\mathbf{\Gamma}(t)=\left(\begin{array}[]{c|c}\mathbb{I}_{2\times 2}&t\cdot\mathbb{I}_{2\times 2}\\ \hline\cr 0_{2\times 2}&\mathbb{I}_{2\times 2}\end{array}\right)\;\delta\mathbf{\Gamma}(0).$$

$$\delta\mathbf{\Gamma}^{\prime}=\left(\begin{array}[]{l}\delta\mathbf{q}-2\left(\delta\mathbf{q}\cdot\mathbf{n}\right)\mathbf{n}\\ \delta\mathbf{p}-2\left(\delta\mathbf{p}\cdot\mathbf{n}\right)\mathbf{n}-2\gamma_{r}\;\tfrac{\left(\delta\mathbf{q}\cdot\mathbf{e}\right)}{\cos\phi}\;\mathbf{e}^{\prime}\\ \end{array}\right)$$

$$\delta \mathbf{\Gamma }(t)=\mathbb{B}\cdot \delta \mathbf{\Gamma }(t_0),$$

$$\displaystyle\mathbb{B}=\left(\begin{array}[]{c|cc}\mathbb{I}_{2\times 2}&\begin{array}[]{c}\rho\sin(\omega t)\\ -\rho\left(\cos(\omega t)-1\right)\end{array}&\begin{array}[]{c}\rho\left(\cos(\omega t)-1\right)\\ \rho\sin(\omega t)\end{array}\\ \hline\cr 0_{2\times 2}&\begin{array}[]{c}\cos(\omega t)\\ \sin(\omega t)\end{array}&\begin{array}[]{c}-\sin(\omega t)\\ \cos(\omega t)\end{array}\end{array}\right)$$

$$\delta {\mathbf{\Gamma }}^{\prime }=$$

$$\omega \frac{\left(\delta \mathbf{q}\cdot \mathbf{n}\right)}{\left(\mathbf{p}\cdot \mathbf{n}\right)}\text{pmqty}undef$$

$$\mathbb{S}=$$

$$|\delta {\mathbf{\Gamma }}_{j+1}(\Delta t)|=\sqrt{(\delta {\mathbf{q}}_j^{\prime }+\delta {\mathbf{p}}_j^{\prime }\Delta t{)}^2+(\delta {\mathbf{p}}_j^{\prime }{)}^2}$$

$$|\delta {\mathbf{\Gamma }}_{j+1}|=z(\Delta t_j)\times a_j\times |\delta {\mathbf{\Gamma }}_j|=z(\Delta t_j)\times |\delta {\mathbf{\Gamma }}_j^{\prime }|.$$

$$|\delta {\mathbf{\Gamma }}_n|$$

$$\log \left(|\delta {\mathbf{\Gamma }}_n|\right)$$

$$\lambda$$

$$\kappa$$

$${\kappa }_{\text{PSB}}=\frac{1-\pi r^2}{2r}.$$

$$z_{\text{PSB}}(t)\approx \sqrt{1+t^2}$$

$$|\delta {\mathbf{\Gamma }}^{\prime }|/|\delta \mathbf{\Gamma }|=a=\sqrt{1+A^2-2A({\mathbf{e}}^{\prime }\cdot \delta \mathbf{p})}$$

$$\frac{|\delta {\mathbf{q}}_{j+1}|}{|\delta {\mathbf{\Gamma }}_{j+1}|}=\frac{\sqrt{\frac{|\delta {\mathbf{q}}_j^{\prime }{|}^2}{|\delta {\mathbf{p}}_j^{\prime }{|}^2}+\Delta t_j^2}}{\sqrt{\frac{|\delta {\mathbf{q}}_j^{\prime }{|}^2}{|\delta {\mathbf{p}}_j^{\prime }{|}^2}+1+\Delta t_j^2}}.$$

$$b(r)=\sqrt{\frac{{\kappa }_{\text{PSB}}^2}{1+{\kappa }_{\text{PSB}}^2}}.$$

$$⟨\log (a)⟩$$

$$=\frac{{\text{csch}}^{-1}\left(\frac{2b}{r}\right)\sqrt{4b^2+r^2}}{r}+\log \left(\frac{b}{r}\right)$$

$${\lambda }_{\text{PSB}}(r)=$$

$$\log \left(\frac{b(r)}{r}\right)+\log \left(\sqrt{1+{\left(\frac{1-\pi r^2}{2r}\right)}^2}\right)).$$

$$|\delta \mathbf{\Gamma }(t)|=$$

$${\left(\delta q_x+\delta p_y\frac{\cos (\omega t)-1}{\omega }+\delta p_x\frac{\sin (\omega t)}{\omega }\right)}^2+$$

$${{\left(\delta q_y-\delta p_x\frac{\cos (\omega t)-1}{\omega }+\delta p_y\frac{\sin (\omega t)}{\omega }\right)}^2]}^{1/2}$$

$$z_{\text{MPSB}}(t)=\sqrt{1+{\left(\frac{1-\cos (\omega t)}{\omega }\right)}^2+{\left(\frac{\sin (\omega t)}{\omega }\right)}^2}.$$

$$V_{\text{REG}}$$

$$V_{\text{TOT}}$$

$$V_{\text{CH}}$$

$$z_{\text{MB}}(\Delta t)=(o_c+s_c\Delta t_c)+f_{\ell }(o_{\ell }+s_{\ell }\Delta t_{\ell })$$

$${\kappa }_{\text{MB}}={\tau }_c+f_{\ell }{\tau }_{\ell }.$$

$${\kappa }_{\text{MB}}(h,w)=\frac{V_{\text{CH}}(h,w)}{2w}$$

Full text

**From the foundations of statistical physics to transport properties of electronic devices, in many areas of physics billiard models are an important tool for understanding complex dynamics. In a billiard model a point particle is moving freely (and frictionless) on a flat (or constantly curved) surface until it hits the boundary of the billiard where it is specularly reflected. Chaotic dynamics in the billiard are characterized by a positive Lyapunov exponent, measuring how initially close trajectories separate exponentially fast. Obtaining its value so far usually requires detailed numerical simulations of the chaotic dynamics. In our paper we assess how well the Lyapunov exponent can be estimated from quite general considerations. Especially we study how parameter changes that vary the phase space structure of the billiard get reflected in the Lyapunov exponent. For example the application of an external magnetic field can force some trajectories in the billiard on closed cyclotron orbits. We show how the mere existence of such orbits varies the Lyapunov exponent of the chaotic dynamics through the phase space volume they occupy . The knowledge of this connection will be helpful to understand physical mechanisms in many systems like the magneto-transport in graphene nanostructures. **

I Introduction

Dynamical billiards are a well-studied class of dynamical systems, having applications in many different fields of physics. Besides playing a prominent role in ergodic theory Sinai (1970); Bunimovich (1974, 1979), billiards are important example systems for understanding quantum chaos Stöckmann and Stein (1990); Barnett and Betcke (2007), with practical applications e.g. in modelling optical microresonators for lasers Stone (2010); Gmachl et al. (1998) and room acoustics Koyanagi et al. (2008). Billiard-models have also been particularly successful in helping to understand transport properties of electronic nanostructures such as quantum dots and antidot super-lattices Datseris et al. (2019); Fleischmann et al. (1992); Weiss et al. (1991); Chernov et al. (1993); Yagi et al. (2015); Maier et al. (2017); Chen et al. (2016); Bird (2003); Beenakker (2005); Ponomarenko et al. (2008).

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) Bunimovich (2001) and the periodic Sinai billiard Sinai (1970) 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. (In this article we will study the Lyapunov exponents of the billiard flow in the physical, continuous time, in contrast to those of the billiard map in a discrete time that counts the number of collisions with the boundary). For a two dimensional billiard the Lyapunov exponents are four numbers $\lambda_{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 $\lambda_{1}=-\lambda_{4},\lambda_{2}=\lambda_{3}=0$. Therefore in the remainder of the text we will be only considering the largest exponent $\lambda\equiv\lambda_{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. Chernov (1997); Chernov and Young (2000); Chernov and Markarian (2006) and references therein.

Quantitative studies of the Lyapunov exponent in actual physical billiards are surprisingly rare however. A computational framework for calculating $\lambda$ in billiard systems was formulated by Dellago, Posch, and Hoover in Refs. Dellago and Posch (1995); Dellago et al. (1996) (to which we refer to as DPH framework in the following, and which we will extend to the dynamics in magnetic fields). Alternative approaches are presented e.g in Refs. Chernov (1991); Garrido (1997). In the literature, Lyapunov exponents have been computed for the PSB on square Gaspard and Baras (1995) and hexagonal Gaspard and Baras (1995); Dellago and Posch (1995) lattices, as well as for for the stadium billiard Dellago and Posch (1995); Posch and Hirschl (2000), which is related to the mushroom billiard. Furthermore, there are results for the magnetic elliptical billiard Meplan et al. (1993) and the inverse magnetic stadium Vörös et al. (2003).

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 Datseris et al. (2019). This is due to the fact that characteristic transport times in the chaotic sea are fundamentally linked to the respective volume 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 $\kappa(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 $\kappa$ at zero magnetic field Datseris et al. (2019), i.e.

[TABLE]

(The value of $\kappa(0)$ is obtained from Eq. (3) below.) For the convenience of the reader we replicate the proof of eq. (1) from Datseris et al. (2019), which uses Kac’s lemma Kac (1947); MacKay (1994); Meiss (1997), in Appendix D.

The Lyapunov exponent in billiards is also linked to mean collision times as the following back-of-the-envelope calculation motivates. Let us study the perturbation growth, i.e. the exponential growth of the phase space distance $|\delta\mathbf{\Gamma}(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 Bunimovich (2018). Assuming an average perturbation growth increase of $C$ between collisions with curved boundaries and an average time $\kappa$ between such collisions, one would expect a perturbation growth of $|\delta\mathbf{\Gamma}(t)|\approx C^{t/\kappa}|\delta\mathbf{\Gamma}(0)|$. This means for the Lyapunov exponent of the ergodic component of phase space (see definition in Eq. (4) below) we expect

[TABLE]

In general, billiards have non-curved boundaries as well as curved ones. The mean collision time $\tau$ between two consecutive collisions with any parts of the billiard boundary is known analytically for any billiard and is given by

[TABLE]

where $|Q|$ is the area of the billiard and $|\partial Q|$ the total length of its boundary Chernov (1997). Since eq. (3) is averaged over the entire boundary of the billiard it includes contributions from both chaotic and regular components (if any). Also notice that a formula similar to eq. (3) exists for billiards of any dimension, see Chernov (1997).

The starting point our work is the observation that in billiards the mean collision time between curved boundaries $\kappa$ (a more precise definition will be given in section II.3) is fundamentally linked to the chaotic phase space volume $V_{\mathrm{CH}}$ by Kac’s lemma Kac (1947); MacKay (1994); Meiss (1997). Therefore also the Lyapunov exponent is linked to the chaotic phase space volume, and the aim of this paper is to explore how far 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 $\lambda$, as well as apply the aforementioned back-of-the-envelope calculation to realistic perturbation growth. Following in sec. III we present our results about the periodic Sinai billiard and the mushroom billiard. We conclude by discussing the generality of our approach, while presenting one additional billiard with mixed phase space, the inverse stadium billiard, which has been studied by Vörös et al. Vörös et al. (2003).

II Lyapunov exponents in billiards

In this section we first give a brief overview of the DPH framework Dellago and Posch (1995); Dellago et al. (1996) for numerically computing $\lambda$, 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 $\delta\mathbf{\Gamma}=(\delta\mathbf{q},\delta\mathbf{p})^{T}$ as

[TABLE]

with $\delta\mathbf{\Gamma}$ evolving according to the evolution equations in tangent space, $\delta\dot{\mathbf{\Gamma}}=J(\mathbf{\Gamma(t))}\cdot\delta\mathbf{\Gamma}$ where $J$ is the Jacobian matrix of the equations of motion. For all $\mathbf{\Gamma}(0)$ inside an ergodic component of phase space the value of $\lambda$ does not depend on the initial condition.

II.1 Without magnetic field

The time evolution in tangent space for a particle moving in a straight line is

[TABLE]

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 Dellago et al. (1996)

[TABLE]

for a collision with a boundary segment of curvature $\gamma_{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 $\gamma_{r}=0$ and for a disc of radius $r$ we have $\gamma_{r}=\pm\frac{1}{r}$, with $-$ for collisions happening from the inside of the disc (as in the MB) and $+$ otherwise (as in the PSB). Here $\mathbf{n}$ denotes the normal vector of the boundary segment at the collision point $\mathbf{q}$, and $\phi$ is the angle of incidence (measured with respect to the normal to the surface). The vectors $\mathbf{e},\mathbf{e}^{\prime}$ are unit vectors orthogonal to the incident and reflected momenta $\mathbf{p}$ and $\mathbf{p}^{\prime}$ respectively (for more details see Dellago et al. (1996)).

Notice that a collision with a straight wall does not change the norm of $\delta\mathbf{\Gamma}$ because both the coordinate and the velocity components are reflected specularly.

II.2 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 magnetic field is uniform, with value $B$ (positive means counter-clockwise rotation). The free evolution of the perturbation vector $\delta\mathbf{\Gamma}(t)$ in the presence of a perpendicular magnetic field is given by

[TABLE]

with the cyclotron frequency $\omega=2B$ and the cyclotron radius $\rho=1/\omega$. As already mentioned in the introduction, 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

[TABLE]

where $\gamma_{r}$ again is the curvature of the wall segment ([math] for a straight wall, $\pm\tfrac{1}{r}$ for a disc).

II.3 The “toy model”

Before finding an approximate expression for the value of $\lambda$ 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}=\sum_{i=0}^{j-1}\Delta t_{i}=t_{j-1}+\Delta t_{j-1}$. There a discontinuous change of the perturbation vector norm happens, so that $|\delta\mathbf{\Gamma}^{\prime}_{j}|=a_{j}|\delta\mathbf{\Gamma}_{j}|$ (in general $a_{j}$ is a function of $\delta\mathbf{\Gamma}_{j}$). The collision event is followed by a time-interval $\Delta 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 $|\delta\mathbf{\Gamma}_{j+1}|$. In the following we will refer to these repeated segments of the growth curve as elementary growth segments, starting with one collision event with a disc and ending just before the next one. In general it is the segment of the perturbation growth curve between successive dispersing or defocussing collisions which are the origin of chaos in billiards Bunimovich (2018). An elementary growth segment thus reflects the perturbation growth during an “effective free path” as defined by Bunimovich Bunimovich (2018).

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 $\Delta t_{j}$ as a function $z(\Delta t_{j})$ of the interval length. The actual precise value $|\delta\mathbf{\Gamma}_{j+1}/\delta\mathbf{\Gamma}^{\prime}_{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)

[TABLE]

(due to the linearity of the equations of motion of the tangent space we can assume a norm of $|\delta\mathbf{\Gamma}_{j}^{\prime}|=1$ in Eq. (20)). Eq. (20) depends on the initial orientations of both the momentum and position deviation vectors and thus is not a function of just $\Delta t$.

We will show, however, by analysing numerical data, that a reasonable approximation can be obtained by assuming that such a function $z(\Delta t)$ exists. Notice that this assumption regards only the existence of $z$. The functional form and its complexity can be completely arbitrary (and in fact in the following we have three different versions of $z$ for the different billiards). For each elementary growth segment we thus write

[TABLE]

We then recursively apply Eq. (21) to get

[TABLE]

(using $|\delta\mathbf{\Gamma}(0)|=1$) and with $T_{n}=\sum^{n-1}_{i=0}\Delta t_{i}$ we use Eq. (4) to write $\lambda=\lim_{n\to\infty}\log\left(|\delta\mathbf{\Gamma}_{n}|\right)/T_{n}$. The quotient of the infinite sums is the same as the ratio of the average over all unit cells (denoted by $\Braket{\cdot}$), i.e.

[TABLE]

Averaging over the unit cells implicitly assumes the ergodicity. Notice also that in some billiards there could be several ergodic chaotic components that are separated from each other. In such cases the above process has to be applied to each component separately, since each component has its own exponent $\lambda$. For the billiards considered here we have found that furthermore $\Braket{\log(z(\Delta t))}=\log(z(\Braket{\Delta t}))=\log(z(\kappa))$ is a good approximation that we will use. This approximation is valid when the standard deviation of $\Delta t$ is small (compared to its mean).

In the remainder of the text we will call Eq. (23) the “toy model”. It is the more detailed version of the back-of-the-envelope calculation given in the introduction. 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.

II.4 Software

All numerical computations presented in this paper were performed with an open source software to simulate billiards, DynamicalBilliards.jl Datseris (2017). In this paper we extend the DPH framework to magnetic fields. We also implemented this extension in the software (which previously only included the non-magnetic case). All code we used for this paper, including all code to replicate the figures we show here, is publicly available on GitHub: https://github.com/Datseris/arXiv_1904.05108.

III Results

III.1 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 $\tau=\kappa$ and the value of $\tau$ is known from Eq. (3)

[TABLE]

An approximation for $z(\Delta t)$ is easily found as well from Eq. (20), namely

[TABLE]

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. $|\delta\mathbf{p}^{\prime}|\gg|\delta\mathbf{q}^{\prime}|$. Numerical calculations show that this approximation is valid for small enough radii (see Fig. 3). The instantaneous change factor $\Braket{\log(a)}$ is rather large for all but very large disc radii. Also as seen in panel (c) and its inset (d), Eq. (26) 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 $\delta\mathbf{p}^{\prime}=\delta\mathbf{p}^{(r)}-\frac{2}{r}\frac{\delta\mathbf{q}\cdot\mathbf{e}}{\cos(\phi)}\mathbf{e}^{\prime}$ with $\delta\mathbf{p}^{(r)}=\delta\mathbf{p}-2\left(\delta\mathbf{p}\cdot\mathbf{n}\right)\mathbf{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=\frac{2}{r}\frac{\delta\mathbf{q}\cdot\mathbf{e}}{\cos(\phi)}$ and carrying out the calculations leads to

[TABLE]

again using the assumption that $|\delta\mathbf{\Gamma}_{j}|=1$ (and recall that $|\delta\mathbf{p}^{(r)}|=|\delta\mathbf{p}|,|\delta\mathbf{q}^{\prime}|=|\delta\mathbf{q}|$).

We now need to average $\Braket{\log(a)}$. We start our approximation by replacing the inner product $(\delta\mathbf{q}\cdot\mathbf{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 $\mathbf{e}$ is a unit vector perpendicular to the particle’s momentum and thus parallel to $\delta\mathbf{q}$, this means that $\Braket{}{(\delta\mathbf{q}\cdot\mathbf{e})}{/|\delta\mathbf{\Gamma}|}=\Braket{}{\delta\mathbf{q}}{/|\delta\mathbf{\Gamma}|}=b(r)$.

Which portion of $|\delta\mathbf{\Gamma}|$ is contained in $|\delta\mathbf{q}|$ is answered based on how perturbations evolve during the free flight part. Starting from the $j$-th collision the perturbations evolve for time $\Delta t_{j}$. Using Eq. (20) (and assuming that the cross-terms $\delta\mathbf{p}\cdot\delta\mathbf{q}$ will drop out in the averaging) we obtain

[TABLE]

To analytically resolve eq. (28) we use the same assumption as in eq. (26), $|\delta\mathbf{p}^{\prime}|\gg|\delta\mathbf{q}^{\prime}|$. We then average, replacing $\Delta t$ by $\kappa$, which leads to

[TABLE]

We discard the term $2A(\mathbf{e^{\prime}}\cdot\delta\mathbf{p})$ in Eq. (27) again assuming that the inner product averages to 0. Now the only variable left to average over is $\phi$, the angle with respect to the normal vector. This is distributed in $[-\pi/2,\pi/2]$ with probability distribution of $\cos(\phi)/2$. Therefore

[TABLE]

with $\text{csch}^{-1}$ the inverse hyperbolic cosecant.

We then put eqs. (25), (26), (29) and (30) into the toy model of Eq. (23) and obtain an analytic approximation for the Lyapunov exponent

[TABLE]

The result is shown in Fig. 3(a), compared with the numerical value of $\lambda$ using the DPH framework as well as with the result of computing the term $\Braket{\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. (31) slightly deviate from the numerical values because the approximation for $b(r)$ in Eq. (29) and thus for $\Braket{\log(a)}$ becomes less accurate.

III.2 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 (other unstable periodic orbits of zero measure are not relevant here). 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 $\kappa_{\text{MPSB}}$ between successive collisions with discs 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

[TABLE]

(where again we assumed $|\delta\mathbf{\Gamma}(0)|=1$). The consequences of Eq. (32) can be seen in Fig. 4(d, e). Using a uni-variate scalar function $z(\Delta 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(\Delta t)$ we simplify Eq. (32) to

[TABLE]

which is also plotted in Fig. 4(d, e).

In the next step we compute $\Braket{\log(a)}$ numerically and use its value in the toy model along with $z_{\text{MPSB}}(\kappa(B))$. We remind that the value of $\kappa$, 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 $\lambda$ very accurately, Fig. 4 nicely shows the impact of phase space restriction on $\lambda$. In our toy model the value of $\lambda$ is composed of five contributions, the first being the denominator $\kappa$. The value of $\Braket{\log(a)}$ itself has two contributions, one again stemming from $\kappa$ (as shown in the previous section) and the other from $B$. The function $z_{\text{MPSB}}$ also has two contributions, one from $B$ and one from $\kappa$. Therefore three out of five contributions to $\lambda$ are inherently linked to the restriction of the chaotic phase space by regular orbits.

III.3 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 these 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 Bunimovich (2001). The tangential circle to these orbits has a radius $\geq 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 Barnett and Betcke (2007) to obtain an analytic expression for the regular phase space volume $V_{\text{REG}}$ of the MB as a function of the billiard parameters

[TABLE]

where $V_{\text{TOT}},V_{\text{CH}}$ is the total and chaotic phase space volume (all lengths are scaled to the cap radius $r$ which is fixed to $r=1$). The parameter dependence of $V_{\text{CH}}$ is illustraded in Fig. 5 c and d. Interestingly, $V_{\text{CH}}$ does not 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 $\Delta 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 Altmann et al. (2005, 2006). Intermittent behaviour in the MB can happen in the stem as well, where orbits stay trapped bouncing between the stem walls.

We therefore have to separate the elementary growth segment into two different “episodes”: the chaotic episode $c$ and the laminar episode $\ell$, where the particle is trapped in the cap. Notice that the second intermittent behaviour, trapping in the stem, does not lead to a new type of dynamics, but is just prolonged chaotic episodes (similarly to a large free flight in the PSB). We show the two episodes in Fig. 5(a). Numeric calculations shown in Fig. 6(e, f) show that each episode has a different average time, $\tau_{c},\tau_{\ell}$ respectively.

During the chaotic episodes the picture is very similar to the PSB. A collision with the cap head gives an instantaneous increase to $|\delta\mathbf{\Gamma}|$, 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 $\Delta t\gtrsim 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_{\ell}$ be the counts of chaotic and laminar episodes up to time $T$. Notice that $n_{\ell}$ 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\to\infty$ we define $f_{\ell}=n_{\ell}/(n_{\ell}+n_{c})$ to be the frequency of the laminar episodes. We then write the function $z$ as

[TABLE]

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 $\Braket{\log(z(\Delta t))}$ simply by replacing $\Delta t$ by its average values $\tau_{c},\tau_{\ell}$ in Eq. (37).

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 $\Braket{\log(a)}$ numerically, we still need a value for $\kappa_{\text{MB}}$, the elementary growth segment average time, to apply our toy model $\lambda=(\Braket{\log(a)}+\log(z_{\text{MB}}(\tau_{c},\tau_{\ell}))/\kappa_{\text{MB}}$. Numerically we can estimate

[TABLE]

However we can estimate $\kappa_{\text{MB}}$ analytically as well, using Kac’s lemma Kac (1947); MacKay (1994); Meiss (1997). The key to this is understanding that the mean elementary growth segment 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

[TABLE]

We compare the analytic formula with the numeric result in Fig. 5(b) and find the expected perfect agreement, since Eq. (39) is exact. In appendix C we also present an analytic approximation for $\tau_{c}$. Since we know $\kappa_{\text{MB}}$ and $\tau_{c}$ analytically, we also know the product $f_{\ell}\tau_{\ell}$ (but we don’t have an expression for $f_{\ell}$ or $\tau_{\ell}$ individually).

We can now use our toy model to compare with $\lambda$, 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 $\tau_{\ell}$ diverges as seen in Fig. 6(e).

As was the case in the MPSB, the average elementary growth segment time $\kappa$ is inversely proportional to $\lambda$ and directly proportional to the chaotic phase space volume $V_{\text{CH}}$. 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 $V_{\textbf{CH}}$ and $\lambda$ increase as $w$ increases. This is simply due to the dependence of $V_{\text{CH}}$ on $w$, as well as the direct dependence of $\kappa_{\text{MB}}$ on $1/w$ (this for example was not the case in Eq. (1) for the MPSB).

IV Discussion

To summarise, we have examined the value of the Lyapunov exponent $\lambda$ in chaotic billiards. We were able to create a conceptually simple model that approximates $\lambda$ 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 eq. (23), which is

[TABLE]

where $a$ the instantaneous change of $|\delta\mathbf{\Gamma}|$ at a collision with a curved boundary and $z(t)$ the continuous change of $|\delta\mathbf{\Gamma}|$ in between collisions with curved boundaries. $\kappa$ is the average elementary growth segment time equal to the mean collision time between curved boundaries. The approximations that lead to the toy model were the following. First we assumed that the chaotic phase space is ergodic and time averages can be replaced by phase space averages and that for $\Delta t,\log(a),\log(z(\Delta t)),$ their averages are finite and greater than 0. We then made the simplifying assumption that the norm of the perturbation vector increases continuously in between successive chaos-inducing collisions (i.e. in each elementary growth segment) as a univariate function of the time interval $z(\Delta t)$.

We used Eq. (23) to find an analytic expression for $\lambda$ in the periodic Sinai billiard (PSB). We have also shown that Eq. (23) can be used to analyse the Lyapunov exponent in the magnetic periodic Sinai billiard (MPSB), and by approximating the numerical curves identified the main contributions. We could follow the same approach for the mushroom billiard (MB), even though the process is complicated in this case by intermittent dynamics. In both billiards with mixed phase space we connected the chaotic phase space volume with $\lambda$ through $\kappa$ and showed that $\lambda$ has a leading contribution given by the inverse of the chaotic phase space volume (for the MPSB we used the chaotic portion instead of volume, because the total phase space volume does not depend on $B$).

To strengthen our point that a prominent contribution to the parameter dependence of the Lyapunov exponent in billiard is given by the inverse chaotic phase space volume, we present one final billiard, called the inverse stadium billiard (ISB) shown in Fig. 7(a), originally studied by Vörös et al. in Ref. Vörös et al. (2003). In this billiard a particle is propagating inside the stadium on straight lines, but after passing the boundary of the stadium it is subjected to a constant magnetic field, which brings the particle back inside the stadium, as depicted in Fig. 7(a). In the limit of infinite magnetic field the ISB recovers the fully chaotic Bunimovich stadium, for finite magnetic fields it has a mixed phase space. Here, we don’t want to analyze the ISB in any detail but only point out that also in this billiard the parameter dependence of its Lyapunov exponent is closely following the inverse of the chaotic phase space volume as shown in Fig. 7(b,c).

We want stress how different the mechanisms are that lead to chaos in the three different billiards. In the MPSB it is dispersing as well as the magnetic field. In MB it is defocusing and in the ISB it is even more involved. Yet in all three cases we find what is suggested by our toy model: the Lyapunov exponent has a leading contribution that is inverse to the chaotic phase space volume.

Because generally chaos in billiards arises via dispersing and defocusing collisions with curved boundary segments Bunimovich (2018), the Lyapunov exponent is necessarily inversely linked with the mean return time to these boundary segments. Furthermore, Kac’s lemma dictates that the mean return time of the chaotic trajectories to these boundaries is directly proportional to the chaotic phase space volume. For this reason we hypothesise that for most chaotic billiards with mixed phase space the Lyapunov exponent has a leading contribution inverse to the chaotic phase space volume.

This should straight-forwardly carry over to higher dimensions as well, since Kac’s lemma, the DPH framework, as well as our toy model, do not depend in any way on the dimensionality of the billiard. So far we didn’t find significant differences between billiards with sharply divided phase space (like the MB and the MPSB) and a fractal phase space structure (like the ISB). To conclude whether there are fundamental differences between sharply-divided and fractal phase spaces one will have to do more research. What we want to point out is that for fractal phase spaces it is much harder to estimate the volume of the chaotic set.

Appendix A Evolution of perturbation vector in a magnetic field

In their paper, Dellago, Posch and Hoover give two main equations to compute the evolution of perturbations in tangent space along a piece-wise smooth flow defined by the autonomous ODE system

[TABLE]

During smooth propagation, the perturbation vector $\delta\mathbf{\Gamma}$ along the trajectory $\mathbf{\Gamma}$ evolves according to

[TABLE]

If smooth propagation is interrupted at discrete times $t_{j}(\mathbf{\Gamma})$ by a discontinuous jump, represented here by a differentiable map $\mathbf{M}_{j}(\mathbf{\Gamma})$, the perturbation vector after the jump is given by

[TABLE]

where $\delta\tau_{c}=t_{j}(\mathbf{\Gamma}+\delta\mathbf{\Gamma})-t_{j}(\mathbf{\Gamma})$. For the case of elastic reflection with an obstacle, $\mathbf{M}_{j}$ can be written as

[TABLE]

Here $\mathbf{n}$ is the normal vector of the obstacle at the collision point, $\mathbb{I}_{2\times 2}$ and $0_{2\times 2}$ are the $2\times 2$ identity and zero matrices respectively, and $\quantity(\mathbf{a}\otimes\mathbf{b})_{jk}=a_{j}b_{k}$ is a second-order tensor.

A.1 Propagation

In magnetic billiards, particles propagate in circular arcs. To get the simplest possible set of equations of motion describing this mode of propagation, it is useful to introduce a phase angle $\theta=\arctan(p_{y}/p_{x})$. For uniform circular motion, this phase angle grows linearly in time as the system rotates with constant angular velocity $\omega$. This can be exploited to determine $\dot{\mathbf{p}}$ using the chain rule

[TABLE]

Using $\lVert\mathbf{p}\rVert=1$, one can easily calculate the explicit relation between $\mathbf{p}$ and $\theta$

[TABLE]

and combine the results of eqs. (44) and (45) to receive the equations of motion

[TABLE]

Using Eq. (41), we can now state the equations of evolution for a perturbation vector $\delta\mathbf{\Gamma}$ using the Jacobian $\mathbf{J}$ of $\mathbf{F}$.

[TABLE]

Using an exponential ansatz, we can compute the general solution of this system and receive a final result of

[TABLE]

where $\rho=1/\omega$ is the cyclotron radius.

A.2 Collisions

The derivation of the collision map for $\delta\mathbf{\Gamma}$ is largely analogous to the process used by DPH to derive their result for non-magnetic billiards. The equations of motion are assumed as stated in Eq. (46). For equation (42), we require the Jacobian matrix of $\mathbf{M}_{j}$, which is

[TABLE]

[TABLE]

By inserting eqs. (46), (43) and (62) into (42), we find

[TABLE]

It is now helpful to continue calculations for $\delta\mathbf{q}^{\prime}$ and $\delta\mathbf{p}^{\prime}$ separately. For the position component $\delta\mathbf{q}^{\prime}$, equation (63) can be written as

[TABLE]

For the momentum component $\delta\mathbf{p}^{\prime}$, we get

[TABLE]

This can be simplified by using the fact that $(\mathbf{b}\otimes\mathbf{a})\mathbf{c}=\left<\mathbf{c},\mathbf{a}\right>\mathbf{b}$ and introducing the quantity $\delta\mathbf{q}_{c}=\delta\mathbf{q}+\delta\tau_{c}\,\mathbf{p}$, which represents the real space difference vector between the collision points of satellite and reference trajectories.

[TABLE]

Using geometric considerations outlined by DPH we can now rewrite the penultimate term to get

[TABLE]

where $\phi$ is the angle of incidence, $\gamma_{R}$ is the local curvature of the obstacle and $\mathbf{e}$ and $\mathbf{e}^{\prime}$ are unit vectors orthogonal to $\mathbf{p}$ and $\mathbf{p}^{\prime}$ respectively.

A.3 Collision delay time

The quantity $\delta\tau_{c}$ in eqs. (42) and (63) can be interpreted as the time delay in between the collisions of the reference trajectory $\mathbf{\Gamma}$ and its satellite $\mathbf{\Gamma}+\delta\mathbf{\Gamma}$.

It can be computed by determining the signed distance from the satellite to its collision point measured along its trajectory at the time $t_{j}(\mathbf{\Gamma})$ of the collision of the reference particle. As we are considering the linearized dynamics of the perturbation and Eq. (42) is valid only to the first order of $\delta\mathbf{\Gamma}$, we will ignore all higher orders of $\delta\mathbf{\Gamma}$ in the subsequent calculations.

Furthermore, we will denote vector components in normal direction of the obstacle by a subscript $n$, i.e. $a_{n}=\left<\mathbf{a},\mathbf{n}\right>$. Similarly, the tangent component of a vector will be denoted by a subscript $t$.

Geometrically, one can immediately derive the following two relations from the two triangles highlighted in Fig. 8a

[TABLE]

where $\mathbf{c}$ is a vector between the obstacle and the cyclotron centre, as shown in fig. 8a. Eliminating the factor $\left<\mathbf{c},\mathbf{n}\right>$ yields

[TABLE]

By construction, $\theta$ cannot exceed $\frac{\pi}{2}$ in absolute value. Therefore, we can safely apply the arcsine function to eq. (70), receiving an expression for $\theta$.

[TABLE]

We can now use this to compute the angle $\alpha$ corresponding to the arc the particle has to travel during $\delta\tau_{c}$.

Using Fig. 8a, we can determine an expression for $\alpha$

[TABLE]

This can be further simplified by expressing $\phi$ in terms of the satellite particle’s momentum (compare fig. 8b).

[TABLE]

We can now linearize Eq. (73) to receive our final result for $\alpha$. As the individual terms are somewhat complicated, but very similar, we will treat them separately. The first-order Taylor expansion of leftmost term is

[TABLE]

where $\zeta$ is of the form $a\cdot{}\delta p_{t}+b\cdot{}\delta p_{n}$ with $a,b$ given by the respective partial derivatives. Expanding the rightmost term in (73) yields the similar result of

[TABLE]

This simplifies the final result for $\alpha$ significantly. Using the anti-symmetry of the arcsine function, we can see that most of (74) and (75) cancel out, leaving only

[TABLE]

Finally, we have to multiply equation (76) with the cyclotron radius $r$, then divide the resulting arclength by $\norm{p}$ to get $\delta\tau_{c}$.

As $\norm{\mathbf{p}}=1$ per convention, we can substitute $\sqrt{1-p_{t}^{2}}=\absolutevalue{p_{n}}$. However, we know that $p_{t}<0$ as the reference particle must have been moving towards the obstacle to collide with it. Therefore, we can further simplify our result to obtain

[TABLE]

This is the same result as for linear propagation in (Dellago et al., 1996, Eq. 18), since higher orders of $\delta\mathbf{\Gamma}$ were neglected.

Combining eqs. (64), (67) and (77), we receive the final result

[TABLE]

Appendix B Mean return time to stem

In this section, we will derive an analytic expression for the mean return time to the stem bottom in a mushroom billiard. Using Kac’s lemma Kac (1947); MacKay (1994); Meiss (1997), which states that for volume-preserving maps, the mean number of iterations $n_{S}$ required to return to a compact subset $S$ of phase space is given by

[TABLE]

where $\mu(\cdot)$ is the volume of a set and $A$ is the subset of phase space accessible to orbits originating in $S$.

To transform the billiard flow into a map, we discretize time in small steps $\Delta t$, implicitly considering the limit of $\Delta t\rightarrow 0$. We now choose the set $S$ of momenta and positions defined by

[TABLE]

where $Y$ is a box of width $w$ and height $\epsilon$ at the bottom of the stem. One should be careful about the choice of $S$. The simplistic approach of choosing the cap semicircle as the returning set (since the elementary growth segments are delimited by collisions with curved boundaries) will not yield the correct result. That is because the mean return time to the cap semicircle inherently includes contributions from both the periodic orbits of the MB as well as the laminar episodes of the elementary growth segments, which we have already shown to correspond to “free flight”-like motion.

The phase space volume of $S$ is $\mu(S)=\pi w\epsilon$ in the limit of $\Delta t\rightarrow 0$. As the chaotic phase space component of mushroom billiards is ergodic Bunimovich (2001), we know that the measure of the subset $A$ of phase space accessible from $S$ is given by the volume of chaotic phase space $V_{\text{CH}}=V_{\text{TOT}}-V_{\text{REG}}$.

Applying Kac’s lemma to get the mean iterations to return to the stem and multiplying by $\Delta t$, we get a result for the mean return time $\kappa_{S}(h,w,r)$ to the set $S$.

[TABLE]

To eliminate $\Delta t$, we divide by $\kappa_{S}$ for $w=2$ to get

[TABLE]

As this equation no longer depends on $\epsilon$, we can now take the limit $\epsilon\rightarrow 0$. The reason to use $w=2$ here is because the MB becomes the stadium billiard for $w=2$ (and thus is fully chaotic with no regular components with measure $>0$).

We can find $\kappa_{S}(h,2)$ because of the ergodicity of the MB for $w=2$. Specifically, it holds that $\kappa_{S}(h,2)=n_{S}(h)\times\tau$ with $n_{S}$ the mean amount of iterations to return to the stem and $\tau$ the mean collision time in the MB given by Eq. (3). We consider the boundary map of the billiard (Birkhoff coordinates), $(\xi,\sin{\phi_{n}})$, with $\xi$ the coordinate along the boundary (i.e. the arc-length) and $\phi_{n}$ the angle of incidence with respect to the normal vector at $\xi$. This coordinate system is a discrete mapping and Kac’s lemma applies directly. Therefore the mean iterations to return to the stem bottom are

[TABLE]

where $\absolutevalue{\partial Q}$ is the perimeter the boundary (the explicit factor of 2 represents the contribution of $\sin\phi_{n}$). Using Eq. (3), we find

[TABLE]

Combining eqs. (83) and (84), we receive an expression for the mean stem return time in the fully ergodic case

[TABLE]

This expression can be simplified by substituting the total phase space volume as defined in eq. (36), yielding

[TABLE]

Inserting this result into eq. (82), we receive

[TABLE]

Appendix C Mean duration of chaotic episodes

A chaotic episode in the MB as defined above consists of the particle travelling from the cap head directly into to the stem and back up to the cap head, without any other collisions inside the cap. To determine the mean duration $\tau_{c}$ of these episodes, we can geometrically determine the average lengths of the trajectories, exploiting that all trajectories are uniquely defined by the angle of incidence $\alpha$ and the distance $\delta$ from the cap center at which they enter the stem.

To simplify the calculation, it is useful to split the trajectory into the mean length of cap transit $\Braket{c}$ and the mean length of stem transit $\Braket{s}$ where $\tau_{c}=2\Braket{c}+2\Braket{s}$.

The stem transit length can be easily computed using simple trigonometry and depends only on the angle $\alpha$

[TABLE]

As the directions of particle momenta are equidistributed, we know that $\alpha$ must be cosine-distributed. We can now integrate over $\alpha$ to obtain the mean $s$, finding

[TABLE]

Determining the cap transit length is more difficult as it depends both on $\alpha$ and $\delta$. Using the law of sines, we can derive

[TABLE]

To get an average result, this expression has to be integrated both over $\delta$ and $\alpha$, again using the fact that $\alpha$ is cosine-distributed. Unfortunately, we were unable to solve the integrals analytically. Therefore, we decided to approximate $\int c\,d\delta$ by a polynomial before performing the second integration, yielding

[TABLE]

Appendix D Proof of Eq. (1)

We are interested in the flow of the PSB but to apply Kac’s lemma we need a discrete system. Thus, we obtain a map of the flow $\Phi^{t}$ by discretizing in time (similarly with appendix B), $f:=\Phi^{\Delta t_{\varepsilon}}$, with $\Delta t_{\varepsilon}\sim\varepsilon$. To prove eq. (1) we will apply Kac’s lemma to a set $\mathcal{S}$.

Let $\mathcal{T}_{\mathcal{S}}(B)=\Delta t_{\varepsilon}\times n_{\mathcal{S}}(B)$ (where $B$ is the magnetic field) denote physical recurrence time instead of map iterations. Let $\mathcal{W}$ be a circle of radius $r+\varepsilon$ concentric to the disk of the PSB and define the phase space subset $\mathcal{S}_{\varepsilon}$ such that

[TABLE]

where $\mathbf{\eta}(\mathbf{x})$ is the vector normal to $\mathcal{W}$. The mean collision time $\kappa$ of the PBS is exactly $\mathcal{T}_{\mathcal{S}}$ in the limit $\varepsilon\to 0$.

To find $\mu(\mathcal{S};B)$, the measure of set $\mathcal{S}$, we first realise that it does not depend on the magnetic field, $\mu(\mathcal{S};B)=\mu(\mathcal{S};0)=\mu(\mathcal{S})$. This is due to the infinitesimal width of $\mathcal{S}$, over which motion can always be approximated by a straight line for all finite magnetic fields values (i.e. equalling the magnetic field free case). Then, using (79) at $B=0$, we have

[TABLE]

because the PSB without magnetic field is fully ergodic and thus $\mu(A;B=0)=\mu(M)=1$ (here $M$ denotes the entire phase space whose measure we set for simplicity to 1). By substitution we get $\mathcal{T}_{\mathcal{S}}(B)=\mu(A;B)\times\mathcal{T}_{\mathcal{S}}(0)$. For small enough magnetic fields all chaotic orbits (and up to measure 0 only those) collide with the disks. In the limit $\varepsilon\to 0$ both $\Delta t_{\varepsilon}$ and $\mu(\mathcal{S}_{\varepsilon})$ go to 0 linearly with $\varepsilon$, therefore their ratio converges, i.e. $\mathcal{T}_{\mathcal{S}}\to\kappa$. Since by definition $\mu(A)=g_{c}$, the portion of chaotic orbits in the PSB (because we set the entire volume to have measure 1), we find that the mean collision time is given by the fraction of chaotic orbits as a function of the magnetic field $B$, times the mean collision time at $B=0$

[TABLE]