Propagation of rays in 2D and 3D waveguides: a stability analysis with Lyapunov and Reversibility fast indicators
PPropagation of rays in 2D and 3D waveguides: a stability analysiswith Lyapunov and Reversibility fast indicators
G. Gradoni, a) F. Panichi, b) and G. Turchetti c) School of Mathematical Sciences and Department of Electrical and Electronic Engineering, University of Nottingham,United Kingdom d)2)
Department of Physics and Astronomy , University of Bologna, Italy Department of Physics and Astronomy, University of Bologna, Italy. INDAM National Group of Mathematical Physics,Italy (Dated: 13 January 2021)
Propagation of rays in 2D and 3D corrugated waveguides is performed in the general framework of stability indicators.The analysis of stability is based on the Lyapunov and Reversibility error. It is found that the error growth follows apower law for regular orbits and an exponential law for chaotic orbits. A relation with the Shannon channel capacity isdevised and an approximate scaling law found for the capacity increase with the corrugation depth.
We investigate the propagation of a ray in a 2D wave guidewhose boundaries are two parallel horizontal lines, witha periodic corrugation on the upper line. The reflectionpoint abscissa on the lower line and the ray horizontalvelocity component after reflection are the phase spacecoordinates and the map connecting two consecutive re-flections is symplectic. The dynamic behaviour is illus-trated by the phase portraits which show that the regionsof chaotic motion increase with the corrugation amplitude.For a 3D wave-guide the 4D map connecting two consec-utive reflections on the lower plane is symplectic, but itsorbit cannot be examined by looking at the intersectionswith a 2D phase plane, since a continuous interpolation ofthe orbits is not available. In this case the fast dynamicindicators allow to perform a stability analysis. For eachpoint of a grid in a 2D phase plane one computes the orbitfor a chosen number of iterations and the correspondingvalue of the fast indicator, which is conveniently visual-ized using a color plot. After the fast Lyapunov indicator,many other indicators have been introduced. Our anal-ysis is based on the Lyapunov error (LE), due to a smallrandom initial displacement and the reversibility error oc-curring when the orbit is reversed in presence of a smalladditive noise (RE) or round off (REM). For integrablemaps the growth of LE and RE follows a power law andfor quasi integrable maps the same growth occurs close toa stable fixed point. More generally the error growth fol-lows a power law for regular orbits and an exponential lawfor chaotic orbits. There is numerical evidence that REMgrows as RE though with large fluctuations. The channelcapacity is related to LE and the dependence of its phasespace average on the corrugation amplitude is considered.These indicators confirm their reliability for the stabilityanalysis of the ray propagation in a 2D and 3D wave guide,providing a measure of the sensitivity of the orbits to ini-tial conditions, noise and round off. a) Electronic mail: [email protected]. b) Electronic mail: [email protected]. c) Electronic mail: [email protected]. d) Also at Maxwell Centre, University of Cambridge, United Kingdom.
I. INTRODUCTION
The equivalence between geometrical optics and mechanicswas established in a variational form by the principles of Fer-mat and Maupertuis. If a ray propagates in a uniform mediumwith a reflecting boundary, then the trajectory is the same as aparticle freely moving and elastically colliding with the sameboundary. As consequence a wave-guide and a billiard areequivalent optical and mechanical systems . Since the veloc-ity of the particle does not change, we can assume it has unitmodulus, just as the ray velocity normalized to the speed oflight. The trajectory is determined by the collision points andthe velocity direction after the collision. The billiards withpolygonal or smooth convex closed boundaries have been in-tensively investigated and the mathematical literature is veryrich, see .We first consider the 2D wave-guide whose boundary is astraight line and a corrugated parallel line. The trajectory isa polygonal line specified by the abscissa x of the reflectionpoints on the straight line and the parallel component v x of thevelocity after reflection. The Fermat principle establishes that,given two points x , x on the lower straight line, the ray, col-liding once with the corrugated line, follows the path of mini-mal length and the collision is just a reflection. In addition theray minimal path length h ( x , x ) is the generating function ofthe area preserving map M connecting two subsequent reflec-tions on the straight line. A similar procedure allows to obtaina symplectic map for the ray propagation in a 3D waveguidemade of an horizontal plane and an upper parallel plane witha periodic corrugation. If the medium between the waveguideboundaries if not uniform the piecewise linear path joiningtwo consecutive reflection points on the lower line or plane isreplaced by the geodesic with respect to the metric n ds , where n is the refraction index.The transition from ordered to chaotic motion in 2D bil-liards and waveguides has been considered ; the transportand diffusion properties have been extensively analysed .The stability properties of the map depend on the corruga-tion amplitude. For the 2D wave guide the phase portrait a r X i v : . [ n li n . C D ] J a n of the corresponding 2D map allows to detect the regions ofregular and chaotic motion. Finite time indicators such asthe Fast Lyapunov Indicator (FLI) have been first pro-posed to analyze the orbital stability. Other short term indi-cators of variational nature have been introduced: the Align-ment Indices (SALI) , the Orthogonal Fast Lyapunov Indi-cator (OFLI) , the Mean Exponential Growth of Nearby Or-bits (MEGNO) , which is an excellent filter of the oscil-lations of FLI, the relative Lyapunov indicator (RLI) , andthe Generalized Lyapunov indicators (GALI), whose asymp-totic behaviour is related the all the Lyapunov exponents .An indicator based on the distribution of the stretching num-bers (SSN) was also proposed . An extensive numericalcomparison of previous indicators for symplectic maps waspresented . Spectral indicators based on the Fourier anal-ysis have also been proposed and extensively used . A ge-ometric chaos indicator based on the Riemann curvature ofthe constant energy manifold and a 0-1 test of chaos have also been proposed. More recently the Lyapunov andreversibility errors have been introduced to measure thesensitivity of the orbits to a small random initial displacementand to a small additive noise along the orbit. The Lyapunoverror (LE) and the reversibility error are simply related sothat the computation of RE does not require expensive Monte-Carlo procedures. The relevant difference of LE withe re-spect to previous variational indicators such as FLI, is that LEdoes not depend on the initial deviation vectors. A full set ofLyapunov error indicators (LEI), whose asyptotic behaviouris related to all the Lyapunov exponents, just as GALI, andwhich independent from the initial deviation vectors, has beenintroduced jointly with reversibility error indicators (REI) .The last REI corresponds to the reversibility entropy, andits asymptotic behaviour is governed by the sum of positiveLyapunov exponents, as the upper-bound to the Kolmogorov-Sinali entropy .The LE and RE grow following a power law, for the regu-lar orbits of quasi integrable symplectic maps as the previousvariational indicators, exponentially for chaotic orbits. The re-versibility error due round off (REM) was first introduced in and its features were examined in . For previous works onthe round effect in the computation of orbits for Hamiltoniansystems see . Previous numerical investigations of Hamil-tonian systems and symplectic maps confirm that LE growslinearly with oscillations, (due to the loss of rotational sym-metry when the coordinates are not normal) for regular orbits,whereas the growth of RE is almost oscillations free. Neglect-ing its large fluctuations REM is comparable with RE, eventhough no rigorous proof is available.For the 3D waveguide no direct inspection of the orbits is pos-sible, since the section of the orbits of the symplectic 4D mapwith a 2D plane would require a continuous interpolation,available only when an interpolating Hamiltonian is known.Normal forms provide the interpolating Hamiltonian for quasiintegrable symplectic maps, but their recursive computation ispossible just for polynomial maps, and in addition the inter-polation is not exact due to the presence of a non integrableremainder. The variational indicators, computed for the orbitsissued from the points of a regular grid in a 2D phase plane for the 4D map, and visualized with a colour plot, allow to deter-mine the stability properties just as for the 2D map. We haveanalyzed only LE and RE for a limited number of iteration.Indeed as stated Froeschlé et al. the variational indicators asFLI computed for short times exhibit some dependency on theinitial conditions of the deviation vectors. Since LE, RE andREM do not depend on the initial deviation vectors our choiceis justified. A careful investigation of the sticky chaos, whichrequires longer orbits, and more extensive numerical explo-ration of the reflection maps based on the full set of invariantindicators LEI and REI and a comparison with the standardvariational indicators, will be the object of a future work.Only a few hundred iterations of the map are required to ob-tain a reliable stability picture, unless one is interested in thedetails of a small region. To obtain the analytic form of thetangent map is rather cumbersome for the 3D waveguide case.The shadow orbit method provides a simple, though less ac-curate, alternative which consists in evaluating the orbits forinitial conditions with small displacements along an orthogo-nal basis, which amounts to replace the partial derivatives withfinite differences. For any initial condition, 4 additional eval-uations of the orbit are required (2 for the 2D map) in orderto obtain LE and RE, whereas just 1 is required to computeREM.We have analyzed a model of 2D waveguide for different val-ues of the corrugation amplitude, showing that the LE, RE,REM provide comparable results, which describe the orbitssensitivity to a small initial random displacement, to a noisealong the orbit and to round off. For a model of 3D waveg-uide the same error plots, for initial conditions on 2D phaseplanes, exhibit a similar behaviour, though the structure isricher with respect to the 2D waveguide, due to the presenceof the Arnold web of resonances The effectiveness of the pro-posed method, already experienced in celestial mechanics and beam dynamics models, is confirmed in these examplesof 2D and 3D waveguides.The paper is organised as follows. In section 1 we presentthe variational derivation of the reflections map for a 2Dwaveguide. In section 3 the extension to the 3D waveguideis outlined. In section 4 our dynamical indicators are definedand their basic properties are illustrated. In section 5 and 6the numerical results of these indicators show how the insta-bility regions grow when the corrugation amplitude increases.In section 7 the link between our indicators and the channelcapacity is established and it is shown how its phase space av-erage increases with the corrugation amplitude. Conclusionsand perspectives are presented in section 8. II. THE 2D WAVEGUIDE
Given two parallel reflecting lines z = z = ( x , z ) plane, the light ray direction after each reflection onthe lower line is the same being specified by the unit vector v = ( x x , v z = (cid:112) − v x ) T , where T denotes the transpose of amatrix or a vector. The time τ between two reflection on theupper and lower line is given by τ v z =
1. We choose as phasespace coordinates ( x , v x ) so that the sequence of reflections ( x n , v xn ) on the lower line is given by v xn + = v xn x n + = x n + τ v xn = x n + v xn (cid:112) − v xn (1)This map is integrable and area preserving. If we let v =( cos θ , sin θ ) T where 0 ≤ θ ≤ π the map becomes θ n + = θ n x n + = x n + θ n (2)and preserves the measure d µ = sin θ dxd θ . If our phasespace is a cylinder T ([ − π , π ])] × [ − , ] rather than the in-finite strip R × [ − , ] , then x is an angle variable and v x anaction variable.The frequency Ω = v x ( − v x ) − / diverges for v x → ± H = − (cid:112) − v x . Since also Ω (cid:48) diverges when v x → ±
1, any perturbation renders the mapchaotic as for an integrable near the the separatrix where Ω vanishes but its derivative diverges. The 2D waveguide is ob-tained by corrugating the upper line according to z = + ε f ( x ) (3)where f ( x ) is a periodic function of period 2 π such that f ( x ) > − ≤ ε ≤
1. The curvilinear abscissa s ( x ) ofa point Q on the corrugated line of coordinates ( x , z ) where z is given by (4) s ( x ) = (cid:90) x (cid:113) + ε f (cid:48) ( x (cid:48) ) dx (cid:48) (4)Consider a ray which starts from P = ( x , ) and reaches cor-rugated line at Q = ( x , z ) . After reflection the ray reaches the x axis at the point P = ( x , ) . Keeping P and P fixed andletting Q vary the path length H of the segments P Q and QP is a function of s ( x ) H ( s ) = h ( x , s ) + h ( x , s ) h ( x , s ) = (cid:113) ( x − x ) + ( + ε f ( x )) h ( x , s ) = (cid:113) ( x − x ) + ( + ε f ( x )) (5)where x = x ( s ) is the inverse of the function s = s ( x ) definedby (4). Referring to Fig. 1 we define ψ and ψ (cid:48) the angleswhich the velocities v and v of the incoming and outgoingray at Q forms with the tangent τ . The angles between thevectors − v and v and the normal ν at Q are π / − ψ and π / − ψ (cid:48) . As a consequence v = ( x − x , + ε f ( x )) T h ( x , s ) v = ( x − x , − − ε f ( x )) T h ( x , s ) τ = ( , ε f (cid:48) ( x )) T (cid:113) + ε f (cid:48) ( x ) ν = ( ε f (cid:48) ( x ) , − ) T (cid:113) + ε f (cid:48) ( x ) (6)We denote with v the velocity of the ray reflected at P andwith θ and θ the angles v and v form with the positive FIG. 1: Geometry of the reflections. On the right frame for n = α = π / − ψ , where ψ is defined on the leftframe. x axis, see (1). The derivatives of h ( x , s ) and of h ( x , s ) aregiven by ∂∂ s h ( x , s ) = x − x + ( + ε ( f ( s )) ε f (cid:48) ( x ) h ( x , s ) (cid:113) + ε f (cid:48) ( x ) == v · τ = cos ψ∂∂ x h ( x , s ) = x − xh ( x , s ) = − v x = − cos θ ∂∂ s h ( x , s ) = x − x + ( + ε ( f ( s )) ε f (cid:48) ( x ) h ( x , s ) (cid:113) + ε f (cid:48) ( x ) == v · τ = − cos ψ (cid:48) ∂∂ x h ( x , s ) = x − xh ( x , s ) = v x = cos θ (7)The stationary point H with respect to s , when x and x arekept fixed, is met for s = s ∗ ( x , x ) where ψ = ψ (cid:48) , whichcorresponds to the the reflection condition. We introducethe function F ( x , x ) equal to H at the stationary point s = s ∗ ( x , x ) and compute its differential F ( x , x ) = h ( x , x , s ∗ ) + h ( x , x , s ∗ ) dF ( x , x ) = − v x dx + v x dx (8)Equation (8) shows that F ( x , x ) is the generating functionof the canonical transformation M from ( x , v x ) to ( x , v x ) .After n iterations the phase space point ( x n , v xn ) is reached andis mapped into ( x n + , v xn + ) by M . In the physical space theray issued from the point P n = ( x n , ) with horizontal velocity v xn is reflected at Q n = ( x n , + ε f ( x n )) and reaches P n + =( x n + , ) where the horizontal velocity after reflection is v xn + see figure 1.In order to obtain the map we first compute the time τ needed for the transit between P n and Q n . Then we determinethe horizontal component of the velocity v xn + of the ray out-going from P n + which is equal to v x ∗ , where the velocity v ∗ of the ray from Q n to P n + is determined by the reflection con-dition v ∗ = v n − ν ( ν · v n ) and the normal of the guide at Q n is given by (6) . Finally x n + − x n is given by v xn τ + v xn + τ (cid:48) where the transit time from τ (cid:48) from Q n to P n + is given by v zn + τ (cid:48) = v zn τ . v zn τ = + ε f ( x n + τ v xn ) , v zn = (cid:113) − v xn v xn + = v xn − ν x ( ν · v n ) == v xn − ε f (cid:48) v xn − ε f (cid:48) v zn + ε f (cid:48) f (cid:48) = f (cid:48) ( x n + τ v xn ) x n + = x n + τ (cid:18) v xn + v xn + v zn v zn + (cid:19) v zn + = (cid:113) − v xn + (9)The map ( x n + , v xn + ) = M ( x n , v xn ) is symplectic since it isimplicitly defined by the generating function F ( x n , x n + ) ac-cording to v xn = − ∂ F / ∂ x n and v xn + = ∂ F / ∂ x n + . The com-putation of the tangent map, is given in the Appendix. III. THE 3D WAVE-GUIDE
In this case we consider two planes z = z = v which can be written inCartesian coordinates according to v = ( v x , v y , v z ) v z = (cid:113) − v x + v y (10)With v we denote the velocity of the ray reflected by the lowerplane z = v z >
0. The time τ between a reflectionon the lower and upper plane is τ v z = ( x , y , v x , v y ) as phase space coordinates, the map between two consecutive reflections reads v xn + = v xn v yn + = v yn x n + = x n + τ v x , n = x n + v xn v z ny n + = y n + τ v y , n = x n + v yn v z n (11)If the map is defined on T ([ − π , π ]) × [ − , ] , then v v , v y are the action variables and the interpolating Hamiltonian is H = − (cid:113) − v x − v y . In the 3D waveguide the plane z = z = z ( x , y ) = + ε f ( x , y ) (12)where f ( x , y ) is a periodic function in x and y with pe-riod 2 π and f > − < ε <
1. As a consequencethe map describing the rays propagation can be defined on T ([ − π , π ]) × [ − , ] . As in the 2D case we consider thesequence of points P , P , . . . , P n , . . . on the z = Q n the point on the upper corrugated plane hitby the ray issued from P n , which after reflection reaches the z = P n + . The sequence points P , P (cid:48) , . . . , P (cid:48) n , . . . on the torus T is simply obtained taking the modulus withrespect to 2 π so that ( x (cid:48) n , y (cid:48) n ) ∈ [ − π , π ] . Letting ν ( x , y ) bethe normal at the corrugated surface at the reflection point Q = ( x , y , z ) = + ε f ( x , y ) , which explicitly reads ν ( x , y ) = ( ε f x , ε f y , − ) T (cid:113) + ε ( f x + f y ) (13)The map ( x n + , y n + , v xn + , v yn + ) = M ( x n , y n , v xn , v yn ) spec-ifying the ray trajectory is given by τ v zn = + ε f ( x n + τ v xn , y n + τ v yn ) v zn = (cid:113) − v xn − v yn v xn + = v xn − ν x ( ν · v n ) ν = ν ( x n + τ v xn , y n + τ v yn ) v yn + = v yn − ν y ( ν · v n ) x n + = x n + τ (cid:18) v xn + v xn + v zn v zn + (cid:19) y n + = y n + τ (cid:18) v yn + v yn + v zn v zn + (cid:19) v zn + = (cid:113) − v xn + − v yn + (14)The first equation implicitly defines the propagation time τ from P n to Q n , the remaining equations define the map M which is symplectic. Indeed starting from the Fermat vari-ational principle it is shown in appendix B that this map isobtained from a generating function F ( x n , y n , x n + , y n + ) . Wedo not quote in this case the expression of the tangent mapsince it is rather involved. IV. DYNAMICAL INDICATORS
We will present a numerical analysis of the 2D and 4D sym-plectic maps whose orbits describe the propagation of a rayin a 2D and 3D periodic waveguide. The stability analysisis intended to discriminate the regions of regular and chaoticmotions using the recently introduced fast stability indicators,denominated Lyapunov and Reversibility error. The first oneis closely related to fast Lyapunov indicator first proposed toanalyze the growth of a small initial displacement. We shallnot present a comparison with other fast indicators, because ithas been performed for other models and since our indicatorsare independent from the initial displacement.
A. Lyapounov error
For a symplectic map in a phase space of dimension 2 d map the Lyapunov error (LE) describes the growth of an initialrandom displacement ε ξ where ξ ∈ R is a random vectorwith zero man and unit variance (cid:104) ξ (cid:105) = (cid:104) ξ ξ T (cid:105) = I (15)in the zero amplitude limit ε →
0. Denoting with (cid:4) n the ran-dom displacement after n iteration and with DM ( x ) the tan-gent map Ξ n ( x ) = lim ε → M n ( x + εξ ) − M n ( x ) ε = L n ( x ) ξ L n ( x ) = DM n ( x ) = D M ( x n − ) L n − ( x ) (16)The square of the Lyapunov error E Ln ( x ) is defined as as thevariance of of the random vector Ξ n or the trace of its covari-ance matrix E Ln ( x ) = (cid:104) Ξ n · Ξ n (cid:105) = Tr ( Σ n ) Σ n = (cid:104) Ξ n Ξ Tn (cid:105) = L n L Tn (17)We might define the error E Ln ( x , η ) for a given initial dis-placement εη with (cid:107) η (cid:107) =
1, when ε →
0. Its logarithm isjust the Fast Lyapunov Indicator. We observe that letting e j be any orthonormal base the sum of E Ln ( x , e j ) extended to allthe vectors of base is equal to E Ln ( x ) . Our definition is inde-pendent from the initial displacement. An expression for theLyapunov error equivalent to (17) is given by E Ln ( x ) = Tr ( L Tn ( x ) L n ( x )) = Tr (cid:16) ( DM n ( x )) T DM n ( x ) (cid:17) (18)The stability analysis is performed by fixing n and observingthe change of E Ln when x varies on a 2 D phase plane (or a 2D manifold). Oseledet theorem states that if x belongs toan ergodic component then ( L Tn L n ) / n has a limit W e Λ W T independent from x , where W is an orthogonal matrix, Λ isdiagonal and its entries λ ≥ λ ≥ . . . ≥ λ d are the Lyapunovexponents with λ d + k = − λ d − k + and λ k ≥ ≤ k ≤ d .As a consequence asymptotically E Ln (cid:39) e n λ + e n λ + . . . + e n λ d (cid:39) e n λ and more rigorously lim n → ∞ ( E Ln ) / n = λ .However we are not interested in the n → ∞ limit but ratheron the dependence on x for a finite (possibly large) value of n ,when x varies in phase space. B. Reversibility error
We consider now the iteration of the map n times followedby the iterations of the inverse map still n times. Inserting asmall additive noise at each forward and backward iteration,the final displacement with respect to the initial condition is astochastic process, whose variance defines the square on thereversibility error (RE). More precisely letting x ε , = x be theinitial conditions, we consider the noisy orbit x ε , k = M ( x ε , k − ) + εξ k Ξ k = lim ε → x ε , k − x k ε (19)for k = , . . . , n . Starting from x ε , n we consider the reversednoisy orbit x ε , n , − k = M − ( x ε , n − k + ) + εξ − k Ξ n , − k = lim ε → x ε , n , − k − x n − k ε (20)The random vector Ξ k satisfies a linear recurrence with initialcondition Ξ = whereas Ξ n , − k satisfies another recurrenceinitialised by Ξ n , = Ξ n , see for an explicit expression.We choose the random vectors ξ k and ξ − k (cid:48) independent for k , k (cid:48) > (cid:104) ξ k ξ Tk (cid:48) (cid:105) = I δ k , k (cid:48) (cid:104) ξ − k ξ T − k (cid:48) (cid:105) = I δ k , k (cid:48) (cid:104) ξ − k ξ Tk (cid:48) (cid:105) = (cid:104) ξ k ξ T − k (cid:48) (cid:105) = Ξ Rn = Ξ n , − n be the stochastic displacement with re-spect to the initial condition x after reversing the orbit and Σ Rn ( x ) the corresponding covariance matrix we define the re-versibility error E R ( x ) according to E Rn ( x ) = (cid:104) Ξ Rn · Ξ Rn (cid:105) = Tr (cid:0) Σ Rn ( x ) (cid:1) Σ Rn ( x ) = (cid:104) Ξ Rn ( Ξ Rn ) T (cid:105) (22)The covariance matrix that Σ Rn ( x ) can be expressed in termsof the tangent maps DM k ( x ) and their inverses DM − k ( x ) . Wedo not quote their expression, which can be found in sinceit can be proved that for a symplectic map M in R d the errorRE is related to LE by the following expression E Rn ( x ) = E L ( x ) + E Ln ( x ) + n − ∑ k = E Lk ( x ) E L = Tr ( I ) = d (23)or by the recurrence E Rn = E Rn − + E Ln + E Ln − initializedby E R = E L = d .We finally define the reversibility error due to to roundoff (REM). Letting M ε ( x ) the map computed with a finiteprecision ε (typically in the 8 bytes representation of reals ε (cid:39) − ) and with M − ε the inverse of M computed withfinite accuracy, so that M − ε (cid:0) M ε ( x ) (cid:1) (cid:54) = x , the modifiedreversibility error (REM) is defined according to E REM n ( x ) = (cid:13)(cid:13) M − n ε (cid:0) M n ε ( x ) (cid:1) − x (cid:13)(cid:13) ε (24)Letting x ε , n = M ε ( x ε , n − ) be the orbit computed with roundoff, the local error (cid:0) M ε ( x ε , n − ) − M ( x ε , n − ) (cid:1) / ε is similar toa random vector, possibly correlated, if the map has a suffi-ciently high computational complexity, see . The differ-ence with respect to the additive noise, we have considered todefine RE, is that we have just a single realization. As a conse-quence even though the behaviour of RE and REM is similar,the last one exhibits large fluctuations, when n varies. On theother side the computation of REM is really trivial since it re-quires a few lines of code, if the inverse map is available. Thisis the case of the reflection map for a waveguide. Indeed toreverse the evolution after n iterations we must simply changethe sign of the horizontal velocity components(s) and iterateagain the same map n times. C. Integrable and quasi-integrable maps
The computation of LE and RE can be analytically per-formed for an integrable map. If the 2D map is a rotation M ( x ) = R ( ω ) x then E Ln = E Rn = n . If the map M islinearly conjugated to a rotation by x = T X then the map innormal coordinates is R ( ω ) X and M ( x ) = TR ( ω ) T − x in thegiven coordinates. The loss of symmetry in the coordinates x induces periodic oscillations ( E Ln = A + ( − A ) cos ( ω n ) with A ≥ E Rn has a linear growth with oscillations. Amap M ( x ) is integrable if there is a coordinate transforma-tion x = T ( X ) such that the map becomes N ( X ) = R ( Ω ( J )) X where the frequency depends on the action J = (cid:107) X (cid:107) . Inthis case M ( x ) = T ◦ N ( Ω ( J )) ◦ T − ( x ) . The Lyapunov andreversibility errors LE and LE, depend on the initial conditionand in the normal coordinates X they grow according to E Ln ( X ) = + ( J Ω (cid:48) ( J ) n ) E Rn ( X ) n + (cid:0) J Ω (cid:48) ( J ) (cid:1) (cid:18) n + n (cid:19) (25)The Lyapunov error LE in the coordinates x exhibits a lineargrowth with oscillations, due to the loss of rotational symme- try as in the linear case and we have E Ln ( x ) = α ( n ) + n J Ω (cid:48) ( J ) α ( n )+ (cid:0) J Ω (cid:48) ( J ) n (cid:1) α ( n ) (26)where J = (cid:107) T − ( x ) (cid:107) and α k are periodic functions of n withfrequency Ω ( J ) . According to equation (23) the square ofthe reversibility error RE grows as n with oscillations whoseamplitude grows as n . If the map is quasi integrable as inthe neighborhood of an elliptic fixed point at x ∗ with non-resonant frequency, then it is conjugated with its normal formup to a remainder which can be made exponentially small with ∼ e − a / r α where r = (cid:107) x − x ∗ (cid:107) . Exponentially small estimatesfor the remainder can also be obtained for the error on the tan-gent map.If the system has an hyperbolic point at x = ( x , p ) T = andit is linear then M ( x ) = ( e − λ x , e λ p ) T and E Ln = ( n λ ) ∼ e n λ . If the system is non linear, in normal coordinates themap is X n + = e − λ ( XP ) X n and P n + = e λ ( XY ) P n with initialcondition X = ( X , P ) T . The Lyapunov error is E Ln ( X ) = e n λ (cid:16) + nXP λ (cid:48) ( XP ) (cid:17) (cid:16) + O (cid:0) e − λ n (cid:1) (cid:17) (27)The square of LE near an elliptic equilibrium point for ageneric symplectic map grows linearly with oscillations, thesquare of RE grows as n with oscillations whose amplitudegrows as n so that asymptotically they become negligible.Near a hyperbolic point LE and RE have the same exponen-tial growth. Finally REM behaves just as RE but exhibits largefluctuations. Right after the break up of the invariant curvessurrounding an elliptic point sticky orbits are found and theindicators have initially a power law growth followed by anexponential growth. We conclude that RE is the smoothest in-dicator, since LE can exhibit oscillations and REM is affectedby large fluctuations. However if we choose a fixed value of n large enough, the portraits all these indicators in a chosenplane of phase space are very similar. V. THE MODEL OF A 2D WAVEGUIDE
We present first a numerical analysis of the orbits for themap of the 2D waveguide where the corrugation is given by f ( x ) = cos ( x ) (28)The phase space coordinates are ( x , v x ) and the map is definedon the cylinder T ([ − π , π ]) × [ − , ] where T (([ − π , π ]) is theinterval [ − π , π ] with identified ends. The map has an ellipticfixed point at x = , v x = f ( x ) has a maximum and isapproximated by f ( x ) = − x + O ( x ) . We expand the mapretaining only the linear terms in x , v x . The first equation in(9) gives τ = + ε neglecting second order terms in x and v x .Since f (cid:48) ( x ) = − x + O ( x ) in the second equation the square of f (cid:48) are neglected and v zn is replaced with 1. Accordingly thesecond and third equations become v xn + = v xn − ε ( x n + τ v xn ) x n + = x n + τ ( v xn + v xn + ) (29)Letting M ( x ) = L x be the linear map we have L = − ε τ τ − ε τ − ε − ε τ (30)so that det L = Tr L = − ε τ . The map is conju-gated to a rotation L = VR ( ω ) V − where sin ( ω / ) = τε = ε ( + ε ) . In this case we have L n = L n . Close to x = z = + ε − ε x / ( x , v x ) correspond to caustics inconfiguration space, namely the plane ( x , z ) where the rayspropagate. Close to x = π the profile is approximated by z = − ε + ε ( x − π ) + O ( x − π ) so that τ = − ε andthe linearized map is ( x n + , v xn + ) T = L ( x n , v xn ) T , where L is given by (29) with ε → − ε . The map L is conjugated to ahyperbolic rotation R H ( α ) where sh ( α / ) = ετ = ε ( − ε ) .Near x = π the wave guide corresponds to an optical systemgiven by a plane and convex mirror .The map defined by (9) was computed by solving the equa-tion for τ with the bisection method initialized by τ = . / (cid:112) − v xn and τ = . / (cid:112) − v xn . The convergence isachieved also when v xn approaches 1 or −
1. The number of it-erations to reach machine accuracy varies between 40 and 60.Newton’s method is faster and machine accuracy is reached inless 10 iterations, but convergence problems are met when v xn approaches 1 or −
1. We have also computed the tangent mapand the explicit expression is written in Appendix A. Writ-ing the implicit equation G ( x n , v xn , τ ) =
0, the solution is notdefined if ∂ G / ∂ τ = ε = . , . N = n varies. Increasing N does not changethe plots significantly. Higher values of N are needed if onewishes to observe details in small regions where the transitionfrom ordered to chaotic orbits occurs. VI. THE 3D MODEL
For the 3D waveguide we choose the corrugated profile ac-cording to f ( x , y ) = cos x cos y =
12 cos ( x + y ) +
12 cos ( x − y ) (31)We remark first that for initial conditions y = v y = x = y and v x = v y theray propagates in the x = y plane as for 2D waveguide. Aftera rotation of π / x (cid:48) , z plane where x (cid:48) = √ x , v (cid:48) x = √ v x , the corrugation function depends onlyon x (cid:48) as f = (cid:2) cos ( √ x (cid:48) ) + (cid:3) and its period is √ π .When the ray does not propagate in a plane the simplest caseto analyze corresponds to the almost vertical propagation of the ray close to to a critical point ( x c , y c ) of the corruga-tion function, where grad f =
0. The fixed point of the map ( x c , y c , v x = , v y = ) is elliptic if f has a maximum at ( x c , y c ) and the corrugated surface near ( x c , y c ) behaves as a concavemirror.If f has a minimum at ( x c , y ) then the surface behaves as aconvex mirror and and the critical point of the map is hyper-bolic. For our model near x = y = f (cid:39) − ( x + y ) . The linearized map becomes v xn + = v xn − ε ( x n + τ v xn ) x n + = x n + τ ( v xn + v xn + ) v yn + = v yn − ε ( y n + τ v yn ) y n + = y n + τ ( v yn + v yn + ) (32)where τ = + ε . Since τ is constant in this approxima-tion the maps in the ( x , v x ) and ( y , v y ) phase planes decou-ple and each of them is area preserving. The linear frequen-cies are equal and given by sin ( ω / ) = ε τ as for the 2Dcase. The degeneracy can be removed choosing for instance f ( x , y ) = cos x + A cos y .Near x = , y = π we have f (cid:39) + ( x + ( y − π ) ) andthe fixed point of the map is hyperbolic. More generally ( , ) , ( ± π , ± π ) , ( ± π , ∓ π ) are elliptic and ( , ± π ) , ( ± π , ) are hyperbolic. We do not examine the projections of individ-ual orbits on 2D phase space planes nor on 3D hyperplanes.Close to the elliptic fixed points single 2D tori are recognis-able in the projections on 3D hyperplanes and even 2D phasespace planes, but the projection of several orbits does not pro-vide any useful information. For this reason we show here theplots of our short term indicators LE, RE and REM of orbitswhose initial points belong to 2D phase planes.The plots are obtained by computing our indicators for ini-tial conditions ( x , v x , y , v y ) in a family of 2D phase planeswhere y and v y / v x are kept fixed. Letting v x = v cos ( φ ) and v y = v sin ( φ ) we compare the errors for fixed valuesof y , φ letting x , v which vary in [ − π , π ] × [ − , ] , wherewe choose a regular grid of N g × N g points. We have fixedthe corrugation amplitude at ε = . y = φ = π / , π / , π /
4. Wehave computed LE by using the shadow orbit to evaluate thetangent map, namely we have replaced the partial derivatives ∂ M i / ∂ x j by finite differences, avoiding the cumbersome an-alytic evaluation. To this and we have chosen four differentinitial conditions x + e j δ for ≤ j ≤ e j are the or-thonormal base vectors in R namely ( e j ) k = δ j k . The tangentmap is approximated by DM n ( x ) e j (cid:39) w j ( n ) ≡ M n ( x + e j δ ) − M n ( x ) δ (33)so that the Lyapunov error becomes E Ln (cid:39) ∑ j = (cid:107) w j ( n ) (cid:107) (34)the discrepancy is of order δ and the choice δ = − wasmade, using double precision accuracy.For φ = x , z plane, but only if the sum in equation (34) runs only up to 2(corresponding to the phase space coordinates to x , v x ). Thisis equivalent to define LE by replacing DM n with its first 2 × y = v y =
0, the components of4D tangent map, not belonging the 2 × ( x , v x ) plane is invariant. To recover agree-ment between REM and RE, when LE is computed accordingto (34), which corresponds to the trace of ( DM n ) T DM n , onecan add a small random displacement of amplitude δ beforereversing the orbit, in order to bring the orbit out of the in-variant plane. When there is no invariant plane the agreementbetween REM and RE is recovered without any random kick.In figure (7) (a)-(c) we show the plot for φ (cid:39) π / y = VII. SHANNON CHANNEL CAPACITY
The possibility of computing fast stability indicators al-lows to establish interesting properties that depend on them.Among others, in information theory it is important to intro-duce the so called Shannon-Hartley channel capacity C = lim t → ∞ t log (cid:18) YX (cid:19) , (35)where the time-dependent relation between input ( X ) andoutput ( Y ) describes the time evolution of the transfer pro-cess taking the message from transmitter to receiver. Equa-tion (40) gives the maximum (ideal) information bit data-rate achievable in a physical system that carry electric sig-nals/electromagnetic waves. In the high-frequency asymp-totics, the propagation of waves through the corrugated chan-nel in Fig. 1 has an isomorphism with classical ray trajecto-ries underlying waves that bounce around the channel, givenby the symplectic map (9). Therefore, one can think that themaximum Lyapunov exponent has an isomorphism with theShannon-Hartley channel capacity. This has been recognisedin to be valid for any dynamical system, hence for a contin-uous time system. The relation between the channel capacityand the Lyapunov exponents of random matrices was first pro-posed by . C = λ = lim t → ∞ t log E L ( t ) = lim t → ∞ t log E R ( t ) (36) For recent works on the subject relating channel capacity toLyapunov exponents and entropy see and Consequently,fast dynamical indicators allow an estimation of the chan-nel capacity. For a map like the waveguide map depending onone parameter ε , given a finite number of iterations n we in-troduce the sequences C Ln ( x , ε ) = n log E Ln ( x , ε ) C Rn ( x , ε ) = n log E Rn ( x , ε ) (37)having the same limit C ( x , ε ) as n → ∞ . The channel capac-ities depend on the initial condition x , however if x belongsto an ergodic component the result is the same for almost anychoice of x . For a parallel wave guide ε = C Ln ( x , ε ) and C Ln ( x , ε ) vanish asymptotically accordingto C Ln ( x , ε ) = n log n + O ( n − ) C Rn ( x , ε ) = n log n + O ( n − ) (38)For a corrugated waveguide the phase plots of C Ln ( x , ε ) and C Rn ( x , ε ) for fixed n of are the same as the plot of E Ln and E Rn in a logarithmic scale shown in the figures 2 and 3 forthe 2D wave guide and 4 to 7 for the 4D waveguide, up to theconstant factor n − .The 2D ray reflection map for ε (cid:28) v x = ± √ ε cos (cid:16) x (cid:17) (39)see figure 8. For ε ∼ . ε . Wehave compared the dependence on n of the channel capaci-ties C Ln ( x , ε ) and C Rn ( x , ε ) for one regular and two chaoticorbits, see figure 9. For the first orbit the channel capacityvanishes as n − log n whereas for the last two chaotic orbits afinite limit is approached.In order to appreciate how the limit for n → ∞ is reached andthe dependence on the corrugated amplitude we compute thefollowing phase space averages C Ln ( ε ) = µ L ( E ) (cid:90) E C Ln ( x , ε ) d x C Rn ( ε ) = µ L ( E ) (cid:90) E C Rn ( x , ε ) d x (40)where E denotes the phase space and µ L ( E ) its volume.We have analyzed the dependence on the corrugation ampli-tude ε of the average channel capacities C Ln ( ε ) and C Rn ( ε ) for n =
100 and n = n =
200 the asymptotic value appears to be reached.We have found the following quadratic fit with ε C ( ε ) = . ε − . ε ≤ ε ≤
12 (41)A similar analysis can be performed for the 4D map by com-puting the phase space average with a Monte Carlo samplingrather than on a regular grid as for the 2D map.
VIII. CONCLUSIONS
We have analyzed the propagation of a ray on 2D and 3Dwaveguide, given by two parallel lines or planes, one of whichhas a periodic corrugation. In both cases the ray reflectionmaps on the uncorrugated line or plane are symplectic. Thestability properties are relevant for the long term propagationof the rays. Indeed in the chaotic regions a diffusion occurspreventing the coherent propagation of a signal. We have an-alyzed the reflection maps using the short time indicators re-cently proposed: the Lyapunov error LE, and the reversibilityerrors RE and REM. The square of LE is the trace of the co-variance matrix of the displacement, after n iterations of themap, induced by a small initial random displacement. Thesquare of RE is the trace of the covariance matrix of the dis-placement from the initial condition, when the map is iteratedforward and backward n times, adding at each step a smallrandom displacement. Replacing the random displacementwith the round off allows to define a modified reversibilityerror REM, which is similar to RE, though affected by largefluctuations due to the absence of averaging. In the 2D casea good qualitative agreement with the phase portrait is found.For a fixed number n of reflection of the rays on the uncorru-gated plane, all the indicators exhibit a similar behaviour witha power law growth for regular orbits and exponential growthfor chaotic orbits. For small corrugations the motion is regularalmost everywhere except in the neighborhood of the separa-trix, occurring when the rays hit a minimum of the corrugatedline or surface, and for rays propagating almost parallel to thewaveguideFor the orbits of the 4D map, describing the ray dynamics ofthe 3D waveguide, LE and RE provide a satisfactory stabilityportrait, just as REM. A stability analysis in this case cannotbe performed by looking at the projection of orbits on a 2Dphase plane. Altogether these indicators quantify separatelythe effect of small initial displacements, additive noise andround off. The limited computational load and the simplicityof implementation of these indicators make them a convenienttool also for a parametric study of 2D and 3D waveguides.We have also considered the relation of our indicators withthe channel capacity C . Indeed the limit of n − log E Ln and n − log E Rn for n → ∞ gives the maximum Lyapunov expo-nent which is isomorphic to the channel capacity. The de-pendence on n for a given initial condition is consistent withtheoretical estimates. The phase space average of the chan-nel capacity is found to rise with the corrugation amplitude ε following a quadratic law. IX. DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openlyavailable in https://github.com/gabrielegradoni/WaveguideStability . ACKNOWLEDGMENTS
GG wishes to acknowledge the financial support of TheRoyal Society under the grant INF | R2 | X. APPENDIX A: TANGENT MAP FOR 2D WAVEGUIDE
The equation which determines τ ( x , v ) is G ( x , v , τ ) = + ε f ( x + τ v ) − τ (cid:112) − v = ∂ τ∂ x = − ∂ G ∂ x ∂ G ∂τ = ε f (cid:48) √ − v − ε v f (cid:48) ∂ τ∂ v = − ∂ G ∂ v ∂ G ∂τ = ε τ f (cid:48) + τ v / √ − v √ − v − ε v f (cid:48) (42)Letting v x = cos θ ≡ v and v z = sin θ ≡ √ − v the recurrencewhich determines the symplectic map M = ( x , v ) is given by G ( x n , v n , τ ) = v n + = v n − ε f (cid:48) v n − ε f (cid:48) (cid:112) − v n + ε f (cid:48) x n + = x n + τ v n + τ (cid:112) − v n (cid:113) − v n + v n + (43)0The tangent map is given by ∂ v n + ∂ x n = ε f (cid:48)(cid:48) ( + ε f (cid:48) ) (cid:20) − ε v n f (cid:48) + (cid:113) − v n ( − ε f (cid:48) ) (cid:21) ×× (cid:18) + v n ∂ τ∂ x n (cid:19) ∂ v n + ∂ v n = − ε f (cid:48) + ε f (cid:48) (cid:32) ε f (cid:48) + v n (cid:112) − v n (cid:33) ++ ε f (cid:48)(cid:48) ( + ε f (cid:48) ) (cid:20) − ε v n f (cid:48) + (cid:113) − v n ( − ε f (cid:48) ) (cid:21) ×× (cid:18) τ + v n ∂ τ∂ v n (cid:19) ∂ x n + ∂ x n = + v n + v n + (cid:112) − v n (cid:113) − v n + ∂ τ∂ x n ++ τ (cid:112) − v n ( − v n + ) / ∂ v n + ∂ x n ∂ x n + ∂ v n = τ + v n + v n + (cid:112) − v n (cid:113) − v n + ∂ τ∂ v n ++ τ (cid:112) − v n ( − v n + ) / ∂ v n + ∂ v n −− τ v n (cid:112) − v n v n + (cid:113) − v n + (44)In all the previous formulae f = f ( x n + τ v n ) with f (cid:48) = f (cid:48) ( x n + τ v n ) and f (cid:48)(cid:48) = f (cid:48)(cid:48) ( x n + τ v n ) .In order to check the limit ε → τ (cid:113) − v n = v n + = v n x n + = x n + v n (cid:112) − v n (45)In this case the tangent map is given by ∂ x n + ∂ x n = ∂ x n + ∂ v n = ( − v n ) / ∂ v n + ∂ x n = ∂ v n + ∂ v n = f ( x ) periodic of period 2 π andscale the coordinates according to x = π x (cid:48) so that the new tangent map is given by ∂ x (cid:48) n + ∂ x (cid:48) n = ∂ x n + ∂ x n ∂ x (cid:48) n + ∂ v n = π ∂ x n + ∂ v n ∂ v n + ∂ x (cid:48) n = π ∂ v n + ∂ x n (47) XI. APPENDIX B. THE 3D WAVE-GUIDE
The equation for the corrugates waveguide is given by equa-tion (12) where f ( x , y ) is a periodic function. The normal vec-tor ν and the tangent vectors τ x , τ y are given by ν ( x , y ) = ( ε f x , ε f y , − ) T (cid:113) + ε ( f x + f y ) τ x = ( , , ε f x ) T (cid:112) + ε f x τ x = ( , , ε f y ) T (cid:113) + ε f y (48)where f x ≡ ∂ f / ∂ x and f y ≡ ∂ f / ∂ y . For an analogy with the2D case we might introduce the curvilinear abscissas s x , s y onthe lines on the corrugated plane having y fixed and x fixedrespectively s x ( x , y ) = (cid:90) x (cid:113) + ε f x ( x (cid:48) , y ) dx (cid:48) s y ( x , y ) = (cid:90) y (cid:113) + ε f y ( x , y (cid:48) ) dy (cid:48) (49)so that, keeping y fixed, ds x = (cid:112) + ε f x ( x , y ) dx and, keep-ing x fixed, ds y = (cid:112) + ε f x ( x , y ) dy . The transformationfrom ( x , y ) to ( s x , s y ) is invertible. Letting P = ( x , y , ) bea point on the plane and Q = ( x , y , z ( x , y )) a point on the cor-rugated plane and P = ( x , y , ) a new point on the plane weconsider a ray whose path is P , Q , P . We denote with v and v the velocity of the ray directed from P to Q and from Q to P . As for the 2 D case we keep P and P fixed letting Q varyand consider the ray path length H which depends on QH ( x , y ) = h ( x , y , x , y ) + h ( x , y , x , y ) h ( x , y , x , y ) = (cid:113) ( x − x ) + ( y − y ) + ( + ε f ( x , y )) (50)The velocities v and v , denoting for brevity h ≡ h ( x , y , x , y ) and h ≡ h ( x , y , x , y ) , are given by v = ( x − x , y − y , + ε f ( x , y )) T h ( x , y , x , y ) v = ( x − x , y − y , + ε f ( x , y ) ) T h ( x , y , x , y ) (51)1We first compute the derivatives of h and h with respect to x and y ∂∂ x h = x − x + ( + ε f ) ε f x h = v · τ x (cid:113) + ε f x ∂∂ y h = y − y + ( + ε f ) ε f y h = v · τ y (cid:113) + ε f y ∂∂ x h = x − x + ( + ε f ) ε f x h = − v · τ x (cid:113) + ε f x ∂∂ y h = y − y + ( + ε f ) ε f y h = − v · τ y (cid:113) + ε f y (52)Choosing h and h a functions of s x , s y rather than x , y onehas ∂ h / ∂ s α = v · τ α and and ∂ h / ∂ s α = − v · τ α for α = x , y . The function H is stationary when v · τ x = v · τ x and v · τ y = v · τ y namely when the projection on the tangentplane of the incoming and outgoing ray velocities are equal.This corresponds to the reflection condition since the normalcomponents of v and v are opposite. The derivatives withrespect to x , y and x , y are given by ∂ h ∂ x = x − xh = − v x ∂ h ∂ y = y − yh = − v y ∂ h ∂ x = x − xh = v x ∂ h ∂ y h = y − yh = − v y (53)Finally letting ( x ∗ , y ∗ ) the point where H is stationary andwhich depends on end points ( x , y ) and ( x , y ) we intro-duce the function F defined as the value of H evaluated at thestationary point and compute its differential F ( x , y , x , y ) = h ( x , y , x ∗ , y ∗ ) + h ( x , y , x ∗ , y ∗ ) dF = − v x dx − v y dy + v x dx + v y dy (54)As a consequence F is the generating function of acanonical transformation M which maps ( x , v x , y , v y ) into ( x , v x , y , v y ) . Denoting with x n = ( x n , v xn , y n , v yn ) T thephase space point reached after n iterations of the symplec-tic map M , the recurrence from x n to x n + is obtained firstby computing the transit time τ from P n = ( x n , y n ) to thepoint Q = ( x ∗ , y ∗ ) on the corrugated surface. Then we no-tice that the horizontal plane projections of the velocities v ∗ = v n − ν ( v n · ν ) and v n + of the ray reflected at Q and P n + = ( x n + , y n + ) are equal. Finally we determine the dis-placement from P n to P n + . The equations defining the recur-rence from x n to x n + are given by (14). (a)(b)(c)(d) FIG. 2: From the top: (a) Phase portrait of the map for the2D waveguide with corrugation z = + ε cos ( x ) and ε = . N =
200 iterations of the map. (c)Plot of RE for N =
200 iterations of the map. (d) Plot ofREM for N =
200 iterations of the map.2 (a)(b)(c)(d)
FIG. 3: From the top: (a) Phase portrait of the map for the2D waveguide with corrugation z = + ε cos ( x ) and ε = . N =
200 iterations of the map. (c)Plot of RE for N =
200 iterations of the map. (d) Plot ofREM for N =
200 iterations of the map. (a)(b)(c)
FIG. 4: From the top: (a) Plot of LE for the 4D reflectionmap of the 3D waveguide with corrugation amplitude ε = .
1. The initial conditions ( x , y , v x , v y ) where v x = v cos φ and v y = v sin φ , are chosen in a 2 D planeobtained by keeping fixed y = π / φ =
0. The error LEis computed for N =
200 iterations of the map. We let ( x , v ) vary on a regular grid of N g × N g points, with N g = [ − π , π ] × [ − , ] . (b)The same plot of RE. (c) The same plot of REM.3 (a)(b)(c) FIG. 5: Same plots as figure 4 for LE, RE, REM for ε = . y = π / φ =
0. The number ofiterations and grid points are the same N = , N g = (a)(b)(c) FIG. 6: Same plots as figure 4 for LE, RE, REM for ε = . y = π / φ =
0. The number ofiterations and grid points are the same N = , N g = (a)(b)(c) FIG. 7: Same plots as figure 3 for LE, RE, REM for ε = . y = φ = π / .
01. The number ofiterations and grid points are the same N = , N g = ε = . ε = . ε = . ε = .
25. The blue and red dots correspond to the initialconditions. Bottom panel: dependence on n of C Ln and C Rn forthe chaotic orbits whose initial points are the first two reddots in left panel (LE pink, RE orange for the first orbit, LEpurple, RE red for the second orbit) and the regular orbitwhose initial point is the last red point (LE blue, RE green) FIG. 10: Growth of the average channel capacities (on aregular grid of 400 points) with ε for two values of n : C Ln cyan for n = n =
200 and C Rn , orange for n = n = A. Bazzani, P. Freguglia, and G. Turchetti, Hamiltonian Analytical Opticsand Simulations of Betatronic Motion by Optical Devices, in
Nonlinear Dy-namics and Collective Effects in Particle Beam Physics , World Scientific,pp. 23-46, (2019). L.A.Bunimovich On the Ergodic Properties of Nowhere Dispersing Bil-liards,
Commun Math Phys.
65 (3) , 295-312, (1979). J. Redmond, S. Tabachnikov
Introducing symplectic billiards S. Tabachnikov billards S. Woo Park
An introduction to dynamical billiards https://math.uchicago.edu/ may/REU2014/REUPapers/Park.pdf D. Holm and G. Kovacic, Homoclinic chaos for ray optics in a fiber ,
Phys-ica D , p. 177 (1991). S. S. Abdullaev and G. M. Zavlaskii, Classical nonlinear dynamics andchaos of rays in problems of wave propagation in inhomogeneous media,
Usp. Fiz. Nauk. , p. 1 (1991). D. Douglas, Chaotic billiard lasers ,
Nature , p. 696 (2010). S. Creagh, Directional Emission from Weakly Eccentric Resonators,
Phys.Rev. Lett. , p. 153901 (2007). G. Tanner, Dynamical energy analysis Determining wave energy distribu-tions in vibro-acoustical structures in the high-frequency regime,
Journalof Sound and Vibration , p. 153901 (2007). E. Leonel, D. da Costa and C. Dettmann, Scaling invariance for the escapeof particles from a periodically corrugated waveguide,
Physics Letters A , 421 (2012). J. de Oliveira, C. Dettmann, D. , da Costa and E. Leonel, Scaling invari-ance of the diffusion coefficient in a family of two-dimensional Hamiltonianmappings,
Phys. Rev. E , p. 062904, (2013). G. Gradoni, J.-H. Yeh, B. Xiao, T. Antonsen, S. Anlage and O. E., Pre-dicting the statistics of wave transport through chaotic cavities by the ran-dom coupling model: A review and recent progress,
Wave Motion , 606(2014). G. Forte, F. Cecconi and A. Vulpiani, Transport and fluctuation-dissipation relations in asymptotic and preasymptotic diffusion acrosschannels with variable section,
Phys. Rev E , p. 062110, (2014).http://denali.phys.uniroma1.it/ ∼ cecconif/MyPapersPDF/gforte_PRE90.pdf F. Cecconi, V. Blakaj, G. Gradoni and A. Vulpiani, Diffusive transport inhighly corrugated channels,
Phys. Lett. A , , pp. 1084-1091, (2018). C. Froeschlé and E. Lega. On the Structure of Symplectic Mappings. TheFast Lyapunov Indicator: a Very Sensitive Tool.
Celestial Mechanics andDynamical Astronomy , , 167–195, 2000. C. Froeschlé, M. Guzzo and E. Lega. Graphical Evolution of the ArnoldWeb: From Order to Chaos.
Science , , 2108–2110, September 2000. C. Skokos. Alignment indices: a new, simple method for determining theordered or chaotic nature of orbits.
Journal of Physics A Mathematical Gen-eral , , 10029–10043, 2001. C. Skokos
The Lyapunov Characteristic Exponents and Their Computa-tion ,J. Souchay and R. Dvorak (eds.), Lecture Notes in Physics, BerlinSpringer Verlag Vol. 790, March 2010. R. Barrio. Theory and Applications of the Orthogonal Fast Lyapunov In-dicator (OFLI and OFLI2) Methods.
Chaos Detection and Predictability , , 55–92, March 2016. P. M. Cincotta and C. Simó. Simple tools to study global dynamics innon-axisymmetric galactic potentials - I.
A&AS , , 205–228, December2000. P. M. Cincotta, C. M. Giordano and C. Simó,
Phase space structure ofmulti-dimensional systems by means of the mean exponential growth factorof nearby orbits , Physica D Nonlinear Phenomena , 151 (August 2003). Z. Sandor, B .Erdi, A. Szell, B. Funk.
The relative Lyapunov indicator:an efficient method of chaos detection , Celestial Mechanics and DynamicalAstronomy , p. 127-138, (2004). https://arxiv.org/pdf/1205.0875.pdf Ch. Skokos, T. Bountis, Ch. Antonopoulos.
Geometrical properties of lo-cal dynamics in Hamiltonian systems: The Generalized Alignment Index(GALI) method , Physica D , p. 3054, (2007). Ch. Skokos, T. Manos.
The smaller (SALI) and the generalized (GALI)alignment indices: Efficient methods of chaos detection , Chaos De-tection and Predictability Springer Lecture Notes in Physics, Edi-tors: Ch. Skokos, Charalampos. G.A. Gottwald, J. Laskar (Eds.).https://arxiv.org/pdf/1412.7401.pdf C. Voglis, G. Contopoulos.
The relative Lyapunov indicator: an efficientmethod of chaos detection Phys. Rev. E
Vol , 372-377 (1998) . C. Voglis, G. Contopoulos, C. Efthymiopoulos. Detection of ordered andchaotic motion using the dynamical spectra Celestial Mechanics and Dy-namical Astronomy
Vol , 211-220 (1999) N. P. Maffione, L. A. Darriba · P. M. Cincotta · C. M. Giordano, compar-ison of different indicators of chaos based on the deviation vectors. Appli-cation to symplectic mappings
Celestial Mechanics and Dynamical Astron-omy (2011) https://arxiv.org/abs/1108.2196 N. P. Maffione, L. A. Darriba · P. M. Cincotta · C. M. Giordano,
Compara-tive study of variational chaos indicators and ODEs numerical integrators , International Journal of Bifurcation and Chaos
World Scientific PublishingCompany (2012) , 1230033 (2012) D. Ruelle,
An inequality for the entropy of differential maps , Bol. Soc. Bras.Mat. , , 83-87, (1978). J. Laskar, C. Froeschlé and A. Celletti. The measure of chaos by the nu-merical analysis of the fundamental frequencies. application to the standardmapping.
Physica D: Nonlinear Phenomena , , 253–269, 05 1992. L. Casetti, C. Clementi, and M. Pettini,
Riemannian theory of Hamiltonianchaos and Lyapunov exponents , Phys. Rev.
E 54 , 5969-5984 (1996). M. Cerruti-Sola, G. Ciraolo, M. Franzosi, M. Pettini.
Riemannian geometryof Hamiltonian chaos: hints for a general theory.
Phys. Rev. E Stat. Nonlin.Soft Matter Phys. , 046205 (2008) G. Gottwald, I. Melbourne,
Testing for chaos in deterministic systems withnoise , Physica D: Nonlinear Phenomena , 100-110 (2005) F. Panichi, L. Ciotti and G. Turchetti,
Fidelity and reversibility in the re-stricted three body problem , Communications in Nonlinear Science andNumerical Simulation , 53 (2015). F. Panichi, K. Go´zdziewski and G. Turchetti,
The reversibility error method(REM): a new, dynamical fast indicator for planetary dynamics , MNRAS , 469 (June 2017). F. Panichi and G. Turchetti, Fast indicators of orbital stability: a survey onLyapunov and Reversibility errors
INTECH G. Turchetti, F. Panichi and K. Gozdziewski A complete set of fast Lya-punov and Reversibility invariant indicators. To be submitted D. Faranda, M. F. Mestre and G. Turchetti, Analysis of Round off Errorswith Reversibility Test as a Dynamical Indicator,
International Journal ofBifurcation and Chaos , , p. 1250215 (September 2012). G. Turchetti, S. Vaienti and F. Zanlungo. Asymptotic distribution ofglobal errors in the numerical computations of dynamical systems.
PhysicaA Statistical Mechanics and its Applications , , 4994–5006, (November2010). E. Hairer, C. Lubich and G. Wanner.
Geometric Numerical Integration:Structure-Preserving Algorithms for Ordinary Differential Equations; 2nded.
Springer, Dordrecht, (2006). G. Turchetti, F. Panichi Birkhoff normal forms and stability indicatorsfor betatronic motion
Proceedings of the NOCE workshop, Arcidosso sept.2017 , S. Di Mitri, S. Chattopadhyay, M. Cornacchia Editors, World Scien-tific (2019) Oseledets, V. I. (1968).
A multiplicative ergodic theorem. characteristicLjapunov, exponents of dynamical systems.
Trudy Moskovskogo Matem-aticheskogo Obshchestva, 19:179-210. G. Friedland and A. Metere (eds.),
Isomorphism between Max-imum Lyapunov Exponent and Shannon?s Channel Capacity ,https://arxiv.org/pdf/1706.08638.pdf, January 2018. Tim Holliday, Andrea Goldsmith, Peter W Glynn
Capacity of Finite StateChannels Based on Lyapunov Exponents of Random Matrices
September2006IEEE Transactions on Information Theory 52(8):3509 - 3532 DOI:10.1109/TIT.2006.878230 (2006) Tim Holliday, Peter Glynn, Andrea Goldsmith
On En-tropy and Lyapunov Exponents for Finite-State Channels J. M. Ebeid