\documentstyle[12pt]{article} \headheight 0cm \headsep 0cm \newlength{\mytopmargin} \newlength{\myleftmargin} \setlength{\mytopmargin}{2.8cm} \setlength{\myleftmargin}{2.5cm} \setlength{\topmargin}{-1in} \setlength{\oddsidemargin}{-1in} \addtolength{\topmargin}{\mytopmargin} \addtolength{\oddsidemargin}{\myleftmargin} \textwidth 16cm \textheight 24cm \setlength{\parindent}{1.5em} \renewcommand{\baselinestretch}{1.4} \begin{document} \title{ A nonlinear equation and its application to nearest neighbor spacings for zeros of the zeta function and eigenvalues of random matrices} \author{P.J. Forrester\thanks{email: matpjf@maths.mu.oz.au; supported by the ARC}\\ Department of Mathematics\\ University of Melbourne\\ Parkville, Victoria 3052\\ Australia \and A.M. Odlyzko\thanks{email: amo@research.att.com} \\AT \& T Bell Laboratories\\ Murray Hill \\ New Jersey 07974\\U.S.A.} \date{} \maketitle \begin{abstract} A nonlinear equation generalizing the $\sigma$ form of the Painlev\'e V equation is used to compute the probability density function for the distance from an eigenvalue of a matrix from the GUE ensemble to the eigenvalue nearest to it. (The classical results concern distribution of the distances between consecutive eigenvalues.) Comparisons are made with the corresponding distribution for zeros of the Riemann zeta function, which are conjectured to behave like eigenvalues of large random GUE matrices. \end{abstract} \noindent {\bf 1. INTRODUCTION} \noindent The so called GUE hypothesis (see e.g.~[1]) states that, in a certain limit, the zeros of the Riemann zeta function on the critical line Re$(z)=1/2$ have the same joint distribution as that of the eigenvalues of a random matrix from the Gaussian Unitary Ensemble (GUE) of large (formally infinite) dimensional random Hermitian matrices. Denoting the zeros by $1/2 + i \gamma_n$, where $n$ labels the zeros sequentially along the critical line, the GUE hypothesis applies in the limit $n \rightarrow \infty$, with each $\gamma_n$ scaled by the mean density of zeros at $1/2 + i \gamma_n$ $(= (1/2\pi) \log(\gamma_n /2 \pi))$, so that the mean spacing between zeros is unity (any finite value will do). We recall (see e.g.~[2]) that a random Hermitian $N \times N$ matrix is said to belong to the GUE if the diagonal elements $x_{jj}$(which must be real) and the upper triangular elements $x_{jk} = u_{jk} + i v_{jk}$ are independently chosen with probability density function (p.d.f.) $$ {1 \over \sqrt{\pi}} e^{-x_{jj}^2} \quad \mbox{and} \quad {2 \over \pi} e^{-2(u_{jk}^2 + v_{jk}^2)} = {2 \over \pi} e^{-2|x_{jk}^2|} \eqno (1.1) $$ respectively. For large $N$ the density of eigenvalues $\rho(\lambda)$ is given by the so called Wigner semi-circle law $$ \rho(\lambda) \sim {\sqrt{2N} \over \pi} \sqrt{1 - {\lambda^2 \over 2N}}. \eqno (1.2) $$ To apply the GUE hypothesis the eigenvalues should therefore be scaled by $\sqrt{2N}/\pi$ before the $N \rightarrow \infty$ is taken, to obtain a mean eigenvalue spacing of unity. One consequence of the GUE hypothesis is that it provides concrete predictions for statistical properties of $\{ \gamma_n \}$, whenever these are known for the GUE random matrices. One such example is the p.d.f., $p(s)$ say, for the spacing between consecutive zeros. In the infinite GUE, scaled so that the mean eigenvalue spacing is $1/\pi$, the corresponding quantity is given in terms of a Fredholm determinant of an integral operator (see e.g.~[2]) by $$ p(2t) = {1 \over 4} {d^2 \over d t^2} \det (1 - K) \eqno (1.3a) $$ where $K$ is the integral operator on the interval $(-t,t)$ with kernel $$ K(x,y) := {\sin (x-y) \over \pi (x-y)}. \eqno (1.3b) $$ Using an eigenvalue expansion of the Fredholm determinant, $p(s)$ can be computed [3] (see also [1]) to give a tabulation or graph of $p(s)$. In a large-scale numerical computation of the non-trivial zeros of the Riemann zeta function by one of the present authors [1], involving over $10^7$ consecutive values of $\gamma_n$ about $n = 10^{20}$ calculated to an accuracy of about six decimal places, the p.d.f.~$p(s)$ has been determined empirically and compared with $p(s)$ as calculated from (1.3). Excellent agreement is found. Similar agreement is found when comparing other empirical statistical distributions with those known exactly for the GUE. In this paper we present the exact calculation of a new statistical quantity for the infinite GUE, which is similar in meaning to $p(s)$. This statistical quantity is the p.d.f., $nn(t)$ say, for the spacing between nearest neighbor eigenvalues. We then test the GUE hypothesis by comparing $nn(t)$ for the infinite GUE with the empirical calculation of $nn(t)$ for the zeros the Riemann zeta function on the critical line from the data of [1] for $\{ \gamma_n \}$. \vspace{.3cm} \noindent {\bf 2. A NON-LINEAR EQUATION} \vspace{.2cm} \noindent {\bf 2.1 Summary of results} \noindent So as to put our calculation of $nn(t)$ in context, we first note that the expression (1.3) for $p(s)$ has been made more explicit by Jimbo et al.~[4], who proved that $$ \det (1 - K) = \exp \int_0^t { \sigma (2 t') \over t'} dt' \eqno (2.1) $$ where $\sigma (s)$ satisfies the $\sigma$ form of the Painlev\'e V equation: $$ (s \sigma'')^2 + 4 (s \sigma' - \sigma)((\sigma')^2 - \sigma + s \sigma') = 0 \eqno (2.2) $$ subject to the boundary condition $$ \sigma (s) \: \sim \: -s/\pi -(s/\pi)^2 \qquad {\rm as \quad s \rightarrow 0}. \eqno (2.3) $$ Subsequent derivations of this result have been given by Its et al.~[5], Mehta [6] and Tracy and Widom [7]. Our expression for $nn(t)$ is given in terms of the solution of a non-linear equation which generalizes (2.2). We have obtained the following results. The p.d.f.~$nn(t)$ for the infinite GUE is given in terms of a Fredholm determinant by $$ nn(t) = - {d \over d t} \det (1 - K_1) \eqno (2.4) $$ where $K_1$ is the integral operator on $(-t,t)$ with kernel $$ K_1(x,y) := { \sqrt{xy} \over 2(x-y)} \Big ( J_{b + 1/2}(x) J_{b - 1/2}(y) - J_{b+1/2}(y) J_{b - 1/2}(x) \Big ) \eqno (2.5) $$ ($J_\alpha (x)$ denotes the Bessel function) and $b=1$. Furthermore $$ \det (1 - K_1) = \exp \int_0^{\pi \rho t} { \sigma_1 (2 t') \over t'} dt' \quad \mbox{ and so} \quad nn(t) = -{ \sigma_1(2 \pi \rho t) \over t} \exp \int_0^{\pi \rho t} {\sigma_1(2t') \over t'} dt' \eqno (2.6) $$ (here the mean eigenvalue spacing is $1/\rho$), where $\sigma_1 (s)$ satisfies the non-linear equation $$ (s \sigma_1'')^2 + 4 (-b^2 + s \sigma_1' - \sigma_1) \Big ( (\sigma_1')^2 + \{ b - (b^2 - s \sigma_1' + \sigma_1)^{1/2} \}^2 \Big ) = 0 \eqno (2.7) $$ with $b=1$, subject to the boundary condition $$ \sigma_1 (s) \: \sim \: -{(s/2)^{2b + 1} \over \Gamma (1/2 + b) \Gamma(3/2 + b)} \quad \mbox{as} \quad s \rightarrow 0 \eqno (2.8) $$ with $b=1$ (the parameter $b$ is included above for later convenience). Note that with $b=0$ (2.7) reduces to (2.2). \vspace{.2cm} \noindent {\bf 2.2 Derivation} \noindent Our derivation of (2.4)-(2.8) uses a recent result of Nagao and Slevin [8] to obtain (2.4), and the theory of Tracy and Widom [9] to obtain (2.7). Nagao and Slevin consider the random matrix ensemble with unitary symmetry defined by the eigenvalue p.d.f. $$ \prod_{j=1}^N |x_j|^{2 b} e^{-x_j^2} \prod_{1 \le j < k \le N} |x_k - x_j|^2, \qquad b > -1/2. \eqno (2.9) $$ They prove that in the thermodynamic limit, with each $x_j$ scaled $x_j \mapsto X_j/\sqrt{2N}$ so that the bulk density is $1/\pi$, the corresponding $n$-particle distribution is given by $$ \rho_n (X_1, \dots, X_n) = \det[ K_1 (X_j, X_k)]_{j,k = 1, \dots, n} \eqno (2.10) $$ where $K_1(x,y)$ is given by (2.5) (for $b \notin Z_{\ge 0}$, $x(y) < 0$, $x(y)$ in the denominator needs to be replaced by $|x|(|y|)$, however below we will only consider the case $b \in Z_{\ge 0}$). It follows by substituting (2.10) into the general formula $$ E_0(-t,t) = 1 + \sum_{n=1}^\infty {(-1)^n \over n!} \int_{-\infty}^\infty dX_1 \dots \int_{-\infty}^\infty dX_n \, \rho_n (X_1, \dots, X_n), $$ where $E_0(-t,t)$ denotes the probability of an interval $(-t,t)$ being free of eigenvalues in the ensemble (2.9), that $$ E_0(-t,t) = \det (1 - K_1). \eqno (2.11) $$ Since in the case $b=1$ (2.9) is precisely the eigenvalue probability density function of the GUE with an eigenvalue fixed at the origin, the result (2.4) follows. The derivation of (2.7) relies on a set of coupled equations for quantities associated with the integral operator $K_1$. To present these equations, the particular quantities must first be defined. \vspace{.2cm} \noindent {\bf Definition 2.1} \noindent Suppose $A$ is an integral operator on the interval $(a_1,a_2)$ with kernel $A(x,y)$: $$ A f = \int_{a_1}^{a_2} A(x,y) f(y) \, dy. $$ We write $$ A \doteq A(x,y) $$ to specify that the kernel of $A$ is $A(x,y)$. Denote by $K$ an integral operator of this type with kernel of the form $$ K(x,y)={\phi(x) \psi(y) - \phi(y) \psi(x) \over x - y} $$ and write $$ (1 - K)^{-1} \doteq \rho (x,y) \qquad K(1 - K)^{-1} \doteq R (x,y). $$ Also let \begin{eqnarray*} Q_k(x) &:= &(1 - K)^{-1} y^k \phi := \int_{a_1}^{a_2} \rho(x,y) y^k \phi(y), \quad q_{kj} = Q_k(a_j) := \lim_{x \rightarrow a_j \atop x \in (a_1,a_2)} Q_k(x) \\ P_k(x)& := &(1 - K)^{-1} y^k \psi := \int_{a_1}^{a_2} \rho(x,y) y^k \psi(y), \quad p_{kj} = P_k(a_j) := \lim_{x \rightarrow a_j \atop x \in (a_1,a_2)} P_k(x) \\ u &:= &\int_{a_1}^{a_2} Q_0(y) \phi (y) \, dy, \quad w := \int_{a_1}^{a_2} P_0(y) \psi (y) \, dy \\ v &:= & \int_{a_1}^{a_2} Q_0(y) \psi (y) \, dy = \int_{a_1}^{a_2} P_0(y) \phi (y) \, dy. \end{eqnarray*} where $k=0,1$ and $j=1,2$. \vspace{.3cm} The coupled equations which imply (2.7) can now be stated. \vspace{.2cm} \noindent {\bf Proposition 2.1} \noindent Consider the kernel (2.5) with $b \in Z_{\ge 0}$ on $(-t,t)$ so that in the setting of Definition 2.1 $a_1 = -t$ and $a_2 = t$. With the notation $q := q_{02}$, $p := p_{02}$, $R(t,t) = R$ we have $$ \begin{array}{ll} \mbox{(i)} \quad tR = 2(-b + u -w)pq + t(p^2 + q^2) + 2(pq)^2 \qquad & \mbox{(ii)} \quad tq' = (-b+u-w)q + tp \\ \mbox{(iii)} \quad tp' = -tq - (-b +u -w)p \qquad & \mbox{(iv)} \quad (tR)' = p^2 + q^2 \\ \mbox{(v)} \quad u' = 2q^2, \quad w' = 2p^2 \qquad & \mbox{(vi)} \quad \Big ( \log (1 - K_1) \Big )' = -2R \end{array} $$ where the dashes denote differentiation with respect to $t$. \vspace{.3cm} The theory of Tracy and Widom [9] allows equations for the quantities of Definition 2.1 to be derived which imply the equations of Proposition 2.1. Before presenting these equations let us show how (2.6) and (2.7) can be derived from the equations of Proposition 2.1. First consider (2.7). We multiply (ii) by $p$, multiply (iii) by $q$, add and use (v) to obtain $$ (pq)' = p^2 - q^2 = {1 \over 2} (w' - u') \eqno (2.12) $$ and consequently $$ pq ={1 \over 2} (w - u). \eqno (2.13) $$ Substituting (2.13) and (iv) in (i) gives $$ tR = - 2b (pq) - 2(pq)^2 +t(tR)' \eqno (2.14) $$ which relates $tR$ to $pq$. On the other hand, another equation relating these two quantities is obtained by squaring (iv) and the first equality in (2.12) and subtracting: $$ ((pq)')^2 - ((tR)')^2 = - 4 (pq)^2. \eqno (2.15) $$ Solving (2.14) for $pq$ (it follows from a small-$t$ expansion that the negative square root is to be taken) and $(pq)'$, substituting in (2.15) and introducing the notation $$ \sigma_1(2t) := - 2t R \eqno (2.16) $$ gives (2.7). The boundary condition (2.8) follows from the fact that $R(s,s) \sim K_1(s,s)$ as $s \rightarrow 0$ and the corresponding behaviour of $K_1(s,s)$ deduced from (7). To derive (2.6) we simply substitute (2.16) in (vi) and integrate. The factor $\pi \rho$ in the upper terminal of (2.6) results from changing the mean eigenvalue spacing from $1/\pi$ to $1/\rho$. \vspace{.3cm} \noindent {\bf 2.3 Theory of Tracy and Widom} \noindent To derive the equations in Proposition 2.1 we use the theory of [9] to first obtain some equations for the quantities in Definition 2.1 considered as functions of the end points $a_1$, $a_2$. These equations are of two types: those which apply independent of the particular functions $\phi$ and $\psi$, and those which are dependent on $\phi$ and $\psi$. The equations of interest which fall into the first category are summarized as follows. \vspace{.2cm} \noindent {\bf Proposition 2.2} \noindent For general values of $\phi$ and $\psi$ in Definition 2.1 we have $$ R(a_j,a_k) = {q_{0j}p_{0k} - p_{0j} q_{0k} \over a_j - a_k} \quad (j \ne k), \qquad {\partial \over \partial a_j} \log \det (1 - K_1) = (-1)^{j-1}R(a_j,a_j) $$ and $$ {\partial u \over \partial a_k} = (-1)^k (q_{0k})^2, \qquad {\partial v \over \partial a_k} = (-1)^k p_{0k} q_{0k}, \qquad {\partial w \over \partial a_k} = (-1)^k (p_{0k})^2, $$ $$ {\partial q_j \over \partial a_k} = (-1)^k R(a_j,a_k) q_k, \quad {\partial p_j \over \partial a_k} = (-1)^k R(a_j,a_k) p_k, \qquad (j \ne k) $$ where $j,k = 1,2$. \vspace{.3cm} To present the second type of equations, note that the kernel $K_1(x,y)$ as given by (2.5) is of the type in Definition 2.1 with $$ \phi (x) = \sqrt{x/2} J_{b+1/2}(x), \quad \psi (x) = \sqrt{x/2} J_{b-1/2}(x). \eqno (2.17) $$ In [9] Tracy and Widom show that equations further to those in Proposition 2.2 exist whenever $\phi$ and $\psi$ satisfy the coupled differential equations $$ m(x) \phi'(x) = A(x) \phi (x) + B(x) \psi (x), \quad m(x) \psi'(x) = -C(x) \phi (x) - A(x) \psi (x). \eqno (2.18) $$ for $m,A,B,C$ polynomials. For the choice (2.6), (2.7) hold with $$ m(x) = x, \quad A(x) = \alpha_0, \quad B(x) = \beta_1 x, \quad C(x) = \gamma_1 x. \eqno (2.19a) $$ with $$ \alpha_0 = -b, \quad \beta_1 = \gamma_1 = 1. \eqno (2.19b) $$ For general $\alpha_0, \beta_1, \gamma_1$ we can read off from the results of [9] additional equations relating the quantities $R(a_j,a_j)$, $q_{kj}$, $p_{kj}$ ($k=0,1, \quad j=1,2$) and $u,v,w$. \vspace{.2cm} \noindent {\bf Proposition 2.3} \noindent Consider the kernel $K$ of Definition 2.1 with $\phi$ and $\psi$ defined by (2.18) with $m,A,B,C$ as in (2.19a). We have $$ \mbox{(a)} \quad q_{1j} = a_j q_{0j} - (v q_{0j} - u p_{0j}) \quad [9,{\rm eq.(2.12)}] \hspace{6.7cm} $$ $$ \mbox{(b)} \quad p_{1j} = a_j p_{0j} - (w q_{0j} - v p_{0j}) \quad [9,{\rm eq.(2.13)}] \hspace{6.5cm} $$ \begin{eqnarray*} \mbox{(c)} \quad a_j {\partial q_{0j} \over \partial a_j} & = & (\alpha_0 + \gamma_1 u) q_{0j} + \beta_1 p_{1j} + \beta_1 v p_{0j} \\ & & - \sum_{k=1 \atop k \ne j}^2 (-1)^k a_k R(a_j,a_k) q_{0k} \quad [9,\mbox{ eq.(2.25) with modification (2.31)}] \end{eqnarray*} \begin{eqnarray*} \mbox{(d)} \quad a_j {\partial p_{0j} \over \partial a_j} & = & -\gamma_{1}q_{1j} + \gamma_1 v q_{0j} - \alpha_0 p_{0j} + \beta_1 w p_{0j} \\ & & - \sum_{k=1 \atop k \ne j}^2 (-1)^k a_k R(a_j,a_k) p_{0k} \quad [9,\mbox{ eq.(2.26) with modification (2.31)}] \vspace{1cm} \end{eqnarray*} \begin{eqnarray*}\lefteqn{ \mbox{(e)} \quad a_jR(a_j,a_j) = (\alpha_0 + \gamma_1 u) q_{0j}p_{0j} + (\beta_1 p_{1j}+ \beta_1 v p_{0j})p_{0j}}\\& & \hspace*{1.4cm} + (\gamma_1 q_{1j} - \gamma_1 v q_{0j}) q_{0j} + (\alpha_0 p_{0j} - \beta_1 w p_{0j}) q_{0j} \\ & & \hspace*{1.4cm} + \sum_{k=1 \atop k \ne j}^2 (-1)^k a_k {(q_{0j} p_{0k} - p_{0j} q_{0k})^2 \over a_j - a_k} \quad [9,\mbox{ eq.(2.27) with modification (2.32)}] \end{eqnarray*} \begin{eqnarray*} \mbox{(f)} \quad {\partial \over \partial a_j}a_jR(a_j,a_j) & = & \beta_1 p_{0j}^2 + \gamma_1 q_{0j}^2 - \sum_{k=1 \atop k \ne j}^2 (-1)^k a_k\Big ( R(a_j,a_k) \Big )^2 \hspace{3.5cm}\\ & & [9,\mbox{ eq.(2.28) with modification (2.32)}], \end{eqnarray*} where $j=1,2$. \vspace{.3cm} To pursue our task of deriving the equations in Proposition 2.1, let us return to the particular case $$ a_1 = -t, \quad a_2 = t, \eqno (2.20) $$ and $\alpha_0$, $\beta_1$, $\gamma_1$ given by (2.8b) so that $\phi$ and $\psi$ are given by (2.6). Since $$ \phi(-x) = (-1)^{b-1}\phi(x), \qquad \psi(-x) = (-1)^{b}\psi(x) \eqno (2.21) $$ (recall $b \in Z_{\ge 0}$) we have $K(x,y)=K(-x,-y)$ and thus $\rho(x,y) = \rho (-x,-y)$, which together with (2.21) implies $$ q_{01} = (-1)^{b-1}q_{02} = (-1)^{b-1}q, \quad p_{01} = (-1)^{b}p_{02} = (-1)^{b}p, \quad v=0. \eqno (2.22) $$ With the aid of (2.22) the equations in Proposition 2.1 can now be deduced in a straightforward manner from the equations in Proposition 2.2 and 2.3. We first use (a) and (b) to eliminate $q_{1j}$ and $p_{1j}$ in (c)-(e) of Proposition 2.3. Equation (i) of Proposition 2.1 now follows by choosing $j=2$ and substituting (2.20) and (2.22). The equations (ii)-(iv) of Proposition 2.1 are deduced from (c),(d) and (f) of Proposition 2.3 respectively. This requires making use of the general formula $$ {d \over dt} f(-t,t) = \Big ( {\partial \over \partial a_2} - {\partial \over \partial a_1} \Big ) f(a_1, a_2) \Big |_{a_2 = -a_1 =t}, \eqno (2.23) $$ using (2.22), and noting from the first equation in Proposition 2.2 with the substitutions (2.20) and (2.22) that $$ R(-t,t) = (-1)^b {pq \over t}. $$ The equations (v) follow immediately from the final line of equations in Proposition 2.2 and (2.22), and the final equation (vi) follows from the second equation in Proposition 2.2 and (2.23). \vspace{.2cm} \noindent {\bf 3. STATISTICS OF THE ZEROS OF THE RIEMANN ZETA FUNCTION} \vspace{.2cm} \noindent {\bf 3.1 Solution of the non-linear equation} \noindent We have computed many terms of the power series expansion of (2.7) about $s=0$ with $b=1$ and subject to (2.8). Substitution of the first seven into (2.4) gives $$ nn(t) ={ 2(\pi \rho t)^2 \over 3} - {4(\pi \rho t)^4 \over 45} +{2 (\pi \rho t)^6 \over 315} - {32 (\pi \rho t)^7 \over 2025 \pi } + O(t^8). \eqno (3.1) $$ Comparison with the analogous expansion of (1.3a) (see e.g.~[2]) shows that $p(t) - (1 /2) nn(t) = O(t^7)$ , which in qualitative terms says that very small spacings between consecutive eigenvalues will most likely be nearest neighbour spacings (the factor of 1/2 accounts for the fact that the nearest neighbour occurs with equal probability to the left or the right). For the large-$t$ expansion, we follow the corresponding analysis of $\sigma(s)$ (see e.g.~[7]) and seek a solution of the form $as^2 + bs + c + d/s + \dots$. Indeed (2.7) has a unique solution of this form, which when integrated according to (2.6) implies $$ \det (1 - K_1) \sim {A \over (\pi \rho t)^{5/4}} \exp \Big ( - (\pi \rho t)^2 / 2 + 2 \pi \rho t + 1 /(4 \pi \rho t) + O(1/t^2) \Big ) \eqno (3.2) $$ As is well known (see e.g.~[7]) this leaves the overall multiplicative constant $A$ unspecified. The large-$t$ expansion of $nn(t)$ is obtained by substituting (3.2) in (2.4). The solution of (2.7) with $b=1$ was computed numerically by first calculating the power series expansion of $\sigma_1 (s)$ up to $O(s^{12})$ and using the corresponding values of $\sigma_1 (1)$ and $\sigma_1' (1)$ as initial data in the Mathematica routine NDSolve. The d.e.~(2.7) was rewritten so that $\sigma''_1$ occured to the first power (the negative square root is to be taken) and it was found necessary to use a high precision setting (AccuracyGoal and PrecisionGoal = 20) to get a stable solution in the interval of interest ($s < 13$). Although the results of Section 2 are only exact in the limit $N \rightarrow \infty$, it is well known that $p(s)$ can be accurately approximated by considering $2 \times 2$ matrices which give the so called Wigner surmise (see e.g.~[2]). This suggests that $nn(t)$ may also be insensitive to the precise dimension of the GUE matrices. Assuming this, to test our exact expression we have compared $nn(t)$ as calculated from (2.6) with $nn(t)$ as determined empirically from 10,000 numerically generalted $15 \times 15$ matrices from the GUE (see Figure 1). The latter calculation was done using Mathematica. For each matrix the eigenvalues were calculated and the nearest neighbour spacing of the middle (8th) eigenvalue was calculated. After scaling the spacings were tested to count how many fell into the intervals $((k-1)/20,k/20), k=1,2,\dots$, and the corresponding empirical value of $nn(t)$ plotted at the points $(k-1/2)/15$. \vspace{.3cm} \noindent{\bf 3.2 Empirical value of $nn(s)$ for Riemann zeta function zeros} \noindent With the p.d.f.~$nn(t)$ for the infinite GUE now evaluated, we can further test the GUE hypothesis by calculating $nn(t)$ for sequences of consecutive zeros $1/2 + i \gamma_n$ of the Riemann zeta function, using the data of [1]. Three sets of $10^6$ consecutive zeros $1/2 + i \gamma_n$ were analyzed, the data sets starting at zero number $N_1 = 1$, $N_2 = 10^6 + 1$ and $N_3 = 10^{20} + 143,782,842$ respectively. The quantity $\delta_n' := \min (\delta_n, \delta_{n-1})$, where $\delta_n := (\gamma_{n+1} - \gamma_n)\rho_n$ with $\rho_n = (1/2 \pi) \log (\gamma_n / 2\pi)$ denoting the smoothed local density of zeros at $1/2 + i\gamma_n$, was calculated and a histogram was constructed for the number of values out of the $10^6$ that were tested that fell into the intervals $((k-1)/20, k/20)$, $k=1,2 \dots$. In Figure 2 the corresponding empirical values of $nn(t)$ at the points $(k-1/2)/20$ are plotted and compared with the value of $nn(t)$ for the infinite GUE. The convergence towards the GUE value as the magnitude of the imaginary part increases is evident. For further comparison the moments $ \langle t^p \rangle := \int_0^\infty t^p nn(t) \, dt \quad (a=1,2,3), $ for $p=1, \dots, 10$ were calculated and compared with the empirical data according to the law of large numbers prediction $ \langle t^p \rangle \approx \langle \delta_n'^p \rangle_a := 10^{-6} \sum_{n = N_a +1}^{N_a + 10^6} \delta_n'^p. $ The results are contained in Table 1. Again the trend is towards convergence to the GUE value. Note in particular the four figure agreement between $\langle t \rangle$ and $\langle \delta_n' \rangle_3$. \pagebreak \begin{description} \item[References] \item[][1] A.M. Odlyzko, The $10^{20}$-th zero of the Riemann zeta function and 70 million of its neighbors, unpublished preprint; On the distribution of spacings between zeros of the zeta function, {\it Math Comput.} {\bf 48} (1987) 273-308 \item[][2] M.L. Mehta, {\it Random Matrices}, 2nd ed. (Academic, New York, 1991) \item[][3] M.L. Mehta and J. des Cloizeaux, The probabilities for several consecutive eigenvalues of a random matrix, {\it Indian J. Pure Appl. Math} {\bf 3} (1972), 329-351. \item[][4] M. Jimbo, T. Miwa, Y. M\^{o}ri and M. Sato, Density matrix of an impenetrable Bose gas and the fifth Painlev\'e transcendent, {\it Physica D} {\bf 1} (1980) 80-158 \item[][5] A.R. Its, A.G. Izergin, V.E. Korepin and N.A. Slavnov, Differential equations for quantum correlation functions, {\it Int. J. Mod. Phys. B} {\bf 4} (1990) 1003-1037 \item[][6] M.L. Mehta, J. Phys.(France), A non-linear differential equation and a Fredholm determinant, {\bf 2} (1992) 1721-1729 \item[][7] C.A. Tracy and H. Widom, Introduction to Random matrices, in {\it Lect. Notes in Physics} {\bf 424} (1993) 407-424 \item[][8] T. Nagao and K. Slevin, Nonuniversal correlations for random matrix ensembles, J. Math. Phys. {\bf 34} (1993) 2075-2085 \item[][9] C.A. Tracy and H. Widom, Fredholm determinants, differential equations and matrix models, {\it Commun. Math. Phys.} {\bf 163} (1994) 33-72 \end{description} \pagebreak \vspace{1cm} \begin{tabular}{|c|c|c|c|c|} \hline $p$ & $\langle t^p$ $\rangle$ & $\langle \delta'^p_n \rangle_1$ & $\langle \delta'^p_n \rangle_2$ & $\langle \delta'^p_n \rangle_3 $\\ \hline \hline 1& 0.725227 &0.731988 &0.730706 &0.725291 \\ \hline 2 &0.603251 &0.606386 &0.605762 &0.602470\\ \hline 3 &0.555775 &0.551262 &0.551956 &0.553540\\ \hline 4 &0.555527 &0.540113 &0.542599 &0.551074\\ \hline 5 &0.594314 &0.563548 &0.568467 &0.586454\\ \hline 6 &0.674002 &0.620786 &0.629172 &0.660788\\ \hline 7 &0.804518 &0.717187 &0.730735 &0.782709\\ \hline 8 &1.00515 &0.864325 &0.885824 &0.969281\\ \hline 9 &1.30870 &1.08177 &1.11583 &1.24935\\ \hline 10 &1.76924 &1.40075 &1.45504 &1.67002\\ \hline \end{tabular} \vspace{.5cm} \noindent {\bf Table 1. } Comparison of the moments of $nn(t)$ for the GUE (second column) and for $10^6$ consecutive zeros of the Riemann zeta function (subsequent columns) on the critical line, starting near zero number 1, $10^6$ and $10^{20}$ respectively. \vspace{2cm} \noindent {\bf Figure captions} \vspace{1cm} \noindent {\bf Figure 1} Comparison of $nn(t)$ for the infinite GUE (solid line) and for 10,000 15 $\times$ 15 computer generated matrices from the GUE (filled circles). \vspace{1cm} \noindent {\bf Figure 2} Comparison of $nn(t)$ for the GUE (solid line) and for $10^6$ consecutive zeros of the Riemann zeta function on the critical line, starting near zero number 1 (open circles), $10^6$ (asterisks) and $10^{20}$ (filled circles) respectively. \end{document}