Periodic Orbit Quantization: How to Make Semiclassical Trace Formulae Convergent ∗ J¨org Main and Gu¨nter Wunner Institut fu¨r Theoretische Physik I, Universita¨t Stuttgart, D-70550 Stuttgart, Germany (February 8, 2008) 1 0 Abstract 0 2 n a Periodic orbit quantization requires an analytic continuation of non-con- J vergent semiclassical trace formulae. We propose two different methods for 2 semiclassical quantization. The first method is based upon the harmonic in- ] version of semiclassical recurrence functions. A band-limited periodic orbit D signalisobtainedbyanalyticalfrequencywindowingoftheperiodicorbitsum. C The frequencies of the periodic orbit signal are the semiclassical eigenvalues, . n and are determined by either linear predictor, Pad´e approximant, or signal i l n diagonalization. The second method is based upon the direct application of [ the Pad´e approximant to the periodic orbit sum. The Pad´e approximant al- 1 lows the resummation of the, typically exponentially, divergent periodic orbit v terms. Both techniques do not depend on the existence of a symbolic dy- 7 0 namics, and can be applied to bound as well as to open systems. Numerical 0 results are presented for two different systems with chaotic and regular clas- 1 0 sical dynamics, viz. the three-disk scattering system and the circle billiard. 1 0 PACS numbers: 05.45. a, 03.65.Sq / − n i l n : v i X r a Typeset using REVTEX ∗Contribution to “Festschrift in honor of Martin Gutzwiller”, eds. A. Inomata et al., to be pub- lished in Foundations of Physics. 1 I. INTRODUCTION Semiclassical periodic orbit quantization is a nontrivial problem for the reason that Gutzwiller’s trace formula [1,2] for chaotic systems and the Berry-Tabor formula [3] for integrable systems do not usually converge in those regions where the physical eigenenergies or resonances arelocated. Various techniques have beendeveloped to circumvent theconver- gence problem of periodic orbit theory. Examples are the cycle expansion technique [4], the Riemann-Siegel type formula and pseudo-orbit expansions [5], surface of section techniques [6], and a quantization rule based on a semiclassical approximation to the spectral staircase [7]. These specific techniques have proven to be very efficient for individual systems with special properties, e.g., the cycle expansion for hyperbolic systems with an existing sym- bolic dynamics. The other methods mentioned have been used for the calculation of bound spectra of specific systems. Recently, an alternative method based upon filter-diagonalization (FD) has been intro- duced for the analytic continuation of the semiclassical trace formula [8,9]. The FD method requires knowledge of the periodic orbits up to a given maximum period (classical action), which depends on the mean density of states. The semiclassical eigenenergies or resonances are obtained by harmonic inversion of the periodic orbit recurrence signal. The FD method can be generally applied to both open and bound systems and has also proven to be a powerful tool, e.g., for the calculation of semiclassical transition matrix elements [10] and the quantization of systems with mixed regular-chaotic phase space [11]. For a review on periodic orbit quantization by harmonic inversion see [12]. Inthispaperwepresenttwodifferenttechniquestomakesemiclassicalperiodicorbitsums convergent. The first method is an advanced version of harmonic inversion adapted to the special structure of periodic orbit signals given as sums of δ functions [13]. The semiclassical signal, in action or time, corresponds to a “spectrum” or response in the frequency domain that is composed of a huge, in principle infinite, number of frequencies. To extract these frequencies andtheir corresponding amplitudes is anontrivial task. Inprevious work [8,9,12] the periodic orbit signal has been harmonically inverted by means of FD [14–16] which is designed for the analysis of time signals given on an equidistant grid. The periodic orbit recurrence signalisrepresented asasumover usually unevenly spaced δ functions. Asmooth signal, from which evenly spaced values can be read off, is obtained by a convolution of this sum with, e.g., a narrow Gaussian function. The disadvantages of this approach are twofold. Firstly, FD acts on this signal more or less like a “black box” and, as such, does not lend itself to a detailed understanding of semiclassical periodic orbit quantization. Secondly, the smoothedsemiclassical signalusually consists ofahugenumber ofdatapoints. Thehandling of such large data sets, together with the smoothing, may lead to significant numerical errors in results for the semiclassical eigenenergies and resonances. Here, we propose alternative techniques for the harmonic inversion of the periodic orbit recurrence signal that avoid these problems. In a first step we create a shortened signal which is constructed from the original signal and designed to be correct only in a window, i.e., a short frequency range of the total band width. Because the original signal is given as a periodic orbit sum of δ functions, this “filtering” can be performed analytically resulting in a band-limited periodic orbit signal with a relatively small number of equidistant grid points. In a second step the frequencies and amplitudes of the band-limited signal are determined from a set of nonlinear equations. 2 To solve the nonlinear system, we introduce three different processing methods, viz. linear predictor (LP), Pad´e approximant (PA), and signal diagonalization (SD). It is important to notethattheseprocessing methodswouldnothaveyielded numerically stablesolutionsifthe signal had not first been band-limited by the windowing (filtering) procedure. Furthermore, thisseparationoftheharmonicinversion procedureintovariousstepsmayelucidate aclearer picture of the periodic orbit quantization method itself, and even provides more accurate results than previous calculations [9,12] using FD. The second method is the direct application of the Pad´e approximant to slowly con- vergent and/or divergent periodic orbit sums [17]. In the former or the latter case, the PA either significantly increases the convergence rate, or analytically continues theexponentially divergent series, respectively. The PA is especially robust for resumming diverging series in many applications in mathematics and theoretical physics [18]. An important example is the summation of the divergent Rayleigh-Schr¨odinger quantum mechanical perturbation series, e.g., for atoms in electric [19] and magnetic [20] fields. In periodic orbit theory the PA has been applied to cycle-expanded Euler products and dynamical zeta functions [21]. It should be noted that the Pad´e approximant is applied in both methods in a completely different context, namely, in the first method as a tool for signal processing [13,22], and in the second for the direct summation of the periodic orbit terms in the semiclassical trace formulae. In Sec. II we present our first method to make semiclassical trace formulae convergent. After briefly reviewing the general idea of periodic orbit quantization by harmonic inversion in Sec. IIA we construct, in Sec. IIB, the band-limited periodic orbit signal which is ana- lyzed, in Sec. IIC, with the help of either LP, PA, or SD. In Sec. III we introduce our second method for semiclassical quantization, viz. the direct application of the Pad´e approximant to the periodic orbit sum. In Sec. IV we present and compare results for the three-disk repeller and the circle billiard as physical examples of relevance. A few concluding remarks are given in Sec. V. II. HARMONIC INVERSION OF PERIODIC ORBIT SIGNALS A. General remarks In order to understand what follows, a brief recapitulation of the basic ideas of periodic orbit quantization by harmonic inversion may be useful. For further details see [12]. Following Gutzwiller [1,2] one can write the semiclassical response function for chaotic systems in the form gsc(E) = gsc(E)+ eiSpo , (1) 0 Apo po X where gsc(E) is a smooth function and S and are the classical actions and weights 0 po Apo (including phase information given by the Maslov index) of the periodic orbit (po) contri- butions. Equation (1) is also valid for integrable systems when the periodic orbit quantities are calculated not with Gutzwiller’s trace formula, but with the Berry-Tabor formula [3] for periodic orbits on rational tori. The eigenenergies and resonances are the poles of the response function. Unfortunately, the semiclassical approximation (1) does not converge in 3 the region of the poles, and hence one is faced with the problem of the analytic continuation of gsc(E) to this region. As in previous work [8,9,12], we will make the (weak) assumption that the classical system has a scaling property, i.e., the shape of periodic orbits is assumed not to depend on a scaling parameter, w, and the classical action scales as S = ws . (2) po po In scaling systems, the fluctuating part of the semiclassical response function, gsc(w) = eiwspo , (3) po A po X can be Fourier transformed readily to yield the semiclassical trace of the propagator 1 +∞ Csc(s) = gsc(w)e−iswdw = δ(s s ) . (4) po po 2π A − Z−∞ po X The signalCsc(s) hasδ spikes at thepositions ofthe classical periods(scaled actions) s = s po ofperiodicorbitsandwithpeakheights(recurrencestrengths) ,i.e.,Csc(s)isGutzwiller’s po A periodic orbit recurrence function. Consider now the quantum mechanical counterparts of gsc(w) and Csc(s) taken as the sums over the poles w of the Green’s function, k d gqm(w) = k , (5) w w +iǫ Xk − k 1 +∞ Cqm(s) = gqm(w)e−iswdw = i d e−iwks , (6) k 2π − Z−∞ k X with d being the residues associated with the eigenvalues. In the case under study, i.e., k density of states spectra, the d are the multiplicities of eigenvalues and are equal to 1 for k non-degenerate states. Semiclassical eigenenergies w and residues d can now, in principle, k k be obtained by adjusting the semiclassical signal, Eq. (4), to the functional form of the quantum signal, Eq. (6), with the w ,d being free, in general complex, frequencies and k k { } amplitudes. This scheme is known as “harmonic inversion”. The numerical procedure of harmonic inversion is a nontrivial task, especially if the number of frequencies in the signal is large (e.g., more than a thousand), or even infinite as is usually the case for periodic orbit quantization. Note that the conventional way to perform the spectral analysis, i.e., the Fourier transform of Eq. (4) will bring us back to analyzing the non-convergent response function gsc(w) in Eq. (3). The periodic orbit signal (4) can be harmonically inverted by application of FD [14–16], which allows one to calculate a finite and relatively small set of frequencies and amplitudes in a given frequency window. The usual implementation of FD requires knowledge of the signal on an equidistant grid. The signal (4) is not a continuous function. However, a smooth signal can be obtained by a convolution of Csc(s) with, e.g., a Gaussian function, 1 Csc(s) = e−(s−spo)2/2σ2 . (7) σ √2πσ Apo po X 4 (σ) As can easily be seen, the convolution results in a damping of the amplitudes, d d = k → k d exp( w2σ2/2). The width σ of the Gaussian function should be chosen sufficiently small k − k to avoid an overly strong damping of amplitudes. To properly sample each Gaussian a dense grid with steps ∆s σ/3 is required. Therefore, the signal (7) analyzed by FD ≈ usually consists of a large number of data points. The numerical treatment of this large data set may suffer from rounding errors and loss of accuracy. Additionally, the “black box” type procedure of harmonic inversion by FD, which intertwines windowing and processing, does not provide any opportunity to gain a deeper understanding of semiclassical periodic orbit quantization. It is therefore desirable to separate the harmonic inversion procedure into two sequential steps: Firstly, the filtering procedure that does not require smoothing and, secondly, a procedure for extracting the frequencies and amplitudes. In Sec. IIB we will construct, by analytic filtering, a band-limited signal which consists of a relatively small number of frequencies. In Sec. IIC we will present methods to extract the frequencies and amplitudes of such band-limited signals. B. Construction of band-limited signals by analytical filtering In general, a frequency filter can beapplied to a given signal by application of the Fourier transform [22–24]. The signal is transformed to the frequency domain, e.g., by application of the fast Fourier transform (FFT) method. The transformed signal is multiplied with a frequency filter function f(w) localized around a central frequency, w . The frequency 0 filter f(w) can be rather general, typical examples are a rectangular window or a Gaussian function. Thefilteredsignalisthenshiftedby w andtransformedbacktothetimedomain 0 − by a second application of FFT. The filtered signal consists of a significantly reduced set of frequencies, and therefore a reduced set of grid points is sufficient for the analysis of the signal. This technique is known as “beam spacing” [23] or “decimation” [22,24] of signals. The special form of the periodic orbit signal (4) as a sum of δ functions allows for an even simpler procedure, viz. analytical filtering. In the following we will apply a rectangular filter, i.e., f(w) = 1 for frequencies w [w ∆w,w + ∆w], and f(w) = 0 outside the 0 0 ∈ − window. The generalization to other types of frequency filters is straightforward. Starting fromthesemiclassical responsefunction(spectrum)gsc(w)inEq.(3),whichisitselfaFourier transform of the “signal” (4), and using a rectangular window we obtain, after evaluating the “second” Fourier transform, the band-limited (bl) periodic orbit signal, 1 w0+∆w Csc(s) = gsc(w)e−is(w−w0)dw bl 2π Zw0−∆w 1 w0+∆w = eisw0−i(s−spo)wdw po 2π A po Zw0−∆w X sin[(s s )∆w] = − po eispow0 . (8) po A π(s s ) Xpo − po The introduction of w into the arguments of the exponential functions in (8) causes a shift 0 of frequencies by w in the frequency domain. Note that Csc(s) is a smooth function and − 0 bl can be easily evaluated on an arbitrary grid of points s < s provided the periodic orbit n max data are known for the set of orbits with classical action s < s . po max 5 Applying now the same filter as used for the semiclassical periodic orbit signal to the quantum one, we obtain the band-limited quantum signal 1 w0+∆w Cqm(s) = gqm(w)e−is(w−w0)dw bl 2π Zw0−∆w K = i d e−i(wk−w0)s , w w < ∆w . (9) k k 0 − | − | k=1 X In contrast to the signal Cqm(s) in Eq. (6), the band-limited quantum signal consists of a finite number of frequencies w , k = 1,...,K, where in practical applications K can be of k the order of (50-200) for an appropriately chosen frequency window, ∆w. The problem of ∼ adjustingtheband-limitedsemiclassical signalinEq.(8)toitsquantummechanicalanalogue in Eq. (9) can now be written as a set of 2K nonlinear equations K Csc(nτ) c = i d e−iwk′nτ , n = 0,1,...,2K 1 , (10) bl ≡ n − k − k=1 X forthe2K unknown variables, viz. theshiftedfrequencies, w′ w w ,andamplitudes, d . k ≡ k− 0 k The band-limited signal now becomes “short”as it can be evaluated on an equidistant grid, s = nτ, with relatively large step width τ π/∆w. It is important to note that the number ≡ of signal points c in Eq. (10) is usually much smaller than a reasonable discretization of the n signal Csc(s) in Eq. (7), which is the starting point for harmonic inversion by FD. Therefore, σ the discrete signal points c Csc(nτ) are called the “band-limited” periodic orbit signal. n ≡ bl Methods to solve the nonlinear system, Eq. (10), are discussed in Sec. IIC below. It should also be noted that the analytical filtering in Eq. (8) is not restricted to periodic orbit signals, but can be applied, in general, to any signal given as a sum of δ functions. An example is the high resolution analysis of quantum spectra [12,25,26], where the density of states is ̺(E) = δ(E E ). n − n P C. Harmonic inversion of band-limited signals In this section we wish to solve the nonlinear set of equations K c = d zn , n = 0,1,...,2K 1 , (11) n k k − k=1 X where z exp( iw′τ) and d are, generally complex, variational parameters. For nota- k ≡ − k k tionalsimplicity we have absorbed the factor of i onthe right-hand side of Eq. (10) into the − d ’s with the understanding that this should be corrected for at the end of the calculation. k We assume that the number of frequencies in the signal is relatively small (K 50 to 200). ∼ Although the system of nonlinear equations is, in general, still ill-conditioned, frequency filtering reduces the number of signal points, and hence the number of equations. Several numerical techniques, that otherwise would be numerically unstable, can now be applied successfully. In the following we employ three different methods, viz. linear predictor (LP), Pad´e approximant (PA), and signal diagonalization (SD). 6 1. Linear Predictor The problem of solving Eq. (11) has already been addressed in the 18th century by Baron de Prony [27], who converted the nonlinear set of equations (11) to a linear algebra problem. Today this method is known as linear predictor (LP). Our method strictly applies theprocedure ofLPexcept with oneessential difference; theoriginalsignalCsc(s) isreplaced with its band-limited counterpart c Csc(nτ). n ≡ bl Equation (11) can be written in matrix form for the signal points c to c , n+1 n+K c zn+1 zn+1 d  n...+1  =  1... ·..·.· K...  ...1  . (12) cn+K   z1n+K ··· zKn+K  dK       From the matrix representation (12) it follows that zn+1 zn+1 −1 c cn = (z1n···zKn) 1... ·..·.· K...   n...+1  = K akcn+k , (13) z1n+K ··· zKn+K   cn+K  kX=1     which means that every signal point c can be “predicted” by a linear combination of the n K subsequent points with a fixed set of coefficients a , k = 1,...,K. The first step of k the LP method is to calculate these coefficients. Writing Eq. (13) in matrix form with n = 0,...,K 1, we obtain the coefficients a as the solution of the linear set of equations, k − c c a c 1 K 1 0 ···  ... ... ...  ...  =  ...  . (14)  cK c2K−1 aK  cK−1   ···          The second step consists in determining the parameters z in Eq. (11). Using Eqs. (13) and k (11) we obtain K K K c = a c = a d zn+k , (15) n k n+k k l l k=1 l=1k=1 X XX and thus K K a zn+l zn d = 0 . (16) " l k − k# k k=1 l=1 X X Equation (16) is satisfied for arbitrary sets of amplitudes d when z is a zero of the poly- k k nomial K a zl 1 = 0 . (17) l − l=1 X The parameters z = exp( iw′τ) and thus the frequencies k − k 7 i w′ = log(z ) (18) k τ k arethereforeobtainedbysearching forthezerosofthepolynomialinEq.(17). Notethatthis is the only nonlinear step of the algorithm, and numerical routines for finding the roots of polynomials are well established. In the third and final step, the amplitudes d are obtained k from the linear set of equations K c = d zn , n = 0,...,K 1 . (19) n k k − k=1 X To summarize, the LP method reduces the nonlinear set of equations (11) for the varia- tional parameters z ,d to two well-known problems, i.e., the solution of two linear sets k k { } of equations (14) and (19) and the root search of a polynomial, Eq. (17), which is a nonlinear but familiar problem. The matrices in Eqs. (14) and (19) are a Toeplitz and Vandermonde matrix, respectively, and special algorithms are known for the fast solution of such linear systems [28]. However, when the matrices are ill-conditioned, conventional LU decomposi- tion of the matrices is numerically more stable, and, furthermore, an iterative improvement of the solution can significantly reduce errors arising from numerical rounding. The roots of polynomials can be found, in principle, by application of Laguerre’s method [28]. However, it turns out that an alternative method, i.e., the diagonalization of the Hessenberg matrix aK−1 aK−2 a1 a0 − aK − aK ··· −aK −aK 1 0 0 0   ··· A = 0 1 0 0 , (20)  ···   ... ... ... ... ...       0 0 1 0   ···    for which the characteristic polynomial P(z) = det[A zI] = 0 is given by Eq. (17) (with − a = 1), is a numerically more robust technique for finding the roots of high degree 0 > − (K 60) polynomials [28]. ∼ 2. Pad´e Approximant As an alternative method for solving the nonlinear system (11) we now propose to apply the method of Pad´e approximants (PA) to our band-limited signal c . Let us assume for n the moment that the signal points c are known up to infinity, n = 0,1,... . Interpreting n ∞ the c ’s as the coefficients of a Maclaurin series in the variable z−1, we can then define the n function g(z) = ∞ c z−n. With Eq. (11) and the sum rule for geometric series we obtain n=0 n P ∞ K ∞ K zd P (z) g(z) c z−n = d (z /z)n = k K . (21) n k k ≡ z z ≡ Q (z) nX=0 kX=1 nX=0 kX=1 − k K The right-hand side of Eq. (21) is a rational function with polynomials of degree K in the numerator and denominator. Evidently, the parameters z = exp( iw′τ) are the poles of k − k g(z), i.e., the zeros of the polynomial Q (z). The parameters d are calculated via the K k residues of the last two terms of (21). We obtain 8 P (z ) K k d = , (22) k z Q′ (z ) k K k with the prime indicating the derivative d/dz. Of course, the assumption that the coeffi- cients c are known up to infinity is not fulfilled and, therefore, the sum on the left-hand n side of Eq. (21) cannot be evaluated in practice. However, the convergence of the sum can be accelerated by application of PA. Indeed, with PA, knowledge of 2K signal points c ,...,c is sufficient for the calculation of the coefficients of the two polynomials 0 2K−1 K K P (z) = b zk and Q (z) = a zk 1 . (23) K k K k − k=1 k=1 X X The coefficients a , k = 1,...,K are obtained as solutions of the linear set of equations k K c = a c , n = 0,...,K 1 , n k n+k − k=1 X which is identical to Eqs. (13) and (14) for LP. Once the a’s are known, the coefficients b k are given by the explicit formula K−k b = a c , k = 1,...,K . (24) k k+m m m=0 X It should be noted that the different derivations of LP and PA yield the same polynomial whose zeros are the z parameters, i.e., the z calculated with both methods exactly agree. k k However, LP and PA do differ in the way the amplitudes, d , are calculated. It is also k important to note that PA is applied here as a method for signal processing, i.e., in a different context to that in Sec. III, where the Pad´e approximant is used for the direct summation of the periodic orbit terms in semiclassical trace formulae. 3. Signal Diagonalization In Refs. [14,16] it has been shown how the problem of solving the nonlinear set of equa- tions (11) can be recast in the form of the generalized eigenvalue problem, UB = z SB . (25) k k k The elements of the K K operator matrix U and overlap matrix S depend trivially upon × the c ’s [16]: n U = c ; S = c ; i,j = 0,...,K 1 . (26) ij i+j+1 ij i+j − Note that the operator matrix U is the same as in the linear system (14), i.e., the matrix form of Eq. (13) of LP. The matrices U and S in Eq. (25) are complex symmetric (i.e., B non-Hermitian), and the eigenvectors are orthogonal with respect to the overlap matrix k S, (Bk S Bk′) = Nkδkk′ , (27) | | 9 where the brackets define a complex symmetric inner product (a b) = (b a), i.e., no complex | | conjugationofeitheraorb.TheoverlapmatrixSisnotusuallypositivedefiniteandtherefore the N ’s are, in general complex, normalization parameters. An eigenvector B cannot be k k normalized for N = 0. The amplitudes d in Eq. (11) are obtained from the eigenvectors k k B via k 2 1 K−1 d = c B . (28) k n k,n Nk "n=0 # X The parameters z in Eq. (11) are given as the eigenvalues of the generalized eigenvalue k problem (25), and are simply related to the frequencies w′ in Eq. (10) via z = exp( iw′τ). k k − k The three methods introduced above (LP, PA and SD) look technically quite different. With LP the coefficients of the characteristic polynomial (17) and the amplitudes d are obtained k by solving two linear sets of equations (14) and (19). Note that the complete set of zeros z k of Eq. (17) is required to solve for the d in Eq. (19). The PA method is even simpler, as only k one linear system, Eq. (14), has to be solved to determine the coefficients of the rational function P (z)/Q (z). Finding the zeros of Eq. (17) provides knowledge about selected K K parameters z , and allows one to calculate the corresponding amplitudes d via Eq. (22). k k The SD method requires the most numerical effort, because the solution of the generalized eigenvalue problem (25) for both the eigenvalues z and eigenvectors B is needed. k k Itisimportanttonotethatthethreetechniques, inspiteoftheirdifferent derivations, are mathematically equivalent and provide the same results for the parameters z ,d , when k k { } the following two conditions are fulfilled: the nonlinear set of equations (11) has a unique solution, when, firstly, the matrices U and S in Eq. (26) have a non-vanishing determinant (detU = 0, detS = 0), and, secondly, the parameters zk are non-degenerate (zk = zk′ for 6 6 6 k = k′). These conditions guaranteetheexistence ofa unique solutionofthelinear equations 6 (14) and (19), the non-singularity of the generalized eigenvalue problem (25), and the non- vanishing of both the derivatives Q′ (z ) in Eq. (22) and the normalization constants N K k k in Eqs. (27) and (28). Equation (11) usually has no solution in the case of degenerate z k parameters, however, degeneracies can be handled with a generalization of the ansatz (11) and modified equations for the calculation of the parameters. Here, we will not further discuss this special case. While the parameters z in Eq. (11) are usually unique, the calculation of the frequencies k w′ via Eq. (18) is not unique, because of the multivalued property of the complex logarithm. k To obtain the “correct” frequencies it is necessary to appropriately adjust the range ∆w of the frequency filter and the step width τ of the band-limited signal (10). From our numericalexperience werecommend thefollowingprocedure. Themostconvenient approach is to choose first the center w of the frequency window and the number K of frequencies 0 within that window. Note that K determines the dimension of the linear systems, and hence the degree of the polynomials which have to be handled numerically, and is therefore directly related to the computational effort required. Frequency windows are selected to be sufficiently narrow to yield values for the rank between K 50 and K 200. The step ≈ ≈ width for the band-limited signal should be chosen as s max τ = , (29) 2K 10