==== Front Commun Math Phys Commun Math Phys Communications in Mathematical Physics 0010-3616 1432-0916 Springer Berlin Heidelberg Berlin/Heidelberg 4692 10.1007/s00220-023-04692-y Article On the Spectral Form Factor for Random Matrices Cipolloni Giorgio gc4233@princeton.edu 1 http://orcid.org/0000-0001-5366-9603 Erdős László lerdos@ist.ac.at 2 Schröder Dominik dschroeder@ethz.ch 3 1 grid.16750.35 0000 0001 2097 5006 Princeton Center for Theoretical Science, Princeton University, Princeton, NJ 08544 USA 2 grid.33565.36 0000000404312247 IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria 3 grid.5801.c 0000 0001 2156 2780 Institute for Theoretical Studies, ETH Zurich, Clausiusstr. 47, 8092 Zurich, Switzerland 23 3 2023 23 3 2023 2023 401 2 16651700 11 10 2021 13 2 2023 © The Author(s) 2023 https://creativecommons.org/licenses/by/4.0/ Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. In the physics literature the spectral form factor (SFF), the squared Fourier transform of the empirical eigenvalue density, is the most common tool to test universality for disordered quantum systems, yet previous mathematical results have been restricted only to two exactly solvable models (Forrester in J Stat Phys 183:33, 2021. 10.1007/s10955-021-02767-5, Commun Math Phys 387:215–235, 2021. 10.1007/s00220-021-04193-w). We rigorously prove the physics prediction on SFF up to an intermediate time scale for a large class of random matrices using a robust method, the multi-resolvent local laws. Beyond Wigner matrices we also consider the monoparametric ensemble and prove that universality of SFF can already be triggered by a single random parameter, supplementing the recently proven Wigner–Dyson universality (Cipolloni et al. in Probab Theory Relat Fields, 2021. 10.1007/s00440-022-01156-7) to larger spectral scales. Remarkably, extensive numerics indicates that our formulas correctly predict the SFF in the entire slope-dip-ramp regime, as customarily called in physics. http://dx.doi.org/10.13039/501100000781 European Research Council 101020331 Erdős László http://dx.doi.org/10.13039/501100012652 ETH Zürich Foundation issue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2023 ==== Body pmcIntroduction Spectral statistics of disordered quantum systems tend to exhibit universal behavior and hence are widely used to study quantum chaos and to identify universality classes. In the chaotic regime, the celebrated Wigner–Dyson–Mehta eigenvalue gap statistics involving the well-known sine-kernel [42] tests this universality on the scale of individual eigenvalue spacing. On this small microscopic scale the universality phenomenon is the most robust and it depends only on the fundamental symmetry type of the model. On larger scales more details of the model influence the spectral statistics, nevertheless several qualitative and also quantitative universal patterns still prevail. The spectral form factor and predictions from physics In the physics literature the standard tool to investigate eigenvalues λ1,λ2,…,λN of a Hermitian N×N matrix (Hamiltonian) H on all scales at once is the spectral form factor (SFF) [40] defined as1.1 SFF(t):=1N2∑i,j=1Neit(λi-λj)=|⟨eitH⟩|2 with a real time parameter t>0, i.e. it is the square of the Fourier transform of the empirical spectral density. Here we denoted the normalized trace of any N×N matrix A by ⟨A⟩=1NTrA. In case of random H, the expectation of SFF(t) is denoted by1.2 K(t):=E[SFF(t)]. For typical disordered Hamiltonians a key feature of SFF(t) is that for larger t (more precisely, in the ramp and plateau regimes, see later) it is strongly dependent on the sample, i.e. the standard deviation of SFF(t) is comparable with K(t). In other words, SFF(t) is not self-averaging [45] despite the large summation in (1.1). The spectral form factor and its expectation K(t) have a very rich physics literature since they contain most physically relevant information about spectral statistics. Quantizations of integrable systems typically result in K(t)∼1/N for all t where N is the dimension of the Hilbert space. Chaotic systems give rise to a linearly growing behavior of K(t) for smaller t (so-called ramp) until it turns into a flat regime, the plateau. The turning point is around the Heisenberg time TH, but the details of the transition depend on the symmetry class of H and on whether the eigenvalues are rescaled to take into account the non-constant density of states (in physics terminology: unfolding the spectrum). For example, in the time irreversible case (GUE symmetry class) the unfolded SFF has a sharp kink, while in the GOE symmetry class the kink is smoothened. The exact formulas can be computed from the Fourier transform of the two point eigenvalue correlation function of the corresponding Gaussian random matrix ensemble, see [42, Eqs. (6.2.17), (7.2.46)], the result is1.3 KGUE(τTH)≈1N×τ,0<τ≤11,τ≥1,KGOE(τTH)≈1N×2τ-τlog(1+2τ),0<τ≤12-τlog2τ+12τ-1,τ≥1, for any fixed τ>0 in the large N limit. Here we expressed the physical time t in units of the Heisenberg time, τ=t/TH, where TH is given by TH=2πρ¯ with ρ¯ being the average density. Choosing the standard normalisation for the independent (up to symmetry) matrix elements,1.4 Ehij=0,E|hij|2=1N, the limiting density of states is the semicircle law ρsc(E)=12π(4-E2)+, so we have N eigenvalues in an interval of size 4, hence ρ¯=N/4 and thus TH=π2N. In particular, in the original t variable1.5 KGUE(t)≈2tπN2,δN≤t≤π2N1N,t≥π2N. Note the lower bound on t: the formula holds in the large N limit in the regime where t≥δN for some fixed δ>0 that is independent of N. The corresponding formulas without unfolding the spectrum (i.e. for the quantity defined in (1.1)) are somewhat different, see e.g. [9, Eq. (4.8)] for the GUE case; they still have a ramp-plateau shape but the kink is smoothened. The ramp-plateau picture and its sensitivity to the symmetry type has been established well beyond the standard mean field random matrix models. In fact, the Bohigas-Giannoni-Schmit conjecture [6] asserts that the formulas (1.3) are universal, i.e. they hold essentially for any chaotic quantum system, depending only on whether the system is without or with time reversal symmetry. The nonrigorous but remarkably effective semiclassical orbit theory [4, 31, 43, 48] based upon Gutzwiller’s trace formula [27] and many follow-up works verified this conjecture for quantizations of a large family of classical chaotic systems, e.g. for certain billiards.Fig. 1 A typical slope-dip-ramp-plateau picture for the spectral form factor of a chaotic system. The figure on log-log scale shows the SFF of a single GUE realisation H of size 500×500, as well as the empirical mean and standard deviation obtained from 500 independent realisations For smaller times, t≪TH, other details of H may become relevant. In particular the drop from K(t=0)=1 to K(t)≪1 for 1≪t≪TH is first dominated by the typical non-analyticity of the density of states at the spectral edges giving rise to the slope regime up to an intermediate minimum point of K(t), called the dip (in the early literature the dip was called correlation hole [40], for a recent overview, see [17]). Figure 1 shows the typical slope-dip-ramp-plateau picture for the GUE ensemble. Formula (1.5) is valid starting from scales t≫N1/2, while K(t) is oscillatorily decreasing for t≲N1/2 with a dip-time tdip∼N1/2. Thus K(t) follows the universal behavior (1.5) only for t≫tdip. In this regime the fluctuation of the SFF is comparable with its expectation, K(t), in fact ⟨eitH⟩ is approximately Gaussian. In contrast, the dominant contribution to the slope regime, t≪tdip, is self-averaging with a relatively negligible fluctuation. However, if the edge effects are properly discounted (e.g. by considering the circular ensemble with uniform spectral density on the unit circle), i.e. the slope regime is entirely removed, then the Gaussian behavior holds for all t≪TH with a universal variance given by (1.5). In more recent works spectral form factors were studied for the celebrated Sachdev-Ye-Kitaev (SYK) model [18, 23, 24, 32, 46] which also exhibits a similar slope-dip-ramp-plateau pattern although the details are still debated in the physics literature and the numerics are much less reliable due to the exponentially large dimensionality of the model. Our results Quite surprisingly, despite its central role in the physics literature on quantum chaos, SFF has not been rigorously investigated in the mathematics literature up to very recently, when Forrester computed the large N limit of K(t) rigorously for the GUE in [21] and the Laguerre Unitary Ensemble (LUE) in [22] in the entire regime t≪N. Both results rely on a remarkable identity from [9] (and its extension to the LUE case) and on previous stimulating work of Okuyama [44]. However, these methods use exact identities and thus are restricted to a few explicitly solvable invariant ensembles. The main goal of the current paper is to investigate SFF beyond these special cases with a robust method, the multi-resolvent local laws. While our approach is valid for quite general ensembles, for definiteness we focus on two models: the standard Wigner ensemble (for both symmetry classes) and the novel monoparametric ensemble introduced recently [25] by Gharibyan, Pattison, Shenker and Wells. The latter consists of matrices of the form Hs:=s1H1+s2H2, where H1 and H2 are typical but fixed realisations of two independent Wigner matrices and s=(s1,s2)∈S1⊂R is a continuous random variable. The normalization s12+s22=1 guarantees that the semicircle law for Hs is independent of s and it also shows that the model has effectively only one random parameter. One may also consider similar ensembles with finitely many parameters (see Remark 2.4) resulting in qualitatively the same behavior but with different power laws, see Table 1. We study the statistics of Hs in the probability space of the single random variable s and probe how much universality still persists with such reduced randomness. We write Es for the expectation wrt. s and EH, StdH for the expectation and standard deviation wrt. H1 and H2. Our main result is to prove a formula for the expectation and standard deviation of SFF for both ensembles up to an intermediate time. While this does not include the ramp regime, it already allows us to draw the following two main conclusions of the paper: The expectation and standard deviation of SFF(t) for Wigner and monoparametric ensembles exhibit the same universal behavior to leading order for 1≪t≪N1/4 if the trivial edge effects are removed. In the monoparametric case it is quite remarkable that already a single real random variable generates universality. For the monoparametric ensemble K(t)=Es[SFF(t)] depends non-trivially on the fixed H1,H2 matrices, but for large t this dependence is a subleading effect whose relative size becomes increasingly negligible as a negative power of t. In particular, while the speed of convergence to universality is much slower for the monoparametric ensemble than for the Wigner case, it is improving for larger t. The second item answers a question raised by the authors of [25] which strongly motivated the current work. In particular, sampling from s does not give a consistent estimator for K(t), but the relative precision of such estimate improves for larger times. We supplement these proofs with an extensive numerics demonstrating that both conclusions hold not only for t≪N1/4 but for the entire ramp regime, i.e. up to t≪TH∼N. Note that recently we have proved [15] that the Wigner–Dyson–Mehta eigenvalue gap universality holds for the monoparametric ensemble, which strongly supports, albeit does not prove, that K(t) in the plateau regime is also universal. We remark that our method applies without difficulty for finite temperatures (expressed by a parameter β>0) and for different-time autocorrelation functions, i.e. for⟨e(-β+it)H⟩⟨e(-β-it′)H⟩ as well, but for the simplicity of the presentation we focus on SFF(t) defined in (1.1), i.e. on β=0 and t=t′. Relations to previous mathematical results Rigorous mathematics for the spectral form factor, even for Wigner matrices or even for GOE, significantly lags behind establishing the compelling physics picture about the slope-dip-ramp-plateau. Given the recently developed tools in random matrix theory, it may appear surprising that they do not directly answer the important questions on SFF. We now briefly explain why. Limitations of the resolvent methods For problems on macroscopic spectral scales (involving the cumulative effect of order N many eigenvalues), and to a large extent also on mesoscopic scales (involving many more than O(1) eigenvalues), the resolvent method is suitable. This method considers the resolvent G(z)=(H-z)-1 of H for a spectral parameter z away from (but typically still close to) the real axis and establishes that in a certain sense G(z) becomes deterministic. This works for η=ℑz≫N-1 (in the bulk spectrum), i.e. on scales just above the eigenvalue spacing (note that the imaginary part of the spectral parameter sets a scale in the spectrum). Such results are called local laws and they can be extended to regular functions f(H) by standard spectral calculus (Helffer–Sjöstrand formula, see (3.3) later). However, the interesting questions about SFF concern a 1/N subleading fluctuation effect beyond the local laws. IndeedTreitH=∑ieitλi is a special case of the well-studied linear eigenvalue statistics, Trf(H)=∑if(λi), with the regular test function f(λ)=eitλ. To leading order it is deterministic and its fluctuation satisfies the central limit theorem (CLT) without the customary N normalisation, i.e.1.6 ∑if(λi)-E∑if(λi)≈N(0,Vf),withE∑if(λi)=N∫Rf(x)ρsc(x)dx+Of(1). is a normal random variable with variance.11.7 Vf=14π2∬-22|f(x)-f(y)x-y|24-xy4-x24-y2dxdy. The computation of higher moments of Trf(H)-ETrf(H) requires a generalization of the local laws to polynomial combinations of several G’s that are called multi-resolvent local laws. Applying (1.6)–(1.7) to f(x)=eitx we obtain, roughly,1.8 SFF(t)=1N2|TreitH|2≈[J1(2t)t+O(VfN2)]2,t≫1, using that∫Rf(x)ρsc(x)dx=∫Reitxρsc(x)dx=J1(2t)t, where J1 is the first Bessel function of the first kind. Note that Vf in (1.7) scales essentially as the H1/2 Sobolev norm of f hence Vf∼t for our f(x)=eitx in the regime t≫1. Therefore the size of the fluctuation term in (1.8) is Vf/N2∼t/N2 and it competes with the deterministic term (J1/t)2∼t-3. The dip time tdip∼N is obtained as the threshold where the fluctuation (the linear ramp function) becomes bigger than the slope function (J1/t)2. This argument, however, is heuristic as it neglects the error terms in (1.6) that also depend on t via f. CLT for linear statistics (1.6) for Wigner matrices H has been proven [1, 3, 13, 26, 28–30, 33, 34, 36, 38, 41, 47, 49] for test functions of the form f(x)=g(Na(x-E)) with some fixed reference point |E|<2, scaling exponent a∈[0,1) and smooth function g with compact support, i.e for macroscopic (a=0) and mesoscopic (00. The only CLT-type result for a special two-scale observable is in [2] where the eigenvalue counting function smoothed on an intermediate scale N-1/3 was considered. Quite remarkably, extensive numerics shows that the formulas (1.6)–(1.7) for f(x)=eitx are in perfect agreement with the expected behavior of K(t) in the entire slope-dip-ramp regime all the way up to t≪N, i.e. the CLT for linear statistics correctly predicts SFF well beyond its proven regime of validity. In the current paper we optimise the error terms specifically for eitx and thus we could cover the regime t≪N5/11 for the variance in (1.6) (corresponding to E[SFF(t)]). Limitations of Dyson Brownian motion techniques For the microscopic scale (i.e. comparable with the eigenvalue spacing, 1/N in the bulk) the resolvent is heavily fluctuating as it strongly depends on single eigenvalues. Local laws cannot access them, but in this regime another approach, the careful analysis of the Dyson Brownian Motion (DBM) becomes applicable. While these two approaches are complementary and apparently cover all scales, the actual methods require additional conditions that seriously restrict their use for SFF. The formulas (1.3) are obtained by computing the Fourier transform of the two point correlation function of the rescaled (unfolded) eigenvalues. Indeed, in the GUE case KGUE(t) in (1.3) is just the Fourier transform of p2(x,y)-1+δ(x-y) in the difference variable x-y, wherep2(x,y):=1-(sin(π(x-y))π(x-y))2, is the two point function, given by the celebrated Wigner–Dyson sine kernel, and KGOE(t) has a similar origin. Wigner–Dyson theory is designed for microscopic scales, i.e. to describe eigenvalue correlations on scales comparable with the local level spacing Δ, this is encoded in the fact that (1.3) holds for any fixed τ>0 in the N→∞ limit (equivalently that (1.5) holds only for t≥δN since Δ∼1/N in the bulk). While this is a very elegant argument supporting (1.3), mathematically it is quite far from a rigorous proof. The mathematical proofs of the sine-kernel universality use test functions that are rapidly decaying beyond scale Δ. The typical statements (so called fixed energy universality [7, 39]) show that for any fixed energy E in the bulk∑i0 the probability of the event is bigger than 1-N-D if N≥N0(D), for some N0(D)>0. Statement of the Main Results Our new results mainly concern the monoparametric ensemble but for comparison reasons we also prove the analogous results for the Wigner ensemble. We start with the two corresponding definitions. Definition 2.1 The Wigner ensemble consists of Hermitian N×N random matrices H with the following properties. The off-diagonal matrix elements below the diagonal are independent, identically distributed (i.i.d) real (β=1) or complex (β=2) random variables; in the latter case we assume that Ehij2=0. The diagonal elements are i.i.d. real random variables with Ehii2=2/(Nβ). Besides the standard normalisation (1.4), we also make the customary moment assumption: for every q∈N there is a constant Cq such that1 E|Nhij|q≤Cq. In the case of Gaussian distributions, it is called the Gaussian Orthogonal or Unitary Ensemble (GOE/GUE), for the real and complex cases, respectively. Remark 2.2 The assumptions Ehij2=0 in the complex case, and Ehii2=2/(βN) are made purely for convenience. All results can easily be generalised beyond this case but we refrain from doing so for notational simplicity. Definition 2.3 The monoparametric ensemble consists of Hermitian N×N random matrices of the form2.1 H=Hs:=s1H1+s2H2, where H1,H2 are independent Wigner matrices satisfying2E|hij(1)|4=E|hij(2)|4 and s=(s1,s2)∈S1 is a random vector, independent of H1,H2. On the distribution of s we assume that it has an square integrable density ρ(s) independent of N. We write Es for the expectation wrt. s and EH, StdH for the expectation and standard deviation wrt. the Wigner matrices H1 and H2. The parameter space S1⊂R2 inherits the usual scalar product and norms from R2, so for s,r∈S1 we have⟨s,r⟩:=s1r1+s2r2,‖s‖p:=(|s1|p+|s2|p)1/p. We also introduce the entrywise product of two vectors:s⊙r:=(s1r1,s2r2). For a fixed s, Hs is just the weighted sum of two Wigner matrices, and, due to the normalisation, itself is just a Wigner matrix. However, the concept of monoparametric ensemble views Hs as a random matrix in the probability space of the single random variable s for a typical but fixed (quenched) realization of H1 and H2. While Wigner matrices have a large (∼N2) number of independent random degrees of freedom, the monoparametric ensemble is generated by one single random variable hence, naively, much less universality properties are expected. Nevertheless, the standard Wigner–Dyson local eigenvalue universality holds [15]. Remark 2.4 In [15] we considered the un-normalized monoparametric model Hs:=H1+sH2, for some real valued random variable s, whose density of states is a rescaled semicircular distribution. In this paper we prefer to work with more homogeneous models since the formulas are somewhat nicer, but our main results also apply to this inhomogeous model with some slightly different exponents in the error terms. One may also consider a different un-normalized ensemble, s1H1+s2H2 with s∈R2 having an absolutely continuous distribution, which is effectively a two parameter model. Similar results also hold for the multi-parametric analogue of (22.2), i.e. s1H1+⋯+skHk for s∈Sk-1, see Remark 2.6 and Sect. 2.4 later. Despite all these options, for definiteness, the main body of this paper concerns the homogenous monoparametric model from Definition 2.3. Central limit theorem for sum of Wigner matrices To understand the effect of the random s, we study the joint statistics of Hs and Hr for two different fixed realisations r, s in the probability space of H1,H2, i.e. we aim at the correlation effects between Hs and Hr. We introduce the short-hand notations2 ⟨f⟩sc:=∫-22f(x)4-x22πdx,⟨f⟩1/sc:=∫-22f(x)1π4-x2dx,κ4:=N2E|h12|4-1-2β. To estimate the error term in the following theorem we introduce a parameter 1≤τ≪N and the weighted norm2.2 ‖f‖τ:=τ2‖f‖∞+τ‖f‖H1+‖f‖H2, where ‖f‖Hk2:=∑j≤k∫R|f(j)|2 is the usual Sobolev norm. For the applications later, the parameter τ will be optimized. Theorem 2.5 For s∈S1 and test functions f∈H2(R) the family of random variables Trf(Hs) is approximately Gaussian of mean2.3 ETrf(Hs)=N⟨f⟩sc+κ4‖s‖44⟨x4-4x2+22f⟩1/sc+1(β=1)[f(2)+f(-2)4-⟨f⟩1/sc2]+O(E1), and fluctuation2.4 E∏i=1p(Trfi(Hsi)-ETrfi(Hsi))=∑P∈Pair([p])∏(i,j)∈Pvsisj(fi,fj)+Op(Ep), for any fixed p∈N, functions f1,⋯,fp∈H2(R), and parameters s1,⋯,sp∈S1, where35 vsr(f,g):=1βπ2∬-22f′(x)g′(y)Vsr(x,y)dxdy+κ42⟨s⊙s,r⊙r⟩⟨(2-x2)f⟩1/sc2Vsr(x,y):=log|1-⟨s,r⟩msc(x)msc(y¯)|-log|1-⟨s,r⟩msc(x)msc(y)|. Here Ep are error terms which for any 1≤τ≪N and any ξ,ϵ>0 may be estimated by42.5 E1:=Nξ‖f‖τN1/2τ1/2,Ep:=Nξ1N1/2τ3/2+NϵN+N-ϵτ2p-11+τ2N1-2ϵ∏i∈[p]‖fi‖τ, for p≥2. Additionally, if s1=⋯=sp, i.e. in the Wigner case, we have the improved bound6 Ep:=1N1/2τ3/2∏i∈[p]‖fi‖τ and the first term of (2.7) for β=2 coincides with (1.7). We note that (2.7) generalizes the standard variance calculation yielding (1.7) to s≠r, see Sect. 3.2.4. Remark 2.6 Theorem 2.5 verbatim holds true also for the multi-parametric models1H1+⋯+skHk upon interpreting ⟨s,r⟩ and ‖s‖p as the Euclidean inner product and p-norm in Rk. Similarly, Theorem 2.5 also applies to the un-normalised case s∈R2 for which on the rhs. of (52.5) the function f has to be replaced by f(‖s‖·) with ‖·‖:=‖·‖2 and vsr from (2.7) has to be replaced by2.6 v~sr(f,g):=‖s‖‖r‖βπ2∬-22f′(‖s‖x)g′(‖r‖y)Vs‖s‖,r‖r‖(x,y)dxdy+κ42⟨s⊙s,r⊙r⟩⟨(2-x2)f(‖s‖x)⟩1/sc⟨(2-x2)f(‖r‖x)⟩1/sc. SFF for Wigner and monoparametric ensemble In this section we specialise Theorem 2.5 to the SFF case. We define the approximate expectation (rescaled by 1/N)2.7 eNs(t):=e(t)+1N[κ4‖s‖441-6t2J0(2t)+κ4‖s‖446t3-4tJ1(2t)-1(β=1)J0(2t)-cos(2t)2]e(t):=J1(2t)t in terms of the Bessel functions Jk of the first kind. We also define the approximate variance8 v±,κsr(t):=vsr(eit·,e±it·)=v±sr(t)+κ4⟨s⊙s,r⊙r⟩J2(2t)2,v±sr(t):=t2βπ2∬-22cos(t(x±y))Vsr(x,y)dxdy, From Theorem 2.5, choosing fi(x)=e±itx and τ=t, and recalling that ⟨e±itHs⟩=N-1Tre±itHs, we readily conclude the following asymptotics for SFF of the Wigner and monoparametric ensemble.Fig. 2 In the first plot we compare the empirical mean (red) and standard deviation (blue) of |⟨eitH⟩|2 obtained from sampling 10,000 independent 100×100 GUE matrices H with our approximation (132.13). In the second plot we similarly compare the empirical mean (red) and variance (blue), with respect to s, obtained from sampling 500 independent scalar random variables s (from the uniform distribution on S1) and 500 independent 100×100 GUE matrix pairs H1,H2, with the prediction (152.15). We also test the precision of the latter GUE-pair sampling by finding the empirical standard deviation (with respect to H1,H2) of the empirical mean of the monoparametric SFF (orange). We observe that for both ensembles our resolvent approximation seems valid for all t0 (possibly N-dependent) we have2.8 EH|⟨eitH⟩|2=Ewig(t)(1+o(1))fort≪N5/11,VarH|⟨eitH⟩|2=Swig(t)2(1+o(1))fort≪N5/17, and we have the asymptotics2.9 Ewig(t):=e(t)2+v-,κee(t)N2≈J1(2t)2t2,1≪t≪N2πtN2,N≪t≪N,Swig(t):=(v+,κee(t)2+v-,κee(t)2N4+2e(t)2v+,κee(t)+v-,κee(t)N2)1/2≈2J1(2t)πtN,1≪t≪N2πtN2,N≪t≪N, where we set e:=(1,0)∈S1. This result shows that Swig(t)≪Ewig(t) in the slope regime, t≪N, and Swig(t)≈Ewig(t) in the ramp regime, N≪t≪N (see the first plot in Fig. 2). In particular, in the ramp regime the SFF is a non-negative random variable whose fluctuations are of the same size as its expectation. Thus the SFF is not self-averaging in the ramp regime, while it is self-averaging in the slope regime but only owing to the dominance of the function e(t) representing the edge effect. If one discounts the edge effect, i.e. artificially removes e(t), then Swig(t)≈Ewig(t) would hold for all 1≪t≪N, demonstrating the universal behavior of SFF in the entire slope-dip-ramp regime. Theorem 2.8 (SFF for the monoparametric ensemble (2.2)). For t>0 (possibly N-dependent) we have10 EHEs|⟨eitHs⟩|2=Ewig(t)(1+o(1))fort≪N3/7EHVars|⟨eitHs⟩|2=(Swig(t)2-Sres(t)2)(1+o(1))fort≪N5/17VarHEs|⟨eitHs⟩|2=Sres(t)2(1+o(1))fort≪N1/4 where the function2.10 Sres(t):=EsEr(v+,κsr(t)2+v-,κsr(t)2N4+2e(t)2v+,κsr(t)+v-,κsr(t)N2) satisfies2.11 Sres(t)∼ψ(t)Nt5/4,1≪t≪Nt3/4N2,N≪t≪N, where ψ(t)∼1 is a positive function with some oscillation. In particular, this result immediately shows the following concentration effect: Corollary 2.9 For 1≪t≪N1/4 it holds that2.12 VarHEs|⟨eitHs⟩|2≲1tVarH|⟨eitH⟩|2, i.e. averaging in s reduces the size of the fluctuation of the SFF by a factor of t-1/4. Note that13 Sres(t)≲t-1/4Swig(t) both in the slope and ramp regimes showing that not only the expectation but also the variance of the SFF for the monoparametric ensemble coincide with those for the Wigner ensemble to leading order, hence they follow the universal pattern (red and blue curves in the second plot in Fig. 2). However, the dependence of Es[SFF(t)] on the fixed Wigner matrix pair (H1,H2) is still present, albeit to a lower order, expressed by the residual standard deviation Sres(t) whose relative size decreases as t-1/4 as t increases (orange curves in Fig. 2). It is quite remarkable that a single random mode s generates almost the entire randomness in the ensemble that is responsible for the universality of SFF. A similar phenomenon was manifested in the Wigner–Dyson universality proven in [15]. Remark 2.10 Based upon extensive numerics (see Fig. 2) we believe that (2.13), (2.15) and (2.18) hold up to any t≪N, i.e. in the entire slope-dip-ramp regime and not only up to some fractional power of N as stated and proved rigorously. The proof for the entire regime t≪N is out of reach with the current technology based upon the multi-resolvent local law Lemma 3.3 whose error term does not trace the expected improvement due to different spectral parameters z1≠z2. We expect that the entire ramp regime t≪N should be accessible by resolvent techniques if a sharp version of Lemma 3.3, tracing the gain from z1≠z2, was available. Remark 2.11 We stated Theorems 2.7 and 2.8 only for the first two moments but the CLT from Theorem 2.5 allows us to compute arbitrary moments E|⟨eitH⟩|2m for the Wigner case and Es|⟨eitHs⟩|2m for the monoparametric case (together with their concentration in the (H1,H2)-space), albeit with worsening error estimates. This would lead to rigorous results of the type (2.13) and (2.15) but for a shorter time scale t≪Nc(m) with some c(m)>0. However, in the spirit of Remark 2.10, we believe that ⟨eitHs⟩ can be approximated for any t≪N, to leading order, by a family of complex Gaussians ξ(t,s) of mean and variance2.13 Eξ(t,s)=e(t),E(ξ(t,s)-e(t))(ξ(t′,s′)-e(t′))=1N2vss′(eit·,eit′·) with vsr from (2.7). Note that (202.20) also specifies the covariance of ξ(t,s) and ξ(t′,s′)¯=ξ(-t′,s′) for different times. The next lemma, to be proven in Sect. 3.2.4, provides explicit asymptotic formulas for v±ss(t), in particular they imply the asymptotics in (2.14) together with e(t)∼t-3/2 (up to some oscillation due to the Bessel function) in the large t regime. Lemma 2.12 For s=r the functions v±ss(t) appearing in (2.12) can be expressed as14 v-ss(t)=t2[J0(2t)2+2J1(2t)2-J0(2t)J2(2t)]=2tπ-1+2sin(4t)16πt+O(t-2)v+ss(t)=-tJ0(2t)J1(2t)=cos(4t)2π-2+sin(4t)16πt+O(t-2). The relation in (2.17) requires a stationary phase calculation that will be done separately in Sect. 5. Implications for sampling Determining the standard deviation of |⟨eitH⟩|2 is important for numerical testing of (2.13). By taking the empirical average EHn of n independent Wigner matrices we may approximate the true expectation EH|⟨eitH⟩|2 at a speed2.14 EHn|⟨eitH⟩|2=EH|⟨eitH⟩|2+Ω(n-1/2StdH|⟨eitH⟩|2)=Ewig(t)+Ω(n-1/2Swig(t)), c.f. the top of Fig. 3. Here Ω(⋯) indicates an oscillatory error term of the given size. In the ramp regime the fluctuation of EHn|⟨eitH⟩|2 thus scales like t/(nN2) using (2.14). In particular, this fluctuation vanishes as the sample size n goes to infinite, hence the statistics via sampling to test (2.13) is consistent.Fig. 3 In the first plot we show the empirical mean of |⟨eitH⟩|2 for k independent GUE matrices H. As expected the standard deviation of the sample average fluctuates within a strip of width n-1/2StdH|⟨eitH⟩|2, in particular the sample average exactly reproduces the mean if n→∞. In the second plot we show the empirical mean of |⟨eitHs⟩|2 for k independently sampled scalar random variables s for a fixed GUE matrix pair H1,H2. We observe that while the sample mean approximates the true mean Es increasingly well as n→∞, the latter is still dependent on the chosen realisation of H1,H2. Thus the empirical mean fluctuates in a strip of width max(n-1/2Swig(t),Sres(t)) around the doubly averaged EHEs|⟨eitHs⟩|2 In contrast, for the monoparametric ensemble, by taking the empirical average of n copies of s we naturally have15 Esn|⟨eitHs⟩|2=Es|⟨eitHs⟩|2+Ω(k-1/2Swig(t)). Replacing the first term by its expectation plus its fluctuation in the H-probability space, we also get2.15 Esn|⟨eitHs⟩|2=EHEs|⟨eitHs⟩|2+Ω(max(n-1/2Swig(t),Sres(t))), where the error term contains both standard deviations and satisfies16 max(n-1/2Swig(t),Sres(t))∼1Ntmax{1n,1t1/4},1≪t≪NtN2max{1n,1t1/4}N≪t≪N, due to (2.15) and (2.17). In particular, both in the slope and in the ramp regimes the size of the fluctuation of Esn|⟨eitHs⟩|2 does not vanish even as the number of samples goes to infinity, n→∞, hence the statistics is not consistent, c.f. the bottom of Fig. 3. However, this lack of consistency, expressed by Sres(t) is still negligible compared with the leading first term in (2.24) by a factor t-1/4≪1 in the large t regime, see (2.19). We recall that mathematically rigorously we can prove all these facts only for t≪N1/4, i.e. well before the dip time, but the numerical tests leave no doubt on their validity in the entire regime 1≪t≪N. Extensions Beside the Wigner ensemble, we formulated our main results on SFF for the normalized monoparametric model in Theorem 2.8. We chose this model for definiteness, but our approach applies to the multi-parametric as well as to the un-normalised models introduced in Remark 2.4. Here we explain the modified results for these natural generalisations. First, for the multi-parametric normalised model, Hs=s1H1+⋯+skHk with k-1 effective parameters s∈Sk-1, Theorem 2.8 holds true verbatim modulo different sizes for the residual standard deviation Sres(t). In fact, we have2.16 Sres(t)≲t-12+14(3-k)+Swig(t), see (5.4) later, hence Sres(t) becomes less relevant compared with Swig(t) for larger k>2, see (2.19). Consequently, the upper bounds on the times of proven validity in (152.15) slightly improve but they still remain below the dip time and we omit the precise formulas. We note that the t-power in (2.26) is not optimal for k≥3. A refined stationary phase estimate could be used to improve the estimate but we refrain from doing so since our primary interest is the mono-parametric model with few degrees of freedom. Second, for the un-normalised model Hs=s1H1+s2H2 with two effective parameters s∈R2, Theorem 2.8 also holds true modulo some minor changes. More precisely, (152.15) becomes17 EHEs|⟨eitHs⟩|2=EsEwig(‖s‖2t)(1+o(1))fort≪N3/7EHVars|⟨eitHs⟩|2=(EsSwig(‖s‖2t)2-S~res(t)2)(1+o(1))fort≪N5/17VarHEs|⟨eitHs⟩|2=S~res(t)2(1+o(1))fort≪N1/7, with S~res obtained from replacing vsr by v~sr from Remark 2.6 in (2.12). For S~res(t) a stationary phase calculation gives the modified2.17 S~res(t)∼ψ(t)Nt7/4,1≪t≪Nt1/4N2,N≪t≪N, assuming that s has an absolutely continuous distribution with a differentiable, compactly supported density ρ on R2 with ρ(0)=0. We will not prove the relation in these formulas in this paper, we only show how to obtain the necessary upper bound on them at the end of Sect. 5. Note that now18 S~res(t)≲t-3/4EsSwig(‖s‖2t), i.e. the fluctuation due to the residual randomness of (H1,H2) after taking the expectation in s remains negligible, in fact it is reduced compared with the normalised case (2.19). As a consequence t1/4 in (2.25) is replaced by t3/4. Analogous results hold for the most general multi-parametric un-normalised model as well as to the mono-parametric inhomogeneous model Hs=H1+sH2, s∈R. We omit their precise formulation, the key point is that the analogue of (2.27) hold in all cases with a residual standard deviation S~res(t) being smaller than the leading term Swig(t) by a polynomial factor in t (e.g. by t-1/2 for Hs=H1+sH2). This guarantees that the universality of SFF holds for all these models. Table 1 summarizes the decay exponents of our main parametric models.Table 1 For our three main parametric models the following table lists the size of the residual fluctuation compared to the fluctuation of the Wigner-SFF Quenched parametric model Randomness s1H1+s2H2 (s1,s2)∈S1 Sres(t)≲t-1/4Swig(t) H1+sH2 s∈R Sres(t)≲t-1/2Swig(t) s1H1+s2H2 (s1,s2)∈R2 Sres(t)≲t-3/4Swig(t) Outline The rest of the paper is organised as follows. In Sect. 3 we outline the resolvent method and explain how via the Helffer–Sjöstrand representation a resolvent CLT implies the CLT for the linear statistics ∑f(λi) of arbitrary test functions f from which our main results Theorems 2.5–2.8 follow. In Sect. 4 we present the proof of the resolvent CLT, while in Sect. 5 we conclude the proof of the asympotics (172.17) via a stationary phase argument. Resolvent Method Let H be a Wigner matrix5 and G(z):=(H-z)-1 its resolvent with a spectral parameter z∈C\R. Define msc(z), the Stieltjes transform of the semicircle law:2.18 m(z)=msc(z):=∫Rρsc(x)x-zdx,ρsc(x):=(4-x2)+2π. The local law for a single resolvent states that the diagonal matrix m(z)·I well approximates the random resolvent G(z) in the following sense (see e.g. [5, 20, 35]):2.19 |⟨(G(z)-m(z))A⟩|≲Nξ‖A‖Nη,⟨x,(G(z)-m(z))y⟩≲Nξ‖x‖‖y‖Nη with η=|ℑz|, for any fixed deterministic matrix A and deterministic vectors x,y. The first bound is called averaged local law, while the second one is the isotropic local law. The bounds (3.2) are understood in very high probability for any fixed ξ>0. The Helffer–Sjöstrand formula20 ⟨f(H)⟩=2π∫C∂z¯fC(z)⟨G(z)⟩d2z, with z=x+iη and d2z:=dηdx, expresses the linear statistics of arbitrary functions as an integral of the resolvent G(z) and the almost-analytic extension2.20 fC(z)=fC(x+iη):=[f(x)+iη∂xf(x)]χ(τη), of f. Here the free parameter τ∈R is chosen such that N-1≪τ-1≲1, and χ a smooth cut-off equal to 1 on [-5,5] and equal to 0 on [-10,10]c. The same τ was used to define the weighted H2-norm (2.4) and eventually we will optimize its value, a procedure that improves the standard error terms in the CLT. By (3.2) it follows that21 ⟨f(H)⟩=2π∫C∂z¯fC(z)m(z)d2z+O(∗)Nξ‖f‖H2N=∫-22ρsc(x)f(x)dx+O(∗)Nξ‖f‖H2N. In order to compute the fluctuation in (3.5) via (3.3) we need to understand the correlation between ⟨G(z)⟩,⟨G(z′)⟩ for two different spectral parameters z,z′ which turns out to be given by2.21 Cov(⟨G(z)⟩,⟨G(z′)⟩)≈1N2⟨G(z)2⟩⟨G(z′)2⟩⟨G(z)G(z′)⟩(1+⟨G(z)G(z′)⟩)⟨G(z)⟩⟨G(z′)⟩, modulo some additional contribution from non-Gaussian fourth cumulant, see (83.8) for the final statement. While G(z)≈m(z), in general it is not true that G(z)G(z′)≈m(z)m(z′) since (3.2) allows only deterministic test matrices multiplying G. Nevertheless G(z)G(z′) is still approximable by a deterministic object:2.22 G(z)G(z′)≈m(z)m(z′)1-m(z)m(z′). Statements of the form (3.7) with an appropriate error term are called multi-resolvent local laws. We will apply this theory to the product of the resolvents Gs of Hs=s1H1+s2H2 for two different parameters s, see the corresponding local law on ⟨GsGr⟩ in (3.11) later. Even though H1 and H2 as well as s and r are independent, the common (H1,H2) ingredients in Hs and Hr introduce a nontrivial correlation between these matrices. We therefore need to extend CLT for resolvents via multi-resolvent local laws to this parametric situation. Resolvent CLT The main technical result of the present paper is the following Central Limit Theorem for product of resolvents of the random matrix Hs:=s1H1+s2H2 with s=(s1,s2)∈S1. Proposition 3.1 Fix ϵ>0, p∈N, s1,⋯,sp∈S1, z1,⋯,zp∈C\R, and define Gi:=(Hsi-zi)-1. Then for any arbitrary small ξ>0 and η∗≥N-1+ϵ it holds2.23 EH∏i∈[p]⟨Gi-EHGi⟩=1Np∑P∈Pair([p])∏(i,j)∈PVij+ONξΨp1L1/2+1Nη∗2+1N2η∗4,Ψp:=∏i∈[p]1N|ηi|. Here ηi:=ℑzi, η∗:=mini|ηi|, L:=mini(Nηiρi), and2.24 Vij:=-2β∂zi∂zjlog(1-⟨si,sj⟩mimj)-⟨si⊙si,sj⊙sj⟩κ4(mi2)′(mj2)′, where mi:=msc(zi), and κ4:=N2E|h12|4-1-2/β. Additionally, for the expectation we have2.25 EH⟨Gi⟩=mi+κ4N‖s‖44mi′mi3+1(β=1)1Nmimi′1-mi2+ONξρi(Nηi)3/2, with ρi:=π-1|ℑmi|. Remark 3.2 For Wigner matrices, i.e. for s1=⋯=sp=(1,0), the error term in (3.8) is given by ΨL-1/2, as a consequence of the fact that the error terms in the first and second line of (3.11) are replaced by (Nη1η2)-1 and (Nη1η22)-1, respectively (see e.g. [16, Remark 3.5]). We point out that similar resolvent CLT have often been used as a basic input to prove CLT for linear eigenvalue statistics of both Hermitian and non-Hermitian matrices down to optimal mesoscopic scales (see e.g. [10, 11, 14, 28–30, 36–38]). The main novelty here is to extend the resolvent CLT to the monoparametric ensemble. Along the proof of Proposition 3.1 we establish the following multi-resolvent local laws. Lemma 3.3 For Gi=Gsi(zi) we have the two- and three-resolvent local laws2.26 |⟨G1G2⟩-m1m21-⟨s1,s2⟩m1m2|≲NξN|η1η2|3/2|⟨G1G22⟩-m1m2′(1-⟨s1,s2⟩m1m2)2|≲NξN|η1||η2|η∗2+1N2|η1η2|3, where mi=msc(zi), with very high probability for any fixed ξ,ϵ>0 and |ℑzi|≥N-1+ϵ. The proofs of Proposition 3.1 and Lemma 3.3 will be presented in Sect. 4. In these proofs we will often use the standard cumulant expansion (see [8, 30, 34] in the random matrix context):2.27 EHhabf(H)=1NEH∂baf(H)+∑k=2R∑q+q′=kκabq+1,q′N(k+1)/2EH∂abq∂baq′f(H)+ΩR. Here ∂ab denotes the directional derivative ∂hab, the first term in the rhs. represents the second order (Gaussian) contribution, while the sum in (3.12) represents the non-Gaussian contribution with κabp,q denoting the joint cumulant of p copies of N1/2hab and q copies of N1/2hab¯. The cumulant expansion is typically truncated at a high (N-independent) level R with an error term ΩR that is negligible. To see this, note that in our applications f will be a product of resolvents at spectral parameters zi with η∗=min|ℑzi|≫1/N hence derivatives of f remain bounded with very high probability by the isotropic local law (3.2) thus the tail of the series (3.12) decays as N-(k+1)/2. Proof of Theorem 2.5 The proof of Theorem 2.5 is divided into three steps: (i) computation of the expectation, (ii) computation of the variance, (iii) proof of Wick Theorem. The expectation is computed in Sect. 3.2.1, while the Wick Theorem and the explicit computation of the variance are proven in Sect. 3.2.2. Expectation Using the bound2.28 |∂z¯fC|≲η|f′′|+τ|χ′|[|f|+iη|f′|], and |⟨Gs-m⟩|≲Nξ(Nη)-1 by (3.2), with m=msc, we conclude that2.29 EH⟨f(Hs)⟩=∫R∫|η|≥η0∂z¯fC(z)EH⟨Gs(z)⟩dηdx+ONξη0‖f‖H2N+Nξη02‖f‖H2, for any N-1≪η0≪τ-1. Note that we chose η0≫N-1 in order to use Proposition 3.1. Plugging (3.10) into (3.14), and using (3.13) to estimate the error term, we get that3.1 EH⟨f(Hs)⟩=∫R∫|η|≥η0∂z¯fC(z)[m+κ4N‖s‖44m′m3+1(β=1)1Nmm′1-m2]dηdx+ONξη0‖f‖H2N+Nξη02‖f‖H2+Nξ‖f‖H2N3/2τ1/2+Nξτ3/2‖f‖∞N3/2+Nξτ1/2‖f‖H1N3/2=∫R∫|η|≥η0∂z¯fC(z)[m+κ4N‖s‖44m′m3+1(β=1)1Nmm′1-m2]dηdx+ONξ‖f‖τN3/2τ1/2, where to go to the last line we chose η0∼N-1+ϵ, for some very small ϵ>0, and we used the norm ‖f‖τ defined in (2.4). Adding back the regime |η|<η0 at the price of a negligible error smaller than the one in (3.15), by explicit computations (exactly as in [13, Section D.1]) in the leading term of (3.15), we conclude3.2 EH⟨f(Hs)⟩=∫-22ρsc(x)f(x)dx+κ42N‖s‖44∫-22x4-4x2+2π4-x2f(x)dx+1(β=1)f(2)+f(-2)4N-12πN∫-22f(x)4-x2dx+ONξ‖f‖τN3/2τ1/2. Second moment and Wick theorem Define3.3 LN(f,s):=N[⟨f(Hs)⟩-EH⟨f(Hs)⟩], then in this section, using Proposition 3.1, we compute the leading order term of EHLN(f1,s1)LN(f2,s2). More precisely, by (3.8) for p=2, and using (3.13) to estimate the error term, it follows that3.4 EHLN(f1,s1)LN(f2,s2)=∬R∬|η1|,|η2|≥η0∂z1¯fC(z1)∂z2¯fC(z2)V12+O(Nξη0(‖f1‖H2‖f2‖∞+‖f2‖H2‖f1‖∞)+Nξ‖f1‖τ‖f2‖τN1/2τ3/2+‖f1‖H2‖f2‖H2N1-ξη0τ1+1Nη02+(‖f1‖H2(τ2‖f2‖∞+τ‖f2‖H1)+‖f2‖H2(τ2‖f1‖∞+τ‖f1‖H1))N1-ξη0τ1+1Nη02+(τ2‖f1‖∞+τ‖f1‖H1)(τ2‖f2‖∞+τ‖f2‖H1)N1+τ2N1-2ϵ)=∬R∬|η1|,|η2|≥N-ϵτ-1∂z1¯fC(z1)∂z2¯fC(z2)V12+ONξ‖f1‖τ‖f2‖τNϵN+N-ϵτ31+τ2N1-2ϵ, where to go to the last line we chose η0∼N-ϵτ-1, for any ϵ>0, and V12 is defined in (3.9). From (3.18), adding back the regimes |ηi|0. We remark that by 1-m2⟨·⟩ in the lhs. of (4.3) we denote the operator acting on matrices R∈CN×N as (1-m2⟨·⟩)[R]=R-m2⟨R⟩. We then start computing:11 EH⟨G-m⟩=-m′mEH⟨s1H1G_+s2H2G_⟩+ONξN2η2ρ, for any small ξ>0, where we used that |1-m2|≳ρ, that m′=m2/(1-m2), and that |⟨G-m⟩|≲Nξ(Nη)-1 by (3.2). Then using cumulant expansion (see (3.12), ignoring the truncation error) we claim (and prove below) that3.11 EH⟨s11H1G_+s21H2G_⟩=EH1N∑k≥2∑ab∑α∈{ab,ba}kκ(1)(ab,α)k!s1∂α(1)+κ(2)(ab,α)k!s2∂α(2)Gba=κ4N‖s‖44m4+ONξρ3/2N3/2η1/2+Nξρ3/2N2η3/2, where κ(i)(ab,α) denotes the joint cumulant of the random variables habi, hα1i,⋯,hαki, and ∂α(i):=∂α1(i)⋯∂αk(i), with i=1,2, where ∂αj(i) denotes the directional derivative in the direction hαji. Here hαji are the entries of Hi. Combining (4.5) with (4.4) we obtain exactly the expansion in (3.10) (recall that here we only present the proof in the complex case, the real case being completely analogous). Proof of the second equality in (4.5) First of all we recall that by (2.1) it follows the bound |κ(i)(ab,α)|≲N-(k+1)/2, with i=1,2. We start with k=2. In this case we can neglect the summation when a=b since it gives a contribution N-3/2. Hence we can assume that a≠b. In this case we have the bounds3.12 N-5/2∑a≠bGab3≲Nξρ3/2N2η3/2,N-5/2∑a≠bGaaGbbGab≲NξN3/2+Nξρ3/2N2η3/2, with very high probability. The first bound in (4.6) follows from the isotropic law in (3.2). The second bound in (4.6) follows by writing G=m+(G-m) and using the isotropic resummation3.13 ∑ab(G-m)aaGab=∑a⟨ea,G1⟩, with ea∈RN the unit vector in the a-direction and 1:=(1,⋯,1)∈RN. For k=3 whenever there are at least two off-diagonal G’s we get a bound N-2η-1ρ. The only way to get only diagonal G’s is that α is one of (ab, ba, ba), (ba, ab, ba), (ba, ba, ab); in this case κ(i)(ab,α)=κ4/N2, with κ4:=κ(i)(ab,ba,ab,ba). For these terms we have (see [13, Lemma 4.2] for the analogous proof for Wigner matrices)3.14 ∂α(i)Gba=-2si3Gaa2Gbb2+ONξρN2η, with very high probability, where the error comes from terms with at least two off-diagonal G’s. Hence we finally conclude that the terms k=3 give a contribution:3.15 -2κ433!‖s‖441N3∑abGaa2Gbb2=κ4N‖s‖44m4+ONξρ3/2N3/2η1/2+NξρN2η. All the terms with k≥4 can be estimated trivially using that |Gab|≲1 with very high probability by (3.2). □ Computation of the variance For the second moment, using (4.3), we compute3.16 EH⟨G1-EHG1⟩⟨G2-EHG2⟩=-EHm1′m1⟨s11H1G1_+s21H2G1_⟩+κ4N‖s1‖44m1′m13⟨G2-EHG2⟩+ONξΨ2L1/2 where si=(s1i,s2i)∈S1 and we used (3.10) to approximate ⟨Gi-EHGi⟩ with ⟨Gi-mi⟩. We made this replacement to use the equation for G-m from (4.3). Then performing cumulant expansion we compute:3.17 -EHm1′m1⟨s11H1G1_+s21H2G1_⟩+κ4N‖s1‖44m1′m13⟨G2-EHG2⟩=⟨s1,s2⟩m1′EH⟨G1G22⟩m1N2-κ4N‖s1‖44m1′m13EH⟨G2-EHG2⟩-m1′m1∑k≥2∑ab∑α∈{ab,ba}kκ(1)(ab,α)k!Ns11∂α(1)+s21κ(2)(ab,α)k!N∂α(2)EH[(G1)ba⟨G2-EHG2⟩]. Using the local law (113.11) we conclude that3.18 m1′m1⟨s1,s2⟩⟨G1G22⟩N2=⟨s1,s2⟩m1′m2′(1-⟨s1,s2⟩m1m2)2N2+ONξN3η1η2η∗2+NξN4|η1η2|3=-1N2∂z1∂z2log(1-⟨s1,s2⟩m1m2)+ONξN3η1η2η∗2+NξN4|η1η2|3, with very high probability. We are now left with the third line of (4.11). The α-derivative in (4.11) may hit either (G1)ba or ⟨G2-E2G2⟩. Define3.19 Φk:=m1′m1∑ab∑α∈{ab,ba}kκ(1)(ab,α)k!Ns11∂α(1)+s21κ(2)(ab,α)k!N∂α(2)EH[(G1)ba⟨G2-EHG2⟩]=∑ab∑αs11κ(1)(ab,α)k!NEHm1′m1∂α1(1)(G1)bak1!∂α2(1)⟨G2-EHG2⟩(k-k1)!+∑ab∑αs21κ(2)(ab,α)k!NEHm1′m1∂α1(2)(G1)bak1!∂α2(2)⟨G2-EHG2⟩(k-k1)!, where k1 denotes the number of derivatives that hit (G1)ba. The summation ∑α indicates the summation over tuples αiki, with i=1,2 and k2:=k-k1. We now claim that3.20 Φk=-1(k=3)(κ4⟨s1⊙s1,s2⊙s2⟩2N2(m12)′(m22)′+κ4N‖s1‖44m1′m13)+ONξΨ2L1/2. Similarly to the proof of [13, Eq. (113)] we readily conclude that the terms in Φk in (4.13) with k=2, or k1 odd and k≥4, or k≥3 and k1 even are bounded by NξΨ2L-1/2. For k=3 and k1=3, analogously to (4.8)–(4.9) we obtain a contribution of3.21 -κ4N‖s1‖44m1′m13+ONξN|η1|L1/2 to (4.14). For k=3 and k1=1 we start computing the action of the α1-derivative on (G1)ba:3.22 ∑α1∂α1(i)(G1)ba=-si1(G1)ba2-si1(G1)aa(G1)bb=-si1m12(1+δab)+ONξρ1N|η1|, with very high probability. Additionally, we have that (see [13, Lemma 4.2] for the analogous proof for Wigner matrices)3.23 ∂ab,ba(i)⟨G2-EHG2⟩=2m2m2′N(si2)2+ONξρ21/2(N|η2|)3/2, with very high probability. We thus conclude that the (k,k1)=(3,1) contribution to (4.14) is3.24 -κ4⟨s1⊙s1,s2⊙s2⟩2N2(m12)′(m22)′+ONξΨ2L1/2, where we used that only the terms with κ4=κ(i)(ab,ba,ab,ba) contribute. This concludes the proof of (3.8) for p=2. Asymptotic Wick Theorem The proof of the Wick Theorem for resolvent is completely analogous to the one for Wigner matrices in [13, Section 4]. The only differences are that along the proof we have to carefully keep track of the si, as we did in Sect. 4.2, since in the Wigner case s1=⋯=sp=(1,0), and that we have to use the three G’s local law in (3.11) with a weaker error term instead of the one in [13, Eq. (45)] to compute the leading order deterministic term (see (4.21)–(4.22) below). Define4.1 YS:=∏i∈S⟨Gi-EHGi⟩, with S⊂N. Similarly to Sect. 4.2 we start computing4.2 EHY[p]=∑i∈[2,p]m1′m1⟨s1,si⟩N2EH⟨G1Gi2⟩Y[p]\{1,i}-κ4N‖s1‖44m1′m13EHY[2,p]-∑k≥2∑ab∑α∈{ab,ba}kκ(1)(ab,α)k!Ns11∂α(1)+s21κ(2)(ab,α)k!N∂α(2)EHm1′m1(G1)baY[2,p]+ONξΨpL1/2. Then proceeding analogously to (4.13)–(4.18) (see also [13, Eqs. (110)-(114)] for the Wigner case) we conclude that4.3 EHY[p]=∑i∈[2,p]m1′m1⟨s1,si⟩N2EH⟨G1Gi2⟩Y[p]\{1,i}-∑i∈[2,p]κ4⟨s1⊙s1,si⊙si⟩2N2(m12)′(mi2)′EHY[1,p]\{1,i}+ONξΨpL1/2. In order to compute the leading deterministic term of ⟨G1Gi2⟩ we use the local law (3.11) and get4.4 EHY[p]=1N2∑i∈[2,p]V1,iEHY[p]\{1,i}+ONξΨp1L1/2+1Nη∗2+1N2η∗4. Finally, proceeding iteratively we conclude (3.8). Multi resolvents local laws The goal of this section is to prove the local laws in (3.11). Starting from (4.3) we get4.5 (1-⟨s1,s2⟩m1m2⟨·⟩)G1G2=m1m2+m1⟨G2-m2⟩-m1(s11H1G1G2_+s21H2G1G2_)+m1⟨s1,s2⟩⟨G1G2⟩(G2-m2)+m1⟨G1-m1⟩G1G2. We estimate |⟨G1G2⟩|≲Nξ(η∗)-1 with very high probability, where η∗:=η1∨η2, using |⟨G1G2⟩|≤⟨|G2|⟩/η1 (in case η∗=η1) and the rigidity of eigenvalues to estimate ⟨|G2|⟩≤Nξ. Then by the single resolvent local law |⟨Gi-mi⟩|≲Nξ(Nηi)-1 from (3.2) we obtain that4.6 (1-⟨s1,s2⟩m1m2)⟨G1G2⟩=m1m2-m1(s11⟨H1G1G2_⟩+s21⟨H2G1G2_⟩)+ONξN|η1||η2|, with very high probability. Finally, using that4.7 |⟨HiG1G2_⟩|≲NξN|η1η2|η∗,i∈[2] with very high probability from an analogous proof to [15, Eq. (5.8)] (see also [12, Eq. (5.10c)]), and that4.8 |1-⟨s1,s2⟩m1m2|≳η∗. we conclude the first local law in (3.11). For the second local law in (3.11) we start writing the equation for G1G22:4.9 G1G22=m1m2′+m1(G22-m2′)-m1(s11H1G1G22_+s21H2G1G22_)+m1⟨s1,s2⟩(⟨G1G2⟩G22+⟨G1G22⟩G2)+m1⟨G1-m1⟩G1G22. Then, using the usual single G local law and the two G’s local law from (3.11), we conclude that4.10 (1-⟨s1,s2⟩m1m2)⟨G1G22⟩=m1m2′+⟨s1,s2⟩m12m2m2′1-⟨s1,s2⟩m1m2-m1(s11H1G1G22_+s21H2G1G22_)+ONξN|η1||η2|η∗. Then, using that4.11 |⟨HiG1G22_⟩|≲NξN|η1η2|η∗2,i∈[2], with very high probability, and (4.26) we conclude (3.11). The proof of (4.29) follows analogously to the one of (4.25). Stationary Phase Calculations The proof of (2.17) is a tedious stationary phase calculation since v±sr(t), the leading part of v±,κsr(t) (see (2.12)), are given in terms of oscillatory integrals for t≫1 being the large parameter. Unlike in the s=r case, no explicit formula similar to (2.21) is available. The main complication is that Vsr(x,y) defined in (2.7) has logarithmic singularities, integrated against a fast oscillatory term from f′g′, so standard stationary phase formulas cannot directly be applied. Nevertheless, a certain number of integration by parts can still be performed before the derivative of the integrand stops being integrable and the leading term can be computed. We will first give a proof of4.12 EsErv-sr(t)∼t then we explain how to modify this argument to obtain4.13 EsErv-sr(t)2∼t3/2, in both cases with a definite large t asymptotics with computable explicit constants. The proof reveals that the corresponding results for EsErv+sr(t) and EsErv+sr(t)2 guarantee only an upper bound with the same behavior4.14 EsErv+sr(t)≲t,EsErv+sr(t)2≲t3/2 depending on the distribution of s on S1, the matching lower bound may not necessarily hold. However, for our main conclusions like (2.19) only an upper bound on Sres(t) is important. All these exponents are valid for the k=2 case, i.e. for Hs=s1H1+s2H2. For the general multivariate model, k≥3, exactly the same proof gives the upper bounds4.15 EsErv±sr(t)≲min{1,t3-k2},EsErv+sr(t)2≲min{1,t5-k2}. The k-dependence of the exponent can directly be related to the tail behavior (5.5) and (5.8) below, so for simplicity we will carry out our main analysis only for k=2. In fact, a more careful analysis yields somewhat better bounds than (5.4), but we will not pursue this improvement here. We introduce a new random variableU:=⟨s,r⟩ then clearly |U|≤1 and since r,s∈Sk have a distribution with an L2 density, it is easy to see that the density ρ∗ of U is bounded by4.16 ρ∗(U)≲(1-U2)k-32. The fact that the main contribution to the lhs. of (5.4) comes from the regime U≈1 is a consequence of the singularity of the logarithm in (2.7) in this regime (see computations below). Indeed, U=cosα where α is the angle between r, s and near U≈±1 we have 1±U≈12α2(1+O(α2)). For example, for k=2 we have4.17 P(1-U=ϵ+dϵ)=dϵϵ(∫S1ρ2(s)ds)(1+O(ϵ)) 4.18 P(1+U=ϵ+dϵ)=dϵϵ(∫S1ρ(s)ρ(s+π)ds)(1+O(ϵ)) in the ϵ≪1 regime. In particular, the bound in (5.5) is actually an asymptotics in the most critical U≈1 regime, while the regime U≈-1 it may happen that the density ρ∗ is much smaller than (5.5) predicts. For symmetric distribution, ρ(s)=ρ(s+π), the two asymptotics are the same. Similar relations hold for k≥3, in which case we have4.19 P(1±U=ϵ+dϵ)≲ϵk-32dϵ with an explicit asymptotics for U≈1. So we will study4.20 R±(t)=t2ℜ∫dUρ∗(U)∬-22dxdyeit(x±y)[log|1-Um(x)m(y¯)|-log|1-Um(x)m(y)|]. Since |m|≤1, as long as |U|≤1-δ for any small fixed δ>0, the arguments of the logarithms are separated away from zero and they allow to perform arbitrary number of integration by parts, each gaining a factor of 1/t. There is a square root singularity of m(x) and m(y) at the spectral edges 2,-2 which still allows one to perform one integration by parts in each variable since m′ is still integrable. Therefore the contribution of the regime |U|≤1-δ to (5.9) is of order t2(1/t)2=O(1), hence negligible compared with the target (5.1). In the sequel we thus focus on the important U≈±1 regimes, in particular every ∫dU integral is understood to be restricted to |U|≥1-δ. Note that m(y)¯=-m(-y), so if U has a symmetric distribution (for example if s∈S1 has a symmetric distribution), then by symmetry we haveR-(t)=-R+(t). For definiteness, we focus on R-(t), the analysis of R+ is analogous. From the explicit form m(x)=12(-x+i4-x2) a simple exercise shows that4.21 |1-Um(x)m(y)¯|2≳(1-U)2+(x-y)2,|1-Um(x)m(y)|2≳(1+U)2+(x+y)2. This shows that the critical regime is U≈1 and x≈y for the first integrand in (5.9) and U≈-1, x≈-y for the second. Again, for definiteness, we focus on the first regime, i.e. on the first log-integrand in (5.9) and establish the following relations for large t and k=2: Lemma 5.1 In the k=2 case we have4.22 t2∫dUρ∗(U)ℜ∬-22eit(x-y)log|1-Um(x)m(y¯)|2dxdy∼t and4.23 t4∫dUρ∗(U)[ℜ∬-22eit(x-y)log|1-Um(x)m(y¯)|2dxdy]2∼t3/2, for t≥1. For t≫1 an analogous asymptotic statement holds with explicitly computable positive constants that depend on the distribution of s. Proof of Lemma 5.1 Introduce the variablesa:=x+y2,b:=x-y2,i.e.x=a+b,y=a-b. Since |x|,|y|≤2 we have4.24 |a|≤2,|b|≤min{|2-a|,|2+a|}. In terms of these variables, we have4.25 |1-Um(x)m(y)¯|2=(1-U+2Ub2b2+d2)2+4U2b2d2(b2+d2)2,d:=12[4-(a+b)2+4-(a-b)2]. Here we also used the identity1-m(x)m(y)¯=2b2b+m(x)-m(y)¯=2bb+id following from the equation -m(x)-1=x+m(x) and similarly for m(y). In the regime (5.13) we have4.26 |b|≤12(4-a2),|b|≤4-a2. Note that by Taylor expansion around a and concavity of the function x→4-x2 in x∈[-2,2], we have4.27 0≤4-a2-d≲b2(4-a2)3/2≤|b|4-a2,aswellas124-a2≤d≤4-a2. We define the function4.28 F=F(U,a,b):=(1-U)2+4U2b24-a2 for |U|≤1, and a, b as in (5.13). We will use F to approximate4.29 M=M(U,a,b):=|1-Um(a+b)m(a-b)¯|2 in the critical regime where |U|≥1-δ and |b|≤δ for some small fixed δ>0. We clearly have5.1 M(U,a,b)≥14F(U,a,b) in the regime (5.13), where |b|≤4-a2≤2d, using (5.16). For the difference function5.2 Δ(U,a,b):=M(U,a,b)-F(U,a,b) an elementary calculation from (5.14)–(5.16) gives5.3 |Δ(U,a,b)|≲b2(4-a2)3/2F in the regime |U|≥1-δ and |b|≤δ. Furthermore, similar estimates hold for the first derivative;5.4 |ddbΔ(U,a,b)]|≲|b|F(4-a2)3/2,|ddaΔ(U,a,b)]|≲b2F(4-a2)5/2≲|b|F(4-a2)3/2, as well as for the second derivatives5.5 |d2db2Δ(U,a,b)]|≲F(4-a2)3/2,|ddaddbΔ(U,a,b)]|≲|b|F(4-a2)5/2≲F(4-a2)3/2. The proof of Lemma 5.1 consists of two parts. First we compute the integral with logF, i.e. we show that5.6 t2∫dUρ∗(U)ℜ∬-22eit(x-y)logF(U,x+y2,x-y2)dxdy∼t with an explicit positive constant factor in the asymptotic regime t≫1. Second, we show that the integrand in (5.11) can indeed be replaced with F up to a negligible error,5.7 |t2∫dUρ∗(U)∬-22eit(x-y)[log|1-Um(x)m(y¯)|2-logF(U,x+y2,x-y2)]dxdy|≲1. Part I. To prove (5.24), we use the a, b variables and the symmetry of F in a to restrict the a integration to 0≤a≤2:5.8 (5.24)=4t2ℜ∫dUρ∗(U)∫02da∫-(2-a)2-adbe2itblogF(U,a,b). Using integration by parts, we have5.9 ∫-(2-a)2-adbe2itblog[(1-U)2+4U2b24-a2]=12it[e2it(2-a)-e-2it(2-a)]log[(1-U)2+4U2(2-a)2+a]-12it4U24-a2∫-(2-a)2-adbe2itb2b(1-U)2+4U2b24-a2. In the boundary terms we can perform one more integration by parts in the a variable when plugged into (5.26). Just focusing on the first boundary term in (5.27), using |U|≤1 we have|12ite4it∫02dae-2italog[(1-U)2+4U2(2-a)2+a]|≲1t2∫02da(1-U)2+U2(2-a)≲|log(1-U)|t2. Since ρ∗(U) is a density bounded by (1-U2)-1/2 in the U≈1 regime from (5.5), the logarithmic singularity is integrable showing that the two boundary terms in (5.27), when plugged into (5.26), give at most an O(1) contribution, negligible compared with the target behavior of order t in (5.1). To compute the main (second) term in the rhs. of (5.27), we first extend the integration limits to infinity and claim that5.10 t2∫dUρ∗(U)|12it∫02da4U24-a2∫2-a∞dbe2itb2b(1-U)2+4U2b24-a2|≲t∫dUρ∗(U)∫02da2-a|∫2-a∞dbe2itb2b(1-U)2+4U2b24-a2| gives a negligible contribution to (5.26) (the lower limit is removed similarly). Indeed, we apply one more integration by parts inside the absolute value in (5.28):|∫2-a∞dbe2itb2b(1-U)2+4U2b24-a2|≲t-1∫2-a∞db(1-U)2+U2b24-a2+t-12-a(1-U)2+(2-a). Its contribution to the rhs of (5.28) is thus bounded by∫dUρ∗(U)∫02da2-a[∫2-a∞db(1-U)2+U2b24-a2+2-a(1-U)2+(2-a)]≲∫dU1-U2[∫02da2-a11-U+2-a+|log(1-U)|]≲1. Summarizing, we just proved that11 (5.24)=-2tℑ∫dUρ∗(U)∫02da4U24-a2∫-∞∞dbe2itb2b(1-U)2+4U2b24-a2+O(1)=tπ∫dUρ∗(U)∫02dae-t4-a2(1-U)/U+O(1)=c0tπ∫dU11-U∫02dae-t4-a2(1-U)/U+O(1)=c0tπ∫0∞e-vvdv∫02da(4-a2)1/4+O(1)=Γ(3/4)2Γ(5/4)c0t+O(1), where in the second line we used residue calculation, in the third line we used thatρ∗(U)=c01-U+O(1) in the regime U≈1 with some positive constant c0>0 depending on the distribution of s (see (5.6)), and finally in the fourth line we used that for large t the main contribution to the integral comes from U≈1 in order to simplify the integrand. This completes the proof of (5.24). Part II. We now prove (5.25). After changing to the a, b variables and considering only the 0≤a≤2 regime for definiteness, we perform an integration by parts in b that gives5.11 (5.25)≲t∫dUρ∗(U)|∫02dae2ita[logM(U,a,b)-logF(U,a,b)]db|+t∫dUρ∗(U)∫02da|∫-(2-a)2-ae2itb∂b[logM(U,a,b)-logF(U,a,b)]db| recalling the definition of M from (5.18). The first term in (5.30) is the boundary term, which is negligible after one more integration by parts using the ∂a derivative estimate from (5.22). In the second term we perform one more integration by parts to obtain12 (5.25)≲t∫dUρ∗(U)|∫02dae2ita∂b[logM(U,a,b)-logF(U,a,b)]db|+∫dUρ∗(U)∫02da∫-(2-a)2-a|∂b2[logM(U,a,b)-logF(U,a,b)]|db, where the first term comes from the boundary. In this term we can perform one more integration by parts in a. The corresponding boundary terms are easily seen to be order one and the main term is analogous to the first term in the rhs of (5.31) just we have the mixed ∂a∂b derivative. Recalling Δ=M-F from (5.20), we use the estimate|∂b2[logM-logF]|≲|∂b2Δ|F+|∂b2F|F2|Δ|+|∂bΔ||∂bM+∂bF|F2+(∂bF)2|Δ|F3 in the situation where M≳F>0 are positive functions (see (5.19)). Similar bound holds for the mixed derivative. Therefore, we can estimate both integrals in (5.31) as follows:5.12 (5.25)≲∫dUρ∗(U)∫02da∫-(2-a)2-a1(4-a2)3/21[(1-U)2+b24-a2]1/2db≲∫dU1-U2∫02da4-a2∫02-adu[(1-U)2+u2]1/2≲∫dU1-U2∫02|logu|+1[(1-U)2+u2]1/2du≲∫dU|log(1-U)|2(1-U)1/2≲1. Here we used the bounds (5.21), (5.22) and (5.23) and that |b|≤2-a≲4-a2 to simplify some estimates. For computing the derivatives of F we used its explicit form (5.17). This completes the proof of (5.25) and thus also the proof of (5.11) in Lemma 5.1. The proof of (5.12) is very similar. We again approximate M=|1-Um(x)m(y)¯|2 by F at the expense of negligible errors. We omit these calculations as they are very similar those for (5.11) and focus only on the main term which is (see the analogous (5.26))5.13 16t4∫dUρ∗(U)[ℜ∫02da∫-(2-a)2-adbe2itblogF(U,a,b)]2. After one integration by parts and neglecting the lower order boundary terms, we have the following analogue of (5.29):5.14 4t2∫dUρ∗(U)[ℜ∫02daU24-a2∫-∞∞dbe2itb2b(1-U)2+4U2b24-a2]2=t2π2∫dUρ∗(U)[∫02dae-t4-a2(1-U)/U]2≈c0t3/2π2∫0∞dvv(∫02dae-4-a2v)2=c0t3/2π2∬02da1da2(4-a12+4-a22)1/2∫0∞e-vvdv∼t3/2 as the leading term. This proves (5.12) and completes the proof of Lemma 5.1. □ We close this section by commenting on the proof of the upper bound in (2.28). Recall from (2.16) that the essential part of S~res(t) in the slope regime is given by EsErv~sr(t) expressed by the oscillatory integrals5.15 R±(t):=t2∬R2ρ(s)ρ(r)dsdr∬-22dxdyeit(‖s‖x±‖r‖y)A(U,x,y) withA(U,x,y):=log|1-Um(x)m(y¯)|-log|1-Um(x)m(y)|, where U=⟨s,r⟩‖s‖‖r‖ is the cosine of the angle between the vectors s,r∈R2. Assuming for the moment that ρ, the density of s, is rotationally symmetric, ρ(s)=ρ(‖s‖) with a slight abuse of notations, we have5.16 R±(t)∼t2∫-11dU1-U2∬-22dxdyA(U,x,y)∫0∞eitxσρ(σ)σdσ∫0∞e±ityσ′ρ(σ′)σ′dσ′∼t∫-11dU1-U2∬-22dxdyρ^(tx)ρ(σ)σ^(±ty)ddxA(U,x,y) performing an integration by parts in x and ignoring lower order boundary term. In the last step we also computed the Fourier transform (we used that ρ(0)=0 to extend ρ to R). The main contribution comes from the regime where A is nearly singular, and considering (5.10), we just focus on the regime U∼1 and x∼y, the singularity from the other logarithmic term is treated analogously. Similarly to the proof of (5.25) we may ignore the edge regime, and effectively we have5.17 |ddxlog|1-Um(x)m(y¯)||≲1(1-U)+|x-y|. Thus we can continue estimating the last line of (5.36)|(5.36)|≲t∫-11dU1-U2∬-22dxdy|ρ^(tx)ρ(σ)σ^(ty)|(1-U)+|x-y|≲t-1/2. Here we used the regularity of ρ, so that the last two factors essentially restrict the integration to the regime |x|,|y|≲1/t. The final inequality is obtained just by scaling. To understand S~res(t) in the ramp regime, we need to compute EsErv~±sr(t)2, i.e. integrals of the following type:5.18 t4∬R2ρ(s)ρ(r)dsdr|∬-22dxdyeit(‖s‖x±‖r‖y)A(U,x,y)|2=t4∫-11dU1-U2∬-22dxdy∬-22dx′dy′A(U,x,y)A(U,x′,y′)¯×∫0∞eit(x-x′)σρ(σ)σdσ∫0∞e±it(y-y′)σ′ρ(σ′)σ′dσ′∼t2∫-11dU1-U2∬-22dxdy∬-22dx′dy′ddxA(U,x,y)ddy′A(U,x′,y′)¯×∫0∞eit(x-x′)σρ(σ)dσ∫0∞e±it(y-y′)σ′ρ(σ′)dσ′∼t2∫-11dU1-U2∬∬-22dxdydx′dy′ρ^(t(x-x′))ρ^(±t(y-y′))ddxA(U,x,y)ddy′A(U,x′,y′)¯. Here we performed two integrations by parts in x and y′ and ignored the boundary terms. Estimating the derivative of A as in (5.37), we can continue|(5.38)|≲t2∫-11dU1-U2∬-22dxdy(1-U)+|x-y|∬-22dx′dy′(1-U)+|x′-y′||ρ^(t(x-x′))ρ^(t(y-y′))|. The last two factors essentially restrict the integration to the regime |x-x′|≲1/t, |y-y′|≲1/t and by scaling we obtain a bound of order t1/2 for |(5.38)|. This completes the sketch of the proof of (2.28) in the radially symmetric case, the general case is analogous but technically more cumbersome and we omit the details. Acknowledgements We are grateful to the authors of [25] for sharing with us their insights and preliminary numerical results. We are especially thankful to Stephen Shenker for very valuable advice over several email communications. Helpful comments on the manuscript from Peter Forrester and from the anonymous referees are also acknowledged. Funding Information Open access funding provided by Institute of Science and Technology (IST Austria). Data Availability All data generated or analysed are included in this published article. 1 Equation (1.7) is for matrices whose second and fourth moments coincide with the ones of GUE, otherwise there are additional terms, see e.g. [13, Theorem 2.4]. 2 We assume equal fourth cumulants merely for notational convenience. Our proof verbatim covers also the more general case. 3 For the applications in this paper, SFF in the regime t≫1, the first term in (2.7) is the only relevant one. 4 The exponent in (2.8) can be optimized depending on τ and f. 5 The resolvent method extends to very general Hermitian matrices possibly with non-centered and correlated entries, see [19], but here we present only the Wigner case for simplicity. László Erdős: Partially supported by ERC Advanced Grant "RMTBeyond" No. 101020331. Dominik Schröder: Supported by Dr. Max Rössler, the Walter Haefner Foundation and the ETH Zürich Foundation. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. ==== Refs References 1. Bai ZD Yao J On the convergence of the spectral empirical process of Wigner matrices Bernoulli 2005 11 1059 1092 10.3150/bj/1137421640 2. Bao Z He Y On Cramér–von Mises statistic for the spectral distribution of random matrices Ann. Appl. Probab. 2022 32 4315 4355 10.1214/22-aap1788 3. Bao, Z., He, Y.: Quantitative CLT for linear eigenvalue statistics of Wigner matrices (2021). arXiv:2103.05402 4. Berry MV Semiclassical theory of spectral rigidity Proc. R. Soc. Lond. Ser. A 1985 400 229 251 10.1098/rspa.1985.0078 5. Bloemendal, A., Erdős, L., Knowles, A., Yau, H.-T., Yin, J.: Isotropic local laws for sample covariance and generalized Wigner matrices. Electron. J. Probab. (33), 53 (2014). 10.1214/ejp.v19-3054 6. Bohigas O Giannoni M-J Schmit C Characterization of chaotic quantum spectra and universality of level fluctuation laws Phys. Rev. Lett. 1984 52 1 4 10.1103/PhysRevLett.52.1 7. Bourgade, P., Erdős, L., Yau, H.-T., Yin, J.: Fixed energy universality for generalized Wigner matrices. Commun. Pure Appl. Math. 69, 1815–1881 (2016). 10.1002/cpa.21624 8. Boutet de Monvel A Khorunzhy A Asymptotic distribution of smoothed eigenvalue density. II. Wigner random matrices Random Oper. Stoch. Equ. 1999 7 149 168 10.1515/rose.1999.7.2.149 9. Brézin E Hikami S Spectral form factor in a random matrix theory Phys. Rev. E (3) 1997 55 4067 4083 10.1103/PhysRevE.55.4067 10. Cipolloni G Fluctuations in the spectrum of non-Hermitian i.i.d. matrices J. Math. Phys. 2022 63 053503 10.1063/5.0089089 11. Cipolloni, G., Erdős, L., Schröder, D.: Mesoscopic central limit theorem for non-Hermitian random matrices (2022). arXiv:2210.12060 12. Cipolloni, G., Erdős, L., Schröder, D.: Central limit theorem for linear eigenvalue statistics of non-Hermitian random matrices. Commun. Pure Appl. Math. (2019). 10.1002/cpa.22028 13. Cipolloni, G., Erdős, L., Schröder, D.: Functional central limit theorems for Wigner matrices (2020). arXiv:2012.13218 14. Cipolloni, G., Erdős, L., Schröder, D.: Functional central limit theorems for Wigner matrices. Ann. Appl. Probab. 33, 447–489 (2023). 10.1214/22-aap1820 15. Cipolloni, G., Erdős, L., Schröder, D.: Quenched universality for deformed Wigner matrices. Probab. Theory Relat. Fields (2021). 10.1007/s00440-022-01156-7 16. Cipolloni, G., Erdős, L., Schröder, D.: Thermalisation for Wigner matrices. J. Funct. Anal. 282, 109394 (2022). 10.1016/j.jfa.2022.109394 17. Cotler J Hunter-Jones N Liu J Yoshida B Chaos, complexity, and random matrices JHEP 2017 1711 2017 048 10.1007/JHEP11(2017)048 18. Cotler JS Gur-Ari G Hanada M Polchinski J Saad P Shenker SH Stanford D Streicher A Tezuka M Black holes and random matrices JHEP 2016 1705 118 19. Erdős L Krüger T Schröder D Random matrices with slow correlation decay Forum Math. Sigma 2019 7 e8 10.1017/fms.2019.2 20. Erdős L Yau H-T Yin J Rigidity of eigenvalues of generalized Wigner matrices Adv. Math. 2012 229 1435 1515 10.1016/j.aim.2011.12.010 21. Forrester PJ Differential identities for the structure function of some random matrix ensembles J. Stat. Phys. 2021 183 33 10.1007/s10955-021-02767-5 22. Forrester PJ Quantifying dip–ramp–plateau for the Laguerre unitary ensemble structure function Commun. Math. Phys. 2021 387 215 235 10.1007/s00220-021-04193-w 23. García-García AM Jia Y Verbaarschot JJM Universality and Thouless energy in the supersymmetric Sachdev–Ye–Kitaev model Phys. Rev. D 2018 97 106003 10.1103/physrevd.97.106003 24. García-García AM Verbaarschot JJM Analytical spectral density of the Sachdev–Ye–Kitaev model at finite N Phys. Rev. D 2017 96 066012 10.1103/PhysRevD.96.066012 25. Gharibyan, H., Pattison, C., Shenker, S., Wells, K.: Work in preparation (2021) 26. Guionnet A Large deviations upper bounds and central limit theorems for non-commutative functionals of Gaussian large random matrices Ann. Inst. H. Poincaré Probab. Stat. 2002 38 341 384 10.1016/S0246-0203(01)01093-7 27. Gutzwiller, M.C.: Chaos in Classical and Quantum Mechanics, vol. 1. Interdisciplinary Applied Mathematics, pp. xiv+432. Springer, New York (1990) 28. He Y Mesoscopic linear statistics of Wigner matrices of mixed symmetry class J. Stat. Phys. 2019 175 932 959 10.1007/s10955-019-02266-8 29. He Y Knowles A Mesoscopic eigenvalue density correlations of Wigner matrices Probab. Theory Relat. Fields 2020 177 147 216 10.1007/s00440-019-00946-w 30. He Y Knowles A Mesoscopic eigenvalue statistics of Wigner matrices Ann. Appl. Probab. 2017 27 1510 1550 10.1214/16-AAP1237 31. Heusler S Müller S Altland A Braun P Haake F Periodic-orbit theory of level correlations Phys. Rev. Lett. 2007 98 044103 10.1103/PhysRevLett.98.044103 17358777 32. Jia Y Verbaarschot JJM Spectral fluctuations in the Sachdev–Ye–Kitaev model J. High Energy Phys. 2020 193 57 10.1007/jhep07(2020)193 33. Johansson K On fluctuations of eigenvalues of random Hermitian matrices Duke Math. J. 1998 91 151 204 10.1215/S0012-7094-98-09108-6 34. Khorunzhy AM Khoruzhenko BA Pastur LA Asymptotic properties of large random matrices with independent entries J. Math. Phys. 1996 37 5033 5060 10.1063/1.531589 35. Knowles A Yin J The isotropic semicircle law and deformation of Wigner matrices Commun. Pure Appl. Math. 2013 66 1663 1750 10.1002/cpa.21450 36. Landon, B., Lopatto, P., Sosoe, P.: Single eigenvalue fluctuations of general Wigner-type matrices (2021). arXiv:2105.01178 37. Landon, B., Sosoe, P.: Almost-optimal bulk regularity conditions in the CLT for Wigner matrices (2022). arXiv:2204.03419 38. Landon B Sosoe P Applications of mesoscopic CLTs in random matrix theory Ann. Appl. Probab. 2020 30 2769 2795 10.1214/20-AAP1572 39. Landon, B., Sosoe, P., Yau, H.-T.: Fixed energy universality for Dyson Brownian motion (2016). arXiv:1609.09011 40. Leviandier L Lombardi M Jost R Pique J Fourier transform: a tool to measure statistical level properties in very complex spectra Phys. Rev. Lett. 1986 56 2449 2452 10.1103/PhysRevLett.56.2449 10032995 41. Lytova A Pastur L Fluctuations of matrix elements of regular functions of Gaussian random matrices J. Stat. Phys. 2009 134 147 159 10.1007/s10955-008-9665-1 42. Mehta, M.L.: Random Matrices, Third, Vol. 142, Pure and Applied Mathematics (Amsterdam), pp. xviii+688. Elsevier, Amsterdam (2004) 43. Müller S Heusler S Braun P Haake F Altland A Semiclassical foundation of universality in quantum chaos Phys. Rev. Lett. 2004 93 014103 10.1103/PhysRevLett.93.014103 44. Okuyama, K.: Spectral form factor and semi-circle law in the time direction. J. High Energy Phys. 161, front matter + 15 (2019). 10.1007/jhep02(2019)161 45. Prange RE The spectral form factor is not self-averaging Phys. Rev. Lett. 1997 78 2280 2283 10.1103/PhysRevLett.78.2280 46. Saad, P., Shenker, S.H., Stanford, D.: A semiclassical ramp in SYK and in gravity (2018). arXiv:1806.06840 47. Shcherbina, M.: Central limit theorem for linear eigenvalue statistics of the Wigner and sample covariance random matrices. Zh. Mat. Fiz. Anal. Geom. 7, 176–192, 197, 199 (2011) 48. Sieber M Richter K Correlations between periodic orbits and their role in spectral statistics Phys. Scr. 2001 T90 128 10.1238/physica.topical.090a00128 49. Sosoe P Wong P Regularity conditions in the CLT for linear eigenvalue statistics of Wigner matrices Adv. Math. 2013 249 37 87 10.1016/j.aim.2013.09.004 50. Watson, G.N.: A Treatise on the Theory of Bessel Functions. Cambridge Mathematical Library. Reprint of the second edition (1944), pp. viii+804. Cambridge University Press, Cambridge (1995)