Nonperturbative theory of power spectrum in complex systems

The power spectrum analysis of spectral fluctuations in complex wave and quantum systems has emerged as a useful tool for studying their internal dynamics. In this paper, we formulate a nonperturbative theory of the power spectrum for complex systems whose eigenspectra – not necessarily of the randommatrix-theory (RMT) type – posses stationary level spacings. Motivated by potential applications in quantum chaology, we apply our formalism to calculate the power spectrum in a tuned circular ensemble of random N × N unitary matrices. In the limit of infinite-dimensional matrices, the exact solution produces a universal, parameter-free formula for the power spectrum, expressed in terms of a fifth Painlevé transcendent. The prediction is expected to hold universally, at not too low frequencies, for a variety of quantum systems with completely chaotic classical dynamics and broken time-reversal symmetry. On the mathematical side, our study brings forward a conjecture for a double integral identity involving a fifth Painlevé transcendent. (Dated: October 16, 2019) ‡ Corresponding author.


Introduction
The power spectrum analysis of stochastic spectra [1] had recently emerged as a powerful tool for studying both system-specific and universal properties of complex wave and quantum systems. In the context of quantum systems, it reveals whether the corresponding classical dynamics is regular or chaotic, or a mixture of both, and encodes a 'degree of chaoticity'. In combination with other long-and shortrange spectral fluctuation measures, it provides an effective way to identify system symmetries, determine a degree of incompleteness of experimentally measured spectra, and get the clues about systems internal dynamics. Yet, the theoretical foundations of the power spectrum analysis of stochastic spectra have not been settled. In this paper, a nonperturbative theory of the power spectrum will be presented.
To set the stage, we review traditional spectral fluctuation measures (Section 1.1), define the power spectrum (Definition 1.1) and briefly discuss its early theoretical and numerical studies as well as the recently reported experimental results (Section 1.2). We then argue (Section 1.3), that a form-factor approximation routinely used for the power spectrum analysis in quantum chaotic systems is not flawless and needs to be revisited.

Short-and long-range measures of spectral fluctuations
Spectral fluctuations of quantum systems reflect the nature -regular or chaoticof their underlying classical dynamics [2,3,4]. In case of fully chaotic classical dynamics, hyperbolicity (exponential sensitivity to initial conditions) and ergodicity (typical classical trajectories fill out available phase space uniformly) make quantum properties of chaotic systems universal [3]. At sufficiently long times t > T * , the single particle dynamics is governed by global symmetries of the system and is accurately described by the random matrix theory (RMT) [5,6]. The emergence of universal statistical laws, anticipated by Bohigas, Giannoni and Schmit [3], has been advocated within a field-theoretic [7,8] and a semiclassical approach [9] which links correlations in quantum spectra to correlations between periodic orbits in the associated classical geodesics. The time scale T * of compromised spectral universality is set by the period T 1 of the shortest closed orbit and the Heisenberg time T H , such that T 1 ≪ T * ≪ T H .
Several statistical measures of level fluctuations have been devised in quantum chaology. Long-range correlations of eigenlevels on the unfolded energy scale [5] can be measured by the variance Σ 2 (L) = var[N(L)] of number of levels N(L) in the interval of length L. The Σ 2 (L) statistics probes the two-level correlations only and exhibit [10] a universal RMT behavior provided the interval L is not too long, 1 ≪ L ≪ T H /T 1 . The logarithmic behavior of the number variance, indicates presence of the long-range repulsion between eigenlevels. Here, β = 1, 2 and 4 denote the Dyson symmetry index [5,6]. For more distant levels, L ≫ T H /T 1 , systemspecific features show up in Σ 2 chaos (L) in the form of quasi-random oscillations with wavelengths being inversely proportional to periods of short closed orbits.
Individual features of quantum chaotic systems become less pronounced in spectral measures that probe the short-range fluctuations as these are largely determined by the long periodic orbits [9]. The distribution of level spacing between (unfolded) consecutive eigenlevels, P (s) = δ(s − E j + E j+1 ) , is the most commonly used shortrange statistics. Here, the angular brackets denote averaging over the position j of the reference eigenlevel or, more generally, averaging over such a narrow energy window that keeps the classical dynamics essentially intact. At small spacings, s ≪ 1, the distribution of level spacings is mostly contributed by the two-point correlations, showing the phenomenon of symmetry-driven level repulsion, P (s) ∝ s β . (In a simple-minded fashion, this result can be read out from the Wigner surmise [5]). As s grows, the spacing distribution becomes increasingly influenced by spectral correlation functions of all orders. In the universal regime (s T H /T * ), these are best accounted for by the RMT machinery which produces parameter-free (but β-dependent) representations of level spacing distributions in terms of Fredholm determinants/Pfaffians and Painlevé transcendents. For quantum chaotic systems with broken time-reversal symmetry (β = 2) -that will be the focus of our study -the level spacing distribution is given by the famous Gaudin-Mehta formula [6,11] P chaos (s) = d 2 ds 2 exp Here, σ 0 (t) is the fifth Painlevé transcendent defined as the solution to the nonlinear equation (ν = 0) subject to the boundary condition σ 0 (t) = −t/2π − (t/2π) 2 + o(t 2 ) as t → 0. The universal RMT laws [Eqs. (1.1) and (1.2)] apply to quantum systems with completely chaotic classical dynamics. Quantum systems whose classical geodesics is completely integrable belong to a different, Berry-Tabor universality class [12], partially shared by the Poisson point process. In particular, level spacings in a generic integrable quantum system exhibit statistics of waiting times between consecutive events in a Poisson process. This leads to the radically different fluctuation laws: the number variance Σ 2 int (L) = L is no longer logarithmic while the level spacing distribution P int (s) = e −s becomes exponential [2], with no signatures of level repulsion whatsoever. Such a selectivity of short-and long-range spectral statistical measures has long been used to uncover underlying classical dynamics of quantum systems. (For a large class of quantum systems with mixed regular-chaotic classical dynamics, the reader is referred to Refs. [13,14,15].)

Power spectrum: Definition and early results
To obtain a more accurate characterization of the quantum chaos, it is advantageous to use spectral statistics which probe the correlations between both nearby and distant eigenlevels. Such a statistical indicator -the power spectrum -has been suggested in Ref. [16]. it is sufficient to consider it in the interval 0 ≤ ω ≤ ω Ny , where ω Ny = π is the Nyquist frequency. In the spirit of the discrete Fourier analysis, one may restrict dimensionless frequencies ω in Eq. (1.4) to a finite set with k = {1, 2, . . . , N/2}, where N is assumed to be an even integer. We shall see that resulting analytical expressions for S N (ω k ) are slightly simpler than those for S N (ω).
Remark 1.2. We notice in passing that similar statistics has previously been used by Odlyzko [17] who analyzed power spectrum of the spacings between zeros of the Riemann zeta function.
Considering Eq. (1.4) through the prism of a semiclassical approach, one readily realizes that, at low frequencies ω ≪ T * /T H , the power spectrum is largely affected by system-specific correlations between very distant eigenlevels (accounted for by short periodic orbits). For higher frequencies, ω T * /T H , the contribution of longer periodic orbits becomes increasingly important and the power spectrum enters the universal regime. Eventually, in the frequency domain T * /T H ≪ ω ≤ ω Ny , long periodic orbits win over and the power spectrum gets shaped by correlations between the nearby levels. Hence, tuning the frequency ω in S M (ω) one may attend to spectral correlation between either adjacent or distant eigenlevels.
Numerical simulations [16] have revealed that the average power spectrum S N (ω k ) discriminates sharply between quantum systems with chaotic and integrable classical dynamics. While this was not completely unexpected, another finding of Ref. [16] came as quite a surprise: numerical data for S N (ω k ) could be fitted by simple powerlaw curves, S N (ω k ) ∼ 1/ω k and S N (ω k ) ∼ 1/ω 2 k , for quantum systems with chaotic and integrable classical dynamics, respectively. In quantum systems with mixed classical dynamics, numerical evidence was presented [18] for the power-law of the form S N (ω k ) ∼ 1/ω α k with the exponent 1 < α < 2 measuring a 'degree of chaoticity'. The power spectrum of interface fluctuations in various growth models belonging to the (1 + 1)-dimensional Kardar-Parisi-Zhang universality class, studied in Ref. [19] both numerically and experimentally, were found to follow the power law with α = 5/3. The power spectrum was also measured in Sinai [20] and perturbed rectangular [21] microwave billiards, microwave networks [22,23] and three-dimensional microwave cavities [24]. For the power spectrum analysis of Fano-Feshbach resonances in an ultracold gas of Erbium atoms [25], the reader is referred to Ref. [26].
For quantum chaotic systems, the universal 1/ω k law for the average power spectrum in the frequency domain T * /T H ω k ≪ 1 can be read out from the existing RMT literature. Indeed, defining a set of discrete Fourier coefficients of level displacements {δε ℓ }, one observes the relation Statistics of the Fourier coefficients {a k } were studied in some detail [27] within the Dyson's Brownian motion model [28]. In particular, it is known that, in the limit k ≪ N , they are independent Gaussian distributed random variables with zero mean and the variance var[a k ] = N/(2π 2 βk). This immediately implies in concert with numerical findings. For larger k (in particular, for k ∼ N ), fluctuation properties of the Fourier coefficients {a k } are unknown. In view of the relation Eq. (1.8), a nonperturbative theory of the power spectrum to be developed in this paper sets up a well-defined framework for addressing statistical properties of discrete Fourier coefficients {a k } introduced in Ref. [27]. An attempt to determine S N (ω k ) for higher frequencies up to ω k = ω Ny was undertaken in Ref. [29] whose authors claimed to express the large-N power spectrum in the entire domain T * /T H ω k ≤ ω Ny in terms of the spectral form-factor [5] of a quantum system, τ ≥ 0. Referring interested reader to Eqs. (3), (8) and (10) of the original paper Ref. [29], here we only quote a small-ω k reduction of their result: (1.11) (Here, the hat-symbol (ˆ) is used to indicate that the power spectrumŜ N (ω k ≪ 1) is the one furnished by the form-factor approximation.) A similar approach was also used in subsequent papers [30,31]. Even though numerical simulations seemed to confirm a theoretical curve derived in Ref. [29], we believe that the status of their heuristic approach needs to be clarified. This will be done in Section 1.3.

Spectra with uncorrelated spacings: Form-factor vs power spectrum
A simple mathematical model of eigenlevel sequences {ε 1 , . . . , ε N } with identically distributed, uncorrelated spacings {s 1 , . . . , s N }, where ℓ-th ordered eigenlevel equals provides an excellent playing ground to analyze validity of the form-factor approximation. Defined by the covariance matrix of spacings of the form cov(s i , s j ) = σ 2 δ ij , such that s i = 1, it allows us to determine exactly both the power spectrum Eq. (1.4) and the form-factor Eq. (1.10).
Power spectrum.-Indeed, realizing that the covariance matrix of ordered eigenlevels equals we derive an exact expression for the power spectrum (N ∈ N) (1.14) Equation (1.14) stays valid in the entire region of frequencies 0 ≤ ω ≤ π. For a set of discrete frequencies ω k = 2πk/N , it reduces to (1.15) Remark 1.3. Notice that Eqs. (1.14) and (1.15) for the power spectrum of eigenlevel sequences with uncorrelated level spacings hold universally. Indeed, both expressions appear to be independent of a particular choice of the level spacings distribution; the level spacing variance σ 2 is the only model-specific parameter. (1.14) with σ 2 = 1. Blue crosses represent the average power spectrum simulated for 10 million sequences of N = 2048 random eigenlevels with uncorrelated, exponentially distributed spacings s i ∼ Exp(1). Inset: a log-log plot for the same graphs.
For illustration purposes, in Fig. 1, we compare the theoretical power spectrum S N (ω), Eq. (1.14), with the average power spectrum simulated for an ensemble of sequences of random eigenlevels with uncorrelated, exponentially distributed level spacings s i ∼ Exp(1). Since the unit mean level spacing s j = 1 is intrinsic to the model, the unfolding procedure is redundant. Perfect agreement between the theoretical and the simulated curves is clearly observed in the entire frequency domain 0 < ω ≤ π.
For further reference, we need to identify three scaling limits of S N (ω) that emerge as N → ∞. In doing so, the power spectrum will be multiplied by ω 2 to get rid of the singularity at ω = 0.
(i) The first -infrared -regime, refers to extremely small frequencies, ω ∼ N −1 . It is described by the double scaling limit where Ω = O(N 0 ). One observes: (1.17) (ii) The second scaling regime describes the power spectrum for intermediately small frequencies ω ∼ N −α with 0 < α < 1. In this case, a double scaling limit becomes trivial: , (1.19) where ω = O(N 0 ). One observes: Spectral form-factor.-For eigenlevel sequences with identically distributed, uncorrelated level spacings, the form-factor K N (τ ) defined by Eq. (1.10) can be calculated exactly, too. Defining the characteristic function of i-th level spacing, where f si (s) is the probability density of the i-th level spacing, we reduce Eq. (1.10) to (1. 22) In Fig. 2, we compare the theoretical form-factor Eq. (1.22) with the average formfactor simulated for an ensemble of sequences of random eigenlevels with uncorrelated, exponentially distributed level spacings as explained below Remark 1.3. The simulation was based on Eqs. (1.10) and (1.12), and involved averaging [32] over ten million realizations. Referring the reader to a figure caption for detailed explanations, we plainly notice a perfect agreement between the simulations and the theoretical result Eq. (1.22).
As N → ∞, three different scaling regimes can be identified for the spectral formfactor. Two of them, arising in specific double scaling limits, appear to be universal.
(i) The first -infrared -regime, refers to extremely short times, τ ∼ N −1 . Assuming existence and convergence of the moment-expansion for the characteristic function Ψ s (τ ), we expand it up to the terms of order N −2 , to derive the infrared double scaling limit for the form factor: where T = O(N 0 ). Notice that this formula holds universally as K (−1) (T ) does not depend on a particular choice of the level spacings distribution; its variance σ 2 is the only model-specific parameter. One observes: we discover that, for intermediately short times, the double scaling limit of the form factor reads where T = O(N 0 ). Hence, in the intermediate double scaling limit, the form-factor exhibits the universal behavior too, as K (−1/2) (T) depends on a particular choice of the level spacings distribution only through its variance σ 2 . One observes: The third scaling regime describes the form-factor for τ = O(N 0 ) fixed as N → ∞. Spotting that in this case the characteristic function Ψ N s (τ ) vanishes exponentially fast, we derive . (1.29) Notably, in the fixed-τ scaling limit, the form-factor is no longer universal as it depends explicitly on the particular distribution of level spacings § through its characteristic function Ψ s (τ ). One observes: (1.30) § For the exponential distribution of level spacings the form-factor in the third scaling regime equals unity, K (0) (τ ) ≡ 1. The three scaling regimes for the form-factor as N → ∞ are illustrated in Fig. 3. The continuity of the entire curve is guaranteed by equality of limits lim T →∞ K (−1) (T ) = lim T→0 K (−1/2) (T) and lim T→∞ K (−1/2) (T) = lim τ →0 K (0) (τ ), see Eqs. (1.25), (1.28) and (1.30). To highlight occurrence of both universal and non-universal τ -domains in the form-factor, the latter is plotted for three different choices of level spacing distributions, s j ∼ Erlang(3, 3), IG(1, 3) and U(0, 2), characterized by the same mean s j = 1 and the variance σ 2 = 1/3: (1.31) The three curves coincide in the universal domains (I) and (II). On the contrary, in the third regime [(III)], the form-factor behavior is non-universal as the three curves evolve differently depending on a particular choice of the level spacing distribution. Yet, all three curves approach unity at infinity.
Implications for the power spectrum.-We now turn to the discussion of a relation between the power spectrum Eq. (1.4) and the form-factor Eq. (1.10). To this end, we shall compare the limiting forms, as N → ∞, of the form-factor, studied both analytically and numerically in the previous subsection, with the limiting behavior of the product ω 2 S N (ω) | ω=2πτ as prompted by the form-factor approximation Eq. (1.11). The latter is plotted in Fig. 3 (1.34) Hence, the two spectral statistics -the form-factor and the power spectrum -cannot be reduced to each other for any finite frequency 0 < ω < π as N → ∞.
Conclusion.-Detailed analytical and numerical analysis, performed for eigenlevel sequences with uncorrelated, identically distributed level spacings, leads us to conclude that the spectral form-factor and the power spectrum are generically two distinct statistical indicators. This motivates us to revisit the problem of calculating the power spectrum for a variety of physically relevant eigenlevel sequences beyond the form-factor approximation. In the rest of the paper, this program, initiated in our previous publication [33], will be pursued for (a) generic eigenlevel sequences possessing stationary level spacings and (b) eigenlevel sequences drawn from a variant of the circular unitary ensemble of random matrices. The latter case is of special interest as its N → ∞ limit belongs to the spectral universality class shared by a large class of quantum systems with completely chaotic classical dynamics and broken time-reversal symmetry.

Main results and discussion
In this Section, we collect and discuss the major concepts and results of our work. Throughout the paper, we shall deal with eigenlevel sequences posessing stationary level spacings as defined below.
Remark 2.2. While stationarity of level spacings is believed to emerge after unfolding procedure in the limit N → ∞, see Ref. [34], it is not uncommon to observe stationarity even for finite eigenlevel sequences. Two paradigmatic examples of finite-N eigenlevel sequences with stationary spacings include (i) a set of uncorrelated identically distributed eigenlevels [35] mimicking quantum systems with integrable classical dynamics and (ii) eigenlevels drawn from the 'tuned' circular ensembles of random matrices appearing in the random matrix theory approach to quantum systems with completely chaotic classical dynamics, see Section 4.

Main results
First result.-For generic eigenlevel sequences, the power spectrum Eq. (1.4) is determined by both diagonal and off-diagonal elements of the covariance matrix δε ℓ δε m . In the important case of eigenlevel sequences with stationary level spacings, the power spectrum can solely be expressed in terms of diagonal elements δε 2 ℓ of the covariance matrix. The Theorem 2.3 below, establishes an exact relation between the power spectrum (see Definition 1.1) and a generating function of variances of ordered eigenvalues. Theorem 2.3 (First master formula). Let N ∈ N and 0 ≤ ω ≤ π. The power spectrum for an eigenlevel sequence {0 ≤ ε 1 ≤ · · · ≤ ε N } with stationary spacings equals where z = e iω , ∆ is the mean level spacing, and For the proof, the reader is referred to Section 3.2.
Second result.-Yet another useful representation -the second master formulaestablishes an exact representation of the power spectrum in terms of a generating function of probabilities E N (ℓ; ǫ) to observe exactly ℓ eigenlevels below the energy ε, (2.5) Here, P N (ǫ 1 , . . . , ǫ N ) is the joint probability density (JPDF) of N unordered eigenlevels taken from a positive definite spectrum; it is assumed to be symmetric under a permutation of its arguments. Such an alternative albeit equivalent representation of the power spectrum will be central to the spectral analysis of quantum chaotic systems.
Theorem 2.4 (Second master formula). Let N ∈ N and 0 ≤ ω ≤ π, and let Φ N (ε; ζ) be the generating function of the probabilities defined in Eq. (2.5). The power spectrum, Definition 1.1, for an eigenlevel sequence with stationary spacings equals where z = 1 − ζ = e iω , ∆ is the mean level spacing, and For the proof, the reader is referred to Section 3.3.
Remark 2.5. Notably, representations Eqs. (2.6) and (2.7) suggest that the power spectrum is determined by spectral correlation functions of all orders. Contrary to the spacing distribution, which is essentially determined by the gap formation probability [5] E N (0; ε), the power spectrum depends on the entire set of probabilities E N (ℓ; ε) with ℓ = 0, 1, . . . , N .
Third result.-To study the power spectrum in quantum systems with broken timereversal symmetry and completely chaotic classical dynamics, let us consider the tuned circular unitary ensemble (TCUE N ). Obtained from the traditional circular unitary ensemble CUE N +1 [5] by conditioning its lowest eigen-angle to stay at zero, the TCUE N is defined by the joint probability density of N eigen-angles {θ 1 , . . . , θ N } of the form Such a seemingly minor tuning of CUE N +1 to TCUE N induces stationarity of level spacings in TCUE N for any N ∈ N, see Corollary 4.3 for the proof. We note in passing that traditional circular unitary ensemble lacks the stationarity property. For the TCUE N , a general Definition 1.1 of the power spectrum can be adjusted as follows.
Definition 2.6. Let {θ 1 ≤ · · · ≤ θ N } be fluctuating ordered eigen-angles drawn from the TCUE N , N ∈ N, with the mean level spacing ∆ and let δθ ℓ δθ m be the covariance matrix of eigen-angle displacements δθ ℓ = θ ℓ − θ ℓ from their mean θ ℓ . A Fourier transform of the covariance matrix is called the power spectrum of the TCUE N . Here, the angular brackets denote average with respect to the JPDF Eq. (2.9).
Theorem 2.7 (Power spectrum in TCUE N ). Let {θ 1 ≤ · · · ≤ θ N } be fluctuating ordered eigen-angles drawn from the TCUE N . Then, for any N ∈ N and all 0 ≤ ω ≤ π, the power spectrum admits exact representation (2.14) Here, z = 1 − ζ = e iω whilst the functionσ N (t; ζ) is a solution to the Painlevé VI equation For the proof of Theorem 2.7, the reader is referred to Section 4.2.
Remark 2.8. Theorem 2.7 provides an exact RMT solution for the power spectrum in the TCUE N . Alternatively, but equivalently, the finite-N power spectrum can be expressed in terms of a Fredholm determinant (Section 4.3), Toeplitz determinant (Section 4.4) and discrete Painlevé V (dP V ) equations (Appendix B). While the Toeplitz representation is beneficial for performing a large-N analysis of the power spectrum, the dP V formulation is particularly useful for efficient numerical evaluation of the power spectrum for relatively large values of N .
Fourth (main) result.-The most remarkable feature of the random matrix theory is its ability to predict universal statistical behavior of quantum systems. In this context, a large-N limit of the power spectrum in the TCUE N is expected to furnish a universal, parameter-free law, S ∞ (ω) = lim N →∞ S N (ω), for the power spectrum. Its functional form is given in the Theorem 2.9 below.
Remark 2.10. As a by-product of this Theorem, we have formulated a conjecture for a double integral identity involving a fifth Painlevé transcendent. A mathematicallyoriented reader is referred to Conjecture 5.9.
Theorem 2.11 (Small-ω expansion). In the notation of Theorem 2.9, the following expansion holds as ω → 0: For the proof of Theorems 2.9 and 2.11, the reader is referred to Section 5.

Discussion
In Figs. 4 and 5, the parameter-free prediction Eq. (2.17) for the power spectrum is confronted with the results of numerical simulations for the large-N circular unitary ensemble CUE N . Two remarks are in order. (i) First, the limiting curve for S ∞ (ω) was approximated by the exact Painlevé VI solution computed for sufficiently large N through its dPV representation worked out in detail in Appendix B. We have verified, by performing numerics for various values of N , that the convergence of dP V representation of S N (ω) to S ∞ (ω) is very fast, so that the N = 10 4 curve provides an  5) lends an independent support to validity of our numerical procedure. (ii) Second, even though the theoretical results used for comparison refer to the TCUE Nrather than the CUE N -ensemble (which differ from each other by the weight function and the way the two are intrinsically unfolded ), the agreement between the TCUE N theory and the CUE N numerics is nearly perfect, which can naturally be attributed to the universality phenomenon emerging as N → ∞. The universal formula for S ∞ (ω), stated in Theorem 2.9, is the central result of the paper. We expect it to hold universally for a wide class of random matrix models belonging to the β = 2 Dyson's symmetry class, as the matrix dimension N → ∞. Expressed in terms of a fifth Painlevé transcendent, the universal law Eq. (2.17) can be viewed as a power spectrum analog of the Gaudin-Mehta formula Eq. (1.2) for the level spacing distribution.
Apart from establishing an explicit form of the universal random-matrix-theory law for S ∞ (ω), our theory reveals two important general aspects of the power spectrum which hold irrespective of a particular model of eigenlevel sequences: (i) similarly to the level spacing distribution, the power spectrum is determined by spectral correlations of all orders; (ii) in distinction to the level spacing distribution, which can solely be expressed in terms of the gap formation probability, the power spectrum is contributed by the entire set of probabilities that a spectral interval of a given length contains The spectra in CUE N and TCUE N ensembles are intrinsically unfolded for any N ∈ N, albeit each in its own way. Indeed, in the CUE N the mean density is a constant [5,6], while in the TCUE N the mean level spacing is a constant, see Corollary 4.3. In the limit N → ∞, the two types of unfolding are expected to become equivalent. exactly ℓ eigenvalues with ℓ ≥ 0. As such, it provides a complementary statistical description of spectral fluctuations in stochastic spectra of various origins.
Considered through the prism of Bohigas-Giannoni-Schmit conjecture, the universal law Eq. (2.17) should hold for a variety of quantum systems with completely chaotic classical dynamics and broken time-reversal symmetry at not too low frequencies T * /T H ω ≤ π, when ergodicity and global symmetries -rather than system specific features -are responsible for shaping system's dynamics.
Potential applicability of our results to the non-trivial zeros of the Riemann zeta function deserves special mention. Indeed, according to the Montgomery-Odlyzko law (see, e.g., Ref. [36]), the zeros of the Riemann zeta function located high enough along the critical line are expected to follow statistical properties of the eigenvalues of large U (N ) matrices. This suggests that the universal law Eq. (2.17) could be detected "experimentally". Extensive, high-precision data accumulated by A. M. Odlyzko for billions of Riemann zeros [37] provide a unique opportunity for a meticulous test of the new universal law.

Power spectrum for eigenlevel sequences with stationary spacings
In this Section, we provide proofs of two master formulae given by Theorem 2.3 and Theorem 2.4.

Stationary spectra
In view of Definition 2.1, we first establish a necessary and sufficient condition for eigenlevel sequences to posses stationarity of level spacings.
Lemma 3.1. For N ∈ N, let {0 ≤ ε 1 ≤ · · · ≤ ε N } 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.

Proof of Theorem 2.3
It follows from Eq. (3.1) of Lemma 3.1 written in the form Here, z k = e iω k and the derivative with respect to z k should be taken as if z k were a continuous variable.

Proof of Theorem 2.4
To prove the Theorem 2.4, we need the following Lemma: For N ∈ N, let {ε 1 ≤ · · · ≤ ε N } be an ordered sequence of eigenlevels supported on the half axis (0, ∞), and let E N (ℓ; ε) be the probability to find exactly ℓ eigenvalues below the energy ε, given by Eq. (2.5), with ℓ = 0, 1, . . . , N . The following relation holds: Here, p ℓ (ε) is the probability density of ℓ-th ordered eigenlevel where p 0 (ε) = p N +1 (ε) = 0 for ε > 0. Equivalently, Proof. Differentiating Eq. (2.5) and having in mind that the probability density of ℓ-th ordered eigenvalue equals Integration by parts performed in the last line is justified provided an average number of eigenlevels N N (ε) in the tail region (ε, ∞) exhibits sufficiently fast decay N N (ε) ∼ ε −(2+δ) with δ > 0 as ε → ∞ ¶. (3.10) Here, z k = e iω k and the derivative with respect to z k should be taken as if z k were a continuous variable.

Power spectrum in the tuned circular unitary ensemble
In this Section, a general framework developed in Section 3 and summed up in Theorem 2.4 will be utilized to determine the power spectrum in the tuned circular ensemble of random matrices, TCUE N , for any N ∈ N. For the definition of TCUE N , the reader is referred to Eqs. (2.9) and (2.10).

Correlations between ordered eigen-angles in TCUE N
The main objective of this subsection is to establish stationarity of spacings between ordered TCUE N eigen-angles. To this end, we prove Lemma 4.1 and Lemma 4.2. The sought stationarity will then be established in Corollary 4.3.
Proof. The proof is based on the circular-symmetry identity between the probability density functions of ℓ-th and (N − ℓ + 1)-th ordered eigenangles in the TCUE N . This relation can formally be derived from the representation Letting ε → ∞, we generate a large-ε expansion of the form is the ℓ-point spectral correlation function. To the first order, the expansion brings Φ N (ε; 1 − z) = z N + z N−1 (1 − z)N N (ε) + . . ., where N N (ε) is the mean spectral density R 1,N (ǫ) integrated over the interval (ε, ∞). Hence, the required decay of N N (ε) at infinity readily follows.
Proof. Indeed, combining Lemma 4.1 taken at q = 1 and Lemma 4.2 taken at q = 1 and m = ℓ − 1, one concludes that the mean spacing is the generating function of the probabilities to find exactly ℓ eigen-angles in the interval (0, ϕ) of the TCUE N spectrum. The JPDF P N (θ 1 , . . . , θ N ) is defined in Eq. (2.9). Substituting Eqs. (4.16) and (2.9) into Eq. (4.15), one derives a multidimensionalintegral representation of the generating function Φ N (ϕ; ζ) in the form satisfying the symmetry relation (4.18) Multidimensional integrals of the CUE-type akin to Eq. (4.17) have been studied in much detail in Ref. [38] whose authors employed the τ -function theory [39] of Painlevé equations. To proceed with evaluation of the generating function of our interest, we introduce a new set of integration variables to write down the generating function Eq. (4.17) in the form (4.20) Its Painlevé VI representation can be read off from Ref. [38] where z ′ k = e iω ′ k . This is essentially Eq. (17) previously announced in our paper Ref. [33].

Power spectrum in TCUE N as a Fredholm determinant
To derive a Fredholm determinant representation of the TCUE N power spectrum, a determinantal structure [5,6] of spectral correlation functions in the TCUE N should be established. This is summarized in Lemma 4.5 below.   where the TCUE N scalar kernel is expressed in terms of the sine-kernel S N +1 (θ) = sin[(N + 1)θ/2)] sin(θ/2) (4.25) of the CUE N +1 ensemble.
Proof. While the determinantal form [Eq. (4.23)] of spectral correlation functions is a universal manifestation of the β = 2 symmetry of the circular ensemble [5,6], a precise form of the two-point scalar kernel κ N (θ, θ ′ ) depends on peculiarities of the TCUE N probability measure encoded in the weight function (z = e iθ ) characterizing the TCUE N measure in Eq. (2.9). For aesthetic reasons, it is convenient to compute a scalar kernel κ N (θ, θ ′ ) in terms of polynomials {ψ j (z)} orthonormal on the unit circle |z| = 1 with respect to the weight function W (z). In such a case, a scalar kernel is given by either of the two representations (w = e iθ ′ ): This being said, we would like to represent the TCUE N scalar kernel in a more suggestive form. To do so, we notice that Szegö-Askey polynomials Eq. (4.31) admit yet another representation ψ ℓ (z) = 2 (ℓ + 1)(ℓ + 2) ℓ+1 j=1 jz j−1 . (4.33) Substituting it further into Eq. (4.28), one obtains: Owing to the representation of the CUE N sine-kernel the above can further be reduced to Calculating derivatives therein, we derive where S N +1 (θ) is the CUE N +1 sine-kernel: Proof. Equation (4.39) is self-evident as the determinant therein is the (ℓ + 1)-point correlation function in the CUE N +1 with one of the eigen-angles conditioned to stay at zero whilst the denominator is the CUE N +1 mean density S N +1 (0) = N + 1.  42) whilst κ N is the TCUE N two-point scalar kernel specified in Lemma 4.5.
A Fredholm determinant representation of the power spectrum is particularly useful for asymptotic analysis of the power spectrum in the deep 'infrared' limit ω ≪ 1 when ζ = 1 − z ≪ 1.

Power spectrum in quantum chaotic systems: Large-N limit
In the limit N → ∞, the exact solution for the TCUE N power spectrum should converge to a universal law. To determine it, we shall perform an asymptotic analysis of the exact solution Eqs. (2.12) and (2.13), stated in Theorem 2.7, with the generating function Φ N (ϕ; ζ) being represented as a Toeplitz determinant specified in Proposition 4.9.

Uniform asymptotics of the Toeplitz determinant
To perform the integral in Eq. (2.12) in the limit N → ∞, uniform asymptotics of the Toeplitz determinant Eq. (4.44) are required in the subtle regime of two merging singularities. In our case, one singularity is of a root type while the other one is of both root and jump types. Relevant uniform asymptotics were recently studied in great detail by Claeys and Krasovsky [46] who used the Riemann-Hilbert technique. Two different, albeit partially overlapping, asymptotic regimes in ϕ can be identified.
Asymptotics at the 'left edge'.-Defining the left edge as the domain 0 ≤ ϕ < ϕ 0 , where ϕ 0 is sufficiently small * , the following asymptotic expansion holds uniformly as N → ∞ (see Theorems 1.5 and 1.8 in Ref. [46]): Hereω = ω/2π is a rescaled frequency so that z = 1 − ζ = e 2iπω . The function σ(s) is the fifth Painlevé transcendent defined as the solution to the nonlinear equation The above holds for 0 ≤ω < 1/2.
Here, Gω is a known function ofω with G(· · ·) being the Barnes' G-function. The above holds for 0 ≤ω < 1/2. The leading term in Eqs. (5.7) and (5.8) is due to Ehrhardt [48]. ♯ Notice that, in distinction to Ref. [46], we kept two reminder terms in Eq. (5.4) -oscillatory and non-oscillatory, even though the latter term is subleading. The reason for this is that the function σ(s) will subsequently appear in the integral Eq. (5.11) which will make the non-oscillatory reminder term dominant.

Asymptotic analysis of the main integral
In doing the large-N asymptotic analysis of our exact solution for the power spectrum [Eqs. (2.12) and (2.14)], we shall encounter a set of integrals where k is a non-negative integer and Φ N (ϕ; ζ) is given by Eq. (4.17). We shall specifically be interested in k = 0 and 1.
Lemma 5.3. In the notation of Eq. (5.12), we have: Proof. To compute the integral Eq. (5.12) at k = 0, we invoke the expansion Eq. (4.15) of Φ N (ϕ; ζ) in terms of probabilities E N (ℓ; ϕ) of observing exactly ℓ eigenangles of TCUE N in the interval (0, ϕ), 14) The integral above can readily be calculated by performing integration by parts: The integral I N,k .-Unfortunately, exact calculation of the same ilk is not readily available for I N,k with k = 1. For this reason we would like to gain an insight from Eq. (5.13) as N → ∞, which, eventually, is the limit we are mostly concerned with. To this end, we extract the leading order behavior of I N,0 (ζ) on the unit circle |z| = |1 − ζ| = 1, The contribution from the bulk of the integration domain appears to be negligible due to strong oscillations e iωϕN of the integrand therein, see Eq. (5.8).
Equipped with these observations, we shall now proceed with an alternative, large-N , analysis of I N,k (ζ) for k = 0 and k = 1, where terms of the same structure (with and without strongly oscillating prefactor) will appear. Aimed at the analysis of the power spectrum [Eq. (2.12)], whose representation contains a very particular z-operator, we shall only be interested in the leading order contributions to both terms. Notably, even though for k = 1 a non-oscillating term is subleading as compared to an oscillating term, we shall argue that its contribution should still be kept.
To proceed with the large-N analysis of I N,k , we first rewrite the integral Eq. (5.12) as a sum of two such that and Here, Φ E N (ϕ; ζ) is an arbitrary integrable function; it will be specified later on. Prompted by the 'edge' and 'bulk' asymptotic expansions of Φ N (ϕ; ζ) [Eqs. (5.2) and (5.8)], we split the integral in Eq. (5.20) into three pieces correspondingly.
To facilitate the asymptotic analysis, we would ideally like to choose Φ E N (ϕ; ζ) in such a way that the contribution of the 'bulk' integral C N,k (ζ) will be dominated by the contributions coming from the 'left-edge' [L (1) N,k (ζ)] and the 'right-edge' [R (1) N,k (ζ)] parts of the integration domain. In fact, the contributions of the left and the right edges are related to each other; an exact relation between the two will be worked out and made explicit later on.
The integral I To get rid of N in the integral over the Painlevé V transcendent, we make the substitution λ = N ϕ to rewrite L N,k (ζ) in the form Further, choosing Ω(N ) to be a slowly growing function, Ω(N ) = ln N , one readily verifies that the third error term in Eq. (5.30) is a dominant one out of the three as 0 <ω < 1/2. Yet, it is smaller as compared to the integral in Eq. (5.30) by a factor Ω(N ) k−1−2ω 2 that tends to zero as N → ∞. Thus, in the leading order, we derive: with k = 0 and 1. Now, let us turn to the 'right-edge' integral R (1) N,k (ζ). Due to the symmetry relation Eq. (4.18) shared by Φ E N (ϕ; ζ) too, we realize that the contributions of the left and the right edges are related to each other: Considering the integral in the r.h.s. of Eq. (5.33) along the lines of the previous analysis, we conclude that the following formula holds as N → ∞: (5.36) The notation → was used here to stress that the r.h.s. contains each leading order contribution of both terms, the oscillating and the non-oscillating, as discussed in the paragraph prior to Eq. (5.19).
The integral I N,k (ζ).-As soon as the function Φ E N (ϕ; ζ) contains a strongly oscillating factor e iωϕN , the integral I (2) N,k (ζ) in Eq. (5.21) can be calculated by the stationary phase method [49]. Since there are no stationary points within the interval (0, 2π), the integral is dominated by contributions L and Proof. Apply the stationary phase method [49] to calculate the integral Eq. (5.21).
The integral I N,k (ζ).-The calculation above implies that the main integral of our interest admits an asymptotic representation with k = 0, 1 and N,k (ζ) negligible. If this is not the case, one should replace Φ E N with someΦ E N by adding to Φ E N the higher-order corrections (up to O(N −2 )) that can be obtained from the full asymptotic expansion of Φ N (ϕ; ζ), see Remark 1.4 of Ref. [45]. Inclusion of these higher-order corrections will reduce the contribution of C On the other hand, the proposed modification of Φ N (ϕ; ζ) will produce corrections to the functions L

Proof of Theorem 2.9
Now we are in position to evaluate the power spectrum as N → ∞. To proceed, we start with the exact, finite-N , representation     for ω = π, one would have to use the Theorem 1.12 of Ref. [46] instead of Theorems 1.5, 1.8 and 1.11 of the same paper.
Since numerical calculations indicate that the power spectrum is continuous at ω = π, we did not study this case analytically.

Proof of Theorem 2.11
Below, the universal law S ∞ (ω) for the power spectrum will be studied in the vicinity of ω = 0. In the language of S N (ω) this corresponds to performing a small-ω expansion after taking the limit N → ∞. Equation (2.17) will be the starting point of our analysis.
Preliminaries.-Being interested in the small-ω behavior of the power spectrum Eq. (2.17), we observe that the functions A(ω) and B(ω), defined by Eqs. (2.18) and (2.19), admit the expansions (5.56) so that the power spectrum, asω → 0, can be written as Here,Λ(ω) denotes a small-ω expansion of the function Small-ω ansatz for the fifth Painlevé transcendent.-To proceed, we postulate the following ansatz for a small-ω expansion of the fifth Painlevé function σ 1 (t): Here, the functions f k (t) with k = 1, 2, . . . satisfy the equations The third order differential equation (5.61) can be solved to bring This representation assumes that the integrals are convergent; integration constants have to be fixed by the boundary conditions Eqs. (5.66), (5.67), (5.68), etc. In particular, we derive: is the Euler's constant, and Here, the function F 3 (t) is known explicitly from Eqs. (5.64), (5.70) and (5.71). We notice that Representation ofΛ(ω) as a partial sum.-Having determined the functions f 1 (t), f 2 (t) and f 3 (t), we now turn to the small-ω analysis of Λ(ω) [Eq. (5.58)]. Expanding the expression in square brackets in small ω, we obtain: The functions G k (λ) can be evaluated explicitly in terms of integrals containing f k (λ) defined in Eq. (5.60). For example, Here, and and, hence, Substituting Eq. (5.74) into Eq. (5.58), we split Λ(ω) into a partial sum Performing the integral, we obtain: Estimate of Λ k (ω).-To treat Λ k (ω) for k ≥ 2, we split it into two parts Here, T is an arbitrary positive number to be taken to infinity in the end. Since a small-ω expansion of A k (ω, T ) is well justified for any finite T , see e.g. Eq. (5.92) below, we conclude that To estimate B k (ω, T ), we refer to Remark 5.7 which implies that G k (λ) = O(λ −1 ) as λ → ∞. Replacing G k (λ) with 1/λ in Eq. (5.88), we perform the integration by parts twice in the resulting integral to conclude that it is of order O(ω −1 ). This entails Since we are interested in calculating Λ(ω) up to the terms O(ω 3 ), see Eq. (5.59), we need to consider A k (ω, T ) and B k+1 (ω, T ) for k ≤ 2 only.
Calculation ofΛ 2 (ω).-A small-ω expansion of A 2 (ω, T ) brings Since G 2 (λ) ∈ R, we even conclude that For this reason, A 2 (ω, T ) does not contribute toΛ 2 (ω). Evaluation of B 2 (ω, T ), given by is more involved. A simplification comes from the fact that, at some point, we shall let T tend to infinity. For this reason, it suffices to consider a large-λ expansion of G 2 (λ) in the integrand. Straightforward calculations bring To determine a small-ω expansion of B 2 (ω, T ), we shall further concentrate on small-ω expansions of its constituents, I 0 (ω, T ), I ±1 (ω, T ) and K ±1 (ω, T ).
(a).-The function I 0 (ω, T ) can be evaluated exactly, Expanding this result aroundω = 0 we derive To analyze a small-ω expansion we proceed in two steps. First, we determine the coefficient α 0 (j, T ) directly from Eq. (5.98) cos[(ω + j)λ] λ 2+2ω 2 (5.106) whose integral term possesses a better convergence whenω approaches zero, as compared to the one given by Eq. (5.98). Differentiating Eq. (5.106) with respect toω and settingω = 0 we derive: This implies the relation (j = 0) As a consequence, we conclude that we follow the same strategy. First, we determine the coefficient κ 0 (j, T ) directly from Eq. (5.99) (5.111) to observe the relation (j = 0) Second, to examine a linear term of a small-ω expansion of K j =0 (ω, T ), we perform integration by parts in Eq. (5.99) in order to improve integral's convergence: Differentiating Eq. (5.113) with respect toω and settingω = 0, we obtain: This implies the relation (j = 0) As a consequence, we conclude that where a 1 , a 2 and a 3 are real coefficients whose explicit values are not required for our analysis; sin(⋆ t) stands to denote sin t and sin(2t), both of which are present in the remainder term. This expansion combined with Eq. (5.73) implies the following large-t behavior of f 3 (t): Here, the coefficients a ′ j ∈ R are real. Now, a large-λ behavior of F 3 (λ) can be read off from Eq. (5.80): where the coefficients a ′′ j ∈ R are real, again. Inspection of Eqs. (5.77), (5.95), (5.96) and (5.121) shows that a large-λ behavior of G 3 (λ) coincides with that of F 3 (λ).
(a).-The first integral, originating from the a ′′ 1 term in Eq. (5.121), admits a small-ω expansion  A straightforward calculation produces a small-ϕ expansion of Φ N (ϕ; ζ) whose several initial terms read: Remark A.1. Since the above procedure is capable of producing the expansion coefficients σ j (N, ζ) of any finite order, it can also be utilized -by virtue of Eq. (2.14) -to generate a small-ϕ expansion of Φ N (ϕ; ζ) up to required accuracy.
B. Generating function Φ N (ϕ; ζ) and discrete Painlevé V equations (dP V ) To avoid intricacies [52] of a numerical evaluation of the six Painlevé functionσ N (t; ζ) appearing in the generating function Eq. (2.14), we opt for an alternative representation of Φ N (ϕ; ζ) in terms of discrete Painlevé V equations. To proceed, we follow Ref. [53] (see also Refs. ( [54,55,56])), to observe that a sequence of U (N ) integrals Here, r N andr N are so-called reflection coefficients appearing in the Szegö theory [57] of orthogonal polynomials on the unit circle. Remarkably, there exists the N -recurrence for reflection coefficients {r N ,r N } as specified in Proposition 4.4 in Ref. [53]; a variation of their Proposition is given below.
The Proposition yields a sought dP V representation of the generating function Φ N (ϕ; ζ), see Eq. where the reflection coefficients {r N ,r N } are determined by equations considered in conjunction with two systems of coupled first order discrete Painlevé equations (dP V ): and Remark B.2. To avoid numerical z-differentiation of Φ N (ϕ; 1 − z) appearing in the formula Eq. (2.12), it is beneficial to produce a similar system of coupled recurrence equations for (∂/∂z)Φ N (ϕ; 1 − z). Since the resulting recurrences are too cumbersome to state them here, we leave their (straightforward) derivation to the inquisitive reader.
Remark B.3. Away from the endpoints ϕ = 0 and ϕ = 2π, the dP V representation opens a way for effective numerical evaluation of both Φ N (ϕ; ζ) and (∂/∂z)Φ N (ϕ; 1−z) for finite N . Since the recurrence procedure tends to accumulate numerical errors, we have used quadruple precision numbers to achieve sufficient precision for very large N (e.g., for N = 10 4 , see Figs. 4 and 5).
Remark B.4. In the vicinity of the endpoints ϕ = 0 and ϕ = 2π, numerical precision of the above recurrence procedure worsens drastically since the recurrence equations start to exhibit a singular behavior. To circumvent this drawback at ϕ = 0, we have used a small-ϕ expansion of Φ N (ϕ; ζ) as described in the Remark A.1. In the vicinity of ϕ = 2π, the symmetry relation Eq. (4.18) combined with a small-ϕ expansion makes the job.