Random Matrix Theory

J. P. Keating
December 19, 2023

1
Contents
1 Overview 5
1.1 Random rotations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Sample covariance matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Systems of linear equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Complex quantum systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 Stability of high-dimensional dynamical systems . . . . . . . . . . . . . . . . . . . . . . 7
1.6 Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7 Complex networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.8 Machine learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.9 Connections with other areas of mathematics . . . . . . . . . . . . . . . . . . . . . . . . 8
1.10 Further reading and remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Examples of random matrix ensembles 10

2.1 Wishart random matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Wigner random matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Gaussian orthogonal random matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Gaussian unitary random matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 Random unitary and orthogonal matrices . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 The semicircle law for Wigner random matrices 13

3.1 Moments of the semicircle law and Catalan numbers . . . . . . . . . . . . . . . . . . . . 16
3.2 Spectral moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 The circular law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4 Gaussian random matrices 24

4.1 Wick’s theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2 The genus expansion for Gaussian random matrices . . . . . . . . . . . . . . . . . . . . . 25

5 The Marchenko-Pastur law for Wishart random matrices 27

5.1 The Marchenko-Pastur Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

6 Stieltjes and R-transforms 33

6.1 The Stieltjes transform and the semicircle law . . . . . . . . . . . . . . . . . . . . . . . . 33
6.2 The R-transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

7 Eigenvalue statistics for GOE, GUE, CUE 38

7.1 Poisson point processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
7.2 Joint eigenvalue density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
7.3 The method of orthogonal polynomials for the GUE . . . . . . . . . . . . . . . . . . . . 41
7.4 Correlation functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
7.5 Scaling at the edge of the spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
7.6 Counting statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
7.7 Other Ensembles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
7.8 Universality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
7.9 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

8 Dyson Brownian motion 57

2
9 Some connections with other areas of mathematics 63
9.1 Longest increasing subsequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
9.2 Stochastic growth models, random tilings, and random paths . . . . . . . . . . . . . . . 64
9.3 Zeros of the Riemann zeta-function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3
Preface
These notes accompany the Part C course on Random Matrix Theory, C7.7. It is a pleasure to thank
Michael Curran, Johannes Forkel, Bhargavi Jonnadula, Isao Sauzedde, and Mo Dick Wong for their
help in preparing them. I would be grateful for comments, suggestions, and for reports of typos and
errors.

J. P. Keating
December 2021

4
1 Overview
Matrices play a central role in mathematics and its applications. Just as one considers the properties
of probability measures defined on a space of real or complex variables, it is therefore also natural to
seek to associate probability measures with spaces of matrices. We shall be interested in the properties
of such random matrices, and in applications of Random Matrix Theory to predict what one should
expect in ‘typical’ situations involving linear algebra.
Historically, questions in Random Matrix Theory (RMT) have arisen in a number of different
mathematical contexts. We shall start by reviewing briefly a selection of these, by way of motivation
and a gentle introduction to some of the key themes that will be explored later in more detail.

1.1 Random rotations

Arguably the first substantial RMT calculation appeared in work of Adolf Hurwitz concerning orthogonal
transformations. Rotations of objects in three dimensions can be parametrized by the Euler angles.
One can therefore consider random rotations in three dimensions by defining a probability measure on
these angles. In the absence of any preferred directions, it is natural to use a probability measure that
is itself invariant under all rotations.
Another way of phrasing this is in terms of orthogonal matrices. Consider r ∈ R3 . Rotating r
corresponds to multiplying it by a 3 × 3 orthogonal matrix1 , i.e. a matrix O satisfying OOT = I, where
I denotes the 3 × 3 identity matrix and, here and hereafter, M T denotes the transpose of a matrix
M . Random rotations can therefore be thought of in terms of a probability measure on the space of
orthogonal matrices. A rotationally invariant probability measure corresponds to a probability measure
that is itself invariant under all orthogonal transformations.
It is an obvious question then how this extends to orthogonal matrices of dimension n; that is, how
one can construct an invariant probability measure on the space of n × n orthogonal matrices. Moreover,
orthogonal transformations form a compact group, so one can ask how this story generalizes to other
related groups, for example to the group of unitary transformations in n dimensions, represented by
n × n unitary matrices U , i.e. matrices satisfying U U † = I, where M † denotes the complex conjugate
of the transpose, M T , of a matrix M . Unitary transformations play a fundamental role in quantum
mechanics, particularly in quantum information theory, and random unitary matrices have important
applications to modelling ‘typical’ quantum transformations.
Hurwitz was the first to investigate invariant probability measures on the orthogonal and unitary
groups. See [6] for an overview of his work and subsequent developments, and [13] for an introduction
to the general theory in wider contexts.

1.2 Sample covariance matrices

Denote by Np (µ, Σ) the multivariate normal distribution of a p-dimensional random vector X =
(X (1) , X (2) , . . . , X (p) )T , with

µ = EX = (EX (1) , EX (2) , . . . , EX (p) )T (1.1)

and
Σij = E(X (i) − µi )(X (j) − µj ) (1.2)
for 1 ≤ i, j ≤ p.
Now let X1 , X2 , . . . , Xn be n i.i.d. p × 1 random vectors with distribution Np (0, Σ). In many
situations in data science, mathematical finance, and statistics one encounters problems of this kind
1 More generally, orthogonal matrices generate linear isometries; the orthogonality of O preserves the dot product

between any two vectors

5
where Σ is not known and where one wants to find it. From the data given, the sample covariance
matrix
n
1X
σn = Xk XkT (1.3)
n
k=1

is a natural estimator for Σ.

(j)
Note that defining the n × p matrix X̃ij = Xi , one also has that

1
σn = X̃ X̃ T . (1.4)
n
Random matrices having this form are known as Wishart random matrices, after John Wishart,
who introduced them in 1928 [20].

1.3 Systems of linear equations

One of the basic problems of numerical linear algebra involves solving the system of linear equations

Ax = b (1.5)

for the n-dimensional vector x, given an m × n matrix A and an m-dimensional vector b, when m and
n are large. One important consideration is: if b is not specified precisely, how does this imprecision
affect the accuracy with which x can be determined? The condition number is the maximum ratio
of the relative error in x to the relative error in b. A natural way to measure the size of these errors
equates the condition number to the ratio of the largest and smallest singular values of A, or in the
case when A is normal2 , the ratio of the maximal and minimal absolute values of the eigenvalues of A.
In 1947 John von Neumann and Herman Goldstine [15] asked what the condition number would be
for a ‘typical’ large matrix and initiated the study of condition numbers of random matrices, taking the
elements of A to be i.i.d. normal3 random variables. This subject has a long and interesting history,
which is beautifully reviewed in [7]. In particular, the distribution of values taken by the condition
number is related to the distribution of the largest and smallest eigenvalues of the associated random
matrices.

1.4 Complex quantum systems

In quantum mechanics, the allowed values of the energy in a closed system – the energy levels – are, in
general, the eigenvalues of a complex Hermitian matrix H = H † , the quantum Hamiltonian. In many
settings this matrix is in fact real and symmetric, so H = H T . It was suggested by Eugene Wigner in
the 1950s that in complex quantum systems the independent entries in the matrix H (Hij with i ≥ j,
say) should be modelled as random variables. The question then is: how are the eigenvalues of random
complex Hermitian or random real symmetric matrices distributed? It turns out that for matrices
whose entries are independent random variables, the eigenvalues are strongly correlated in a distinctive
and mathematically interesting way. These characteristic correlations are indeed seen in the energy
level statistics of typical quantum systems, ranging from atomic nuclei to superconducting systems.
The fact that large perfectly ordered systems conduct electricity while disordered systems do not is
because in the latter case the quantum Hamiltonians behave like random matrices.
In fact this is not an exclusively quantum phenomenon: it is observed in all wave theories, including
acoustics, optics, elasticity, electromagnetism, etc. Moreover, it is observed on the widest range of
scales, ranging from the description of sub-nuclear physics in terms of quantum chromodynamics to the
structure of the cosmic microwave background.
2 i.e. if AA† = A† A.
3 i.e. Gaussian – not to be confused with the use of ‘normal’ in the sense of the previous footnote.

6
1.5 Stability of high-dimensional dynamical systems
Let x(t) be an n-dimensional vector satisfying

dx
= −λIx (1.6)
dt
where λ is a positive constant and I is the n × n identity matrix. This system has a stable fixed point
at x = 0. Clearly in (1.6) the different components of x are uncoupled (because I is diagonal). The
question is: if one now introduces random coupling between these components, does this make the fixed
point more or less stable? Specifically, what is the stability of the system
dx
= −λIx + Ax (1.7)
dt
where A is a random matrix?
This question was first raised by Robert May (Lord May of Oxford) in 1972 in a famous paper [11]
on population dynamics. In that context, it is called the May model. Then the components of x
represent the differences of the populations of various species from some equilibrium values. If the
species do not interact, then it is assumed that the equilibrium values are stable. It had long been
assumed that adding random interactions would make the equilibrium populations more stable when n
is taken to be large. May asked whether this was indeed the case and analysed the situation using the
simple model system (1.7).
Obviously the question of stability is related to knowing the expected size of the eigenvalues of A
when n is large.
This question arises in many other contexts as well, for example in the dynamics of neural models
and on complex networks.

1.6 Principal Component Analysis

Let X1 , X2 , . . . , Xn be p × 1 vectors, where p is to be treated as large with respect to the number of
vectors, n. Let Y be another p × 1 vector. Which of X1 , X2 , . . . , Xn is closest to Y ? Obviously one
can check each of the distances ||Y − Xi ||2 , but doing so is costly. This is an important problem in
high-dimensional data analysis.
The idea of Principal Component Analysis is to project onto a subspace in which the vectors
X1 , X2 , . . . , Xn show maximum variability. So one could look for the unit vector u ∈ Rp such that
uT X1 , uT X2 , . . . , uT Xn have maximum variability. For simplicity, let us assume that X1 + X2 + · · · +
Xn = 0. Then the variance of the set uT X1 , uT X2 , . . . , uT Xn is
 
Xn Xn
1 1
(uT Xj )2 = uT  Xj XjT  u = uT σn u. (1.8)
n j=1 n j=1

The maximum value uT σn u can take is the largest eigenvalue4 of σn , and this occurs when u is the
corresponding eigenvector.
Clearly this can be extended so that one can project onto the subspace spanned by the eigenvectors
corresponding to the m largest eigenvalues. How large should m be taken? The idea is to compare the
eigenvalues to those of the random matrix Σ defined by (1.2), taken to be a null model, and to keep
those that differ significantly and so represent non-random features of the data. One therefore needs to
know how the eigenvalues of a random Wishart matrix are distributed.
Returning to the problem of comparing Y with X1 , X2 , . . . , Xn , it is obviously natural to seek to
do this in the subspace where the vectors X1 , X2 , . . . , Xn exhibit maximum variability, because this
is where their differences are largest. This can be thought of as enacting ‘feature selection’ in many
4 The eigenvalues of σn are real because the matrix is clearly real and symmetric.

7
applications. The reduction in the dimension of the space where the comparison is made increases the
efficiency of the search considerably.
For an overview of this and similar applications (e.g. to mathematical finance, image analysis, etc)
see [8] and [4].

1.7 Complex networks

Many problems in data science and mathematics are related to properties of networks. These can be
modelled by graphs, that is by ordered pairs G = (V, E) comprising a set of vertices V and edges E.
The structure of the graph may be represented in terms of the adjacency matrix A. For a simple graph
this is a square matrix of dimension |V | such that Aij is equal to 1 when there is an edge from vertex i
to vertex j, and zero otherwise5 .
The statistical properties of complex graphs are then modelled by random square matrices where
the elements are taken from {0, 1}, each with probability p, in some cases with constraints (e.g. in the
case of regular graphs on the total number of ones in each row and column).
In the physics literature, these network models, where often the non-zero entries of the adjacency
matrix are augmented with a ‘hopping probability’, are usually called tight binding models. One well
known example is the Anderson tight binding model used to explore quantum localization and its
implications for electrical conductivity.

1.8 Machine learning

In machine learning one frequently wishes to minimise functions, which are often highly complex, in
an extremely high-dimensional space. One technique is to use stochastic gradient descent. This raises
the important question as to how easy it is to explore effectively random surfaces in high dimensional
problems; are there many local minima and saddles where one can get stuck? And at a saddle, how
many downward directions are there typically? This is the subject of intense study. One way of
modelling the problem is to take the Hessian at a point on the surface to be a random symmetric matrix.
Understanding the structure of the critical points then reduces to understanding the distribution of the
eigenvalues of this random matrix, for example how many eigenvalues are expected to be positive, and
how many negative. See, for example, [5].

1.9 Connections with other areas of mathematics

In the above examples it is hopefully clear that matrices play an important role, and that it is reasonable
that one might wish to use random matrices as a statistical model. What is more surprising is that the
mathematical structures one finds in random matrix theory are more general than one might expect, in
that they arise in contexts with no obvious connection with linear algebra.
One example concerns the length of the longest increasing subsequence in random permutations. Let
Sn be the group of permutations of 1, 2, . . . , n. If π ∈ Sn , π(i1 ), . . . , π(ik ) is an increasing subsequence
in π if i1 < i2 < · · · < ik and π(i1 ) < π(i2 ) < · · · < π(ik ). Let ln (π) be the length of the longest
increasing subsequence. For example, if n = 5 and π is the permutation 5 1 3 2 4, then the longest
increasing subsequences are 1 2 4 and 1 3 4, and ln (π) = 3. Equip Sn with uniform distribution,

#{π ∈ Sn : ln ≤ m}
P(ln ≤ m) = . (1.9)
n!
What is the asymptotics of this distribution as n → ∞? Remarkably, it is the same as the asymptotics
of the distribution of the largest eigenvalue of a random complex Hermitian n × n matrix, despite the
fact that there is no matrix in the problem.
5 This matrix is related to the discrete Laplacian acting on the vertex set.

8
 A second example concerns the Riemann zeta-function
X∞
1
ζ(s) = s
. (1.10)
n=1
n

The sum converges in Res > 1, and the zeta function has an analytic continuation to all s (except for
s = 1, where it has a pole). The Riemann Hypothesis asserts that all of the zeros of ζ(s) with non-zero
imaginary parts lie on the line Res = 1/2. How are the zeros distributed along the line? Like random
points, or are their positions correlated? It turns out that they are correlated, and all of the evidence
we have suggests that the correlations between them are the same as those between the eigenvalues of
large complex Hermitian random matrices, despite the fact that again there is no obvious matrix in the
problem.
Moreover, the same correlations have been found in the distances between parked cars in London,
and the gaps between bus arrivals in the Mexican city of Cuernavaca.
Finally, many of the calculations in random matrix theory have deep connections with ideas and
results in other areas of mathematics, including combinatorics and topology.

1.10 Further reading and remarks

In addition to these notes, students may wish to consult the following for further details and extensions
of the material covered.

• M. Potters & J.-P. Bouchad, A First Course in Random Matrix Theory for Physicists, Engineers
and Data Scientists (Cambridge University Press)
• G. Livan, M. Novaes & P. Vivo, Introduction to Random Matrices (Springer Briefs in Mathematical
Physics)
• E.S. Meckes, The Random Matrix Theory of the Classical Compact Groups (Cambridge University
Press)
• G.W. Anderson, A. Guionnet & O. Zeitouni, An Introduction to Random Matrices (Cambridge
Studies in Advanced Mathematics)

• M.L. Mehta, Random Matrices (Elsevier, Pure and Applied Mathematics Series)
• G. Akemann, J. Baik & P. Di Francesco, The Oxford Handbook of Random Matrix Theory
(Oxford University Press)
• Z. Bai & J.W. Silverstein, Spectral Analysis of Large Dimensional Random Matrices (Springer).

• T. Tao, Topics in Random Matrix Theory, available online from

[Link]

It is worth remarking that random matrix theory is mathematically an extremely broad subject.
Papers and books on the subject vary considerably in style and intended readership, from analysis
and probability theory to engineering, theoretical physics, and data science. They therefore assume
different background knowledge and differ in levels of abstraction and notions of rigour. These notes
will attempt to chart a middle course, focusing on the main ideas and using a mixture of techniques,
hopefully to give a sense of the various viewpoints that have influenced the development of the subject.

9
2 Examples of random matrix ensembles
Motivated by the examples described in the previous section, we start by defining explicitly the random
matrix ensembles that will be studied in this course6 . Here an ensemble is a space of matrices endowed
with a probability measure.
We start with some general notation. We say that an n × n matrix M is symmetric if Mij is real
∗
and Mij = Mji , and that it is Hermitian if Mij is complex and Mij = Mji , where ∗ denotes complex
conjugation. We denote Sn = {n × n symmetric matrices} and Hn = {n × n Hermitian matrices}.
The Lebesgue measure on Sn is by definition the product of the Lebesgue measures on the linearly
independent entries of M :
Y Yn
dM = dMij dMii . (2.1)
1≤i<j≤n i=1

Similarly on Hn :
Y n
Y
dM = d(ReMij )d(ImMij ) dMii . (2.2)
1≤i<j≤n i=1

2.1 Wishart random matrices

Let p ≤ n. A real Wishart matrix with covariance matrix Σ ∈ Sp is a p-dimensional random symmetric
positive definite matrix M of the form
M = XX T (2.3)
where X is an p × n matrix with columns that are zero-mean, independent, Gaussian vectors with
covariance matrix Σ.
The probability density function of M was computed by Wishart in 1928 to be7
 
1 (n−p−1)/2 1 −1
f (M ) = np/2 (det M ) exp − Tr(Σ M ) (2.4)
2 Γp ( n2 )(det Σ)n/2 2

with respect to Lebesgue measure on the cone of symmetric positive definite matrices. Here

Yp  
n n j−1
Γp ( ) = π p(p−1)/4 Γ − . (2.5)
2 j=1
2 2

2.2 Wigner random matrices

A Wigner matrix is a random matrix M ∈ Sn or M ∈ Hn such that (i) Mij with i < j are i.i.d. random
variables (real or complex) with EMij = 0 and E|Mij |2 = 1 , (ii) Mii are i.i.d. real random variables
(with a distribution not necessarily the same as for the off-diagonal elements) with EMii = 0 and
EMii2 < ∞. In some circumstances one rescales the matrix elements by a power of n. This then rescales
the eigenvalues. For example, we shall see that some important results emerge when one rescales the
matrix elements, and hence the eigenvalues, by n−1/2 , and others emerge when one scales them by n1/2 .

2.3 Gaussian orthogonal random matrices

The Gaussian Orthogonal Ensemble, or GOE, is the particular class of Wigner random matrices for
which the matrix elements are independent real normal random variables. Specifically, (i) Mij with
i < j are i.i.d. random variables with distribution N (0, 1) (i.e. Gaussian with zero mean and variance
6 There are many important ensembles that we will not have the time to investigate. Those defined here are a

representative set and are amongst the most widely studied.

7 The precise form of the density is rarely used.

10
1), and (ii) Mii are i.i.d. real random variables with distribution N (0, 2) (i.e. Gaussian with zero mean
and variance 2).
Hence, the probability measure on the entries of M ∈ Sn is
n
Y Y
1 1 2 1 1 2 1 1 2
√ e− 4 Mii √ e− 2 Mij dM = n/2 n(n+1)/4
e− 4 TrM dM. (2.6)
i=1
4π 1≤i<j≤n
2π 2 (2π)

The following is one of the most important properties of the GOE. Let O be a non-random n × n
orthogonal matrix (i.e. OOT = I) and let M be an n × n GOE matrix. Then the distribution of M is
the same as that of OM OT ; that is, the GOE is invariant under conjugation by all orthogonal matrices8 .
To prove this, note first that

Tr[OM OT ]2 = Tr[OM 2 OT ] = Tr[M 2 OT O] = TrM 2 ; (2.7)

second, it follows from calculating the Jacobian of the transformation that

d[OM OT ] = dM. (2.8)

To see this directly, note that the Hilbert-Schmidt norm of M is precisely TrM 2 , which we have just
seen is invariant under conjugation by O. The map M 7→ OM OT is therefore an isometry and so its
Jacobian determinant is 1. Alternatively one can compute the Jacobian by hand and confirm this by
explicit calculation.
The only ensemble of symmetric random matrices whose entries are independent (up to the symmetry
constraint) and whose distribution is invariant under conjugation by all orthogonal matrices is the
GOE – so the orthogonal invariance is a special property of the Gaussian nature of the matrices. One
can easily construct matrix ensembles that are invariant under conjugation by all orthogonal matrices
but in which the matrix elements are not statistically independent.

2.4 Gaussian unitary random matrices

The Gaussian Unitary Ensemble, or GUE, is the particular class of Wigner random matrices for which
the matrix elements are independent complex normal random variables. Specifically, (i) Mij with i < j
are i.i.d. random variables with real and imaginary parts that are independently Gaussian, each with
zero mean and variance 1/2, and (ii) Mii are i.i.d. real random variables with distribution N (0, 1)
(i.e. Gaussian with zero mean and variance 1).
In this case the probability measure on the entries of M ∈ Hn is
n
Y Y
1 1 2 1 2 1 2 1 1 2
√ e− 2 Mii √ e−(ReMij ) √ e−(ImMij ) dM = n/2 n2 /2 e− 2 TrM dM (2.9)
i=1
2π π π 2 π
1≤i<j≤n

and it is invariant under conjugation by all non-random unitary matrices U ; i.e. for U an n × n unitary
matrix (i.e. U U † = I), M and U M U † have the same distribution9 . As in the case of the GOE, to prove
this, note that
Tr[U M U † ]2 = Tr[U M 2 U † ] = Tr[M 2 U † U ] = TrM 2 , (2.10)
and that it follows from calculating the Jacobian of the transformation that

d[U M U † ] = dM. (2.11)

Again, this may be viewed as a consequence of the fact that TrM 2 is the Hilbert-Schmidt norm of M ,
and so conjugation by a unitary matrix is an isometry.
8 This is the reason for the name Gaussian Orthogonal Ensemble.
9 This is the reason for the name Gaussian Unitary Ensemble.

11
 The only ensemble of Hermitian random matrices whose entries are independent (up to the symmetry
constraint) and whose distribution is invariant under unitary conjugation is the GUE. One can easily
construct matrix ensembles that are invariant under unitary conjugation but in which the matrix
elements are not statistically independent.
Note that the GOE and GUE measures may, up to a normalization constant, be written in the
form10 β 2
e− 4 TrM dM (2.12)
where for the GOE β = 1 and for the GUE β = 2.
It is also worth noting that the GOE matrix elements are real and the GUE matrix entries are
complex. There is a third Gaussian ensemble, which we will not have time to discuss, for which the
matrix elements are quaternions. This is called the Gaussian Symplectic Ensemble and is invariant
under conjugation by all symplectic matrices.

2.5 Random unitary and orthogonal matrices

The circular unitary ensemble, or CUE, is the space of n×n unitary matrices endowed with a probability
measure that is invariant under all unitary transformations (i.e. under left and right multiplication by
all unitary matrices). In this case the matrices represent elements of the compact Lie group U (n) and
the invariant measure is called Haar measure. This measure does not have a simple expression in terms
of the matrix elements, as in the case of Wigner matrices.
Similarly, one can consider the space of n × n orthogonal matrices endowed with a probability
measure that is invariant under all orthogonal transformations (i.e. under left and right multiplication
by all orthogonal matrices). In this case the matrices represent elements of the compact Lie group O(n)
and the invariant measure is the Haar measure on this group. Again, this measure does not have a
simple expression in terms of the matrix elements.
For a careful introduction to Haar measure on compact Lie groups, see, for example, [13] and [14].
The latter reference discusses how to generate random unitary and orthogonal matrices numerically.

10 Note that dM has a different meaning in the two cases.

12
3 The semicircle law for Wigner random matrices
Our goal in this section is to take the first steps in understanding how the eigenvalues of Wigner
random matrices are distributed. As before, let M denote an n × n matrix with real or complex entries
such that for 1 ≤ i ≤ j ≤ n the entries Mij are independent, M = M † in the case when the entries are
complex, and M = M T when the entries are real. Clearly M depends on n, but we will not record this
explicitly, unless it is necessary to do so. We shall assume11 that EMij = 0 and

lim max |nE[|Mij |2 ] − 1| = 0. (3.1)

n→∞ 1≤i,j≤n

We denote the eigenvalues of M by λ1 , λ2 , . . . , λn . It is important to note that in this case we are

taking the mean-square size of the matrix entries of M to be asymptotically 1/n, rather than 1. This
rescaling means that the eigenvalues are effectively rescaled by n−1/2 compared to those of matrices
where the mean-square size of the matrix elements is independent of n.

Exercise Recall the proof that for complex Hermitian matrices (and so as a special case real
symmetric matrices) the eigenvalues are all real.

We will show here the following.

Theorem 1. For bounded continuous functions f (x) and for M satisfying EMij = 0 and (3.1)
n Z
1X
lim f (λi ) = f (x)dσ(x) (3.2)
n→∞ n
i=1

where dσ(x) is given by (

1
√
2π 4 − x2 dx if |x| ≤ 2
dσ(x) = . (3.3)
0 if |x| > 2
and where the convergence holds in expectation and almost surely.
This is known as Wigner’s semicircle law. Various versions of the proof of this result yield different
notions of convergence to the limit: Wigner’s original approach [19] gave convergence in expectation,
but this can be extended to convergence in probability and almost sure convergence, etc. The different
proofs rely on different assumptions about the moments of the matrix elements Mij . We will explain
the proof of convergence in expectation first, and then the extension to almost sure convergence.
Before proceeding we first illustrate the semicircle law by showing the results of some numerical
experiments. These experiments involved generating random matrices from the GOE and the GUE,
plotting histograms of the eigenvalues, and comparing with the semicircle law (3.3). For ease of
visualisation, the eigenvalues have each been divided by 2, so now the support of the rescaled semicircle
is |x| ≤ 1 and it has area π/2. Figure 1 shows the result for a single GOE matrix of dimension 2000,
and Figure 2 shows the result of averaging over 100 GOE matrices of the same dimension. Similarly
Figure 3 shows the result for a single GUE matrix of dimension 2000, and Figure 4 shows the result of
averaging over 100 GUE matrices of the same dimension.
The strategy for proving the semicircle law is as follows.

• The main idea is to compare the moments of the semicircle law (3.3) with the corresponding
spectral moments, i.e. to consider the case when f (x) = xk .
• The computation of the spectral moments maps onto a problem in combinatorics, leading to a
proof that they coincide with the moments of (3.3) in the limit when n → ∞.
11 These assumptions can be relaxed in various ways but are sufficient for our purposes here. Following the steps in the

proof may be easier if one has in mind the weaker condition that E[|Mij |2 ] = 1/n; the difference will not be important to
us.

13
 0.6

0.5

0.4

���� ��
0.3

0.2

0.1

0.0
-1.0 -0.5 0.0 0.5 1.0

Figure 1: The spectral density of a single GOE matrix of dimension 2000 compared to the Wigner
semicircle law.

0.6

0.5

0.4

���� ��
0.3

0.2

0.1

0.0
-1.0 -0.5 0.0 0.5 1.0

Figure 2: The spectral density obtained by averaging over 100 GOE matrices of dimension 2000
compared to the Wigner semicircle law.

14
 0.6

0.5

0.4

���� ��
0.3

0.2

0.1

0.0
-1.0 -0.5 0.0 0.5 1.0

Figure 3: The spectral density of a single GUE matrix of dimension 2000 compared to the Wigner
semicircle law.

0.6

0.5

0.4

���� ��
0.3

0.2

0.1

0.0
-1.0 -0.5 0.0 0.5 1.0

Figure 4: The spectral density obtained by averaging over 100 GUE matrices of dimension 2000
compared to the Wigner semicircle law.

15
 • We first show this for convergence in expectation.
• A variance calculation then implies almost sure convergence.

• This proves (3.3) when f (x) is a polynomial function of x.

• The proof of (3.3) for bounded continuous functions follows from the fact that such functions can
be approximated by polynomials. We shall discuss this at the end.

3.1 Moments of the semicircle law and Catalan numbers

The moments of the semicircle law (3.3) are
Z 2 p
1
αk = xk 4 − x2 dx (3.4)
2π −2

Clearly, by symmetry, αk = 0 when k is odd. Setting k = 2m we then have

Z 2 p
1
α2m = x2m 4 − x2 dx
2π −2
Z
22m+1 π
= (cos θ)2m sin2 θdθ
π 0
Z π Z π 
22m+1
= (cos θ)2m dθ − (cos θ)2m+2 dθ
π
Z π 0 0
Z 
2 −iθ 2m 1 π iθ
= iθ
(e + e ) dθ − (e + e−iθ )2m+2 dθ
π 0 4 0
   
2m 1 2m + 2
=2 −
m 4 m+1
 
1 2m
= , (3.5)
m+1 m

where, in passing from the fourth to the fifth line, we used the binomial expansion and the fact that
only one term in this expansion gives a non-zero contribution when integrated.
The numbers  
1 2m
Cm = (3.6)
m+1 m
are known as the Catalan numbers. They satisfy
   
2m 2m
Cm = − , (3.7)
m m+1

the recurrence relation

m
X
Cm+1 = Ci Cm−i , (3.8)
i=0

with C0 = 1, and they are all integers12 . The sequence of Catalan numbers begins 1, 1, 2, 5, 14, 42,
132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440, 9694845, . . . .
Note that it follows from (3.6) that Cm ≤ 4m .
The Catalan numbers play an important role in combinatorics because they count various interesting
objects, ranging from lattice paths to geometric constructions in polygons. See, for example, [10, 17].
We now describe briefly one such connection.
12 This is straightforward to see from (3.7) or from (3.8).

16
 For k ∈ N, the set Dk of Dyck paths is the set of paths of length 2k starting at 0 and ending at 0
with increments ±1 and which remain non-negative13 . Then |Dk | = Ck .
To prove this, let A denote the set of walks of length 2k with increments ±1 which start at 0 and
end at -2, and similarly let B denote the set of all walks of length 2k with increments ±1 which start
at 0 and end at 0. Then |Dk | = |B| − |A|, which may be seen in the following way. Let C denote the
set of walks of length 2k with increments ±1 which start at 0 and end at 0, but which hit -1 at some
intermediate point. Obviously |Dk | = |B| − |C|. However, |C| = |A| because the respective paths are in
bijection: for a walk in C, let j be the last visit to −1; reflecting the portion of the path after j about
−1 gives a unique path

terminating at −2 with the same set of visits to −1.

Clearly |B| = 2kk , because B consists of k steps +1 and k steps −1 in some order, and |A| = k+1 ,
2k

because A consists of k − 1 steps +1 and k + 1 steps −1. Hence

   
2k 2k
|Dk | = − = Ck . (3.9)
k k+1

Alternatively, consider a Dyck path of length 2(k + 1). Let 2(m + 1) denote the first ‘time’ this
path hits 0. We can break the path into two pieces: the section up to time 2(m + 1) and the section
after this time. The path up to time m is itself a Dyck path, as is the path afterwards. If we denote
|Dn | = dn , then clearly
Xk
dk+1 = dm dk−m , (3.10)
m=0

because there are dm dk−m paths which first hit 0 at time 2(m + 1), and d0 = 1. Hence dk satisfies the
recurrence relation (3.8).

3.2 Spectral moments

We begin by focusing on polynomial functions f (x) in (3.2). This means that we need only consider
f (x) = xk . The sum on the left hand side of (3.2), which we refer to as the kth spectral moment, is
then simply
n
1X k 1
λi = TrM k . (3.11)
n i=1 n
Our goal now is to prove the following theorem.
Theorem 2. (Wigner 1958) Assume that for all k ∈ N
√
Bk := sup sup E[| nMij |k ] < ∞. (3.12)
n∈N (i,j)∈{1,...,n}2

Then (
1 0 if k is odd
lim TrM k = (3.13)
n→∞ n Ck otherwise
2

where the convergence holds in expectation and almost surely.

The main idea underlying the proof of this theorem is to expand TrM k in terms of products
of elements of M , to use the independence of these to select those products that make a non-zero
contribution to the expectation value of the trace, and then from these to identify the products of
elements that make the largest contribution when n → ∞.
We begin by proving Theorem 2 for convergence in expectation. When this is done we explain how
to extend the proof to almost sure convergence.
13 Note that D is in obvious bijection with the "latex bracketing problem", that is with the number of ways of writing
k
2k brackets { and } starting with { and in such a way that the brackets correctly match, if we identify { with a step of
+1 and } with a step of −1.

17
 √ √
First, we set µij = nMij , and denote µ = nM . Then
  n
X
1 k 1
E TrM = E[µi1 i2 µi2 i3 . . . µik i1 ]. (3.14)
n i
n1+k/2
1 ,i2 ,...,ik =1

Note that for all sets of indices, E[µi1 i2 µi2 i3 . . . µik i1 ] is uniformly bounded by Bk . Each set of indices
corresponds to a sequence i1 i2 . . . ik i1 , which can be thought of as a closed connected path on the set
of vertices {i1 , i2 , . . . , ik } with edges corresponding to pairs of consecutive indices ij ij+1 and with the
convention that ik+1 = i1 . We denote the graph with this set of vertices and edges G(i1 ,i2 ,...,ik ) . Clearly
the graph is connected. One can think of the paths as closed walks on this graph. Crucially, because
the matrix entries are independent, centred random variables, E[µi1 i2 µi2 i3 . . . µik i1 ] = 0 for all paths
except those for which every edge is traversed at least twice, possibly in reverse. There can therefore
be at most k/2 unique edges and hence at most k/2 + 1 distinct vertices in paths giving a non-zero
contribution to the expectation value.
For a sequence i1 i2 . . . ik i1 we define its weight w to be the number of distinct indices. The terms
making a non-zero contribution to (3.14) have w ≤ k/2 + 1. Two sequences i1 i2 . . . ik i1 and i01 i02 . . . i0k i01
are termed equivalent if there exists a bijection on the set {1, 2, . . . , n} mapping each ij to i0j . Sequences
that are equivalent have essentially the same weight and their contributions to (3.14) are asymptotically
the same. Note that the total number of distinct equivalent classes depends on k but is independent of
n, because each class has a representative where all i1 i2 . . . ik are in {1, 2, . . . , k}.
The next step is to show that those paths with weight w < k/2 + 1 make a contribution to (3.14)
that vanishes in the limit as n → ∞. As noted above, E[µi1 i2 µi2 i3 . . . µik i1 ] is uniformly bounded by
Bk . For a given index sequence i1 i2 . . . ik i1 of weight w, the number of sequences equivalent to it
is n(n − 1)(n − 2) . . . (n − w + 1) ≤ nw . The contribution of the terms in this equivalence class to
(3.14) is therefore ≤ Bk nw−k/2−1 and if w < k/2 + 1 this tends to zero as n → ∞. Therefore, because
the number of distinct equivalent classes is independent of n, the total contribution from paths with
w < k/2 + 1 tends to zero in the limit, and when k is even it does so at least as fast as 1/n.
It follows that the only sequences making a non-zero contribution to (3.14) in the limit n → ∞ are
those for which w = k/2 + 1. When k is odd this equation has no solution, because w is an integer.
Hence when k is odd the left hand side of (3.14) tends to zero in the limit n → ∞.
Consider now the case when k is even and w = k/2 + 1. In this case the graphs we need to consider
are connected, have k/2 + 1 vertices and k/2 distinct edges. They are therefore trees.
To see this note that for any finite connected graph G = (V, E) consisting of a set of vertices V and
a set of edges E,
|V | ≤ |E| + 1, (3.15)
with equality if and only if G is a tree. If G is a tree, it follows easily by induction that |E| = |V | − 1.
If G is not a tree, let T = (V, E 0 ) denote a spanning tree of G; i.e. a tree which has the same vertices as
G and for which E 0 ⊆ E. Then |E| ≥ |E 0 | = |V | − 1. If |E| = |V | − 1 then we would have |E| = |E 0 |,
and so G = T .
When k is even and w = k/2 + 1, the sequence i1 i2 . . . ik i1 therefore corresponds to a closed path on
a tree that traverses each edge exactly twice, once in each direction. This means that for any 1 ≤ j ≤ k
ij+1 6= ij . Hence in (3.14) only off-diagonal matrix entries contribute, and these are precisely paired, so
each appears squared. Therefore, using the fact that14 E|µij |2 ∼ 1 for i = 6 j, and that the off-diagonal
matrix entries are independent, we have

E[µi1 i2 µi2 i3 . . . µik i1 ] ∼ 1. (3.16)

Recall that for a given index sequence i1 i2 . . . ik i1 of weight w, the number of sequences equivalent to it
is n(n − 1)(n − 2) . . . (n − w + 1), and that we now have w = k/2 + 1, so that this number grows like
n1+k/2 , as n → ∞ which exactly cancels the factor multiplying the expectation in (3.14).
14 When f (n)/g(n) → 1 as n → ∞ we write f (n) ∼ g(n).

18
 The problem has therefore reduced to that of counting the number of sequences i1 i2 . . . ik i1 corre-
sponding to paths of length k on the tree with k/2 + 1 vertices and k/2 distinct edges, each traversed
twice. Such paths are called non-crossing. We can perform this count in the following way. As a
given path is traversed, record at each edge on it whether that edge has been traversed before, or
not. If it is being traversed for the first time, term it open, if it is being traversed for the second
(and last) time, term in closed. For each non-crossing path of length k we associate a sequence, called
the path sequence, whose jth entry is the number of open edges minus the number of closed edges
in the path i1 i2 . . . ij+1 . Now we reach the main point: this sequence starts with a 1 and ends with
a 0, and successive terms differ by ±1. For example, for the path corresponding to the set of edges
({1, 2}, {2, 3}, {3, 2}, {2, 4}, {4, 5}, {5, 4}, {4, 6}, {6, 4}, {4, 2}, {2, 1}), the edges are respectively open,
open, closed, open, open, closed, open, closed, closed, closed, and the associated sequence is 1, 2, 1, 2,
3, 2, 3, 2, 1, 0. If we attach a label 0 at the start, then these sequences are clearly in bijection with the
set of Dyck paths of length k.
Collecting together what we have shown so far, when k is even
1
lim ETrM k = #{path sequences of length k}
n→∞ n
= #{Dyck paths of length k}
= Ck/2 , (3.17)

and when k is odd, the limit is zero. This proves Theorem 2 for convergence in expectation.
We next indicate how to prove almost sure convergence to the same limit. The strategy here is,
having computed the expectation value of the moments, now to compute the variance. We will see that
the variance vanishes as n → ∞ sufficiently fast that a sequence of matrices with increasing size will
converge to the mean almost surely. This is proved by appealing to the Borel-Cantelli lemma. It is an
example of a phenomenon known as concentration of measure that plays an important role in RMT.
The variance is    2  2
1 k 1 k 1 k
var TrM = E TrM − E TrM (3.18)
n n n
and so
  n
X n
X
1 k 1
var TrM = 2+k E[µi1 i2 µi2 i3 . . . µik i1 µj1 j2 µj2 j3 . . . µjk j1 ]
n n i1 ,i2 ,...,ik =1 j1 ,j2 ,...,jk =1

− E[µi1 i2 µi2 i3 . . . µik i1 ]E[µj1 j2 µj2 j3 . . . µjk j1 ]. (3.19)

As in the calculation of the expectation value of the moments above, one can analyse this sum in
terms of paths on a graph. In the case of the first term on the right hand side of (3.19), the vertices of
the graph in question are labelled by {i1 , i2 , . . . , ik , j1 , j2 , . . . , jk }. For this graph to contribute to the
sum, it must be connected15 , otherwise the first term factorises and is exactly cancelled by the second
term. As before, each edge must appear at least twice and must be traversed an even number of times.
From the arguments above, the terms making a non-zero contribution to (3.19) have w ≤ k + 1.
The fact that the summand is uniformly bounded implies that the variance tends to zero at least as
fast as 1/n. However, we shall need a better bound, and we can find one as follows. For the reasons
explained in the previous calculation, the largest contribution is expected to come from paths with
w = k + 1, which is when the graph is a tree and each edge appears exactly twice, traversed once in
each direction. However, in this case no such terms exist. To see this note that for a pair {i1 , i2 , . . . , ik }
and {j1 , j2 , . . . , jk } to make a non-zero contribution, they must have at least one edge in common;
otherwise the first term in (3.19) factorizes and is cancelled by the second term. However, if each edge
is to appear exactly twice in the i-path or the j-path, this is impossible unless the graph has a closed
cycle; and this is a contradiction, because the graph is a tree.
15 and so i1 = j1

19
 It follows that there are no contributing paths with w = k + 1, and so in fact the variance tends to
zero as n → ∞ at least as fast as 1/n2 .
We now apply Chebychev’s inequality, which asserts that for a random variable X with zero mean
and variance σ 2
σ2
P(|X| ≥ δ) ≤ 2 (3.20)
δ
for any δ > 016 . In our case, the fact that the variance tends to zero at least as fast as 1/n2 then means
that for any δ > 0  
1 k 1 k αk
P TrM − E TrM > δ ≤ 2 2 (3.21)
n n δ n
where αk is a uniform bound on the summand in (3.19).
The Borel-Cantelli lemma asserts that if E1 , E2 , . . . is a sequence of events in some probability
space, and if the sum of their probabilities is finite, i.e. if
∞
X
P(En ) < ∞, (3.22)
n=1

then
P(lim sup En ) = 0 (3.23)
n→∞

i.e. the probability that infinitely many of them occur is 0.

We can apply the Borel-Cantelli lemma to (3.21) precisely because the sum of the right-hand side
over n converges17 , giving that the convergence in Theorem 2 holds almost surely; that is,
1 1
lim TrM k − E TrM k = 0 (3.24)
n→∞ n n
almost surely.
We have therefore established Theorem 2 for polynomials f (x). It remains to show that the
result extends to all bounded continuous functions f (x). One can do this using the Weierstrass
approximation theorem, which asserts that in a bounded domain f (x) can be approximated by a
polynomial. Specifically, for any B and δ > 0, one can find a polynomial pδ (x) such that
δ
sup |f (x) − pδ (x)| < . (3.25)
|x|≤B 6

We now introduce some additional notation. Let

(
1 if z ≤ 0
Θ(z) = (3.26)
0 if z > 0.

We define the empirical spectral distribution function to be

n
1X
µM (x) := Θ(λj − x). (3.27)
n j=1

So this is the fraction of eigenvalues of M less than or equal to x. Note that

n Z
1X ∞
φ(λj ) = φ(x)dµM (x). (3.28)
n j=1 −∞

16 Chebychev’s inequality follows from Markov’s inequality: P(|X|2 > (δσ)2 ) ≤ E[|X|2 ]/(δσ)2 .
17 This explains why the straightforward bound, 1/n, was not good enough for our purposes and why we needed the
stronger bound (3.21) .

20
We also define
n
1X
µ̄(x) := E Θ(λj − x). (3.29)
n j=1

We start with the following result which we will need.

Lemma 3. Let p ∈ N and  > 0. For B > 4
Z !
p
lim sup P |x| dµM (x) >  =0 (3.30)
n→∞ |x|>B

Proof. By Markov’s inequality

Z ! Z !
p 1 p
P |x| dµM (x) >  ≤ E |x| dµM (x) (3.31)
|x|>B  |x|>B

and Z Z
1
|x|p dµM (x) ≤ x2p dµM (x), (3.32)
|x|>B Bp
which follows by multiplying the integrand by xp /B p .
Hence Z ! Z 
p 1 2p 1 ETrM 2p
P |x| dµM (x) >  ≤ E x dµM (x) = . (3.33)
|x|>B B p B p n
We have shown above that the right hand side converges to Cp /B p , and so, for n sufficiently large,
using the fact that Cp ≤ 4p ,
Z !  p
p 1 4
P |x| dµM (x) >  ≤ . (3.34)
|x|>B  B

In the range of integration, |x|p is strictly increasing, and so the left hand side increases as p increases.
However, the right hand side decreases as p increases. This is only possible if the limsup of the sequence
on the left is actually 0, as was to be shown.

We return now to proving that Theorem 2 extends to all bounded continuous functions f (x), making
use of this Lemma. Let f (x) be a (fixed) bounded continuous function, fix B > 4 and δ > 0. By the
triangle inequality,
Z Z Z Z
f (x)dµM (x) − f (x)dσ(x) ≤ f (x)dµM (x) − pδ (x)dµM (x)
Z Z
+ pδ (x)dµM (x) − pδ (x)dσ(x)
Z Z
+ pδ (x)dσ(x) − f (x)dσ(x) , (3.35)

where σ(x) is the semicircle law. Therefore,

Z Z  Z Z 
P f (x)dµM (x) − f (x)dσ(x) > δ ≤ P f (x)dµM (x) − pδ (x)dµM (x) > δ/3
Z Z 
+P pδ (x)dµM (x) − pδ (x)dσ(x) > δ/3
Z Z 
+P pδ (x)dσ(x) − f (x)dσ(x) > δ/3 . (3.36)

21
By design of (3.25), the last term is identically zero. In the first term we can use the triangle inequality
Z Z Z Z
f (x)dµM (x) − pδ (x)dµM (x) ≤ |f (x) − pδ (x)|dµM (x) + |f (x) − pδ (x)|dµM (x).
|x|≤B |x|>B
(3.37)
Hence
Z Z  Z !
P f (x)dµM (x) − pδ (x)dµM (x) > δ/3 ≤ P |f (x) − pδ (x)|dµM (x) > δ/6
|x|≤B
Z !
+P |f (x) − pδ (x)|dµM (x) > δ/6 . (3.38)
|x|>B

Again, it follows from (3.25) that the first term is identically zero.
We are therefore left with
Z Z  Z Z 
P f (x)dµM (x) − f (x)dσ(x) > δ ≤ P pδ (x)dµM (x) − pδ (x)dσ(x) > δ/3
Z !
+P |f (x) − pδ (x)|dµM (x) > δ/6 . (3.39)
|x|>B

We have already shown that the first term tends to zero as n → ∞. In the second term, because f (x)
is bounded and pδ (x) is a polynomial, we have |f (x) − pδ (x)| ≤ c|x|k for some positive constant c and
some integer k. Therefore
Z ! Z !
P |f (x) − pδ (x)|dµM (x) > δ/6 ≤P c|x|k dµM (x) > δ/6 . (3.40)
|x|>B |x|>B

It follows from Lemma 3 that as n → ∞, the lim sup of the righthand side is 0, concluding the proof of
Theorem 2.
We note finally that establishing Theorem 2 implies that the measure dµM (x) almost surely converges
weakly to the measure dσ(x) when n → ∞.

3.3 The circular law

The semicircle law holds for either real-symmetric or complex-Hermitian matrices, and in both cases the
eigenvalues are real. It is natural to seek to extend it to non-symmetric matrices, when the eigenvalues
are complex numbers. This turns out to be harder, so we simply state the result.
Let X be a centred complex random variable with mean 0 and variance 1, and let M be the random
n × n matrix with entries being i.i.d. copies of X. Let λ1 , λ2 , . . . , λn be the eigenvalues of √1n M . Define
the empirical spectral distribution function µM by
1
µM (s, t) = #{k ≤ n|Re(λk ) ≤ s; Im(λk ) ≤ t}. (3.41)
n
We shall compare this to the uniform distribution on the unit disk in the complex plane
1
ν(s, t) = meas{z ∈ C : |z| ≤ 1; Re(z) ≤ s; Im(z) ≤ t}. (3.42)
π
The weak circular law is said to hold for X if for any fixed s, t µM converges to ν in probability.
The strong circular law is said to hold for X if, with probability 1, µM converges uniformly to ν over
the unit disk as n → ∞. Both laws hold under the general assumptions on X stated.

22
3.4 Applications
The semicircle law has applications in settings where real-symmetric or complex-Hermitian matrices
arise and can be modelled by random matrices. As noted in the introduction, many problems in
population dynamics, modelling random neural networks, etc, reduce to the following mathematical
question.
Let x(t) be an n-dimensional vector satisfying

dx
= −λIx (3.43)
dt
where λ is a positive constant and I is the n × n identity matrix. This system has a stable fixed point
at x = 0. Clearly in (3.43) the different components of x are uncoupled (because I is diagonal). The
question is: if one now introduces random couplings between these components, does this make the
fixed point more or less stable? Specifically, what is the stability of the system
dx
= −λIx + Ax (3.44)
dt
where A is a random matrix. What, for example, should one expect if |λ| is large compared to the
typical size of the elements of A?
In population dynamics, this is called the May model. Then the components of x represent
the differences of the populations of various species from some equilibrium values and A represents
interactions between the different species.
Now, if A is a real-symmetric random matrix in which the mean-square size of the matrix entries
(i.e.√the interactions)
√ is 1, it is a consequence of the semicircle law that the eigenvalues of A lie between
−2 n and 2 n with increasing probability as n grows. We expect to find eigenvalues over the whole
of this range, but very few outside it. And any outside it will not be far from it. So when n > λ2 /4
we expect to see a transition in the dynamics from exponential stability to exponential instability. In
large random neural networks, this question of stability is related to a transition in the dynamics from
ordered to ‘chaotic’.
If A is not symmetric, so the matrix entries are independent random variables, the circular law
implies a similar transition occurs for large enough systems. In this case the stability also depends on
the eigenvectors, which no longer need be (and typically are not) orthogonal, so the details are more
complicated, but the qualitative picture is the same.
Now, one can question whether it is reasonable to model the interactions as i.i.d. random variables
with the same variance. In population dynamics, some species pairs interact far more strongly than
others. For example, lions interact strongly with zebras, wildebeests, and impalas (and the interaction
is clearly not symmetric!), but barely at all with butterflies or fish. There are food chains and these
give rise to complex networks. Building this topological structure into population models is a significant
area of research involving random matrices with structured correlations between the matrix entries
and where the matrices may be rather sparse. It is also an important area of research in neuroscience,
where again the connectivity matrix for neural pathways is both highly structured and also stochastic.
Finally, we remark that this discussion relates to linear stability. Extending it to nonlinear stability
analysis is currently an active area of research.

23
4 Gaussian random matrices
In the previous section we calculated the spectral moments for general Wigner random matrices. We
now look at the special case when the matrix elements have a Gaussian distribution. In this case we
can go more deeply into the formula embodied in Theorem 2. Recall that Wigner’s theorem relates the
moments to the Catalan numbers in the limit as n → ∞, and that this is proved by linking them to
counts of certain non-crossing paths on graphs. In this limit they therefore have an interpretation in
terms of counting topologically defined objects. We shall see that for Gaussian random matrices this
may be viewed as part of a bigger picture.
To start with, we need to explain a formula that will prove useful in analysing the expectation value
on the right-hand side of (3.14).

4.1 Wick’s theorem

Let X be a real standard normal random variable, i.e. one with mean 0 and variance 1, so that
Z b
1 2
P(a ≤ X ≤ b) = √ e−y /2
dy. (4.1)
2π a

Hence Z (
∞
n 1 n −y 2 /2 0 if n is odd
E[X ] = √ y e dy = (4.2)
2π −∞ (n − 1)!! if n is even
where, when m is odd, m!! = m(m − 2)(m − 4) . . . 1. Therefore E[X n ] is equal to the number of ways
of splitting n objects into disjoint pairs. For example E[X 4 ] = 3, corresponding to the 3 ways to split
the numbers 1, 2, 3, 4 into disjoint pairs (i.e in this case [(1, 2)(3, 4)], [(1, 3)(2, 4)], [(1, 4)(2, 3)]).
We label the different ways of splitting n objects into disjoint pairs by π and the set of all pairings
by P2 (n), so #P2 (n) = (n − 1)!! when n is even and zero otherwise.
The following is then a straightforward application of (4.2) due to Isserlis in 1918 and Wick in 1950
(and many others).
Theorem 4. Let Y1 , . . . , Yp be independent standard Gaussian random variables and consider x1 , . . . , xn ∈
{Y1 , . . . , Yp }. Then X Y
E[x1 x2 . . . xn ] = E[xi xj ]. (4.3)
π∈P2 (n) (i,j)∈π

Importantly, the covariances in the Wick formula can be computed straightforwardly from (4.2):
(
1 if xi = xj
E[xi xj ] = (4.4)
0 if xi 6= xj .

Note that because it is linear in the variables x1 , . . . , xn , this√formula applies in exactly the same way
to standard complex Gaussian random variables Z = (X +iY )/ 2, where X and Y are independent real
standard Gaussian random variables. In this case E[Z] = E[Z 2 ] = 0 and E[|Z|2 ] = 1, and if Z1 , . . . , Zp
are independent complex standard Gaussian random variables then for z1 , . . . , zn ∈ {Z1 , . . . , Zp },
E[z1 z2 . . . zk ] is equal to the number of parings which connect zi with its complex conjugate.
Historically, Isserlis introduced this formula in the context of statistical analysis. It was later
re-introduced by Wick to count Feynman diagrams in Quantum Field Theory and it tends now to be
associated with his name.

24
4.2 The genus expansion for Gaussian random matrices
Recall that the spectral moments are given by
  X n
1 k 1
E TrM = E[µi1 i2 µi2 i3 . . . µik i1 ]. (4.5)
n i ,i ,...,i
n1+k/2
1 2 k =1

In the previous section we analysed this for general Wigner matrices and found that when n → ∞ the
leading order asymptotics may be calculated by counting non-crossing paths. If we now specialize to
complex Hermitian Gausssian random matrices, for example, we can apply Wick’s theorem to analyse
this expression much more explicitly.
We have from Wick’s theorem that
X Y
E[µi1 i2 µi2 i3 . . . µik i1 ] = E[µia ,ia+1 µib ,ib+1 ] (4.6)
π∈P2 (k) (a,b)∈π

where ik+1 = i1 . Hence,

  X n
X Y
1 k 1
E TrM = δia ,ib+1 δia+1 ,ib , (4.7)
n n1+k/2
π∈P2 (k) i1 ,i2 ,...,ik =1 (a,b)∈π

where δi,j denotes the Kronecker δ-symbol.18

We can think of π as a particular kind of permutation, namely one that factorizes into permutations
between the elements being paired, so permutations made up of cycles of length 2. Hence (a, b) ∈ π
means that π(a) = b, or equivalently π(b) = a. Therefore

Y k
Y
δia ,ib+1 δia+1 ,ib = δia ,iπ(a)+1 . (4.8)
(a,b)∈π a=1

If we introduce the shift notation γ(a) = a + 1 mod k on the set {1, 2, . . . , k}, so that γ ∈ Sk is a
permutation with cycle (1, 2, . . . , k), then we have that
  X n
X k
Y
1 1
E TrM k = δia ,iγπ(a) . (4.9)
n n1+k/2 a=1
π∈P2 (k) i1 ,i2 ,...,ik =1

One can think of the indices {i1 , i2 , . . . , ik } as a function i : {1, . . . , k} → {1, . . . , n}. Then
k
(
Y 1 if i is constant on the cycles of γπ
δia ,iγπ(a) = (4.10)
a=1
0 otherwise.

We therefore have that

  X
1 k 1
E TrM = 1+k/2 #{i : {1, . . . , k} → {1, . . . , n} : i is constant on the cycles of γπ}.
n n
π∈P2 (k)
(4.11)
The count in the summand is straightforward: we simply need to choose one value for each of the cycles
in γπ (with repeats allowed). For any permutation σ ∈ Sk , let #(σ) denote the number of cycles in σ.
Then we have  
1 X
E TrM k = n#(γπ)−k/2−1 . (4.12)
n
π∈P2 (k)

18 Recall that for standard complex normal random variables, we have to pair a random variable in Wick’s theorem

with its complex conjugate, and that the matrices we are dealing with are Hermitian.

25
 Equation (4.12) is an exact formula for the spectral moments of M . To understand its structure,
note first that if k is odd then P2 (k) is empty and so the corresponding spectral moment is identically
zero. We therefore now set k = 2m, so that
  X
1 2m
E TrM = n#(γπ)−m−1 . (4.13)
n
π∈P2 (2m)

This sum is known as the genus expansion, because there is an interpretation of #(γπ) that combines
geometric and topological ideas, which we now explain.
We start by drawing a 2m-gon and labelling its vertices in cyclic order v1 , v2 , . . . , v2m . We can then
label its edges by the vertices they connect; so for i = 1, . . . , 2m − 1, ei = vi vi+1 , and e2m = v2m v1 .
A pairing π ∈ P2 (2m) corresponds to an association between pairs of edges, which can then be glued
together to form a compact surface. If this is done so that when vi vi+1 is glued to vj vj+1 , vi is glued
to vj+1 , and vi+1 to vj , then the surface that results is orientable. So, for example, if π(1) = 3 one
identifies e1 and e3 by gluing v1 to v4 and v2 to v3 . We term this the "tail-to-head" convention.
Consider now the surface S̃π obtained by this gluing procedure. The number of distinct vertices in
the image Gπ of the 2m-gon in S̃π is precisely #(γπ). To see this, note that ei is glued to eπ(i) , and so
vi is glued to vγπ(i) for each i ∈ {1, . . . , 2m}. The edge eγπ(i) is glued to eπγπ(i) and so vγπ(i) , which is
now the tail of the edge in question, gets glued to vγπγπ(i) etc. Continuing on this way, we see that vi
ends up being identified with precisely those vj for which j = (γπ)l (i) for some l ∈ N. Therefore the
cycles of γπ count the number of distinct vertices after gluing.
The Euler characteristic of S̃π , χ(S̃π ), is an even integer which may be defined as follows: if G is
any embedded polygonal complex in S̃π , then
χ(S̃π ) = (number of vertices in G) − (number of edges in G) + (number of faces of G). (4.14)
Now, any orientable compact surface is homeomorphic to a g-holed torus for some g ≥ 0 (the g = 0
case is the sphere). The topological invariant g is known as the genus of the surface. It is a theorem of
Cauchy that the Euler characteristic is 2 minus twice the genus, and so
χ(S̃π ) = 2 − 2g(S̃π ). (4.15)
The embedded complex in S̃π that we constructed above has one face, m edges, and #(γπ) vertices,
and so
2 − 2g(S̃π ) = χ(S̃π ) = #(γπ) − m + 1. (4.16)
We therefore have, finally, that
  X X
1
E TrM 2m = n−2g(S̃π ) = τg (m)n−2g (4.17)
n
π∈P2 (2m) g≥0

where
τg (m) = #{genus-g surfaces obtained by gluing together pairs of edges of a 2m-gon}. (4.18)
2m
This remarkable formula therefore implies that averaging TrM over the GUE leads to a method for
counting the number of genus-g surfaces obtained by gluing together pairs of edges of a 2m-gon!
We see immediately that since g ≥ 0,
 
1
lim E TrM 2m = τ0 (m) = #{spheres obtained by gluing together pairs of edges of a 2m-gon}
n→∞ n
(4.19)
and so we see that by Wigner’s theorem
#{spheres obtained by gluing together pairs of edges of a 2m-gon} = Cm (4.20)
where Cm is the m-th Catalan number.
This is the simplest of a wide and important variety of examples where random matrix theory
connects with enumerative topology.

26
5 The Marchenko-Pastur law for Wishart random matrices
The Marchenko-Pastur19 law is the analogue for Wishart random matrices of the Wigner semicircle law
for Hermitian matrices.

5.1 The Marchenko-Pastur Theorem

Let X be a p × n matrix with entries Xij that are i.i.d. real random variables with EXij = 0 and
2
EXij = 1. Denote by σn the p × p matrix

1
σn = XX T ∈ Rp×p (5.1)
n
and denote by 0 ≤ λ1 ≤ λ2 ≤ · · · ≤ λp the eigenvalues of σn .
Defining the random spectral distribution by
p
1X
µσ (x) := Θ(λj − x), (5.2)
p j=1

The Marchenko-Pastur law asserts that

Theorem 5. For σn and µσ as defined above, when p → ∞ and n → ∞ such that p/n → γ ∈ (0, 1],
µσ → µ in expectation and almost surely, where µ is the deterministic measure satisfying
( p
1
dµ (a+ − x)(x − a− ) if a− ≤ x ≤ a+
= 2πγx (5.3)
dx 0 otherwise,
√ √
where a− = (1 − γ)2 and a+ = (1 + γ)2 . When γ > 1, one needs to add (1 − γ −1 )δ(x) to the
right-hand side of (5.3), where δ(x) denotes a Dirac delta-function/mass at the origin.

The additional factor (1 − γ −1 )δ(x) which needs to be added to the right-hand side of (5.3) when
γ > 1 is explained by the fact that, since the rank of σn is the smaller of p and n, there are then
approximately p − n zero eigenvalues which will contribute a mass of (1 − γ −1 ) at 0 in the limiting
measure.
Note that when γ = 1, a− = 0 and a+ = 4, and that then (5.3) coincides with the semicircle law
under the mapping x → x2 .
We illustrate the Marchenko-Pastur law by showing the results of numerical experiments. These
involved generating random matrices σn , plotting histograms of the eigenvalues, and comparing with
(5.3). Figure 5 shows the result when p = 103 and n = 104 , compared to (5.3) with γ = 0.1. Similarly
Figure 6 shows the result when p = 1000 and n = 2000, compared to (5.3) with γ = 0.5.
One can prove the Marchenko-Pastur formula in a similar way to the proof of the semicircle law
already described. The main steps are:

• proving that the spectral moments converge in expectation to the moments of the limiting
deterministic Marchenko-Pastur distribution when n → ∞;
• showing that the variance vanishes in the limit when n → ∞, and then using the Borel-Cantelli
lemma to conclude that the moments converge almost surely;
• demonstrating that the distribution is determined by its moments.
19 Sometimes written Marčenko-Pastur.

27
1.2

1.0

0.8

0.6

0.4

0.2

0.0
0.6 0.8 1.0 1.2 1.4 1.6

Figure 5: Eigenvalue density when p = 103 and n = 104 , compared to (5.3) with γ = 0.1.

28
0.8

0.6

0.4

0.2

0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0

Figure 6: Eigenvalue density when p = 1000 and n = 2000, compared to (5.3) with γ = 0.5.

We will describe the proof of the first step: the convergence of the spectral moments to those of
the limiting deterministic Marchenko-Pastur distribution in expectation. The second and third steps
follow very similar lines to the corresponding ones in the proof of the semicircle law and so we shall not
repeat them here. We focus on the case when γ ≤ 1.
√
The moments of the Marchenko-Pastur distribution are, recalling that a− = (1 − γ)2 and
√ 2
a+ = (1 + γ) , given by
Z a+
1 p
α̃k = xk (a+ − x)(x − a− )dx
a− 2πγx
X γ r k k − 1
k−1
= . (5.4)
r=0
r+1 r r

To see this note that a− + a+ = 2(1 + γ) and that a− a+ = (1 − γ)2 , and hence
Z 2 p
1 √
α̃k = ( γy + 1 + γ)k−1 4 − y 2 dy. (5.5)
2π −2
√
Expanding ( γy + 1 + γ)k−1 binomially in powers of y, and using the formula for the moments of the
semicircle law and Vandermonde’s identity
  Xk   
m+n m n
= (5.6)
k r=0
r k−r

then leads immediately to the result.

29
 The goal is to show that
p  k
1X k 1 1
E λ = E Tr XX T (5.7)
p i=1 i p n

converges to (5.4) in the appropriate limit.

As in the case of Wigner’s semicircle theorem, we proceed by expanding the trace in terms of the
matrix entries and using the independence of these to identify the products that make the largest
contributions. So,
 k
1 1 1 X
E Tr XX T = E [Xi1 j1 Xi2 j1 Xi2 j2 Xi3 j2 . . . Xik jk Xi1 jk ] (5.8)
p n pnk
IJ

where I ∈ [p]k and J ∈ [n]k .

It is helpful to picture this in terms of cycles of length 2k on a directed bipartite graph with one set
of vertices labelled by i1 , i2 , . . . , ik and the second set by j1 , j2 , . . . , jk , with directed edges from i1 to
j1 , j1 to i2 , i2 to j2 , j2 to i3 etc, and with finally an edge from jk to i1 .
Since EXij = 0 and the Xij are independent, the summand in (5.8) vanishes unless each edge
appears at least twice in a cycle. Hence there can be at most k edges and therefore k + 1 vertices in the
skeleton of the graph; i.e. in the subgraph consisting of the vertices and edges visited, with multiplicities
removed. The number of vertices corresponds to the weight of the contributing cycle.
We consider first the case when the weight w < k + 1. Let the cycle visit lI i-vertices and lJ
j-vertices (so w = lI + lJ ). The total number of ways of choosing these sites is bounded from above by
a constant (i.e. a factor independent of n) times plI nlJ . Hence the contribution of all such terms in
(5.8) is bounded from above by a constant (i.e. a factor independent of n) times plI nlJ /pnk , which
tends to zero as n → ∞ if p ∼ γn.
As in the proof of the semicircle law, we therefore need to focus on the case when there are exactly k
edges and k + 1 vertices. The graph then is a double tree, namely a tree where vertices can be labelled
either by i or j and these are visited alternately on the cycle, which passes through each edge exactly
twice. So, if we start at an i-vertex, other i-vertices can be reached if and only if the number of steps
taken is even, and j-vertices reached if and only if the number of steps taken is odd.
We call the vertices numbered in order of appearance, up to an overall shift in the numbers assigned,
the shape of the double tree. The question then is: how many double trees are there with a given
shape? The answer is that we need to choose r + 1 i-vertices from [p] and k − r j-vertices from [n], and
there are
p(p − 1) . . . (p − r)n(n − 1) . . . (n − k + r + 1) (5.9)
ways of doing this. Note that if we set p ∼ γn, then as n → ∞

p(p − 1) . . . (p − r)n(n − 1) . . . (n − k + r + 1) = pnk γ r (1 + O(1/n)), (5.10)

where f (x) = O(g(x)), with g(x) strictly positive for sufficiently large values of x, means |f (x)| ≤ Cg(x)
for some constant C > 0 and for all sufficiently large x.
Hence if G(r, k) denotes the number of different shapes of double trees with r + 1 i-vertices and
k − r j-vertices, we have
p k−1
1X k X r
E λi = γ G(r, k)(1 + O(1/n)). (5.11)
p i=1 r=0

We therefore need to show that   

1 k k−1
G(r, k) = (5.12)
r+1 r r
if we are to demonstrate that the spectral moments converge to (5.4).

30
 In order to do this we establish a bijection with a set of paths of 2k steps, or equivalently a set of
sequences, much as we did in the case of the semicircle law. In this case the paths (or sequences) are
generated as follows. Take (
{−1, 0} if m is odd
sm ∈ (5.13)
{0, 1} if m is even,
and take s2k = 0. We then consider all sequences of steps sm such that: (i) for any t = 1, 2, . . . , 2k
t
X
sm ≥ 0, (5.14)
m=1

that is, the path never crosses 0; (ii) #{m : sm = 1} = #{m : sm = −1} = r, that is, there are exactly
r positive steps and r negative steps; and hence (iii)
2k
X
sm = 0. (5.15)
m=1

To see that this set of sequences is in one-to-one correspondence with the set of double trees that
we wish to count, we start by constructing a mapping from sequences of this kind to double trees in
the following way. If m is odd and sm = 0 go from the i-vertex to a new j-vertex, and if sm = −1 go
back to the adjacent, previously-constructed j-vertex. If m is even and sm = 1, go from a j-vertex to a
new i-vertex, and if sm = 0 go back to an adjacent, previously-constructed i-vertex.
In the reverse direction, consider a path on the double tree as described above, starting at an
i-vertex. After an even number 2m of edge crossings one sits at an i-vertex. If this vertex is new, that
is, if it is being visited for the first time, then set s2m = +1, and otherwise set s2m = 0. After an odd
number 2m − 1 of edge crossings one must sit at a j-vertex, having just left an i-vertex. If this edge
crossing is the last exit from the i-vertex in question, set s2m−1 = −1, and otherwise set s2m−1 = 0. It
is easy to see that this construction gives sequences satisfying properties (ii) and (iii) set out above.
We need, therefore, to check (i). We do this by establishing a contradiction. Let us suppose we have
P2t−1 P2t−2
a first t such that m=1 sm = −1. This requires m=1 sm = 0 and s2t−1 = −1. Now, we could use
the sequence up to 2t − 2 to construct a double tree, and having s2t−1 = −1 would then mean not
establishing a new vertex, but instead going back to one previously created. This however, contradicts
already having built a double tree.
In order to evaluate G(r, k), we therefore need to count the number of sequences {sm }2k m=1 . We
begin by counting sequences that do not necessarily satisfy constraint (i). Recall that s2k = 0. Since
s2k 6= +1 , we have to count the number of ways of allocating r numbers +1 to k − 1 positions, and r
numbers −1 to
k
positions. The total number of sequences not necessarily satisfying constraint (i) is
therefore k−1 r
k
r .
We now count the number of sequences that fail to satisfy (i). For each of these there must exist
P2t−1
a first t such that m=1 sm = −1. Given this, construct a new sequence {s0m }2k m=1 in the following
way. Set s0m = sm for m = 1, 2, . . . , 2t − 1, and s02k = 0. For t ≤ j ≤ k − 1 put (s02j , s02j+1 ) =
(1, −1) if (s2j , s2j+1 ) = (1, −1), (s02j , s02j+1 ) = (0, 0) if (s2j , s2j+1 ) = (0, 0), (s02j , s02j+1 ) = (1, 0) if
(s2j , s2j+1 ) = (0, −1), and (s02j , s02j+1 ) = (0, −1) if (s2j , s2j+1 ) = (1, 0). The set of sequences {sm }2k m=1
that fail to satisfy (i) is in bijection with the set of sequences {s0m }2k m=1 , but the latter can be counted
straightforwardly to be the number of ways
ofk choosing

r − 1 out of k − 1 +1s, and r + 1 out of k −1s.
The number of ways of doing this is k−1 r−1 r+1 . Hence the number of sequences satisfying all of the
conditions, and hence G(r, k), is
        
k−1 k k−1 k 1 k k−1
− = . (5.16)
r r r−1 r+1 r+1 r r

31
5.2 Applications
Many applications of mathematics involve matrices that are rectangular, not square. In these cases
one is often interested in the singular values. The Marchenko-Pastur law describes the distribution of
singular values when the matrix entries are taken at random.
In data science, we have already seen that the empirical correlations between data points can be
represented as a matrix. These can be analysed through the spectrum of the sample correlation matrix
and so the Marchenko-Pastur law gives a ‘null hypothesis’ against which to test the spectrum, enabling
one to identify the directions (eigenvectors) associated with statistically significant correlations, in the
same way that one might test one-dimensional data against a normal distribution.
An example of this is in Principal Component Analysis, where one wants to identify the effective
dimensionality of the subspace representing the main variation in the data. A rule of thumb for how
many of the largest eigenvalues should be retained is to focus on those lying outside the support of the
Marchenko-Pastur distribution, or those near the edge if these is an unexpected density of them, as
these are most likely to represent non-random features.
For an overview covering a wide range of applications see [8], and for specific applications to
mathematical finance (where this is as an important technique), see [3, 9].
Finally, in many algorithms in numerical linear algebra, the rate of convergence is determined by
the largest singular values, and the Marchenko-Pastur law determines where these are expected to lie
(at the edge of the support).

32
6 Stieltjes and R-transforms
6.1 The Stieltjes transform and the semicircle law
Let µ be a non-negative finite measure on R.
Definition. The Stieltjes transform of µ is, for z ∈ C\R
Z
1
gµ (z) = dµ(x). (6.1)
R x − z
R
If all of the moments of µ, mk = R xk dµ(x), are bounded20 , then for sufficiently large z
Z ∞
1 1 X mk
gµ (z) = − dµ(x) = − (6.2)
R z(1 − x/z) z zk
k=0

and so gµ (z) can be thought of as the moment generating function of µ.

Note that when Imz = 6 0, it follows immediately from the fact that |(x − z)|−1 ≤ |Im(z)|−1 that
−1
|gµM (z)| ≤ |Im(z)| when µ is a probability measure. Moreover, when Imz > 0, gµ (z) is continuous
and analytic.
In some areas of mathematics and physics, −gµ (z) is referred to as the Green function.
When the measure is the empirical spectral distribution function µM (x) defined in (3.27),
1X 1 1 1
gµM (z) = = Tr . (6.3)
n j λj − z n M − zI

Note for future reference that it follows from this formula that if Imz > 0, then ImgµM (z) > 0. Also,
if µM is compactly supported, then as z → ∞
1
gµM (z) ∼ − . (6.4)
z
When µ is the measure corresponding to the semicircle law (3.3),
Z 2 Z Z 2π
1 1 p 1 π 2 2 1 2
gσ (z) = 4 − x dx =
2 sin ydy = sin2 ydy. (6.5)
2π −2 x − z π 0 2 cos y − z 2π 0 2 cos y − z

Setting w = eiy (and so w−1 = e−iy ) then gives

I
1 (w2 − 1)2
gσ (z) = − dw. (6.6)
4πi |w|=1 w2 (w2 + 1 − zw)
√
The integrand here has poles at w = 0, (z ± z 2 − 4)/2 (but only two lie inside the unit circle).
Evaluating the integral using Cauchy’s theorem gives
√
−z + z 2 − 4
gσ (z) = (6.7)
2
where the branch cut associated with the square root is chosen so that Imgσ (z) > 0 when Imz > 0.
Observe that gσ (z) satisfies the quadratic functional equation
1 1
gσ (z) = − − gσ (z)2 (6.8)
z z
with (6.7) being the solution satisfying Imgσ (z) > 0 when Imz > 0.
20 or, more generally, if mk = O(C k ) for some constant C > 0 when k → ∞

33
 (We remark in passing that for the Marchenko-Pastur distribution, the formula analogous to (6.8) is

γzg(z)2 + (z − 1 + γ)g(z) + 1 = 0, (6.9)

and the Stieltjes transform of the Marchenko-Pastur distribution is the solution of this equation
satisfying Img(z) > 0 when Imz > 0.)
Importantly, the Stieltjes transform can be inverted in the following way. Let us define a sub-
probability measure to be a measure µ on R such that µ(R) ≤ 1. For such a measure consider an
interval [a, b] such that µ{a} = µ{b} = 0. Then
Z b
1
µ[a, b] = lim Imgµ (x + iη)dx. (6.10)
η→0 a π
The proof of this inversion formula follows from interchanging the order of integration in
Z b Z b Z ∞
1 1 η
Imgµ (x + iη)dx = dµ(λ)dx. (6.11)
a π −∞ (λ − x) + η
2 2
a π

This gives
Z b Z ∞
1
Imgµ (x + iη)dx = V (λ)dµ(λ) (6.12)
a π −∞

where
      
1 b−λ a−λ 1 η(b − a)
V (λ) = arctan − arctan = arctan . (6.13)
π η η π (b − λ)(a − λ) + η 2

As η → 0, V (λ) → 1 for λ ∈ (a, b), V (a) = V (b) → 1/2, and V (λ) → 0 otherwise. Since V (λ) is
uniformly bounded and, for any η, V (λ) ∼ η/λ2 as λ → ∞, it can be majorized uniformly in η by a
positive integrable function. Hence it follows from the dominated convergence theorem that
Z ∞
V (λ)dµ(λ) → µ[a, b]. (6.14)
−∞

Note that this means that if two sub-probability measures µ and ν have gµ = gν , then µ = ν.
The Stieltjes transform provides an alternative route to proving both the Wigner semicircle law and
the Marchenko-Pastur law that avoids the combinatorics and path-counting of the proofs described
earlier in these notes; instead the Stieltjes transform approach is more rooted in analysis and probability
theory. We sketch here the proof of the semicircle law for a specific class of Wigner random matrices –
those with Gaussian entries. The proof can be extended to Wigner random matrices in general, but
this is an important class of examples and specialising will allow us to introduce new techniques that
provide some shortcuts.
We start with
1 1 1
gµM (z) = Tr := TrGM (z) (6.15)
n M − zI n
where GM is termed the resolvent of M . For this calculation we will specialise to the case where M is
a real symmetric matrix with entries that are i.i.d. Gaussian random variables with EMij = 0 and
2
EMij = 1/n. (Note the scaling of the variance, which is the one √ that leads to the semicircle law when
n → ∞ and which corresponds to scaling the eigenvalues with 1/ n with respect to when the variance
is 1.) Clearly
1
EgµM (z) = ETrGM (z). (6.16)
n
We now use the fact, which follows straightforwardly from the definition of the resolvent, that

GM +A (z) − GM (z) = −GM +A (z)AGM (z) (6.17)

34
and so, setting A = −M
1 1
GM (z) = − I + M GM (z). (6.18)
z z
Another consequence of (6.17) , which follows from setting, for a fixed pair i and j, Aij = Aji = 
and Auv = 0 otherwise and then taking the limit as  → 0, is that
∂(GM )uv
= −(GM )ui (GM )jv − (GM )uj (GM )iv . (6.19)
∂Mij

We note now that if x is a normal random variable with mean zero and variance σ 2 , then for any
function f that is differentiable and grows no faster than a polynomial, it follows from integrating by
parts that Exf (x) = σ 2 Ef 0 (x).
Combining the last results together, we then have that
1 1 1
ETrGM (z) = − + ETr(M GM )
n z nz
1 1 X
=− + E (Mij (GM )ji )
z nz ij
 
1 1 X ∂(GM )ji
=− + 2 E
z n z ij ∂Mij
1 1 X
=− + 2 E [−(GM )ji (GM )ij − (GM )jj (GM )ii ]
z n z ij
" 2 #
1 1  2
 1 1
= − − 2 E TrGM − E TrGM . (6.20)
z n z z n

Now note that

n
1  2
 1 1X 1
2
TrGM =
n z n|z| n i=1 |λi − z|2
1
≤
n|z||Imz|2
→0 (6.21)

when n → ∞.
Hence we have that
1 1

EgµM (z) = − − E gµM (z)2 + E(z, n) (6.22)
z z
where for any z, E(z, n) → 0 as n → ∞.
2
We now want to establish that we can replace E gµM (z)2 with (EgµM (z)) , with an error that
tends to zero as n → ∞. Note that gµM (z) depends on n, because M is an n × n matrix. In the
following calculation the n-dependence will be important, so we record that dependence explicitly,
augmenting the notation to read gµMn (z). The strategy will be to show that replacing Mn by an
(n − 1) × (n − 1) minor, Mn−1 , does not change the Stieltjes transform appreciably when n is large –
i.e. that gµMn−1 (z) is close to gµMn (z). This will allow us to use an estimate to establish that gµMn (z)


lies sufficiently close to E gµMn (z) for us to prove what we need.
(n)
Let z = a + ib, with b > 0. We denote the eigenvalues of Mn by {λj }nj=1 and the eigenvalues of
(n−1) n−1
Mn−1 by {λj }j=1 . Now, Cauchy’s Interlace Theorem21 implies that the eigenvalues of Mn and
21 Cauchy’s Interlace Theorem may be stated as follows. Let A be an n × n Hermitian matrix with eigenvalues

α1 ≤ α2 ≤ · · · ≤ αn , and let B be an (n − 1) × (n − 1) principal submatrix of A with eigenvalues β1 ≤ β2 ≤ · · · ≤ βn−1 .

Then α1 ≤ β1 ≤ α2 ≤ β2 ≤ · · · ≤ αn−1 ≤ βn−1 ≤ αn .

35
 (n) (n−1) (n)
those of Mn−1 interlace, that is λj ≤ λj ≤ λj+1 for j = 1, . . . , n − 1. Hence the difference
n−1
X n
X
b b
(n−1)
− (n)
(6.23)
j=1 (λj − a)2 + b2 j=1 (λj − a)2 + b2
(n) (n−1) (n) (n−1) (n−1) (n)
is bounded in n , because (x−a)b2 +b2 has finite total variation and λ1 , λ1 , λ2 , λ2 , . . . , λn−1 , λn
forms a partition. The same conclusion holds if one replaces the numerators in the summands with
(n−1) (n)
λj − a and λj − a respectively. It follows that
 
1
gµMn (z) = gµMn−1 (z) + O . (6.24)
n
gµMn (z) is therefore said to be stable in n. Note as well that the right-hand side of this equation
depends only on the top left (n − 1) × (n − 1) minor of M and is independent of its nth row and nth
column. We can therefore apply the following inequality, due to McDiarmid.
Lemma 6 (McDiarmid’s inequality). Let x1 , . . . , xn be independent random variables taking values in
ranges R1 , . . . , Rn and let F : R1 × · · · × Rn → C be a function having bounded differences, so that,
there exist constants c1 , . . . , cn such that for all i,

|F (x1 , . . . , xn ) − F (x1 , . . . , xi−1 , x0i , xi+1 , . . . , xn )| ≤ ci (6.25)

for all xi , x0i ∈ Ri . Then for any κ > 0

P (|F (x1 , . . . , xn ) − EF (x1 , . . . , xn )| ≥ κσ) ≤ C exp(−cκ2 ) (6.26)

Pn
for some C > 0 and c > 0, independent of F and n, where σ 2 = i=1 c2i .
We can apply this inequality in the following way. First note that (6.24) still applies if we resample the
nth row and nth column. Denoting the Stieltjes transform of the resampled matrix by gµM 0 (z), it follows
from (6.24) that gµM (z) = gµM 0 (z) + O(1/n). Moreover, this applies to resampling any row and column
with the same label. Therefore, applying McDiarmid’s inequality with xj = (Mj,j , Mj,j+1 , . . . , Mj,n )
gives √
P(|gµMn (z) − E(gµMn (z))| ≥ κ/ n) ≤ C exp(−cκ2 ). (6.27)


This means that as n → ∞ gµMn (z) lies sufficiently close to E gµMn (z) that we can replace

2
E gµM (z)2 with (EgµM (z)) in (6.22). Hence we have that
1 1 2
EgµM (z) = − − (EgµM (z)) + E 0 (z, n) (6.28)
z z
where E 0 (z, n) → 0 as n → ∞. Hence EgµM (z) converges to the solution of
1 1
s(z) = − − s(z)2 (6.29)
z z
which we have already seen is the Stieltjes transform of the semicircle distribution, once one specifies
that Ims(z) > 0 when Imz > 0. Therefore, because the Stieltjes transform is invertible (as a continuous
operator), we see that the empirical spectral measure converges weakly to the semicicle distribution
when n → ∞ in expectation.
Finally, taking κ = n1/3 in (6.27) and applying the Borel-Cantelli lemma, we see that almost surely

gµMn (z) − E[gµMn (z)] ≤ O(n−1/6 ) (6.30)

for all large n. Therefore gµMn (z) − E[gµMn (z)] → 0 almost surely for all z with Imz > 0.
Note that one can view (6.27) and (6.30) as examples of concentration of measure in the context of
the Stieltjes transform.

36
6.2 The R-transform
As noted previously, −gµ (z) is sometimes referred to as the Green function. The function Bµ (z) that is
the formal inverse of the Green function, i.e. the function which satisfies

Bµ (−gµ (z)) = z, (6.31)

is sometimes referred to as the Blue function.

The R-transform of µ is then defined to be
1
Rµ (z) = Bµ (z) − . (6.32)
z
The R-transform of the semicircle distribution is easily seen (e.g. from the quadratic equation
satisfied by the Stieltjes transform of the semicircle distribution) to be

Rsc (z) = z (6.33)

and similarly for the Marchenko-Pastur distribution

1
RM P (z) = . (6.34)
1 − γz
The fact that the R-transforms of the semicircle distribution and the Marchenko-Pastur distributions
take such elementary forms is in indication that mathematically they are rather natural and simple, in
a way that may be obscured by other representations. Indeed, the semicircle distribution can be viewed
as playing a similar role to that played by the normal distribution for commuting random variables.
The R-transform plays an important role in analysing the spectral density of sums of random
matrices which satisfy certain natural conditions, but we will not examine this further here.

37
7 Eigenvalue statistics for GOE, GUE, CUE
In this section we will focus on matrices from the GOE, the GUE, and the CUE. We start with the
GOE and the GUE and turn to the CUE at the end of the section.
We have seen that GOE and GUE matrices satisfy the semicircle law, so that if we take the
eigenvalues of a GOE or a GUE matrix√(in which the matrix elements have a variance independent of
the matrix size n) and scale them by 1/ n then the support of the spectrum is almost surely [−2, 2] in
the limit n → ∞ with a density in that interval given by the semicircle distribution.
The question we wish to address now is: how are the eigenvalues distributed on the scale of their
mean spacing? This is a much smaller scale than that for which the semicircle law holds. It might be
assumed that if the matrix elements are i.i.d. random variables, then the eigenvalues will also behave
like i.i.d. random variables on this scale; that is, that they behave like a Poisson point process. We
shall see, however, that this is not the case, and that instead the eigenvalues are strongly correlated in
an interesting and significant way. This is an important feature of random matrices.

7.1 Poisson point processes

A Poisson point process on R with parameter λ > 0 corresponds to a set of random points such that
the probability that the number of points in an interval (a, b], N (a, b), is exactly equal to k is

[λ(b − a)]k
P(N (a, b) = k) = exp[−λ(b − a)]. (7.1)
k!
Moreover, the numbers of points in any finite collection of disjoint intervals are independent of each
other.
The probability density PPoiss (s) that the spacing between two nearest-neighbouring points is s may
be computed as follows. Let us fix a point at one end of an interval of length s. PPoiss (s)ds is then the
probability that the interval is empty, and that there is a point in the next infinitesimal interval of
length ds. This is
PPoiss (s)ds = exp[−λs]λds. (7.2)
Clearly Z ∞
PPoiss (s)ds = 1. (7.3)
0
If now we choose to set the mean spacing to be 1, i.e.
Z ∞
sPPoiss (s)ds = 1, (7.4)
0

this corresponds to setting λ = 1, in which case PPoiss (s) = e−s . Hence the probability to find
nearest-neighbouring random points a distance s apart is maximal at s = 0 and decreases as s increases.
We shall later want to compare this to how the eigenvalues of GOE and GUE matrices are distributed.

7.2 Joint eigenvalue density

We have seen earlier that the GOE and GUE measures can be expressed, up to a normalization constant,
in the form β 2
e− 4 Tr(M ) dM (7.5)
where for the GOE β = 1, for the GUE β = 2, and dM is the Lebesgue measure on the appropriate
space of matrix entries. Our aim now is to re-express this in terms of the Lebesgue measure on the
space of the eigenvalues. Specifically, denoting the eigenvalues of a matrix M by λ1 , λ2 , . . . , λn , define
a class function of M to be a function, f (M ) = f (λ1 , . . . , λn ), that depends only on the eigenvalues
and which is symmetric with respect to permutations of its arguments. The following theorem gives a
formula for averaging class functions with respect to the GUE and the GOE.

38
Theorem 7. The average of f (λ1 , . . . , λn ) over the GUE and GOE can be expressed as the multiple
integral
Z ∞ Z ∞ Y β P 2
E[f (λ1 , λ2 , . . . , λn )] = c(β)
n · · · f (λ1 , . . . , λn ) |λi − λj |β e− 4 i λi dλ1 . . . dλn (7.6)
−∞ −∞ i<j

(1) (2)
where β = 1 for the GOE, β = 2 for the GUE, and cn and cn are normalization constants.
We remark that the normalization constants can be evaluated to be:
n
Y
1
(β)
= (2π)n/2 β −n/2−βn(n−1)/4 [Γ(1 + β/2)]−n Γ(1 + βj/2) (7.7)
cn j=1

although we shall not make explicit use of this fact and do not include the proof here.
In the case of the GUE, the idea behind the proof of Theorem 7 is based on the observation that
any Hermitian matrix M can be diagonalized by a unitary transformation; that is, we can write
M = U DU † , where D is diagonal and U is unitary. The idea is then to change variables from the
matrix elements of M to the matrix elements of U and D. Because the GUE measure is invariant
under unitary transformations, the matrix elements of U can be integrated out trivially, leaving just
the matrix elements of D, which are the eigenvalues λ1 , . . . , λn .
We first note that when the eigenvalues of M are distinct, the decomposition M = U DU † is
unique up to (i) permuting the eigenvalues, and (ii) multiplying U by a diagonal matrix with entries
eiθ1 , eiθ2 , . . . , eiθn , and that this holds almost surely in the GUE (and the GOE). To see this, observe
first that the joint distribution of the entries of M has a smooth density. Since the eigenvalues are
continuous functions of the matrix entries, and because the event of two eigenvalues coinciding has
codimension ≥ 1 in the n-dimensional variety of eigenvalues, it follows that the eigenvalues are almost
surely distinct 22 . So we can ignore this situation.
Let y1 , . . . , yn , w1 , . . . , wn(n−1)
 be the local parameters on (Rn /Sn ) × (U (n)/Tn ). Our goal now is
∂Mij ∂Mij 2
to compute the Jacobian det ∂yα , ∂wβ . To this end we write M as an element of Rn , setting
 
M11 M22 M
φ(M ) = √ , √ , . . . , √nn , ReM12 , ImM12 , ReM13 , ImM13 , . . . , ReMn−1,n , ImMn−1,n . (7.8)
2 2 2
It is easily seen that
1 1
|φ(M )|2 = Tr[M 2 ] = Tr[(U M U † )2 ] = |φ(U M U † )|2 (7.9)
2 2
for any n × n unitary matrix U , because the trace is invariant under unitary conjugation. Consequently,
2 2
the linear transformation τU : Rn → Rn , τU (y) = φ(U † φ−1 (y)U ) is isometric; that is det τU = 1.
Hence     
∂Mij ∂Mij ∂Mij ∂Mij
det , = det τU , . (7.10)
∂yα ∂wβ ∂yα ∂wβ
Clearly we can also think of τU as acting on matrices X according to τU (X) = τU (φ(X)). Now,
   
∂M ∂D † ∂D
τU = τU U U = . (7.11)
∂yi ∂yi ∂yi
22 An alternative proof of this uses the fact that the zero set of any polynomial in k variables has zero Lebesgue measure

in Rk . There are degenerate eigenvalues if and only if the characteristic polynomial of M , det(M − xI), and its derivative
have a common zero. A necessary and sufficient condition for this is that the discriminant i<j (λi − λj )2 vanishes. The
Q
result then follows from the fact that discriminant is itself a polynomial in the entries of M .

39
Note that acting with φ on the right-hand side gives a vector with all n2 coordinates except the ith
equal to zero. We also have that
   
∂M ∂U ∂U † ∂U ∂U †
τU = τU DU † + U D = U† D+D U. (7.12)
∂wβ ∂wβ ∂wβ ∂wβ ∂wβ

Since U † U = I, we have, differentiating with respect to wβ ,

∂U ∂U †
U† + U = 0. (7.13)
∂wβ ∂wβ

Therefore, setting Tβ = U † ∂w
∂U
β  
∂M
τU = Tβ D − DTβ . (7.14)
∂wβ
Bearing in mind that D is diagonal,
  
∂M
τU = (Tβ )ij (yj − yi ). (7.15)
∂wβ ij
 
∂Mij ∂Mij
Therefore, the matrix τU ∂yα , ∂wβ has the following form
 
In 0n 0n 0n ...
0n Re(T1 )12 (y2 − y1 ) Im(T1 )12 (y2 − y1 ) Re(T1 )13 (y3 − y1 ) . . .
 
0n Re(T2 )12 (y2 − y1 ) Im(T2 )12 (y2 − y1 ) Re(T2 )13 (y3 − y1 ) . . . (7.16)
 
.. .. .. .. ..
. . . . .

where 0n is the n × n zero matrix.

Hence, computing the determinant, we have that
 
  Re(T1 )12 Im(T1 )12 ...
∂Mij ∂Mij Y  
det , = (yj − yi )2 det Re(T2 )12 Im(T2 )12 ... . (7.17)
∂yα ∂wβ i<j
.. .. .. ..
. . . .
Q
The point is that the problem has now factorized and we are left with a product of i<j (yj − yi )2 and
the determinant, which is a function only of the w-variables. Hence in computing E[f (λ1 , λ2 , . . . , λn )],
where the expectation is computed with respect to the GUE, we can integrate with respect to all of the
w-variables,
Q obtaining some constant, and we are left only with the integrals over y1 , . . . , yn , picking up
a factor i<j (yj − yi )2 . This proves the theorem for the GUE.
The proof for the GOE follows an identical path, except that in this case the matrix M is diagonalized
by an orthogonal transformation, so M = ODOT . The number of independent random variables in M
(e.g. the matrix elements on and above the diagonal) is then n(n + 1)/2. These can be replaced by the
diagonal elements of D (the eigenvalues), y1 , . . . , yn and n(n − 1)/2 
other parameters
 w1 , . . . , wn(n−1)/2 .
∂M ∂M
One again needs to compute the Jacobian of the transformation det ∂yαij , ∂wβij . The calculation goes
through essentially as above (but with orthogonal matrices replacing the unitary matrices), and one
ends up with a Jabobian which is similar to (7.16) but missing the columns involving Im(Tβ ) because
the orthogonal matrices can be taken to be real valued. Therefore Q each factor (yj − yi ) appears only
once and the Jacobian in this case factorizes as the product of i<j |yj − yi | and a function only of
the w-variables. Hence again in computing E[f (λ1 , λ2 , . . . , λn )], where the expectation is computed
now with respect to the GOE, we can integrate with respect to all of the w-variables, Q obtaining some
constant, and we are left only with the integrals over y1 , . . . , yn , picking up a factor i<j |yj − yi |. This
proves the theorem for the GOE.

40
 We turn next to the Circular Unitary Ensemble of random unitary matrices M . In this case too
the matrices in question can be diagonalized by a unitary transformation: M = U DU † , where D is
diagonal with entries eiθ1 , eiθ2 , . . . , eiθn . The ensemble itself is, by construction, invariant under unitary
transformations. A calculation that is essentially the same as for the GUE then gives, for a class
function f (θ1 , . . . , θn ), that
Theorem 8. The average of f (θ1 , . . . , θn ) over the CUE can be expressed as the multiple integral
Z 2π Z 2π Y
1
E[f (θ1 , . . . , θn )] = n
··· f (θ1 , . . . , θn ) |eiθj − eiθk |2 dθ1 . . . dθn . (7.18)
(2π) n! 0 0 j<k

This is known as the Weyl integration formula.

Finally, for the Wishart Ensemble, the corresponding result is23
Theorem 9. The average of f (λ1 , . . . , λp ) over the Wishart Ensemble can be expressed as the multiple
integral
Z ∞ Z ∞ Y Y (n−p−1)/2 1
E[f (λ1 , λ2 , . . . , λp )] = c(Wishart)
n,p ··· f (λ1 , . . . , λp ) |λi − λj | λk e− 2 λk dλ1 . . . dλp
0 0 i<j k
(7.19)
(Wishart)
where cn,p is a normalization constant.
One of the most important consequences of Theorem 7, Theorem 8 and Theorem 9 is that the
eigenvalues of random matrices repel each other. It is clear from Theorem 7 that the probability that
λi − λj ∈ (s, s + ds) vanishes like sβ as s → 0; so it vanishes quadratically for the GUE and linearly for
the GOE. The same repulsion phenomenon may be seen in Theorem 8 and Theorem 9. Note that this
is unlike the Poisson point process, where the probability to find nearest-neighbouring random points a
distance s apart is maximal at s = 0.

7.3 The method of orthogonal polynomials for the GUE

Q
For the GUE, the integrand in Theorem 7 involves i<j |λi − λj |2 . It is a key observation that the latter
factor can be written as a determinant. Specifically, it is the square of the Vandermonde determinant:
2
1 1 ... 1
λ1 λ2 ... λn
Y
|λi − λj |2 = λ21 λ22 ... λ2n . (7.20)
i<j
.. .. .. ..
. . . .
λn−1
1 λn−1
2 ... λn−1
n

(To see this, note that the determinant is a homogeneous polynomial in the variables λ1Q , λ2 , . . . , λn of
degree n(n − 1)/2. It vanishes if any pair of the variables are set equal, so is divisible by i<j (λj − λi ).
This is also a homogeneous polynomial of degree n(n − 1)/2. Therefore the two polynomials are equal up
to a constant multiplier. That this constant has to be one may be seen by comparing any corresponding
terms on the two sides; for example
Q the product of terms on the diagonal of the matrix corresponds to
multiplying the first terms in i<j (λj − λi ).)
It is a second key observation that by combining rows in the determinant we also have
2
p0 (λ1 ) p0 (λ2 ) ... p0 (λn )
p1 (λ1 ) p1 (λ2 ) ... p1 (λn )
Y
|λi − λj |2 = p2 (λ1 ) p2 (λ2 ) ... p2 (λn ) (7.21)
i<j
.. .. .. ..
. . . .
pn−1 (λ1 ) pn−1 (λ2 ) . . . pn−1 (λn )
23 This
Q
is the formula for real Wishart matrices. For complex Wishart matrices the factor i<j |λi − λj | is squared.

41
where pk is a monic polynomial of degree k. Hence the integrand in Theorem 7 is the function we wish
to average, f (λ1 , λ2 , . . . , λn ), multiplied by
1 2 1 2 1 2 2
p0 (λ1 )e− 4 λ1 p0 (λ2 )e− 4 λ2 ... p0 (λn )e− 4 λn
1 2 1 2 1 2
p1 (λ1 )e− 4 λ1 p1 (λ2 )e− 4 λ2 ... p1 (λn )e− 4 λn
1 2 1 2 1 2
p2 (λ1 )e− 4 λ1 p2 (λ2 )e− 4 λ2 ... p2 (λn )e− 4 λn . (7.22)
.. .. .. ..
. . . .
1 2 1 2 1 2
pn−1 (λ1 )e− 4 λ1 pn−1 (λ2 )e− 4 λ2 ... pn−1 (λn )e− 4 λn

This is true for any monic polynomials pk . There is, however, one choice that considerably simplifies
the calculations we shall want to perform: we henceforth take {pk (x)} to be orthogonal with respect to
2
the measure e−x /2 dx. Specifically, setting
1 2 1 2 1 2 2
p0 (λ1 )e− 4 λ1 p0 (λ2 )e− 4 λ2 ... p0 (λn )e− 4 λn
1 2 1 2 1 2
p1 (λ1 )e− 4 λ1 p1 (λ2 )e− 4 λ2 ... p1 (λn )e− 4 λn
1 2 1 2 1 2
ρ(λ1 , . . . , λn ) = c(2)
n
p2 (λ1 )e− 4 λ1 p2 (λ2 )e− 4 λ2 ... p2 (λn )e− 4 λn (7.23)
.. .. .. ..
. . . .
1 2 1 2 1 2
pn−1 (λ1 )e− 4 λ1 pn−1 (λ2 )e− 4 λ2 ... pn−1 (λn )e− 4 λn

we have !
n−1
Y

ρ(λ1 , . . . , λn ) = c(2)
n κ2i det AAT (7.24)
i=1

where A is the n × n matrix with elements

1 1 2
Aij = pj−1 (λi )e− 4 λi (7.25)
κj−1

and the κi are constants satisfying

Z ∞
2
pi (x)pj (x)e−x /2
dx = κ2i δij . (7.26)
−∞

This clearly holds for all j < i if Z ∞

2
xk pi (x)e−x /2
dx = 0 (7.27)
−∞
2
for all k < i. We can arrange this by taking pi (x)e−x /2 to be an exact ith derivative of some function
that vanishes as |x| → ∞, and the obvious choice is to set

2 di −x2 /2
pi (x) = (−1)i ex /2
e . (7.28)
dxi
These are known as the Hermite polynomials and are usually denoted Hi (x)24 . We then also have
Z ∞ Z ∞ √
−x2 /2 2
2
κi = 2
pi (x)e dx = xi pi (x)e−x /2 dx = i! 2π. (7.29)
−∞ −∞

We shall use pk (x) and Hk (x) interchangeably, since they are equal.
2
24 In some references the Hermite polynomials are defined to be orthogonal with respect to the measure e−x /2 dx, and
2
in others with respect to the measure e−x dx. We choose the former convention. The latter convention is more common
in the physics literature.

42
 The normalized Hermite polynomials κ−1
k Hk (x) form an orthonormal set of functions with respect
−x2 /2
to the weight e : Z ∞
−x2 /2
κ−1 −1
j Hj (x)κk Hk (x)e dx = δjk . (7.30)
−∞

They form an orthonormal basis of the (complex) Hilbert space of functions satisfying
Z ∞
2
|f (x)|2 e−x /2 dx < ∞ (7.31)
−∞

in which the inner product is taken to be

Z ∞
2
hf, gi = f (x)g(x)e−x /2
dx. (7.32)
−∞

2
The functions φk (x) := κ−1
k Hk (x)e
−x /4
are called the Hermite functions. They satisfy the
differential equation
d2 Ψ 1 2 1
− + x Ψ = (k + )Ψ (7.33)
dx2 4 2
which is the Schrödinger equation of the simple harmonic oscillator.
We note now that the entries of the matrix AT A appearing in (7.24) can be written
n
X n−1
X


AAT ij
= Aik Ajk = φk (λi )φk (λj ). (7.34)
k=1 k=0

The function
n−1
X n−1
1 X 1 2 2
Kn (x, y) = φk (x)φk (y) = √ pk (x)pk (y)e−(x +y )/4 (7.35)
k=0
2π k=0 k!
is known as the Christoffel-Darboux kernel in the theory of orthogonal polynomials. It is the kernel of
the projection onto the span of φ1 , . . . , φn−1 .
We note for later use the following general properties of the Hermite polynomials.
Lemma 10. The Hermite polynomials satisfy

pk+1 (x) = xpk (x) − p0k (x), (7.36)

xpk (x) = pk+1 (x) + kpk−1 (x), (7.37)

p00k (x) − xp0k (x) = −kpk (x), (7.38)
and for x 6= y
n−1
X pk (x)pk (y) pn (x)pn−1 (y) − pn−1 (x)pn (y)
= . (7.39)
k! (n − 1)!(x − y)
k=0

Proof. (7.36) follows from differentiating (7.28). (7.37) follows from expanding xpk (x) in terms of
the basis formed by all of the Hermite polynomials and using (7.27). (7.38) follows by differentiating
(7.36) and then using (7.37). (7.39) follows from multiplying Kn by x − y and then using (7.37); it is
a special case of a general formula known as the Christoffel-Darboux Theorem, which goes back to
Christoffel in 1858 (for more on this see, for example, [16]).
Note that, applying L’Hôpital’s rule, it follows from (7.39) that
n−1
X p2k (x) p0 (x)pn−1 (x) − p0n−1 (x)pn (x)
= n . (7.40)
k! (n − 1)!
k=0

43
 We have from (7.39) that for x 6= y
√ φn (x)φn−1 (y) − φn−1 (x)φn (y)
Kn (x, y) = n (7.41)
x−y
and that √

Kn (x, x) = n φ0n (x)φn−1 (x) − φ0n−1 (x)φn (x) . (7.42)
In terms of the Christoffel-Darboux kernel, we have that the right-hand side of (7.24) can be written
as
const. det (Kn (λi , λj )) . (7.43)
It follows from the orthonormality of the Hermite polynomials that
Z ∞
Kn (x, y)Kn (y, z)dy = Kn (x, z). (7.44)
−∞

This allows us to apply the following general lemma

Lemma 11 (Gaudin’s lemma). Let
Jn = (Jij )1≤i,j≤n = (k(xi , xj ))1≤i,j≤n (7.45)
where k(x, y) is a kernel satisfying
Z ∞
k(x, y)k(y, z)dy = k(x, z). (7.46)
−∞

Then Z ∞
det(Jn )dxn = (r − n + 1) det(Jn−1 ) (7.47)
−∞
R∞
where r = −∞
k(x, x)dx.
Proof. We note first that expanding the determinant gives
Z ∞ Z ∞ X
det(Jn )dxn = sgn(σ)k(x1 , xσ(1) ) . . . k(xn , xσ(n) )dxn
−∞ −∞ σ∈S
n
Z ∞ Xn X
= sgn(σ)k(x1 , xσ(1) ) . . . k(xn , xm )dxn . (7.48)
−∞ m=1 σ:σ(n)=m

In the term m = n, σ runs over all permutations in Sn−1 and so the integral evaluates to r det(Jn−1 ),
because the sign of σ as a permutation in Sn−1 is the same as it is in Sn .
When m < n, let j = σ −1 (n), and take σ̂ ∈ Sn−1 to be given by
(
σ(i), if i 6= j
σ̂(i) = (7.49)
m, if i = j.
The map {σ ∈ Sn : σ(n) = m} −→ Sn−1 , σ 7→ σ̂, is straightforwardly seen to be a bijection, and
sgn(σ̂) = −sgn(σ), because the two permutations differ by a transposition (mn). Therefore
Z ∞ X
sgn(σ)k(x1 , xσ(1) ) . . . k(xn , xm )dxn
−∞ σ:σ(n)=m
Z ∞ X
= sgn(σ)k(x1 , xσ(1) ) . . . k(xn−1 , xσ(n−1) )k(xj , xn )k(xn , xm )dxn
−∞ σ:σ(n)=m
X
= −sgn(σ̂)k(x1 , xσ̂(1) ) . . . k(xn−1 , xσ̂(n−1) )
σ̂∈Sn−1

= − det Jn−1 . (7.50)

44
Finally, using the fact that this last equality holds independently of the value of m, the equation in the
statement of the theorem is established.

We now return to Theorem 7. Note that in the notation introduced above, we have for the GUE
that Z ∞ Z ∞
E[f (λ1 , λ2 , . . . , λn )] = ··· f (λ1 , . . . , λn )ρ(λ1 , . . . , λn )dλ1 . . . dλn (7.51)
−∞ −∞

with
ρ(λ1 , . . . , λn ) = c̃n det (Kn (λi , λj )1≤i,j≤n ) (7.52)
where c̃n denotes the constant.
In many situations we shall see that we will want to average functions f that depend only on k < n
eigenvalues. In this case we need to evaluate the integral of ρ over the remaining n − k variables. This
can be achieved using (Gaudin’s) lemma 11.
Lemma 12. For the GUE we have
Z ∞ Z ∞
(n − k)!
··· ρ(λ1 , . . . , λn )dλk+1 . . . dλn = det (Kn (λi , λj )1≤i,j≤k ) . (7.53)
−∞ −∞ n!

Proof. Applying Gaudin’s lemma, and using the fact that r = n for this kernel, because the Hermite
functions are orthonormal, we have immediately that
Z ∞
ρ(λ1 , . . . , λn )dλn = c̃n det (Kn (λi , λj )1≤i,j≤n−1 ) . (7.54)
−∞

Applying Gaudin’s lemma a second time now gives

Z ∞Z ∞
ρ(λ1 , . . . , λn )dλn−1 dλn = 2c̃n det (Kn (λi , λj )1≤i,j≤n−2 ) . (7.55)
−∞ −∞

Repeated applications of Gaudin’s lemma therefore give

Z ∞ Z ∞
··· ρ(λ1 , . . . , λn )dλk+1 . . . dλn = (n − k)!c̃n det (Kn (λi , λj )1≤i,j≤k ) . (7.56)
−∞ −∞

In the case when k = 0, this reduces to

Z ∞ Z ∞
··· ρ(λ1 , . . . , λn )dλ1 . . . dλn = n!c̃n (7.57)
−∞ −∞

and so we see that c̃n = 1/n!.

It is worth pointing out that this is a rather remarkable result: it asserts that integrating an n × n
determinant, involving a kernel evaluated at all pairs drawn from n variables, over n − k of those
variables, gives a k × k determinant involving the same kernel. Point processes with this property are
called determinantal point processes.

7.4 Correlation functions

Definition. The k-point correlation function is defined to be
Z ∞ Z ∞
n!
Rk (x1 , . . . , xk ) := ··· ρ(x1 , . . . , xn )dxk+1 . . . dxn (7.58)
(n − k)! −∞ −∞
= det (Kn (xi , xj )1≤i,j≤k ) . (7.59)

45
 These correlation functions have the following interpretation. For any Borel set B,
Z
R1 (x)dx = E[#{eigenvalues in B}] (7.60)
B

and Z Z
R2 (x, y)dxdy = E[#{ordered pairs of eigenvalues in B}] (7.61)
B B
etc. √
Clearly if we rescale the eigenvalues by 1/ n and normalise the 1-point correlation function by 1/n,
we expect to recover the semicircle law from the 1-point correlation function; that is, we expect
1 √ 1 √ √ 1 p
lim √ R1 ( nx) = lim √ Kn ( nx, nx) = 4 − x2 (7.62)
n→∞ n n→∞ n 2π

for x ∈ [−2, 2], and for the limit to give 0 if |x| > 2. This can be shown to be true using asymptotic
formulae for the Hermite polynomials when n → ∞, leading to yet another proof of the semicircle law
for the GUE.
Specifically, let us define √
1
Ξn+p (x) = n 4 φn+p (x n) (7.63)
for p = −2, −1, 0. It follows from the properties of the Hermite polynomials listed above that
r
1 √ √ 2 n−1
√ Kn ( nx, nx) = Ξn−1 (x) − Ξn−2 (x)Ξn (x). (7.64)
n n

The asymptotic formulae we need25 are as follows. First, setting x = 2 cos φ with 0 < φ < π, when
n→∞      
1 1 1 π
Ξn+p (x) ∼ √ cos n φ − sin(2φ) + p + φ− (7.65)
π sin φ 2 2 4
uniformly for φ in a compact subset of (0, π). Second, setting |x| = 2 cosh φ with 0 < φ,

e(p+1/2)φ − n (e2φ +1−2φ)

Ξn+p (x) ∼ √ e 2 (7.66)
2π sinh φ

uniformly for φ in a compact subset of (0, ∞)26 . Substituting (7.65) into (7.64) gives the semicircle law
in |x| < 2. Since the semicircle law has a total mass of one in (−2, 2), there can be no limiting mass at
x = 2 or in |x| > 2. When |x| > 2 this also follows by substituting (7.66) into (7.64).
Note that (7.42) gives that
1 √ √ √ √ √
√ R1 ( nx) = φ0n ( nx)φn−1 ( nx) − φ0n−1 ( nx)φn ( nx) (7.67)
n

which gives an exact formula for the expectation of the empirical spectral density of a GUE matrix in
terms of the Hermite functions. We illustrate this by showing in Figure 7 the results of a numerical
experiment. This involved generating 10,000 random matrices from the GUE with n = 10, plotting a
histogram of the eigenvalues, and comparing with the exact formula (7.67). For ease of visualisation,
the eigenvalues have each been divided by 2, so the support of the rescaled semicircle that emerges
when n → ∞ is |x| ≤ 1.
√The semicircle law describes the limiting eigenvalue distribution when the eigenvalues are scaled by
1/ n. There are then n scaled eigenvalues lying between -2 and 2, and so their mean separation is
25 These go under the general name Plancherel-Rotach asymptotic formule and can be established via a WKB analysis

of (7.33).
26 Note that when φ > 0, e2φ + 1 − 2φ > 0.

46
0.6

0.5

0.4

0.3

0.2

0.1

0.0
-1.0 -0.5 0.0 0.5 1.0

Figure 7: The spectral density of 10,000 GUE matrices of dimension 10 compared to (7.67).

of order 1/n. We now want to look at correlations on the scale of the mean eigenvalue spacing. We
therefore need to scale the √ √ so that they have constant mean spacing. There are n eigenvalues
eigenvalues √
lying roughly between −2 n and √ 2 n. They therefore have a mean separation that scales
√ like 1/ n
and so we need to rescale them by n to achieve a constant mean spacing. Setting yi = nλi , we need
to consider the n → ∞ asymptotics of
h √ √ ir h ir
det n−1/2 Kn (yi / n, yj / n) := det K̂n (yi , yj ) (7.68)
i,j=1 i,j=1

where, when yi 6= yj ,
√ √ √ √
φn (yi / n)φn−1 (yj / n) − φn−1 (yi / n)φn (yj / n)
K̂n (yi , yj ) = , (7.69)
yi − yj
and √ √ √ √
K̂n (y, y) = φ0n (y/ n)φn−1 (y/ n) − φ0n−1 (y/ n)φn (y/ n). (7.70)
The important result that we now establish is the following one.
Theorem 13. For any fixed x, y
sin(π(x − y))
lim K̂n (x, y) = . (7.71)
n→∞ π(x − y)
This is an immediate consequence of the following lemma.
Lemma 14. For any fixed x and any fixed l
√
lim n1/4 φn−l (x/ n) − π −1/2 cos (x − (n − l)π/2) = 0. (7.72)
n→∞

47
 Proof. Let us set m = n − l. Note that
m Z ∞
m x2 /2 dm −x2 /2 m x2 /2 d 2 dz
(−1) e m
e = (−1) e m
e−ixz−z /2 √
dx dx −∞ 2π
Z
1 x2 /2 ∞ 2
=√ e (iz)m e−ixz−z /2 dz (7.73)
2π −∞

and so Z ∞
√ 2 √
n−z 2 /2
n1/4 φm (x/ n) = (2π)−3/4 (m!)−1/2 ex /4n n1/4 (iz)m e−ixz/ dz. (7.74)
−∞

Using Stirling’s formula to evaluate the prefactor and Laplace’s method27 to evaluate the integral when
n → ∞ proves the lemma. (See, for example, Chapter VIII in [16].)

Substituting the asymptotic estimate in the lemma into the formula for K̂n , and using (7.36) and
(7.37) to evaluate the asymptotics of the derivative of φ, then leads to the formula in the theorem.
The theorem means that in the limit when n → ∞ the k-point correlation function of the scaled
eigenvalues at the centre of the spectrum is simply given by the k × k determinant involving the kernel
sin(π(x − y))
. (7.75)
π(x − y)
This is known as the sine kernel. For example, the pair correlation of the scaled eigenvalues is simply
 2
sin(π(x − y))
1− . (7.76)
π(x − y)

So, for example, for f (x, y) such that the sum and integral converge,

X Z ∞ Z ∞  2 !
√ √ sin(π(x − y))
lim f ( nλi , nλj ) = f (x, y) 1 − dxdy. (7.77)
n→∞
ij −∞ −∞ π(x − y)

Note that the quadratic repulsion manifest in Theorem 7 is clear here in (7.76), which, if we set
x − y = w, vanishes like w2 as w → 0. We also remark that for the Poisson point process, the
corresponding formula is obtained by replacing the two-point correlation function (7.76) by 1, because
the points are, by definition, uncorrelated.
We have shown that at the centre of the spectrum (i.e. at the centre of the semicircle) the limiting
correlations are determined by the sine kernel. In fact, the sine kernel can be shown to describe the
local correlations at any point in (−2, 2) provided one normalizes the mean spacing to be unity using
the semicircle density.

7.5 Scaling at the edge of the spectrum

Having seen how the sine kernel describes correlations in the bulk of the spectrum, we now consider
how the eigenvalues are distributed near to the edges. By the semicircle law, we expect to have, when 
is small, that Z 2 p
√ n 2 3/2
#{λi : λi / n ≥ 2 − } ≈ 4 − x2 dx ≈ n . (7.78)
2π 2− 3π
We might therefore expect to get a nice limit if we set  ∝ n−2/3 , i.e. if we rescale the distance
between the eigenvalues and the edge by n−2/3 . If we want to find a limiting √ distribution for the
eigenvalues close to the upper edge of the spectrum, this suggests setting λi / n = 2 + αi n−2/3 , or
27 One can also use the method of steepest descent; there are two saddle points and the contributions from each have

to be summed.

48
 √
λi = 2 n + αi n−1/6 and looking at correlations between the values taken by α1 , . . . , αn . This involves
studying the asymptotics of
h √ √ ir
det n−1/6 Kn (2 n + αi n−1/6 , 2 n + αj n−1/6 ) . (7.79)
ij=1

Defining √ √
K̃n (x, y) = n−1/6 Kn (2 n + xn−1/6 , 2 n + yn−1/6 ) (7.80)
we have the following theorem describing the n → ∞ asymptotics.
Theorem 15. For any fixed x, y

Ai(x)Ai0 (y) − Ai0 (x)Ai(y)

lim K̃n (x, y) = (7.81)
n→∞ (x − y)

where Z ∞
1
Ai(x) = cos(t3 /3 + xt)dt (7.82)
π 0
is the Airy function.
This is an immediate consequence of the following lemma
Lemma 16. For any fixed x
√
lim n1/12 φn (2 n + x/n1/6 ) − Ai(x) = 0. (7.83)
n→∞

We only sketch the proof of this lemma. This starts with (7.74). One again seeks to evaluate the
integral by the method of steepest descent. The difference is that in this case the two saddle points lie
so close together that one cannot treat them separately. Instead of expanding to quadratic order, one
has to expand to cubic order. This leads to the Airy function, rather than the usual Gaussian integral.
The Airy function is plotted in Figure 8.
Substituting the asymptotic estimate in the lemma into the formula for K̃n in terms of φ, and using
(7.36) and (7.37) to evaluate the asymptotics of the derivative of φ, then leads to the formula in the
theorem.
The kernel representing the n → ∞ limit of K̃n is known as the Airy kernel. Far away from the
edge on the inside of the semicircle, it reduces to the sine kernel. Away from the edge on the outside of
the semicircle, it decays exponentially.

7.6 Counting statistics

Having understood the structure of spectral correlation functions in the bulk of the spectrum and at
the edge, we now turn our attention to counting statistics. Let us start by looking at an interval I in
the bulk of the spectrum. Denote by ν(M, I) the number of eigenvalues of a GUE matrix M lying in I.
We shall focus on
Dm (I) = P(ν(M, I) = m). (7.84)
Let II (x) denote the characteristic function on I, so
(
1 if x ∈ I
II (x) = (7.85)
0 if x ∈
/ I.

Then, for example,

n
Y
D0 (I) = E (1 − II (λi )). (7.86)
i=1

49
 0.4

0.2

-15 -10 -5 5 10

-0.2

-0.4

Figure 8: The Airy function Ai(x).

We define
n
Y
GI (t) = E (1 − tII (λi )), (7.87)
i=1

so that D0 (I) = GI (1). It is easy to see that −G0I (1) = D1 (I), and that more generally
(−1)m (m)
Dm (I) = GI (1). (7.88)
m!
GI (t) is therefore the generating function for Dm (I).
To calculate GI (t), we expand the product over the eigenvalues to get28
X X
GI (t) = 1 − tE II (λi ) + t2 E II (λi )II (λj ) − . . .
i i≤j
Z Z Z2 Z Z Z
t t3
= 1 − t R1 (x1 )dx1 + R2 (x1 , x2 )dx1 dx2 − R3 (x1 , x2 , x3 )dx1 dx2 dx3 + . . .
I 2 I I 3! I I I
X∞ Z Z
(−t)k
=1+ · · · det (Kn (xi , xj ))1≤i,j≤k dx1 . . . dxk . (7.89)
k! I I
k=1

The last equation can be viewed as defining the Fredholm determinant of the operator with kernel
Kn (x, y) acting on L2 (I):
∞
X Z Z
(−t)k
det (I − tKn (x, y)) := 1 + · · · det (Kn (xi , xj ))1≤i,j≤k dx1 . . . dxk . (7.90)
k! I I
k=1
28 Note the the expansion is valid because the rank of Kn is at most n.

50
 In the theory of Fredholm
R determinants29 an operator Rwith kernel K(x, y) acts on a function f
according to (Kf )(x) = K(x, y)f (y)dy. When Tr(K) = K(x, x)dx exists, then K is called trace
class and one can calculate the Fredholm determinant in a number of (equivalent) ways:
∞
!
X tm m
det (I − tKn (x, y)) = exp − TrK (7.91)
m=1
m

where K m is to be understood as an m-fold convolution of K. Alternatively, denoting the eigenvalues

of K by βj , Y
det (I − tKn (x, y)) = (1 − tβj ). (7.92)
R P
We also have that then K(x, x)dx = j βj
The probability Pn (s)ds that the spacing between adjacent (i.e. nearest neighbouring) eigenvalues
lies between s and s+ds can be calculated in terms of the Fredholm determinant as follows. Consider the
interval I to have length s. The probability that this interval contains no eigenvalues is F (s) = D0 (I).
Now extend the interval by ds at one end. The probability that this extended interval has no eigenvalues
in it is F (s + ds). Hence F (s) − F (s + ds) is the probability that there is no eigenvalue in the original
interval, but at least one in the extended interval. Now extending the interval at the other end, one
sees that
d2 F d2 GI(s) (t)
Pn (s) = 2
= . (7.93)
ds ds2 t=1

Hence
√ one can calculate Pn (s) from the Fredholm determinant. When the eigenvalues are scaled by
n, then when n → ∞ the limiting spacing distribution is the second derivative of the Fredholm
determinant of the sine kernel30 . This distribution is often called the Gaudin-Mehta distribution
We illustrate this by showing in Figure 9 the results of a numerical experiment. This involved
generating 10,000 random matrices from the GUE with n = 10, plotting a histogram of the normalized
distances between neighbouring pairs of eigenvalues, and comparing with the Fredholm determinant of
the sine kernel.
The analysis of Dm (I) extends immediately to the edge of the spectrum, with, in this case, the
kernel in (7.89) and (7.90) being the Airy kernel. Rather than the spacing distribution, the natural
question in this case is: what is the distribution of the largest eigenvalue of a GUE matrix? Let λmax
denote the largest eigenvalue and set

Fn(max) (t) := P(λmax ≤ t). (7.94)

Then from the scaling at the edge we have that

√
F (max) (s) = lim Fn(max) (2 n + s/n1/6 ) (7.95)
n→∞

exists. Clearly P(λmax < t) is simply the probability that the interval (t, ∞) contains no eigenvalues,
and so
F (max) (s) = det(I − KAiry ) (7.96)
where
Ai(x)Ai0 (y) − Ai0 (x)Ai(y)
KAiry = (7.97)
(x − y)
acting on L2 ((s, ∞)).
29 We limit ourselves to stating results in the theory of Fredholm determinants without proof, and without setting out
the theory systematically. For an introductory exposition see, forexample, [12].
30 The Fredholm determinant of the sine kernel is equal to exp
R πs σ(x;t)
0 x
dx where σ is the solution of the Painlevé
V equation (xσ 00 )2 + 4(xσ 0 − σ)(xσ 0 − σ + (σ 0 )2 ) = 0 with σ(x; t) ∼ − πt x as x → 0, but showing this is beyond the scope
of the present course.

51
 0.75
density

0.50

0.25

0.00

0 1 2 3
spacing

Figure 9: The spacing distribution between neighbouring eigenvalues of 10,000 GUE matrices of
dimension n = 10 compared to the the Fredholm determinant of the sine kernel. (Plot kindly provided
by Johannes Forkel.)

52
0.4

0.3

0.2

0.1

0.0
-5 -4 -3 -2 -1 0 1

√
Figure 10: (λmax − 2 n)n1/6 for 10,000 GUE matrices of size n = 1000.

In this case too F (max) (s) can be expressed in terms of a solution of a second order nonlinear
o.d.e. which is another of the Painlevé equations. Specifically,
 Z ∞ 
(max) 2
F (s) = exp − (x − s)q(x) dx (7.98)
s

where q is a solution of the Painlevé II equation q 00 (x) − xq(x) + 2q(x)3 = 0 with q(x) ∼ Ai(x) as
x → ∞. This is called the Tracy-Widom distribution.
We illustrate this by showing in Figure 10 the results of a numerical experiment, which involved
generating
√ 10,000 random matrices from the GUE with n = 1000, plotting a histogram of (λmax −
2 n)n1/6 , and comparing with the Tracy-Widom distribution.

7.7 Other Ensembles

The method of orthogonal polynomials extends straightforwardly to the CUE. Indeed, in this case it is
considerably simpler, because the functions that are orthogonal with respect to the uniform weight
on the unit circle are simply the Fourier functions eipθ – all that is involved is the manipulation of
Fourier series and the asymptotic analysis is elementary. So, for the Vandermonde factor in the Weyl
integration formula
Y Y
Y −iθj

|eiθj − eiθk |2 = eiθj − eiθk e − e−iθk
j<k j<k j<k
ipθj −ipθj
= det[e ] det[e ]. (7.99)

53
In these determinants j = 1, 2, . . . , n and p = 0, 1, . . . , n − 1, or, equivalently, p = − 12 (n − 1), − 12 (n −
1) + 1, . . . , 12 (n − 1). Hence
Y h i
|eiθj − eiθk |2 = (2π)n det S̃n (θj − θk ) (7.100)
j,k=1,...,n
j<k

with
1
2 (n−1)
X
1 1 sin(nθ/2)
S̃n (θ) = eipθ = . (7.101)
2π 2π sin(θ/2)
p=− 12 (n−1)

Using its representation as a Fourier series

Z 2π
S̃n (θj − θk )S̃n (θk − θl )dθk = S̃n (θj − θl ). (7.102)
0

We can therefore apply Gaudin’s lemma, but now with integrals over [0, 2π). Trivially,
Z 2π
S̃n (0)dθ = n (7.103)
0

and so Gaudin’s lemma yields, exactly as for the GUE, that the k-point correlation function is
h i
Rk (θ1 , . . . , θk ) = det S̃n (θi − θj ) . (7.104)
i,j=1,...,k

This is again a k × k determinant, but now with a simpler kernel. It is easy to see that
n
R1 (θ) = S̃n (0) = (7.105)
2π
which is the analogue of the semicircle law in this setting – the eigenvalues have a constant average
density around the unit circle – and that
 n 2  2
R2 (θ1 , θ2 ) = − S̃n (θ1 − θ2 ) . (7.106)
2π
Now crucially, observe that if we scale the eigenvalues to have unit mean spacing, this corresponds
n
to setting φk = 2π θk . Making this substitution in S̃n and rescaling by 2π/n, we see that

1 sin(π(φj − φk ) sin(π(φj − φk )
lim S̃n (2π(φj − φk )/n) = lim = (7.107)
n→∞ n n→∞ n sin(π(φj − φk )/n) π(φj − φk ))

i.e. the limit gives precisely the sine kernel we found for the GUE. Hence, even though the statistics are
different when n is finite, in the limit when n → ∞ the local statistics of the GUE and CUE are the
same.
We illustrate this by showing in Figure 11 the results of a numerical experiment. This involved
generating 200 random matrices from the CUE with n = 200, plotting a histogram of the distances
between all pairs of eigenvalues, and comparing with the 2 × 2 determinant involving the sine kernel.
The method of orthogonal polynomials works as well for the GOE, but is a little more complicated in
that case. It works for complex Wishart matrices, when the orthogonal polynomials are the generalized
Laguerre polynomials. In the case of real Wishart matrices, as for the GOE, the method also applies,
but it is a little more complicated.

54
 1.00

0.75
density

0.50

0.25

0.00

0 1 2 3 4
distances

Figure 11: The two-point correlation function of the eigenvalues of 200 CUE matrices of dimension
200 compared to the 2 × 2 determinant involving the sine kernel. (Plot kindly provided by Johannes
Forkel.)

55
7.8 Universality
Given that the GUE and the CUE exhibit the same local statistics on the scale of the mean eigenvalue
spacing in the limit n → ∞ one might ask whether other random matrix ensembles also fall into the
same class. Does the result depend on having simple measures, for example? Or does it rely on having
nice formulae for the measure in terms of the eigenvalues, expressed via Vandermonde determinants?
The answer is that in the large-matrix limit the results are independent of the measure, provided we
have matrices that are essentially complex and Hermitian (unitary matrices can be thought of as the
exponentials of complex Hermitian matrices). So complex Hermitian Wigner random matrices all have
the same limit, once some mild conditions are satisfied, no matter what the measure one applies to the
matrix elements, and for most measures one doesn’t have nice expressions in terms of the eigenvalues.
Similarly, real-symmetric matrices form a different universality class, characterised by linear repulsion
between the eigenvalues as opposed to quadratic repulsion. The GOE falls into this class, along with
other real-symmetric Wigner matrices, irrespective of the measure.
Proving universality has been a major theme of research in Random Matrix Theory over the past
15 years.

7.9 Applications
The local spectral statistics of random matrices in the bulk of the spectrum have many applications.
For example, they describe fluctuation statistics in complex quantum systems, including quantum
chaotic systems, atomic spectra and nuclear spectra, and in other complex wave problems, such as
lasers, elastic vibrations, and acoustics etc. They also provide the main method for modelling statistical
properties of quantum dynamics in complex systems, such as in conductivity through disordered media.
Remarkably, GUE spectral statistics also appear to describe correlations between the positions of
parked cars in London, and arrival times of buses in the Mexican city of Cuernavaca.
Knowing the distribution of largest and smallest eigenvalues is important in many contexts in
numerical linear algebra in determining bounds on the efficiency and convergence of algorithms. In
Wishart random matrices the ratio of the extremal eigenvalues determines how errors in the input data
get magnified in solving systems of linear equations, as measured by the condition numbers of the
matrices appearing.
For more on these applications, see the respective chapters in The Oxford Handbook of Random
Matrix Theory [1].

56
8 Dyson Brownian motion
We can express the GOE and GUE measures in terms of the eigenvalues, as embodied in Theorem 7, in
the form
−βW (λ1 ,...,λn )
c(β)
n e (8.1)
where
n
1X 2 X
W (λ1 , . . . , λn ) = λ − ln |λi − λj |. (8.2)
4 j=1 j i<j

One can think of this as a Boltzmann weight, with W (λ1 , . . . , λn ) representing the potential energy
associated with particles located at positions λ1 , . . . , λn and with β playing the role of the inverse
temperature. Hence computing expectations with respect to the RMT measures in this case is the same
as determining the equilibrium thermodynamics of particles moving in one dimension with interactions
described by this energy function.
We now explore the dynamics of the eigenvalues if we allow the matrix elements of a real-symmetric
or complex-Hermitian matrix to vary stochastically with time t. Specifically, let us consider the matrix
M to be a function of t satisfying, for t ≥ 0, the matrix-valued stochastic differential equation
1 1
dM (t) = √ dB(t) − M (t)dt (8.3)
n 2

with initial data M0 . Here in the real-symmetric√case, when β = 1, B(t) is an n × n symmetric matrix
such that Bij (t), for 1 ≤ i < j ≤ n, and Bii (t)/ 2, for 1 ≤ i ≤ n, are independent standard Brownian
motions; and
√ in the complex-Hermitian
√ case, when β = 2, B(t) is an√n × n complex Hermitian matrix
such that 2ReBij (t) and 2ImBij (t), for 1 ≤ i < j ≤ n, and Bii (t)/ 2, for 1 ≤ i ≤ n, are independent
standard Brownian motions. This defines what is known as a matrix-valued Ornstein-Uhlenbeck process.
In terms of the matrix elements of M , this differential equation reads
1 1
dMij (t) = √ dBij (t) − Mij (t)dt (8.4)
n 2

where Bij (t) has variance t in the complex-Hermitian case, while in the real-symmetric case Bij (t) has
variance t when i 6= j and variance 2t when i = j.
The question we seek to address is: if the matrix M evolves in time according to (8.3), what is
the equation satisfied by the eigenvalues λ1 (t), λ2 (t), . . . , λn (t)? First, one might wonder about the
labelling of the eigenvalues if they become degenerate or exchange positions. In fact, if we assume
that the eigenvalues at t = 0 are non-degenerate and if we label them in order of increasing size,
λ1 (0) < λ2 < (0) < · · · < λn (0), then as time evolves it can be shown that the eigenvalues remain
simple and are continuous functions of t, and so the labelling makes sense and is preserved for all t > 0.
The equation governing the time evolution of the eigenvalues might be expected to involve the
eigenvalues and eigenvectors vj (t) ∈ Rn . However, remarkably, this turns out not to be the case: the
eigenvalues satisfy an autonomous system of stochastic differential equations that do not involve the
eigenvectors. This system of SDEs is the following one.
 
√
2 λ i 1 X 1
dλi = √ dB̃i + − +  dt (8.5)
βn 2 n λi − λj
j6=i

where B̃i is a set of real-valued independent standard Brownian motions. The solution of this equation
is known as Dyson Brownian motion with parameter β, after Freeman Dyson, who introduced it in a
seminal paper in 1962.
As noted above, it can be proved (see, for example, section 4.3 in [2]) that there exists a unique
(strong) solution of this equation in the space of continuous functions. Moreover, if the initial conditions

57
are such that all of the eigenvalues are simple at t = 0, then this remains true for all t > 0 – i.e. none
of the eigenvalue processes intersect.
Dyson’s theorem can be stated as follows.
Theorem 17. Let M (t) solve (8.3) in a strong sense. Then its eigenvalue process satisfies (8.5).
Proof.. We denote the eigenvalues of an n × n real-symmetric or complex-Hermitian matrix M by
(λm )1≤m≤n and the corresponding eigenvectors by (vm )1≤m≤n , so that
M vm = λm vm . (8.6)
Note that the eigenvectors are orthonormal, i.e.
†
vm vp = δm,p . (8.7)
Differentiating these equations with respect to Mij yields
∂M ∂vm ∂λm ∂vm
vm + M = vm + λm , (8.8)
∂Mij ∂Mij ∂Mij ∂Mij
†
∂vm † ∂vp
vp + vm =0 (8.9)
∂Mij ∂Mij
for p 6= m, and
†
∂vm
vm = 0. (8.10)
∂Mij
Hence
∂λm † ∂M
= vm vm (8.11)
∂Mij ∂Mij
and for p 6= m
∂M ∂vm ∂vm
vp† vm + vp† M = λm vp† , (8.12)
∂Mij ∂Mij ∂Mij
from which we have that
∂M ∂vm ∂vm
vp† vm + λp vp† = λm vp† . (8.13)
∂Mij ∂Mij ∂Mij
Hence
∂vm X  ∂vm  X 1

∂M

= vp† vp = vp† vm vp . (8.14)
∂Mij ∂Mij λm − λp ∂Mij
p6=m p6=m

In order to keep the notation as simple as possible, we now focus on the real-symmetric case, when
β = 1. The complex-Hermitian calculation follows exactly the same lines. In the real-symmetric case,
(i)
denoting the ith component of vm by vm , (8.11) becomes
∂λm (i) (j)
= vm vm (2 − δij ), (8.15)
∂Mij
and (8.14) becomes
(k)
∂vm X vp(i) vm
(j) (j) (i)
+ vp vm (1 − δij ) (k)
= vp . (8.16)
∂Mij λm − λp
p6=m

Differentiating the first of these equations and then using the second gives
" #
(i) (k)
∂ 2 λm ∂vm (k) (i) ∂vm
= (2 − δik ) v + vm
∂Mlj ∂Mik ∂Mlj m ∂Mlj
X 1 h i
= (2 − δik ) (vp(j) vm(l)
+ vp(l) vm
(j)
(1 − δjl ))(vp(i) vm
(k)
+ vp(k) vm
(i)
) . (8.17)
λm − λp
p6=m

58
 We are now in a position to use the assumption that the matrix elements evolve as an Ornstein-
Uhlenbeck process (8.4), which means that in a time interval of length dt their square scales like dt.
We therefore expand dλm to quadratic order in the dMij 31 . We then get32
X ∂λm 1 X X ∂ 2 λm
dλm = dMik + (dMik )(dMlj ). (8.18)
∂Mik 2 ∂Mlj ∂Mik
i≤k i≤k j≤l

Now, using (8.4), (8.15), (8.17), and the fact that

1
(dMik )(dMlj ) = δil δkj (1 + δik )dt (8.19)
n
we have, using the fact that B is symmetric, that
X  
(i) (k) 1 Mik
dλm = vm vm √ dBik − dt
n 2
i,k
1 XX 1 h i
+ |vp(i) |2 |vm
(k) 2 (i) 2 (k) 2
| + |vm | |vp | dt. (8.20)
2n λm − λp
i,k p6=m

Hence, using the fact that vm is the eigenvector of M with eigenvalue λm , and so
X X X
(i) (k) (i)
vm vm Mik = vm (M vm )(i) = λm (i) 2
|vm | = λm , (8.21)
i,k i i

we have finally that

1 X (i) (k) 1 1 X 1
dλm = √ vm vm dBik − λm dt + dt. (8.22)
n 2 n λm − λp
i,k p6=m

The key point now is that, looking at the first term in (8.22), we can define
X
(i) (k)
B̃m := vm vm Bik . (8.23)
i,k

This is itself a real Gaussian process satisfying EdB̃m = 0 and

XX 0 0
(i) (k)
EdB̃m dB̃p = E vm vm dBik vp(i ) vp(k ) dBi0 k0
i,k i0 ,k0
X
(i) (i) (k) (k)
=2 vm vp vm vp dt
i,k

= 2δmp dt (8.24)
√
We see therefore that B̃m = 2Bm in distributional sense, where (Bm )nm=1 is the standard Brownian
motion in Rn . This completes the proof in the real-symmetric case. A similar calculation along the
same lines gives the corresponding result when β = 2.

It is clear from (8.3) that the Ornstein-Uhlenbeck process starting at t = 0 from

√ some fixed matrix
M0 is equal in distribution at any later fixed time t to the sum of e−t/2 M0 , and 1 − e−t G, where G is
31 This corresponds to applying Ito’s formula from stochastic calculus.
32 There is a technical subtlety that we are ignoring in this step: one needs to prove that the singularities coming from
the factors (λi − λj )−1 do not overwhelm the other terms. To do this one needs to show that under the dynamics the
eigenvalues do not get too close. That this is true is a consequence of the eigenvalue repulsion in the GOE and GUE. We
do not give the argument here; see, for example, Section 4.3.1 of [2] for the proof that this subtlety can indeed be ignored.

59
Figure 12: A simulation of GOE Dyson Brownian motion when n = 10. (Plot kindly provided by
Johannes Forkel.)

a GOE matrix when β = 1 or a GUE matrix when β = 2. Hence for sufficiently large t the GOE/GUE
is the equilibrium solution which the dynamics reaches, no matter what M0 is taken to be. Put another
way, the GOE/GUE measure is invariant under Dyson Brownian motion, and is the attractor for the
dynamics. It was conjectured by Dyson that this equilibrium is actually reached very quickly, on time
scales of the order of 1/n, and this has subsequently been proved. This fact plays a key role in the
proof of the universality of the spectral statistics.
We illustrate this by showing in Figures 12, 13 and 14 the results of a numerical simulation of the
Ornstein-Uhlenbeck process with, respectively, n = 10, n = 20 and n = 50 for GOE matrices. The
invariance of the GOE under this process is illustrated in Figures 15, 16 and 17, where M0 is taken to
be a GOE matrix in each case, again with, respectively, n = 10, n = 20 and n = 50.

60
Figure 13: A simulation of GOE Dyson Brownian motion when n = 20. (Plot kindly provided by
Johannes Forkel.)

Figure 14: A simulation of GOE Dyson Brownian motion when n = 50. (Plot kindly provided by
Johannes Forkel.)

61
Figure 15: A simulation of GOE Dyson Brownian motion when n = 10, starting from a GOE matrix.
(Plot kindly provided by Johannes Forkel.)

Figure 16: A simulation of GOE Dyson Brownian motion when n = 20, starting from a GOE matrix.
(Plot kindly provided by Johannes Forkel.)

62
Figure 17: A simulation of GOE Dyson Brownian motion when n = 50, starting from a GOE matrix.
(Plot kindly provided by Johannes Forkel.)

9 Some connections with other areas of mathematics

We have already mentioned some of the applications of Random Matrix Theory, for example to data
science, quantum mechanics, mathematical finance, population dynamics etc. There are also many
beautiful and deep connections with other problems in mathematics, and we mention a few examples
here.

9.1 Longest increasing subsequences

As discussed in the introduction, let Sn be the group of permutations of 1, 2, . . . , n. If π ∈ Sn ,
π(i1 ), . . . , π(ik ) is an increasing subsequence in π if i1 < i2 < · · · < ik and π(i1 ) < π(i2 ) < · · · < π(ik ).
Let ln (π) be the length of the longest increasing subsequence. For example, if n = 5 and π is the
permutation 5 1 3 2 4, then the longest increasing subsequences are 1 2 4 and 1 3 4, and ln (π) = 3.
Equip Sn with uniform distribution,

#{π ∈ Sn : ln (π) ≤ m}
P(ln ≤ m) = . (9.1)
n!
Then it was proved by Baik, Deift and Johansson in 1998 that
 √ 
ln − 2 n
lim P ≤ s = F (max) (s) = det(I − KAiry ) (9.2)
n→∞ n1/6

where
Ai(x)Ai0 (y) − Ai0 (x)Ai(y)
KAiry = (9.3)
(x − y)

63
acting on L2 ((s, ∞)). Or equivalently,
 √   Z ∞ 
ln − 2 n 2
lim P ≤ s = exp − (x − s)q(x) dx (9.4)
n→∞ n1/6 s

where q is a solution of the Painlevé II equation q 00 (x) − xq(x) + 2q(x)3 = 0 with q(x) ∼ Ai(x) as
x → ∞.
That is, ln has the same limiting distribution as the largest eigenvalue of a random GUE matrix.

9.2 Stochastic growth models, random tilings, and random paths

There has been a considerable focus in the past few years on the geometrical properties of structures
grown or generated by simple random processes. One example involves random tilings of Aztec
diamonds, and a second involves randomly grown surfaces, such as by a random deposition of atoms.
These are generated by simple probabilistic rules. We do not give the details, but show an example of
a random tiling in Figure 18. The outside can be thought of as a ‘frozen’ phase, and the interior a
‘liquid’ phase. The boundary of the liquid phase is a random curve whose fluctuations are again the
same as those found in the largest eigenvalue of a GUE matrix. This models laboratory experiments
involving growing interfaces of liquid-crystal turbulence, where the fluctuations of the interface are
found to match those of the largest eigenvalue of a GUE matrix rather well. For more details (and a
movie) see [18]. For more about random growth models and random tilings, see the article by P. Ferrari
and H. Spohn in [1].
One finds the same again in models of random paths. For example, the Hammersley process refers
to the following stochastic model. In the unit square mark in points uniformly at random according to
a Poisson point process with intensity α. We call a path from (0, 0) to (1, 1) through these random
points up/right if the points it passes through have coordinates xk ≤ xk+1 and yk ≤ yk+1 for each k.
Let L(α) denote the maximum number of points on such a path. Then for all s ∈ R,
 √ 
L(α) − 2 α
lim P ≤ s = F (max) (s) = det(I − KAiry ). (9.5)
α→∞ α1/6

9.3 Zeros of the Riemann zeta-function

The Riemann zeta-function, ζ(s) is defined when Res > 1 by

X∞
1 Y 1
−1
ζ(s) = = 1 − (9.6)
n=1
ns p
ps

where the product runs over all primes p. It has an analytic continuation to the rest of the complex
plane, except for a pole at s = 1. It is important because it encodes information about the distribution
of the primes. The zeta function has trivial zeros at s = −2n for n ∈ N, and infinitely many other
nontrivial zeros. The Riemann Hypothesis asserts that the nontrivial zeros all lie on the line Res = 1/2;
that is, they are all of the form 1/2 + itn with tn ∈ R.
Let us assume that the Riemann Hypothesis is true and so the numbers tn are all real. The theory
of the zeta function then implies that
 
T T
#{n : 0 < tn ≤ T } = log + O(log T ) (9.7)
2π 2π

Therefore, setting  
tn |tn |
wn = log (9.8)
2π 2π

64
Figure 18: A random tiling, generated by P. Ferrari.

65
 Z 1
fˆ(ø ) = f (x)e2ºixø dx (2.5)
°1

has support in (°1, 1). Numerical evidence suggests that Montgomery’s conjecture is correct
[58], see Figure 1.

R[x]
1

0.8

0.6

0.4

0.2

x
0.5 1 1.5 2 2.5 3

Figure 1: Two-point correlation function of the scaled Riemann zeros plotted using 106
zeros Figure
around 19: the 1020 thcorrelation
The 2-point
20
zero, computed
function for by A. Odlyzko,
the Riemann and compared
zeros, computed withforRzeros
by A. Odlyzko 2 (x) from
near to the 10 th, compared to the GUE 2-point correlation function (9.11).
(2.4).

we have that the mean density of the numbers wn is asymptotically 1, in that

Shortly after Montgomery completed 1 this work he was introduced to Freeman Dyson.
lim #{n : 0 < wn ≤ W } = 1 (9.9)
After hearing about Montgomery’s W →∞latest
W results Dyson [34] recognized in R2 (x) the two-
point correlation function
Consider now of eigenvalues
the pair correlation function of random
of the unitary
scaled zeros matrices
wn , defined, defined,
assuming forexists,
the limit a suitable
test function
by f (x, y), as: X Z ∞
1
lim f (wn − wm ) = f (x)R2 (x)dx (9.10)
N →∞ N −∞
1≤n<m≤N
2
It is a theorem due to Montgomery in 1973 that, assuming the Riemann Hypothesis, for functions f (x)
whose Fourier transform has support in (−1, 1) the limit in the left hand-side exists and then
 2
sin πx
R2 (x) = 1 − (9.11)
πx

which is precisely the GUE/CUE 2-point correlation function. Montgomery conjectured that this
remains true for all functions f for which the sums converge. This has since been extended to all
k-point correlation functions, where the result and conjecture lead to the k × k determinant of the
sine kernel. It is therefore currently the belief that that all local statistics of the zeros coincide with
those of the GUE/CUE, and this is supported by extensive numerical computations – see for example
Figure 19. Assuming this is true, many interesting properties of the zeta function and the primes can
be calculated using random matrix theory. See the chapter on connections between number theory and
random matrix theory in [1]. This remains a highly active area of research.

References
[1] G. Akemann, J. Baik & P. Di Francesco, The Oxford Handbook of Random Matrix Theory (Oxford
University Press).

66
[2] G.W. Anderson, A. Guionnet & O. Zeitouni, An Introduction to Random Matrices (Cambridge
Studies in Advanced Mathematics).
[3] J.P. Bouchaud & M. Potters, Financial Applications of Random Matrix Theory: a short review,
in The Oxford Handbook of Random Matrix Theory, edited by Gernot Akemann, Jinho Baik, &
Philippe Di Francesco (OUP), ;also available from arXiv:0910.1205
[4] Z. Bai & J.W. Silverstein, Spectral Analysis of Large Dimensional Random Matrices (Springer).
[5] A. Choromanska, M. Henaff, M. Mathieu, G. Ben Arous & Y. LeCun, The loss surfaces of multilayer
networks, Artificial intelligence and statistics, 192–204 (2015).
[6] P. Diaconis & P.J. Forrester, Hurwitz and the origins of random matrix theory in mathematics,
Random Matrices: Theory and Applications, 06, 1730001 (2017); arXiv:1512.09229.
[7] A. Edelman & N. Raj Rao, Random matrix theory, Acta Numerica 14, 233–297 (2005).

[8] I.M. Johnstone, High dimensional statistical inference and random matrices, Proceedings of the
International Congress of Mathematicians, Madrid, Spain, 2006; arXiv:math/0611589.
[9] L. Laloux, P. Cizeau, M. Potters & J.P. Bouchaud, Random matrix theory and financial correlations,
available from [Link]
[10] T. Koshy, Catalan Numbers with Applications, Oxford University Press (2008).

[11] R.M. May, Will a large complex system be stable?, Nature 238, 413-414 (1972).
[12] H.P. McKean, Fredholm Determinants, Central European Journal of Mathematics 9, 205–243
(2011).
[13] E.S. Meckes, The Random Matrix Theory of the Classical Compact Groups,
Cambridge Tracts in Mathematics 218 (CUP, 2019); available online from
[Link]
[14] F. Mezzadri, How to generate random matrices from the classical compact groups, Notices of the
American Mathematical Society 54, 592–604 (2007); arXiv: math-ph/0609050.

[15] J. von Neumann & H. H. Goldstine, Numerical inverting of matrices of high order, Bull. Amer.
Math. Soc. 53, 1021–1099 (1947).
[16] G. Szegö, Orthogonal Polynomials, 4th ed. Providence, RI: Amer. Math. Soc., 42–44, (1975).
[17] R.P. Stanley, Catalan numbers, Cambridge University Press (2015).

[18] K.A. Takeuchi, M. Sano, T. Sasamoto, & H. Spohn, H., Growing interfaces uncover universal
fluctuations behind scale invariance, Scientific Reports, 1: 34 (2011); arXiv:1108.2118
[19] E.P. Wigner, On the distribution of the roots of certain symmetric matrices, Annals of Mathematics
67, 325–327 (1958).
[20] J. Wishart, The Generalised Product Moment Distribution in Samples from a Normal Multivariate
Population, Biometrika 20A, 32–52 (1928).

67