Eigenvalue fluctuations of Gaussian random matrices

# Replica-symmetric approach to the typical eigenvalue fluctuations of Gaussian random matrices

Fernando L. Metz Instituto de Física, Universidade Federal do Rio Grande do Sul, 91501-970 Porto Alegre, Brazil
Departamento de Física, Universidade Federal de Santa Maria, 97105-900 Santa Maria, Brazil
London Mathematical Laboratory, 14 Buckingham Street, London WC2N 6DF, United Kingdom
July 26, 2019
###### Abstract

We discuss an approach to compute the first and second moments of the number of eigenvalues that lie in an arbitrary interval of the real line for Gaussian random matrices. The method combines the standard replica-symmetric theory with a perturbative expansion of the saddle-point action up to (), leading to the correct logarithmic scaling of the variance as well as to an analytical expression for the correction to the average . Standard results for the number variance at the local scaling regime are recovered in the limit of a vanishing interval. The limitations of the replica-symmetric method are unveiled by comparing our results with those derived through exact methods. The present work represents an important step to study the fluctuations of in non-invariant random matrix ensembles, where the joint distribution of eigenvalues is not known.

## 1 Introduction

In the last decades, random matrices have established itself as a fundamental tool to model complex disordered systems, with many applications in different branches of science [BookMehta, oxford]. The primary aim in random matrix theory is the study of the eigenvalue statistics, which is accessed by computing suitably chosen observables. Important observables are the moments of the number of eigenvalues that lie between and . For , the random variable is the so-called index, i.e., the number of eigenvalues contained in the unbounded interval . A great deal of attention has been devoted to the index of random matrices, especially due to its ability to probe the stability properties of the energy landscape characterising disordered systems [Cavagna2000, wales2003energy]. Another prominent observable, derived from and studied originally by Dyson [Dyson2], is the variance of the number of eigenvalues within the bounded interval (), also known as the number variance. Recently, this observable has been applied to study the statistics of non-interacting fermions confined in a harmonic trap [Marino2014, Marino2016], and the ergodicity of the eigenfunctions in a mean-field model for the Anderson localisation transition [MetzLett, Metz2017].

A powerful tool to study the fluctuations of is the Coulomb gas approach [Dyson1], whose starting point is an exact correspondence between the joint distribution of eigenvalues and the partition function of a two-dimensional Coulomb gas confined to a line. The Coulomb gas technique has been employed to derive analytical results for typical and atypical fluctuations of in the case of rotationally invariant random matrices (RIRM), including Gaussian [Marino2014, Marino2016, Majumdar1, Majumdar2, Perez2014b, Isaac2016], Wishart [Vivo, Perez2015] and Cauchy [Majumdar3] random matrices. Restricting ourselves to typical fluctuations of , one of the central outcomes for RIRM is the logarithmic increase of the variance as a function of the matrix dimension , due to the strong correlations among the eigenvalues. In spite of the success of the Coulomb gas method, its application is strictly limited to RIRM, where the joint distribution of eigenvalues is analytically known and, consequently, the analogy with the Coulomb gas partition function can be readily established.

Recently, a novel technique to study the fluctuations of has been developed [MetzLett, Metz2015]. This method is based on the replica-symmetric theory of disordered systems and its main advantage lies in its broad range of possible applications, which goes beyond the realm of RIRM. In fact, the replica approach allows to derive analytical results for typical and atypical fluctuations of in the case of random matrices where the joint distribution of eigenvalues is not even known. The most typical example in this sense is the adjacency matrix of sparse random graphs [MetzLett, Metz2015], where the eigenvalues behave as uncorrelated random variables and the leading term of the variance scales as for large . However, it is unclear whether the formalism of [MetzLett, Metz2015] is able to grasp the variance behaviour of random matrices with strong correlated eigenvalues. In the case of a positive answer, this would pave the way to inspect the fluctuations of in random matrices that are not rotationally invariant, but in which the eigenvalues are strongly correlated. The ensembles of Gaussian random matrices are the ideal testing ground for this matter, since one can make detailed comparisons with exact results and identify eventual limitations of the replica method.

In this work we show how the replica approach, as developed in [Metz2015, MetzLett], has to be further adapted in order to derive the correct logarithmic scaling for the GOE and the GUE ensembles of random matrices [Dyson1, BookMehta]. Following the approach of [MetzPar], the central idea here consists in writing the characteristic function of as a saddle-point integral, in which the action is computed perturbatively around its limit. We show that the leading term () is correctly recovered only when the fluctuations due to finite are taken into account. However, the present approach does not yield the exact expression for the contribution to the variance of , due to our assumption of replica-symmetry for the order-parameter. In the case of the number variance, our analytical expression converges, in the limit , to standard results of random matrix theory [BookMehta], valid in the regime where the spectral window is measured in units of the mean level spacing. This result supports the central claim of reference [Metz2017], where the limit of the number variance is employed to study the ergodic nature of the eigenfunctions in the Anderson model on a regular random graph. As a by-product, our method yields the correction to the intensive average , whose exactness is tested against numerical diagonalisation and previous analytical results. While the contribution to is exact for the GOE ensemble, it fails in the case of the GUE ensemble due to our assumption of replica symmetry. The present approach can be also employed to compute the moments of characteristic polynomials of non-invariant random matrices, where the key quantity under study is similar to eq. (7). The moments of characteristic polynomials of invariant random matrices have been largely studied [Brezin, Brezin1, Borodin], in view of their connection with the distribution of zeros of the Riemann zeta function [keating].

The paper is organised as follows. In the next section we show how the computation of the characteristic function of can be pursued using the replica approach. Section 3 explains the basic steps of the replica method, including the perturbative calculation of the action up to . We make the replica-symmetric assumption for the order-parameter and derive an analytical expression for the characteristic function in section 4. Section 5 derives the analytical results for the index, the number variance and the fluctuations of in an arbitrary interval. The last section summarises our work and discusses the impact of replica symmetry in our results.

## 2 The number of eigenvalues inside an interval

We study here two different ensembles of Gaussian random matrices with real eigenvalues , defined through the probability distribution . If the elements of the ensemble are real symmetric matrices, we have the Gaussian orthogonal ensemble (GOE) [BookMehta]

 P(\boldmathM)=NGOEexp(−N4\Tr\boldmathM2). (1)

If the ensemble is composed of complex-Hermitian random matrices, we have the Gaussian unitary ensemble (GUE) [BookMehta]

 P(\boldmathM)=NGUEexp[−N2\Tr(\boldmathM†\boldmathM)], (2)

where represents the conjugate transpose of a matrix. The explicit form of the normalization factors in eqs. (1) and (2) are not important in our computation.

The number of eigenvalues lying between and is given by

 IN(a,b)=N∑α=1[Θ(b−λα)−Θ(a−λα)]a

with the Heaviside step function. The statistics of is encoded in the characteristic function

 GN(μ)=⟨exp[iμIN(a,b)]⟩, (4)

where is the ensemble average over the random matrix elements. In particular, the first two moments of read

 ⟨IN(a,b)⟩=−i∂GN(μ)∂μ∣∣μ=0,⟨I2N(a,b)⟩=−∂2GN(μ)∂μ2∣∣μ=0. (5)

In order to calculate the ensemble average in eq. (4), one has to write in terms of the random matrix . By following [Cavagna2000, Metz2015] and representing as the discontinuity of the principal value of the complex logarithm along the negative real axis, may be written as the limit

 IN(a,b)=12πilimη→0+ln[Z(zb)Z(z∗a)Z(z∗b)Z(za)],za=a+iη,zb=b+iη, (6)

with . Here is an arbitrary complex number, denotes its complex-conjugate, and is the identity matrix. Since the imaginary part of the principal logarithm is bounded in , the right hand side of eq. (6) is not extensive and, consequently, unfit to count the number of eigenvalues within for single realizations of . This issue comes from our naive derivation of eq. (6), where we assume that the principal complex logarithm fulfills the same standard properties as those valid for the logarithm of real numbers. In spite of that, this is a necessary step to apply the replica method. We will see that, after calculating the ensemble average and introducing an appropriate order-parameter, the problem factorises over sites and the extensivity of the moments of is restored. This procedure is heuristic, but it yields correct results for the moments of for large . Thus, eq. (6) completely encodes the statistics of , even though it is unsuitable to obtain for single instances of .

Inserting eq. (6) in eq. (4), we find

 GN(μ)=limη→0+⟨[Z(zb)Z(z∗a)]μ2π[Z(z∗b)Z(za)]−μ2π⟩. (7)

At this point we invoke the replica method in the form

 GN(μ)=limη→0+limn±→±μ2πGN(n±,η), (8)

in which we have introduced the function

 GN(n±,η)=⟨[Z(zb)Z(z∗a)]n+[Z(z∗b)Z(za)]n−⟩ (9)

for finite . The idea consists in assuming that are positive integers, which allows to calculate for through a saddle-point integration. After we have derived the behaviour of for , we take the replica limit and reconstruct the original function in the limit . This is nothing more than the general strategy of the replica approach [BookParisi]. The only difference lies in the fact that we continue the arbitrary positive integers to purely imaginary numbers.

## 3 Replica method and finite size corrections

Our aim in this section is to calculate for large but finite . Firstly, we have to recast eq. (9) into an exponential form, which is suitable to perform the ensemble average. This is achieved by representing the functions and as Gaussian integrals. For instance, the function is written as the multidimensional Gaussian integral [negele]

 Z(za)=1det(\boldmathM−\boldmathINza)=1(2π)N∫(N∏i=1dϕidϕ∗i)exp[−iN∑ij=1ϕ∗i(\boldmathM−%\boldmath$I$Nza)ijϕj], (10)

where are complex integration variables. This representation of is appropriate to deal with both the GOE ensemble and the GUE ensemble in the same framework. Introducing analogous identities for , and , eq. (9) can be compactly written as

 GN(n±,η)=∫(N∏i=1d\boldmathφid\boldmathφ†i)exp(iN∑i=1\boldmathφ†i.\boldmathZ% \boldmathφi)⟨exp[−iN∑ij=1Mij(\boldmathφ†i.\boldmathA% \boldmathφj)]⟩, (11)

where the following block matrices have been introduced

 \boldmathA=⎛⎜ ⎜ ⎜ ⎜⎝\boldmathI+0000\boldmathI−0000−\boldmathI+0000−\boldmathI−⎞⎟ ⎟ ⎟ ⎟⎠,\boldmathZ=⎛⎜ ⎜ ⎜ ⎜⎝zb\boldmathI+0000za\boldmathI−0000−z∗a\boldmathI+0000−z∗b\boldmathI−⎞⎟ ⎟ ⎟ ⎟⎠,

with () denoting the () identity matrix. The integration variables in eq. (11) are -dimensional complex vectors in the replica space, defined according to

 \boldmathφi=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝ϕ1,iψ1,iϕ2,iψ2,i⎞⎟ ⎟ ⎟ ⎟ ⎟⎠i=1,…,N,

where and have dimension , while and have dimension . Each one of the vectors comes from the Gaussian integral representation of a single function in eq. (9). The off-diagonal blocks of and , filled with zeros, have suitable dimensions, such that their product with is well-defined.

The ensemble average in eq. (11) is easily performed for the GOE and the GUE ensembles by using eqs. (1) and (2), leading to

 GN(n±,η)=∫(N∏i=1d\boldmathφid\boldmathφ†i)exp[iN∑i=1\boldmathφ†i.\boldmathZ% \boldmathφi−12NN∑ij=1Kβ(\boldmathφi,\boldmathφ†i,\boldmathφ% j,\boldmathφ†j)], (12)

with the kernel

 Kβ(\boldmathφi,\boldmathφ†i,\boldmathφj,\boldmathφ†j)=|\boldmathφ†i.\boldmathA\boldmathφj|2+(2−β)Re(\boldmathφ†i.\boldmathA\boldmathφj)2, (13)

which depends on the random matrix ensemble through the Dyson index : for the GOE ensemble and for the GUE ensemble. Following the standard approach to decouple the sites in eq. (12), we introduce the order-parameter

 ρ(\boldmathφ,\boldmathφ†)=1NN∑i=1δ(\boldmathφ−\boldmathφi)δ(\boldmathφ†−\boldmathφ†i) (14)

through a functional Dirac delta, yielding

 GN(n±,η) = ∫DρD^ρexp[i∫d% \boldmathφd\boldmathφ†ρ(\boldmathφ,\boldmathφ†)^ρ(\boldmathφ,\boldmathφ†)] (15) × exp[−N2∫d\boldmathφ1d% \boldmathφ†1d\boldmathφ2d% \boldmathφ†2ρ(\boldmathφ1,% \boldmathφ†1)Kβ(\boldmathφ1,\boldmathφ†1,\boldmathφ2,% \boldmathφ†2)ρ(\boldmathφ2,% \boldmathφ†2)] ×

where is the functional integration measure over and its conjugate order-parameter . After performing the Gaussian integral over and introducing a new integration variable

 ^ρ(\boldmathφ,\boldmathφ†)=−iN∫d\boldmathφ1d\boldmathφ†1Kβ(\boldmathφ,\boldmathφ†,% \boldmathφ1,\boldmathφ†1)Φ(% \boldmathφ1,\boldmathφ†1), (16)

we obtain the compact expression

 GN(n±,η)=√detKβ∫DΦexp(NSN[Φ]), (17)

with the action:

 SN[Φ] = 12∫d\boldmathφ1d\boldmathφ†1d\boldmathφ2d\boldmathφ% †2Φ(\boldmathφ1,\boldmathφ†1)Kβ(\boldmathφ1,\boldmathφ% †1,\boldmathφ2,\boldmathφ†2)Φ(\boldmathφ2,\boldmathφ†2) (18) + ln[∫d\boldmathφd\boldmathφ†exp(i\boldmathφ†.\boldmathZ% \boldmathφ−∫d\boldmathφ1d\boldmath% φ†1Kβ(\boldmathφ,\boldmathφ†,\boldmathφ1,\boldmathφ†1)Φ(\boldmathφ1,\boldmathφ†1))].

The functional integration measure in eq. (17) can be intuitively written as

 DΦ=∏\boldmathφ\boldmathφ†(−i)√N2πdΦ(\boldmathφ,% \boldmathφ†), (19)

in which the product runs over all possible arguments of the function .

The next step consists in solving the integral in eq. (17) using the saddle-point method. In the limit , this integral is dominated by the value that extremizes the action , i.e.

 δSN[Φ]δΦ(\boldmathφ,\boldmathφ†)∣∣∣Φ=Φ0=0, Φ0(\boldmathφ,\boldmathφ†)=exp[i\boldmathφ†.% \boldmathZ\boldmathφ−∫d\boldmathφ1d\boldmathφ†1Kβ(\boldmathφ,\boldmathφ†,\boldmathφ1,% \boldmathφ†1)Φ0(\boldmathφ1,%\boldmath$φ$†1)]∫d\boldmathφ2d\boldmathφ†2exp[i\boldmathφ% †2.\boldmathZ\boldmathφ2−∫d% \boldmathφ1d\boldmathφ†1Kβ(%\boldmath$φ$2,\boldmathφ†2,% \boldmathφ1,\boldmathφ†1)Φ0(%\boldmath$φ$1,\boldmathφ†1)]. (20)

We are not interested only on the behaviour of strictly in the limit , but also on the first perturbative correction due to large but finite . We will see that such correction yields precisely the correct logarithmic scaling of the variance of with respect to . It is important to note that, up to now, we have not done any approximation regarding the system size dependence. In other words, eqs. (17) and (18) are exact for finite . As we can see from eq. (18), the only source of finite size fluctuations in the action comes from fluctuations of the order-parameter around its limit . Let us assume that such fluctuations are of the form [MetzPar]

 Φ(\boldmathφ,\boldmathφ†)=Φ0(\boldmathφ,\boldmathφ†)−1√Nχ(\boldmathφ,\boldmathφ†). (21)

After expanding around up to , we substitute eq. (21) and rewrite eq. (17) as follows

 GN(n±,η)=√detKβexp(NS0[Φ0]) ×∫Dχexp[12∫d% \boldmathφ1d\boldmathφ†1d% \boldmathφ2d\boldmathφ†2χ(% \boldmathφ1,\boldmathφ†1)H(% \boldmathφ1,\boldmathφ†1,% \boldmathφ2,\boldmathφ†2)χ(% \boldmathφ2,\boldmathφ†2)], (22)

in which

 H(\boldmathφ1,\boldmathφ†1,\boldmathφ2,\boldmathφ†2)=δ2SNδΦ(\boldmathφ1,% \boldmathφ†1)δΦ(\boldmathφ2,\boldmathφ†2)∣∣∣Φ=Φ0. (23)

The leading contribution to the action, defined as , is formally given by eq. (18) with replaced by , while the functional integration measure in eq. (22) reads

 Dχ=∏\boldmathφ,\boldmathφ†i√12πdχ(\boldmathφ,\boldmath% φ†). (24)

By computing explicitly the derivatives in eq. (23) and using eq. (20), the matrix of second derivatives can be expressed as

 \boldmathH=\boldmathKβ+\boldmathKβ\boldmathT, (25)

with the elements of defined according to

 T(\boldmathφ1,\boldmathφ†1,\boldmathφ2,\boldmathφ†2) = Kβ(\boldmathφ1,\boldmathφ% †1,\boldmathφ2,\boldmathφ†2)Φ0(\boldmathφ1,\boldmathφ†1) (26) − Φ0(\boldmathφ1,\boldmathφ†1)∫d\boldmathφd\boldmathφ†Kβ(\boldmathφ2,\boldmathφ†2,\boldmathφ,\boldmathφ†)Φ0(\boldmathφ,\boldmathφ†).

Now we are ready to obtain an expression for when is large but finite. After integrating over the Gaussian fluctuations in eq. (22), we substitute eq. (25) and derive

 GN(n±,η)=exp(NS0[Φ0]+12∞∑ℓ=1(−1)ℓℓTr\boldmathTℓ). (27)

The above equation is determined essentially by the behaviour of , which is obtained from the solution of eq. (20). In the next section we show the outcome for when is symmetric with respect to the interchange of the replica indexes.

## 4 The replica symmetric characteristic function

The next step in our calculation of the characteristic function consists in performing the replica limit of eq. (27). Thus, we need to understand how depends on , which is ultimately determined by the solutions of eq. (20). In order to proceed further, we follow previous works [Kuhn2008, MetzPar, Metz2015, MetzLett] and we make the following Gaussian ansatz for the order-parameter

 Φ0(\boldmathφ,\boldmathφ†)=det\boldmathC(2πi)2(n++n−)exp(−% \boldmathφ†.\boldmathC\boldmathφ), (28)

which is parametrised by the diagonal block matrix

 \boldmathC=⎛⎜ ⎜ ⎜ ⎜⎝\boldmathI+Δ∗b0000\boldmathI−Δ∗a0000\boldmathI+Δa0000\boldmathI−Δb⎞⎟ ⎟ ⎟ ⎟⎠,ReΔa>0ReΔb>0,

given in terms of the complex parameters and . The conditions and ensure the convergence of any Gaussian integrals over . The off-diagonal blocks in have suitable dimensions such that the standard matrix operations involving and are well-defined. The above assumption for is symmetric with respect to the permutation of replicas inside each group and . We do not consider here solutions of eq. (20) that break replica symmetry.

Let us explore the consequences of the replica symmetric (RS) assumption. Substituting eq. (28) in eq. (20), considering the explicit form of the kernel (see eq. (13)), and noting that

 ∫d\boldmathφ1d\boldmathφ†1|\boldmathφ†.\boldmathA% \boldmathφ1|2Φ0(\boldmathφ1,% \boldmathφ†1)=\boldmathφ†.% \boldmathC−1\boldmathφ, ∫d\boldmathφ1d\boldmathφ†1Re(\boldmathφ†.\boldmathA\boldmathφ1)2Φ0(\boldmathφ1,\boldmathφ†1)=0, (29)

we conclude that the RS form of solves the self-consistent eq. (20) provided the parameters and fulfill the quadratic equations

 Δ2a−iz∗aΔa−1=0,Δ2b−iz∗bΔb−1=0. (30)

Now we are in a position to derive the explicit dependency of with respect to . Inserting the RS assumption for in eq. (18), the leading contribution to the action is derived

 S0(n±)=12⎡⎣n+(Δ∗b)2+n−(Δ∗a)2+n+Δ2a+n−Δ2b⎤⎦−n+ln(ΔaΔ∗b)−n−ln(Δ∗aΔb). (31)

The second contribution appearing in eq. (27) involves an infinite series, so that we have to evaluate the RS form of the coefficients . Plugging eq. (28) in eq. (26) and performing the Gaussian integrals, we have derived the following expression

 Tr\boldmathTℓ = (2−β)⎡⎣n+(Δ∗b)2ℓ+n−(Δ∗a)2ℓ+n+Δ2ℓa+n−Δ2ℓb⎤⎦ (32) + 2β⎡⎣n+(Δ∗b)ℓ+n−(Δ∗a)ℓ+(−1)ℓn+Δℓa+(−1)ℓn−Δℓb⎤⎦2,

in which we have used the fact that the Dyson index is limited to the values or . Finally, eqs. (31) and (32) are substituted in eq. (27), the limit is taken, and the following expression for the characteristic function is obtained

 (33)

where and are, respectively, the mean and the variance of for finite

 ⟨IN⟩ηN = i4π[1Δ2b−1(Δ∗b)2−1Δ2a+1(Δ∗a)2]−i2πln(ΔbΔ∗aΔ∗bΔa) − i(2−β)4πN∞∑ℓ=1(−1)ℓℓ[1Δ2ℓa−1(Δ∗a)2ℓ−1Δ2ℓb+1(Δ∗b)2ℓ], ⟨I2N⟩η − ⟨IN⟩2η=−12βπ2∞∑ℓ=1(−1)ℓℓ[(−1)ℓΔℓa−1(Δ∗a)ℓ−(−1)ℓΔℓb+1(Δ∗b)ℓ]2. (35)

Note that the contribution of in the exponent of eq. (27) depends linearly on , only providing the mean value . If one wants to extract any information about the typical fluctuations around , one has to take into account the finite size fluctuations of the order-parameter . This is in contrast with some models of sparse random matrices, where the calculation of the leading term of the action is enough to obtain the linear scaling of the variance of with [Metz2015]. As we will see below, here the system size dependence of manifests itself as the divergence of the infinite series present in eq. (35). The finite-size fluctuations of also yield the correction to appearing in eq. (4).

## 5 The mean and the variance of In

In this section we derive explicit analytical results from eqs. (4) and (35). The solutions of eqs. (30) read

 Δa=12(iz∗a±√4−(z∗a)2),Δb=12(iz∗b±√4−(z∗b)2). (36)

We consider below the behaviour of and in the limit for specific observables depending on the values of and . For this purpose, it is instrumental to recognise that the eigenvalues of the GOE and the GUE random matrix ensembles are distributed, for , according to the Wigner semicircle law with support in [BookMehta].

### 5.1 The index

The first observable considered here is the index, i.e., the number of eigenvalues smaller than a certain threshold . In this case we set and , which leads to the following solutions for

 Δa = |Δa|e−iπ2, (37) Δb = eiθb, (38)

with the argument:

 θb=arctan(b√4−b2). (39)

After plugging eqs. (37) and (38) in eq. (4), we can sum the convergent series and obtain

 ⟨IN⟩N=12+12πsin(2θb)+1πθb+(2−β)2NC(θb), (40)

where

 (41)

The leading term of eq. (40) agrees with an exact result [Perez2014b]. The coefficient can be derived from the integral , where is the correction to the average spectral density. In the case of (GOE), eq. (40) is in full agreement with references [verbaar1984, dhesi1990, Itoi1997, Kalisch2002], where is computed exactly using replicas [verbaar1984, dhesi1990, Itoi1997] and supersymmetry [Kalisch2002]. In figure 1 we also compare the analytical result for with numerical diagonalisation results for real symmetric matrices drawn from the GOE ensemble. The agreement between our analytical expression and numerical diagonalisation is very good until the band edge is approached, where finite size fluctuations become stronger and a discrepancy between theory and numerics is evident. The results shown on the inset exhibit the convergence of the numerical results to the analytical formula as increases.

For (GUE), the correction in eq. (40) is zero. This is in contrast to the exact expression for obtained through supersymmetry [Kalisch2002], replicas [Kamenev] and the large expansion of orthogonal polynomials [Kamenev]. These methods show that is an oscillatory function of with maxima. For large but finite , the integration of over an interval within the support of the spectral density yields a negligible contribution to the correction of . This has been confirmed by computing through numerical diagonalisation and then subtracting the leading term of eq. (40), which yields an outcome orders of magnitude smaller than . The present approach is not able to recover the correct contribution to due to our replica symmetric assumption for the order parameter. This is evident from the calculation of the spectral density using the fermionic replica method [Kamenev], in which the oscillatory correction is derived by including saddle-points that break replica symmetry.

Let us derive an expression for the variance of the index. Inserting the above expressions for and in eq. (35) and performing the sum of the convergent series, we can write the formal expression

 ⟨I2N⟩−⟨IN⟩2=1βπ2[∞∑ℓ=11ℓ+12ln(2+2cos(2θb))]. (42)

We have isolated in eq. (42) the divergent contribution to the variance in the form of the harmonic series. This is not surprising, since the leading term of should indeed diverge for . The question here is how scales with . In order to extract this behaviour in the replica framework, the authors of [Cavagna2000] keep the regularizer finite until the end of the calculation, and then assume there is a functional relation between and . Here a different strategy has been pursued, where the limit has been performed before considering the convergence of the series in eqs. (4) and (35). This approach gives rise to the divergent contribution in eq. (42), which is naturally interpreted as the leading term . Thus, in order to understand how the variance scales with , we have to study how the harmonic series diverges. For large , the partial summation behaves as

 N∑ℓ=11ℓ=lnN+γ+O(1/N), (43)

where is the Euler-Mascheroni constant. Consequently, we conclude that the variance behaves for as follows

 ⟨I2N⟩−⟨IN⟩2=1βπ2[lnN+γ+12ln(2+2cos(2θb))]. (44)

The above equation recovers exactly the leading behaviour of the index variance for [Majumdar1, Majumdar2, Isaac2016, Perez2014b]. However, eq. (44) fails in reproducing the correction to . For and , the term in eq. (44) is given by , which is only part of the exact result for the GUE ensemble [Majumdar2]. For , the term of eq. (44) does not agree as well with the available result of [Cavagna2000], obtained from a fitting of numerical diagonalisation results. As we shall discuss below, this inaccuracy comes from the replica-symmetric assumption for the order-parameter.

### 5.2 The number of eigenvalues in a symmetric interval

For the second observable we set and , with . In this case the random variable quantifies the number of eigenvalues within . The solutions for and in the limit read

 Δa=Δ∗b, (45) Δb=eiθL,θL=arctan(L√4−L2). (46)

Inserting the above forms in eq. (4) and summing the series we obtain

 ⟨IN⟩N=1πsin(2θL)+2πθL+(2−β)NC(θL), (47)

with defined in eq. (41).

The leading term of eq. (47) agrees with the exact result [Marino2014]. For , the correction in the above equation is absent, due to the replica symmetric ansatz for the order parameter. The contribution to can be derived by integrating the correction to the spectral density over . In the replica framework, the latter quantity is exactly calculated only when replica symmetry breaking is taken into account [Kamenev], i.e., the situation here is completely analogous to the case of the index discussed above. For , our result for the correction in eq. (47) can be derived from references [verbaar1984, dhesi1990, Itoi1997, Kalisch2002], where the contribution to the spectral density is computed exactly using replicas [verbaar1984, dhesi1990, Itoi1997] and supersymmetry [Kalisch2002]. In figure 2 we compare the analytical behaviour of with numerical diagonalisation of real symmetric matrices drawn from the GOE ensemble. Similarly to the correction to the average index, figure 2 illustrates the convergence of the numerical diagonalisation results to the analytical formula of for increasing . However, the variance of for Gaussian random matrices displays an abrupt change of behaviour as we reach the scaling regime [Marino2014, Marino2016], which indicates that our formula for the coefficient might in fact breakdown sufficiently close to . A detailed analysis of the behaviour close to the band edges will not be pursued here.

The variance of the number of eigenvalues within is the so-called number variance [BookMehta]. The substitution of eqs. (45) and (46) in eq. (35) reads

 ⟨I2N⟩−⟨IN⟩2=2βπ2[∞∑ℓ=11ℓ+12ln(sin2(2θL))], (48)

where the divergent contribution appears once more as a harmonic series. Following the reasoning of the previous subsection, we conclude that the number variance for