{'fraction_non_alphanumeric': 0.0448592922819727, 'fraction_numerical': 0.01870795685228675, 'mean_word_length': 4.281479924322052, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fractal Weyl law for chaotic microcavities: Fresnel's laws imply multifractal scattering

We demonstrate that the harmonic inversion technique is a powerful tool to analyze the spectral properties of optical microcavities. As an interesting example we study the statistical properties of complex frequencies of the fully chaotic microstadium. We show that the conjectured fractal Weyl law for open chaotic systems [W. T. Lu, S. Sridhar, and M. Zworski, Phys. Rev. Lett. 91, 154101 (2003)] is valid for dielectric microcavities only if the concept of the chaotic repeller is extended to a multifractal by incorporating Fresnel's laws.


I. INTRODUCTION
Optical microcavities are expected to be key components in future photonic applications, such as ultralow threshold lasers [1,2], single-photon emitters [3,4], and correlated photon-pair emitters [5]. Microdisk cavities with non-circular cross-sectional shape have recently attracted considerable interest in the quantum chaos community since the internal ray dynamics is non-integrable [6,7]. In fact, microdisks can be considered as open billiards. In a billiard a point-like particle moves freely in a two-dimensional domain with elastic reflections at the boundary [8]. Depending on the shape of the boundary the system shows a variety of dynamical behaviors ranging from integrable to fully chaotic [9]. Light rays in a microcavity behave similarly as they are totally reflected at the boundary as long as the angle of incidence χ (measured from the normal) is larger than the critical angle for total internal reflection χ c = arcsin (1/n), where n is the index of refraction. If, however, χ < χ c then the light ray escapes refractively according to Snell's and Fresnel's laws.
An important topic of quantum chaos is the analysis of the statistical properties of energy levels in quantum systems whose classical dynamics is fully or partially chaotic [10]. In the last years, the focus has shifted from closed to open systems. Quantum eigenenergies of open systems (resonance frequencies in the case of microcavities) are complex valued with the imaginary part being related to the lifetime of the state. Of particular interest is the fractal Weyl law for open chaotic systems [11,12,13,14,15,16,17]. This is an extension of the well-known Weyl's formula for bounded systems which states that the number of levels with wave number k n ≤ k is asymptotically N (k) ∼ k 2 for the particular case of a two-dimensional system which scales with the energy such as quantum billiards. The fractal Weyl law for an open chaotic system (which again scales with the energy) having complex wave numbers k n can be written as N (k) = {k n : Im(k n ) > −C, Re(k n ) ≤ k} ∼ k α . (1) The cutoff constant C > 0 removes fast decaying states. It is conjectured that the non-integer exponent is where D is the fractal dimension of the chaotic repeller of the open chaotic system [14]; d = D − 1 is the dimension of the repeller in a properly chosen surface of section [18]. The chaotic repeller is the set of points in phase space that never leave the system both in forward and in backward time dynamics, see, e.g., Ref. [18]. An essential prerequisite for a statistical analysis is a sufficient amount of resonance data. Resonance spectra can either be measured experimentally or computed by exact quantum or semiclassical methods. For some systems, e.g., atoms in external electric and magnetic fields the complex S-matrix poles can be obtained directly by diagonalization of non-Hermitian matrices using complex dilated basis sets [19,20]. However, in many theoretical calculations of, e.g., open billiards [21,22,23,24], optical microcavities [25], or in the cycle-expanded Gutzwiller-Voros ζ-function [26] the resonances cannot be obtained directly by matrix diagonalization but rather by a numerically expensive root search where the matrix M(E) depends on the complex valued energy E (or on the complex wave number k for billiards and microcavities). The calculation of a real valued spectrum σ(E) as the superposition of Lorentzian line shapes is much less expensive, and it would be a great advantage to be able to extract the complex S-matrix poles from a real spectrum σ(E). In principle, this can be achieved by fitting the real valued spectrum to a sum of Lorentzians. Unfortunately, the established methods for fitting Lorentzians are restricted to isolated or weakly overlapping resonances, whereas typical open systems in physically interesting regions are characterized by strongly overlapping resonances. In this paper we will apply the harmonic inversion method to obtain the resonances of open dielectric microcavities. We will show that the conjectured fractal Weyl law for open chaotic systems is valid for microcavities only if the concept of the chaotic repeller is extended to a multifractal by incorporating Fresnel's laws. The paper is organized as follows. In Sec. II we present the numerical technique for the extraction of the resonance poles from a real spectrum σ(E). In Sec. III we introduce the dielectric microstadium. The high-index and low-index microstadium are investigated in detail in Secs. IV and V. Conclusion are given in Sec. VI.

II. THE HARMONIC INVERSION TECHNIQUE
Here, we will apply methods for high-resolution signal processing to extract the resonance positions, widths, and amplitudes from a real spectrum σ(E). In a first step the spectrum is Fourier transformed with appropriate frequency filters to obtain a band-limited time signal. This time signal is analyzed, in the second step, by a high-resolution harmonic inversion method. The harmonic inversion method for signal processing has been introduced by Wall and Neuhauser [27] and then refined and improved by Mandelshtam and Taylor [28,29]. It has proved to be a powerful tool for both the highresolution analysis of quantum spectra [30] and the semiclassical quantization of non-integrable systems [31,32,33,34]. For a review on the use of the harmonic inversion method in semiclassical physics see Ref. [35]. Technically, the harmonic inversion problem can be recast as a set of nonlinear equations, which can be solved by either linear predictor, Padé approximant, or direct signal diagonalization [36]. Although details of the procedure have already been published we provide here, for the convenience of the reader and to make the paper self-contained, a brief description of the harmonic inversion method extended to the extraction of resonances from experimental or theoretical spectra given as superpositions of Lorentzians, i.e., spectra of the form where the {Ω k , Γ k , A k } are the frequency positions, widths, and amplitudes of the resonances. The frequency spectrum (4) can be interpreted as the Fourier transform of a time signal with complex frequencies ω k = Ω k − i 2 Γ k and complex values d k simply related to the amplitudes A k . The challenge is to determine very broad resonances deep in the complex plane and to resolve individual resonance poles in regions with strongly overlapping resonances. Basically, the harmonic inversion algorithm is split into the following two steps.
In a first step, a frequency window [ω 0 −∆ω/2, ω 0 +∆ω/2] is chosen. A band-limited signal C bl (t) with a finite number of about 50 to 200 frequencies is obtained by the application of a windowed discrete Fourier transform to the frequency spectrum (4). The band-limited signal is evaluated on an equidistant grid with time step τ = 2π/∆ω and reads with n = 0, 1, . . . , N bl − 1. The limits j 1 and j 2 of the sum in Eq. (6) are taken so that the equidistant grid points ω j of the digitized spectrum (4) cover the selected frequency filter The frequencies in the exponent of Eq. (6) are shifted by −ω 0 to reduce the phase oscillations of the band-limited signal. The number N bl of data points c n of the band-limited signal is restricted by the number of grid points in the frequency window, i.e., N bl ≤ j 2 − j 1 + 1. We choose N bl = 2K, where K is an upper bound of the number of resonances in the window (see below). In a second step, the (shifted) frequencies ω ′ k = ω k − ω 0 of the band-limited signal (6) are obtained by solving the nonlinear set of equations where z k ≡ exp(−iω ′ k τ ) and d k are in general complex variational parameters. As the number of frequencies in the bandlimited signal is relatively small (K ∼ 50 − 200), several variants of Prony's method [37], which otherwise would be numerically unstable, can now be applied successfully, e.g., linear predictor, Padé approximant, or direct signal diagonalization [36]. Here, we briefly elaborate on the Padé approximant.
Let us assume for the moment that the signal points c n are known up to infinity, n = 0, 1, . . . ∞. Interpreting the c n 's as the coefficients of a Maclaurin series in the variable z −1 , we can define the function g(z) = ∞ n=0 c n z −n . With Eq. (7) and the sum rule for geometric series we obtain The right-hand side of Eq. (8) is a rational function with polynomials of degree K in the numerator and denominator. Evidently, the parameters z k = exp (−iω ′ k τ ) are the poles of g(z), i.e., the zeros of the polynomial Q K (z). The parameters d k are calculated via the residues of the last two terms of Eq. (8). Application of the residue calculus yields with the prime indicating the derivative d/dz. Of course, the assumption that the coefficients c n are known up to infinity is not fulfilled and, therefore, the sum over all c n z −n in Eq. (8) cannot be evaluated in practice. However, the convergence of the sum can be accelerated by use of the Padé approximant. Indeed, knowledge of 2K signal points c 0 , . . . , c 2K−1 is sufficient for the calculation of the coefficients of the two polynomials The coefficients a k with k = 1, . . . , K are obtained as solutions of the linear set of equations Once the a k are known, the coefficients b k are given by the explicit formula Note that the Padé approximant does not only accelerate the convergence of the sum over all c n z −n in Eq. (8) but also yields an analytic continuation for z values, where the sum is not absolutely convergent. The parameters z k = exp(−iω ′ k τ ) and thus the frequencies are obtained by searching for the zeros of the polynomial Q K (z) in Eq. (10). This is the only nonlinear step of the algorithm. The roots of polynomials can be found, in principle, by application of Laguerre's method [38]. However, it turns out that an alternative method, i.e., the diagonalization of the Hessenberg matrix with a k the coefficients of the polynomial Q K (z) in Eq. (10) is a numerically more robust technique for finding the roots of high degree (K 60) polynomials [38]. When the harmonic inversion method is applied to a continuous spectrum of a scattering system, given as a superposition of Lorentzian line shapes, the resulting ω k are the complex poles of the S-matrix.

III. THE MICROSTADIUM
Below a certain cutoff frequency ω cutoff , microdisk-like cavities can be regarded as two-dimensional systems. In this case Maxwell's equations reduce to a two-dimensional scalar mode equation [39] − with effective index of refraction n, frequency ω and the speed of light in vacuum c. The mode equation (15) is valid for both transverse magnetic (TM) and transverse electric (TE) polarization. We focus on TM polarization with the electric field E(x, y, t) = (0, 0, ψ(x, y)e −iωt ) perpendicular to the cavity plane. The wave function ψ and its normal derivative are continuous across the boundary of the cavity. At infinity, outgoing wave conditions are imposed. Note that we ignore in the following a slight frequency dependence of n.
There are only a few reports in the literature on the spectral statistics of microcavities. In a recent work the resonance width distribution in the circular disk has been studied [40]. The integrability of this system allows to easily gather sufficient resonance data. In Ref. [41] the resonance width distribution of a circular cavity with strong surface roughness has been analyzed. However, an additional approximation was necessary in order to get a sufficient amount of resonances.
In our approach, no such approximation is needed. Moreover, it is not restricted to a particular geometry. As example, we consider the dielectric microstadium as illustrated in Fig. 1. The cross-sectional area is bounded by two semicircles with radius R > 0 and two straight lines of length 2L > 0. The corresponding stadium billiard is a paradigm for fully chaotic systems [42,43]. We choose the case L = R which is the most chaotic one [44]. The spatial structure of optical modes in microstadiums has been extensively studied in the context of ray-wave correspondence in open systems [45,46,47,48,49]. Our aim is to study the spectral properties of modes in such a kind of cavity.

IV. HIGH-INDEX MICROSTADIUM
We consider first a GaAs microstadium with a high refractive index n = 3.3, the experimental realization of which has been reported in Refs. [50,51].

A. Spectral analysis of resonances
The modes (15) can be numerically computed with the boundary element method (BEM) using a root search algorithm in the complex frequency plane [25]. Even though the BEM is very efficient, it does not allow for obtaining a sufficient amount of resonances for a statistical analysis. Our improved strategy is therefore based on the BEM and the harmonic inversion technique. It consists of two steps. In a first step, we compute spectral data on the real frequency axis. For convenience, we choose the total scattering cross section σ as function of the normalized frequency Ω = ωR/c = kR using the BEM. For a definition of the scattering cross section we refer to Ref. [25]. We are able to cover the interval Ω ∈ [0, 25] where we found 3772 resonances. The upper bound Ω = 25 corresponds to R = 3.2 µm assuming a typical wavelength of 800 nm in GaAs devices [50,51]. To keep the picture clear, Fig. 2 shows the calculated scattering cross section in a smaller frequency interval. In a second step we analyze the spectral data on the real frequency axis with the help of the harmonic inversion technique. Fig. 3 shows the resulting resonances in the complex plane with Im(Ω) ≥ −0.06. The imaginary part is related to the quality factor via Q = −Re(Ω)/[2Im(Ω)]. The limitations of the harmonic inversion method are as follows: Very weak and broad resonances, i.e., resonances with small amplitudes A k and large widths Γ k in Eq. (4) may be overlooked. Also, in rare cases one true resonance can be fitted by a sum of two Lorentzians with nearly identical resonance positions and widths. We have two ways to verify the harmonic inversion technique. First, we can check exemplarily resonances with the BEM by using root searching in the complex plane. We always find good agreement. Second, from the extracted resonances we can reconstruct the scattering cross section. Fig. 2 reveals good agreement with the original scattering cross section.
The statistics of imaginary parts Im(Ω) = −R/(2cτ ), where τ is the lifetime of the given mode, is depicted in Fig. 4. We see a clear maximum around Im(Ω) ≈ −0.014. Figure 5 shows the number of modes with cutoff C = 0.06 according to Eq. (1) versus normalized frequency Ω = kR in a log-log plot. The data can be extremely well fitted by a straight line with slope α ≈ 1.98. Because of the well pronounced maximum in the statistics of imaginary parts in Fig. 4 we expect only a small change of the exponent α when the cutoff C is varied. Indeed, this expectation is confirmed by the numerics as we find α ∈ [1.96, 2.02] for C ∈ [0.03, 0.1]. As an exponent larger than 2 for two-dimensional open systems is rather unphysical, we believe that the regime [2, 2.02] is solely due to the numerical uncertainty.

B. The chaotic repeller
The chaotic repeller is the set of points in phase space that never leave the system both in forward and in backward time dynamics [18]. The stable (unstable) manifold of a chaotic repeller is the set of points that converges to the repeller in forward (backward) time evolution. The unstable manifold therefore describes the route of escape from the chaotic system. For the case of chaotic microcavities it has been demonstrated that the unstable manifold can play an important role for the far-field emission pattern [48,52,53,54,55].
Here, we discuss the chaotic repeller of the microstadium restricted to the Poincaré surface of section (SOS) [18]. It is a plot of the intersection points of a set of ray trajectories with The reduced phase space is parametrized by the arclength coordinate along the circumference of the cavity's boundary, s, and the tangential momentum, p = sin χ, where χ is the angle of incidence measured from the surface normal; see Fig. 1.
To compute the repeller we portion the phase space above the upper critical line for total internal reflection, p > | sin χ c | = 1/n, into L × L = 1024 × 1024 boxes. For each box we would like to determine whether it belongs to the repeller or not. Obviously, this depends on the chosen point within the box. To soften this arbitrariness the following averaging procedure is used. In each box m = 400 trajectories are started and propagated both forward and backward in time. Given a cutoff length L cutoff = 20.5R, we assign an intensity I = 1 to rays which never enter the leaky region |p| < 1/n, otherwise I = 0. The weight w ∈ [0, 1] associated with each box is then Figure 6 shows the chaotic repeller determined in this way. The magnifications in Fig. 7 demonstrate the fractal nature of the repeller. In order to compute the fractal dimension of the chaotic repeller, we apply the thermodynamic formalism for fractal measures; see, e.g., Ref. [56]. As for the determination of the chaotic repeller we divide the region (x, y) = [s/s max , (p − 1/n)/(1 − 1/n)] ∈ [0, 1] × [0, 1] into L × L = 1024 × 1024 boxes with size ε × ε where ε = 1/L. The weights w kl = w(x k , y l ) from Eq. (16) are normalized such that kl w kl = 1. We write the coarse-grained partition function as  Fig. 7(a). The fractal structure of the chaotic repeller can be clearly seen.
where q is a real-valued parameter. The scaling behaviour of the partition function defines the generalized or Rényi dimensions d(q). The dimension d(q) can be determined numerically for given q = 1 by computing the partition function Γ as in Eq. (17) for increasing box sizes ε M = 2 −M with M = 10, 9, 8, 7, 6, 5. The weight of each resulting box is then the sum of the weights of the corresponding smaller boxes. Finally, the dimension d is extracted by fitting the function (18). The dimension d(0) is equal to the Minkowski dimension, also called capacity or box counting dimension. In the limit q → 1, the Rényi dimension becomes the information dimension Similar to the general case of d(q) with q = 1, d(1) can be determined by fitting d(1) ln ε M to kl w kl ln w kl as function of the box size ε M = 2 −M . The dimension d(2) is called correlation dimension since it is related to correlation functions as, e.g., in Ref. [57].
In the case of a fractal set with uniform measure, i.e., a point belongs to the set (I = 1) or not (I = 0), all dimensions d(q) have the same value. In the more general case of fractal sets with nonuniform measure (I ∈ IR + ) the dimensions d(q) differ. In such a case one speaks about a multifractal. The chaotic repeller is strictly speaking a set with uniform measure. For numerical reasons, however, we have seen that it is more convenient to consider a real-valued I(s, p). That this approximation is justified can be seen by the fact that the numerically computed fractal dimensions d(0), d (1), and d (2) are all approximately 1.68. The predicted exponent according to Eq. (2) and d = 1.68 is then around 1.84. Unexpectedly, this is not close to the Weyl exponent α ∈ [1.96, 2.02] obtained from the counting of resonances. Hence, it seems that the conjectured fractal Weyl law fails for the dielectric microstadium.
In order to see the origin of this failure, let us consider the escape rate γ of the chaotic repeller. Points in the vicinity of the repeller escape exponentially in time with rate γ.
To determine this rate, we remark that the numerically computed repeller in Fig. 6 is an approximation to the real repeller, in that numerically computed phase space points are never exactly located on the repeller. We can exploit this fact to determine γ by evolving the numerically computed repeller in time. Its spatial structure does not change under time evolution, only the total intensity decays exponentially. The corresponding decay rate is nothing else than the escape rate γ. We find γ ≈ 0.056c/R. This translates into Im(Ω) = −γR/(2c) ≈ −0.028 which is significantly larger than the mean escape rate of the quasi-bound modes in Fig. 4. We conclude that the chaotic repeller as defined above also fails to describe the mean lifetimes of the quasi-bound modes.
In the following we show that both problems stem from the fact that the conventional chaotic repeller ignores the partial escape at dielectric interfaces according to Fresnel's laws. Rays entering the leaky region are not completely transmitted to the exterior of the cavity. For TM polarization the reflection coefficient is with the angle of incidence χ and the angle of transmission χ t related by Snell's law n sin χ = sin χ t . We can include Fresnel's laws in the concept of the chaotic repeller by allowing for a real-valued intensity I assigned to each ray. A similar approach has been used to describe far-field emission pattern based on unstable manifolds [48,55]. Initially, we choose I = 1 for a given ray, and whenever the ray enters the leaky region at a phase space point (s, p) we multiply the intensity I by the reflection coefficient R TM (s, p). Finally, the ray's intensity is 0 ≤ I ≤ 1. The chaotic repeller defined in this way is a set with nonuniform measure I(s, p); it is a multifractal.
Using the extended concept of the chaotic repeller we compute the escape rate to be γ ≈ 0.034c/R. This gives Im(Ω) ≈ −0.017 which is in much better agreement with the distribution of imaginary parts of the quasi-bound modes depicted in Fig. 4. Figure 8 shows the chaotic repeller including Fresnel's laws. The fractal dimensions turn out to be d(0) ≈ 1.986, d(1) ≈ 1.913, and d(2) ≈ 1.877. The predicted exponent for the Weyl law is therefore 1.99, 1.96, and 1.94, correspondingly. The values resulting from the box-counting dimension and the information dimension are in the interval of the Weyl exponent α ∈ [1.96, 2.02] calculated from the resonances. For the high-index microstadium we therefore have confirmed the conjectured fractal Weyl law using the extended concept of the chaotic repeller.

V. LOW-INDEX MICROSTADIUM
Now we proceed to the microstadium with low refractive index n = 1.5. Such a cavity has been fabricated recently by using a PMMA polymer matrix [47,58]. Here, we can compute the resonances in the interval Re(Ω) ∈ [0, 75] as can be observed in Fig. 9. For example, Re(Ω) = 75 corresponds to R = 7.2 µm for λ = 600 nm [47,58]. The comparison with Fig. 3 reveals that the resonances of the low-index stadium are much deeper in the complex plane, which reflects the stronger degree of openness due to the larger leaky region | sin χ c | = 1/n. This can also be seen from the statistics of imaginary parts depicted in Fig. 10. Here, we observe a clear maximum around Im(Ω) ≈ −0.15.  The chaotic repeller including Fresnel's laws is shown in Fig. 12. The escape rate for phase space points near the chaotic repeller is determined to be γ ≈ 0.26c/R. This translates into Im(Ω) ≈ −0.13 in good agreement with the distribution of the quasi-bound modes in Fig. 10. Note that Fresnel's laws do not change much here since the reflectivity in the leaky region is in general small due to the low refractive index. The numerically computed dimensions are d(0) ≈ 1.512 and d(1) ≈ d(2) ≈ 1.593. Hence, the predicted exponent is 1.76 and 1.78 which is fully consistent with the fractal Weyl law with exponent α ∈ [1.68, 1.88].

VI. CONCLUSION
We have demonstrated that the harmonic inversion technique allows to compute very efficiently a huge amount of complex resonances in dielectric microcavities. This is of high practical value for the statistical analysis of resonances in the context of quantum chaos and complex scattering. The approach has been illustrated for the chaotic microstadium made of GaAs and polymers. We demonstrate that the fractal Weyl law can be applied to the class of dielectric cavities by extending the conventional concept of the chaotic repeller by including partial escape of light according to Fresnel's laws. We expect that our finding is of high relevance for the understanding of ray-wave correspondence in chaotic microcavities. More-over, it can have implications on multi-mode lasing in such cavities [59].

VII. ACKNOWLEDGEMENT
We would like to thank E. Bogomolny, M. Hentschel, and S. Shinohara for discussions. Financial support from the DFG research group 760 is acknowledged.