A nonlinear equation and its application to

nearest neighbor spacings for zeros of the zeta

function and eigenvalues of random matrices

P.J. Forrester^{∗} Department of Mathematics

University of Melbourne Parkville, Victoria 3052 Australia

A.M. Odlyzko^{†} AT&T Bell Laboratories Murray Hill New Jersey 07974 U.S.A.

Abstract

A nonlinear equation generalizing the σ form of the Painlevé 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.

1. INTRODUCTION

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γ_{n}, where n labels the zeros sequentially along the critical line, the GUE hypothesis applies in the limit n → ∞, with each γ_{n} scaled by the mean density of zeros at 1/2 + iγ_{n} (= (1/2π) log(γ_{n}/2π)), so that the mean spacing between zeros is unity (any finite value will do).

∗email: matpjf@maths.mu.oz.au; supported by the ARC †email: amo@research.att.com

We recall (see e.g. [2]) that a random Hermitian N × 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} + iv_{jk} are independently chosen with probability density function (p.d.f.)

\frac{1}{\sqrt{π}} e^{−x_{jj}^{2}} and \frac{2}{π} e^{−2(u_{jk}^{2}+v_{jk}^{2})} = \frac{2}{π} e^{−2|x_{jk}^{2}|}  (1.1)

respectively. For large N the density of eigenvalues ρ(λ) is given by the so called Wigner semi-circle law

ρ(λ) ∼ \frac{\sqrt{2N}}{π} \sqrt{1 − \frac{λ^{2}}{2N}} .  (1.2)

To apply the GUE hypothesis the eigenvalues should therefore be scaled by \sqrt{2N}/π before the N → ∞ 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 {γ_{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/π, the corresponding quantity is given in terms of a Fredholm determinant of an integral operator (see e.g. [2]) by

p(2t) = \frac{1}{4} \frac{d^{2}}{dt^{2}} det(1 − K)  (1.3a)

where K is the integral operator on the interval (−t, t) with kernel

K(x, y) := \frac{sin(x − y)}{π(x − y)}.  (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 γ_{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 {γ_{n}}.

2. A NON-LINEAR EQUATION

2.1 Summary of results

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 ∫_{0}^{t} \frac{σ(2t′)}{t′}dt′  (2.1)

where σ(s) satisfies the σ form of the Painlevé V equation:

(sσ′′)^{2} + 4(sσ′ − σ)((σ′)^{2} − σ + sσ′) = 0  (2.2)

subject to the boundary condition

σ(s) ∼ −s/π − (s/π)^{2} as s → 0.  (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) = −\frac{d}{dt} det(1 − K_{1})  (2.4)

where K_{1} is the integral operator on (−t, t) with kernel

K_{1}(x, y) := \frac{\sqrt{xy}}{2(x − y)}(J_{b+1/2}(x)J_{b−1/2}(y) − J_{b+1/2}(y)J_{b−1/2}(x))  (2.5)

(J_{α}(x) denotes the Bessel function) and b = 1. Furthermore

det(1 − K_{1}) = exp ∫_{0}^{πρt} \frac{σ_{1}(2t′)}{t′}dt′ and so nn(t) = −\frac{σ_{1}(2πρt)}{t} exp ∫_{0}^{πρt} \frac{σ_{1}(2t′)}{t′}dt′

(2.6)

(here the mean eigenvalue spacing is 1/ρ), where σ_{1}(s) satisfies the non-linear equation

(sσ′′_{1})^{2} + 4(−b^{2} + sσ′_{1} − σ_{1})((σ′_{1})^{2} + {b − (b^{2} − sσ′_{1} + σ_{1})^{1/2}}^{2}) = 0  (2.7)

with b = 1, subject to the boundary condition

σ_{1}(s) ∼ −\frac{(s/2)^{2b+1}}{Γ(1/2 + b)Γ(3/2 + b)} ass → 0  (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).

2.2 Derivation

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.

∏_{j=1}^{N} |x_{j}|^{2b}e^{−x_{j}^{2}} ∏_{1≤j<k≤N} |x_{k} − x_{j}|^{2},  b > −1/2.  (2.9)

They prove that in the thermodynamic limit, with each x_{j} scaled x_{j} ↦ X_{j}/\sqrt{2N} so that the bulk density is 1/π, the corresponding n-particle distribution is given by

ρ_{n}(X_{1}, …, X_{n}) = det[K_{1}(X_{j}, X_{k})]_{j,k=1,...,n}  (2.10)

where K_{1}(x, y) is given by (2.5) (for b ∉ Z_{≥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 ∈ Z_{≥0}).

It follows by substituting (2.10) into the general formula

E_{0}(−t, t) = 1 + ∑_{n=1}^{∞} \frac{(−1)^{n}}{n!} ∫_{−∞}^{∞} dX_{1} … ∫_{−∞}^{∞} dX_{n} ρ_{n}(X_{1}, …, X_{n}),  (2.10)

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}).  (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.

Definition 2.1 Suppose A is an integral operator on the interval (a_{1}, a_{2}) with kernel A(x, y):

Af = ∫_{a_{1}}^{a_{2}} A(x, y)f(y) dy.

We write

A ≐ 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) = \frac{ϕ(x)ψ(y) − ϕ(y)ψ(x)}{x − y}

and write

(1 − K)^{−1} ≐ ρ(x, y) K(1 − K)^{−1} ≐ R(x, y).

Also let

Q_{k}(x) := (1 − K)^{−1}y^{k}ϕ := ∫_{a_{1}}^{a_{2}} ρ(x, y)y^{k}ϕ(y), q_{kj} = Q_{k}(a_{j}) := lim_{{x→a_{j} \atop x∈(a_{1},a_{2})}} Q_{k}(x)
P_{k}(x) := (1 − K)^{−1}y^{k}ψ := ∫_{a_{1}}^{a_{2}} ρ(x, y)y^{k}ψ(y), p_{kj} = P_{k}(a_{j}) := lim_{{x→a_{j} \atop x∈(a_{1},a_{2})}} P_{k}(x)
u := ∫_{a_{1}}^{a_{2}} Q_{0}(y)ϕ(y) dy, w := ∫_{a_{1}}^{a_{2}} P_{0}(y)ψ(y) dy
v := ∫_{a_{1}}^{a_{2}} Q_{0}(y)ψ(y) dy = ∫_{a_{1}}^{a_{2}} P_{0}(y)ϕ(y) dy.

where k = 0, 1 and j = 1, 2.

The coupled equations which imply (2.7) can now be stated. Proposition 2.1 Consider the kernel (2.5) with b ∈ Z_{≥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

(i) tR = 2(−b + u − w)pq + t(p^{2} + q^{2}) + 2(pq)^{2}  (ii) tq′ = (−b + u − w)q + tp
(iii) tp′ = −tq − (−b + u − w)p  (iv) (tR)′ = p^{2} + q^{2}
(v) u′ = 2q^{2}, w′ = 2p^{2}  (vi) ( log(1 − K_{1}))′ = −2R

where the dashes denote differentiation with respect to t.

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} = \frac{1}{2}(w′ − u′)  (2.12)

and consequently

pq = \frac{1}{2}(w − u).  (2.13)

Substituting (3.2) and (iv) in (i) gives

tR = −2b(pq) − 2(pq)^{2} + t(tR)′  (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}.  (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

σ_{1}(2t) := −2tR  (2.16)

gives (2.7). The boundary condition (2.8) follows from the fact that R(s, s) ∼ K_{1}(s, s) as s → 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 πρ in the upper terminal of (2.6) results from changing the mean eigenvalue spacing from 1/π to 1/ρ.

2.3 Theory of Tracy and Widom

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 ϕ and ψ, and those which are dependent on ϕ and ψ.

The equations of interest which fall into the first category are summarized as follows. Proposition 2.3 For general values of ϕ and ψ in Definition 2.1 we have

R(a_{j}, a_{k}) = \frac{q_{0j}p_{0k} − p_{0j}q_{0k}}{a_{j} − a_{k}} (j ≠ k), \frac{∂}{∂a_{j}} log det(1 − K_{1}) = (−1)^{j−1}R(a_{j}, a_{j})

and

\frac{∂u}{∂a_{k}} = (−1)^{k}(q_{0k})^{2}, \frac{∂v}{∂a_{k}} = (−1)^{k}p_{0k}q_{0k}, \frac{∂w}{∂a_{k}} = (−1)^{k}(p_{0k})^{2}.

where j, k = 1, 2.

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

ϕ(x) = \sqrt{x/2}J_{b+1/2}(x), ψ(x) = \sqrt{x/2}J_{b−1/2}(x).  (2.17)

In [9] Tracy and Widom show that equations further to those in Proposition 2.3 exist whenever ϕ and ψ satisfy the coupled differential equations

m(x)ϕ′(x) = A(x)ϕ(x) + B(x)ψ(x), m(x)ψ′(x) = −C(x)ϕ(x) − A(x)ψ(x). (2.18)

for m, A, B, C polynomials. For the choice (2.6), (2.7) hold with

m(x) = x, A(x) = α_{0}, B(x) = β_{1}x, C(x) = γ_{1}x.  (2.19a)

with

α_{0} = −b, β_{1} = γ_{1} = 1.  (2.19b)

For general α_{0}, β_{1}, γ_{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, j = 1, 2) and u, v, w. Proposition 2.4 Consider the kernel K of Definition 2.1 with ϕ and ψ defined by (2.18) with m, A, B, C as in (2.19a). We have

(a) q_{1j} = a_{j}q_{0j} − (vq_{0j} − up_{0j}) [9, eq.(2.12)]

(b) p_{1j} = a_{j}p_{0j} − (wq_{0j} − vp_{0j}) [9, eq.(2.13)]

(c) a_{j} \frac{∂q_{0j}}{∂a_{j}} = (α_{0} + γ_{1}u)q_{0j} + β_{1}p_{1j} + β_{1}vp_{0j}

− ∑_{{k=1 \atop k≠i}}^{2}(−1)^{k}a_{k}R(a_{j}, a_{k})q_{0k} [9, eq.(2.25) with modification (2.31)]

(d) a_{j} \frac{∂p_{0j}}{∂a_{j}} = −γ_{1i} + γ_{1}vq_{0j} − α_{0}p_{0j} + β_{1}wp_{1j}

− ∑_{{k=1 \atop k≠i}}^{2}(−1)^{k}a_{k}R(a_{j}, a_{k})p_{0k} [9, eq.(2.26) with modification (2.31)]

(e) a_{j}R(a_{j}, a_{j}) = (α_{0} + γ_{1}u)q_{0j}p_{0j} + (β_{1}p_{1j} + β_{1}vp_{0j})p_{0j}
+(γ_{1}q_{1j} − γ_{1}vq_{0j})q_{0j} + (α_{0}p_{0j} − β_{1}wp_{0j})q_{0j}

+ ∑_{{k=1 \atop k≠i}}^{2}(−1)^{k}a_{k} \frac{(q_{0j}p_{0k} − p_{0j}q_{0k})^{2}}{a_{j} − a_{k}} [9, eq.(2.27) with modification (2.32)]

(f) \frac{∂}{∂a_{j}} a_{j}R(a_{j}, a_{j}) = β_{1}p_{0j}^{2} + γ_{1}q_{0j}^{2} − ∑_{{k=1 \atop k≠i}}^{2}(−1)^{k}a_{k}R(a_{j}, a_{k})

[9, eq.(2.28) with modification (2.32)],

where j = 1, 2.

To pursue our task of deriving the equations in Proposition 2.1, let us return to the particular case

a_{1} = −t, a_{2} = t,  (2.20)

and α_{0}, β_{1}, γ_{1} given by (2.8b) so that ϕ and ψ are given by (2.6). Since

ϕ(−x) = (−1)^{b−1}ϕ(x),  ψ(−x) = (−1)^{b}ψ(x)  (2.21)

(recall b ∈ Z_{≥0}) we have K(x, y) = K(−x, −y) and thus ρ(x, y) = ρ(−x, −y), which together with (2.21) implies

q_{01} = (−1)^{b−1}q_{02} = (−1)^{b−1}q, p_{01} = (−1)^{b}p_{02} = (−1)^{b}p, v = 0.  (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.3 and 2.4.

We first use (a) and (b) to eliminate q_{1j} and p_{1j} in (c)-(e) of Proposition 2.4. 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.4 respectively. This requires making use of the general formula

\frac{d}{dt}f(−t, t) = (\frac{∂}{∂a_{2}} − \frac{∂}{∂a_{1}} )f(a_{1}, a_{2})|_{a_{2}=−a_{1}=t},  (2.23)

using (2.22), and noting from the first equation in Proposition 2.3 with the substitutions (2.20) and (2.22) that

R(−t, t) = (−1)^{b} \frac{pq}{t} .

The equations (v) follow immediately from the final line of equations in Proposition 2.3 and (2.22), and the final equation (vi) follows from the second equation in Proposition 2.3 and (2.23).

3. STATISTICS OF THE ZEROS OF THE RIE-

MANN ZETA FUNCTION

3.1 Solution of the non-linear equation

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) = \frac{2(πρt)^{2}}{3} − \frac{4(πρt)^{4}}{45} + \frac{2(πρt)^{6}}{315} − \frac{32(πρt)^{7}}{2025π} + O(t^{8}).  (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 σ(s) (see e.g. [7]) and seek a solution of the form as^{2} + bs + c + d/s + . . .. Indeed (2.7) has a unique solution of this form, which when integrated according to (2.6) implies

det(1 − K_{1}) ∼ \frac{A}{(πρt)^{5/4}} exp ( − (πρt)^{2}/2 + 2πρt + 1/(4πρt) + O(1/t^{2}))  (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 σ_{1}(s) up to O(s^{12}) and using the corresponding values of σ_{1}(1) and σ′_{1}(1) as initial data in the Mathematica routine NDSolve. The d.e. (2.7) was rewritten so that σ′′_{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 → ∞, it is well known that p(s) can be accurately approximated by considering 2 × 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×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, …, and the corresponding empirical value of nn(t) plotted at the points (k − 1/2)/15.

3.2 Empirical value of nn(s) for Riemann zeta function zeros

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γ_{n} of the Riemann zeta function, using the data of [1]. Three sets of 10^{6} consecutive zeros 1/2 + iγ_{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 δ′_{n} := min(δ_{n}, δ_{n−1}), where δ_{n} := (γ_{n+1} − γ_{n})ρ_{n} with ρ_{n} = (1/2π) log(γ_{n}/2π) denoting the smoothed local density of zeros at 1/2 + iγ_{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 . . .. In Figure 2 the corresponding empirical values of nn(t) at the points (k − 1/2)/20 are plotted as a histogram 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 ⟨t^{p}⟩ := ∫_{0}^{∞} t^{p}nn(t) dt (a = 1, 2, 3), for p = 1, …, 10 were calculated and compared with the empirical data according to the law of large numbers prediction ⟨t^{p}⟩ ≈ ⟨δ_{n}^{′p}⟩_{a} := 10^{−6} ∑_{n=N_{a}+1}^{N_{a}+10^{6}} δ_{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 ⟨t⟩ and ⟨δ′_{n}⟩_{3}.

Table 1

\begin{matrix}p & ⟨t^{p} ⟩ & ⟨δ_{n}^{′p}⟩_{1} & ⟨δ_{n}^{′p}⟩_{2} & ⟨δ_{n}^{′p}⟩_{3} \\ 1 & 0.725227 & 0.731988 & 0.730706 & 0.725291 \\ 2 & 0.603251 & 0.606386 & 0.605762 & 0.602470 \\ 3 & 0.555775 & 0.551262 & 0.551956 & 0.553540 \\ 4 & 0.555527 & 0.540113 & 0.542599 & 0.551074 \\ 5 & 0.594314 & 0.563548 & 0.568467 & 0.586454 \\ 6 & 0.674002 & 0.620786 & 0.629172 & 0.660788 \\ 7 & 0.804518 & 0.717187 & 0.730735 & 0.782709 \\ 8 & 1.00515 & 0.864325 & 0.885824 & 0.969281 \\ 9 & 1.30870 & 1.08177 & 1.11583 & 1.24935 \\ 10 & 1.76924 & 1.40075 & 1.45504 & 1.67002 \end{matrix}

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.

Figure captions

Figure 1 Comparison of nn(t) for the infinite GUE (curved line) and for 10,000 15 × 15 computer generated matrices from the GUE (histogram).

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.

References

[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, Math Comput. 48 (1987) 273-308

[2] M.L. Mehta, Random Matrices, 2nd ed. (Academic, New York, 1991)

[3] M.L. Mehta and J. des Cloizeaux, The probabilities for several consecutive eigenvalues of a random matrix, Indian J. Pure Appl. Math 3 (1972), 329-351.

[4] M. Jimbo, T. Miwa, Y. Môri and M. Sato, Density matrix of an impenetrable Bose gas and the fifth Painlevé transcendent, Physica D 1 (1980) 80-158

[5] A.R. Its, A.G. Izergin, V.E. Korepin and N.A. Slavnov, Differential equations for quantum correlation functions, Int. J. Mod. Phys. B 4 (1990) 1003-1037

[6] M.L. Mehta, J. Phys.(France), A non-linear differential equation and a Fredholm determinant, 2 (1992) 1721-1729

[7] C.A. Tracy and H. Widom, Introduction to Random matrices, in Lect. Notes in Physics 424 (1993) 407-424

[8] T. Nagao and K. Slevin, Nonuniversal correlations for random matrix ensembles, J. Math. Phys. 34 (1993) 2075-2085

[9] C.A. Tracy and H. Widom, Fredholm determinants, differential equations and matrix models, Commun. Math. Phys. 163 (1994) 33-72