Power spectrum and form factor in random diagonal matrices and integrable billiards

Triggered by a controversy surrounding a universal behaviour of the power spectrum in quantum systems exhibiting regular classical dynamics, we focus on a model of random diagonal matrices (RDM), often associated with the Poisson spectral universality class, and examine how the power spectrum and the form factor get affected by two-sided truncations of RDM spectra. Having developed a nonperturbative description of both statistics, we perform their detailed asymptotic analysis to demonstrate explicitly how a traditional assumption (lying at the heart of the controversy) -- that the power spectrum is merely determined by the spectral form factor -- breaks down for truncated spectra. This observation has important consequences as we further argue that bounded quantum systems with integrable classical dynamics are described by heavily truncated rather than complete RDM spectra. High-precision numerical simulations of semicircular and irrational rectangular billiards lend independent support to these conclusions.


Contents 1 Introduction 3
2 Power spectrum and spectral form factor for truncated spectra with stationary spacings 6 2.1 Random spectra with stationary and infinitely-stationary spacings . . .

Introduction
Since the second half of the XX century, a variety of statistical indicators have been devised to study fluctuations in spectra of bounded quantum systems, at both shortand long-range energy scales. For high-lying energy levels, these indicators, or spectral fluctuation measures, were found to exhibit an identical behavior in different quantum systems, thereby revealing the phenomenon of spectral universality [1]. To observe it, an influence of system-specific mean level density has to be eliminated from sequences of measured energy levels {λ 1 , . . . , λ N }. This is achieved by means of the unfolding procedure [2] represented by the map ε ℓ =ˆλ ℓ −∞ dλ̺ N (λ), ℓ = 1, . . . , N, (1.1) where ̺ N (λ) is the mean level density of the original (e.g., measured) spectrum. It is the unfolded energy levels {ε 1 , . . . , ε N } that may obey universal statistical laws as N → ∞.
There is a broad consensus, based on a vast amount of experimental, numerical and theoretical evidence [3,4,5], that a universal statistical behavior of quantum systems correlates with the nature -chaotic or regular -of their underlying classical dynamics. This gives rise to two paradigmatic universality classes in quantum chaology.
The Wigner-Dyson universality class accommodates generic quantum systems which are fully chaotic in the classical limit. According to the Bohigas-Giannoni-Schmit (BGS) conjecture [6], spectral fluctuations of highly excited energy levels in such systems -exhibiting long-range correlations and local repulsion -are governed by global symmetries rather than by system specialties and are accurately described by the random matrix theory [2,7]. Mostly built on numerical indications [8,9,10,6], the BGS conjecture has later been advocated within the field-theoretic [11] and semiclassical [12] approaches. On the contrary, generic quantum systems whose classical dynamics is completely integrable belong to a different -Poisson -universality class as was first conjectured by Berry and Tabor [13]. Spectral fluctuations therein are radically different from those in Wigner-Dyson spectra, with energy levels being completely uncorrelated.
Since the most commonly used spectral fluctuation measures [2,7] -be it the distribution of level spacings [14] or that of the ratio of two consecutive level spacings [15], the number variance [16] or spectral rigidity [17] -probe spectral correlations on either local or global scales, a combination of both short-and long-range statistical indicators is often required to get a reliable information on the degree or type of chaoticity of a quantum system in question. This remark is even more pertinent, if not crucial, for analysis of quantum systems with mixed classical dynamics [18] and for an accurate interpretation of experimental data in which the completeness of energy spectra is rather a rare situation [19,20].
In this context, we would like to draw reader's attention to an alternative measure of spectral fluctuations -the power spectrum of eigenlevel sequences -which is capable of probing the correlations between both nearby and distant eigenlevels. Introduced long ago by Odlyzko in his famous study [21] of the spacing distribution between nontrivial zeros of the Riemann zeta function, and re-invented fifteen years later by Relaño and collaborators [22], it considers a sequence of discrete energy levels in a quantum system as a discrete time series, where indices of ordered eigenlevels play the role of discrete time.
Definition 1.1. Let {ε 1 ≤ . . . ≤ ε N } be a sequence of ordered unfolded eigenlevels, N ∈ N, with the mean level spacing ∆ and let δε ℓ δε m be the covariance matrix of level displacements δε ℓ = ε ℓ − ε ℓ from their mean ε ℓ = ℓ∆. A Fourier transform of the covariance matrix is called the power spectrum of a sequence. Here, the angular brackets stand for an average over an ensemble of eigenlevel sequences.
Remark 1.3. It readily follows from Definition 1.1 that, by tuning the frequency ω in the power spectrum, one may attend to spectral correlations between either adjacent or distant eigenlevels. Indeed, at 'large' frequencies, ω = O(N 0 ) yet below the Nyquist frequency, the distant eigenlevels barely contribute to the power spectrum; characterized by large values of |ℓ − m|, they produce strongly oscillating terms in Eq. (1.2) which effectively cancel each other. As the result, S N (ω) is mainly shaped by correlations between the nearby levels. At low frequencies, ω ≪ 1, these oscillations are by far less pronounced thus making a contribution of distant eigenlevels increasingly important.
Although the power spectrum, as an alternative spectral statistics, was suggested quite some time ago [21,22], its nonperturbative theory [23,24] was formulated only recently. (For an earlier, heuristic approach exploiting a form factor approximation, see Ref. [25].) In particular, in Ref. [23], the power spectrum of energy level fluctuations in fully chaotic quantum structures with broken time-reversal symmetry was determined by employing a random-matrix-theory approach. A universal, parameter-free prediction for the power spectrum expressed in terms of a fifth Painlevé transcendent was the main result of that study. The ideas underlying a randommatrix-theory calculation of Ref. [23] were further developed in Ref. [24], where the power spectrum was studied in a more general setting which is applicable to a much broader class of quantum systems whose eigenspectra -not necessarily of the randommatrix-theory type -possess stationary level spacings.
In this note, we utilize a nonperturbative approach of Ref. [24] to determine the power spectrum for generic quantum systems with integrable classical dynamics ‡. We argue that an adequate description of such systems, modelled by random diagonal matrices [27], must take into account the effect of eigenspectra truncations which are unavoidable in any realistic experimental and numerical setup.
The paper is organized as follows. In Section 2, we introduce the notions of random spectra with stationary [24] and infinitely-stationary level spacings, and meticulously examine how truncations of such spectra affect the power spectrum and the form factor of a quantum system. Our analysis reveals that a much sought-after information about individual fluctuations of eigenlevels in truncated sequences gets hindered by their collective fluctuations. This observation leads us to formulate meaningful definitions ‡ Recent report [26] makes our previous [24] and present study potentially relevant to the problem of universal correlations in spectra of quantum many-body localized systems. Notice that the main analytical result of Ref. [26], see their Eqs. (8) and (9), is a particular case of the spectral form factor derived in Ref. [24] for a model of eigenlevel sequences with identically distributed, uncorrelated level spacings, see Eqs. (1.21) and (1.22) therein. Indeed, setting the characteristic function Ψs(τ ) of level spacings of Ref. [24] to Ψs(τ ) = (1−2iπτ ) −1 describing exponentially fluctuating spacings, one obtains the result of Ref. [26].  Figure 1. Comparison of the regularized power spectrum for quantum irrational rectangular billiards (red solid line, numerical simulation, Section 4.2) with that for complete (blue solid line) and heavily truncated (black dashed line) spectrum of random diagonal matrices (theory, Section 3.3). The power spectrum in rectangular billiards is seen to follow predictions for heavily truncated random diagonal matrices. Oscillatory behavior of two upper curves is a fingerprint of inevitable finite-size effects originating from a fixed length L of 'measured' spectral sequences; in all cases, sequences with L = 2048 were used. Inset: a log-log plot for the power spectra shown up to the Nyquist frequency ω Ny = π; the same color legend is used.
of the power spectrum (Definition 2.7) and of the form factor (Definition 2.11) for truncated eigenlevel sequences. Theorem 2.9 for the power spectrum of truncated spectra is the main result of Section 2. Our treatment of both statistics here is not restricted to a particular spectral model or universality class; yet, it will later be used to address the issue of validity of the form factor approximation in the context of quantum systems with integrable classical dynamics.
In Section 3, a general framework of Section 2 is applied to a model of random diagonal matrices (RDM) often used to mimic spectral fluctuations in quantum systems exhibiting regular classical dynamics, see e.g. Ref. [27].
In Section 3.1, we introduce a procedure of spectral unfolding in the finite-N random diagonal matrices, show that their ordered unfolded eigenlevels are beta-distributed (Proposition 3.1), and prove infinite-stationarity of level spacings on the unfolded scale (Proposition 3.2). These results hold true for RDM of arbitrary size N .
In Section 3.2, we evaluate the power spectrum of complete and truncated eigenlevel sequences drawn from an unfolded RDM spectrum. Both exact (Theorem 3.7, finite-N ) and asymptotic (Corollaries 3.10 and 3.12, N → ∞) expressions for the power spectrum constitute the main results of Section 3.2. Their implications are further discussed in the context of a quantitative description of the universal spectral fluctuations in generic quantum systems with regular classical geodesics. In paticular, we argue that heavily truncated -rather than complete -spectra of random diagonal matrices should be employed to adequately determine the power spectrum in realistic quantum systems. This leads us to re-confirm our previous result [23,24] S int (ω) = 1 2 for the power spectrum of generic 'integrable' quantum systems, away from the origin ω = 0 (here ω 0 > 0 is an arbitrary small real number). Notice the factor 1 / 2 which differs from the factor 1 / 4 erroneously claimed in the earlier study [25] which used a so-called form factor approximation. We stress that the parameter free prediction Eq. (1.4) is a universal limit of a more general expression for the power spectrum [Eq. (3.21) of Corollary 3.10] which accounts for finite-size effects arising due to a finite length L of 'measured' spectral sequences. For readers' convenience, these conclusions are visualized in Fig. 1, where the power spectrum of quantum irrational rectangular billiards is seen to follow the ' 1 / 2 ' (rather than ' 1 / 4 ') law claimed in Eq. (1.4). In Section 3.3, we focus on a nonperturbative analysis of the spectral form factor for complete and truncated eigenlevel sequences sampled from an unfolded, finite-N RDM spectrum. The exact solution for the spectral form factor (Theorem 3.14) is further analysed asymptotically (Theorem 3.19) to pave the way for an in-depth discussion of a status of the form factor approximation. Theorem 3.14 and Theorem 3.19 are the main results of Section 3.3.
Such a discussion is held in Section 3.4, where we compare a large-N behavior of the spectral form factor and the power spectrum in various time and frequency scaling limits, in an attempt to scrutinize a heuristic relation Eq. (3.60) claimed [25] to link the power spectrum to the spectral form factor. While we show that this relation does hold true for complete spectra of infinite-dimensional random diagonal matrices up to ω = O(N 0 ) as soon as ω ≪ 1, it becomes invalid for truncated spectra (of the length L < N ). In the latter case, the validity of the form factor approximation is restricted to very low frequencies ω = O(L −1 ); it breaks down for ω = O(L −1/2 ) and certainly gives no access to the bulk of the power spectrum.
In Section 4, we present a thorough numerical study of the power spectrum and the spectral form factor for quantum (i) semicircular and (ii) irrational rectangular billiards. The numerics lends unequivocal support to our theory which associates generic quantum systems possessing integrable classical dynamics with heavily truncated (rather than complete) spectra of random diagonal matrices.
Appendix A contains proofs of various asymptotic formulae for the Kummer confluent hypergeometric function 1 F 1 of imaginary argument, required for the large-N analysis of the RDM spectral form factor, Section 3.3. To the best of our knowledge, explicit asymptotic results obtained in Appendix A are unavailable in the existing literature on special functions and thus can be of independent interest.
2. Power spectrum and spectral form factor for truncated spectra with stationary spacings

Random spectra with stationary and infinitely-stationary spacings
To set up the scene, we first collect the most useful definitions and statements regarding random spectra with stationary and infinitely-stationary level spacings [24]. Definition 2.1 (Stationarity). Consider an ordered sequence of eigenlevels {0 ≤ ε 1 ≤ . . . ≤ ε N } with N ∈ N. Let {s 1 , . . . , s N } be the sequence of spacings between consecutive eigenlevels such that s ℓ = ε ℓ − ε ℓ−1 with ℓ = 1, . . . , N and ε 0 = 0. The sequence of level spacings is said to be stationary if (i) the average spacing is independent of ℓ = 1, . . . , N and (ii) the covariance matrix of spacings is of the Toeplitz type: for all ℓ, m = 1, . . . , N . Here, I n is a function defined for non-negative integers n.
A necessary and sufficient condition for eigenlevel sequences to possess stationary level spacings, proven in Section 3 of Ref. [24], is formulated in Lemma 2.2.
be an ordered sequence of unfolded eigenlevels such that ε 1 = ∆. An associated sequence of spacings between consecutive eigenlevels is stationary if and only if for ℓ > m and both q = 1 and q = 2.
Remark 2.3. Examples of finite-N eigenlevel sequences with stationary spacings include (i) random spectra with independent, identically distributed (i.i.d.) level spacings ¶; (ii) unfolded spectra in ensembles of random diagonal matrices (see Section 3); (iii) spectra of 'tuned' circular ensembles of random matrices, see Section 4 in Ref. [24]. Interestingly, spectral fluctuations in all aforementioned models appear to satisfy the correlation constraint Eq. (2.3) also for q > 2. This leads us to define the notion of random spectra with infinitly-stationary level spacings.
Definition 2.4 (Infinite stationarity). An ordered sequence of eigenlevels {0 ≤ ε 1 ≤ . . . ≤ ε N } with N ∈ N is said to possess an infinitly-stationary sequence of level spacings if the relation holds for ℓ > m and all q ∈ N.

Power spectrum for truncated eigenlevel sequences with stationary spacings
For complete (untruncated) sequences of eigenlevels with stationary level spacings, the power spectrum can be represented as a differential operator acting on a generating function of eigenlevel variances, see Theorem 2.3 in Ref. [24]: Theorem 2.5. Let N ∈ N and 0 ≤ ω ≤ π. The power spectrum for an eigenlevel sequence {0 ≤ ε 1 ≤ . . . ≤ ε N } with stationary spacings equals

5)
where z = e iω , ∆ is the mean level spacing, and is the variance of ℓ-th ordered eigenlevel.
Remark 2.6. Theorem 2.5, as it stands, cannot directly be applied to describe the power spectrum in a realistic quantum system with unbounded spectrum for two reasons. First, in an experiment one normally measures a truncated spectrum represented by a subsequence + {ε M+1 ≤ . . . ≤ ε M+L } of, say, L eigenlevels obtained from a complete spectrum {0 ≤ ε 1 ≤ . . . ≤ ε N } of the length N by skipping both M low-lying eigenlevels {ε 1 ≤ . . . ≤ ε M } and N − M − L high-lying eigenlevels {ε M+L+1 ≤ . . . ≤ ε N }; here, M + L ≤ N . Second, obviously, the N → ∞ limit should be implemented. While the latter problem is merely technical, the former one raises the question of how the original Definition 1.1 of the power spectrum should be modified to account for truncated eigenlevel sequences in a meaningful way. ¶ As the ℓ-th ordered eigenlevel is a sum of ℓ i.i.d. random variables, ε ℓ = ℓ j=1 s j , both the l.h.s. and r.h.s. in Eq. (2.3) represent the q-th moment of a sum of (ℓ − m) i.i.d. random variables, q ∈ N. (Infinite) stationarity follows immediately. + The eigenlevels are assumed to be unfolded. Unfortunately, there is a serious drawback to this proposal. In order to reveal it, we assume stationarity of eigenlevel spacings [Eq. (2.3) at q = 2] which brings the relation where δε ℓ = ε ℓ −ℓ∆. Substituting Eq. (2.8) into the trial definition Eq. (2.7), we reduce the latter tõ Both terms in Eq. (2.9) can be simplified by observing that M+L ℓ,m=M+1 and M+L ℓ,m=M+1 Combining Eqs. (2.9), (2.10) and (2.11), we manage to rewrite Eq. (2.7) in a more suggestive form * : var[ε ℓ ]z ℓ + δS N ;M,L (ω), (2.12) where Comparing the first term in Eq. (2.12) with the one in Eq. (2.5), one realizes that it does not depend on the number M of dropped low-lying eigenlevels, effectively describing the power spectrum of the first L eigenlevels as if M = 0. The second term in Eq. (2.12), specified explicitly in Eq. (2.13), does depend on M . Nullifying at M = 0, it represents a parasitic signal, which may dramatically mask a contribution of the first term in Eq. (2.12). This is best seen in the model of eigenlevel sequences with identically distributed uncorrelated spacings [24], characterized by the variance var[ε ℓ ] = σ 2 ℓ linearly depending on the number of ℓ-th ordered eigenlevel (here σ 2 is the variance of level spacings). Setting the mean level spacing to unity, ∆ = 1, we observe that the second term in Eq.
Remark 2.8. Definition 2.7 is essentially the power spectrum of a 'collectively tuned' Notice that stationarity of the complete spectrum implies stationarity of a 'collectively tuned' subsequence, in the sense of Definition 2.1.

16)
where z = e iω , ∆ is the mean level spacing, and is the variance of ℓ-th ordered eigenlevel.
Proof. Simple rewriting of Eq. (2.15) yields:  Performing summation over m, one obtains after some algebra: , from now on we shall use the notation S N ;L (ω) instead of S N ;M,L (ω).

Spectral form factor for truncated eigenlevel sequences with stationary spacings
Prompted by the above discussion, we define the form factor of truncated spectrum as follows: The form factor of the truncated spectrum equals Remark 2.12. Definition 2.11 is essentially the form factor of a 'collectively tuned' Proof. Since the infinite stationarity of level spacings, Definition 2.4, implies equality of moments (ε ℓ − ε m ) q and ε q ℓ−m for all q ∈ N, one is allowed to simultaneously shift eigenvalues' indices in the differences (ε ℓ − ε m ) appearing under the sign of averages . . . in Eq. (2.21) by the same integer M : , from now on, we shall use the notation K N ;L (τ ) instead of K N ;M,L (τ ).

Random diagonal matrices: Spectral statistics of truncated eigenlevel sequences
Random diagonal matrices (RDM) are routinely used to model spectral fluctuations in quantum systems exhibiting regular classical dynamics. Yet, much to our astonishment, we failed to find a comprehensive and transparent treatment of this paradigmatic model in the literature. In Section 3.1 we partially bridge this gap: the issues of unfolding, ordering, infinite stationarity and correlations of individual eigenlevels are of our primary interest. These results, combined with a general framework of Section 2, are utilized in Sections 3.2 and 3.3, where the power spectrum and the spectral form factor are calculated nonperturbatively for both complete and truncated spectra of the finite-N RDM-model. The two statistics are further analyzed in various scaling limits to lay the ground for Section 3.4, where we discuss a status of the form factor approximation in the RDM context.

Unfolding, ordering and (infinite) stationarity of level spacings
Unfolding.-Consider an ensemble of N × N random diagonal matrices of the form H = diag(λ 1 , . . . , λ N ), where λ j ∈ R are independent identically distributed random variables, each characterized by the probability density function f λj (λ) = f (λ) for all j = 1, . . . , N . The mean spectral density of H explicitly depends on the fluctuation properties of individual eigenlevels and, thus, is system-specific. Here, the angular brackets denote an appropriate ensemble averaging.
To study the universal properties of random diagonal matrices, one needs to rescale their fluctuating eigenlevels in such a way that the average spacing between two adjacent levels becomes a constant. This can be achieved through the 'unfolding procedure' [2] which introduces a new set of eigenlevels according to the map where j = 1, . . . , N and ε j ∈ (0, N ). The above definition ensures that all ε j 's are independent, uniformly distributed random variables, ε j ∼ U(0, N ). Indeed, the probability density function of ε j equals The latter will now be denoted Here, 1 X is the indicator function; it equals unity if X is true; otherwise it brings zero.
Ordering induces correlations.-Aimed at studying statistical properties of truncated eigenlevel sequences introduced in Section 2, we further concentrate on ordered unfolded eigenlevels §. (i) the probability density g ℓ (ε) of the ℓ-th ordered eigenlevel equals is a univariate beta-variable; (ii) the joint probability density g ℓm (ε, ε ′ ) of the ℓ-th and m-th ordered eigenlevels (ℓ < m) equals is a bivariate beta-variable [28].
Proof. Given the joint probability density of all N unfolded (not ordered!) eigenvalues and and (ℓ > m) .

Power spectrum of complete and truncated sequences
Now that the spadework has all been done, an exact calculation of the power spectrum becomes starightforward.
Complete spectrum.-Due to stationarity of level spacings proven in Proposition 3.2, the power spectrum of a complete eigenvalue sequence can be determined by virtue of Theorem 2.5. Since both the mean can be read off from Eqs. (3.9) and (3.10), the power spectrum of a complete sequence can be calculated from Eq. (2.5) after substituting Remark 3.3. The first, linear in ℓ, term in Eq. (3.13) essentially coincides with the variance of ℓ-th ordered eigenvalue in a simple model of spectral sequences with identically distributed uncorrelated spacings considered in Section 1.3 of Ref. [24]. The second, negative (and quadratic in ℓ), term in Eq. (3.13) is a fingerprint of weak albeit long-range negative correlations between the spacings which appear to be an intrinsic feature ♯ of spectral fluctuations in random diagonal matrices. Hence, as soon as the complete spectra are concerned, the power spectrum in a model of random diagonal matrices should differ from the one in a model of spectral sequences with uncorrelated spacings. Explicit calculations confirm this expectation.
and S (c) Proof. Make use of Theorem 2.5, Proposition 3.2 and substitute Eq. (3.13) into Eq. (2.5) to complete the proof.   Truncated spectrum.-Following the discussion in Section 2.2, we adopt Definition 2.7 of the power spectrum for truncated subsequence {ε M+1 ≤ . . . ≤ ε M+L } of the complete spectrum {0 ≤ ε 1 ≤ . . . ≤ ε N }. The power spectrum of such a truncated subsequence is determined as follows.
be a complete sequence of unfolded eigenlevels, N ∈ N, generated by the map Eq. (3.2) out of ordered eigenvalues {λ 1 ≤ . . . ≤ λ N } of a random diagonal matrix, and let {ε M+1 ≤ . . . ≤ ε M+L } be a spectral subsequence of the length L ∈ N obtained from the complete sequence by and S (c) Proof.
Its amplitude can be significantly higher than the power spectrum S N ;L (ω) itself; see also Eq. (2.14) and a discussion around it. L (ω) deriving from weak long-range correlations between level spacings, bears a weight factor (L + 1)/(N + 2). Since it becomes vanishingly small for subsequences of finite length L as N → ∞, the power spectrum of RDM sequences of finite length reduces to the power spectrum of finite eigenvalue sequences with uncorrelated spacings, see Section 1.3 of Ref. [24].
Remark 3.11. Equation (3.21) coincides with the power spectrum calculated in Ref. [24] for random spectra with uncorrelated level spacings (with unit mean and variance). As L → ∞, two scaling limits of the power spectrum S ∞;L (ω) have been identified in Ref. [24] -in the infrared frequency domain ω = O(L −1 ) and in the bulk ω = O(L 0 ). Corresponding limiting laws can alternatively be recovered from Eqs. (3.23) and (3.24) of Corollary 3.12 upon formal setting α L = 0 therein.
In Fig. 2, our theoretical predictions for regularized ‡ power spectrum are confronted with the results of numerical simulations performed for complete and truncated spectra of random diagonal matrices. Referring the reader to a figure caption for detailed explanations, we plainly notice a perfect agreement between the simulations and the theoretical results.
Motivated by a forthcoming discussion of validity of the form factor approximation in Section 3.4, below we are going to concentrate on two scaling limits of the power spectrum for truncated eigenlevel sequences of random diagonal matrices as their dimensionality N → ∞.  . Regularized power spectrum ω 2 S N;L (ω) as a function of frequency ω for truncated RDM eigenlevel sequences characterized by various truncation parameters α L . (i) Solid curves represent regularized power spectra computed numerically by sampling L = 2048 consecutive eigenlevels out of complete spectra of random diagonal matrices whose N = L/α L independent eigenvalues were uniformly distributed U(0, N ). Averaged over M = 10 6 realizations, the green, purple, yellow, red and blue curves mark the regularized power spectrum for truncation parameters α L = 0.1, 0.2, 0.5, 0.8, and 1, respectively. (ii) Dashed lines on the top of colored ones display theoretical curves from Theorem 3.7. (iii) Dotted line represents a theoretical curve for the regularized power spectrum for eigenlevel sequences of the length L = 2048 with independent and exponentially distributed level spacings (Corollary 3.10 and Remark 3.11). Its comparison with a simulated curve can be found in Fig. 1 of Ref. [24]. Inset: the same graphs for higher frequencies.
L(N ) → ∞ as N → ∞ and the two limits exist, then the regularized power spectrum ω 2 S N ;L(N ) (ω) admits the following scaling limits: and where ω 0 > 0 can be arbitrary small but fixed. Furthermore, any other scaling limit of S N ;L(N ) is trivial, i.e. for ω =Ω/L(N ) δ with the scaling exponent 0 < δ < 1, we have: As far as a complete spectrum of random diagonal matrices is concerned, the limiting laws for the power spectrum are furnished by Corollary 3.12 taken at α L = 1. In this case, the power spectrum in the bulk and near the origin is described by and  .27) have been derived for random diagonal matrices (which are widely believed to mimic spectral statistics in quantum systems with integrable classical dynamics), they do not describe the power spectrum in realistic quantum systems. Indeed, in any conceivable experimental setup, only a finite number L = O(N 0 ) of eigenlevels -out of infinitely many (N → ∞) -can be measured. This situation is described by Corollary 3.10 which brings, as L → ∞, different power spectrum laws in the bulk  .27)] derived for complete RDM spectra. In particular, comparison of the leading terms in Eqs. (3.26) and (3.28) reveals that, at any finite frequency (that is, away from the origin ω = 0), the power spectrum of infinite-dimensional random diagonal matrices is twice smaller than the power spectrum of a spectral subsequence of a large ¶ but finite length. The origin of such a reduction was highlighted in Remark 3.9.
This discussion leads us to conclude that the bulk of the power spectrum in generic quantum systems with completely regular classical geodesics is described by Eq. (3.28) (notice the factor 1 / 2 ) rather than by Eq. (3.26) (notice the factor 1 / 4 ) previously claimed in Refs. [25]. This conclusion is in concert with our previous findings [23,24]. Numerical studies of the power spectrum in rectangular and semicircular billiards, to be presented in Section 4, lend independent support to the universal laws Eqs. (3.28) and (3.29).

Spectral form factor of complete and truncated sequences
Motivated by recent discussions [23,24] of the status of a form factor approximation in the theories of the power spectrum of bounded quantum systems, we now turn to a nonperturbative analysis of the form factor in the model of random diagonal matrices; both complete and truncated eigenspectra will be considered.
Complete spectrum.-To determine the form factor for a complete spectrum of random diagonal matrices, we adopt Definition 2.11: Here, summations run over a set {ε 1 , . . . , ε N } of either ordered or unordered unfolded eigenlevels. Truncated spectrum.-Infinite stationarity of level spacings in ensemble of random diagonal matrices makes it possible to determine the spectral form factor exactly also for truncated sequences. This being said, a calculation of the form factor for truncated spectrum is way more involved since eigenlevel ordering becomes essential.
where 1 F 1 (a; b; z) is the Kummer confluent hypergeometric function.
Proof. By virtue of Lemma 2.13 and Proposition 3.2, the form factor K N ;L (τ ) of the truncated subsequence is defined by Eq. (2.22), where the averaging runs over an ensemble of random diagonal matrices with the mean level spacing ∆ = N/(N + 1), see Eq. (3.11). To facilitate the calculation of K N ;L (τ ) we spot that it can be written in the form 34) where 3.35) and The latter equality is a direct consequence of the infinite stationarity of level spacings in ensemble of random diagonal matrices, see Proposition 3.2.
Having expressed the form factor in term of a single function F N ;L (τ ), we now turn to its evaluation. To proceed, we notice that it can be expressed in terms of the partial mean spectral density of the first L ordered eigenvalues where g ℓ (ε) is the mean density of the ℓ-th ordered eigenlevel, see Eq. (3.5), so that The remaining integral is nothing but a characteristic function for the random variable Y ∼ Beta 1 (L, N − L). This observation brings Comparison with simulated curve can be found in Fig. 2 of Ref. [24]. Inset: the same graphs for higher τ .  In Section 3.2 (see Remark 3.9 and Corollary 3.10) we have shown that the power spectrum of spectral subsequences of finite length L drawn from infinite-dimensional random diagonal matrices coincides with the one for finite eigenvalue sequences with uncorrelated spacings addressed in Section 1.3 of Ref. [24]. Not surprisingly, such a correspondence between two spectral models holds for their spectral form factors, too, see Remark 3.17 below. Figure 3 demonstrates a perfect agreement between the analytical prediction of Theorem 3.14 and the form factor computed out of simulated RDM spectra.
(1 + 4π 2 τ 2 ) L/2 .  To address the issue of validity of the form factor approximation raised in our previous papers [23,24], we need to determine an asymptotic behavior of the spectral form factor K N ;L (τ ) for truncated unfolded spectra of random diagonal matrices as N → ∞. Various scaling limits of the form factor are established in Theorem 3.19 below. These large-N results, combined with a similar analysis of the power spectrum presented in Section 3.2, will set the scene for discussion in Section 3.4. The proof of Theorem 3.19 relies heavily on various asymptotic properties of the Kummer confluent hypergeometric function. Their detailed study is presented in Appendix A. where L(N ) + M (N ) ≤ N . Then, the form factor K N ;L (τ ) admits the following scaling limits: and (3.50) Furthermore, any other scaling limit of K N ;L is trivial, i.e. for τ = T/L(N ) δ with fixed scaling exponent δ, we have (3.51) Proof. (i) To prove Eq. (3.49), we start with the exact Eq. (3.33) of Theorem 3.14 and make use of Proposition A.3 to deduce the asymptotic expansion of Eq. (A.36) proven in Proposition A.6 to observe that as N → ∞. Notice that r.h.s. in Eqs. (iii-b) For 1 / 2 < δ < 1, Proposition A.6 (see also Eq. (3.53)) implies that   α L (τ ), describing the regimes (I), (II) and (III), correspondingly, are plotted vs variables T = Lτ , T = L 1/2 τ and τ , each running over the entire real half-line compactified using the transformation (0, ∞) = tan((0, π/2)). Solid blue, red, yellow, purple and green curves correspond to the form factor for truncated RDM sequences with truncation parameters α L = 1, 0.75, 0.5, 0.25, 0. The dashed grey lines display the limiting curve for the regularized power spectrum lim N→∞ ω 2 S N;L (ω) with 0 ≤ ω = 2πτ ≤ π (that is, 0 ≤ τ ≤ 1 / 2 ). In the scaling regimes (I), (II) and (III), the curve is described by Eq.

Form factor approximation: Power spectrum vs form factor
Equipped with the nonperturbative results for the power spectrum and the form factor for the exactly solvable model of random diagonal matrices, we are in position to examine a status of the so-called form factor approximation, introduced heuristically in Ref. [25]. Assuming that the power spectrum of a quantum system is merely determined by its two-point spectral correlations, the authors of Ref. [25] claimed to relate the power spectrum S(ω) of a quantum system to its spectral form factor K(τ ) in the entire frequency domain 0 ≤ ω ≤ π. Referring interested readers to Eqs. (3), (8) and (10) of the original paper Ref. [25], we only quote a small-ω reduction of their result: see also Ref. [30].
correspondingly. This immediately implies that the heuristic equality Eq.
sin Ω Ω (3.63) and  (3.65) and the form factor   As the two limits Eqs. (3.67) and (3.68) differ from each other for all 0 ≤ α L < 1, we conclude that, in the domain of finite frequencies ω ≪ 1 and times τ ≪ 1, the form factor approximation is broken down, too, for any truncated eigenlevel sequence. The heavier the truncation the higher discrepancy is; the discrepancy factor 2 − α L reaches its maximum 2 for heavily truncated sequences (α L = 0) describing spectral fluctuations in realistic quantum systems with integrable classical dynamics, see Discussion in Section 3.2. The above consideration is visualized in Fig. 4, where we compare the N → ∞ limiting curves for the regularized power spectrum and for the form factor, in all three scaling regimes. While both spectral indicators are on the top of each other in the first scaling regime (characterized by extremely low frequencies ω = O(L −1 ) and short times τ = O(L −1 )), they start to deviate from each other already at ω = O(L −1/2 ) and τ = O(L −1/2 ), except for α L = 1 (the case of complete spectrum). The deviations keep growing with further increase of ω and τ . Figure 5 demonstrates that a similar behavior is also observed for RDM spectra of a large but finite size.
To conclude, our analysis shows that the heuristic relation Eq. (3.60), considered in the context of infinite-dimensional random diagonal matrices, holds true for their complete spectra up to ω = O(L 0 ) as long as ω ≪ 1. However, for truncated spectra , This includes the case α L = 0 which corresponds to a model of random spectra with uncorrelated spacings discussed in Ref. [24]. the form factor approximation breaks down much earlier, at the very beginning of the scaling regime (II), see Fig. 4. Hence, the form factor approximation is restricted to very low frequencies ω = o(L −1/2 ) which give no access to the bulk of the power spectrum.

Numerical study of power spectrum in quantum systems with integrable classical dynamics
In numerical studies of spectral observables of individual quantum systems, an energy averaging, rather than an ensemble averaging, is usually involved. Specifically, an observable of interest is computed along a spectral segment of the length L and then averaged over various locations of such segments in the spectrum [31]. In the context of our study, such an energy averaging is straightforward to define for both the power spectrum and the form factor Here, {E 1 ≤ E 2 ≤ . . .} is a set of deterministic, ordered and unfolded eigenlevels of an individual quantum system; angular brackets (. . .) M0 stand for energy averaging performed by sampling spectral segments of the length L at various locations M 0 along unbounded spectrum. Following a detailed discussion of collective fluctuations in Section 2, both Eqs. (4.1) and (4.2) deal with tuned spectral segments. To capture the universal aspects of spectral statistics expected to emerge for high-lying eigenlevels, the number of skipped eigenlevels M 0 should be large enough (see below).

Semicircular billiards
To study the power spectrum in a circular billiard of the radius R, we start with the Schrödinger equation Owing to the identity κ n,−m = κ n,m , the energy spectrum Eq. (4.5) is double degenerate for all m = 0. To remove such a degeneracy, we follow the idea of Ref. [32] and consider a semicircular billiard comprised by the domain D = {0 ≤ r ≤ R; 0 ≤ θ ≤ π}. For the Dirichlet boundary conditions, the eigenstates of D χ n,m (r, θ) = sin(mθ) J m κ n,m r R , n, m = 1, 2, . . . , (4.6) are antisymmetric with respect to the reflection across the line θ = 0 and the eigenspectrum is given by Eq. (4.5) with n, m = 1, 2, . . .. The large-E asymptotics of spectral counting functionN (E), describing the number of eigenlevels below energy E, is given by the Weyl law [33] where A = πR 2 /2 is the semicircle area, L = (2 + π)R is its circumference, and C is a curvature-and-corner's contribution whose precise value is known but irrelevant for our purposes. In what follows we set R = 2 √ 2 to ensure that the unit mean spacing ∆ = (dN/dE) −1 = 1 as E → ∞.
In concert with a discussion in the beginning of Section 4, the following computational protocol will be used to compute the form factor and power spectrum in the semicircular billiard.

Protocol 1.
To generate an ensemble E = {E 1 , . . . , E Q } of Q eigenlevel sequences (Q ≫ 1 is required for high-accuracy computation of the form factor and the power spectrum), we produce a sufficiently long, ordered, eigenlevel sequence {E 1 ≤ E 2 ≤ . . .} for a semicircular billiard of the radius R = 2 √ 2. To attain the universal spectral regime, we choose a high-lying threshold energy E * > 10 9 (measured in the units of mean level spacing ∆ = 1) and count the number M (E * ) of eigenvalues below the threshold E * ; all of them will be skipped. The first truncated spectral subsequence E 1 of the length L belonging to the ensemble E is set to be E 1 = {E M(E * )+1 − E M(E * ) , . . . , E M(E * )+L − E M(E * ) }, notice the tuning! To produce the remaining (Q − 1) spectral subsequences of the ensemble E, we generate a set of (Q − 1) independent discrete random variables {n 1 , n 2 , . . . , n Q−1 }, with n k denoting a random gap between the last eigenvalue of the k-th and the first eigenvalue of the (k + 1)-th truncated subsequence of the same length L. The j-th member of the ensemble E is comprised by truncated subsequence of the form The averaging procedure is then performed over the ensemble E = {E 1 , . . . , E Q } of Q ≫ 1 spectral sequences, each of length L.
This protocol was applied for numerical studies of the power spectrum and the spectral form factor in a semicircular billiard of the radius R = 2 √ 2 at energies larger than the threshold energy E * = 10 11 . The spectral data were gathered by generating Q = 10 5 spectral sequences of length L = 1024 each, separated by independent, uniformly distributed random gaps {n 1 , n 2 , . . . , n Q−1 }, with ‡ n j ∼ U (50, 150). The results are presented in Fig. 6.
(i) Spectral form factor and power spectrum computed numerically for semicircular billiards show excellent agreement with the theoretical predictions derived for heavily truncated spectra of infinite-dimensional random diagonal matrices (Theorem 3.16 and Corollary 3.10). In accordance with the discussion in Sections 3.2 and 3.3, this sector (α L = 0) of the theory of RDM spectra is equivalent to a spectral model of finite eigenlevel sequences with independent, exponentially distributed level spacings, see also Ref. [24].
(ii) For comparison, we have also plotted theoretical curves for the form factor and the power spectrum calculated for complete spectra of random diagonal matrices (see Propositions 3.4 and 3.13). It appears that the two have nothing in common with the corresponding statistics in semicircular billiards thus invalidating a direct use of random diagonal matrices to description of spectral fluctuations in realistic quantum systems with integrable classical dynamics.
(iii) As a side remark, we note that while the power spectrum in truncated and complete RDM spectra are nowhere the same, their form factors approach unity in both models as τ grows; see Fig. 4 for further details.

Irrational rectangular billiards
In the case of a rectangular billiard with sides a and b, we start with the Schrödinger equation The spectral counting functionN (E), describing the number of eigenlevels below energy E, is asymptotically given by the Weyl law Eq. (4.7), where A = ab is the billiard area, L = 2(a + b) is its circumference, and C is a corner's contribution.
In what follows we set A = 4π to ensure the unit mean spacing ∆ = (dN/dE) −1 = 1 as E → ∞. To focus on generic, universal aspects of spectral statistics, we choose the squared aspect ratio to be an irrational number hereby suppressing contributions of non-universal accidental spectral degeneracies [13,34]. This leaves us with the energy levels parametrized as follows: To compute the form factor and the power spectrum with high accuracy, one has to create an ensemble of eigenlevel sequences and appropriately perform ensemble averaging. Contrary to the numerical study of a semicircular billiard, where only one protocol was devised, this can be achieved in several ways.
Below we discuss three different computational protocols particularly suitable for numerical studies of rectangular billiards. In Protocol 1, statistical ensemble arises as a result of samplings of a sufficiently long deterministic eigenlevel sequence 'measured' for an individual billiard, see the discussion of 'energy averaging' in the beginning of Section 4. The Protocols 2 and 3 are inapplicable to studies of an individual billiard; instead, the two are tailor-made for statistical analysis of an ensemble of billiards with different aspects ratios but asymptotically equivalent mean spacings. Choosing the ratios at random introduces a true randomness similar to the one encountered in our analysis of ensembles of random diagonal matrices. Computation of the form factor and the power spectrum according to Protocol 2 is in one-to-one correspondence with the tuning procedure discussed and implemented in Section 2. The Protocol 3, describing a somewhat heuristic way of dealing with collective fluctuations, appears to be more effective numerically-wise. ≤ . . .} produced for a fixed value of the squared aspect ratio that we set to the golden ratio α 0 = (1 + √ 5)/2. Similarly to the Protocol 1 for the semicircular billiard, we choose a high-lying threshold energy E * > 10 9 to attain the universal spectral regime and count the number M (E * ; α 0 ) of eigenvalues below the threshold E * ; all of them will be skipped. The first truncated spectral subsequence E ′ 1 of the length L belonging to the ensemble E ′ is chosen to be notice the tuning! Remaining (Q − 1) spectral subsequences of the ensemble E ′ will be produced by generating a set of (Q − 1) independent random gaps {n 1 , n 2 , . . . , n Q−1 }, where n k denotes a random gap between the last eigenvalue of the k-th and the first eigenvalue of the (k + 1)-th truncated subsequence of the same length L. The j-th member of the ensemble E ′ is then comprised by truncated subsequence of the form The averaging procedure is then performed over the ensemble M+L } taking particular care that there are no spectral degeneracies therein, within the machine precision. Finally, we tune this subsequence, as discussed in Section 2.2 for the power spectrum (Definition 2.7 and Remark 2.8) and in Section 2.3 for the form factor (Definition 2.11, Lemma 2.13 and Remark 2.10) to obtain the j-th member of the spectral ensemble, The averaging procedure is then performed over the ensemble E ′′ = {E ′′ 1 , . . . , E ′′ Q } of Q ≫ 1 eigenlevel sequences, each of length L.

Protocol 3.
There exists an alternative way to get rid of collective fluctuations in the truncated sequence of eigenlevels, as described below. To build an appropriate ensemble E ′′′ = {E ′′′ 1 , . . . , E ′′′ Q } of Q eigenlevel sequences, we again generate Q aspect ratios {α 1 , . . . , α Q } at random, see Protocol 2. For each α j , we compute a sufficiently long, ordered, eigenlevel sequence {E ≤ . . .} by virtue of Eq. (4.13). However, in distinction to Protocol 2, a truncated sequence of L ordered eigenlevels will be constructed in a slightly different manner. Specifically, we choose a high-lying threshold energy E * > 10 9 (measured in units of the mean level spacing ∆ = 1) and count the number of eigenvalues below the threshold E * ; say, there will be M (E * ; α j ) of such eigenvalues; all of them will be discarded. The number M (E * ; α j ) is a random variable as is the aspect ratio α j . The associated truncated sequence of L ordered eigenlevels comprising the j-th member of the spectral ensemble reads: We verify, within the machine precision, that E ′′′ j is free of spectral degeneracies. It should be noticed that spectral subsequence E ′′′ j does not exhibit collective fluctuations since its lowest eigenvalue -the first above the threshold energy -may only fluctuate on the scale of the mean level spacing. The averaging procedure is then performed over the ensemble E ′′′ = {E ′′′ 1 , . . . , E ′′′ Q } of Q ≫ 1 eigenlevel sequences, each of length L.
The Protocol 3 was applied to study the power spectrum and the spectral form factor in an ensemble of Q = 10 6 rectangular billiards with randomly chosen aspect ratios α j ∈ (1, 2) at energies larger than the threshold energy E * = 10 13 . For each billiard realization, a set of L = 2048 ordered eigenlevels was recorded for further statistical processing. Remarkably, the power spectrum + and the form factor computed numerically for an ensemble of rectangular billiards followed closely the theoretical curves Eqs. (3.21) and (3.47), respectively. In fact, a correspondence between simulations and theoretical predictions for irrational rectangular billiards is even better than for semicircular billiards. This can be attributed to a higher threshold energy E * in the former case. Hence, the discussion and conclusions presented in the end of Section 4.1 equally apply to rectangular billiards. Numerical implementation of Protocols 1 and 2 produced very similar curves for the two statistics, thus validating all three computational protocols.

Appendices
A. Asymptotics of the Kummer confluent hypergeometric function For any x ∈ C and Re b > Re a > 0, the Kummer confluent hypergeometric function 1 F 1 (a, b; ix) admits the integral representation Lemma A.1. For all b ≥ a ≥ 0 and x ∈ R, the following bound holds: Proof. For any x ∈ R and 0 < a < b we have Since 1 F 1 (0; b; ix) = 1 and 1 F 1 (a; a; ix) = e ix , the bound also holds for a = 0 and a = b.
Proposition A.2. Let a > 0, x ∈ R and j ∈ R. Then, as λ → ∞, the expansion holds uniformly for any bounded a, x and j.
Proof. As soon as 1 F 1 (a; b; i0) = 1, we shall only consider x = 0. Furthermore, it is sufficient to treat the case x > 0 only since the domains x > 0 and x < 0 are related to each other by virtue of the symmetry relation I(a; b; −x) = I * (a; b; x).
Owing to the integral representation Eq. (A.1), we shall focus on the asymptotic analysis of I(a, λ + j, λx), Eq. (A.2), as λ → ∞. This will be done by applying the Laplace method [35] after deforming the integration contour t ∈ (0, 1) in the integral into the integration path γ lying in the complex plane and coinciding, at least partially, with a path γ 0 of steepest descent determined by equation The latter brings an explicit parameterization of the path γ 0 in the form z = 1 − re −iφ , where In particular, if x = 0, the real line will be a path of steepest descent. For x > 0, a path of steepest descent is described by the curve where φ 0 = min(x, π/2). Depending on the values of x, Eq. (A.10) brings out three different types of the deformed integration contour γ.
(i) If 0 < x < π/2, the path of steepest descent γ 0 is a finite-length curve starting at the origin z = 0 and finishing at the endpoint z = 1, see the left panel in Fig. 7. Hence, the integration path γ in Eq. (A.7) is identical to γ 0 .
(ii) If x = π/2, the path of steepest descent γ 0 is a finite-length curve starting at the origin z = 0 and finishing at the point z = 1 + 2i/π thus not reaching the endpoint z = 1. To close the integration contour, we need to add the straight (steepest descent) path γ 1 connecting the points z = 1 + 2i/π and z = 1, see the middle panel in Fig. 7. Hence, the integration path γ in Eq. (A.7) is γ = γ 0 ∪ γ 1 .
(iii) Finally, if x > π/2, the path of steepest descent γ 0 is a curve of infinite length starting at the origin z = 0 and approaching complex infinity at z = π/(2x) + i∞. To build the integration contour γ, we need to return from complex infinity along the straight (steepest descent) path γ 2 connecting the points z = 1 + i∞ and z = 1, and connecting the curves γ 0 and γ 2 with an arc γ ∩ , see the right panel in Fig. 7. Hence, the integration path γ in Eq.
Notice that for x ≥ π/2 the integration over additional paths γ 1 and γ 2 will bring exponentially small corrections to the dominating contribution coming from the contour γ 0 as λ → ∞. This holds true since for any z belonging to γ 1 , γ 2 we observe the inequality Re h 0 (x; z) > h 0 (x; 0) = 0. The integration over γ ∩ brings no contribution even for finite λ.
Hence, we are left to consider a contribution I γ0 (a, λ+ j, λx) to the integral Eq. (A.7) coming from a path of steepest descent γ 0 . Making use of Eq. (A.10), we obtain: and Since the function h(φ), appearing in the exponent of the integral Eq. (A.11), reaches its minimum at the origin φ = 0 and h(0) = 0, the main contribution to the integral I γ0 (a, λ + j, λx) comes from a neighborhood of φ = 0. To determine the dominating contribution to I γ0 (a, λ + j, λx), we make use of Theorem 1 from Chapter II in Ref. [35] to write down a λ → ∞ asymptotic expansion where the coefficients c k can be related to the coefficients a k and b k arising in the expansion of the functions h and f at φ = 0, The leading order term in Eq. (A.14) originates from k = 0 for which c 0 = b 0 /a a 0 . Given that 16) we eventually derive the asymptotic expansion of I γ0 (a, λ + j, λx) in the form as λ → ∞. This essentially finishes the proof after noticing that the pre-factor in Eq. (A.1) exhibits an asymptotic behavior as λ → ∞ while a and j are kept fixed. Finally, since for any x ∈ R the saddle point φ * of h 0 , occuring at φ * = 1 + i/x, stays away from the origin φ = 0 at which the global minimum of h 0 is achieved along the path γ, we conclude that the asymptotic expansion holds uniformly in x. Furthermore, since the integrand is regular in the vicinity of the origin φ = 0 for any j and a > 0, the uniformity in j and a is also secured. Proposition A.3. Let 0 < b < 1, x ∈ R and j ∈ R. Then, as λ → ∞, the asymptotic expansion holds uniformly for any bounded x and j.
Proof. With Eq. (A.1) in mind, we shall focus on the asymptotic analysis of I(bλ + j, λ+j, λx), Eq. (A.2), as λ → ∞. To proceed, we start with the integral representation Real-valuedness of the exponent h b (t) allows us to apply the Laplace method, see Theorem 1 from Chapter II in Ref. [35]. A saddle point t * of the large-λ integral Eq. (A.20) is seen to be located at t * = b at which h b (t) reaches its minimum for t ∈ [0, 1]. In view of Theorem 1 of Ref. [35], expansions written in a vicinity of t * = b are of particular (technical) interest. There, and To meet requirements of the aforementioned Theorem 1, the point t * at which the exponent h b (t) reaches its minimum should coincide with a lower limit of the integral Eq. (A.20). To achieve this, one has to split the integral at the saddle point t * into two integrals over domains (0, t * ) and (t * , 1). It appears that their asymptotic expansions, as λ → ∞, partially cancel each other furnishing a series I(bλ + j, λ + j, x) = 2e −λh b (b) Keeping the k = 0 and k = 1 contributions in Eq. (A.27), we derive Combined with the large-λ expansion of the pre-factor in Eq. (A.1), . be a function of z ∈ C with 0 < b < 1 and x ≥ 0 possessing two saddle points and let γ b be a path in the complex plane such that Imh b (x; z) = Imh b (x; z − ) along γ b . Then, for x small enough, such a path connects the points z = 0 and z = 1 passing through the saddle point z − but not through z + .
Proof. Due to a forthcoming application of this lemma to an asymptotic evaluation of the integral Eq. (A.43), we shall call γ b a path of steepest descent throughout the proof. Let us parameterize the path γ b passing through the saddle point z − as where χ b (φ) ∈ R is a solution to the equation Im h b (x; z b (φ)) = Im h b (x; z − ) whose explicit form reads Since Im h b (x; z − ) is a continuous function of x in a vicinity of x = 0, then so is its solution χ b (x; φ). As x → 0, the path γ b approaches the real line, becoming a segment (0, 1) at x = 0. This is consistent with the behavior of z − as a function of x as x → 0 such that z − = b at x = 0.
To prove that z = 0 belongs to the path of steepest descent for sufficiently small x, we must show that χ(φ) → 0 as φ → 0. Let us assume that this is not the case, that is χ(φ) → χ * = 0. Then the φ → 0 limit of Eq. (A.34) yields 35) which should hold in a vicinity of x = 0. In particular, at x = 0, it brings χ * = 0 which contradicts our observation that γ b becomes a segment (0, 1) along the real line at x = 0. (Essentially, this is secured by continuity of χ b (x; φ).) Hence, z = 0 does belong to γ b . A similar reasoning can be applied to show that the point z = 1 belongs to γ b as well.
Finally, spotting that |z + | → ∞ as x → 0, we conclude that z + stays away from γ b . The path of steepest descent γ b is illustrated in the left panel of Fig. 8.
Remark A.5. The topology of the steepest descent path depends on the value of x. For x small enough, this path was discussed in Lemma A.4. As x grows, one approaches a critical situation in which a path of steepest descent goes through both saddle points z − and z + as illustrated in the middle panel of Fig. 8. A critical value of x is determined by the equation Im h b (x; z − ) = Im h b (x; z + ). Further increase of x drives a steepest descent path connecting the points z = 0 and z − to infinity (see red curve in the right panel of Fig. 8). To reach the endpoint z = 1, one has to choose another steepest descent path originating at infinity and approaching z = 1 through z + as depicted by a blue curve in the right panel of Fig. 8.
Proposition A.6. Let 0 < b < 1, x ∈ R, j ∈ R and x 0 be sufficiently small. Then, as λ → ∞, the asymptotic expansion  Remark A.7. There are strong indications that the asymptotic result Eq. (A.36) holds true irrespective of the smallness of x provided that b = 1/2. According to Remark A.5, an asymptotic evaluation of the integral Eq. (A.43) for larger values of x necessitates a deformation of the integration contour which has to pass through both saddle points z − and z + , see the middle and right panels in Fig. 8. Interestingly, even though both saddle points belong to the integration path, numerics provides an evidence that the integral Eq. (A.43) is still dominated by the contribution of z − . This is not the case for b = 1/2 when both saddle points become equally important.