Chaos Detection by Fast Dynamic Indicators in Reflecting Billiards

The propagation of electromagnetic waves in a closed domain with a reflecting boundary amounts, in the eikonal approximation, to the propagation of rays in a billiard. If the inner medium is uniform, then the symplectic reflection map provides the polygonal rays’ paths. The linear response theory is used to analyze the stability of any trajectory. The Lyapunov and reversibility error invariant indicators provide an estimate of the sensitivity to a small initial random deviation and to a small random deviation at any reflection, respectively. A family of chaotic billiards is considered to test the chaos detection effectiveness of the above indicators.


Introduction
Electromagnetic cavities exhibit wave chaos that can be predicted by a semi-classical analysis and random matrix theory; see [1] for a two-dimensional open electromagnetic cavity, ref. [2] for closed three-dimensional cavities and [3] for a review.These powerful prediction tools hold true and can be extended to include coupling through antennas and wave-guides, as well as interconnected cavities [4].Furthermore, the semi-classical treatment of electromagnetic cavities has been tackled extensively in [5,6].The experimental verification of wave chaos in microwave billiards has been addressed by several international research groups over the last few decades.This extensive work led to the observation of a rich phenomenology originating from complex wave dynamics [7], including: Wave-function scars [8], chaotic dynamics in superconducting billiards [9], time-reversal symmetry breaking [10], nodal domains in rough billiards [11], time-invariance violation at and around exceptional points [12] as well as electromagnetic reverberation [13].Concerning optical resonators, it has been shown that chaotic dynamics are a key mechanism in asymmetric geometries [14][15][16].This effect is important in the generation of lasers [17] with directional emissions [18] of eccentric cavities.Chaotic cavities have been employed in a plethora of electromagnetic engineering applications, including time reversal energy focusing [19] and electromagnetic compatibility (EMC) testing [20].The equations for the propagation of acoustic and electromagnetic waves in resonant cavities are similar and coincide when the medium within the cavity is uniform.Usually, the cavities have a cylindrical or spherical symmetry and the resonant modes can be computed in a closed mathematical form.When the eikonal approximation is applicable, for sufficiently short wavelengths, the symplectic reflection map is integrable [21].Deformations of the boundary cause the loss of integrability and the emergence of chaotic behavior in the particles or rays trajectories [22].Even if the dynamics of integrable billiards are still an exciting research topic that continuously unveils novel phenomena, a growing interest has been developing within the scientific community for the properties of chaotic billiards [23][24][25][26][27].
For a sphere, the trajectory of any ray develops in a plane, which is determined by the initial ray position and velocity.The intersection of the sphere with the invariant plane is a circle, where the 2D area preserving the reflection map is integrable.
The symplectic 4D reflection map on a sphere is no longer integrable if the sphere is deformed.Thus, as opposed to a 4D flow, a first integral is no longer available and the 2D Poincaré map cannot be computed, preventing the analysis of the dynamical structures in phase plane portraits.For a cylinder, the trajectory of any ray develops in a plane, which intersects it on an ellipse, where the 2D reflection map is integrable.When the cylinder is deformed, the orbits no longer belong to a plane, the 4D symplectic map is not integrable and the phase portraits of the reflection map cannot be drawn [28].This is important also in wave propagation within deformed cavities as the evolution of quantum states show footprints of classical trajectories [29].
The Lyapunov error (E n (x)) and the reversibility error (E R n (x)) are dynamic indicators, based on the linear response theory, and allow the testing of the sensitivity of the orbits to small random deviations (see, for instance, [30]).For any fixed number of iterations, this sensitivity can be compared on a set of initial conditions chosen in a phase plane (see, for instance, [31]).
In Section 2, we present a brief introduction of the fast indicator used, their properties and relations.In Appendix B, we provide the definition of Lyapunov and reversibility errors.For a single orbit, the dependence on n can be investigated and the limit n → ∞, limit of n −1 log E n (x), provides the maximum Lyapunov exponent.E R n (x) has been shown to be very sensitive to multi-dimensional problems, such as chaos detection in planetary systems (see, for instance, [32]) and in 2D and 3D waveguides (see [33]).In Section 3, we present a numerical analysis of the 2D reflection map in a convex domain, given by a deformed circle.We compare the phase portraits with the color plots of the Lyapunov and reversibility errors computed in a regular grid of phase space.
The mathematical description of the convex billiard and its new parametrization are presented in Appendix A.
To conclude, we consider the transport of particles and rays within the billiard.Given a source of particles or rays within the billiard, the time evolution of the probability density of particles or the energy density of rays is analyzed.Such a density can hardly be determined analytically, even for integrable billiards, and a numerical strategy is presented.
The dynamic indicators could be evaluated at time t for the inner points of the billiard after performing an average with respect to the initial ray or particle direction.A detailed analysis of 2D and 3D billiards can be worked out by using the algorithms described here.

Lyapunov and Reversibility Error Indicators
The dynamic stability of the reflection map can be analyzed with the Lyapunov error indicator (i.e., E n (x)) or with the Reversibility error indicator (i.e., E R n (x)).The former measures the sensitivity of the initial conditions to small random deviations, and the latter measures the sensitivity to small deviations at each reflection.To implement the Reversibility error indicator, we iterate both the randomly perturbed and exact maps' n times, and compute the distance of the final point of the perturbed orbit with respect to the one obtained from the exact orbit.Letting be the amplitude of the random perturbations, we consider the linear response given by the → 0 limit.Given a symplectic or measure preserving map M(x), defined as a compact manifold of R 2d , we denote with DM(x) the tangent map where (DM) ij (x) = ∂ M i (x)/∂x j .The orbit x n is obtained by iterating n times the map M. We introduce the matrix L n obtained by taking the products of the tangent map along the orbit: For the initial condition y = x + ξ, where ξ is a random vector with a zero mean and unit covariance matrix ξ ξ T = I, we compute the orbit y n = M(y n−1 ) where y 0 = y.We compare the perturbed and the reference orbit by defining the linear response Ξ n according to The square of E n (x) is defined as the trace of the covariance matrix Σ 2 n (x) of the random vector Ξ n (x): The matrix L n L T n has the same eigenvalues as the Lyapunov matrix L T n L n but their eigenvectors are different.The invariants I (k) n (x) of the Lyapunov matrix are equal to the sum of the ( d k ) products of the eigenvalues and can be computed with the Faddeev-Leverrier algorithm.From a geometrical viewpoint, the invariants are equal to the sum of the squared volumes of the ( d k ) parallelotopes whose edges are L n (x)e j , and e j are any orthonormal base vectors.According to the Oseledec theorem, by writing the Lyapunov matrix as L T n (x)L n (x) = W n (x) e 2n Λ n (x) W n (x), the diagonal matrix Λ n (x) and the eigenvectors matrix W n (x) have a limit, for n → ∞ is equal to Λ(x) and W(x), respectively.As a consequence, letting Λ(x) = diag λ 1 (x), . . ., λ 2d (x) , we have If the map M is symplectic, the exponents are pairwise opposite λ d+j = −λ d−j+1 .If the map is integrable, then the first d eigenvalues of the Lyapunov matrix grow as n 2 , while the last d decrease as n −2 , so that all the Lyapunov exponents vanish.For a generic map, the power law growth or exponential growth with n of the first d invariants guarantee the classification of the phase space regions of regular and chaotic evolution.Below, we briefly introduce the reversibility error.For a complete description of the indicators and their properties, see Appendix B.

Reversibility Error
We consider the Backward-Forward process (BF) as one implementation of the Reversibility error indicator.In this case, we iterate n times the map with a random perturbation of amplitude first, and then the inverse map.The recurrence is given by and the linear response is defined by where L −1 k is the inverse of the matrix L k defined by (1).In the absence of perturbation, we are back to the initial condition y 2n ≡ x 2n = x.The covariance matrix of the random vector Ξ BF 2n is given by The square of the BF reversibility error E BF n is defined as the trace of Σ 2 BF n .The asymptotic limit of the invariants BF n (x) is the same as (4) for k ≤ d.The log of the last invariant BF n (x) is the Gibbs entropy of the process with the covariance matrix Σ 2 BF n (x) and its asymptotic behavior is the same as the Kolmogorov-Sinai entropy.
If the map is symplectic, then the matrices L n L T n and L T n L n are symplectic.Therefore, the trace of L T n L n and its inverse are equal.In this case, the BF reversibility error is simply related to the Lyapunov error E n as Previously, the BF process has been considered with the noise applied to both the B and F iterations.The covariance matrix of the linear response in in this case is given by and proof of (9) can be found in [3] Section 2.3.The asymptotic behavior of the invariants of the reversibility error covariance matrix Σ 2 nBF is determined by the positive Lyapunov exponents.Indeed, the limit of (2t) −1 log(I (k) nBF ) for t → ∞ is the sum of the positive exponents among the first k.As a consequence, 1  2 log(I nBF ), which corresponds to the Gibbs entropy of the BF random process where , is the sum of all the positive Lyapunov exponents (the first d for a symplectic map ), just as the Kolmogorov-Sinai entropy.Notice that the difference with the previous definition and ( 7) is negligible for a large n.Such a definition was initially proposed to compare the reversibility error due to a small random displacement with respect to the reversibility error due to the round off.The Reversibility Error indicator due to round off, denoted by REM, is a numerical implementation of E R n (x) and numerical equivalence has been proven for simple maps [34].Denoting with M the map evaluated with round off and M −1 the inverse map evaluated with the round off, the REM indicator is then The round off is a pseudo-random process in which just one realization is available.When different values for n or x are used, it is evident that REM BF n exhibits significant fluctuations with respect to E BF n , which is the result of an averaging process over the random displacements.No higher invariants can be defined for the reversibility error due to round off.We can also consider the specular implementation, namely the Forward-Backward REM FB , which consists of iterating n times first the perturbed inverse map M −1 and then the map M. For an autonomous symplectic map however, all the invariants are the same as for the BF process.Moreover, for a Hamiltonian flow, equivalence of the FB and BF invariants is a consequence of the time reversal invariance.
The linear approximation given by y n = x n + Ξ n can be considered; however, since the remainder of order 2 can generically be neglected only for a number of iterations, small with respect to log(1/ ), it is of no practical use.The linear response being based on the → 0 limit is valid for any number n of iterations.

Numerical Results for the 2D Billiard
In this section, we introduce numerical results for a 2D billiard; the reflection map for a convex billiard defined by (A12) with f (x) = x 2 /3 is presented.The billiard analyzed in this section has a closed boundary for ≤ 1/ √ 3 ∼ 0.577 as one can easily show.Indeed, letting V(x) = x 2 /2 + x 3 /3, the maximum occurs for x = −1/ and we require V(−1/ ) = 1/(6 2 ) ≥ 1/2 to have a closed boundary so that ≤ 1/ √ 3. We have restricted our analysis to ≤ /1/2, observing that for ≤ 0.1, the map is almost integrable; for = 0.2, the area of chaotic regions is significant and already for = 0.4, a large fraction of the phase space is chaotic.We have used the coordinates (φ, p) in the phase space even though the map preserves the measure but is not area preserving.We have compared the phase portraits with the Lyapunov error E n (x), the reversibility error E BF n and the round-off induced reversibility error REM BF n , computed according to (A28), (A38) and (A41), respectively.In Figure 1, we show the billiard for: It is important to notice that the indicators give additional information with regard to the phase portrait, since they provide a quantitative measurement of the chaos for an orbit.Additionally, for 4D billiards Poincare' sections are not available and the indicators are the only methods to investigate the stability of the phase space.
In Figures 2, we compare the phase portrait against the E n (x) for the billiard with = 0.1.The correspondence between the phase portrait and E n (x) color plot is good and in both cases, a thin layer of chaotic orbits is seen at the boundary of the main chains of islands.In the interior of the chains of islands, the error tends to be zero because the de-tuning (derivative of the frequency with respect to the action) is low.By approaching the separatrix, the de-tuning grows, diverging on the separatrix itself, where the error growth with n, changing from a power law to an exponential one.
In Figure 3, E R n (x) and REM BF n (x) show a similar behavior.REM is similar to the error induced by the small random displacements, but exhibits higher fluctuations because the averaging over the random process is missing.It is worth noticing that E R n 2 constitutes the sum along the orbit of (E n ) 2 , provided (s, p) are used as canonical coordinates.The Lyaponov error oscillates with n in the regions of quasi-integrable motion and for fixed n, oscillations are observed in phase space.These oscillations can be eliminated by using the MEGNO average (see, for instance, Ref. [35]) and they disappear for E R n (x).In Figures 4 and 5, the Poincare section, E R n (x), E n (x) and REM BF n (x) are shown for a larger perturbation = 0.2.The larger deformation with respect to the integrable billiard increases the area of the chaotic orbits.For = 0.3, almost one half of the unit phase space area-in the (φ 0 /2pi, p 0 ) initial coordinates used in the figures-is filled with chaotic orbits and this fraction increases approaching 1, when tends to the limit value = 1/ √ 3.For a 2D map, the phase space portraits provide the required information on the orbits stability; however, one advantage of the proposed indicators is that they provide a quantitative value to discriminate between regular and chaotic orbits.To explore small details in the phase space, one can analyze a smaller region and increase the number of iterations.Moreover, the indicators become the unique stability analysis tool when the dimensionality of the problem is increased.This is the case for the 3D billiard which leads to a 4D reflection map.In this case, the 2D phase portraits are no longer available.

Conclusions
We have considered the motion of particles on a 2D billiard with a convex reflecting boundary using a parametrization different from the one proposed by [36].Our implementation can be easily extended to a 3D billiard.The computation of the arc length s on the boundary requires a numerical integration, which can be avoided choosing the phase φ of the position vector r rather then the curvilinear abscissa s, even though, in this case, the map is only measure preserving.The stability of the orbits has been analyzed using the Lyapunov and Reversibility error fast indicators.Both indicators are invariant with respect to the choice of the initial deviation and of the orthogonal reference frame.We have shown that the logarithm of the second reversibility error invariant is the Gibbs entropy of the random vector defining the deviation from reversibility.The numerical results of the Lyapunov and reversibility error indicators, presented for a selected family of billiards, confirm the reliability of the proposed methods to explore the sensitivity of ray propagation to small random perturbations.
The minimum condition is satisfied for s = s 1 , where ψ = ψ 1 and setting v 1 = v(s 1 ) as the reflection condition ψ 1 = ψ 0 is fulfilled.Indeed, the angles that the rays form with the normal are θ 0 = π/2 − ψ 0 and θ 1 = π/2 − ψ 1 and the standard condition θ 1 = θ 0 is recovered.In order to determine the reflection map, we notice that p 0 = cos ψ = v 0 • τ 0 and p 1 = cos ψ 1 = v 1 • τ 1 are the momenta conjugated to s 0 and s 1 .To this end, we consider h(s 0 , s 1 ) as the generating function of the map.The moments p 0 and p 1 conjugated to s 0 and s 1 are given by As a consequence, the map M from x 0 = (s 0 , p 0 ) T to x 1 = (s 1 , p 1 ) T generated by h(s 0 , s 1 ) is area preserving and s, p are canonical coordinates.

Appendix A.1. Berry's
A parametrization r = r(φ) which allows the analytical computation of the curvilinear abscissa s = s(φ) was proposed in [36].The basic idea is to use ρ = ρ(φ) to represent the curve where ρ is the curvature radius and φ is the angle which the tangent forms with one of the coordinates' axes.We choose the origin on the curve at the point where x = x min and we move clock-wise.In the case of a circle of radius R, we have s = Rφ and the tangent and the inner normal are given by The functions r(φ) and ρ(φ) are periodic with period 2π and can be expanded in a Fourier series.Since r(0) = r(2π), Equation (A7) implies that the integral of ρ cos φ and ρ sin φ in [0, 2π] vanishes.The simplest deformation of the circle is given by but any trigonometric polynomial in cos(kφ) and sin(kφ) with k ≥ 2 might be considered.In this case, the curve parametrization becomes In order to find the map, we first determine the phase φ n+1 corresponding to the new reflection at r n+1 = r(φ n+1 ).To this end, we notice that the direction of the ray emerging at r n is given by the vector v n , forming an angle ψ n with τ n which in turn forms an angle φ n with the y axis.As a consequence, v n forms an angle φ n + ψ n with the y axis.Recalling that , the ratio of the x and y components is given by This equation implicitly determines φ n+1 as a function of φ n and ψ n .Next, we compute s n+1 = s(φ n+1 ) and r n+1 = r(φ n+1 ) via (A9).The angles entering the reflections at r n and r n+1 are shown in Figure A1.The last step is to determine ψ n+1 ; to this end, we consider the triangle whose vertices are P n , P n+1 and C n where the half lines directed along the normals n n and n n+1 intersect.Figure A1 shows the displacement vectors with respect to the origin being r n , r n+1 , r c .Taking into account the reflection condition, the angles of this triangle are π/2 − ψ n , π/2 − ψ n+1 and φ n+1 − φ n from which we obtain

. Alternative Parametrization and Reflection Map
We propose another equation for a deformed circle given by with > 0 and sufficiently small so that the billiard is convex.A parametrization r = r(φ) can be easily obtained but requires the solution of an implicit equation.This parametrization was already obtained in [38] and re-derived here.The advantage is that the present form can be easily extended to a deformed sphere.Indeed, at a given point, the curvatures of a surface depend on the two lines traced on it and therefore, on the corresponding tangents which are usually not orthogonal.It is not trivial to represent the curves by assigning the principal curvatures because this implies the knowledge of two families of curves that are mutually orthogonal.The parametrization we propose for the billiard boundary, defined by (A12), is the following: In this case, r(φ) is the distance from the origin and we denote by ρ(φ), the radius of curvature.Setting ṙ = dr/dφ, the curvilinear abscissa s(φ) is given by The tangent vector τ , the inner normal vector n, and the radius of curvature ρ are given by The expressions of τ and n are still relatively simple, unlike the radius of curvature, which never enters in the computation of the reflection map.Notice that φ is no longer the phase of the vector τ with respect to the y axis.Denoting with φ τ the phase of τ , which we choose to be 0 when τ = (1, 0) and increasing clock-wise as φ, we have As a consequence, given τ , its phase is given by When the equation of the billiard boundary is given by (A13), we determine r = r(φ) and ṙ = ṙ(φ) by solving the equations iteratively staring with r = 1 and ṙ = 0.At first order, r(φ) = 1 − f (− cos φ) + O( 2 ) .The solution r = r(φ) is obtained iteratively and the convergence is very fast for << 1.
The iterations compute the fixed point stops when machine accuracy is reached.Then, s = s(φ) is obtained by numerical integration according to (A14).
To compute the reflection map, we start from a point r 0 with ray velocity v 0 .At the point r 0 , the tangent and inner normal vectors are τ 0 and n 0 .We define the phase of a vector b with respect to a vector a as the angle between a and b counted clock-wise and its range is [0, 2π].To the initial position, we associate, according to (A13), the angle φ 0 defined as the phase of −r 0 with respect to the x axis.To the initial tangent vector, we associate, according to (A18), the angle φ τ 0 defined as the phase of τ 0 with respect to the y axis.If y 0 = 0 and x 0 < 0, we have φ 0 = φ τ 0 = 0.The phase of v 0 with respect to τ 0 is ψ 0 and its range is [0, π] whereas the range of φ 0 and φ τ 0 is [0, 2π].As a consequence, the phase of v 0 with respect to the y axis is φ τ 0 + ψ 0 mod 2π.The tangential component of the velocity is p 0 = v 0 • τ 0 ≡ cos ψ 0 and if p 0 > 0, the reflection is clock-wise; if p 0 < 0, it is counter clock-wise.Alternatively, we can start with (φ 0 , p 0 ) and determine the vectors r 0 , v 0 .
After n iterations, the point r n on the boundary is reached and the velocity of the outgoing ray is v n .The angles φ n , φ τ n , ψ n and the momentum p n = cos ψ n are determined and eventually the curvilinear abscissa s n is obtained with a numerical integration.each reflection.In the latter case, it is convenient to perturb the orbit during the first n iterations and integrate the unperturbed inverse map.
These indicators are computed according to the linear response theory, namely when the limit of the perturbation tends to 0. For this reason, the linear approximation is valid only for a limited number of iterations n.The number n of iterations up to which the linear approximation is valid depends on x and and uniform estimates are typically n ∼ log(1/ ).The linear response is valid for any iteration number N since it involves the limit to zero of the noise amplitude .Given a symplectic or measure preserving map M(x), defined in R 2d , or a compact manifold of dimension 2d, we denote with x n the point of the orbit obtained by iterating the map n times, starting from an initial condition x, and by L n (x), its tangent map.Denoting the composition of the map with M Consider a nearby initial condition y = x + ξ where ξ is a random vector with a zero mean and unit covariance matrix ξ ξ T = I.The orbit y n = M(y n−1 ) is compared with the unperturbed orbit x n = M(x n−1 ) and the linear response is given by a random vector Ξ n defined by The square of E n (x) is defined as the trace of the covariance matrix: The matrix L n L T n has the same eigenvalues of the Lyapunov matrix L T n L n but the eigenvectors are different.The invariants products of the eigenvalues.From a geometric viewpoint, the invariant I (k) n is the sum of the squared volumes of the ( d k ) parallelotopes, whose edges are all the k distinct choices among the vectors L n e 1 , . . ., L n e 2d , having denoted with e j the base vectors (e j ) i = δ ij .From the polar decomposition L n = R n e n Λ(n) W T n , where R n and W n are orthogonal matrices, and the fact that Λ(n) is a real diagonal matrix whose elements are λ If the matrix is symplectic, letting r T = (q, p) with q = (q 1 , . . ., q d ) and p = (p 1 , . . ., p d ), 2 and the next invariants can be recursively computed with the Faddeev-Leverrier formula, involving the traces of the powers of L T n L n .The last invariant is given by the determinant.
The square of the E R n (x) is defined as the first invariant, namely, the trace of the covariance matrix Σ BF n : The next invariants of Σ BF n 2 can be computed.If the map is symplectic, then L n L T n and L T n L n are symplectic matrices.Therefore, the traces of L T n L n and its inverse are equal and the RF reversibility error is given by a quadratic average of the E n (x).
Previously, the BF process has been considered with the noise applied to both the backward and forward iterations.In this case, the covariance matrix is given by A proof of (A40) can be found in [3] Section 2.3.This definition, whose difference with the previous one (A38) is negligible for large n, was initially proposed to compare E R n (x), due to a small random displacement, against REM.Denoting with M the map and with M −1 its inverse, both computed by including the round-off, then the REM is defined by The round-off is a pseudo-random process of which just one realization is available.It is evident that, as opposed to E BF n (x), REM BF n , it exhibits significant fluctuations when n or x vary.No higher invariants can be defined for the reversibility error due to round-off.
The specular process FB, in which we first iterate n times the perturbed inverse map M −1 and then the map M, can be considered.The previous formulae hold where L n (x) = DM n (x) is replaced with L −n (x) = DM −n (x) = DM n (M −n )) = L n (x −n ) −1 and therefore, the covariance matrix for the FB process is This formula can be proved just as the formula (A40) for the BF covariance matrix.One should take into account exchanging F and B amounts to change the map M with M −1 so that L T n L n is replaced by L T −n L −n since L n = DM •n and L −n = DM • −n ≡ D(M −1 ) •n .For a symplectic map, the BF and FB covariance matrix are equivalent, provided that one exchanges the sign on the momentum, whereas for a dissipative map, they are intrinsically different.In this case, their different behavior reflects the break-up of the time reversal invariance for the unperturbed system.As a consequence, for a symplectic map, the FB reversibility error is related to the Lyapunov error by The corresponding round-off induced reversibility error REM is defined by iterating n times the inverse map M −1 first, and then the map M .If x belongs to an ergodic component, then the asymptotic behavior of E n (x) is the same for any x except for a subset of zero measure.As a consequence, the asymptotic of E BF n (x) and E FB n (x) is the same almost everywhere.

Figure 1 .
Figure 1.Left panel: billiard boundary defined by (A12) with f (x) = x 2 /3 and =0, 0.1, 0.2, 0.4, 0.5 (green, cyan, purple, grey, blue).Right panel: rays trajectory for = 0.5 and n = 10 (purple line) and the reversed ray trajectory (red line).At the initial point, the tangent and inner normal vectors are shown (black lines).We have analyzed the orbits for the billiard, comparing them with the color map of the Lyapunov error E n (x) , reversibility error E BF n (x) and the round-off induced reversibility error REM BF n (x), where x = (φ, p) T is chosen on a regular grid of the N g × N g points.A logarithmic color map is used to show the results.By defining the tangent vector τ (s) = ∂r(s)/∂s and the ray velocity as v(s) , then the range of the phase φ of the position vector is [0, 2π], and the range of the momentum p = τ • v is [−1, 1].Since the orbits are symmetric, by changing p into −p, the chosen range for p is then in the range [0, 1].For a deeper discussion on the properties of the billiard, see Appendix A.It is important to notice that the indicators give additional information with regard to the phase portrait, since they provide a quantitative measurement of the chaos for an orbit.Additionally, for 4D billiards Poincare' sections are not available and the indicators are the only methods to investigate the stability of the phase space.In Figures 2, we compare the phase portrait against the E n (x) for the billiard with = 0.1.The correspondence between the phase portrait and E n (x) color plot is good and in both cases, a thin layer of chaotic orbits is seen at the boundary of the main chains of islands.In the interior of the chains of islands, the error tends to be zero because the de-tuning (derivative of the frequency with respect to the action) is low.By approaching the separatrix, the de-tuning grows, diverging on the separatrix itself, where the error growth with n, changing from a power law to an exponential one.In Figure3, E R n (x) and REM BF n (x) show a similar behavior.REM is similar to the error induced by the small random displacements, but exhibits higher fluctuations because the averaging over the random process is missing.It is worth noticing that E R

Figure 2 .
Figure 2. Left panel: phase portrait of the billiard with = 0.1.Each orbit is computed for n = 1000 and the initial points of each orbit are the red dots.Right panel: Lyapunov error E n (x) where φ 0 /(2π) and p are chosen in a regular grid of the unit square with N g × N g points and N g = 200.The iterations number is n = 200 and the results are shown in a logarithmic color scale.

Figure 3 .
Figure 3. Left panel: billiard with = 0.1 reversibility BF error E BF n (x) where φ 0 /(2π) and p are chosen in a regular grid of the unit square with N g × N g points and N g = 200 and the iterations number is n = 200.Right panel: round-off induced BF reversibility error REM BF n (x).

Figure 4 .
Figure 4. Left panel: phase portrait of the billiard with = 0.2.Each orbit is computed for n = 1000 and the initial points of each orbit are the red dots.Right panel: Lyapunov error E n (x) where φ 0 /(2π) and p are chosen in a regular grid of the unit square with N g × N g points and N g = 200.The iterations number is n = 200 and the results are shown in a logarithmic color scale.

Figure 5 .
Figure 5. Left panel: billiard with = 0.2 reversibility BF error E BF n (x) where φ 0 /(2π) and p are chosen in a regular grid of the unit square with N g × N g points and N g = 200 and the iterations number is n = 200.Right panel: round-off induced reversibility error REM BF n (x).

(
Figure A1.Angles and vectors defining the ray propagation.