Level Spacings in Random Matrix Theory and Coulomb Blockade Peaks in Quantum Dots

# Level Spacings in Random Matrix Theory and Coulomb Blockade Peaks in Quantum Dots

## Abstract

We obtain analytic formulae for the spacing between conductance peaks in the Coulomb blockade regime, based on the universal Hamiltonian model of quantum dots. New random matrix theory results are developed in order to treat correlations between two and three consecutive spacings in the energy level spectrum. These are generalizations of the Wigner surmise for the probability distribution of single level spacing. The analytic formulae are shown to be in good agreement with numerical evaluation.

###### pacs:
73.23.Hk,05.45.Mt,73.21.La,02.10.Yn

## I Introduction

The Coulomb blockade of electrons has been a remarkable tool for probing fundamental physics in nanoscale systems review (); IhnBook (); glazmanreview (); oregreview (). In a semiconductor quantum dot (QD), weakly coupled to leads via tunneling, the effect is manifested as a series of spikes in the conductance of the device as a function of a gate voltage which controls the number of electrons on the QD. In the quantum Coulomb blockade regime (defined below), the spacing between the spikes is determined by the ground state energy of the QD. Experimentally, the spacing shows mesoscopic fluctuations, and there have been several detailed studies of the statistical distribution of the peak spacings sivan (); patel1 (); patel2 (); simmel (); luscher01 (); ihn02 (); fuhrer03 (); ong ().

The distribution of peak spacings results from the interplay of electron-electron interaction and randomness. A typical semiconductor QD in the experiments of Ref. patel1, consists of a droplet of several hundred electrons confined to a two dimensional region of a few tenths of a micron in size. In the absence of electron-electron interaction, the motion of the electrons would be randomized by impurity scattering or by chaotic scattering from the boundaries of the QD; hence, it is expected BohigasLesHouches (); EfetovBook () that the single particle energy levels in the absence of interaction would be described by random matrix theory (RMT) mehta (). In a real QD, however, it is necessary to take the additional effects of interactions into account. This is indicated by the failure of (essentially) non-interacting models to account for the observed distribution of Coulomb blockade peak spacings sivan (); patel1 ().

The interplay of electron-electron interaction and randomness is a notoriously hard problem. Nonetheless, due to the finite size of the QD, it has been possible to make significant progress. In particular, it is now believed that the energy levels of a weakly interacting QD are described statistically by a “universal Hamiltonian” (see Refs. glazmanreview, and oregreview, , and references therein). According to this model, each QD is characterized by a set of single particle orbitals with single particle energies and wavefunctions distributed according to the appropriate ensemble of RMT. Given the occupation numbers of these single particle orbitals, the non-interacting contribution to the energy of a QD follows directly. The interaction contribution is determined entirely by the net charge on the QD (the “charging energy”) and its total spin (the “exchange energy”). Thus the universal Hamiltonian is

 H=∑iσϵi^niσ+e22Cδ^N2−J^S2, (1)

where 0 or 1 is the occupation number of orbital with spin , is the number of excess electrons on the QD, and is the total spin-operator. denotes the capacitance of the QD, and is the exchange constant. Strictly, there is a fourth term in the universal Hamiltonian corresponding to interaction in the Cooper channel. However, this term vanishes when time reversal symmetry is broken by the application of a magnetic field, and, even when present, it is significantly smaller than the others and may be neglected glazmanreview ().

It should be emphasized that the universal Hamiltonian does not help calculate the energy levels of any particular QD. Rather it provides a description of the universal statistical features of the levels of all QDs that belong to the universality class under consideration.

The form of the universal Hamiltonian is dictated by the few symmetries that remain to a random system such as a QD kurland (). For example, Eq. (1) applies provided the QD has spin-rotation invariance (no spin-orbit or spin-flip scattering). If in addition time reversal symmetry is intact, the single particle levels are distributed according to the orthogonal ensemble of RMT; if time reversal symmetry is broken (by a magnetic field, for example), the unitary ensemble applies.

If the exchange constant is set to zero, the universal Hamiltonian reduces to the old constant interaction model which has proven unsuccessful in accounting for the observed peak spacing distribution. It is instructive to compare the physics of the two models. According to the constant interaction model, each filled QD level is doubly occupied in the ground state, except for the top level, which would be singly occupied if the total number of electrons in the dot were odd. Thus the total spin of the ground state is zero or one-half depending on whether the number of electrons is even or odd. Within the universal Hamiltonian model, higher spin ground states are possible: the higher single particle energy of these states is offset by the exchange energy which favors parallel alignment of spins brouwer (); baranger (). Indeed if the exchange constant is above a certain threshold, the QD should have a ferromagnetic ground state that is fully spin-polarized kurland (); AndKam98 (). Here we restrict our attention to the paramagnetic regime in which the ground state may have a substantial spin but is not fully polarized.

It is also instructive to compare the physics of QDs to that of atoms. Atoms can be approximately understood by considering that the electrons occupy a set of self-consistent hydrogen-like orbitals. The precise filling of the incomplete shells (which determines the total spin of the atomic ground state) is then controlled by the interplay of electron-electron interaction and spin-orbit coupling encoded in Hund’s rules. The universal Hamiltonian provides a comparable description of the energy levels of weakly interacting QDs. The key difference between atoms and QDs is that for atoms the confining potential is spherically symmetric whereas for QDs it is essentially random.

The universal Hamiltonian is found to give a good quantitative account of the experiments of Ref. patel1, when finite temperature effects and small non-universal corrections are included ong (); ullmo01 (); gonzalo (); UsajB02 (); AlhassidM02 (). Here, we obtain analytic formulae for the distribution of conductance peak spacings in the low temperature and large QD limit. Though the distribution has been evaluated numerically before ullmo01 (); gonzalo (); UsajB02 (); AlhassidM02 (), analytic expressions should prove helpful in the analysis of experimental data.

The quest to develop such analytic expressions raises some interesting problems in RMT. Within the constant interaction model, the distribution of peak spacings, , is controlled by the distribution of the spacing between consecutive energy levels, denoted . Level spacing distributions are notoriously hard to calculate but by now this distribution is well understood mehta (). The exact analytic expression for is far too complicated to be useful in practice. Fortunately, the celebrated Wigner surmise provides a simple and remarkably accurate approximation to the exact result. Within the universal Hamiltonian model, however, the peak spacing distribution is controlled by the joint probability distribution of the spacings between several consecutive levels, which we shall denote as , , etc. As for the single spacing distribution, the exact expressions for the joint spacing distributions are too cumbersome for practical use.

We have, therefore, developed approximate expressions analogous to Wigner’s surmise for the joint level spacing distributions and used these to obtain formulae for the conductance peak spacing distribution. In addition, we have developed new numerical techniques for evaluating the exact expressions which are much more efficient than previously known methods. The approximate analytic expressions are in excellent agreement with the numerics.

To our knowledge there has been no previous investigation of these joint spacing distributions in a physics context Orsay99 (). Number theorists have studied a quantity dubbed the nearest neighbor spacing distribution in connection with the zeros of Riemann’s zeta function ForrOdlyzko (). That distribution is closely related to , and, indeed, the zeta function work, as well as other known results, provide a useful test of the accuracy of the approximate formulae we develop.

The rest of the paper is organized as follows. First, we start with the main results of interest from a QD point of view: Section II gives our analytic results for the peak spacing distribution. Then, in Section III we discuss the RMT quantities on which the peak spacing distribution is based. The approximate expressions for the joint distribution of several consecutive levels are developed here. Section IV presents our exact results for these RMT quantities and compares them with the previously given surmises. Finally, in Section V we summarize and conclude.

## Ii Peak Spacing Distribution

The peaks in the conductance occur at values of the gate voltage for which the number of electrons on the dot changes by one. In the quantum Coulomb blockade regime, the spacing between the and peaks is

 δ2≡EN+1+EN−1−2EN, (2)

where denotes the ground state energy of the QD with electrons. The sequence of peak spacings is sometimes called the addition spectrum of the QD. For later use it is convenient to define the shifted peak spacing

 s=δ2−e2C. (3)

It is instructive to first consider the constant interaction model. For even, the transition involves adding an electron to level , while in the transition an electron is added to level . Thus, for an odd-even-odd transition (for brevity, an even spacing hereafter)

 s=Δ, (4)

where is the spacing between the two levels. Similarly, for an even-odd-even transition (an odd spacing hereafter)

 s=0. (5)

Thus in the constant interaction model, the peak spacing distribution is bimodal: the even spacings are distributed in the same way as the single particle level spacings while the odd spacing distribution is a delta function spike at .

Let us now consider the case of the universal Hamiltonian model, Eq. (1). For even , the state with the lowest single particle levels doubly occupied has total spin . It is therefore an eigenstate of the universal Hamiltonian although, as we shall see, it is not necessarily the ground state. Next, consider the four states obtained by promoting an electron from the highest occupied level to the lowest unoccupied one (the occupations here are stated relative to the constant interaction ground state). The single-particle energy of these states is greater by , and they can be combined into a spin triplet and singlet. We have thus identified five eigenstates of the universal Hamiltonian: the constant interaction ground state with energy , the degenerate triplet states with energy and the singlet state with energy . Proceeding in this way we can construct all the eigenstates of the universal Hamiltonian from the excitations of the constant interaction model.

For , the constant interaction ground state is the true ground state, but for the triplet states become lower in energy. In principle a state with still higher spin might have even lower energy: such states are costlier in terms of single particle energy but have large (negative) exchange contributions to their energy. One might expect that in the paramagnetic limit (the mean level spacing) the single particle energy would dominate, and a high spin state would be unlikely. Indeed for as large as we have verified that, in systems without time-reversal symmetry, approximately of the ground states for even have or . Similarly, and are the dominant spin states for odd. Hence, in the following we will consider only the competition between the two lowest spin states.

Now consider the peak spacing for an even spacing. For simplicity, first assume that and or ; in this case the spacing is given by

 s=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Δ−32JforΔ>2J52J−ΔforΔ<2J. (6)

Note that . The cumulative probability that the spacing is smaller than a certain value is

 F(s)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩0fors<12J∫32J+s52J−sP1(Δ)dΔfor12J

Here denotes the single particle level spacing distribution, for which both the exact result and an extremely accurate approximation—Wigner’s surmise—are known in RMT.mehta () By differentiation we can convert the cumulative distribution to the peak spacing probability distribution

 (8)

Note that there is a sharp jump in the distribution at the left edge of its support , in contrast to the smooth behavior in the constant interaction model. Also is continuous at but kinked in systems with time-reversal symmetry [see Eq. (12) below].

The odd spacing is similar but more tedious to analyze because there are four cases to consider: the QD may have either spin zero or spin one in both initial and final states. We denote the spacing between the level and the level above it and between the level and the level above that . Straightforward analysis then shows that

 s=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩32JforΔ1>2J~{}and~{}Δ2>2JΔ2−12JforΔ1>2J~{}and~{}Δ2<2JΔ1−12JforΔ1<2J~{}and~{}Δ2>2JΔ1+Δ2−52JforΔ1<2J~{}and~{}Δ2<2J. (9)

Note . Provided and (always true in the constant interaction limit of ) all odd peak spacings have the same shifted value independent of the values of the level spacings. However if these conditions are not met the spacing does depend on the values of and . A short calculation reveals

 Πodd(s)=δ(s−32J)∫∞2J∫∞2JP2(Δ1,Δ2)dΔ1dΔ2+C(s). (10)

Here is the joint probability distribution of two consecutive level spacings, for which we will derive an accurate approximation analogous to the Wigner surmise in Section III. The continuum distribution is given by

 C(s)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩0fors<−52Jors>32J∫s+52J0P2(μ,~μ)dvfor−52J

with and .

Eqs. (8), (10), and (11) are the expressions for the Coulomb blockade peak spacing distributions when it is assumed that the QD cannot have a ground state spin greater than one. These expressions apply whether time-reversal symmetry is intact or broken: the only difference in the two cases is that it is necessary to use the level spacing distributions , for the Gaussian orthogonal ensemble (GOE) when time-reversal symmetry is intact and for the Gaussian unitary ensemble (GUE) when it is broken.

Unfortunately, this very simple approach is not enough to describe some important features of the peak spacing distribution for values of . States with have to be included. Table 1 shows all possible spin transitions, the corresponding peak spacing, and the conditions on the level spacings for the states involved to be the ground states. We used the spin selection rule upon the addition of an electron. By a simple inspection of Table 1 we note several new features that appear in this analysis: (i) The even peak spacing distribution now depends not only on and but also on the joint probability distribution of three consecutive level spacings, —an accurate approximation to this quantity is derived in Section III; (ii) even transitions with and odd transitions with are now possible; and (iii) a discontinuity appears in the odd distribution at .

An explicit expression for the peak spacing distribution is also possible in this case. We present it in Appendix A as the calculation is rather lengthy. Here, we simply show the results for the GOE in Figure 1. Transitions and both give the same contribution, so only one is displayed; similarly for and .

A comparison of the peak spacing distribution, calculated in the approximation that allows no spin state higher than , to a numerical simulation ullmo01 () is shown in Figures 2 and 3 for the GOE and GUE, respectively. The numerical calculation involved the diagonalization of random matrices of size . Each random matrix corresponds to a different realization of the single particle levels in Eq. (1), and the peak spacing distribution is inferred by histogramming the calculated peak spacings for these “virtual quantum dots”. As can be seen from the figures, the agreement between the simulation and the analytic formula is excellent.

We have found that it is not necessary to consider higher spins for as large as . The analytic expression offers considerable advantage over the numerical calculation if, for example, it is necessary to compare the universal Hamiltonian model to an experimental distribution with the exchange constant as fitting parameter.

## Iii Level Spacing Distributions: Simple Approximations

In this Section we develop accurate approximations to the joint probability density of consecutive level spacings. These expressions are needed to complete the expressions for the Coulomb blockade peak spacing distribution discussed above.

First it is helpful to recall Wigner’s work mehta (). We are interested in , the spacing between two consecutive levels near the center of the spectrum of a large random matrix. The exact result, derived by a lengthy and intricate calculation mehta (), involves an infinite product of eigenvalues of prolate spheroidal functions which are difficult to evaluate numerically. It is therefore of limited practical value.

Now, following Wigner, consider a real symmetric Gaussian random matrix. Such a matrix has two eigenvalues, and an elementary calculation shows that the spacing between them is distributed according to

 Pgoe1WS(Δ)=π2Δexp(−π4Δ2). (12)

The scale is chosen so that the mean level spacing is unity,

 ⟨Δ⟩=∫∞0ΔPgoe1WS(Δ)dΔ=1. (13)

Eq. (12) is the celebrated Wigner surmise. For the GUE, the surmise is similarly derived by considering a complex hermitian Gaussian random matrix, yielding

 Pgue1WS(Δ)=32π2Δ2exp(−4πΔ2). (14)

The Wigner surmise is known to be an excellent approximation to the level spacing distribution for the orthogonal ensemble where it has been checked against numerical calculations and rigorous bounds (see page 157 of Ref. mehta, ). For the unitary case, to our knowledge, the surmise has not been subject to comparably rigorous scrutiny, but it is generally believed to be very accurate guhr () (see also Fig. 8 below where the unitary Wigner surmise is compared to our numerical calculation of the level spacing distribution).

Let us turn to the joint probability distribution of two consecutive level spacings, . The exact result may be derived using the theory of Fredholm or Toeplitz determinants, as in Section IV.4 below. However, a more useful approximation may be derived by considering, in the spirit of the Wigner surmise, a real symmetric Gaussian random matrix (for the GOE). Using the standard joint probability distribution of the levels of a random matrix, we obtain

 Pgoe2;3×3(Δ1,Δ2) = 4√23πa52Δ1Δ2(Δ1+Δ2) × exp(−2a3[Δ21+Δ22+Δ1Δ2])

with . The distribution has been normalized and the energy scale chosen so that . We conjecture that Eq. (III) is an accurate approximation to the true joint spacing distribution for large GOE matrices.

A similar result for the GUE may be derived by considering complex hermitian Gaussian random matrices. We obtain

 Pgue2;3×3(Δ1,Δ2) = 4b4π√3Δ21Δ22(Δ1+Δ2)2 × exp(−2b3[Δ21+Δ22+Δ1Δ2])

with , which we conjecture is a good approximation to the true joint spacing distribution for the GUE.

There is no a priori justification for either Wigner’s surmise or for our analogous conjectures. To test our conjectures we subjected them to many checks, some of which we describe here. For definiteness we focus here on the unitary ensemble; but we have done similar tests in the orthogonal case.

First, we note the relationship

 Pgue1(Δ1)=∫∞0Pgue2(Δ1,Δ2)dΔ2 (17)

between the exact single level spacing distribution and the exact joint spacing distribution . Performing the indicated integral using our surmise Eq. (III), we obtain

 Pgue1;3×3(Δ1) = ∫∞0Pgue2;3×3(Δ1,Δ2)dΔ2 = 9b322√2πΔ21exp(−b2Δ21)h(Δ1√b6)

where

 h(x)=−x3e−x2+32xe−x2+√π(34−x2+x4)Erfc(x) (19)

and is the complementary error function.

Fig. 4 shows a comparison between the Wigner surmise and the approximation Eq. (III). The agreement is excellent over the entire region where the distribution has substantial weight and also for small spacings. In the tails, the relative error is larger but, of course, very small on an absolute scale. To be precise of the weight of lies in the range . Over this range, and disagree by less than .

Another quantity related to is the nearest neighbor spacing introduced by Forrester and Odlyzko in connection with studies of the zeros of the zeta function:ForrOdlyzko () the nearest neighbor is the closer of the level just above and the level just below a given level. The distribution of the nearest neighbor spacing is, therefore,

 Pguenn(Δ) = ∫∞0dΔ1∫Δ10δ(Δ−Δ2)Pgue2(Δ1,Δ2)dΔ2 + ∫∞0dΔ1∫∞Δ1δ(Δ−Δ1)Pgue2(Δ1,Δ2)dΔ2.

If we substitute the approximation (III) into Eq. (LABEL:eq:nearestneighbour), we obtain the approximation to the nearest neighbor distribution, where

 gnn(x) = 272√2πx2e−12x2[18x+84x3 (21) +e9x2π12(3−4x2+4x4)Erfc(3x)]

Fig. 5 shows a comparison between this approximation to and the exact numerical computation in Section IV. Again the agreement is very good except in the tails which have negligible weight. We have also confirmed that these two results are in good agreement with the exact numerical computation of using the method of Forrester and Odlyzko.ForrOdlyzko ()

The covariance of consecutive spacings is also known mehta (). The exact result for is in the orthogonal ensemble and in the unitary ensemble; within our approximation for they are found to be and , respectively.

As a final check we can compare Eq. (III) directly to the exact numerical computation of discussed in Section IV. Again the agreement is reasonable (see Fig. 9 below).

Next, we turn to the joint probability distribution of three consecutive spacings. To this end, we consider a real symmetric Gaussian random matrix. From the standard expression for the distribution of eigenvalues for a random matrix mehta (), we obtain the normalized distribution

 Pgoe3;4×4(Δ1,Δ2,Δ3)=8√πΔ1Δ2Δ3(Δ1+Δ2)(Δ2+Δ3)(Δ1+Δ2+Δ3)exp[−14{2(Δ1+Δ2)2+2(Δ2+Δ3)2+(Δ1+Δ3)2}].

Using this distribution, we find the average of the middle spacing, . A difficulty that now arises is that due to end effects in our extremely finite sized matrix, . We resolve this problem by the standard procedure of “unfolding” BohigasLesHouches (), introducing the rescaled distribution

 Pgoe3;4×4(Δ1,Δ2,Δ3)=f′2fPgoe3;4×4(f′Δ1,fΔ2,f′Δ3). (22)

We conjecture that is an accurate approximation to the true joint spacing distribution for the GOE.

We have subjected this conjecture to tests similar to those applied to the distribution above. For example, we can integrate out the end spacings and to obtain an approximation to the single level spacing distribution or we can integrate out the lower ones and to obtain a second approximation. The two approximations are compared to each other and to the nearly exact Wigner surmise in Fig. 6 and are found to be accurate outside the tails. On the basis of tests such as these we conclude that our conjectured distribution Eq. (22) is sufficiently accurate for our purposes.

The corresponding distribution for the unitary ensemble can be similarly derived by consideration of a complex hermitian matrix. We obtain the normalized distribution

 Pgue3;4×4(Δ1,Δ2,Δ3)=κΔ21Δ22Δ23(Δ1+Δ2)2(Δ2+Δ3)2(Δ1+Δ2+Δ3)2exp[−14{2(Δ1+Δ2)2+2(Δ2+Δ3)2+(Δ1+Δ3)2}], (23)

where . In this case we find whereas . Unfolding the distribution (23), we obtain

 Pgue3;4×4(Δ1,Δ2,Δ3)=h′2hPgue3;4×4(h′Δ1,hΔ2,h′Δ3), (24)

which we find to be an accurate approximation to the true joint spacing distribution for the GUE.

In summary, we have obtained simple expressions for the joint spacing distribution of two consecutive spacings, Eqs. (III) and (III), and for three consecutive spacings, Eqs. (22) and (24), for the orthogonal and unitary ensembles, respectively. Together with the Wigner surmise, Eqs. (12) and (14), these distributions are the random matrix theory tools needed to compute the Coulomb blockade peak spacing distribution in the paramagnetic limit. Using these tools, we obtained the peak spacing distributions already given in Section II.

## Iv Level Spacing Distributions: Exact Results

In this Section we derive a number of exact results in random matrix theory. These results allow us to certify the accuracy of the more useful simple approximations developed in the previous section. Some of the results, notably the new exact calculation of and the analysis of its asymptotic behavior, are also of intrinsic interest in random matrix theory.

The problem under consideration is the following: We are given a set of energy levels and , their joint probability distribution. According to random matrix theory the energy level distribution is

 P(x1,…,xn)=N∏1≤i

where the exponent , , or corresponds respectively to the orthogonal, unitary, and symplectic ensembles of random matrix theory. is a normalization constant, and is an energy scale set by the mean spacing between levels.

¿From the joint probability distribution we wish to calculate quantities such as the consecutive level spacing distribution

 P1(t) ≡ K1n(n−1)∫outdx3…∫outdxn (26) P(x1→−t2,x2→t2,x3,…,xn),

where

 ∫outdx≡∫−t/2−∞dx+∫∞t/2dx. (27)

This is the probability that one level is at , another at , and none of the others are in between. We may place any of the levels at and any of the remaining levels at ; this is the origin of the combinatoric prefactor in (26). The normalization ensures that . Generally we wish to calculate the level spacing distribution in the limit . In addition to the level spacing, similar quantities of interest here are a more primitive quantity that we shall call the spacing determinant (defined below), the nearest neighbor spacing, and the joint probability distribution of two consecutive spacings.

The essential difficulty in random matrix calculations is the highly correlated nature of the probability distribution Eq. (25). Level spacing calculations are particularly difficult due to the piecewise continuous nature of the integration domain. Nonetheless, exact formal expressions for these quantities have been derived. For example, the spacing determinant can be expressed as the Fredholm determinant of a particular integral operator , the limit of a certain Toeplitz determinant, and the solution to an ordinary non-linear second-order differential equation. Explicit evaluation of the spacing determinant is then traditionally carried out numerically from the Fredholm determinant or the non-linear differential equation. In Section IV.1 we find that direct evaluation of the Toeplitz determinant for large but finite provides a third efficient numerical method of calculating the spacing determinant.

In Sections IV.2, IV.3, and IV.4 we derive new expressions for the consecutive level spacing , the nearest neighbor spacing , and the joint distribution of two consecutive spacings , expressing each of these quantities as a Toeplitz determinant. These quantities may then be numerically evaluated by direct calculation of the determinants. By contrast, in the conventional approach they are expressed in terms of the eigenvalues and eigenfunctions of the kernel , which are harder to handle numerically. The Toeplitz representation also makes it easy to obtain the large asymptotic behavior: the previously unknown asymptotic behavior of is given in Section IV.5. For simplicity, we concentrate upon the unitary ensemble, ; the other ensembles will be studied in future work.

Toeplitz determinants arise in many contexts, among them signal processing numerical () and statistical mechanics fisher (). Consequently, a great deal is known or conjectured about their asymptotic behavior and efficient numerical methods exist for their evaluation; that is the chief virtue of the Toeplitz expressions for , , and derived here.

### iv.1 Level Spacing Determinant

The level spacing determinant is the probability that a band of width is entirely void of energy levels,

 E(t)≡∫outdx1…∫outdxnP(x1,…,xn), (28)

where is defined in Eq. (27). Evidently , whereas . By straightforward differentiationfootnote1 ()

 F(t)≡−dEdt =n∫outdx2…∫outdxnP(x1→−t2,x2,…,xn).

is the probability that one level is at the edge of the band while the others lies outside it. Evidently, while where is the mean density of levels, or the one point correlation, defined as

 R1(x)=n∫∞−∞dx2…∫∞−∞dxnP(x1→x,x2…xn). (30)

Differentiating once again, we find that the consecutive level spacing distribution, Eq. (26), is given by

 P1(t)=K1d2Edt2. (31)

In deriving this equation we have assumed that the joint probability distribution vanishes if two levels coincide, and have also neglected a bulk term in which the derivative acts on the integrand in (IV.1). The error is found to vanish as . Eq. (31) shows that in principle is determined by . In practice, is usually computed numerically, and hence Eq. (31) is not the best way to determine . However, we may use this identity to show that the normalization constant is .

It is now useful to briefly recount how the exact formal expressions for are obtained. It is more convenient for this purpose to use Dyson’s circular unitary ensemble rather than the Gaussian unitary ensemble [Eq. (25) with ]. The circular unitary ensemble describes angles distributed around the unit circle according to the normalized distribution

 P(θ1,θ2,…,θn)=1n!1(2π)n∏1≤r

It is well known that (up to an irrelevant scale factor that determines the mean level spacing) the local statistical properties of these angles are identical to those of energy levels governed by the GUE.

In the circular ensemble, is the probability that all angles lie outside the arc . To compute it is helpful to rewrite

 P(θ1,…,θn)=1n!∑P,Q(−1)P(−1)Q ×φP(1)(θ1)φ∗Q(1)(θ1)…φP(n)(θn)φ∗Q(n)(θn).

Here , . and represent permutations of the integers and is the parity of the permutation. This representation follows from (32) by means of the Vandermonde identity

 ∏1≤i

where is the matrix with elements .

If we define

 grs≡∫2π−t/2t/2dθφr(θ)φ∗s(θ)=∫2π−t/2t/2dθ2πexp[i(r−s)θ], (35)

it follows

 E(t) = 1n!∑PQ(−1)P(−1)QgP(1),Q(1)…gP(n),Q(n) (36) = 1n!∑PR(−1)RgP(1),R[P(1)]…gP(n),R[P(n)] = 1n!∑PR(−1)Rg1,R(1)…gn,R(n).

To obtain the second line, we have introduced and used . The third line follows from a simple rearrangement of the terms in the summand. The sum is now trivial, and we find

 E(t)=detg (37)

where is an matrix with matrix elements given by (35).

The mean spacing between angles is (show by computing the one point correlation ). Hence we introduce ,

 Δ≡n2πt, (38)

as a measure of the band-width in units of the mean level spacing. Making this change of variables, we may finally write

 E(Δ)=limn→∞detg (39)

with given by Eqs. (35) and (38).

An Toeplitz matrix has elements

 Trs=∫2π0dθ2πf(θ)exp[i(r−s)θ] (40)

that are controlled by a single function . Comparing Eqs. (35) and (40), we conclude that the matrix is of the Toeplitz form with for and zero otherwise. Hence Eq. (39) is the representation of the level spacing determinant in terms of a Toeplitz determinant.

To motivate the Fredholm representation, consider the eigenvectors of ,

 n∑s=1grsψ(α)s=λ(α)ψ(α)r. (41)

Here labels distinct eigenvectors and their eigenvalues . In the limit, we set , , and . Explicitly evaluating in (35) and taking the limit, we obtain the integral eigenvalue equation

 ∫10dyK(x,y)ψ(α)(y)=λ(α)ψ(α)(x) (42)

with kernel

 K(x,y)=δ(x−y)−1π(x−y)sin[πΔ(x−y)]. (43)

Since the determinant of a matrix is the product of its eigenvalues, we may write

 E(Δ)=∞∏α=1λ(α)=detK. (44)

This is the representation of as a Fredholm determinant, a virtue of which is that has been explicitly taken.

Evaluation of the eigenvalues is facilitated by a remarkable connection between the integral kernel and the prolate spheroidal functions of classical mathematical physics mehta (). It is found that the prolate spheroidal differential operator and the kernel commute and therefore have the same eigenfunctions. In practice it is easier to determine the eigenfunctions of and then to compute by applying to these eigenfunctions.

Finally, we note that a third expression for was derived by Jimbo et al., expressing it as the solution to a second-order non-linear differential equation jimbo ().

Each of these representations has proved useful in the past. The asymptotic behavior of for large [and hence of via Eq. (31)] was derived by des Cloizeaux and Mehta by asymptotic analysis of the differential operator and independently by Dyson by means of an ingenious application of inverse scattering theory to the analysis of the kernel mehta (). It can also be obtained from the Toeplitz representation using an asymptotic formula due to Widom widom (); mehta ().

The first numerical computation of was based on evaluation of the Fredholm determinant, taking advantage of the connection to prolate spheroidal functions mehta (). Numerical solution of the non-linear differential equation was subsequently found to be more efficient jimbo ().

Here we find that direct evaluation of the Toeplitz determinant Eq. (39) for large (but finite) also provides an accurate way to calculate . Fig. 7 shows a plot of calculated in this way with ; the inset shows the convergence of as a function of for . The results for and differ only in the sixth significant figure. The convergence is slower for larger , but it is clear that for all values of for which the distribution has significant weight, and to the accuracy needed for applications to quantum dots, Toeplitz determinants with equal to a few hundred should suffice. Fig. 7 is the main new result of this subsection.

### iv.2 Consecutive Level Spacing

In this sub-section we express the consecutive level spacing (sometimes simply called the level spacing) as a Toeplitz determinant. We then show that this expression can be used to compute with good precision with modest numerical effort.

For simplicity let us suppose that there are angles on the circle. is the probability that one of these equals , another equals and the others all lie outside the range . Thus we must consider

 P(θ1,…,θn,θ