Spectral statistics of the uni-modular ensemble

TL;DR

This paper studies the spectral properties of large Hermitian matrices with elements uniformly distributed on the unit circle, employing three methods to analyze spectral moments and fluctuations, and compares their effectiveness.

Contribution

It introduces three complementary analytical methods to study the spectral statistics of the uni-modular ensemble and compares their advantages, providing new insights into spectral fluctuations.

Findings

Derived expressions for 1/N corrections to spectral moments

Analyzed fluctuations around the mean spectral moments

Compared the effectiveness of supersymmetric, combinatorial, and Brownian motion methods

Abstract

We investigate the spectral statistics of Hermitian matrices in which the elements are chosen uniformly from U (1), called the uni-modular ensemble (UME), in the limit of large matrix size. Using three complimentary methods; a supersymmetric integration method, a combinatorial graph-theoretical analysis and a Brownian motion approach, we are able to derive expressions for 1/N corrections to the mean spectral moments and also analyse the fluctuations about this mean. By addressing the same ensemble from three different point of view, we can critically compare their relative advantages and derive some new results.

Equations

$$M_{\mu \nu }=(1-{\delta }_{\mu \nu })\exp \{i{\varphi }_{\mu \nu }\}\mathrm{and}{\varphi }_{\mu \nu }=-{\varphi }_{\nu \mu },$$

$$P(\mathcal{H})=\frac{1}{Z_N}\exp \left\{-\frac{1}{2}\mathrm{Tr}({\mathcal{H}}^2)\right\}$$

$${\tau }_k(\mathcal{H})=\int d\lambda {\lambda }^k\left(\frac{1}{N}\sum _{i=1}^N\delta (\lambda -{\lambda }_i(\mathcal{H}))\right)=\frac{1}{N}\mathrm{Tr}({\mathcal{H}}^k).$$

$$m_k:=⟨{\tau }_{2k}\left(\frac{\mathcal{H}}{\sqrt{N}}\right)⟩=⟨\frac{1}{N^{k+1}}\text{Tr}({\mathcal{H}}^{2k})⟩.$$

$$(k+1)m_k=(4k-2)m_{k-1}+\frac{(k-1)(2k-1)(2k-3)}{N^2}{m}_{k-2}.$$

$$m_k=C_k\left(1+\frac{1}{N^2}\frac{(k-1)k(k+1)}{12}+\mathcal{O}(N^{-4})\right).$$

$$L_f(\mathcal{H}):=\mathrm{Tr}[f(\mathcal{H})]-⟨\mathrm{Tr}[f(\mathcal{H})]⟩,$$

$$⟨H_{\mu \nu }{H}_{{\nu }^{\prime }{\mu }^{\prime }}⟩={\delta }_{\mu {\mu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}.$$

$$⟨M_{\mu \nu }{M}_{{\nu }^{\prime }{\mu }^{\prime }}⟩={\delta }_{\mu {\mu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}(1-{\delta }_{\mu \nu }).$$

$$⟨\mathrm{Tr}({\mathcal{H}}^n)⟩=⟨\mathrm{Tr}({\mathcal{M}}^n)⟩\mathrm{for}\mathrm{all}\mathrm{integer}n\ge 1.$$

$$⟨\mathrm{Tr}({\mathcal{H}}^{2k})⟩=C_k{N}^{k+1}+\dots .$$

$$⟨M_{\mu \nu }{M}_{{\mu }^{\prime }{\nu }^{\prime }}{M}_{{\mu }^{\prime \prime }{\nu }^{\prime \prime }}{M}_{{\mu }^{\prime \prime \prime }{\nu }^{\prime \prime \prime }}{⟩}_4=[{\delta }_{\mu {\nu }^{\prime }}{\delta }_{\nu {\mu }^{\prime }}(1-{\delta }_{\mu \nu })]$$

$$\times [{\delta }_{{\mu }^{\prime \prime }{\nu }^{\prime \prime \prime }}{\delta }_{{\nu }^{\prime \prime }{\mu }^{\prime \prime \prime }}(1-{\delta }_{{\mu }^{\prime \prime }{\nu }^{\prime \prime }})]{\delta }_{\mu {\mu }^{\prime \prime }}{\delta }_{\nu {\nu }^{\prime \prime }}+\dots .$$

$$⟨\mathrm{Tr}({\mathcal{M}}^{2k})⟩=C_k{N}^{k+1}+\dots .$$

$$⟨\mathrm{Tr}({\mathcal{H}}^{n_1})\times \dots \times \mathrm{Tr}({\mathcal{H}}^{n_k})⟩=⟨\mathrm{Tr}({\mathcal{M}}^{n_1})\times \dots \times \mathrm{Tr}({\mathcal{M}}^{n_k})⟩.$$

$$\rho (E)=-\frac{1}{N\pi }\underset{\epsilon \to 0}{\lim }\Im \mathrm{Tr}[G^{(r)}(E)].$$

$$\mathrm{Tr}[G^{(r)}(E)]=\sum _{n=0}^{\infty }\frac{1}{(E^{+}{)}^{n+1}}\mathrm{Tr}[{\mathcal{M}}^n].$$

$$⟨G^{(r,N)}(E)⟩=\frac{1}{E}+\frac{1}{E}⟨\mathcal{M}(⟨G^{(r,N)}(E)⟩)\mathcal{M}⟩⟨G^{(r,N)}(E)⟩.$$

$$⟨G^{(r,N)}(E)⟩=\frac{1}{E}+\frac{1}{E}\mathrm{Tr}(⟨G^{(r,N)}(E)⟩)⟨G^{(r,N)}(E)⟩.$$

$$\mathrm{Tr}(⟨G^{(r,N)}(E)⟩)=\frac{E}{2}-i\sqrt{N}\sqrt{1-\frac{E^2}{4N}}.$$

$$\langle G^{(r,N)}(E)\rangle_{\mu\nu}=\frac{1}{N}\bigg{(}\frac{E}{2}-i\sqrt{N}\sqrt{1-\frac{E^{2}}{4N}}\bigg{)}\delta_{\mu\nu}\ .$$

$$\delta G_{\mu \nu }^{(r,bin)}(E)=-\frac{1}{E}⟨G^{(r,N)}(E){⟩}_{\mu \mu }⟨G^{(r,N)}(E){⟩}_{\mu \nu }.$$

$$\delta G_{\mu \nu }^{(r,four)}(E)$$

$$⟨\mathrm{Tr}{G}^{(r)}(E+{\epsilon }_1)\times \dots \times \mathrm{Tr}{G}^{(r)}(E+{\epsilon }_P)$$

$$\times \mathrm{Tr}{G}^{(a)}(E-{\tilde{\epsilon }}_1)\times \dots \times \mathrm{Tr}{G}^{(a)}(E-{\tilde{\epsilon }}_Q)⟩.$$

$${\rm Tr}\frac{1}{E^{+}-{\cal M}}=\frac{1}{2}\frac{\partial}{\partial j}{\cal G}(j)\bigg{|}_{j=0}\ {\rm where}\ {\cal G}(j)=\frac{\det(E^{+}-{\cal M}+j)}{\det(E^{+}-{\cal M}-j)}\ .$$

$$\psi =(s_1,\dots ,s_N,{\chi }_1,\dots ,{\chi }_N{)}^T$$

$$\mathcal{C}=(E^{+}{1}^N-\mathcal{M})1^s-j{\sigma }_3{1}^N.$$

$$\mathcal{G}(j)=\int \mathrm{d}({\psi }^{\ast },\psi )\exp \{(i/2)\tilde{\psi }\mathcal{C}\psi \}.$$

$$⟨(\tilde{\psi }\mathcal{M}{1}^s\psi {)}^2⟩$$

$$⟨(\tilde{\psi }\mathcal{M}{1}^s\psi {)}^4⟩$$

$$\displaystyle\bigg{\langle}\exp\{-(i/2)(\tilde{\psi}{\cal M}1^{s}\psi)\}\bigg{\rangle}$$

Full text

Spectral statistics of the uni-modular ensemble Christopher H. Joyner1, Uzy Smilansky2 and Hans A. Weidenmüller3

1 Introduction

In this work we investigate the spectral statistics of the uni-modular ensemble of random matrices (UME) in the limit of large matrix dimension $N$. The ensemble is defined as the set of Hermitian matrices ${\cal M}=\{M_{\mu\nu}\}$, where $\mu,\nu=1,2,\ldots,N$, with elements of the form

[TABLE]

and where, except for the symmetry relation, the phases $\phi_{\mu\nu}$ are uncorrelated real random variables distributed uniformly in the interval $[0,2\pi)$. In contrast, the Gaussian unitary ensemble (GUE) is given by the set of Hermitian matrices $\mathcal{H}=\{H_{\mu\nu}\}$ of dimension $N$ endowed with the probability distribution

[TABLE]

where $Z_{N}$ is a normalization factor. The UME serves as a paradigmatic example of a Wigner ensemble, i.e. a set of random Hermitian matrices with independently distributed elements that do not follow a Gaussian distribution (see e.g. [1] for details). Wigner [2, 3] was the first to show that many spectral properties of Wigner ensembles coincide with those of the GUE in the limit of large $N$. Little is known, however, about the $1/N$ deviations of particular Wigner ensembles from the universal limit.

We address this question for the UME. We have chosen that ensemble because, as we shall show, there exist exact relations that allow us to determine the spectral properties of the UME even though the ensemble is not unitarily invariant. That allows us to compute the $1/N$ corrections to the universal limit. To our knowledge there has been little previous work on uni-modular ensembles, except for the work by Sodin and Feldheim, who investigate the fluctuations of spectral moments in the UME and the minimum eigenvalue of unimodular covariance matrices [4, 5, 6], and that of Lakshminarayan, Puchala and Zyczkowski [7], who obtain exact expressions for the first four spectral moments of unimodular covariance matrices111In contrast to the present article, the term ‘unimodular ensemble’ is used in [7] for the non-Hermitian counterpart of (1). and provide a conjecture for all moments.

The article is structured as follows: In the remainder of the present Section we briefly discuss the distribution of the spectral moments in the context of the GUE and related ensembles. We then use a general argument to show that in the large $N$ limit the probability distribution of these moments in the UME approaches that of the GUE and discuss the difficulties in proceeding with corrections. In Section 2 we obtain corrections to the mean spectral density via a supersymmetric approach and point to the difficulties of going to higher order. In Section 3 we show how these difficulties may be addressed in some cases using a graph-theoretical and combinatorial techniques. Finally, in Section 4 we use a Brownian motion approach, similar in spirit to Dyson’s original approach [8], with a theorem of Meckes [9] utilising Stein’s method, to show that the fluctuations of the spectral moments are Gaussian in the large $N$ limit, and provide rates of convergence.

1.1 Spectral moments in random-matrix theory

The spectral moments $\tau_{k}(\mathcal{H})$ of a Hermitian matrix $\mathcal{H}$ are given by the moments of the empirical spectral density,

[TABLE]

Here $\lambda_{i}$ with $i=1,2,\ldots,N$ are the eigenvalues of $\mathcal{H}$. For Wigner matrices, the mean values of these moments vanish for odd $k$ and, for even $k=2\nu$, converge in the limit of large matrix dimension to $C_{\nu}N^{k}$ where $C_{\nu}=(2\nu)!/(\nu!(\nu+1)!)$ are the Catalan numbers [2, 3]. In addition, the variances converge sufficiently rapidly to conclude that the density converges weakly, almost surely, to Wigner’s semi-circle law $\sigma(\lambda)=\frac{1}{(2\pi)}\sqrt{4-\lambda^{2}}$ (see e.g. [1] for instance). For the GUE, with the average defined by (2) and the mean values of the moments (indicated by angular brackets) by

[TABLE]

Harer and Zagier [10] discovered the three-term recurrence relation

[TABLE]

That recurrence relation immediately leads to the following correction to Wigner’s leading term,

[TABLE]

Recurrence relations similar to (5) were found for the GOE and GSE by Ledoux [11] and for ensembles characterized by the index $\beta$ by Witte and Forrester [12]. For Gaussian, Laguerre and Jacobi $\beta$ ensembles, exact expressions for the moments were given by Mezzadri, Reynolds and Winn in terms of Jack polynomials [13]. These, however, do not seem to lend themselves to asymptotic expansions in $1/N$. The systematic approach to the $1/N$ expansion has also been addressed in the context of RMT distributions of the mean delay time. These involve random matrix ensembles with exact expressions for the joint probability density functions of the eigenvalues (see e.g. [14, 15, 16] and references therein). We are not aware of attempts to go beyond leading order in other matrix ensembles.

Fluctuations of the moments (for an arbitrary ensemble defined in analogy to (3)) are often discussed in terms of the so-called linear-statistic

[TABLE]

where ${\rm Tr}[f(\mathcal{H})]:=\sum_{i=1}^{N}f(\lambda_{i}(\mathcal{H}))$. We note that if $f$ is a polynomial then ${\rm Tr}[f(\mathcal{H})]$ is simply a weighted sum over the moments. The distribution of $L_{f}(\mathcal{H})$ and related quantities were first analysed by Jonson in the case of Wishart matrices [17], by Johansson in the case of unitarily invariant matrices [18], and by a number of authors in the case of Wigner matrices [19, 20]. In all cases one observes convergence to a Gaussian distribution, with a universal variance, in the limit of large matrix size.

There exists a large number of papers - too many for detailed referencing - that prove universality of $L_{f}(\mathcal{H})$ for various types of random matrices. We only emphasise some results that are particularly relevant to this article. Using similar techniques to those presented here, Sodin has shown for the UME that the moments of $L_{f}({\cal M})$ are Gaussian (see e.g. [6] and references therein) but does not discuss rates of convergence. Chatterjee has previously used Stein’s method along with estimations of Poincaré inequalities to provide bounds on the total-variation distance between a Gaussian and $L_{f}(\mathcal{H})$ for appropriate random matrix ensembles [21]. Finally Cabanal-Duvillard has used a Brownian-motion approach to derive similar results for the GUE [22], using the eigenvalue motion directly.

1.2 Moments of the UME

We show that for $N\to\infty$, the distribution of the spectral moments of the UME coincides with that of the GUE. We do so by showing that in the limit, all moments and all products of moments of the UME have the same values as for the GUE. The latter, defined in (2), consists of Hermitean matrices ${\cal H}$ with elements $H_{\mu\nu}=H^{*}_{\nu\mu}$ that are Gaussian-distributed zero-centred random variables with second moments

[TABLE]

The normalization of the matrix elements in (8) implies that the support of the spectral density is $(-2\sqrt{N},+2\sqrt{N})$. The elements of the UME have zero average and second moments

[TABLE]

We note that for $\mu\neq\nu$, $|M_{\mu\nu}|^{2}=1$ without averaging.

We first show that to leading order in $1/N$ we have

[TABLE]

The Gaussian distribution of $H_{\mu\nu}$ implies $\langle{\rm Tr}({\cal H}^{n})\rangle=0$ for odd values of $n$. For $n=2k$ even, the trace is calculated using Wick contraction. Contributions of leading order in $1/N$ arise only from a subset of all Wick contraction patterns (“nested contractions”) where contraction lines connecting pairs of contracted matrix elements do not intersect. The result is

[TABLE]

The Catalan numbers count the number of nested contractions in the $n^{th}$ moment, with $n=2k$. The dots indicate terms of order $N^{l}$ with $l\leq k$.

For the UME, we obviously have $\langle{\rm Tr}({\cal M}^{n})\rangle=0$ for $n$ odd. For ${\rm Tr}({\cal M}^{2})$, (9) yields $N^{2}-N$. That differs from the GUE result ${\rm Tr}({\cal H}^{2})=N^{2}$ by a term of order $N^{-1}$. That term is due to the last Kronecker symbol in (9). Higher even moments of ${\cal M}$ receive contributions not only from the pairwise correlations displayed in (9), but also from correlations of order $4,6,...$. We demonstrate the existence of such correlations for the case of order four. The average of the term $M_{\mu\nu}M_{\mu^{\prime}\nu^{\prime}}M_{\mu^{\prime\prime}\nu^{\prime\prime}}M_{\mu^{\prime\prime\prime}\nu^{\prime\prime\prime}}$ vanishes unless the indices are pairwise equal but differ within each factor ${\cal M}$. For the correlation of order four (as opposed to a product of correlations of order two) that gives

[TABLE]

The factors in straight brackets impose the conditions $\phi_{\mu\nu}=-\phi_{\mu^{\prime}\nu^{\prime}}$ and $\phi_{\mu^{\prime\prime}\nu^{\prime\prime}}=-\phi_{\mu^{\prime\prime\prime}\nu^{\prime\prime\prime}}$. The last two Kronecker deltas yield $\phi_{\mu\nu}=\phi_{\mu^{\prime\prime}\nu^{\prime\prime}}$, a condition that would be absent for a product of two correlations of order two. The dots indicate terms obtained by a permutation of the indices. Direct calculation yields ${\rm Tr}({\cal M}^{4})=2N^{3}-3N^{2}+N$. That differs from the GUE result ${\rm Tr}({\cal H}^{4})=2N^{3}+N$ by terms of order $1/N$. The difference is due to the last Kronecker delta in (9) and to the two Kronecker deltas in (12). Each Kronecker delta reduces the number of independent summations in the expression for the trace and, thus, produces terms of order $1/N$. Correlations of higher order than four exist and carry additional Kronecker symbols beyond the ones in (12). Therefore, with increasing $n$ the expressions for ${\rm Tr}({\cal M}^{2n})$ and for ${\rm Tr}({\cal H}^{2n})$ become ever more different. The differences are confined, however, to terms of order $1/N$ or smaller. The term of leading order in ${\rm Tr}({\cal M}^{2n})$ is obtained by taking account only of binary correlations and by omitting the last Kronecker delta in (9). Hence,

[TABLE]

The terms indicated by dots are of lower order in $N$. They differ from the terms indicated in the same manner in (11). Comparison of (11) and (13) shows that all moments of the GUE and of the UME become asymptotically ($N\to\infty$) equal, and that (10) holds.

To see in which sense (10) applies we consider next-order corrections. We start with corrections due to last Kronecker delta in (9). In (13), terms of order $N^{k}$ arise from all nested contributions involving that term once. The additional Kronecker delta can be affixed to each one of the pairwise contractions ($k$ pairs). The number of contributions is, therefore, equal to $kC_{k}$ and the total contribution is given by $kC_{k}N^{k}$. In comparison with the result (13) that term is of order $k/N$. For fixed $k$ the contribution vanishes for $N\to\infty$. It does not vanish, however, for fixed $N$ and $k\to\infty$. An analogous conclusion holds for the contribution of correlations of higher order to the right-hand side of (13). We conclude that (10) is an asymptotic relation. It establishes the identity of the $k^{\rm th}$ moments for fixed $k$ and $N\to\infty$.

We turn to the average of products of moments and show that for all positive integer $k,n_{1},n_{2},\ldots,n_{k}$ and to leading order in $N$ we have

[TABLE]

The left-hand side of (14) is evaluated by calculating all Wick contractions of pairs of matrix elements of ${\cal H}$. That rule comprises pairs of matrix elements occurring under the same trace and pairs that occur in different traces. Only nested contributions contribute to the leading order in $N$. The right-hand side of (14) is evaluated using the binary correlation of (9) as well as all higher-order correlations as exemplified in (12). These likewise comprise sets of matrix elements that occur either under the same trace or under two or more different traces. As in the case of (10) we use the fact that the leading-order terms in $N$ are obtained by suppressing the minimum number of independent summations over matrix indices. That rules out all higher-order correlations and leaves us with the binary correlations of (9). For the terms of leading order in $N$ we suppress the last Kronecker delta in that equation. As a result we find that the leading-order terms in $N$ of the right-hand side of (14) are obtained by calculating all Wick contractions of matrix elements of ${\cal M}$ (occurring either under the same trace or as arguments of different traces). For each Wick-contracted pair the rule is the same as for the GUE in (8). The rules for calculating the right-hand side of (14) being the same as for the left-hand side, the results are the same, too, and (14) is seen to hold to leading order in $N$. Again, that is an asymptotic result. It holds for fixed $n_{1},n_{2},\ldots,n_{k}$ and $N\to\infty$.

To leading order in $1/N$, these results imply the equality of the mean spectral density of the GUE and the UME and also the convergence in distribution of the $\tau_{k}$ and by extension of the $L_{f}({\cal M})$ for polynomial functions $f$. They do not, however, allow us to obtain corrections to this density or rates of convergence for the distributions. These aspects are explored in subsequent sections.

1.3 Mean spectral density

The empirical density $\rho(E)=\frac{1}{N}\sum_{i=1}^{N}\delta(E-\lambda_{i}(\mathcal{H}))$ (see Eqn. (3)), normalised so that $\int dE\rho(E)=1$, can be written in terms of the retarded Green’s function $G^{(r)}(E)=(E^{+}-{\cal M})^{-1}$, where $E^{+}=E+i\varepsilon$ with $\varepsilon$ infinitesimal and positive, as

[TABLE]

We expand the retarded Green’s function for the UME as

[TABLE]

We first calculate $\langle G^{(r)}(E)\rangle$ to leading order in $1/N$ and then consider the sub-leading contributions. Using the expansion (16) and taking into account only nested contributions to the average, we obtain the Pastur equation [23]

[TABLE]

The upper index $N$ stands for the leading-order term. We use the binary correlator (9), suppress the last Kronecker delta, and obtain

[TABLE]

We take the trace of (18) and solve the resulting quadratic equation for ${\rm Tr}(\langle G^{(r,N)}(E)\rangle)$. That gives

[TABLE]

The range of the spectral density is $(-2\sqrt{N},+2\sqrt{N})$. For the full Green function we find

[TABLE]

To leading order the spectral density is the same as for the GUE, as expected.

Correction terms of order $1/N$ to (18) arise when either the last Kronecker delta in (9) or the fourfold correlation (12) are taken into account once. Non-nested contributions and higher-order correlations do not contribute to that order. Using in (17) the last Kronecker delta in (9) we obtain

[TABLE]

The additional contribution in (17) due to the fourfold correlation term (12) is

[TABLE]

Equation (20) shows that $G_{\mu\mu}(E)$ is of order $1/N$. Therefore, both contributions (21) and (22) are of order $1/N$ compared to the leading contribution in (18). Adding the results (21) and (22) we obtain as a $(1/N)$-correction to the spectral density of the UME a polynomial of fourth order in $\langle G^{(r,N)}(E)\rangle$. That correction is completely different from the $1/N$ oscillations of the spectral density displayed in later sections of the paper. The reason is that the Pastur equation is valid only asymptotically. It is derived with the help of the same asymptotic expansion as used for the averaged traces in Eq. (10). That approach cannot be used for a systematic evaluation of terms of next order in $1/N$.

In Section 2 and Section 3 we go beyond leading order by using a supersymmetry approach first developed in [28, 29] and a graph theoretic approach adapted from $d$-regular graphs [30, 31].

2 Supersymmetry

Equation (14) suggests that all level correlation functions for the UME coincide to leading order in $1/N$ with those of the GUE. The argument goes as follows. The $(P,Q)$ level correlation function for the UME is defined as

[TABLE]

Here $G^{(r)}(E)$ and $G^{(a)}(E)$ are the retarded and the advanced Green functions for the UMA, respectively. The increments $\varepsilon_{p}$, $p=1,\ldots,P$ and $\tilde{\varepsilon}_{q}$, $q=1,\ldots,Q$ are of the order of the mean level spacing. The $(P,Q)$ level correlation function for the GUE has the same form except for the replacement ${\cal M}\to{\cal H}$ in each of the Green’s functions.

We use the expansion (16) for ${\rm Tr}G^{(r)}(E)$ and proceed correspondingly for ${\rm Tr}G^{(a)}(E)$. Each term in the resulting expansion of the correlation function (23) contains an ensemble average over products of traces of powers of ${\cal M}$ that has the form of the right-hand side of (14). We proceed analogously for the level correlator of the GUE, expanding the Green’s functions in powers of ${\cal H}$. Each term in the resulting series is obtained from the corresponding term of the UME by the formal replacement ${\cal M}\to{\cal H}$. That same replacement converts the ensemble average over products of traces of powers of ${\cal M}$ into the ensemble average over products of traces of powers of ${\cal H}$. With (14) showing that these averages are equal to leading order in $N$ we conclude that all $(P,Q)$ level correlationfunctions of the UME coincide with those of the GUE in that order.

The argument lacks stringency, however. It is based upon a perturbative expansion of the correlation functions. In contrast to the spectral density, all correlation functions possess a zero mode. The two-point function, for instance, has a zero mode at $\varepsilon_{1}=0=\tilde{\varepsilon}_{1}$ and thus, perturbatively, a singularity. That is why we turn to the supersymmetry approach where the zero mode is treated exactly.

The one-point function is written as

[TABLE]

The generating function ${\cal G}(j)$ is written as a superintegral. The $2N$-dimensional supervector

[TABLE]

contains the commuting complex variables $s_{k}$ and the anticommuting variabless $\chi_{k}$, $k=1,\ldots,N$ with $\int\chi_{k}{\rm d}\chi_{k}=(2\pi)^{-1/2}=\int\chi^{*}_{k}{\rm d}\chi^{*}_{k}$ for all $k$. The integration measure is the flat Berezinian ${\rm d}(\psi^{*},\psi)=\prod_{k}{\rm d}\Re(s_{k}){\rm d}\Im(s_{k}){\rm d}\chi^{*}_{k}{\rm d}\chi_{k}$. In the $2N$-dimensional superspace (the direct product of the $N$-dimensional ordinary space with indices $k=1,\ldots,N$ and the two-dimensional superspace with indices $s=1,2$) we define

[TABLE]

Here $1^{s}$ and $\sigma_{3}$ are the unit matrix and the third Pauli spin matrix, respecively, in two-dimensional superspace while $1^{N}$ is the unit matrix in ordinary space. With $\tilde{\psi}=(\psi^{*})^{T}$ we have

[TABLE]

The ensemble average of ${\cal G}$ is defined as an average over the independent phases $\phi_{kl}$ with $k<l$.

To average ${\cal G}$ we calculate the expectation value of $\exp\{-(i/2)(\tilde{\psi}{\cal M}1^{s}\psi)\}$. We first consider the moments of $(\tilde{\psi}{\cal M}1^{s}\psi)$. All odd moments vanish. For the second moment we use (9). The fourth moment is the sum of the binary and the fourfold correlations given in (9) and (12), respectively. Thus,

[TABLE]

That gives

[TABLE]

The first two terms in the exponent strongly suggest how the series continues although we have not checked that. Writing

[TABLE]

we observe that the second term on the right-hand side of (29) is of order $1/N$ compared to the first one. The same is true of the second term on the right-hand side of the second of (27) in comparison with the first one. And the same statement holds a fortiori for higher-order correlations. To leading order in $1/N$ we, therefore, have

[TABLE]

For the GUE we have correspondingly

[TABLE]

The equality of the right-hand sides of (30) and (31) implies that to leading order in $1/N$ the spectral densities of the UME and of the GUE are identical.

We turn to the $(P,Q)$ level correlation function of the UME as defined in (23). We skip the construction of the generating function for the correlation function (23) because it runs completely parallel to that for the GUE given in Ref. [24]. Suffice it to say that the result is similar in form to (26), with the following differences. The vectors $\psi$ and $\tilde{\psi}$ and the matrix ${\cal C}$ now have dimension $2N(P+Q)$, the matrix ${\cal C}$ contains the matrix ${\cal M}$ in block-diagonal form $(P+Q)$ times, the scalar $j$ becomes a matrix of dimension $(P+Q)$, the vector $\psi$ ($\tilde{\psi}$) is multiplied from the left (right) by the matrix ${\bf L}^{1/2}$ where ${\bf L}={\bf 1}$ in the retarded and ${\bf L}=-{\bf 1}$ in the advanced sector, and the energy increments $\varepsilon_{1},\ldots,\varepsilon_{P}$ and $\tilde{\varepsilon}_{1},\ldots,\tilde{\varepsilon}_{Q}$ appear in the exponent. Evaluating the expectation value of $\exp\{-(i/2)\tilde{\psi}{\cal M}\psi\}$ by using (27) and dropping terms of higher order in $1/N$ we obtain exactly the same form for the averaged generating function as for the GUE. That implies that to leading order in $1/N$, all $(P,Q)$ level correlation functions for the UME coincide with those of the GUE.

There are two possibilities to go beyond these results. First, the average over the phases $\phi_{\mu\nu}$ can be done exactly using the color-flavor transformation [25]. For quantum graphs, that transformation was used in Refs. [26, 27, 24]. In the present context the color-flavor transformation seems uneconomical, however. The reason is seen by considering the generating function (26) for the spectral density. In (26) the ensemble average has to be taken by integrating over the real phase angles $\phi_{\mu\nu}$ with $\mu<\nu$ and $\mu,\nu=1,\ldots,N$. The color-flavor transformation performs these integrations at the expense of introducing for every $\phi_{\mu\nu}$ a pair of supermatrices $(Z_{\mu\nu},\tilde{Z}_{\nu\mu})$ with $\mu<\nu$. That increases the number of integration variables by a factor eight. We have, therefore, not followed that possibility. The second, more attractive possibility is to examine correction terms of order $1/N$ within the Hubbard-Stratonovich approximation to the supersymmetry approach. We do that, confining ourselves to the spectral density.

2.1 Terms of Order $1/N$

The papers by Kalisch and Braak [28] and by Shamis [29] show how corrections of order $1/N$ to the spectral density of the GUE can be worked out. Shamis makes heavy use of the unitary invariance of the GOE. Such invariance is not shared by the UME. That is why we follow the work of Kalisch and Braak. These authors use the Hubbard-Stratonovich transformation, perform the integration over the two remaining anticommuting variables exactly, and then use the saddle-point approximation. We apply that method in our more general context.

We approximate the ensemble average of $\exp\{-(i/2)\tilde{\psi}{\cal M}\psi\}$ by keeping only correction terms of order $1/N$. That is done by dropping in (28) the last term and by writing the second term as in (29). Thus,

[TABLE]

We eliminate the first quartic term by a single Hubbard-Stratonovich transformation and the second one by $N$ such transformations, one for each term $(\tilde{\psi}_{k}\psi_{k})^{2}$. We define the two-by-two supermatrices

[TABLE]

so that

[TABLE]

We use

[TABLE]

The supermatrices $\sigma$ and $\tau_{k}$ all have dimension two. We insert that into (33) and carry out the Gaussian integrals over the original integration variables. That gives

[TABLE]

The indices indicate that the supertraces extend only over the superindices $s$. We remark in parentheses that (36) shows the limitations of the supersymmetry approach. For the GUE we deal with a single supermatrix $\sigma$. The $1/N$ correction to the GUE introduces $N$ additional supermatrices $\tau_{k}$. Corrections of higher order in $1/N$ lead to ever more complex integrals, bringing the method to its limit.

2.1.1 GUE

We are guided by the GUE case. We consider only terms up to first order in $j$ and indicate that fact by using the sign $\approx$ in Eqs. (37) through (40),

[TABLE]

We write

[TABLE]

Here $s_{B},s_{F}$ are real commuting and $\alpha,\alpha^{*}$ anticommuting variables. We define $a=E^{+}-s_{B}$, $b=E-is_{F}$ and carry out the integrals over the anticommuting variables. Then

[TABLE]

We rescale $E\to x=E/\sqrt{N}$ and with it $s_{B}\to q=s_{B}/\sqrt{N}$, $s_{F}\to p=s_{F}/\sqrt{N}$, $a\to c=a/\sqrt{N}=x^{+}-q$, $b\to d=b/\sqrt{N}=x-ip$, and $j\to j^{\prime}=j\sqrt{N}$. The last operation assures that the average level density is normalized to unity. Choosing $j^{\prime}=j/\sqrt{N}$ would yield an average level density normalized to $N$. Then

[TABLE]

Using Eq. (24) we have performed the integrals over $p$ and $q$ analytically. The resulting expression in Hermite polynomials agrees for every $N$ with the standard result.

For the saddle-point approximation we define the effective action

[TABLE]

It is the sum of the effective actions ${\cal A}_{q}$ and ${\cal A}_{p}$ for the variables $q$ and $p$. Therefore, the saddle points for $p$ and for $q$ are unrelated. The saddle points for $q$ are $q_{0}=(x/2)\pm i\sqrt{1-x^{2}/4}$. Because of the singularity of the integrand at $x^{+}$ we admit only the the solution in the lower half plane so that

[TABLE]

As $x^{+}$ moves from $-2^{+}$ to $+2^{+}$ the saddle point $q_{0}$ moves from $-1$ on a semicircle into the lower half plane, reaches the value $-i$ for $x=0$, and continues on the semicircle to $+1$ for $x=2$. Without crossing the singularity, we can shift the path of integration for all values of $x$ with $|x|\leq 2$ so that it runs parallel to the real axis and passes through $q_{0}$. We write $q=q_{0}+t_{B}$ and expand ${\cal A}_{q}$ around $q_{0}$ in powers of $t_{B}$ up to second order,

[TABLE]

The saddle points for $p$ are

[TABLE]

As $x$ increases from $-2$ to zero, the two saddle points (that are degenerate at $x=-2$ with value $+i$) move on semicircles in opposite directions, reaching the real axis at $-1$ and at $+1$ for $x=0$ and continue into the lower half plane, reaching the degenerate point $-i$ at $x=+2$. For all values of $|x|\leq 2$ the two saddles lie on a straight line parallel to the real axis. We shift the path of integration so that it runs along that line. For ${\cal A}_{p}$ we write $p=p_{0}+t_{F}$ and expand ${\cal A}_{p}$ in powers of $t_{F}$ up to second order,

[TABLE]

We evaluate $\langle{\cal G}(j^{\prime})\rangle$ at the two saddle points defined by $ip_{0}=q_{0}$ and $ip_{0}=q^{*}_{0}$. We do so by expanding up to and including terms of order $1/N$. Upon integration over $t_{B}$ and $t_{F}$ that gives for $ip_{0}=q_{0}$ the result $\langle{\cal G}(j^{\prime})^{-}\rangle=1+2j^{\prime}q_{0}$. Equation (24) and the fact that the imaginary part of the retarded Green function equals $-\pi\delta(E-{\cal H})$ then yields for the spectral density

[TABLE]

That is the asymptotic expression ($N\to\infty$). For the first saddle point, it holds up to and including terms of order $1/N$.

For $ip_{0}=q^{*}_{0}$ the leading-order contribution is of order $1/N$ and given by

[TABLE]

The $(1/N)$-correction to the level density (46) is [28, 29]

[TABLE]

Characteristic features are the rapid oscillations with frequency $1/N$ and the singularities at the end points $x=\pm 2$ of the spectrum. For finite $N$ the exact expression for the spectral density is non-singular for all values of the energy. The singularities occur only in the $1/N$ expansion.

2.2 UME

We start from (36). In the calculation of the normalization integral ${\cal G}(0)$, we proceed as in Section 2.1.1, see also the calculation of the source terms discussed below. The effective action ${\cal A}$ is defined as the contribution of leading order in $1/N$ to the negative exponent of the integrand in ${\cal G}(0)$. We mention without proof that the effective action turns out to be equal to the effective action for the GUE in Eq. (41). The saddle points are the same. When we calculate the leading-order contribution of the two saddle points to ${\cal G}(0)$ we find that the result is identical to the GUE expression in (40). Therefore, ${\cal G}^{-}(0)=1$ and ${\cal G}^{+}(0)=0$.

We turn to the source terms. For $\sigma$ we use the parametrization (38). For $\tau_{k}$ with $k=1,2,\ldots,N$ we write

[TABLE]

Here $t_{kB},t_{kF}$ are real commuting and $\gamma_{k},\gamma^{*}_{k}$ are anticommuting variables. Then

[TABLE]

We define

[TABLE]

In (36) we keep terms of first oder in $j$. We have

[TABLE]

We use the identity

[TABLE]

which is valid for an arbitrary supermatrix of dimension two. The prefactor in expression (52) becomes

[TABLE]

A change of integration variables shows that that factor is the same for every value of $k$. Therefore

[TABLE]

We integrate explicitly over all anticommuting variables. That gives

[TABLE]

Integration over the commuting variables $\tau_{k,B}$ and $\tau_{kF}$ shows that the leading terms in powers of $1/N$ are

[TABLE]

We rescale $E\to x=E/\sqrt{N}$ and with it $s_{B}\to q=s_{B}/\sqrt{N}$, $s_{F}\to p=s_{F}/\sqrt{N}$, and $j\to j^{\prime}=j\sqrt{N}$. Then

[TABLE]

That expression agrees with the source terms in Eq. (40) showing that the spectral density and its oscillations are in leading order the same for the UE and for the GUE.

We have reported in Section 2.1.1 that the supersymmetry approach, when evaluated exactly, yields the correct finite-$N$ expression for the spectral density of the GUE in terms of Hermite polynomials. The approximation used in (32) for the UME defines a new random-matrix ensemble. That ensemble is identical to the GUE except that the diagonal elements vanish. Therefore, we expect that exact expressions for the spectral density can be obtained also in this case from the supersymmetry approach.

Starting from the approximate form of the action given in Eq. (32), we have performed the integrals as done in the GOE case. We have carried the $1/N$ expansion beyond the leading-order terms. We have used the result to construct an expansion of the spectral density of the UME in terms of Hermite polynomials as was done for the GUE. We have failed to obtain meaningful results. We ascribe that to the fact that for technical reasons we have not been able to include all terms in the $1/N$ expansion. For a truncated expansion, the spectral density is expected to be singular. We have not succeeded in separating these singularities from the oscillatory behaviour of the spectral density.

3 Graph-theoretical approach

3.1 Background and definitions

We start with a few definitions of graph-theoretical concepts which are helpful for the subsequent discussion.

In a simple graph with $N$ vertices, any two different vertices are connected by at most one edge; no edge begins and ends in the same vertex. In a complete graph, every two different vertices are connected by an edge. The elements $A_{i,j}=A_{j,i}$ of the symmetric vertex adjacency matrix $A$ of dimension $N$ equal $1$ ([math]) if the vertices $i$ and $j$ are connected (not connected, respectively). For a simple graph, $A_{i,i}=0$. A directed edge $e=(j,i)$ connects the vertices $j,i$ and has direction $i\rightarrow j$. The vertex $i$ is the origin of $e$ and $j$ is its terminus: $i=o(e),\ j=\tau(e)$. The direction of the edge $\hat{e}$ is opposite to that of $e$. The number of directed edges equals $\sum_{i,j}A_{i,j}$. The matrix ${\cal B}$ with elements $B_{e^{\prime},e}=\delta_{o(e^{\prime}),\tau(e)}$ describes the way the vertices are connected in the space of directed edges. For a complete graph, the matrix ${\cal B}$ has dimension $N(N-1)$.

A walk of length $t$ on the graph is defined by a list of directed edges $e_{1},e_{2},\cdots,e_{t}$ where $B_{e_{i},e_{i+1}}\neq 0$. For a $t$-periodic walk $B_{e_{t},e_{1}}\neq 0$. A cycle is the set of all periodic walks that differ only by a cyclic permutation of their edges. A cycle is primitive if the edge list is not a repetition of a shorter list. Writing $J_{e^{\prime},e}=\delta_{e^{\prime},{\hat{e}}}$ we define the Hashimoto matrix ${\cal Y}={\cal B}-J$ which connects only directed edges that are not reversed to each other. Thus, while ${\rm Tr}({\cal B}^{t})$ counts the number of $t$-periodic walks on the graph, ${\rm Tr}({\cal Y}^{t})$ counts the number of $t$-periodic walks where back-tracking is not allowed.

To use these definitions for the UME we write the phases $\phi_{\mu\nu}$ of (1) as $\phi_{e}$ with $e=(\mu,\nu)$. Following [30, 31], we include these phases in the definitions of the matrices ${\cal B}$ and ${\cal Y}$. Denoting the set of phases of the matrix ${\cal M}$ by $\Phi$ we define magnetic edge connectivity matrix ${\cal B}(\Phi)$ and the magnetic Hashimoto matrix ${\cal Y}(\Phi)$ as

[TABLE]

The term “magnetic” relates the phases to a (fictitious) magnetic field. Contributions to ${\rm Tr}[{\cal Y}(\Phi)^{n}]$ arise from the set $\Omega_{n}$ of all $n$-periodic non-backtracking walks on the graph. In the magnetic case each walk $\omega$ contributes a phase $\Phi_{w}=\sum_{e\in w}\phi_{e}$ so that

[TABLE]

An important identity due to Bass, generalized by Bartholdi and extended to the magnetic case in Ref.[30], connects the spectra of the UME matrix ${\cal M}$ and of ${\cal Y}(\Phi)$. The Bass identity for complete magnetic graphs is valid for any $\eta\in\mathbb{C}$. With $I^{(k)}$ the identity matrix of dimension $k$, it reads

[TABLE]

The identity shows that but for a factor $\eta^{N}(\eta^{2}-1)^{\frac{N(N-3)}{2}}$, the characteristic polynomial of ${\cal Y}(\Phi)$ is proportional to the characteristic polynomial of ${\cal M}$ evaluated at $\frac{\eta^{2}+(N-2)}{\eta}$. The spectra of the two matrices are, therefore, related. Let $\sigma({\cal M})\doteq\{\lambda_{k}\}_{k=1}^{N}$ ($\sigma({\cal Y}(\Phi))\doteq\{\eta_{r}\}_{r=1}^{N(N-1)}$) be the spectrum of ${\cal M}$ (of ${\cal Y}(\Phi)$, respectively). The factor on the right-hand side of (61) vanishes at $\eta=\pm 1$ with multiplicity $N(N-3)/2$ and at the $2N$ eigenvalues of ($I^{(N)}(\eta^{2}+(N-2))-\eta M$). These can be expressed in terms of the $\lambda_{k}$,

[TABLE]

From the left-hand side of (61) we see that the nontrivial part of the spectrum of ${\cal Y}(\Phi)$ consists of two sets of points. The first set is confined to the circle of radius $\sqrt{N-2}$ in the complex plane. It corresponds to the spectral points of ${\cal M}$ that lie in the interval $[-2\sqrt{N-2},2\sqrt{N-2}]$. We write $\sigma^{R}({\cal M})=\{\lambda_{k}:|\lambda_{k}|\leq 2\sqrt{N-2}\}$. The second set consists of real pairs $(\eta_{+},\eta_{-})$ whose product is $(N-2)$. These correspond to the spectral points in $\sigma^{NR}({\cal M})=\sigma({\cal M})-\sigma^{R}({\cal M})$. Matrices ${\cal M}$ for which the entire spectrum belongs to $\sigma^{R}({\cal M})$ are referred to as “Ramanujan” matrices - a term which we borrow freely from an analogous situation in the spectra of $d$-regular graphs. Conversely, matrices for which at least one of the spectral points does not lie in $\sigma^{R}({\cal M})$ are called “non-Ramanujan”.

For convenience we scale the UME matrices in the following manner;

[TABLE]

so that $\epsilon_{k}=\lambda_{k}/(2\sqrt{N-2})$ are the eigenvalues222We note that with this scaling, the semi-circle density (the limit of the mean spectral density for large $N$) is defined on the interval $[-1,1]$, in contrast to the previous sections where different scaling led to the interval $[-2,2]$. Using the connection of the two spectra given by the Bass identity (61), we can write the normalized trace of ${\cal Y}(\Phi)^{n}$ as

[TABLE]

Here $T_{n}(x)$ are the Chebyshev polynomials of the first kind, given by

[TABLE]

Since we sum over eigenvalues in both $\sigma^{\rm R}$ and $\sigma^{\rm NR}$, equation (66) is valid for every matrix ${\cal M}$ in the UME, independently of whether it is Ramanujan or not. Equation (66) has been derived in an alternative manner by Sodin (see e.g. Section 4.2.6 in [6]).

3.1.1 Expected value of $y_{2n}(\Phi)$

Using the expressions (60) and (66) we see that the ensemble average of $y_{2n}(\Phi)$ is given by

[TABLE]

The expectation value $\left\langle\exp\{i\Phi_{w}\}\right\rangle$ vanishes unless the total phase $\Phi_{w}$ of the non-backtracking walk $w$ obeys $\Phi_{w}=0$. That condition is met only if each edge is traversed both forwards and backwards the same number of times. The argument implies $\left\langle y_{n}(\Phi)\right\rangle=0$ for all odd $n$.

We display $y_{2n}$ pictorially in terms of subgraphs as in Figure 1. Vertices and edges that contribute to $y_{2n}$ are depicted as dots and as bonds, respectively. We characterize the topology of each subgraph in terms of the Betti number $\beta$ (not to be confused with the matrix index $\beta$ used in Section 1). Informally, $\beta$ counts the number of two-dimensional holes in the planar representation of a graph. For instance, the subgraphs displayed in Figure 1 both have $\beta=2$. For $n\ll N$ the dominant contributions to (68) come from walks in which each edge is traversed only once in each direction. Trees ($\beta=0$) and single loops ($\beta=1$) cannot occur as neither allows backtracking. Thus, the dominant contributions come from walks on subgraphs with $\beta=2$. There are two such subgraphs333Other forms of subgraphs with two loops do exist. These do not contribute to the same order, however, because their edges must be traversed more than twice. called type I and type II and shown in Figure 1. To describe these we introduce the following notation. The total number of vertices on a subgraph is $v$. Vertices located at the intersections are called junctions. The remaining vertices are called simple vertices. The latter are arranged in the form of linear chains that begin and terminate in a junction. Each such chain is called a branch. In subgraphs of types I and II there are at most three branches. The number of simple vertices on a branch is denoted by letters $p,q,r$. The relation between $n$ (the length of the walk), $v$ and $p,q,r$ depends on the topology of the subgraph (see Fig. 1).

• •

Type I:  There is one junction linked to four edges and there are two branches carrying $p$ and $q$ vertices, respectively. Each branch forms a loop connected to the junction. Then $n=v+1,\ v=p+q+1$ but $p,q\geq 2$ to avoid backscattering. Note that due to the constraints we must have $n\geq 6$.

• •

Type II:  There are two junctions linked to three edges each and three branches containing $p,q,r$ vertices, respectively. Then $n=v+1,\ v=p+q+r+2$. We may have either $p,q,r\geq 1$, or $p,q\geq 1$ and $r=0$ (cyclic). Note that even though a periodic walk of length $2n$ traverses the entire subgraph with each edge traversed in both directions, the walk of length $n$ is not a periodic walk, in contrast to walks in type I. Note that due to the constraints we must have $n\geq 6$ if $p,q,r\geq 1$ or $n\geq 5$ if one of $p,q$ or $r$ is zero.

We calculate the contributions to (68) from type I and type II subgraphs. For type I the contribution is determined by the total number $\left\langle{\rm Tr}[{\cal Y}(\Phi)^{2n}]\right\rangle_{\rm I}$ of walks $w=(e_{1},\ldots,e_{2t})$ that trace out the subgraph such that every edge is traversed once in each direction. Let $a$ denote the vertex at the junction, $p^{+}=(\mu_{1},\ldots,\mu_{p})$ the ordered vertices on the $p$-branch, $p^{-}=(\mu_{p},\ldots,\mu_{1})$ the reversed order, and analogously for $q$. With $e_{1}=(a,\mu_{1})$ fixed, there are two possible traversals of the subgraph given by $w=ap^{+}aq^{+}ap^{-}aq^{-}$ and $w=ap^{+}aq^{-}ap^{-}aq^{+}$. There are $2n$ possible directed edges from where to start, giving a total of $4n$ walks for each labelled subgraph. With the total number of vertices given by the dimension $N$ of ${\cal M}$, there are $N!/(N-v)!$ ways of choosing the vertices. We have to sum over all values of $(p,q)$ subject to the constraints $p+q=v-1$ and $p,q\geq 2$. To avoid overcounting (due to reversing of the orientations of walks in each branch and the exchange of the two branches) we have to divide by $s_{\rm I}=2\times 2\times 2=8$. That gives

[TABLE]

which is valid for $n\geq 6$.

To count the number $\left\langle{\rm Tr}[{\cal Y}(\Phi)^{2n}]\right\rangle_{\rm II}$ of walks $w$ that traverse subgraphs of type II, we denote the two vertices at the junctions by $a$ and $b$ and the branches by $p^{\pm},q^{\pm},r^{\pm}$ similarly as before. For every labelled subgraph, there are two possible traversals $ap^{+}bq^{-}ar^{+}bp^{-}aq^{+}br^{-}$ and $ap^{+}br^{-}aq^{+}bp^{-}ar^{+}bq^{-}$ and $2n$ directed edges from where to start, giving a total of $4n$ traversals. There are $N!/(N-v)!$ ways of choosing the vertices. In addition, we have to sum over $p,q,r$ subject to the constraints (either $p,q,r\geq 1$, or $p,q\geq 1$ and $r=0$ (cyclic)). To avoid overcounting we must divide by $s_{\rm II}=3!\times 2=12$. This comes from 3! ways of exchanging the three branches and from the reflection of the subgraph about its centre (i.e. exchanging $a\leftrightarrow b$ and relabelling $(p^{+},q^{+},r^{+})\leftrightarrow(p^{-},q^{-},r^{-})$). Thus, for $n\geq 6$ and $v=n-1$,

[TABLE]

For $n=5$ we only keep the second term on the right-hand side and so

[TABLE]

For $n=1,\ldots,4$ we have $\left\langle{\rm Tr}[{\cal Y}(\Phi)^{2n}]\right\rangle=0$ since there do not exist any non-backtracking paths of length 8 or less in which all the edges are traversed the same number of times in both directions. We may now combine the expressions (69), (70) and (71) for $n\geq 5$. So, using that $v=n-1$ and $N!/(N-v)!=N^{n-1}-\mathcal{O}(N^{n-2})$ for $N\gg v$ we get

[TABLE]

The relation (68) between $\left\langle{\rm Tr}[{\cal Y}(\Phi)^{2n}]\right\rangle$ and $\left\langle y_{2n}(\Phi)\right\rangle$ then leads to

[TABLE]

Therefore, from Eq. (66) the expectation value of the trace of the Chebyshev polynomial is

[TABLE]

whereas for $n<5$ the relation (66) gives $\left\langle{\rm Tr}[T_{2n}({\cal W})]\right\rangle=-\frac{N}{2}\frac{(N-3)}{(N-2)^{n}}$.

In comparison to (74) the equivalent expectation for the GUE can be obtained by combining the result (6) with the form of the Chebyshev polynomial (67) to obtain

[TABLE]

We highlight that for the GUE there is no constant term in the $1/N$ expansion. One may see immediately why this is the case from the form of the moments (6), which after multiplying by a factor of $N$ do not contain any constant term444This is not the case for other Gaussian $\beta$-ensembles however.. Whereas for $n=2$ the UME satisfies $\left\langle\Tr[T_{2n}({\cal W})]\right\rangle\to-\frac{1}{2}$ as $N\to\infty$. We could remove this constant term for $n=2$ by changing the scaling of ${\cal W}$ from $(N-2)^{-1/2}$, which arises naturally in (66). However, choosing another scaling of the form $(N-c)^{-1/2}$ for some constant $c$ will induce an order $\mathcal{O}(1/N)$ correction (from the leading term, which is 0 for all $n>1$) which will result in $\left\langle\Tr[T_{2n}({\cal W})]\right\rangle$ having constant terms for other values of $n$.

3.1.2 Correlations of $y_{n}(\Phi)$.

From expression (60) the covariance of the traces of powers of ${\cal Y}(\Phi)$ is given by

[TABLE]

Here $\Omega_{n,m}(0):=\{(w,w^{\prime})\in(\Omega_{n},\Omega_{m}):\Phi_{w}=-\Phi_{w^{\prime}}\neq 0\ \forall\ \Phi\}$ denotes the set of pairs of non back-tracking walks that have nonvanishing opposite phases for every ${\cal M}\in$ UME.

Variance of $y_{n}(\Phi)$

The dominant contribution to $|\Omega_{n,n}(0)|$ comes from pairs of walks $w,w^{\prime}$ that reside on a subgraph with a single loop in which every edge is traversed exactly once by $w$ and in the opposite direction by $w^{\prime}$. Hence there are $v=n$ vertices on the subgraph. The path $w$ can start from each edge and in either direction, giving $2n$ possible starting positions and then there are $N!/(N-v)!$ ways of labelling the vertices. That, however, overcounts by a factor $s_{\beta=1}=2n$ since all starting position of $w$ can also be obtained by relabelling the vertices. For every $w$ we have $n$ possible walks $w^{\prime}$ which gives, for $n\geq 3$

[TABLE]

The next-to-leading-order contribution to $|\Omega_{n,n}(0)|$ comes from pairs of non-backtracking walks on subgraphs containing two loops ($\beta=2$). As shown in the previous section there exist two types of subgraphs, for which we require $w$ to traverse each edge of these subgraphs precisely once and $w^{\prime}$ to traverse the same subgraph in the opposite direction. That is not possible for type II subgraphs as these do not support non-backtracking walks in which the loops are traversed only once. Thus leaving subgraphs of type I.

To obtain the contribution from the type I subgraphs (where $v=n-1$), we see that if a walk $w$ starts at a particular edge then there are two possible traversals of the subgraph, which in the notation of the previous section, are given by $w=ap^{+}aq^{+}$ and $w=ap^{+}aq^{-}$. There are $n$ possible starting edges for the walk $w$ and hence a further $n$ possible choices for $w^{\prime}$. We may also relabel the vertices, which gives a factor of $N!/(N-v)!$ but must gain mitigate for the overcounting. Thus we must divide through by a factor of 4, coming from the possibility of swapping the two loops and the two different possible traversals of the subgraph by $w$. So altogether, for $n\geq 6$

[TABLE]

Combining (77) and (78) leads to

[TABLE]

Therefore, using that $(N-2)^{-n}=N^{-n}-2N^{-n-1}+\ldots$ we have from (66) that

[TABLE]

This in turn implies that

[TABLE]

The leading term in this expression coincides with the result of Johansson [18] for the GUE. The manner in which it has been derived (i.e., counting the number of non-backtracking walks on single loops) is the same as in [6] (see also [1, 36, 37]). In fact, this approach is capable of showing that all the joint moments of the type $\left\langle{\rm Tr}[T_{1}({\cal W})]^{a_{1}}{\rm Tr}[T_{2}({\cal W})]^{a_{2}}\ldots{\rm Tr}[T_{k}({\cal W})]^{a_{k}}\right\rangle$ with $a_{i}>0$, coincide in the large $N$ limit with averages of the form $\left\langle Z_{1}^{a_{1}}Z_{2}^{a_{2}}\ldots Z_{k}^{a_{k}}\right\rangle=\left\langle Z_{1}^{a_{1}}\right\rangle\left\langle Z_{2}^{a_{2}}\right\rangle\ldots\left\langle Z_{k}^{a_{k}}\right\rangle$, where the $Z_{n}$ are independent and identically distributed Gaussian random variables with zero mean and variance $\sigma_{n}^{2}=n/4$.

Covariance of $y_{n}(\Phi)$ and $y_{m}(\Phi)$

The covariance is obtained by setting $n\neq m$ (we take $n>m$ without loss of generality) in (76) and computing the number of pairs of non-backtracking walks $|\Omega_{n,m}(0)|$ of length $n$ and $m$ which retrace each other. Obviously ${\rm Cov}({\rm Tr}[{\cal Y}(\Phi)^{n}],{\rm Tr}[{\cal Y}(\Phi)^{m}])=0$ if $(n-m)$ is odd since then the phases $\Phi_{w}$ and $\Phi_{w^{\prime}}$ cannot be equal in general.

For even $(n-m)$ the leading contribution comes from subgraphs of type I in which $w\in\Omega_{n}$ traverses one of the loops twice in opposite directions whilst $w^{\prime}\in\Omega_{m}$ only traverses the other loop once. For example $w=ap^{+}aq^{+}ap^{-}$ and $w^{\prime}=aq^{-}$. Thus one loop contains $m$ edges, the other $(n-m)/2$ edges and there are a total of $v=(n+m-2)/2$ vertices. The number of such pairs $(w,w^{\prime})$ is given by noting that the $w$ walk has $2n$ possible ways of traversing the subgraph, given by the 2 possible orders of traversing the $(n-m)/2$ loop and the $n$ possible starting edges. Then for each $w$ there are $m$ ways of choosing the $w^{\prime}$ walk. Relabelling the vertices also gives a factor of $N!/(N-v)!$ however we must then divide by a factor of $2$ to account for reversing the orientation of the loop that is traversed twice by $w$. Altogether this gives

[TABLE]

which in turn leads to

[TABLE]

3.2 Spectral density and two-point correlation function

3.2.1 Trace formula for the spectral density

Given any matrix ${\cal M}\in{\rm UME}$, the density of the scaled eigenvalues $\epsilon_{k}$ (see (65))

[TABLE]

is a distribution which we study by applying it to test functions (observables) that are analytic on the entire real line. We restrict our attention to this space of functions because the maximum (scaled) spectral radius (achieved by setting $M_{\mu\mu}=1$ for all $\mu\neq\nu$) for the UME is given by $[(N-1)/(2\sqrt{N-2})]\thicksim\sqrt{N}$.

Let $f(x)$ be an allowed test function. It can be expanded in terms of Chebyshev polynomials:

[TABLE]

The series converges on the entire real line, since it is a rearrangement of the Taylor expansion. We note, however, that the coefficients $f_{m}$ are derived from the restriction of $f(\epsilon)$ to the interval $[-1,1]$.

Recalling (66) and the fact that $y_{0}(\Phi)=(N-1),y_{1}(\Phi)=y_{2}(\Phi)=0$ we can write

[TABLE]

The sum over $y_{n}(\Phi)f_{n}$ converges because $|y_{n}|<(2\sqrt{N-1})^{n}$. Therefore, inserting the expression (82) for the coefficients $f_{n}$, using the absolute convergence of the series to exchange summation and integration, and noting that $T_{0}(\epsilon)=1$ leads to

[TABLE]

where going to the final line we have made use of the identity

[TABLE]

For appropriate test functions $f(x)$ and for any ${\cal M}\in$ UME (3.2.1) is an exact, absolutely convergent trace formula and provides the correct manner in which one can apply the formal trace formulae derived in [30].

The absolute convergence of (3.2.1) permits computing the ensemble average term by term. Due to (73), the infinite sum over $\left\langle y_{n}(\Phi)\right\rangle f_{n}$ is then of lower order in $N$ than the leading term, the integral, which does not depend on any periodic walk information. Therefore, in the limit of large $N$, the expression within the curly brackets in the integrand can be interpreted as the mean spectral density and indeed, in this limit, it converges to the semi-circle density.

3.2.2 Mean spectral density for finite $N$

The term appearing in the curly brackets in (3.2.1),

[TABLE]

can be identified as the mean spectral density. It includes $1/N$ corrections to the semi-circle law. However, the integration domain is the interval $[-1,1]$ so that the possible contributions for finite $N$ due to the non-Ramanujan part of the spectrum can only come from the periodic orbit sum in (3.2.1). In Figure 2 we compare the numerical mean density for $N=10$ with the expression given in (84), as well as the semi-circle density. None accounts for the oscillations that are due to the contributions from the periodic orbit sum.

To proceed further we introduce a family of $\delta$-like functions:

[TABLE]

where $c$ is a numerical constant, and ${N^{\star}}$ a large but finite integer. We consider $\delta_{N^{\star}}(x;\xi)$ a function of $x$ that depends on the parameter $\xi$, with $\xi$ restricted so as to have a finite distance from the end points of the interval $[-1,1]$. The function $\delta_{N^{\star}}(x;\xi)$ is concentrated about the point $x=\xi$ where it takes its maximum value $(N^{\star}+1)/(2(1+\frac{\pi}{2}\sqrt{1-\xi^{2}}))$ with full width at half maximum $[2\sqrt{3}/N^{\star}]\sqrt{1-\xi^{2}}$. The integral over the domain $[-1,1]$ is unity up to a correction of order $\mathcal{O}(1/N^{\star})$. For $x$ values sufficiently far from $\xi$, $\delta_{N^{\star}}(x;\xi)$ oscillates about [math] with mean amplitude of order $1$. The function $\delta_{N^{\star}}(x;\xi)$ is a polynomial in $x$ and therefore it belongs to the class of test functions relevant to our discussion. Because of these properties the function

[TABLE]

provides a smooth mean spectral density. The smoothing is done over spectral intervals of order $1/N^{\star}$ that lie within the domain $|\xi|<1-\frac{c}{{N^{\star}}}$. The coefficients $f_{m}$ for $\bar{\rho}_{N^{\star}}(\xi)$ are obviously proportional to $T_{m}(\xi)$ for $m\leq{N^{\star}}$ and zero otherwise. Thus the oscillatory part of the mean spectral density is

[TABLE]

Our numerical results (Figure 2) suggest that in the interval $[-1,1]$, the oscillatory part of the mean spectral density possesses $N$ oscillations about the mean, and this indicates that the parameter $N^{\star}$ should be at least $N$. Thus, in order to use (87) to match the data, one needs to compute $\left\langle y_{2m}(\Phi)\right\rangle$ at least for all $m\leq N$. Unfortunately the combinatorial computation (73), which provides an estimate for $\left\langle y_{2m}(\Phi)\right\rangle$, is valid only for $m\leq\sqrt{N}$.

3.2.3 Spectral form factor

For the discussion in the present Section it is convenient to map the spectral interval $-1\leq\epsilon\leq 1$ onto the unit circle by $\theta=2\arccos(\epsilon)$. We restrict the attention to Ramanujan matrices. Evaluating, formally, $N^{-1}{\rm Tr}[\delta(\epsilon-{\cal W})]$ in (3.2.1) gives us the spectral density, which may be written as

[TABLE]

The mean term here is given by the transformation of (84) and the oscillatory term by the sum over the $y_{n}(\Phi)$ in (3.2.1)

[TABLE]

The spectral two-point correlation function is defined as

[TABLE]

Using Eq. (3.2.3) we can write this as

[TABLE]

where we have substituted $\rho(\phi)$ for $\left\langle\rho_{{\cal M}}(\theta)\right\rangle$. Now,

[TABLE]

[TABLE]

[TABLE]

To leading order,

[TABLE]

The spectral form factor is defined as

[TABLE]

Substituting and collecting the terms, we get

[TABLE]

The spectral information used to calculate the form factor (92) consists of the entire set of values within the asymptotic support. The spectral density is not constant. Therefore, Eq. (92) cannot be compared directly with the GUE “local” form factor. The latter is defined for the unfolded spectrum with constant mean spectral density. For a comparison one has to transform the GUE result using a convolution integral introduced previously in [31]. The resulting GUE expression for $K_{2}(t,N)$ is a function of the scaled variable $\tau=t/N$ and is given by

[TABLE]

We observe that the leading term in (92) is $\tau$-inherited from the GUE expression. The next-order terms consist of odd powers of $\tau^{\frac{1}{2}}$. The numerical data and the expression for $K_{2}^{\rm GUE}(t/N)$ above are compared in Fig. 3 for $N=20$. The agreement is very good.

4 Brownian motion approach

In the previous section we obtained expressions for the average spectral moments of the UME in a basis of Chebyshev polynomials (see Eqn. (74)). By inverting the relation (67) one may obtain the monomials from this basis and hence obtain expressions for the average value of the expressions $\frac{1}{N}\left\langle\Tr[{\cal W}^{2n}]\right\rangle$, whose leading terms correspond to the moments of the semicircle distribution $\frac{2}{\pi}\sqrt{1-\lambda^{2}}$. The result (74) therefore allows us to obtain the deterministic deviations for averages $\frac{1}{N}\left\langle\Tr[f({\cal W})]\right\rangle$ of polynomial test functions over the UME from the average over the semicircle distribution, i.e.

[TABLE]

In the present section we turn to analysing the behaviour of fluctuations of the traces of Chebyshev polynomials about their mean. In particular we will show, using a Brownian motion approach, that in the large $N$ limit the random variables $\Tr[T_{n}({\cal W})]$ behave like independent Gaussian random variables (as was briefly highlighted at the end of Section 3.1.2) and provide bounds on the rates of this convergence utilising results derived from Stein’s method [9]. This is encapsulated in our main result of this section, which is Theorem 3.

Let us therefore first define our centred traces of Chebyshev polynomials by removing their mean

[TABLE]

Here $\Lambda_{n}:=\Omega_{n}\setminus\{w\in\Omega_{n}:\Phi_{w}=0\ \forall\ \Phi\in{\rm UME}\}$ is the set of periodic non-backtracking walks of length $n$ whose phases are not identically 0 for every member of the UME.

As was mentioned in the introduction, there are many examples of random matrix ensembles in which the linear-statistic $L_{f}(\mathcal{H})$ converges in distribution to a Gaussian random variable with a universal variance (which we do not state here but may be found in [18, 37] for example) as $N$ becomes large. Many works have sought to establish regularity conditions for these linear statistics, i.e. finding the class of test functions for which $L_{f}(\mathcal{H})$ converges to a Gaussian random variable in the large $N$ limit (see e.g. [33] and references therein).

In the case of polynomial functions $f$ of order $k$ this is equivalent to showing that the joint distribution of the first $k$ traces of Chebyshev polynomials $\{F_{n}(\mathcal{H})\}_{n=1}^{k}$ converge in distribution to independent Gaussian distributed random variables. We are going to show that this result is also true of the UME, except that, as noted in Section 3.2.1, we have $\Tr[{\cal Y}(\Phi)]=\Tr[{\cal Y}(\Phi)^{2}]=0$ and therefore $F_{0}({\cal M})=F_{1}({\cal M})=F_{2}({\cal M})\equiv 0$. For this reason it only makes sense to investigate $F_{n}({\cal M})$ for $n\geq 3$. One may consider this to be a substantial deviation from the GUE, however it is shown in [32] that one may always scale the first two moments in the Gaussian $\beta$-ensembles in such a way that they may be considered independently of all other moments.

For the UME (and more generally for Wigner matrices) the convergence of $F_{n}({\cal M})$ to independent Gaussians (see Theorem 1 below) is well known (see e.g. [6] for instance) and can be obtained by showing the convergence of all moments via combinatorial methods.

Theorem 1**.**

For ${\cal M}$ distributed according to the UME and $k$ fixed555Although we are not aware of a specific result for the UME, it has been shown for real Wigner matrices by Sinai and Soshnikov [20] that the traces $\Tr[\mathcal{H}^{k}]$ are still Gaussian distributed for growing $k$ (such that $k<N^{\frac{1}{2}}$) and one should expect a similar outcome here.

[TABLE]

in distribution as $N\to\infty$, where $Z_{n}$ are iid Gaussian variables with mean zero and variance $\left\langle Z_{n}^{2}\right\rangle=\frac{n}{4}$.

However the combinatorial approach does not easily lend itself to obtaining rates of convergence. To the best of our knowledge, the only result that discusses rates of convergence for Wigner matrices (including the UME) is by Chatterjee [21], in which he opts instead for an analytical approach that combines Stein’s method with estimates of second order Poincaré inequalities. In contrast, our result in Theorem 3 shows that one may still obtain similar convergence rates using combinatorial procedures and, moreover, one only requires estimates on finite moments of the matrix elements to achieve this.

In the present section we follow an Brownian motion approach pioneered by Dyson in his seminal article [8], which outlined an alternative description for the canonical Gaussian ensembles. His insight was that one could replace the (static) Gaussian matrix entries with independent (dynamical) Ornstein-Uhlenbeck processes such that the equilibration in the large time limit corresponds to the original Gaussian ensemble. The motion in the entries then induces a corresponding stochastic motion in the eigenvalues, now known as Dyson Brownian motion, whose stationary distribution is then given by the joint probability density function for the eigenvalues of the Gaussian ensemble.

In the setting of the UME (See Eq. (1)) we define our motion such that each of the $N(N-1)/2$ phases $\phi_{\mu\nu}$ is an independent standard Brownian motion on the torus $([0,2\pi))^{N(N-1)/2}$. Thus, in a small time $\delta s$ each phase moves an amount $\delta\phi_{\mu\nu}:=\phi_{\mu\nu}^{\prime}-\phi_{\mu\nu}$, which is characterised by its drift and diffusion, given by the two respective moments

[TABLE]

and all higher moments are of order $\mathcal{O}(\delta s^{2})$. Here $\rho(\Phi\to\Phi^{\prime};\delta s)=\prod_{\mu<\nu}\rho(\phi_{\mu\nu}\to\phi_{\mu\nu}^{\prime};s)$ denotes the probability of $\Phi=\{\phi_{\mu\nu}\}_{\mu<\nu}$ moving to $\Phi^{\prime}=\Phi+\delta\Phi$ in a time $\delta s$.

The above formulation of the motion is equivalent to considering the probability distribution $P(\Phi;s)$ of finding the particles at position $\Phi$ at time $s$, subject to some initial distribution $P(\Phi;0)$ at time 0. $P(\Phi;s)$ satisfies the following Fokker-Planck equation

[TABLE]

with periodic boundary conditions on the torus. Note that if $P(\Phi;0)=\delta(\Phi-\Phi^{\prime})=\prod_{\mu<\nu}\delta(\phi_{\mu\nu}-\phi^{\prime}_{\mu\nu})$, i.e. the particles are conditioned to be at position $\Phi$ at time 0, then $P(\Phi^{\prime};s)=\rho(\Phi\to\Phi^{\prime};s)$ is simply the transition probability. One may solve this equation explicitly (see e.g. [38, 39]), although we shall not need the exact form of the solution here. In the limit of large times the probability distribution satisfies

[TABLE]

which corresponds to the stationary distribution of the UME - the solution obtained once the left hand side of (97) is set to zero and appropriately normalised.

If the process is started from equilibrium, i.e. $P(\Phi;0)=(2\pi)^{-N(N-1)/2}$, then for all times $s$ we have $P(\Phi^{\prime};s)=P(\Phi;0)$, where $\Phi$ and $\Phi^{\prime}$ are related by $\rho(\Phi\to\Phi^{\prime};s)$. In the probability literature, the pair $(\Phi,\Phi^{\prime})$ is then termed an exchangeable pair, since both $\Phi$ and $\Phi^{\prime}$ have the same distribution (see e.g. [9]).

With the same philosophy as used by Dyson [8], a change $\delta\Phi$ in the matrix elements induces a corresponding change in our traces of Chebyshev polynomials of $\delta F_{n}:=F_{n}({\cal M}(\Phi^{\prime}))-F_{n}({\cal M}(\Phi))$, where $\Phi^{\prime}=\Phi+\delta\Phi$, in time $\delta s$. The key difference in the UME in comparison to the GUE, is that the motion is not invariant under unitary transformations666See [34] for further discussion on this point in the case of Bernoulli matrices. and consequently the evolution of $F_{n}({\cal M})$ is not closed (i.e. cannot be described entirely in terms of the $F_{n}({\cal M})$ themselves). Nevertheless we can estimate the remainder terms using the combinatorial procedures outlined in the previous section.

For large $N$ we will show in the following subsections that the drift and diffusion coefficients associated to the motion of $F_{n}({\cal M})$ approximately (in the probabilistic sense of Theorem 2 and Theorem 3 below) satisfy

[TABLE]

[TABLE]

with cross terms $\lim_{\delta s\to 0}\mathbb{E}[\delta F_{n}\delta F_{m}|\Phi]/\delta s\approx 0$ and higher terms identically zero. This means the process behaves in a similar way to an Ornstein-Uhlenbeck (OU) process generated by the following Fokker-Planck equation

[TABLE]

where $X=(X_{3},\ldots,X_{n})$. Moreover, since $\lim_{s\to\infty}P(\Phi;s)=(2\pi)^{-N(N-1)/2}$ is the stationary distribution of the UME, the associated stationary distribution of $F({\cal M})$ will, in turn, be approximately equal to the stationary solution of (100) given by

[TABLE]

Theorem 2 below (due to [Meckes [9]) makes this notion precise and for reasons of brevity we shall not repeat their proof here. However we would like to highlight that the result follows from the Taylor expansion of an appropriate observable $h(F({\cal M}))$, with $\delta h:=h(F({\cal M}^{\prime}))-h(F({\cal M}))$

[TABLE]

Dividing through by a factor $\delta s$ and then taking the limit $\delta s\to 0$ we see that from (98) and (99) the combination of the first two terms are (approximately) equal to $\mathcal{A}h(F({\cal M}))$, where

[TABLE]

Notice this is the adjoint operator to that in the right hand side of (100). Stein’s Lemma ensures that if for any twice-differentiable test function $h$ we have $\left\langle\mathcal{A}h(Z)\right\rangle=0$ then $Z=(Z_{3},\ldots,Z_{k})$ must be a multi-dimensional Gaussian random variable with $\left\langle Z_{n}\right\rangle=0$ and $\left\langle Z_{n}Z_{m}\right\rangle=n\delta_{nm}/4$. In our case the average over the UME is not identically 0 but $\left\langle\mathcal{A}h(F({\cal M})\right\rangle\approx 0$, which suggests that $F({\cal M})$ is approximately Gaussian. Stein’s method (originally developed in order to provide an alternative proof of the CLT [35]) then allows one to estimate the distance between $F({\cal M})$ and the Gaussian variable $Z$ in a suitable metric by bounding the variance of $\mathcal{A}h(F({\cal M})$.

Definition 1** (Wasserstein distance).**

Let us denote $\mathcal{L}:=\{f:\mathbb{R}^{k}\to\mathbb{R}:|f(x)-f(y)|\leq\|x-y\|\}$ to be the set of all Lipschitz continuous functions and $X,Y$ be two $k$-dimensional random variables, then the Wasserstein distance between $X$ and $Y$ is

[TABLE]

The Wasserstein distance provides a particular way to measure the distance between two probability distributions and it often emerges as a natural distance when utilising Stein’s method. Moreover, if one has sequence of random variables $X_{N}$ in which $d_{\rm W}(X_{N},Y)\to 0$ then this implies that $X_{N}\to Y$ in distribution.

The following theorem, utilising Stein’s method, is due to Meckes. Note that their results are more general than the following statement but we adapt it to our setting for purposes of clarity.

Theorem 2** (Meckes [9]).**

Let ${\cal M}$ and ${\cal M}^{\prime}$ be two random matrices with the same probability distribution (they are an exchangeable pair) and related via some transition probability $\rho({\cal M}\to{\cal M}^{\prime};s)$. Let $X\equiv X({\cal M})=(X_{3}({\cal M}),\ldots,X_{k}({\cal M}))$ and $X^{\prime}\equiv X({\cal M}^{\prime})=(X_{3}({\cal M}),\ldots,X_{k}({\cal M}^{\prime}))$ be two $k-2$ dimensional random variables dependent of ${\cal M}$ and ${\cal M}^{\prime}$. If

[TABLE] 2. 2.

[TABLE] 3. 3.

[TABLE]

for all $n,m,l=3,\ldots,k$ and $R_{n}({\cal M})$ and $R_{nm}({\cal M})$ are (potentially) random variables depending on ${\cal M}$, then

[TABLE]

with $Z$ the multi-dimensional Gaussian random variable stated in Theorem 1.

We also comment that an alternative version of this theorem by Döbler and Stolz [40] has also been used by Webb [41] to estimate the Wasserstein distance between the traces in the circular $\beta$-ensembles and Gaussian random variables.

Theorem 3**.**

Let ${\cal M}$ be distributed according to the UME (see Equation (1)) and let $F({\cal M})=(F_{3}({\cal M}),\ldots,F_{k}({\cal M}))$ be defined as in (4). Then, for $k$ fixed, we have

[TABLE]

with $Z$ that of Theorem 1.

Proof.

In the following subsections we will show the remainders for our drift (98) and diffusion (99) terms satisfy $\left\langle|R_{n}({\cal M})|\right\rangle=\mathcal{O}(N^{-1})$ and $\left\langle|R_{nm}({\cal M})|\right\rangle=\mathcal{O}(N^{-\frac{1}{2}})$ respectively. Incorporating these estimates into the Wasserstein distance (101) in Theorem 2 then gives the result. ∎

We remark that a convergence rate of order $\mathcal{O}(N^{-\frac{1}{2}})$ is also found in [21] for certain classes of Wigner matrices using the total-variation metric.

4.1 Drift term

Before proceeding, let us first introduce the following notation. Let $E:=\{e=(\mu,\nu):\mu<\nu\}$ be the set of directed edges on our graph such that $e(\sigma)=(\mu,\nu)$ if $\sigma=+$ and $(\nu,\mu)$ if $\sigma=-$. Then for a non-backtracking walk $w=(e_{1}(\sigma_{1}),e_{2}(\sigma_{2}),\ldots,e_{n}(\sigma_{n}))$ the total phase can be written as

[TABLE]

where $\kappa^{(w)}_{e}=\#\{e(+)\in w\}-\#\{e(-)\in w\}$ counts the net number of traversals of the edge $e$ by $w$ (it may be positive, negative or zero). This means, using the properties of the motion (95) and (96), we have

[TABLE]

Therefore, using the form of $F_{n}({\cal M})$ from (4) we have

[TABLE]

where, writing $x_{w}:=\sum_{e\in w}(\kappa^{(w)}_{e})^{2}$ for simplicity, the remainder is given by

[TABLE]

and $\Lambda^{\prime}_{n}:=\{w\in\Lambda_{n}:x_{w}\neq n\}$. In particular, this excludes those $n$-periodic non-backtracking walks in which every edge is only traversed once. To estimate the value of this remainder we must compute

[TABLE]

By averaging over the phases we find the main contributions to (103) will come from pairs of walks in which $w=w^{\prime}$. These, however, cannot be walks in which all the edges are distinct (as this would imply $x_{w}=n$). Therefore, the main contribution if from those $w$ containing precisely one edge that is traversed twice (once in each direction) and the remaining edges are connected to this edge by two loops in which every edge is traversed once. Taking $w^{\prime}$ to be the same walk but in the opposite direction means we have, using the notation from Section 3.1.1, $v=n-\beta+1=n-2$ vertices (note that we have effectively $\beta=3$ loops since traversing an edge twice can be viewed as creating an additional loop). Following the arguments of Section 3.1.1 and Section 3.1.2 this means $\sum_{w,w^{\prime}\in\Lambda^{\prime}_{n}}x_{w}x_{w^{\prime}}\left\langle\exp\{i(\Phi_{w}-\Phi_{w}^{\prime})\}\right\rangle=\mathcal{O}(N^{n-2})$ and thus $\left\langle|R_{n}({\cal M})|\right\rangle=\mathcal{O}(N^{-1})$.

4.2 Diffusion term

We now show that the remainder for our diffusion term satisfies $\left\langle|R_{nm}({\cal M})|\right\rangle=\mathcal{O}(N^{-1/2})$. To begin, using the notations above, we note that

[TABLE]

Thus, writing $x_{w,w^{\prime}}=\sum_{e}\kappa^{(w)}_{e}\kappa^{(w^{\prime})}_{e}$ and taking the complex conjugate of $F_{n}({\cal M})$ in the following (since it is real), we have for the diffusion term

[TABLE]

where

[TABLE]

In order to obtain an estimate for $\left\langle|R_{nm}({\cal M})|\right\rangle$ we treat the cases $n=m$ and $n\neq m$ separately.

For $n=m$ we have

[TABLE]

Expanding out the brackets inside the square root above gives

[TABLE]

[TABLE]

Now, we first note that $x_{w,w^{\prime}}=0$ if $w$ and $w^{\prime}$ do not share an edge. Therefore the main contribution to the first sum in the above comes from pairs of non-backtracking walks $w_{1}=w_{2}$ and $w_{3}=w_{4}$ that reside on disconnect subgraphs comprised of a single loop. As was determined in (77), and using that $x_{w_{1},w_{1}}=n$ for such walks, means the first summation gives a contribution $n^{4}+\mathcal{O}(N^{-1})$ and the second summation $n^{2}+\mathcal{O}(N^{-1})$. Therefore, taking the square root, we have $\left\langle|R_{nn}({\cal M})|\right\rangle=\mathcal{O}(N^{-1/2})$.

For $n\neq m$ we have

[TABLE]

Due to the presence of $x_{w_{1},w_{2}}$ and $x_{w_{3},w_{4}}$ it must be the case that for $\left\langle\exp\{i(\Phi_{w_{1}}-\Phi_{w_{2}}+\Phi_{w_{3}}-\Phi_{w_{4}})\}\right\rangle x_{w_{1},w_{2}}x_{w_{3},w_{4}}$ to be non-zero $w_{1}$ and $w_{2}$ must share at least one edge, and similarly for $w_{3}$ and $w_{4}$. Therefore, since $w_{1}$ and $w_{2}$ are of different lengths, those subgraphs supporting the walks $(w_{1},\ldots,w_{4})$ in which every edge is traversed twice must have at least three loops. The number of vertices in such a situation is therefore $v=n+m-2$ and hence the contribution inside the square root is of order $\mathcal{O}(N^{-2})$, meaning $\left\langle|R_{nm}({\cal M})|\right\rangle=\mathcal{O}(N^{-1})$.

4.3 Remainder term

In order to complete the necessary conditions for Theorem 3 requires verifying that our motion satisfies Part (3) in Theorem 1. This is easily shown, since

[TABLE]

This comes from noting that the main contribution to $\mathbb{E}[(\delta F_{n}\delta F_{m}\delta F_{l})^{2}|{\cal M}]$ will come from terms in which the same edges appears six times. Therefore, taking this conditional expectation we find $\mathbb{E}[\delta\phi_{e}^{6}|\Phi]=\mathcal{O}(\delta s^{3})$, because $\delta\phi_{e}$ behaves like a Gaussian in the small time limit with variance $\sigma^{2}\propto\delta s$. Thus, dividing through by $\delta s$ and taking the limit we get the result

[TABLE]

as required.

5 Conclusions

In this paper we have studied the spectral density and moments of the UME for large matrix dimension $N$. The UME is a paradigmatic example of a Wigner ensemble where the distribution of the matrix elements is not Gaussian. The spectral properties of the UME agree with those of the GUE in the limit $N\to\infty$. We have focused attention on contributions of next order(s) in an asymptotic expansion in $1/N$. These are of interest in their own right, possibly showing deviations from universality. Their study also casts light on the power of the approaches used for their study. We have used three very different approaches: The supersymmetry approach, the graph-theoretical approach, and the Brownian-motion approach.

Using the supersymmetry approach and the saddle-point approximation we have shown that the leading $1/N$ corrections to the spectral density of the UME are the same as for the GUE and account for the oscillations that feature so prominently in the numerical data. We have not been able to push the supersymmetry approach beyond that point. The graph-theoretical approach has yielded leading $1/N$ corrections to the mean values of Chebyshev moments and their covariances. These differ from GUE values. That information did not suffice, however, to account for the oscillations in the spectral density. Progress might be possible via more complicated combinatorial computations. These would go beyond the scope of the present paper. Finally we used a Brownian motion approach combined with Stein’s method to analyse the fluctuations of the spectral moments. In particular the combinatorial ideas outlined in Section 3 allowed us to estimate the error terms and provide bounds on the rate of convergence to a multi-dimensional Gaussian in the large $N$ limit.

One issue worth mentioning, that has not been addressed here, is that due to gauge invariance the spectrum of any matrix in the UME depends only on the magnetic fluxes $\Phi_{c}=\sum_{i}\phi_{v_{i},v_{i+1}}$ on the fundamental cycles $c$. The number of such fundamental cycles is $\beta=\frac{1}{2}N(N-1)-N+1$, which is less than the number of independent phases by $N-1$. We did not attempt to make use of the gauge invariance or the freedom in choosing the independent cycles.

Acknowledgements

CHJ and US would like to thank Professor Sasha Sodin for many discussions and illuminating comments. CHJ is also grateful to the Leverhulme Trust (ECF-2014-448) for financial support.