Fractal Weyl law for three-dimensional chaotic hard-sphere scattering systems
FFractal Weyl law for three-dimensional chaotic hard-sphere scattering systems
Alexander Ebersp¨acher, J¨org Main, and G¨unter Wunner
Institut f¨ur Theoretische Physik 1, Universit¨at Stuttgart, 70550 Stuttgart, Germany (Dated: October 26, 2018)The fractal Weyl law connects the asymptotic level number with the fractal dimension of thechaotic repeller. We provide the first test for the fractal Weyl law for a three-dimensional openscattering system. For the four-sphere billiard, we investigate the chaotic repeller and discussthe semiclassical quantization of the system by the method of cycle expansion with symmetrydecomposition. We test the fractal Weyl law for various symmetry subspaces and sphere-to-sphereseparations.
PACS numbers: 05.45.Mt, 03.65.Sq, 42.25.p
I. INTRODUCTION
The asymptotic eigenvalue distribution of partial dif-ferential equations such as the Schr¨odinger equation forfree particles or the one-dimensional Helmholtz equationfor sound waves has been of interest as early as from 1912on when Hermann Weyl and Richard Courant first stud-ied the problem [1, 2]. They found expressions for theasymptotic level number N ( k ) in closed systems to beproportional to k d , with d the spatial dimension of thesystem. The so called “Weyl law”, which has been wellknown from then on [3], states that for closed quantumsystems, every accessible Planck cell in phase space is oc-cupied by one quantum state. A generalization to chaotic open systems, where complex resonances k n = ¯ k n − Γ n / k n and lifetimes Γ n replace real eigen-values k , has been proposed in the 1990s [4, 5]. Thenumber of resonances N ( k ) = { k n : Re( k n ) ≤ k ; Im( k n ) > − C } (1)inside a rectangle in the complex plane defined by theenergy k and the strip width C is conjectured to be pro-portional to k α with the exponent α = D + 12 (2)being related to the non-integer fractal dimension D ofthe chaotic repeller . The number α takes the role ofthe effective number of degrees of freedom. The repelleris the set of all classical trajectories that stay trappedfor t → + ∞ or t → −∞ . Considering only the stablemanifold W s of trapped trajectories for t → + ∞ in asuitable Poincar´e surface of section, the fractal Weyl lawreads [6] N ( k ) ∝ k d H +1 , (3)with d H the Hausdorff dimension of the chaotic repeller’sstable manifold.The fractal Weyl law (3) has been investigated for var-ious two-dimensional systems, e.g. a triple Gaussian po-tential [7], the three-disk billiard [6], an optical microsta-dium resonator [8] and a modified H´enon-Heiles potential FIG. 1. The four-sphere billiard. Four spheres of equal radius R are located on the vertices of an equilateral tetrahedron(indicated by bars) with edge length d . Shown is the case of d/R = 6. [9] as well as for quantum maps, e.g. the kicked rotator[10, 11]. The systems under consideration so far have allbeen at most two-dimensional. We provide a first inves-tigation of a three-dimensional system in this paper.The problem under consideration is scattering in thefour-sphere billiard. This system is characterized by fourspheres of the same radius R located on the vertices of anequilateral tetrahedron with edge length d as visualized inFig. 1. The relevant configuration parameter is the ratio d/R . Different from e.g. the two-dimensional three-diskbilliard, the four-sphere billiard is an open system evenfor the case of touching spheres at d/R = 2.The paper is organized as follows. The fractal dimen-sion of the repeller is computed in Sec. II. Periodic-orbittheory and cycle expansion methods are applied in Sec.III to obtain semiclassical resonance spectra. The clas-sical escape rate is determined in Sec. IV. Results whichallow for a test of the fractal Weyl law are presented inSec. V. Concluding remarks are given in Sec. VI. II. GAUGING THE REPELLER
As the fractal dimension D of the chaotic repeller en-ters into the fractal Weyl law (3), it is crucial to deter-mine D accurately. The fractal structures in the four- a r X i v : . [ n li n . C D ] A ug sphere billiard have already been studied experimentallyand theoretically [12–15], however, all investigations sofar have been limited to small values of the configurationparameter d/R (cid:46) .
5. In this paper the fractal dimen-sion of the repeller is determined accurately for a widerange of the parameter d/R . A. Fractal repeller
The stable manifold W s of the chaotic repeller is afractal in phase space. The time spent in the scatter-ing system can be measured by the time-delay function T counting the number of reflections experienced by atrajectory. Choosing initial conditions in a plane parallelto the plane spanned by three of the spheres, it is possi-ble to iterate trajectories entering the scattering systemsuch that the fourth sphere is visited first. For large ra-tios d/R , the boundary to the region that contains thoseinitial conditions is the projection of the fourth sphereonto the plane from which the trajectories are iterated.As there are three distinct possibilities to visit the nextsphere, there are three regions of higher values of T in-side this circle. Repeating this line of argument for anyregion belonging to a given visitation sequence, the struc-ture of the fractal repeller can be understood. Figure 2illustrates the fractal structure.In principle, the Hausdorff dimension d H can be calcu-lated from box-counting. This procedure, however, is notsuited for the billiards under consideration as it requiresiteration of a vast amount of initial conditions on a grid.The regions of high values of T that exhibit fractal prop-erties may not be resolved with acceptable computationaleffort. A method better suited to billiards is introducedbelow. B. Estimating d H through Hausdorff sums Finite numerical precision and finite computing timeavailable prevent the determination of initial conditionsthat lead to trapped orbits. Instead, it is possible to es-timate the Hausdorff dimension from regions of finite T .We introduce the auxiliary quantity A ( i ) n that denotes thearea of the n -th region of initial conditions with T ≥ i reflections in the surface of section, and define the quan-tities K ( i ) ( s ) := (cid:88) n (cid:16) A ( i ) n (cid:17) s (4)which will be called Hausdorff sums below. These sumshave the following properties [16, 17]:lim i →∞ K ( i ) ( s ) = ∞ for 0 ≤ s < d H const. > s = d H d H < s < ∞ (5) (a) -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2x-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 y (b) -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02x 0.105 0.11 0.115 0.12 0.125 0.13 0.135 0.14 0.145 y FIG. 2. (a) Time-delay functions in the range 2 ≤ T ≤ d/R = 3 and (b) a magnification thereof in the range 2 ≤ T ≤
7. The time-delay functions are drawn as functions of thecoordinates x, y in the surface of section. The colors indicatethe value of T and for clear identification some regions areexplicitly labeled with T . Only the boundaries of individualregions are drawn; initial conditions chosen within the regionsexperience the same number of reflections. The self-similaritysuggests that the stable manifold is a fractal set. The property for s = d H stems from the fact that theHausdorff sums by definition are smooth functions of thevariable s . This allows one to estimate the Hausdorffdimension d H by intersecting K ( i ) ( s ) for different i [14,15].Existing methods of estimating the repeller’s fractal di-mension have been confined to a narrow range of the pa-rameter d/R , in particular to the case of almost-touchingspheres [14, 15]. The algorithm discussed in the followingallows for calculations in a wider range of d/R . To esti-mate the stable manifold W s in the surface of initial con-ditions, it will be necessary to accurately find boundariesof regions of a given visitation sequence of the spheresand reflection count. Once both a point inside the regionin question and a point outside are known, interval bi-section may be used to compute the boundary point onthe line connecting both points. The bisection conditionuses the visitation sequence and the reflection count T ,i.e. the length of the symbolic code.
1. Finding regions
The algorithm used in this paper relies on the structureof the time-delay functions discussed in Sec. II A. In thecalculations, the following assumptions are made: • Regions of a specific order of visits with the scat-terers described by the symbolic code are non-overlapping. • Within a region of a given order of visits, there areexactly three more regions each corresponding toadditional visits at one of the three other spheres.Both conditions may be violated for d/R close to 2, i.e.the case of almost touching spheres. We find the assump-tions to be fulfilled for regions with T ≥ d/R (cid:38) . z -axissuch that the uppermost sphere is visited first. Let usassume that regions with T min ≤ T ≤ T max are sufficientfor an estimation of d H . Under this assumption, it ispossible to find regions approximating the repeller withthe following procedure.In the a first step, the projection of the sphere visitedfirst onto the surface of section is determined. For small d/R , this region is a circle, for larger d/R , projectionsof the other spheres may be cut out of the circle. Thisis done by randomly choosing a point in the surface ofsection and an interval bisection between this point andpoints equally distributed on a large circle fully contain-ing the projection of sphere 4. By assumption, insidethis region there are three other regions with T = 2. Thecorresponding visitation sequences differ in the secondcharacter. Once a point inside each of the regions with T = 2 is found, in a second step, polygonal chains form-ing boundaries to each of the new distinct regions arecalculated. This procedure is iterated until all desiredregions corresponding to T min ≤ T ≤ T max have beenfound.
2. Areas from polygonal chains
One possibility to store the boundaries is by keepinga polygonal chain. As in this procedure the number N regions of regions grows exponentially with N regions =3 T − , this way of data storage is memory-expensive.However, the areas enclosed by the polygons are easily TABLE I. Numerical values of the Hausdorff dimensions d H ofthe stable manifold W s for various configuration parameters d/R . All decimal digits are significant. d/R . d H . . . . . . . calculated using numerical quadrature of the area givenby A = 12 π (cid:90) r ( ϕ ) d ϕ , (6)with r ( ϕ ) the distance of the boundary point from the“midpoint” of the region. A fairly low number of sup-porting points has proved to be sufficient for very highprecision. All calculations have been performed with 101supporting points.
3. Areas from ellipses
An alternative to the memory-expensive storage ofpolygonal chains is to approximate the boundary by anellipse described by the polynomial a x + a y + a xy + a x + a y = 1 . (7)For five or more known points ( x i , y i ) of the boundarythe coefficients a , . . . , a can be determined from a lin-ear least-squares fit. The semi-major axes a, b , the centershift ( x , y ) as well as the rotation angle ϕ of the ellipsescan be easily extracted from the coefficients of the poly-nomial in Eq. (7).Fitted ellipses allow to improve accuracy as it is nowpossible to shift the regions’ “midpoints” used in the con-struction of the polygonal chains to the midpoint of theellipses. All bisections for the polygonal chain boundariesare repeated in such a way that all lines connecting theellipse’s midpoints and the boundary points intersect atidentical angles. This will be beneficial for the quadra-ture of the areas entering into the Hausdorff sums.Once all desired boundaries have been calculated, thefractal dimension d H can be estimated. To build theHausdorff sums (4), the areas enclosed in the individ-ual regions have to be known. From the ellipses fitted tothe boundaries, the area A is trivially given by A = πab , (8)where a and b are the semi-major axes.Calculations have been performed for d/R = 2 . d/R = 10. Sample plots for intersected Hausdorff sums K ( i ) ( s ) are shown in Fig. 3. Results for the Hausdorffdimension d H are compiled in Fig. 4 and in Table I. Thecalculations using polygonal chains agree up to four dec-imal digits with the calculations using fitted ellipses. -0.9-0.8-0.7-0.6-0.5-0.4-0.3 0.23 0.232 0.234 0.236 0.238 0.24 Loga r i t h m o f H au s do r ff s u m s s(a) T=4T=5T=6T=7T=8T=9T=10 -0.9-0.8-0.7-0.6-0.5-0.4-0.3 0.23 0.232 0.234 0.236 0.238 0.24 Loga r i t h m o f H au s do r ff s u m s s(b) T=4T=5T=6T=7T=8T=9T=10 FIG. 3. Intersected Hausdorff sums K ( i ) ( s ) for various re-flection numbers i = T calculated from (a) a polygonal chainwith 101 supporting points and (b) ellipses fitted to polygonalchains for d/R = 6. As can clearly be seen, the intersectionpoints for all shown curves agree perfectly. For this reason,the Hausdorff dimensions can be determined to a precision ofat least 4 significant digits. Figure 4 clearly shows that with decreasing d/R theintersection of the stable manifold W s with the Poincar´esurface of section fills the plane denser. The repeller’s di-mension d H thus increases as the tetrahedron gets packedmore densely.In summary, the method presented here establishesa fast and very precise method of gauging the repeller.Though the assumptions are quite strong, they hold overa wide range of the ratio d/R . III. SEMICLASSICAL RESONANCES
Studying a billiard system in a purely quantum me-chanical fashion turns out to be intricate since – althoughparticles move freely in between the scatterers – it is ademanding task to find wavefunctions that vanish on allscatterer’s boundaries simultaneously. For attempts on N -sphere scattering systems in three dimensions, see [18];for the two-dimensional three-disk scattering system see[19]. Results and comparisons of methods for the four-sphere scattering system are presented in [20]. The tech- H au s do r ff d i m en s i on d H d/R FIG. 4. Hausdorff dimension d H of the stable manifold W s asfunction of the ratio d/R . All data points have been obtainedby intersecting Hausdorff sums as demonstrated in Fig. 3. niques of semiclassical quantization presented below arebetter suited for billiard systems.As one of the great achievements of semiclassicalphysics, Gutzwiller’s trace formula provides a mean ofquantizing a system via periodic orbits [21]. Unfortu-nately, the trace formula is plagued by serious conver-gence problems. In chaotic systems, the number of pe-riodic orbits typically grows exponentially with length l ,and this growth usually cannot be compensated by thedecrease of the amplitude factors. A method of improv-ing the convergence ideally suited for billiard systems isbased on the Gutzwiller-Voros zeta function.The logarithmic derivative of the function Z ( k ) = (cid:89) n ( k − k n ) (9)with the quantized wavenumbers k n of a billiard systemyields the density of states g ( k ) = − π Im dd k ln Z ( k ) = − π Im (cid:88) n k − k n + i (cid:15) . (10)Voros [22] proposed a semiclassical formulation of Z ( k )which reads for billiard systems Z GV ( k ) = exp (cid:32) − (cid:88) p ∞ (cid:88) r =1 r ( − rn p e i rl p k (cid:112) | det( M pr − ) | (cid:33) , (11)where l p is the length of a primitive periodic orbit and n p the number of reflections on hard wall boundaries. Theindex r counts the number of repetitions of a primitiveperiodic orbit. The reduced monodromy matrix M p pro-vides information on the linear evolution of a small devi-ation from an initial condition belonging to a periodic or-bit over one period [23]. The eigenvalues of M p quantifythe stability of the periodic orbit. Due to the symplec-tic structure of Hamiltonian mechanics, the eigenvaluescome in tuples λ , 1 /λ , λ ∗ , and 1 /λ ∗ . TABLE II. Character table for the group T d [20]. The grouphas five irreducible representations: the one-dimensional rep-resentations A and A , the two-dimensional representation E and the two three-dimensional representations T and T . T d E C C S σ d A A − − E − T − − T − − For systems with symbolic dynamics such as bil-liards, the method of cycle expansion [24–27] has provedto be especially successful. A cycle expansion of theGutzwiller-Voros zeta function Z GV ( k ) in Eq. (11) isachieved by replacing ( − rn p in Eq. (11) by the term( − z ) rn p depending on the book-keeping variable z , ex-panding Z GV as a power series in z and then truncatingthe series. The highest power of z equals the maximumcycle length n max contributing to the cycle expansion.After truncation, z has to be set to z = 1. The cycle-expanded zeta function has better convergence behaviorover the trace formulas as individual terms tend to can-cel. A. The symmetry group T d The four-sphere billiard has discrete tetragonal sym-metry [20]. The associated symmetry group T d containsall symmetry operations that leave a regular tetrahedroninvariant. In particular, there are the identity operation E , 4 rotations C by 2 π/ C by 4 π/ C by π around theaxes intersecting the middle points of opposing edges,6 reflections σ d at planes perpendicular to the tetrahe-dron’s edges and also containing another vertex; and,furthermore, 3 permutations of the vertices S which canbe written as a combination of a rotation C by π/ σ h at the plane perpendicular to the mainrotation axis, i.e. the axes of C . Finally, the symmetrygroup T d also contains 3 distinct three-times repeatedrotary reflections S .The character table of the symmetry group is givenin Table II. The symmetry group can be decomposedinto 5 invariant subspaces, i.e. the representation matri-ces D of T d can be decomposed into block-diagonal formwhere the diagonal elements contain matrices represent-ing the group’s elements. The representation is called“irreducible” if no further decomposition is possible.From the character table, statements on the wavefunc-tions ψ can be made. For the one-dimensional representa-tions A and A , the effect of the symmetry transforma-tions is described by a multiplication with the character χ . For the totally symmetric A subspace, all charactersare equal to 1, therefore the wavefunctions have the fullsymmetry of T d , i.e. ψ is not altered by any symmetrytransformation. In the A subspace, the wavefunctionchanges sign under the reflection σ d and permutation S .For representations of higher dimension, the effect of thesymmetry transformations cannot be described in such asimple way.We note that repeated application of symmetry trans-forms may be identical with other elements of the sym-metry group, e.g., σ d = C = C = S = E and S = C . These identities will be useful for the symme-try decomposition of zeta functions discussed below. B. Symbolic dynamics and periodic orbits in thefour-sphere scattering system
The method of cycle expansion requires all periodicorbits up to a given maximum cycle length n r . In specialbilliards it is convenient to assign a symbolic code to eachperiodic orbit. For the four-sphere scattering system thesymbolic dynamics and periodic orbits have already beenintroduced in Ref. [20]. For the convenience of the readerwe briefly recapitulate the central ideas and properties ofthe orbits.In the four-sphere scattering system, periodic orbitsare determined by the periodic sequence in which the fourscatterers are visited. For most cases, any sequence cor-responds to one periodic orbit. If this one-to-one corre-spondence is not given, i.e. if some orbits become unphys-ical because they penetrate the scatterers, one speaks of pruning . This is the case for small center-to-center sep-arations. In the four-sphere scattering system, the sym-bolic dynamics has been shown to be pruned for config-urations with d/R < .
1. Periodic orbits in the fundamental domain
In systems with discrete symmetries such as the four-sphere billiard, which is invariant under all symmetryoperations of the tetrahedron group T d , whole classes oforbits are equivalent to each other, e.g., the 6 orbits whichare scattered back and forth between two spheres can bemapped onto each other using symmetry operations of T d . Furthermore, cyclic permutation of the sequence ofspheres leaves the orbits invariant. For these reasons, itis appropriate to use the symmetry properties to intro-duce the following short notation [20]. First, define theplane of reflection as the plane that contains the centersof the last three distinct spheres visited. Then, instead oflabeling all spheres individually, the label 0 will be used ifthe orbit visits the last sphere once more, 1 will indicatea visit at the third other sphere in the same plane of re-flection, whereas the label 2 will be used for a visit at thefourth sphere outside the plane of reflection. With thisnomenclature, all symbolic codes containing the charac- TABLE III. Primitive periodic orbits up to cycle length n ˜ p =2 for d/R = 4. The reduced symbolic code ˜ p as well as thesymmetry h ˜ p of each cycle is given. Furthermore, the real andimaginary parts of the stability eigenvalues λ ( i )˜ p are tabulated.All numbers have been rounded to five decimal digits. Theshortest cycle, labeled by 0, which visits two spheres in turnshas ambiguous symmetry. Both the rotation about π , C , aswell as the reflection about the plane perpendicular to the lineconnecting the sphere’s center, σ d , map this particular orbitonto itself.˜ p h ˜ p l ˜ p Re λ (1)˜ p Im λ (1)˜ p Re λ (2)˜ p Im λ (2)˜ p σ d , C C S σ d C S ter 2 are three-dimensional, whereas orbits correspondingto sequences of 0 and 1 are two-dimensional. The newlabeling reduces the number of characters in the alphabetto three, i.e. the code is ternary. This reduction corre-sponds to a reduction of the full physical phase space M to the fundamental domain ˜ M from which the wholephase space can be reconstructed by applying the symme-try group’s elements. Note that the symmetry reducedorbits are in general shorter than the corresponding phys-ical ones. Only the symmetry reduced orbits that havethe identity operation E as maximum symmetry have thesame length as the corresponding physical orbits. The re-duced orbits of symmetry classes σ d and C yield twiceas long physical orbits, C orbits are three times longerthan in the physical space, and, finally, S orbits havequadruple length.
2. Finding periodic orbits
Periodic orbits of the four-sphere scattering system arecalculated by varying, for a given symbolic code, the re-flection points on the spheres until the length of the or-bit takes its minimum value. Details of the numericalperiodic orbit search are described in [20] and the com-putation of the monodromy matrix is explained in Refs.[28, 29]. Table III lists the first few periodic orbits in thefundamental domain as well as their properties for theratio d/R = 4. The table also gives the maximum sym-metry operation that leaves the orbit invariant. In thefundamental domain, this operation corresponds to theoperation that maps the endpoint of the orbit in the fun-damental domain onto the starting point. Note that, forexample, the cycle t which visits two spheres in turnsis periodic in the fundamental domain, but not in fullphysical space. In the full domain, the start and end-point of the 0-cycle are not identical; application of therotation C respectively the reflection σ d yields back the full periodic orbit. C. Discrete symmetries and cycle expansion
In systems with discrete symmetries the full physicalspectrum can be decomposed into spectra belonging todifferent representations of the symmetry group. Thediscrete symmetries lead to symmetry factorized zetafunctions, which allow for the computation of quantumspectra belonging to a specific symmetry subspace. Thesymmetry decomposition of zeta functions has been elab-orated by Cvitanovi´c and Eckhardt [30] and exampleshave been given for the symmetry groups of various two-dimensional N -disk pinball models. Here, we present ex-plicit results for the tetrahedron group T d of the three-dimensional four-sphere scattering system. As in [30] wefirst discuss the symmetry decomposition of the dynam-ical zeta function Z ( k ) = (cid:89) p (1 − t p ( k )) (12)which is obtained from the Gutzwiller-Voros zeta func-tion (11) with the approximation | det( M pr − ) | − / ≈ | λ (1) p λ (2) p | − r/ = e − ( u (1) p + u (2) p ) r/ (13)and the definition t p ( k ) = e i l p k − ( u (1) p + u (2) p ) / , (14)and then generalize the results for the symmetry decom-position of the Gutzwiller-Voros zeta function.In quantum mechanics, the full Hilbert space H of theproblem factorizes into subspaces belonging to certainirreducible representations of the symmetry group, i.e. H = H A ⊗ H A ⊗ H E ⊗ H T ⊗ H T , (15)for the four-sphere scattering system. In [30] it is pointedout that zeta functions can be factorized in a similar way.The fundamental domain of phase space is sufficient forall computations, as the whole phase space M can beobtained from the fundamental domain ˜ M by M = (cid:88) h ∈ G h ˜ M , (16)where G is the symmetry group. Evaluating traces oftransfer operators in the fundamental domain ˜ M , thissymmetry reduction results in [26, 30](1 − t p ) m p = det ( − D ( h ˜ p ) t ˜ p ) , (17)with m p the multiplicity of a primitive cycle p .These expressions could be evaluated using a certainexplicit representation D α ( h ) of the group’s symmetryoperations h . However, this is a computationally ratherdemanding endeavor. Instead, the determinants can beexpressed in terms of traces χ which can be read off fromthe symmetry group’s character table (see Table II). Forexample, the expansion of det( − D ( h ) t ) for dimension d = 3 readsdet( − D ( h ) t ) = 1 − χ ( h ) t + 12 (cid:0) χ ( h ) − χ (cid:0) h (cid:1)(cid:1) t + 16 (cid:0) χ ( h ) − χ ( h ) χ (cid:0) h (cid:1) + 2 χ (cid:0) h (cid:1)(cid:1) t , (18)where the trace of D ( h ) is as usual denoted by χ ( h ).Carrying out this procedure explicitly, one obtains thefactorizations given in Table IV. Thus, the zeta functionin Eq. (12) can be rewritten in a symmetry reduced ver-sion Z α = (cid:89) ˜ p (1 − D α ( h ˜ p ) t ˜ p ) (19)for the subspace α . The zeta function now depends onlyon the fundamental cycles ˜ p . By this procedure, a fac-torization Z ( k ) = (cid:89) α Z α ( k ) d α (20)is achieved. The zeta function Z factorizes into zeta func-tions belonging to certain irreducible representations α ofthe symmetry group. The dimensions d α of the represen-tations enter into the full zeta function – and with them,the quantum multiplicities of resonances belonging to acertain subspace.
1. Assigning weight factors
The method of cycle expansion expands the zeta func-tion Z into a truncated series in which all cycles up toa certain cutoff length enter [24–27]. However, besidesthe primitive cycles, also multiple traversals contribute.Therefore, it needs to be clarified how repeated revolu-tions can be taken into account. Let us assume the prim-itive fundamental cycles ˜ p are known. Then, the contri-bution of an r -times repeated revolution to the symmetryreduced zeta function (19) is given by polynomials suchas (1 − z r t ˜ p,r ) , (21)where a dummy variable z has been introduced. The cy-cle weights t ˜ p,r have the form of the terms in (11) andare thus easily calculable from the cycle weight t ˜ p of theprimitive fundamental cycle. By using the factorizationsgiven in Table IV, it is possible to determine the weightfactor w ˜ p,r ( k ) as the sum of all roots z ir of the polyno-mials given in the table, w ˜ p,r = (cid:88) i z ir . (22) If this is possible, a way to use the ˜ p for repetitions aswell has been found. As an example, for the contributionof the r -times repeated cycle ˜ p to the A spectrum, weneed to solve (1 − z r t ˜ p,r ) = 0 , (23)which is true for z r = 1. Thus, in the A subspace, allweight factors are w ˜ p,r = 1. By this choice, the symmetryfactorization is retained. As another example, considerthe E subspace for cycles with C symmetry. Here, so-lutions to the equation (cid:16) z r t ˜ p,r + z r t ˜ p,r (cid:17) = 0 (24)are needed. A factorization is given by (cid:16) − e π i r/ t ˜ p,r (cid:17) (cid:16) − e − π i r/ t ˜ p,r (cid:17) = 0 , (25)where the exponentials are the roots z i . Evaluatingthe sum z r + z r , we find the weight factors w ˜ p,r = − , − , , − , − . . . for r = 1 , , . . . . A short notationfor this sequence is given by w ˜ p,r = 2 cos(2 πr/ w ˜ p,r given in Table Vare determined.
2. Ambiguous symmetry
The shortest cycle labeled by 0 in the four-sphere sys-tem has ambiguous symmetry. It is possible to map thiscycle onto itself by both the rotation C and the reflec-tion σ d . This ambiguity requires special care in the sym-metry decomposition. This is particularly important asthe 0-cycle is one of the fundamental cycles that act asbuilding block for longer cycles in the sense of cycle ex-pansion. The group theoretical weight of the 0-cycle canbe written as [30] h = C + σ d . (26)The symmetry factorization can thus be not one of thosegiven in Table IV. However, it is possible to use a fac-torization that contains factors in such a way that thefactorization is at most the greatest common divisor ofthe factors given for C and σ d in Table IV. The factor-izations and weight factors w ,r are given in Table VI.With that factorization the product (20) of all zeta func-tions belonging to the irreducible representations of thesymmetry groups coincides with the cycle expansion (12)using all orbits in the full domain.The final form of the Gutzwiller-Voros zeta functionwe use for our calculations is Z GV; α ( k ) = exp − (cid:88) ˜ p ∞ (cid:88) r =1 r w ˜ p,r ; α ( − z ) rn ˜ p e i rl ˜ p k (cid:112) | det( M ˜ pr − ) | , (27) TABLE IV. Symmetry factorization of the zeta function Z for all five irreducible representations of the group T d . The tableentries give the contribution of each fundamental cycle ˜ p to the Euler product Z = (cid:81) ˜ p (1 − t ˜ p ). This factorization allows thecomputation of quantum spectra for each symmetry subspace. E C C S σ d A (1 − t ˜ p ) (1 − t ˜ p ) (1 − t ˜ p ) (1 − t ˜ p ) (1 − t ˜ p ) A (1 − t ˜ p ) (1 − t ˜ p ) (1 − t ˜ p ) (1 + t ˜ p ) (1 + t ˜ p ) E (1 − t ˜ p ) (1 + t ˜ p + t p ) (1 − t ˜ p ) (1 − t ˜ p )(1 + t ˜ p ) (1 − t ˜ p )(1 + t ˜ p ) T (1 − t ˜ p ) (1 − t ˜ p )(1 + t ˜ p + t p ) (1 − t ˜ p )(1 + t ˜ p ) (1 − t ˜ p )(1 + t p ) (1 − t ˜ p )(1 + t ˜ p ) T (1 − t ˜ p ) (1 − t ˜ p )(1 + t ˜ p + t p ) (1 − t ˜ p )(1 + t ˜ p ) (1 + t ˜ p )(1 + t p ) (1 + t ˜ p )(1 − t ˜ p ) TABLE V. Weight factors w ˜ p,r for r traversals of the primitivecycle ˜ p . These factors allow for symmetry factorizations withrepetitions of primitive fundamental cycles. E C C S σ d A A − r ( − r E πr − r − r T πr − r πr − r T πr − r ( − r + 2 cos πr − r TABLE VI. Factorizations of Z ( k ) and weight factors w ,r for the fundamental cycle 0 with ambiguous symmetry classes C , σ d in all subspaces α of the symmetry group T d . C , σ d w ,r A (1 − t ) 1 A E (1 − t ) 1 T (1 + t ) ( − r T (1 + t )(1 − t ) 1 + ( − r with ˜ p the primitive symmetry reduced cycles, r the num-ber of repetitions of it and α the symmetry subspace.A symmetry reduced version of the cycle expansion isobtained by expanding Eq. (27) into a power series in z which is truncated at a maximum cycle length n max .Then, z has to be set to z = 1. D. Harmonic inversion method
The Gutzwiller-Voros zeta function Z ( k ) in Eq. (27)contains all energy eigenvalues k as complex zeros, and,in principle, it is possible to obtain spectra by a numeri-cal root search. This method has been successfully usedfor billiards, see e.g. [27]. However, the root search incycle expansions of high order is numerically expensive.For statistical purposes it is important not to miss anyresonances in the strip of the complex plane under consid-eration. Therefore, a dense grid of initial root guesses hasto be used for the root search. Consequently, many reso-nances will be found several times. Thus, the problem is to distinguish for every new root whether a new distinctresonance has been found or if the new zero has alreadybeen computed. As the number of resonances enters intothe fractal Weyl law (3) through the counting functions N ( k ), it is crucial to count individual resonances onlyonce.An alternative to the computation of zeros is based onthe harmonic inversion method for high-resolution spec-tral analysis [31–33]. When Eq. (10) is evaluated alonga line of real-valued k or a shifted line k + i δ with real k and δ , we obtain a spectrum g ( k ) = (cid:88) n π Γ n / δ (cid:0) k − ¯ k n (cid:1) + (Γ n / δ ) (28)which is a superposition of resonances with a Lorentzianshape. For negative shifts δ the Lorentzians are locatedat the positions ¯ k n , but with reduced widths Γ n + 2 δ .The basic idea is now to reformulate, via a Fourier trans-form, the problem of extracting eigenvalues as a signalprocessing task. Details of the method are given in [8].The procedure of calculating quantum spectra is sum-marized as follows: First, the spectrum g ( k ) is calculatedas a superposition of Lorentzians. We use the cycle-expanded zeta function Z ( k ) for this purpose. The quan-tity g ( k ) = − π Im dd k ln Z ( k ) = − π Im Z (cid:48) ( k ) Z ( k ) (29)is evaluated along lines parallel to the real axis with dif-ferent shifts δ . Thus, the shifts that allow for betterresults in harmonic inversion enter into the cycle expan-sion. Then, harmonic inversion is used to obtain theLorentzians’ parameters ¯ k n and Γ n for spectra calculatedwith different shifts. In the next step, the spectra arefiltered via the amplitudes. The quantity ˜ g ( k ) given in(29) should give resonances with an amplitude of A n = 1.True resonances with amplitudes A n ≈ δ arejoined such that the individual strips do not overlap. IV. CLASSICAL ESCAPE RATE
The classical escape rate γ can be interpreted descrip-tively as follows [34]: presume the scattering system un-der consideration is located in a box much larger thanthe system itself. Conducting N scattering experimentswith the same incident energy k , but different incidentdirections, one finds that the number N t of trajectoriesthat are inside the box after the time t has passed decaysexponentially as N t ∝ N e − γ t . (30)The relation of the escape rate and the imaginary partof the quantum resonances can be understood from thecorrespondence principle. The number of classical trajec-tories inside the box corresponds to the quantum proba-bility density (cid:104) ψ | ψ (cid:105) . The decay of this probability, (cid:104) ψ | ψ (cid:105) ∝ e − Γ t , (31)relates to the decay (30) of the number of classical tra-jectories inside the box. Thus, in the classical limitIm k n = − Γ n → − γ γ is found to be the largest real zero of a dynamical zetafunction Z ( s ) = (cid:89) ˜ p (1 − t ˜ p ( s )) , (33)with ˜ p the primitive periodic cycles and t ˜ p the cycleweights. For a three-dimensional system, t ˜ p ( s ) = e − l ˜ p s (cid:12)(cid:12)(cid:12) λ (1)˜ p λ (2)˜ p (cid:12)(cid:12)(cid:12) . (34)The quantities λ ( i )˜ p are the leading stability eigenvalues.A generalization to a zeta function for three dimensionsand multiple traversals r of the primitive cycle ˜ p is givenby Z ( s ) = exp − (cid:88) ˜ p (cid:88) r r e − rl ˜ p s (cid:12)(cid:12)(cid:12) λ (1)˜ p λ (2)˜ p (cid:12)(cid:12)(cid:12) r , (35)with λ ( i )˜ p the leading two stability eigenvalues of ˜ p . Thiszeta function can be cycle-expanded as described in Sec.III. Results for the escape rate γ at various configura-tions d/R are given in Table VII. V. RESULTS
The fractal Weyl law has been put to test for billiardsystems before. In [6], the 3-disk billiard has been stud-ied. To make our own results comparable to those givenin [6], we carry out a similiar discussion.
TABLE VII. Classical escape rates γ ( n )0 in order n of the cycleexpansion for various values of the configuration parameter d/R . d/R γ (1)0 γ (2)0 γ (3)0 γ (4)0 . . . . . . . . . . . . . . . . A. Defining a scale for the strip widths
For the 3-disk system discussed in [6], the strip widths C have been chosen in relation to the classical escaperate γ . For large values k → ∞ , the imaginary partof quantum resonances converges to Im k = − γ / C to˜ C := Cγ / , (36)which defines a universal scale independent of the sym-metry subspace and the ratio d/R . Similar to [6], weevaluate the fractal Weyl law for scaled strip widths˜ C ∈ [1; 1 . B. Counting resonances
We have computed spectra for various values of d/R in all symmetry subspaces. Generally, we find the bestconvergence behavior of cycle expansions for large valuesof d/R . Furthermore, the one-dimensional group repre-sentations A and A yield the largest number of con-verged resonances. The two-dimensional representation E and the three-dimensional T , T representations con-verge not as well in cycle expansion since the shadowingof individual cycles is less efficient for the weight factorsof these subspaces. For A , where all weight factors areequal to 1, the best convergence is observed.It is important to note that for the tests of the fractalWeyl law we have used only converged resonances. Forexample, the A resonances at d/R = 10 are sufficientlyconverged in the region Re k < k > − .
45 sothat counting functions N ( k ) obtained in orders 10 and11 of the cycle expansion fully agree. Thus, comparisonswith the fractal Weyl law are not affected by the orderof the cycle expansion.The fractal Weyl law (3) suggests that the countingfunctions N ( k ) obey a power law, N ( k ) ∝ k α , (37)thus, in a logarithmic plot of N ( k ), straight lines of slope α = 1 + d H are expected. A sample spectrum and cor-responding counting functions for d/R = 10 are shownin Fig. 5. This figure is generic in structure, i.e. we have0 -0.45-0.4-0.35-0.3-0.25-0.2 0 1000 2000 3000 4000 5000 I m k Re k- γ /2(a) N ( k ) k(b) C = 0.28C = 0.34C = 0.35C = 0.36C = 0.37C = 0.45const ⋅ x FIG. 5. (a) Spectrum and (b) counting functions N ( k ) for A resonances calculated from cycle expansion in order 11 forthe ratio d/R = 10. Several counting functions for differentstrips C are shown. The curves can be used as “raw data”to fit power laws. In this way, the exponent α in the fractalWeyl law can be obtained and compared to the classical cal-culations. More than 50 000 resonances have been used in theanalysis. found similar behavior of N ( k ) in other subspaces andfor other ratios d/R as well. Thus, a brief discussion ofthese features will be given in the following.We first note that the strip width C has to be suffi-ciently large, since otherwise, counting would not involvereasonably large numbers of resonances. In the spectrumshown in Fig. 5(a), we have found converged resonancesin the relevant strip ˜ C ∈ [1; 1 .
6] for Re k (cid:46) C = 0 . ⇔ ˜ C = 0 .
97, it isevident from Fig. 5(b) that counting involves only a lim-ited number of resonances. Larger strip widths involvemore resonances in the counting. However, choosing thestrip width too large, the counting may also involve res-onances that may not have been converged. Taking theasymptotic behavior of the resonances’ imaginary partsinto account, choosing rescaled strip widths in the inter-val ˜ C ∈ [1 .
0; 1 .
6] turns out to be a reasonable choice.Figure 5(b) reveals that the counting functions N ( k )deviate from power laws that led to straight lines in the α γ (a) k ∈ [0; 2000] α γ (b) k ∈ [0; 2000] FIG. 6. Exponents α obtained from least-squares fits of apower law to measured counting functions for (a) A reso-nances and (b) A resonances calculated for d/R = 6 in order13 of the cycle expansion. The power law has been fitted tothe interval k ∈ [0; 2000]. The vertical dotted line gives theclassical exponent α = 1 + d H = 1 . plot. From this observation one infers that the exponent α will clearly depend on the k -range one fits to. We follow[6] and choose the largest interval converged resonanceshave been computed for. C. Putting the fractal Weyl law to test
To provide several tests for the fractal Weyl law, wewill provide plots showing the exponents α obtained fromleast-squares fitting a power function N ( k ) ∝ k α to thecounting functions calculated from the spectra for vari-ous subspaces, fitting ranges [0; k max ], configuration pa-rameters d/R and strip widths ˜ C . We have performedleast-squares fits to match the power function (37) to themeasured N ( k ).
1. Configuration d/R = 6
For d/R = 6, we obtain the exponents shown in Fig. 6.For the A subspace, we find a very good agreement for1 α γ (a) k ∈ [0; 2000]k ∈ [0; 2500]k ∈ [0; 3000] α γ (b) k ∈ [0; 2000]k ∈ [0; 2500]k ∈ [0; 3000] FIG. 7. Exponents α obtained for d/R = 8 from least-squaresfits of a power law to measured counting functions for (a) A resonances and (b) A resonances. The power law has beenfitted to several k -intervals. The vertical dotted line gives theclassical exponent α = 1 + d H = 1 . moderate values of ˜ C < .
4. The relative error in this˜ C -interval is less than 2 percent. However, in the A subspace, all computed exponents are too large by about8 percent for the same ˜ C -interval. One possible reason isthat the k -range used for fitting is too small.
2. Configuration d/R = 8
Performing the same procedure for a configuration pa-rameter of d/R = 8, we obtain the plots shown in Fig. 7.Both plots reveal a clear tendency to obey the fractalWeyl law within a smaller error range when the range of k values used for the fit increases. However, for reasonsof convergence, longer spectra have not been used. Wenote that for ˜ C < . k ∈ [0; 3000], the error is lessthan 7 percent for the A resonances. The exponentsobtained from A resonances are larger than the expectedexponent. For ˜ C = 1 .
3, the relative error is about 15percent. α γ (a) k ∈ [0; 2000]k ∈ [0; 3000]k ∈ [0; 4000]k ∈ [0; 5000]k ∈ [0; 6000] α γ (b) k ∈ [0; 2000]k ∈ [0; 3000]k ∈ [0; 4000]k ∈ [0; 5000] FIG. 8. Exponents α obtained for d/R = 10 from least-squares fits of a power law to measured counting functionsfor (a) A resonances and (b) A resonances. Results for sev-eral k -intervals are shown. The vertical dotted line gives theclassical exponent α = 1 + d H = 1 .
3. Configuration d/R = 10
Finally, the system configuration given by d/R = 10has been studied. Results are shown in Fig. 8.The counting functions for the A subspace once moretend to give the expected exponent when the k -intervalused for the fit increases. For k ∈ [0; 6000], the error isless than 3 percent. Again, the A spectra yield expo-nents that are too large. Possibly the k -range investi-gated here is not large enough to exhibit the asymptoticbehavior clearly. D. Discussion
While the classical calculations for the fractal dimen-sion d H are accurate to at least four significant digits,the agreement of the exponents α is at best 2 to 7 per-cent for the A spectra. However, we note that for thetwo-dimensional three-disk billiard, the errors in the ex-ponents are about 5 to 10 percent [6]. Therefore, weconclude that the fractal Weyl law for the four-spherescattering system is confirmed with roughly the same ac-2curacy as for the three-disk billiard.A very large k -range seems to be necessary for a properinvestigation. This tendency is also visible for the A spectra. The exponents obtained from the A spectraare too large. However, using larger k -intervals, the ex-ponents seem to approach the correct value for large stripwidths ˜ C . Possibly, if larger spectra were available, theexpected exponents could be obtained. Unfortunately,we are limited by the convergence of the cycle expansionswe use. The higher-dimensional symmetry subspaces E , T and T could not be used to put the fractal Weyl totest law since the spectra did not contain enough con-verged resonances. VI. SUMMARY AND OUTLOOK
This paper provides a test of the fractal Weyl law fora three-dimensional scattering system. The four-spherebilliard was investigated both classically and quantummechanically.In Sec. II, we have developed a fast and very precisemethod to gauge the repeller. We found estimates for theHausdorff dimension d H with a relative accuracy of 10 − .Although the algorithm is based on strong assumptions, it works over a wider range of the configuration parame-ter d/R than existing methods.In Sec. III, we have discussed the methods of semiclas-sical quantization. We have applied the method of cycleexpansion to the four-sphere billiard. Furthermore, forthe first time, the method of symmetry decompositionwas demonstrated for the Gutzwiller-Voros zeta functionof the system.We have given results in Sec. V. We have providedtests of the fractal Weyl law for various configurations ofthe system. Although we have found the counting func-tions N ( k ) to deviate from power functions, we couldconfirm the fractal Weyl law for the A resonances of thefour-sphere scattering within a small error range. Forthose spectra we did not find a convincing agreementof calculated level numbers N ( k ) with the prediction N ( k ) ∝ k d H for, there is hope that larger spectra wouldlead to the expected exponent. We also assume that thedeviations from pure power laws are due to the fact thatthe energy range under consideration is too small.As an outlook, the physical origin of the modulationsin the counting functions N ( k ) will have to be investi-gated. Moreover, it is desirable to study further three-dimensional scattering systems to find out to what extentthe results found for the four-sphere billiard carry overand are generic. [1] H. Weyl, Mathematische Annalen , 441 (1912).[2] R. Courant, Mathematische Zeitschrift , 1 (1920).[3] H. P. Baltes and E. R. Hilf, Spectra of Finite Systems (B.I. Wissenschaftsverlag, 1976).[4] J. Sj¨ostrand, Duke Math. J. , 1 (1990).[5] M. Zworski, Invent. Math. (1999).[6] W. T. Lu, S. Sridhar, and M. Zworski, Phys. Rev. Lett. , 154101 (2003).[7] K. K. Lin and M. Zworski, Chem. Phys. Lett. , 201(2002).[8] J. Wiersig and J. Main, Phys. Rev. E , 036205 (2008).[9] J. A. Ramilowski, S. D. Prado, F. Borondo, and D. Far-relly, Phys. Rev. E , 055201 (2009).[10] D. L. Shepelyansky, Phys. Rev. E , 015202 (2008).[11] M. Kopp and H. Schomerus, Phys. Rev. E , 026208(Feb 2010).[12] D. Sweet, E. Ott, and J. A. Yorke, Nature , 315(1999).[13] D. Sweet, B. W. Zeff, E. Ott, and D. P. Lathrop, PhysicaD , 207 (2001).[14] Q. Chen, M. Ding, and E. Ott, Phys. Lett. A , 93(1990).[15] A. E. Motter and P. S. Letelier, Phys. Lett. A , 18(2000).[16] K. J. Falconer, The geometry of fractal sets , CambridgeTracts in Mathematics, Vol. 85 (Cambridge UniversityPress, Cambridge, 1985).[17] G. M. Zaslavsky,
Hamiltonian Chaos & Fractional Dy-namics (Oxford University Press, Oxford, 2005).[18] M. Henseler, A. Wirzba, and T. Guhr, Annals of Physics , 286 (1997).[19] P. Gaspard and S. A. Rice, J. Chem. Phys. , 2255(1989).[20] J. Main, E. Atilgan, H. S. Taylor, and G. Wunner, Phys.Rev. E , 056227 (2004).[21] M. C. Gutzwiller, Chaos in Classical and Quantum Me-chanics (Springer, New York, Berlin, Heidelberg, 1990).[22] A. Voros, J. Phys. A , 685 (1988).[23] M. Brack and R. K. Bhaduri, Semiclassical Physics (Addison-Wesley, Reading, MA, 2003).[24] P. Cvitanovi´c, Phys. Rev. Lett. , 2729 (1988).[25] P. Cvitanovi´c and B. Eckhardt, Phys. Rev. Lett. , 823(1989).[26] R. Artuso, E. Aurell, and P. Cvitanovi´c, Nonlinearity ,325 (1990).[27] A. Wirzba, Phys. Rep. , 1 (1999).[28] H. Primack and U. Smilansky, Phys. Rep. , 1 (2000).[29] M. Sieber, Nonlinearity , 1607 (1998).[30] P. Cvitanovi´c and B. Eckhardt, Nonlinearity , 277(1993).[31] M. R. Wall and D. Neuhauser, J. Chem. Phys. , 8011(1995).[32] J. Main, Phys. Rep. , 233 (1999).[33] Dˇz. Belki´c, P. A. Dando, J. Main, and H. S. Taylor, J.Chem. Phys. , 6542 (2000).[34] P. Gaspard and S. A. Rice, J. Chem. Phys. , 2225(1989).[35] K. Pance, W. Lu, and S. Sridhar, Phys. Rev. Lett.85