Chaotic properties of the truncated elliptical billiards

Chaotic properties of the truncated elliptical billiards

Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323 Contents lists available at ScienceDirect Commun Nonlinear Sci Numer Simulat journal homepage: ...

2MB Sizes 5 Downloads 110 Views

Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Chaotic properties of the truncated elliptical billiards V. Lopac a,*, A. Šimic´ b a b

Department of Physics, Faculty of Chemical Engineering and Technology, University of Zagreb, Croatia Department of Physics and Applied Mathematics, Faculty of Science, University of Navarra, Pamplona, Spain

a r t i c l e

i n f o

Article history: Received 28 January 2010 Accepted 16 March 2010 Available online 19 March 2010 Keywords: Chaotic billiards Truncated elliptical billiard Elliptical stadium billiard Poincaré sections Box-counting method Orbit stability Chaotic fraction Resonant cavities Optical resonator Hamiltonian system

a b s t r a c t Chaotic properties of symmetrical two-dimensional stadium-like billiards with piecewise flat and elliptical segments are studied numerically and analytically. Their sensitivity to small variations of the shape parameters can be usefully applied for optimal construction of the dielectric and polymer optical microresonators. For the two-parameter truncated elliptical billiards (TEB) the existence and linear stability of several periodic orbits are investigated in the full parameter space. Poincaré plots are computed and used for evaluation of the chaotic fraction of the phase space by means of the box-counting method. A highly chaotic behavior prevails in the region of elongated elliptical arcs, where most of the existing orbits are either neutral or unstable. In the parameter space, the upper limit of the fully chaotic behavior isffi reached when the relation between the two shape parameters becomes pffiffiffiffiffiffiffiffiffiffiffiffiffi d ¼ 1  c2 , corresponding to truncated circles. Above this limit, for flattened elliptical arcs, mixed dynamics with numerous stable elliptic islands is found. In both regions parabolic orbits are present, many of them identical to orbits within an ellipse. These properties of the TEB differ remarkably from the behavior in the elliptical stadium billiards (ESB), where the chaotic region in the parameter space was strictly bounded from both sides. In order to follow the transition from TEB to ESB, a generalisation to a novel three-parameter family (GTEB) of stadium-like boundary shapes with elliptical arcs is proposed. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Two-dimensional planar billiards are nonlinear Hamiltonian systems with rich and interesting dynamical properties. In the classical mechanics, a point particle, moving with constant velocity within a closed boundary and exhibiting specular reflections on the walls, can have regular, mixed or fully chaotic dynamics. This behavior depends strongly on the shape of the billiard boundary, and the resulting dynamical properties are extremely sensitive to small variations of the shape parameters. A detailed exposition of the fundamental properties of classical chaotic billiards, including a list of relevant references, is given in the book by Chernov and Markarian [1]. With exception of some special classes of polygonal billiards, the only known boundary shape with integrable classical dynamics is ellipse, with the circular billiard included [2]. Complete classical solutions for the family of confocal elliptical billiards, including also the full elliptical shapes, have been explored in [3]. The fully chaotic billiards are the Sinai billiard [4], the Bunimovich stadium billiard [5] and the cardioid (Robnik) billiard [6–8]. The Sinai billiard is dispersing, predictably leading to chaotic dynamics. However, the stadium billiard, with a rectangle inserted between the two half-circles, is also chaotic. It has been found that the defocusing mechanism responsible for chaos in the circular stadium is present also in a number of other shapes and different billiard shapes with fully chaotic dynamics have been constructed [9–12]. Interpretations of mechanisms responsible for creating chaoticity in billiards, explaining how and why the absolutely focusing * Corresponding author. Tel.: +385 1 45 97 106; fax: +385 1 45 97 135. E-mail address: [email protected] (V. Lopac). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.03.007

310

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

components may generate hyperbolicity in billiards, containing the most recent results, are reviewed in [13]. Much attention has been paid to numerical investigations of billiards with mixed dynamics, especially of those with boundaries consisting of circular arcs, but also with elliptical, hyperbolic and parabolic segments [2,14–20]. A two-dimensional planar billiard can also be approached as a quantum system, since the corresponding Schrödinger equation is mathematically equivalent to the Helmholtz equation for the vibrations of a two-dimensional membrane. The quantal energy spectra and wave functions, reflecting the chaotic properties of the corresponding classical billiards, are equally sensitive to the shape variations [21–26]. Importance of the billiard shapes has been stressed also in the recent investigations of the graphene quantum dots [27]. All these properties place billiards among the most attractive types of dynamical models in the study of ergodic and mixing properties of Hamiltonian systems [28,2]. The recent revival of general interest for billiards is mostly due to the fact that there are many physical systems whose study can be reduced to a billiard-type problem. The problem of stationary electromagnetic waves within the waveguides and optical cavities in the limit of small wavelengths can be approximated by the ray propagation mechanism, which is mathematically equivalent to the problem of a point particle moving freely in a billiard. In applications to laser technology, the consideration of chaotic dynamics was crucial for increasing the output power of the dielectrical, semiconducting and polymer microlasers [29]. In experiments with the asymmetric resonant cavities, the critical element in technological improvements has been provided either by distorting slightly the circular shape of the optical resonators [30,31] or by considering the long and narrow stadium and quasi-stadium shapes [32–34]. The fine adjustments of the resonator shapes described above can be achieved by using the two-parameter planar boundary with two symmetrically placed elliptical arcs at the opposite sides of a rectangle. In our previous work we investigated the two-parameter elliptical stadium billiard (ESB), first explored by Donnay [10], which is also a special case of the mushroom billiard [11,35–37]. The upper limit of the fully chaotic region in the parameter space for ESB was found by Wojtkovski [9]. Our numerical investigation in the full parameter space has shown that the ESB are fully chaotic for a bounded set of shape parameters and have mixed behavior in other parameter regions [17,38]. This confirmed the results of analysis [39–42] identifying the emergence of the stable pantographic orbits as the mechanism generating the lower limit of chaos in the elliptical stadium billiards. In this work, we extend our investigations to a larger class of stadium-like billiards. We first examine the dynamical properties of the truncated elliptical billiards (TEB), another family of billiards with elliptical arcs joined by flat segments. We follow the procedure used in [17] for elliptical stadia (ESB) and compare the results for these two billiard families, as a possible basis for estimating the directional output power from the corresponding dielectrical cavities, in dependence on the resonator shapes [43]. The TEB is defined as a point particle in the two-parameter planar domain, constructed by truncating an ellipse on opposite sides (Fig. 1). A symmetrical stadium-like shape thus obtained consists of a rectangle with two elliptical arcs added at its opposite ends. The corresponding billiard has been analysed by Del Magno [44] who, investigating a restricted part of the parameter space and applying the method of invariant cones, determined the region of hyperbolic behavior and presented an estimate of the region where such billiard could be ergodic. Although chaotic dynamics of the truncated elliptical billiard was not analysed elsewhere, the shape itself has been introduced in different contexts: as a particular cross-section of the three-dimensional cavity, for example in description of the barrel billiard introduced by Bunimovich [45,46], or as a geometrical model for the mathematical description of the human heart in [47]. The chaos in billiards obtained by various ways of truncating the circle has been considered in [48–50]. In our approach, the two billiard families, TEB and ESB, have in common the two shape parameters d and c, which determine half of the width and the height, respectively, of the central rectangle. They both include as limits the circular, elliptical, rectangular and quadratic shapes. However, whereas the ESB comprises the Bunimovich (circular) stadium, this shape does not exist within the TEB. For computational reasons, the horizontal diameters of these billiards were normalised to 2, in the same way as in our previous work on symmetrical lemon-shaped billiards with the parabolic, generalised power-law, hyperbolic and elliptical arcs [18–20]. In this paper, we investigate numerically and analytically the TEB in the full parameter space, described by two parameters d and c. In Section 2, we define the TEB boundary and describe its geometrical properties. In Section 3, the existence and stability conditions of selected periodic orbits are discussed and illustrated by Poincaré plots and orbit diagrams. In Section 4, the Poincaré sections are used to compute, with the box-counting method, the chaotic fraction of the phase space, for a large set of shape parameters. The region in which the billiard seems to be fully chaotic is identified, and for the remaining part of the parameter space the numerical estimates of the chaotic fraction, in dependence on the shape, are presented. The results are shown in the parameter-space diagram and compared with the same type of diagram for the elliptical stadium billiard. In Section 5, we propose a possible generalisation (GTEB) of the truncated elliptical billiard, providing a transition between the two types of the stadium-like elliptical billiards ESB and TEB. In Section 6, we summarise the obtained results and propose further investigations. At the end, in a brief Appendix A, we discuss the box-counting method and analyse the dependence of results on relevant parameters.

2. Geometrical properties of the truncated elliptical billiard In our parametrisation the TEB is defined in the ðx; yÞ plane by means of the two shape parameters d and c, satisfying conditions 0 6 d 6 1 and 0 < c < 1. The billiard boundary is described as

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

311

 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 Fig. 1. Three types of the  TEBpshape: ffiffiffiffiffiffiffiffiffiffiffiffiffiffi(a) circular shape, with d ¼ 1  c ; (b) the shape with flattened elliptical arcs d > 1  c ; (c) the shape with elongated elliptical arcs d < 1  c2 ; (d) notation for an impact point T on the billiard boundary, the point S of intersection of the incoming path with the x-axis, and the corresponding angles a, b, /, /0 and h.

yðxÞ ¼

8 c; > > <

if 0 6 jxj < d

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi > c ffi > 1  x2 ; if d 6 jxj 6 1 :  pffiffiffiffiffiffiffi 2

ð1Þ

1d

The horizontal diameter is normalised to 2, so that the horizontal semiaxis of the ellipse is 1. The possible height 2c of the billiard extends from 0 to 1. The horizontal length of the central rectangle is 2d and the vertical semiaxis of the ellipse is .pffiffiffiffiffiffiffiffiffiffiffiffiffi ffi c 1  d2 . In special cases, for d ¼ 0 thepshape ffiffiffiffiffiffiffiffiffiffiffiffiffiis ffi a full ellipse, for d ¼ 0 and c ¼ 1 a full circle. For d ¼ 1 the shape is rectangular and for d ¼ c ¼ 1 a square. which For d ffi¼ 1  c2 one obtains a set of truncated circle billiards,p pffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiseparates two distinct billiard classes, one for d < 1  c2 with elongated elliptical arcs and the other with d > 1  c2 and flattened elliptical arcs. Fig. 1(a)–(c) shows three typical shapes of the truncated elliptical billiard, with circular, flattened and elongated elliptical arcs, respectively. Fig. 1(d) shows the notation concerning the impact point on the boundary and the angles corresponding to the incoming and outgoing paths. pffiffiffiffiffiffiffiffiffiffiffiffiffiffi The coordinates of the focal points are, for d < 1  c2 ,

2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 1  d2  c2 5 F 4 ;0 1  d2 and, for d >

ð2Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  c2 ,

2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3 c2 þ d2  15 F 40;  1  d2

ð3Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 Theypcontain ffiffiffiffiffiffiffiffiffiffiffiffiffiffi the important term s ¼ c þ d  1 which is negative for d < 1  c , positive for d > 1  c , and zero for 2 d ¼ 1  c (circular arcs). In Fig. 2, presenting the structure of the ðc; dÞ parameter space, the circular limit is shown as the thick circular line denoted by letter (a).

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

312

Fig. 2. Diagram of the two-dimensional parameter space (c, d). The letters denoting the lines are explained in the text. The thick black line corresponds to the truncated circles, the limit at which the stable diametral orbits emerge. The parabolic bow-tie orbit exists everywhere within the blue shaded region. The hour-glass orbit exists and is parabolic everywhere within the red shaded region. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)

For d < jxj 6 1 the curvature radius is

 R¼

  3=2 1  d2  c2 1  x2 þ c2   c 1  d2

ð4Þ

For 0 6 jxj < d the boundary is flat and the curvature radius is R ¼ 1, but for jxj ¼ d has a discontinuity and drops to h i3=2   2  c 1  d2 . At the endpoints of the horizontal axis of the ellipse ðjxj ¼ 1Þ the curvature radius Rd ¼ 1  d2 þ c2 d2

  is R1 ¼ c2 1  d2 , which reduces to R1 ¼ 1 for circular arcs. For full ellipses with d ¼ 0, at x ¼ 0 the curvature radius is R0 ¼ 1=c. In the same way as in [17], the symbols h, / and /0 , respectively, denote the angles which the normal to the boundary, the incoming path and the outcoming path close with the x-axis (Fig. 1(d)). The angle between the incoming (or outcoming) path and the normal to the boundary is b ¼ ð/0  /Þ=2. The angle between the tangent to the boundary at the point Tðx; yÞ of impact and the incoming (or outcoming) path is a ¼ ðp=2Þ  b. The angles h, / and /0 satisfy the relation ð2 tan hÞ=ð1  tan2 hÞ ¼ ðtan / þ tan /0 Þ=ð1  tan / tan /0 Þ. Denoting the slope of the incoming path as a ¼ tan / and taking into account that the slope of the normal to the boundary at the point Tðx; yÞ is given by C ¼ tan h ¼ 1=y0 ðx; yÞ, one obtains the slope of the outcoming path a0 = tan /0 , which can be written as

a0 ¼

2C  að1  C2 Þ

ð5Þ

2aC þ 1  C2

After inserting the dependence of a and a0 , deduced from the geometrical properties of the orbits, on impact coordinates, the expression (5) becomes the basis for finding the existence criteria for particular periodic orbits [17]. Since in our further analysis we are dealing only with orbits which are symmetrical with respect to both x- and y-axis, we refer to the impact points Tðx; yÞ in the first quadrant ðx P 0 and y P 0Þ, with no loss of generality for the obtained results. The points PðX; V x Þ in the Poincaré sections are obtained by plotting the x-component of the velocity V x ¼ cos / versus X, the x-coordinate of the intersection point S of the incoming path with the x-axis (Fig. 1(d)), as explained in [17–20]. The Poincaré diagrams obtained in this way are area preserving. It should be stressed that the number of impact points Tðx; yÞ is equal to the number of iterations for a special orbit, but the corresponding number of points PðX; V x Þ in the Poincaré section is usually smaller, since many of the orbit segments do not cross the x-axis. The linear stability of a periodic orbit is assured if the absolute value of the trace of the stability matrix M is smaller than 2, thus if

2 < Tr M < 2

ð6Þ

The orbits satisfying Eq. (6) are elliptic, those with jTr Mj ¼ 2 are neutral (parabolic) and those with jTr Mj > 2 hyperbolic [2,17]. The stability matrix of the closed orbit of period N can be written as M ¼ M 12 M 23  M N1 , where the 2  2 matrix Mik for two subsequent impact points T i and T k , connected by a rectilinear chord of the length qik , is

0 B M ik ¼ @

qik sin ai  sin a þ Ri sin a k

qik

 Ri R þ k

k

sin ak Ri

þ

sin ai Rk

 sin aqi iksin a ak  sin sin ai

þR

k

1 k

qik

sin ai

C A:

ð7Þ

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

313

3. Classical dynamics of the truncated elliptical billiard Both in ESB and in TEB, the full parameter space ð0 6 d 6 1 and 0 < c < 1Þ is strictly divided into two parts by a line denoting the emergence of the stable horizontal two-bounce diametral orbit. For TEB this limit is given by the equation

c2 þ d2 ¼ 1

ð8Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi The corresponding curve d ¼ 1  c2 is shown in Fig. 2 denoted by letter (a). It divides the parameter space into two parts. A similar division exists in the ESB, and in [17] we discussed separately the orbits within each of the two regions. For the TEB another distinction is much more important, since the periodic orbits in this billiard fall into two categories. First are the orbits with some of the impact points on the flat parts of the boundary, whereas the second group consists of orbits which have impact points only on the curved (elliptical) parts of the boundary. Those orbits are identical to the orbits within the ellipse, which is integrable and in which, with the exception of the two diametral orbits along the two main axis, all remaining orbits are parabolic. This significantly affects the properties of orbits in the truncated elliptical billiards. In this paper, we analyse numerically the dynamics of TEB in the full parameter space. We also examine analytically the existence and stability conditions for several typical orbits, similarly as in the description of the ESB [17]. In Fig. 3(a)–(d), we show ffi pffiffiffiffiffiffiffiffiffiffiffiffiffi thepffiffiffiffiffiffiffiffiffiffiffiffiffi Poincaré sections for d ¼ 0:19 and different values of c < 1  d2 ¼ 0:982. Similar results for d ¼ 0:60 and c ffi 6 1  d2 ¼ 0:80 are shown in Fig. 4(a)–(d). The elliptical arcs for these shapes have an elongated form. These pictures reveal a highly chaotic behavior. There are no visible fixed points or elliptic islands. However, flights of points typical for neutral orbits can be discerned. This is remarkably different from the results for the elliptical stadium billiards [17], where in the corresponding parameter region there were many fixed points and elliptic islands due to stable pantographic and some other types of orbits. The Poincaré sections for the part of the parameter space with flattened elliptical arcs, which has not been investigated previously, are shown for d ¼ 0:19 and several values of c > 0:982 in Fig. 3(e)–(h), and for values d ¼ 0:60 and c > 0:80 in Fig. 4(e)–(h). In this parameter region, many fixed points and stable islands are noticed, especially the large elliptic island corresponding to the horizontal two-bounce orbit, but also here the traces of the parabolic orbits can be discerned. Several typical symmetrical low-period orbits are shown in Fig. 5, together with the corresponding fixed points in the Poincaré diagrams. 3.1. Orbits with impact points on the flat parts of the boundary 3.1.1. Diamond orbit The diamond orbit of period four, shown in Fig. 5(e), exists for any parameter choice. It has two bouncing points at the ends of the horizontal semiaxis, and the other two on the flat parts on the boundary. To assess its linear stability, one should calculate the stability matrix M ffi¼ ðM 01 M 10 Þ2 . The angles contained in the matrix are given as sin a0 ¼ c=q and sin a1 ¼ 1=q, pffiffiffiffiffiffiffiffiffiffiffiffiffi where the chord is q ¼ 1 þ c2 . The curvature radius at x ¼ 1 is given by Eq. (4). This leads to the trace

Fig. 3. Poincaré plots, showing the coordinate 1 < x < 1 and the slope 1 < V x < 1 at the intersections points of the incoming path with the x-axis, for ffi pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi ffi d ¼ 0:19 and various c. For d < 1  c2 the points densely fill the phase space. The large bands, visible for d > 1  c2 , correspond to the degenerated diametral horizontal bouncing mode.

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

314

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fig. 4. Poincaré plots for d ¼ 0:60 and various c. For d < 1  c2 the points densely fill the phase space. For d > 1  c2 the large bands corresponding to the degenerated horizontal bouncing modes are visible. With increasing c the stable diamond orbit and multidiamond orbits sequentially emerge.

Fig. 5. Typical lowest order periodic orbits in TEB: (a) bow-tie orbit; (b) rectangular orbit; (c) horizontal two-bounce orbit; (d) tilted two-bounce orbit; (e) diamond orbit; (f) multidiamond orbit with n ¼ 2; (g) hour-glass orbit; and (h) the 8-shaped orbit. Below each orbit the corresponding fixed points in the Poincaré diagrams are shown.

" # 2 2q2 1 1 R

Tr M ¼ 2 2

and to the condition q2 < R or 1 þ c2 < c2

1 d > pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ c2

ð9Þ

  1  d2 . Thus the stable diamond orbit appears if

ð10Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   or c > 1=d2  1. This limit is shown in Fig. 2 as the line denoted by letter (b). The only exception is for pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d ¼ ð2 þ c2 Þ=ð2 þ 2c2 Þ, when R ¼ 2q2 and Tr M ¼ 2, thus the diamond orbit is neutral.

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

315

3.1.2. Multidiamond orbits The multidiamond orbit of order n is the orbit of period 2 þ 2n. It has two bouncing points at the ends of the horizontal axis and 2n bouncing points on the flat parts of the boundary (see Figs. 4(h) and 5(f)). Such orbit exists if d > 1  ð1=nÞ. As explained for p a similar case ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi in [17], the chord q in Eq. (9) should be replaced by L ¼ nqn , where, for the truncated elliptical billiard, qn ¼ ð1=n2 Þ þ c2 . The trace of the stability matrix is then

" # 2 2q2n n2 Tr M ¼ 2 2 1 1 R

ð11Þ

with R given by Eq. (4). The resulting condition for the stability of the multidiamond orbit is

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d>

1

c2 1 þ c2 n2

ð12Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     ffi or c > 1  d2 1  n2 1  d2 . The limiting curves (12) are plotted in Fig. 2. As already explained, the line with n ¼ 1 (diamond orbit) is denoted by letter (b), and the lines with n ¼ 2; n ¼ 3 and n ¼ 4 are denoted by letters (c), (d) and (e), respectively. Above them lie the lines for n > 4. For c ! 1 the minimal values of d above which the multidiamond orbits appear are

lim d ¼

c!1

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2 n

ð13Þ

The emergence of multidiamond orbits can be followed by observing the Poincaré sections for a set of shapes with d ¼ c (Fig. 6). The values of this parameter, for which an orbit of new n appears, are given as intersections of the straight line d ¼ c with curves (12) and obey the equation n2 d4  ðn2  2Þd2  1 ¼ 0. For the diamond orbit ðn ¼ 1Þ this equation reads r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . pffiffiffi 5  1 2 ¼ 0:78615. For the same type of boundary, the stable d4 þ d2  1 ¼ 0 and the orbit appears for d ¼ c ¼ .pffiffiffi two-bounce orbit, dividing elongated from flattened shapes, appears at d ¼ c ¼ 1 2 ¼ 0:707107. Besides the diamond and multidiamond orbits, in Figs. 4 and 6 one discerns the presence of the 8-shaped stable orbit, shown in Fig. 5(h). 3.2. Orbits with impact points only on the elliptical parts of the boundary 3.2.1. Horizontal diametral two-bounce orbits The horizontal two-bounce p orbit (Fig. obviously ffiffiffiffiffiffiffiffiffiffiffiffiffiffi exists for all combinations of d and c. The stability condition [2,17] ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 5(c)) p 2 or c > q=ð2RÞ < 1 takes the form d > 1  c 1  d2 , so that bifurcations giving birth to stable diametral orbits appear at pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 the values d ¼ 1  c , corresponding to circular arcs. The corresponding line in the parameter space diagram in Fig. 2 is denoted by letter (a). In the Poincaré diagrams, the corresponding stable islands are visible as two large bands near V x ¼ 1.

.pffiffiffi .pffiffiffi Fig. 6. Poincaré plots for a set of different values d ¼ c. For d ¼ c < 1 2 the points densely fill the phase space. For d ¼ c > 1 2 the narrow bands corresponding to the horizontal bouncing modes are visible. Besides, at d ¼ c ¼ 0:85 four islands appear corresponding to the diamond orbit. For d ¼ c ¼ 0:90 new six narrow islands appear due to the multidiamond orbit with n ¼ 2. Two symmetrical triangle-shaped islands for d ¼ c ¼ 0:90 originate from the ‘‘8-shaped” orbit. For higher values of c successive multidiamond orbits of higher order appear.

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

316

3.2.2. Tilted diametral two-bounce orbits For a tilted two-bounce orbit (Fig. 5(d)) a ¼ a0 and a ¼ C, so that, according to Eq. (5), it bounces at the point Tðx; yÞ on the billiard boundary with derivative y0 if yy0 þ x ¼ 0. This is realised for any d < x < 1, provided that c2 þ d2 ¼ 1, thus only for the truncated circle. Since in this case the chord is q ¼ 2 and the radius is R ¼ 1; q=ð2RÞ ¼ 1, and these orbits are neutral. 3.2.3. Rectangular orbit Further we investigate properties of the rectangular orbit shown in Fig. 5(b). For this orbit a ¼ 0; a0 ¼ 1 and C ¼ 1, so that according to Eq. (5), this orbit exists if the derivative on the boundary is y0 ¼ 1. Corresponding solutions for the impact qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi    1  d2 þ c2 and y ¼ c2 1  d2 1  d2 þ c2 . The condition d < x < 1 leads to the requirement points are x ¼ 1  d2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2 þ 4  c d< 2

ð14Þ

 

or c < 1  d2 d. This limit is shown in Fig. 2, denoted by letter (f). If we denote the points with positive x by 1 and the points on the negative side by 1, the deviation matrix can be computed as

M ¼ ðM 11 M11 Þ2

ð15Þ

.pffiffiffi using the matrix (7) with a = a0 and R = R0 , where the angle by sin a ¼ 1 2. The chords are h pffiffiffi a is given     

 3=2 q  q1;1 ¼ 2x and q0  q1;1 ¼ 2y, and the curvature radius is R ¼ 2 2c2 1  d2 . If we define 1  d2 1  d2 þ c2



q

q0

R sin a R sin a



q R sin a

þ

q0



R sin a

ð16Þ

the trace of the deviation matrix becomes

h i Tr M ¼ 2 2ð2U þ 1Þ2  1

ð17Þ

For the rectangular orbit one obtains U ¼ 0 and Tr M ¼ 2, leading to the conclusion that this orbit is neutral for all allowed shapes, both flattened and elongated. 3.2.4. The bow-tie orbit We investigate the existence and stability of the bow-tie orbit (the lowest pantographic orbit), shown in Fig. 5(a). This orbit exists if the coordinates x and y of the impact point and the derivative y0 of the boundary at this point satisfy Eq. (5), which now reads 2yy0 þ x  xy02 ¼ 0, giving as solution the coordinates of the point of impact qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi    1  d2 1  d2  c2 , possible only if d < 1  c2 . The expression 1  d2  c2 and y ¼ c2 x¼ 1  d2  2c2 d < x < 1, requiring that this point should lie on the elliptical part of the boundary, leads to the existence condition

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  c2  c4 þ 4c2 d< 2

ð18Þ

 .pffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi or c < 1  d2 2  d . This upper limit for the existence of the bow-tie orbit is shown in Fig. 2 and is denoted with the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi      1  d2 2 1  d2  c2 . The chords letter (g). The corresponding angle a needed in the matrix (7) is given by sin a ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi      ffi are q  q1;1 ¼ 2y ¼ 2c2 1  d2 . The curvature ra1  d2 1  d2  c2 and q0  q1;1 ¼ 2 x2 þ y2 ¼ 2 1  d2  c2  pffiffiffi   dius at the impact point is obtained by substituting the impact coordinate x into (4) and reads R ¼ 2c2 2 1  d2 . Denoting the impact points with positive x by 1 and the points where it is negative by 1, the deviation matrix is given by Eq. (15). If we define U as in Eq. (16), the trace of the deviation matrix is given again by Eq. (17). The left-hand side of the stability condition (6) is valid automatically, but the right-hand side is fulfilled only if 1 < U < 0. By using the given expressions for sin a, R, q and q0 and substituting them into (16), one obtains U ¼ 1 and Tr M ¼ 2 for all allowed cases. The conclusion is that the bow-tie orbit is neutral for all parameter values below the existence limit given in Eq. (18). 3.2.5. The hour-glass orbit The hour-glass orbit (Fig. 5(g)) looks like the bow-tie orbit rotated by p=2. It exists if the coordinates x and y of the impact point and the derivative y0 of the boundary at this point satisfy the equation 2xy0 þ yy02  y ¼ 0, giving as solution the coorqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi     2     d þ c2  1 and y ¼ c 2d2  2 þ c2 1  d2 d2 þ c2  1 . The condition dinates of the impact point x ¼ 1  d2 d < x < 1 that this point should lie on the elliptical part of the boundary leads to the requirement

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2 þ c4 þ 4 1
317

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi

c2

ð19Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi  ffi or 2 1  d2 < c < 1  d4 d. These limits are shown in Fig. 2 and are denoted by letters (h) and (i). If we denote the points with positive y by 1 and the points on the negative side by 1, the deviation matrix can again be qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

   calculated from Eq. (15). The angle a needed in the calculation is given as sin a ¼ c2 2 d2 þ c2  1 . The curvature radius qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  

  c2 . The chords are q  q1;1 ¼ 2x ¼ 2 1  d2 d2 þ c2  1 and q0  q1;1 ¼ 2 at this point is given as R ¼ 8 1  d2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi   x2 þ y2 ¼ 2 d2 þ c2  1 1  d2 . If we define U as in Eq. (16), the trace of the deviation matrix is again given by Eq. (17). After substituting the calculated values of R, q, q0 and sin a, we obtain U ¼ 1 for all allowed shapes, thus Tr M = 2, and conclude that the hour-glass orbit is neutral. This means that in the TEB there is no stable hour-glass orbit, at variance with the ESB, where such an orbit having interesting properties was stable in a large region of the parameter space [17]. 3.2.6. Comparison of different descriptions of the parameter space The elongated truncated elliptical billiards were discussed in [44], where the boundary shapes were described with the help of two parameters h and a, related to our parameters d and c as

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi h ¼ 1  d2 ;



pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  d2

c

ð20Þ

In [44] the truncated elliptical billiards have been examined in a restricted space,ffi described as .pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi region of the parameter a > 1 and h < 1. In our parameters it corresponds to the region d < 1  c2 . A smaller region h < 1 1 þ a2 has been rigorously proved to be ergodic [44]. Written with our parameters, it obeys the conditions

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi c4 þ 4  c < d < 1  c2 2

ð21Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi  

or 1  d2 d < c < 1  d2 . The corresponding part of the parameter space in Fig. 2 is denoted with A. The comparison with our results shows that the lower limit in (21) is identical with the limit (14) for which the parabolic rectangular orbits emerge. The  .pffiffiffi hyperbolic behavior has been identified in the region h < min 1=a; 1 2 . In our parameters, this corresponds to the quasipffiffiffiffiffiffiffiffiffiffiffiffi triangular region in the parameter space, consisting of subregions A and B in Fig. 2. It is delimited by curves d ¼ 1  c, denoted . pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi in Fig. 2 by letter (j), d ¼ 1  c2 denoted as (a), and by the straight line d ¼ 1 2 denoted by letter (k). 4. The box-counting numerical analysis of the chaotic fraction of the phase space In this section, we return to the question of limits within which the TEB is chaotic. Finding the method able to determine exactly the chaotic volume fraction in the phase space has been considered one of the open problems in the theory of chaotic behavior in dynamical systems [51]. We wish to test the extent of chaos in the truncated elliptical billiard for a large number of shape parameters. Although analytical methods have been successfully applied to certain billiards, the numerical computation is often used to find the degree of chaoticity for large regions of the parameter space. One of possible methods is the computation of the chaotic fraction qclass of the classical phase space, which may be a significant indicator of the bifurcations and sharp transitions from integrable or mixed to chaotic behavior. The classical chaotic fraction can be calculated from the Poincaré surfaces of section, by means of the box-counting method. One usually assumes that the phase plane of a mixed Hamiltonian system consists of fixed points, some islands of stability surrounding these points and the rest of the phase space filled with chaotic trajectories. If one removes the islands and fixed points, the chaotic part measures the degree of chaoticity. However, as pointed out by Zaslavsky [28], the obtained result may not be correct because of the slow convergence of the method, due to the stickiness of some orbits. There are orbits which for a long time remain in the same region of the phase space and then suddenly start spreading. However, the method was successfully used in the investigation of a large family of billiards with circular arcs [15]. The statistical background for this method was investigated by Prosen and Robnik [52,53], with the conclusion that the box-counting method is fastest and the most reliable among the three numerical methods used for determining the value of the chaotic fraction [54]. In the quantum approach to chaotic systems, the conjecture that the corresponding chaotic fraction qquan can be obtained from the statistical analysis of the energy level densities has been largely accepted and applied to billiards [21,22]. In our investigation of the ESB, the method was applied to calculate the extent of the fully chaotic behavior and the results were in agreement with the analytic calculation [17,40]. Here we compute the chaotic component in the TEB using the boxcounting method for a large part of the parameter space, taking 0:00 6 d 6 1:00 and 0:01 6 c 6 2:20 with intervals of 0.01, thus for totally 100  220=22 000 shapes. The result is shown in Fig. 7. The corresponding results for the ESB for the same range of shape parameters are given in Fig. 8. The number of initial conditions is n1 ¼ 100 and the number of points in the Poincaré surfaces of section for each choice of initial values is n2 ¼ 5000. The third parameter used in the calculation for Figs.

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

318

Fig. 7. Diagram showing the degree of chaoticity qclass of the TEB, in dependence on the shape parameters, with 0 < c < 2:2 and 0 6 d 6 1. Black points denote shapes with qclass ¼ 1:00, red with 0:90 6 qclass < 1:00, green with 0:80 6 qclass < 0:90, blue with 0:70 6 qclass < 0:80, yellow with 0:60 6 qclass < 0:70 and grey with 0:00 6 qclass < 0:60. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)

Fig. 8. Diagram showing the degree of chaoticity qclass of the elliptical stadium billiards (ESB), in dependence on the shape parameters, with 0 < c < 2:2 and 0 6 d 6 1. The meaning of the black and colored points is the same as in Fig. 7. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)

7 and 8 is the number of squares (boxes) in the first quadrant, determined by the parameter n3 ¼ 100, and is n3  n3 ¼ 10 000. In Appendix A, we discuss in more detail the method of calculation, give the explanations for the choice of parameter values and present diagrams showing how the parameter choice affects the results. In Figs. 7 and 8, we are plotting the points representing the pairs of shape parameters in the ðd; cÞ plane. Points are plotted in different colors, depending on the corresponding value of qclass . The full chaos, corresponding to qclass ¼ 1:00, is depicted by black points. Colored points denote shapes within intervals between 0 and 0.999. Fig. 7 confirms that for the TEB, in the region below the onset of the stable two-bounce horizontal orbit, dynamics is practically completely chaotic. The chaotic part is sharply separated from the mixed part by the function c2 þ d2 ¼ 1 and suggests that the chaotic behavior extends over a region larger than the domain examined in [44]. This is in strong contrast with the behavior of the elliptical stadium billiard, for which the similar diagram is shown in Fig. 8. Due to the emergence of stable pantographic orbits, for the ESB the region of chaos was strictly bounded also from the lower side [40]. To examine the possible mechanism leading to these differences, we assume that the TEB and the ESB are two extremes and we search for a possible transition between them. 5. Generalised truncated stadium-like elliptical billiards In this section, we propose a new large class of stadium-like billiards which we call generalised truncated elliptical billiards (GTEB). Such a billiard depends on three shape parameters d, c and j. The allowed values of the shape parameters are

0 6 d 6 1;

0 < c < 1; 1 6 j 6 1

ð22Þ

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

319

For the limiting values of j we obtain the two billiard families considered before: for j = 1 GTEB reduces to the ESB, and for j = 1 GTEB becomes the TEB. The new GTEB boundary is obtained by adding elliptical arcs symmetrically at the two opposite ends of a rectangle with sides 2d and 2c. Elliptical arcs are cut out from the two identical but displaced ellipses by two horizontal straight lines at y ¼ c (Fig. 9). The two ellipses have centers at the points with Y  ¼ 0 and

1j X  ¼ d 2

ð23Þ

The distance between the two centers is D ¼ dð1  jÞ. The horizontal and vertical semiaxis are given, respectively, as

1j Ax ¼ 1  d 2

ð24Þ

 c 1j Ay ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  d 2 ð1  dÞð1 þ jdÞ

ð25Þ

and

The equation describing the two ellipses reads

2 2 x  X y þ ¼1 Ay Ax

ð26Þ

The horizontal diameter of the billiard is 2. For j ¼ 1 and c ¼ 1  d the GTEB becomes the Bunimovich stadium billiard. In Fig. 10, the Poincaré sections are shown for d ¼ 0:2 and c ¼ 0:6 with j assuming different values between 1 and 1. The four islands typical for the bow-tie orbit are present for all j except for j = 1 (TEB), where this orbit becomes parabolic and the island reduces to a characteristical flight of points. The limits separating chaotic from mixed behavior are determined by the onset of the stable horizontal 2-bounce orbit which is stable if R1 > 1. Since the curvature radius at jxj ¼ 1 is

R1 ¼

c2 ½2  dð1  jÞ 2ð1  dÞð1 þ djÞ

ð27Þ

the limit separating the chaotic from the mixed region is given as



sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1  dÞð1 þ jdÞ 2  dð1  jÞ

ð28Þ

In Fig. 11, the chaotic fraction qclass is shown for the special case d ¼ c for different j, in dependence on c. It is noticed that in the case j = 1 (ESB) there is a narrow, strictly limited region of full chaos, outside of which the values of qclass are low. For

Fig. 9. Construction and parameters of the generalised elliptical truncated stadium-like billiards (GTEB).

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

320

Fig. 10. Poincaré plots for the generalised truncated elliptical stadium-like billiards (GTEB), for d ¼ 0:2 and c ¼ 0:6 with different values of j: (a) j = 1.0; (b) j = 0.6; (c) j = 0.3; (d) j = 0; (e) j = 0.3; (f) j = 0.6; (g) j = 0.95; and (h) j = 1.0.

Fig. 11. Dependence of qclass on c for a set of shapes with d ¼ c and the choice of j as in Fig. 10.

j = 1 (TEB) the fully chaotic region is much larger and extends practically over all values d and c below the chaotic limit. Between these two limits, the regions of full chaos ðqclass ¼ 1Þ are shorter and limited, but there are many shapes with chaotic parameter close to 1 (between 0.90 and 0.999). This corresponds to a selection of extremely narrow islands in the Poincaré plots, noticed in Fig. 10(g).

6. Discussion and conclusions In conclusion, our investigation of the elliptical stadium-like billiards has revealed a rich variety of integrable, mixed and chaotic behavior, resulting from the character of the two elliptical arcs and their mutual position. This strong dependence on parameters d and c is confirmed for the TEB, but is even more enhanced when a third shape parameter j is added, providing a link between the two extreme cases, TEB and the previously treated ESB. The analysis, based on the box-counting method for calculating the chaotic fraction of the phase space, shows that among all considered shapes the TEB, created by cutting a single ellipse with two parallel straight lines, has exceptional properties, suggesting strongly that it is chaotic in a region sig-

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

321

Fig. 12. Dependence of qclass on the number n1 of initial conditions, number n2 of points in the Poincaré section for each pair of initial conditions and n3 determining the number n3  n3 of boxes in the first quadrant. The total number of points included in the box-counting procedure is n1  n2 .

Fig. 13. (a) Poincaré section for TEB with the shape parameters d ¼ 0:6 and c ¼ 1:3 and the corresponding empty and filled boxes with (b) n3 ¼ 100; (c) n3 ¼ 150; and (d) n3 ¼ 200.

nificantly larger than is the region examined formerly. Notable is the presence of many neutral orbits in this region. It is not surprising, since many of these orbits actually can be identified as orbits in an ellipse. However, neutral orbits are found also

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323

322

among the orbits with impact points on the flat segments of the boundary, increasing the chaotic regions in the parameter space. The stable islands appear only for flattened arcs, and are mostly due to the two-bounce horizontal orbit and to the diamond and multidiamond orbits. Our main result can be expressed as a conjecture that in the TEB there are no visible stable orbits with fixed points and corresponding islands in the region of elongated elliptical arcs, i.e. below the circular limit ffi pffiffiffiffiffiffiffiffiffiffiffiffiffi d < 1  c2 . Our investigations can be applied in designing the semiconducting optical devices and in the technology of microwave and acoustic resonant cavities, whose properties sensitively depend on the shape. With this purpose in mind, we propose the further analysis of generalised truncated stadium-like billiards with elliptical arcs and an extension of the present investigation to different types of open billiards. Acknowledgments The authors thank A. Bäcker, N. Pavin, T. Prosen, M. Robnik and T. Tanaka for discussions and useful comments, and V. Dananic´ and D. Radic´ for help with numerical procedures. This work was supported by MZOS of Croatia (Research Grant 125-1252971-2868). Appendix A. The box-counting method for computing the chaotic fraction Here we describe some details of the box-counting method used to estimate the degree of chaoticity for a particular billiard shape. The method is based on the calculation of the Poincaré surfaces of section. They contain a certain number of points dispersed within the ðX; V x Þ plane, where both coordinates can take values between 1 and 1. The Poincaré section for a chosen pair of shape parameters is computed starting with n1 sets of initial conditions and iterating each orbit for n2 intersections with the x-axis, thus obtaining n1  n2 points in the Poincaré diagram. The number n2 of intersection points, given in advance as an input parameter, is in most cases smaller than the number of bouncing points on the boundary, but in the input also the number of bouncing points can be limited. In Figs. 3, 4, 6 and 10, the computed Poincaré sections were shown, illustrating the structure of the phase space. In those computations it was important that both chaotic and quasiperiodic orbits were included. Therefore, the total number n1  n2 of points ðX; V x Þ had to be quite large, and we have set it to 100 000. The initial positions and velocities were taken from a list containing a random series of data. The question is whether it is appropriate to take the same or even greater numbers for the calculation of the chaotic fraction. Detailed testing has shown [43] that with the values n1 ¼ 100; n2 ¼ 5000 and n3 ¼ 100, used in our present calculation, one obtains reliable results. The computation proceeds as follows. The (+, +) quadrant of the phase plane is first divided into a grid of n3  n3 squares (boxes), then one counts the boxes which contain one or more points ðX; V x Þ and calculates the ratio of the number thus obtained to the total number of boxes. The obtained ratio is denoted as qclass . In this way also certain points belonging to invariant curves within the regular islands are included. Since our main aim is to examine precisely the onset of full chaos and to find only a qualitative estimate of the qclass outside the chaotic region, this method gives satisfactory results. It could seem questionable why one allows such a large number of orbits, when a single chaotic orbit should be sufficient. However, if only one orbit is included, it might be periodic or quasiperiodic and not appropriate for calculating the chaotic fraction. With a small number of orbits, two, four or ten, at least one among them might be chaotic, but this cannot be assured in advance, what is extremely important when one calculates, as in our case, 22 000 different shapes for each of the two billiard families. There are two important issues concerning this calculation: first, one looks for the shapes with total chaos, where all boxes are full. The obtained result will be more reliable if many orbits are included, covering together the white microregions in the Poincaré sections due to the fractality and stickiness of some orbits and those caused by the bouncing modes between the flat parts of the boundary or by other types of neutral orbits. On the other side, a larger number of periodic or quasiperiodic orbits leads to a larger deviation from the exact values for nonchaotic and mixed regions. In Fig. 12, the dependence of results for qclass on n1 and n3 is shown for the shape with d ¼ 0:6 and changing c, with n3 ¼ 100; n3 ¼ 150 and n3 ¼ 200. In Fig. 13, the filled and the empty boxes are shown, with the Poincaré section, for d ¼ 0:6 and c ¼ 1:5. Observing the inserted enlarged details in Fig. 12, one sees that the value n1 ¼ 100 is the most acceptable for the calculation of the chaotic part. Outside the chaotic region the results for different n3 differ in absolute values, but exhibit identical variations with c. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

Chernov N, Markarian R. Chaotic billiards. In: American Mathematical Society, mathematical surveys and monographs, vol. 127; 2006. Berry M. Regularity and chaos in classical mechanics, illustrated by three deformations of a circular billiard. Eur J Phys 1981;2:91–102. Bandres MA, Gutiérrez-Vega C. Classical solutions for a free particle in a confocal elliptic billiard. Am J Phys 2004;72:810–7. Sinai YaG. Dynamical systems with elastic reflections: ergodic properties of dispersing billiards. Russ Math Surv 1970;25:137–89. Bunimovich LA. On ergodic properties of certain billiards. Funct Anal Appl 1974;8:254–5. Robnik M. Classical dynamics of a family of billiards with analytic boundaries. J Phys A 1983;16:3971–86. Dullin HR, Bäcker A. About ergodicity in the family of limaçon billiards. Nonlinearity 2001;14:1673–87. Markarian R. New ergodic billiards: exact results. Nonlinearity 1993;6:819–41. Wojtkovski M. Principles for the design of billiards with nonvanishing Lyapunov exponents. Commun Math Phys 1986;105:391–414. Donnay VJ. Using integrability to produce chaos: billiards with positive entropy. Commun Math Phys 1991;141:225–57.

V. Lopac, A. Šimic´ / Commun Nonlinear Sci Numer Simulat 16 (2011) 309–323 [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54]

323

Bunimovich LA. Mushrooms and other billiards with divided phase space. Chaos 2001;11:802–8. Bunimovich LA. On the ergodic properties of nowhere dispersing billiards. Commun Math Phys 1979;65:295–312. Bunimovich LA. Criterion of absolute focusing for focusing component of billiards. Regular Chaotic Dynam 2009;14:42–8. Benettin G, Strelcyn J-M. Numerical experiments on the free motion of a point mass moving in a plane convex region: stochastic transition and entropy. Phys Rev A 1978;17:773–85. Dullin HR, Richter PH, Wittek A. A two-parameter study of the extent of chaos in a billiard system. Chaos 1996;6:43–58. Oliveira DFM, Leonel ED. On the dynamical properties of an elliptical–oval billiard with static boundary. Commun Nonlinear Sci Numer Simulat 2010;15:1092–102. Lopac V, Mrkonjic´ I, Pavin N, Radic´ D. Chaotic dynamics of the elliptical stadium billiard in the full parameter space. Physica D 2006;217:88–101. Lopac V, Mrkonjic´ I, Radic´ D. Chaotic dynamics and orbit stability in the parabolic oval billiard. Phys Rev E 2002;66(5):035202. Lopac V, Mrkonjic´ I, Radic´ D. Chaotic behavior in lemon-shaped billiards with elliptical and hyperbolic boundary arcs. Phys Rev E 2001;64(8):016214. Lopac V, Mrkonjic´ I, Radic´ D. Classical and quantum chaos in the generalized parabolic lemon-shaped billiard. Phys Rev E 1999;59:303–11. Bohigas O, Giannoni M-J, Schmit C. Characterization of chaotic quantum spectra and universality of level fluctuation laws. Phys Rev Lett 1984;52:1–4. Berry MV, Robnik M. Semiclassical level spacings when regular and chaotic orbits coexist. J Phys A 1984;17:2413–21. Gutzwiller MC. Chaos in classical and quantum mechanics. New York: Springer; 1990. Makino H, Harayama T, Aizawa Y. Effects of bifurcations on energy level statistics for oval billiards. Phys Rev E 1999;59:4026–35. Casati G, Prosen T. Mixing property of triangular billiards. Phys Rev Lett 1999;83:4729–32. de Menezes DD, Jar e Silva M, de Aguiar FM. Numerical experiments on quantum chaotic billiards. Chaos 2007;17(10):023116. Raedt HD, Katsnelson MI. Electron energy level statistics in graphene quantum dots. JETP Lett 2008;9:607–10. Zaslavsky GM. Physics of chaos in Hamiltonian systems. London: Imperial College Press; 2007. Gmachl C, Capasso F, Nöckels JU, Stone AD, Sivco DL, Cho AY. High-power directional emission from microlasers with chaotic resonators. Science 1998;280:1556–64. Tureci HE, Schwefel HGL, Douglas Stone A, Narimanov EE. Gaussian-optical approach to stable periodic orbit resonances of partially chaotic dielectric micro-cavities. Opt Exp 2002;10:752–61. Hentschel M, Richter K. Quantum chaos in optical systems: the annular billiard. Phys Rev E 2002;66(13):056297. Harayama T, Davis P, Ikeda KS. Stable oscillations of a spatially chaotic wave function in a microstadium laser. Phys Rev Lett 2003;90(4):063901. Lebental M, Djellali N, Arnaud C, Lauret J-S, Zyss J, Dubertrand R, et al. Inferring periodic orbits from spectra of simply shaped microlasers. Phys Rev A 2007;76(13):023830. Tanaka T, Hentschel M, Fukushima T, Harayama T. Classical phase space revealed by coherent light. Phys Rev Lett 2007;98(4):033902. Altmann EG, Motter AE, Kantz H. Stickiness in mushroom billiards. Chaos 2005;15(6):033105. Lansel S, Porter Mason A, Bunimovich LA. One-particle and few-particle billiards. Chaos 2006;16(12):013129. Ketzmerick R, Bäcker A, Löck S, Robnik M, Vidmar G, Höhmann R, et al. Dynamical tunneling in mushroom billiards. Phys Rev Lett 2008;100(4):174103. Lopac V, Movre I, Mrkonjic´ I, Radic´ D. Chaotic properties of the elliptical stadium billiard. Prog Theor Phys Suppl 2003;150:371–5. Del Magno G, Markarian R. Bernoulli elliptical stadia. Commun Math Phys 2003;233:211–30. Canale E, Markarian R, Oliffson Kamphorst S, Pinto de Carvalho S. A lower bound for chaos on the elliptical stadium. Physica D 1998;115:189–202. Markarian R, Oliffson Kamphorst S, Pinto de Carvalho S. Chaotic properties of the elliptical stadium. Commun Math Phys 1996;174:661–79. Oliffson Kamphorst S, Pinto de Carvalho S. Elliptic islands on the elliptical stadium. Discr Cont Dynam Syst 2001;7:663–74. Šimic´ A. B.Sc. Thesis, Faculty of Science, University of Zagreb, Zagreb; 2009. Del Magno G. Ergodicity of class of truncated elliptical billiards. Nonlinearity 2001;14:1761–86. Zaslavsky GM, Strauss HR. Billiard in a barrel. Chaos 1992;2:469–72. Bunimovich L, Casati G, Guarneri. Chaotic focusing billiards in higher dimensions. Phys Rev Lett 1996;77:2941–4. Mohr MB, Seemann G, Sachse FB, Dössel O. Deformation simulation in an elastomechanical ventricular model. Comput Cardiol 2004;31:777–80. Heller EJ, Tomsovic S. Postmodern quantum mechanics. Phys Today 1993;46(7):38–46. Ree S, Reichl LE. Classical and quantum chaos in a circular billiard with a straight cut. Phys Rev E 1999;60:1607–15. Michel C, Tascu S, Doya V, Legrand O, Mortessagne F. Gain-controlled wave chaos in a chaotic optical fibre. J Eur Opt Soc – Rapid Publ 2009;4:09020. Bunimovich LA. Relative volume of Kolmogorov–Arnold–Moser tori and uniform distribution, stickiness and nonstickiness in Hamiltonian systems. Nonlinearity 2008;21:T13–7. Prosen T, Robnik M. General Poissonian model of diffusion in chaotic components. J Phys A: Math Gen 1998;31:L345–53. Robnik M. Randomness of classical chaotic motion. Progr Theor Phys Suppl 2003;150:229–42. Robnik M, Dobnikar J, Rapisarda A, Prosen T, Petkovšek M. New universal aspects of diffusion in strongly chaotic systems. J Phys A: Math Gen 1997;30:L803–13.