Introduction to Random Matrices - Theory and Practice

Preface

This is a book for absolute beginners. If you have heard about random matrix theory, commonly denoted RMT, but you do not know what that is, then welcome!, this is the place for you. Our aim is to provide a truly accessible introductory account of RMT for physicists and mathematicians at the beginning of their research career. We tried to write the sort of text we would have loved to read when we were beginning Ph.D. students ourselves.

Our book is structured with light and short chapters, and the style is informal. The calculations we found most instructive are spelt out in full. Particular attention is paid to the numerical verification of most analytical results. The reader will find the symbol [$\spadesuit$ test.m] next to every calculation/procedure for which a numerical verification is provided in the associated file test.m located at . We strongly believe that theory without practice is of very little use: in this respect, our book differs from most available textbooks on this subject (not so many, after all).

Almost every chapter contains question boxes, where we try to anticipate and minimize possible points of confusion. Also, we include To know more sections at the end of most chapters, where we collect curiosities, material for extra readings and little gems - carefully (and arbitrarily!) cherrypicked from the gigantic literature on RMT out there.

Our book covers standard material - classical ensembles, orthogonal polynomial techniques, spectral densities and spacings - but also more advanced and modern topics - replica approach and free probability - that are not normally included in elementary accounts on RMT.

Due to space limitations, we have deliberately left out ensembles with complex eigenvalues, and many other interesting topics. Our book is not encyclopedic, nor is it meant as a surrogate or a summary of other excellent existing books. What we are sure about is that any seriously interested reader, who is willing to dedicate some of their time to read and understand this book till the end, will next be able to read and understand any other source (articles, books, reviews, tutorials) on RMT, without feeling overwhelmed or put off by incomprehensible jargon and endless series of “It can be trivially shown that….”.

So, what is a random matrix? Well, it is just a matrix whose elements are random variables. No big deal. So why all the fuss about it? Because they are extremely useful! Just think in how many ways random variables are useful: if someone throws a thousand (fair) coins, you can make a rather confident prediction that the number of tails will not be too far from $500$. Ok, maybe this is not really that useful, but it shows that sometimes it is far more efficient to forego detailed analysis of individual situations and turn to statistical descriptions.

This is what statistical mechanics does, after all: it abandons the deterministic (predictive) laws of mechanics, and replaces them with a probability distribution on the space of possible microscopic states of your systems, from which detailed statistical predictions at large scales can be made.

This is what RMT is about, but instead of replacing deterministic numbers with random numbers, it replaces deterministic matrices with random matrices. Any time you need a matrix which is too complicated to study, you can try replacing it with a random matrix and calculate averages (and other statistical properties).

A number of possible applications come immediately to mind. For example, the Hamiltonian of a quantum system, such as a heavy nucleus, is a (complicated) matrix. This was indeed one of the first applications of RMT, developed by Wigner. Rotations are matrices; the metric of a manifold is a matrix; the $S$-matrix describing the scattering of waves is a matrix; financial data can be arranged in matrices; matrices are everywhere. In fact, there are many other applications, some rather surprising, which do not come immediately to mind but which have proved very fruitful.

We do not provide a detailed historical account of how RMT developed, nor do we dwell too much on specific applications. The emphasis is on concepts, computations, tricks of the trade: all you needed to know (but were afraid to ask) to start a hopefully long and satisfactory career as a researcher in this field.

It is a pleasure to thank here all the people who have somehow contributed to our knowledge of RMT. We would like to mention in particular Gernot Akemann, Giulio Biroli, Eugene Bogomolny, Zdzisław Burda, Giovanni Cicuta, Fabio D. Cunden, Paolo Facchi, Davide Facoetti, Giuseppe Florio, Yan V. Fyodorov, Olivier Giraud, Claude Godreche, Eytan Katzav, Jon Keating, Reimer Kühn, Satya N. Majumdar, Anna Maltsev, Ricardo Marino, Francesco Mezzadri, Maciej Nowak, Yasser Roudi, Dmitry Savin, Antonello Scardicchio, Gregory Schehr, Nick Simm, Peter Sollich, Christophe Texier, Pierfrancesco Urbani, Dario Villamaina, and many others.

This book is dedicated to the fond memory of Oriol Bohigas.

The final publication is available at Springer via .

Chapter 1 Getting Started

Let us start with a quick warm-up. We now produce a $N\times N$ matrix $H$ whose entries are independently sampled from a Gaussian probability density function (pdf)You may already want to give up on this book. Alternatively, you can brush up your knowledge about random variables in Section 1.1. with mean and variance $1$. One such matrix for $N=6$ might look like this:

Some of the entries are positive, some are negative, none is very far from . There is no symmetry in the matrix at this stage, $H_{ij}\neq H_{ji}$.

Any time we try, we end up with a different matrix: we call all these matrices samples or instances of our ensemble. The $N$ eigenvalues are in general complex numbers (try to compute them for $H$!).

To get real eigenvalues, the first thing to do is to symmetrize our matrix. Recall that a real symmetric matrix has $N$ real eigenvalues. We will not deal much with ensembles with complex eigenvalues in this book…but we will deal a lot with matrices with complex entries (and real eigenvalues)..

Try the following symmetrization $H_{s}=(H+H^{T})/2$, where $(\cdot)^{T}$ denotes the transpose of the matrix. Now the symmetric sample $H_{s}$ looks like this:

Congratulations! You have produced your first random matrix drawn from the so-called GOE (Gaussian Orthogonal Ensemble)… a classic - more on this name later.

You can now do several things: for example, you can make the entries complex or quaternionic instead of real. In order to have real eigenvalues, the corresponding matrices need to be hermitian and self-dual respectivelyHermitian matrices have real elements on the diagonal, and complex conjugate off-diagonal entries. Quaternion self-dual matrices are $2N\times 2N$ constructed as A=[X Y; -conj(Y) conj(X)]; A=(A+A’)/2, where X and Y are complex matrices, while conj denotes complex conjugation of all entries. - better have a look at one example of the former, for $N$ as small as $N=2$

You have just met the Gaussian Unitary (GUE) and Gaussian Symplectic (GSE) ensembles, respectively - and are surely already wondering who invented these names.

We will deal with this jargon later. Just remember: the Gaussian Orthogonal Ensemble does not contain orthogonal matrices - but real symmetric matrices instead (and similarly for the others).

Although single instances can sometimes be also useful, exploring the statistical properties of an ensemble typically requires collecting data from multiple samples. We can indeed now generate $T$ such matrices, collect the $N$ (real) eigenvalues for each of them, and then produce a normalized histogram of the full set of $N\times T$ eigenvalues. With the code [$\spadesuit$ Gaussian_Ensembles_Density.m], you may get a plot like Fig. 1.1 for $T=50000$ and $N=8$.

Roughly half of the eigenvalues collected in total are positive, and half negative - this is evident from the symmetry of the histograms. These histograms are concentrated (significantly nonzero) over the region of the real axis enclosed by (for $N=8$)

You can directly jump to the end of Chapter 5 to see what these histograms look like for big matrices.

Question. Can I compute analytically the shape of these histograms? And what happens if $N$ becomes very large? $\blacktriangleright$ Yes, you can. In Chapters 10 and 12, we will set up a formalism to compute exactly these shapes for any finite $N$. In Chapter 5, instead, we will see that for large $N$ the histograms approach a limiting shape, called Wigner’s semicircle law.

Attributed to Giancarlo Rota is the statement that a random variable $X$ is neither random, nor is a variableIn the following we may use both upper and lower case to denote a random variable..

Whatever it is, it can take values in a discrete alphabet (like the outcome of tossing a die, $\{1,2,3,4,5,6\}$) or on an interval $\sigma$ (possibly unbounded) of the real line. For the latter case, we say that $\rho(x)$ is the probability density functionFor example, for the GOE matrix (1.2) the diagonal entries were sampled from the Gaussian (or normal) pdf $\rho(x)=\exp(-x^{2}/2)/\sqrt{2\pi}$. We will denote the normal pdf with mean $\mu$ and variance $\sigma^{2}$ as $N(\mu,\sigma^{2})$ in the following. (pdf) of $X$ if $\int_{a}^{b}dx\rho(x)$ is the probability that $X$ takes value in the interval $(a,b)\subseteq\sigma$.

A die will not blow up and disintegrate in the air. One of the six numbers will eventually come up. So the sum of probabilities of the outcomes should be $1$ ($=100\%$). People call this property normalization, which for continuous variables just means $\int_{\sigma}dx\rho(x)=1$.

The cumulative distribution function $F(x)$ is the probability that $X$ is smaller or equal to $x$, $F(x)=\int_{-\infty}^{x}dy\ \rho(y)$. Clearly, $F(x)\to 0$ as $x\to-\infty$ and $F(x)\to 1$ as $x\to+\infty$.

If we have two (continuous) random variables $X_{1}$ and $X_{2}$, they must be described by a joint probability density function (jpdf) $\rho(x_{1},x_{2})$. Then, the quantity $\int_{a}^{b}dx_{1}\int_{c}^{d}dx_{2}\rho(x_{1},x_{2})$ gives the probability that the first variable $X_{1}$ is in the interval $(a,b)$ and the other $X_{2}$ is, simultaneously, in the interval $(c,d)$.

When the jpdf is factorized, i.e. is the product of two density functions, $\rho(x_{1},x_{2})=\rho_{1}(x_{1})\rho_{2}(x_{2})$, the variables are said to be independent, otherwise they are dependent. When, in addition, we also have $\rho_{1}(x)=\rho_{2}(x)$, the random variables are called i.i.d. (independent and identically distributed). In any case, $\rho(x_{1})=\int\rho(x_{1},x_{2})dx_{2}$ is the marginal pdf of $X_{1}$ when considered independently of $X_{2}$.

The above discussion can be generalized to an arbitrary number $N$ of random variables. Given the jpdf $\rho(x_{1},\ldots,x_{N})$, the quantity $\rho(x_{1},\ldots,x_{N})dx_{1}\cdots dx_{N}$ is the probability that we find the first variable in the interval $[x_{1},x_{1}+dx_{1}]$, the second in the interval $[x_{2},x_{2}+dx_{2}]$, etc. The marginal pdf $\rho(x)$ that the first variable will be in the interval $[x,x+dx]$ (ignoring the others) can be computed as

Question. What is the jpdf $\rho[H]$ of the $N^{2}$ entries $\{H_{11},\ldots,H_{NN}\}$ of the matrix $H$ in (1.1)? $\blacktriangleright$ The entries in $H$ are independent Gaussian variables, hence the jpdf is factorized as $\rho[H]\equiv\rho(H_{11},\ldots,H_{NN})=\prod_{i,j=1}^{N}\left[\exp\left(-H_{ij}^{2}/2\right)/\sqrt{2\pi}\right]$.

If a set of random variables is a function of another one, $x_{i}=x_{i}(\bm{y})$, there is a relation between the jpdf of the two sets

where $J$ is the Jacobian of the transformation, given by $J(\bm{x}\to\bm{y})=\det\left(\frac{\partial x_{i}}{\partial y_{j}}\right)$. We will use this property in Chapter 6.

Chapter 2 Value the eigenvalue

In this Chapter, we start discussing the eigenvalues of random matrices.

Consider a $2\times 2$ GOE matrix H_{s}=\left(\begin{array}[]{cc}x_{1}&x_{3}\\ x_{3}&x_{2}\\ \end{array}\right), with $x_{1},x_{2}\sim N(0,1)$ and $x_{3}\sim N(0,1/2)$. What is the pdf $p(s)$ of the spacing $s=\lambda_{2}-\lambda_{1}$ between its two eigenvalues ($\lambda_{2}>\lambda_{1}$)?

The two eigenvalues are random variables, given in terms of the entries by the roots of the characteristic polynomial

therefore $\lambda_{1,2}=\left(x_{1}+x_{2}\pm\sqrt{(x_{1}-x_{2})^{2}+4x_{3}^{2}}\right)/2$ and $s=\sqrt{(x_{1}-x_{2})^{2}+4x_{3}^{2}}$.

Note that we used $\cos^{2}\theta+\sin^{2}\theta=1$ to achieve this very simple result: however, we could only enjoy this massive simplification because the variance of the off-diagonal elements was $1/2$ of the variance of diagonal elements - try to redo the calculation assuming a different ratio. Observe also that this pdf is correctly normalized, $\int_{0}^{\infty}ds\ p(s)=1$.

It is often convenient to rescale this pdf and define $\bar{p}(s)=\langle s\rangle p\left(\langle s\rangle s\right)$, where $\langle s\rangle=\int_{0}^{\infty}dsp(s)s$ is the mean level spacing. Upon this rescaling, $\int_{0}^{\infty}\bar{p}(s)ds=\int_{0}^{\infty}s\bar{p}(s)ds=1$. For the GOE as above, show that $\bar{p}(s)=(\pi s/2)\exp(-\pi s^{2}/4)$, which is called Wigner’s surmiseWhy is it defined a ’surmise’? After all, it is the result of an exact calculation! The story goes as follows: at a conference on Neutron Physics by Time-of-Flight, held at the Oak Ridge National Laboratory in 1956, people asked a question about the possible shape of the distribution of the spacings of energy levels in a heavy nucleus. E. P. Wigner, who was in the audience, walked up to the blackboard and guessed (= surmised) the answer given above., whose plot is shown in Fig. 2.1.

In spite of its simplicity, this is actually a quite deep result: it tells us that the probability of sampling two eigenvalues ’very close’ to each other ($s\to 0$) is very small: it is as if each eigenvalue ’felt’ the presence of the other and tried to avoid it (but not too much)! A bit like birds perching on an electric wire, or parked cars on a street: not too close, not too far apart. If this metaphor does not win you over, check this out .

2 Eigenvalues as correlated random variables

In the previous Chapter, we met the $N$ real eigenvalues $\{x_{1},\ldots,x_{N}\}$ of a random matrix $H$. These eigenvalues are random variables described by a jpdfWe will use the same symbol $\rho$ for both the jpdf of the entries in the upper triangle and of the eigenvalues. $\rho(x_{1},\ldots,x_{N})$.

Question. What does the jpdf of eigenvalues $\rho(x_{1},\ldots,x_{N})$ of a random matrix ensemble look like? $\blacktriangleright$ We will give it in Eq. (2.15) for the Gaussian ensemble. Not for every ensemble the jpdf of eigenvalues is known.

The important (generic) feature is that the $\{x_{i}\}$’s are not independent: their jpdf does not in general factorize. The most striking incarnation of this property is the so-called level repulsion (as in Wigner’s surmise): the eigenvalues of random matrices generically repel each other, while independent variables do not - as we show in the following section.

3 Compare with the spacings between i.i.d.’s

It is useful at this stage to consider the statistics of gaps between adjacent i.i.d. random variables. In this case, we will not see any repulsion.

Consider i.i.d. real random variables $\{X_{1},\ldots,X_{N}\}$ drawn from a parent pdf $p_{X}(x)$ defined over a support $\sigma$. The corresponding cdf is $F(x)$. The labelling is purely conventional, and we do not assume that the variables are sorted in any order.

We wish to compute the conditional probability density function $p_{N}(s|X_{j}=x)$ that, given that one of the random variables $X_{j}$ takes a value around $x$, there is another random variable $X_{k}$ ($k\neq j$) around the position $x+s$, and no other variables lie in between. In other word, a gap of size $s$ exists between two random variables, one of which sits around $x$.

The reasoning goes as follows: one of the variables sits around $x$ already, so we have $N-1$ variables left to play with. One of these should sit around $x+s$, and the pdf for this event is $p_{X}(x+s)$. The remaining $N-2$ variables need to sit either to the left of $x$ - and this happens with probability $F(x)$ - or to the right of $x+s$ - and this happens with probability $1-F(x+s)$.

Now, the probability of a gap $s$ between two adjacent particles, conditioned on the position $x$ of one variable, but irrespective of which variable this is is obtained by the law of total probability

where one uses the fact that the variables are i.i.d. and thus the probability that the particle $X_{j}$ lies around $x$ is the same for every particle, and given by $p_{X}(x)$.

To obtain the probability of a gap $s$ between any two adjacent random variables, no longer conditioned on the position of one of the variables, we should simply integrate over $x$

As an exercise, let us verify that $p_{N}(s)$ is correctly normalized, namely $\int_{0}^{\infty}ds\ p_{N}(s)=1$. We have

Changing variables $F(x+s)=u$ in the $s$-integral, and using $F(+\infty)=1$ and $du=F^{\prime}(x+s)ds=p_{X}(x+s)ds$, we get

Setting now $F(x)=v$ and using $dv=F^{\prime}(x)dx=p_{X}(x)dx$, we have

As there are $N$ variables, it makes sense to perform the ’local’ change of variables $s=\hat{s}/(Np_{X}(x))$ and consider the limit $N\to\infty$. The reason for choosing the scaling factor $Np_{X}(x)$ is that their typical spacing around the point $x$ will be precisely of order $\sim 1/(Np_{X}(x))$: increasing $N$, more and more variables need to occupy roughly the same space, therefore their typical spacing goes down. The same happens locally around points $x$ where there is a higher chance to find variables, i.e. for a higher $p_{X}(x)$.

which for large $N$ and $\hat{s}\sim\mathcal{O}(1)$, can be approximated as

the exponential law for the spacing of a Poisson process. From this, one deduces easily that i.i.d. variables do not repel, but rather attract: the probability of vanishing gaps, $\hat{s}\to 0$, does not vanish, as in the case of RMT eigenvalues!

4 Jpdf of eigenvalues of Gaussian matrices

The jpdf of eigenvalues of a $N\times N$ Gaussian matrix is given byThis jpdf goes back to the prehistory of RMT. It is an immediate consequence of Theorem 2 in , a 1939 statistics paper published in the journal Annals of Eugenics (a rather scary title, isn’t it?). In its full glory, it appeared explicitly for the first time in .

This jpdf corresponds exactly to eigenvaluesFor $\beta=4$, each matrix has $2N$ eigenvalues that are two-fold degenerate. generated according to the algorithm in Chapter 1Quite often, however, you find in the literature a Gaussian weight including extra factors, such as $\exp(-(\beta/2)\sum_{i}x_{i}^{2})$ or $\exp(-(N/2)\sum_{i}x_{i}^{2})$. One then needs to be very careful when comparing theoretical results (obtained with such conventions) to numerical simulations - in particular, a rescaling of the numerical eigenvalues by $\sqrt{\beta}$ or $\sqrt{N}$ before histogramming is essential in these two modified scenarios., and provided in the code [$\spadesuit$ Gaussian_Ensembles_Density.m].

Where does (2.15) come from? Let us postpone the proof for a while and draw some conclusions by just staring at it for a few minutes.

The Gaussian factor $e^{-\frac{1}{2}\sum_{i=1}^{N}x_{i}^{2}}$ kills any configuration of eigenvalues $\{\bm{x}\}$ where some $x_{j}$’s are “big” (far from zero, in absolute value): the eigenvalues do not like to stay too far from the origin. On the other hand, the term $\prod_{j<k}|x_{j}-x_{k}|$ kills configurations where two eigenvalues get “too close” to each other.

The “repulsion” factor $\prod_{j<k}|x_{j}-x_{k}|$ has another effect: it makes the eigenvalues strongly non-independent! Every eigenvalue feels the presence of all the others, and the jpdf (2.15) does not factorize at all. Hence, the classical tools for independent random variables are of little use here. We will use (2.15) in the next Chapter to deduce Wigner’s semicircle law in a few simple steps.

This interplay between confinement and repulsion is the physical mechanism at the heart of many results in RMT.

As a final remark, go back to the spacing pdf in Eq. (2.5), which was obtained for $N=2$ and $\beta=1$ (a $2\times 2$ GOE matrix). Armed with (2.15) one may redo the calculation as

Try to compute this integral, and recover Eq. (2.5).

Chapter 3 Classified Material

In this Chapter, we continue setting up the formalism and provide a simple classification of matrix models.

Question From the jpdf of eigenvalues $\rho(x_{1},\ldots,x_{N})$, how do I compute the shape of the histograms of the $N\times T$ eigenvalues as in Fig. 1.1, for $T$ sufficiently large? $\blacktriangleright$ To cut a long story short, all you have to do is to take the marginal $\rho(x)=\int\cdots\int dx_{2}\cdots dx_{N}\rho(x,x_{2},\ldots,x_{N})\ ,$ (3.1) and this function will reproduce the histogram profile you are after for any finite $N$. Note that $\rho(x)$ is correctly normalized to $1$, as your histogram is.

Take a single, fixed matrix $H$ with real eigenvalues - no randomness in here - and perform the following task: define a counting function $n(x)$ such that $\int_{a}^{b}n(x^{\prime})dx^{\prime}$ gives the fraction of eigenvalues $x_{i}$ between $a$ and $b$.

The way to define it is to setAs we know, the Dirac delta function (or rather distribution) $\delta(x)$ is basically an extremely peaked function at the point $x=0$, like the limit of a Gaussian pdf as its variance goes to zero, $\delta(x)=\lim_{\epsilon\to 0^{+}}\frac{1}{2\sqrt{\pi\epsilon}}e^{-x^{2}/(4\epsilon)}$.

the (normalized) sum of a set of “spikes” at the location $x_{i}$ of each eigenvalue. Using the following property of the delta function

we can show that indeed (3.2) does the job properlyCompute $N\int_{a}^{b}n(x)dx=\sum_{i=1}^{N}\int_{a}^{b}\delta(x-x_{i})dx=\sum_{i=1}^{N}\chi_{[a,b]}(x_{i})\ ,$ (3.4) where the indicator function $\chi_{[a,b]}(z)$ is equal to $1$ if $z\in(a,b)$ and otherwise. This is by definition the number of eigenvalues between $a$ and $b$, as it should..

If $H$ is now a random matrix, the function $n(x)$ becomes a random measure on the real line - a function of $x$ that changes from one realization of $H$ to another. The average of it over the set of random eigenvalues $\{x_{1},\ldots,x_{N}\}$ becomes interesting nowWe use again the shorthand $d\bm{x}=\prod_{j=1}^{N}dx_{j}$.

where $\rho(x)=\int\cdots\int dx_{2}\cdots dx_{N}\rho(x,x_{2},\ldots,x_{N})$ is the marginal density of $\rho$. Try to prove the last equality in (3.5) using the properties of delta function, and the fact that $\rho(x_{1},\ldots,x_{N})$ is symmetric upon the exchange $x_{i}\to x_{j}$. This is indeed the case for the Gaussian jpdf (2.15) and will remain generally true.

The quantity $\langle n(x)\rangle=\rho(x)$ has many names: most often, it is called the (average) spectral density. Fig. 3.1 helps you visualize how $T=4$ sets of $N=8$ randomly located “spikes” conspire to produce the continuous shape $\rho(x)=\langle n(x)\rangle$.

Question. What is the meaning of the unexpected rescaling factor $\sqrt{\beta N}$? $\blacktriangleright$ This means that the histograms of eigenvalues for larger and larger $N$ become concentrated over the interval $[-\sqrt{2\beta N},\sqrt{2\beta N}]$, in agreement with our numerical findings in Fig. 1.1. The points $\pm\sqrt{2\beta N}$ are called (spectral) edges. Note that: 1. The edges are growing with $\sqrt{N}$ - bigger matrices have a wider range of eigenvalues, can you explain why? To get histograms that do not become wider and wider with $N$, we need to divide each eigenvalue by $\sqrt{\beta N}$ before histogramming. This is what we do in Fig. 3.2, using the very same eigenvalues collected to produce Fig. 1.1. You can see that the histograms for different $\beta$s nicely collapse on top of each other, reproducing an almost perfect semielliptical shape between $-\sqrt{2}$ and $\sqrt{2}$. 2. The edges are at $\pm\sqrt{2\beta N}$ for the jpdf $\rho(x_{1},\ldots,x_{N})$ given in (2.15). If you put ad hoc extra factors in the exponential, like $\exp(-(\beta/2)\sum_{i}x_{i}^{2})$ or $\exp(-(N/2)\sum_{i}x_{i}^{2})$, as you sometimes find in the literature, this is tantamount to rescaling the eigenvalues by an appropriate factor. For example, for the choice $\exp(-(N/2)\sum_{i}x_{i}^{2})$, the edges are fixed - they do not grow with $N$ - at $\pm\sqrt{2\beta}$. 3. The edges of the semicircle are called soft: for large but finite $N$, there is always a nonzero probability of sampling eigenvalues exceeding the edge points. For example, for a GOE matrix $10\times 10$, you have a tiny but nonzero probability to sample eigenvalues larger than $\sqrt{2\beta N}\approx 4.47...$. Other ensembles have spectral densities with hard edges - this means impenetrable walls, which the eigenvalues can never cross.

2 Layman’s classification

We deal here with ensembles of square matrices with real eigenvalues (the entries can be real, complex or quaternionic random variables). Can we classify these ensembles according to simple features?

A useful scheme (covering several scenarios encountered in real life) is the following (see Fig. 3.3):

Independent entries: the first group on the left gathers matrix models whose entries are independent random variables - modulo the symmetry requirements. Random matrices of this kind are usually called Wigner matrices.

Examples: in this category, you may find adjacency matrices of random graphs , or matrices with independent power-law entries (so-called Lévy matrices ), and power-law banded matrices among others. Take a moment to download and read these papers - remember the following sentence, found on Richard Feynman’s blackboard at the time of his death: “Know how to solve every problem that has been solved”.

Rotational invariance: the second group on the right is characterized by the so-called rotational invariance. In essence, this property means that any two matrices that are related via a similarity transformation$U$ is orthogonal/unitary/symplectic if $H$ is real symmetric/complex hermitian/quaternion self-dual, respectively. You surely have noticed that this is precisely the origin of the names given to the ensembles: Orthogonal, Unitary and Symplectic. ${H}^{\prime}={U}{H}{U}^{-1}$ occur in the ensemble with the same probability

This requires the following two conditions:

$\rho[H]=\rho[{U}{H}{U}^{-1}]$. This means that the jpdf of the entries retains the same functional form before and after the transformation. This imposes a severe constraint on the allowable functional forms thanks to Weyl’s lemma , which states that $\rho[H]$ can only be a function of the traces of the first $N$ powers of $H$,

$dH_{11}\cdots dH_{NN}=dH^{\prime}_{11}\cdots dH^{\prime}_{NN}$, i.e. the flat Lebesgue measure is invariant under conjugation by $U$. This is a classical result.

The rotational invariance property in essence means that the eigenvectors are not that important, as we can rotate our matrices as freely as we wish, and still leave their statistical weight unchanged.

Examples: you may find in this category the Wishart-Laguerre (Chapter 13) and Jacobi classical ensembles, the so-called “weakly-confined” ensembles and many others. The same advice (“download-and-study”) applies here.

What about the intersection between the two classes? It turns out that it contains only the Gaussian ensembleIn its three incarnations: GOE, GUE and GSE..

This is a consequence of a theorem by Porter and Rosenzweig . And is bad news, isn’t it? We have to make a choice: if we insist that the ensemble has independent entries, then eigenvectors do matter. If we require a high level of rotational symmetry, then the entries get necessarily correlated. No free lunch (beyond the Gaussian)!

Question. I can see that the Gaussian ensemble has independent entries. But I do not easily see that it has this “rotational invariance”. $\blacktriangleright$ This can be seen from the jpdf of entries in the upper triangle (1.7). Show that you can rewrite this jpdf as $\rho[H_{s}]\propto\exp\left(-\frac{1}{2}\mbox{Tr}(H_{s}^{2})\right)\ ,$ (3.9) where $\mbox{Tr}(\cdot)$ is the matrix trace (the sum of diagonal element). For example, for the $2\times 2$ real symmetric matrix H_{s}=\left(\begin{array}[]{cc}a&b\\ b&c\\ \end{array}\right), the trace of $H_{s}$ is $a+c$, and the trace of $H_{s}^{2}$ is $a^{2}+c^{2}+\mathbf{2}b^{2}$. You can actually rewrite (1.7) as (3.9) only thanks to that factor $2$…check this! Now, from (3.9), the rotational invariance property is much easier to see: for a similarity transformation ${H_{s}}^{\prime}={U}{H_{s}}{U}^{-1}$, one has $\mbox{Tr}({H_{s}}^{\prime 2})=\mbox{Tr}(H_{s}^{2})$ (cyclic property of the trace).

3 To know more…

Anything worth mentioning beyond the above classification? One important class is represented by the biorthogonal ensembles: these are non-invariant, with non-independent entries, and yet their jpdf of eigenvalues is known in terms of the product of two determinants. Check these papers out for further information.

We suggest the following paper about “histogramming without histogramming”. Solid maths and an insightful and unconventional perspective on RMT spectra.

For a proof of the Porter-Rosenzweig theorem in the simplified $2\times 2$ case, as well as for a nice and pedagogical introduction to the Gaussian ensembles, we highly recommend the review .

For the mathematically oriented reader, who is looking for more formal classifications of random matrix models, we recommend the mini-review and references therein.

Chapter 4 The fluid semicircle

In this Chapter, we set up a statistical mechanics formalism to compute Wigner’s semicircle law for Gaussian matrices. You will learn here the so-called “Coulomb gas technique”.

The Coulomb gas (or fluid) technique is usually attributed to Dyson . Actually, a few years before, Wigner had already used it for the derivation of the semicircle law .

Take the jpdf for the Gaussian ensemble (2.15)

and rescale the eigenvalues as $x_{i}\to x_{i}\sqrt{\beta N}$.

The normalization constant now reads (set $C_{N,\beta}=(\sqrt{\beta N})^{N+\beta N(N-1)/2}$)

The factor $1/2$ in front of the logarithmic term is due to the symmetrization from $i<j$ to $i\neq j$. Stare at (4.2) intensely.

We have just exponentiated the product $\prod_{j<k}$, and obtained a canonical partition functionWe are integrating the Gibbs-Boltzmann weight $e^{-\beta N^{2}\mathcal{V}[\bm{x}]}$ over all possible positions of the particles.!

The Gibbs-Boltzmann weight $e^{-\beta N^{2}\mathcal{V}[\bm{x}]}$ corresponds to a thermodynamical fluid of particles with positions $\{x_{1},\ldots,x_{N}\}$ on a line, in equilibrium at “inverse temperature” $\beta$ under the effect of competing interactions: a quadratic (single-particle) potential (see fig. 4.1), and a repulsive (all-to-all) logarithmic term. The fluid is “static”, as there is no kinetic term in $\mathcal{V}[\bm{x}]$.

The presence of the pre-factor $\beta N^{2}$ shows - at least formally - that the limit $N\to\infty$ is a simultaneous thermodynamic and zero-temperature limit. A standard thermodynamic argument tells us how to find the equilibrium positions at zero temperature of the particles (eigenvalues) under such interactions: all we need to do is to minimize the free energy $F=-(1/\beta)\ln\mathcal{Z}_{N,\beta}$ of this system. The calculation greatly simplifies in the limit $N\to\infty$.

Question. Why is this called a “Coulomb” gas? $\blacktriangleright$ Because we have a logarithmic interaction among charged particles. More precisely, we have a 2D “fluid” of charges constrained to a line. We know that in 2D the electrostatic potential generated by a point charge is proportional to the logarithm of the distance from it - while in 3D, this potential is inversely proportional to the distance, and in 1D is proportional to the distance. Therefore, a 2D charged fluid confined to a line is not quite the same as a 1D fluid! A simple way to see this is by using Gauss’s law, with a single charge $q$ sitting at the origin on a 2D plane. If we enclose the charge in a $1$-sphere $S$ (i.e. a circle), then we must have $\int_{S}\bm{E}\cdot\bm{n}\propto q$, where $\bm{n}$ is the normal vector to the circle. If you assume that the electric field $\bm{E}$ is rotationally symmetric, i.e. $\bm{E}=E(r)\hat{\bm{r}}$, this turns into $E(r)2\pi r\propto q$, implying that $E(r)\propto q/r$. Integrating a field that goes like $1/r$ gives you a logarithmic potential.

2 Do it yourself (before lunch)

So, our goal is to find the free energy $F=-(1/\beta)\ln\mathcal{Z}_{N,\beta}$ for a large number of particles $N\to\infty$. As in many branches of physics, “larger is easier”.

We now provide a “continuum” description of the fluid, based on the following steps.

Define first a normalized one-point counting function

Instead of directly summing - or rather integrating - over all configurations of eigenvalues $\{x_{1},\ldots,x_{N}\}$, which in stat-mech we would call microstates of our fluid, we first fix a certain one-point profile $n(x)$ (non-negative, smooth and normalized).

This coarse-graining procedure can be put on slightly cleaner grounds introducing the following representation of unity as a functional integral

which enforces the definition (4.4). The functional integral runs (so to speak) over all possible normalized, non-negative and smooth functions $n(x)$. See for more details on functional integrations.

Inserting this representation of unity inside the multiple integral (4.2) and exchanging the order of integrations, we end up with

Using the identitiesProve them inserting the definition of $n(x)$ into the integrals and using properties of the delta function.

we can rewrite the two terms in the energy (4.3) as

where $\Delta(x)$ is a position-dependent short-distance cutoff. What does this mean?

4. $\mathcal{V}[\bm{x}]\to\mathcal{V}[n(x)]$

Note that in (4.9) and (4.10) the sums over eigenvalues $\{x_{1},\ldots,x_{N}\}$ have been expressed through the counting function $n(x)$, which - with a slight abuse of notation - will denote from now on its smooth limit as $N\to\infty$.

5. Evaluate the integral $I_{N}[n(x)]$ for large $N$

It is quite easy to give a physical interpretation of this multiple integral. It is basically counting how many microstates - microscopic configurations of the fluid charges - are compatible with a given macrostate - the density profile $n(x)$. We know from standard statistical mechanics arguments that the logarithm of this number should be proportional to the entropy of the fluid. Let us see how.

Introducing a ’functional’ analogue of the standard integral representation for the delta function , we can write

This type of integrals is music to the statistical physicist’s ears! It is of the form $\int d(\cdot)\exp[\Lambda f(\cdot)]$, with $\Lambda\equiv N$ a very large parameter. Hence it can be evaluated with a Laplace (or saddle-point) approximation .

Finding the critical point of the action $S[\hat{n}(x)|n(x)]$

to leading order in $N$. As expected, the term inside square brackets has precisely the form of the Shannon entropy of the density $n(x)$.

Look back again at (4.12). The short-distance cutoff $\Delta(x)$ is yet to be fixed.

A standard, physically motivated argument - going back to Dyson for charges on a ring - posits that $\Delta(x)$ - the so-called self-energy term - should be taken of the form

as the higher the density of particles around $x$, the smaller the average distance between themWe have already met a similar argument in section 2.3.. Also, $N$ charges spread over a distance of $\mathcal{O}(1)$ have a mean spacing $\sim\mathcal{O}(1/N)$, and this justifies the $1/N$ factor. This argument, however plausible, does not seem to have been made rigorous yet, though. Note, in particular, that the constant $c$ in (4.19) cannot be fixed by this simple heuristic argument. While conceptually quite important (see e.g. ), this missing bit will prove rather inconsequential in the following.

Combining (4.11), (4.12), (4.18) and (4.19), the partition function eventually reads

Note that the term $(\beta/2)N\ln N$ is essentially independent of the potential, and can be absorbed into the overall normalization constant. The $\mathcal{O}(N)$ contribution is composed by i) the self-energy term, ii) the entropic term, and iii) a contribution coming from the unknown constant $c$ in (4.19).

8. Flash-forward: cross-check with finite-$N$ result

Inserting the semicircle law into (4.21) and (4.22) - and evaluating the corresponding integrals - we obtain

Therefore, the partition function (4.20) reads for large $N$

The constants $a_{\beta}$ and $b_{\beta}$ are given as follows:

Can we check that this result is plausible?

Note that for $\beta=2$, the partition function $\mathcal{Z}_{N,\beta=2}$ from (2.16) has a particularly simple expression at finite $N$,

where $G(x)$ is a Barnes G-functionThe Barnes G-function is defined via the recursion $G(z+1)=\Gamma(z)G(z)$, with $G(1)=1$.. Hence, if everything was done correctly, the large-$N$ asymptotics of (4.29) should precisely match the large-$N$ behavior (4.25).

Using known asymptotics of the Barnes G-function, we deduce that

which coincides (up to the term $N\ln N$ included) with the asymptotics of $\mathcal{Z}_{N,\beta}$ in (4.25) once $\beta$ is set to $2$.

This check should convince you that the “mean-field” approach - based on a continuum description of the charged fluid of eigenvalues - is indeed capable of capturing the first three terms of the free energy, and only fails at the level of $\mathcal{O}(N)$ contributions - as the renormalized self-energy term cannot be precisely determined by a simple-minded scaling argument.

Let us recap what we have done so far. The normalization constant $\mathcal{Z}_{N,\beta}$ of the Gaussian model has been re-interpreted as the canonical partition function of a 2D static fluid of charged particles confined on a line, in equilibrium at inverse temperature $\beta$. For a large number of particles, among all possible configurations, the fluid will choose the one that minimizes its free energy, i.e. the logarithm of this partition function.

The partition function has been written as a functional integral over the space of normalized counting functions $n(x)$, see (4.20). For large $N$, it lends itself to a saddle-point evaluation, which will be carried out in the next Chapter.

Chapter 5 Saddle-point-of-view

Let us continue the study of the Coulomb gas method for large random matrices.

Earlier we showed that the partition function for the Gaussian model could be represented as

Quite interestingly, the leading term in the exponential is of order $\sim\mathcal{O}(N^{2})$ and not of $\sim\mathcal{O}(N)$ as in standard short-range models. As a consequence of the all-to-all coupling between the charged particles, the free energy per particle is dominated by the “energetic” component at the expenses of the “entropic” part (sub-leading for large $N$).

A saddle-point evaluation yieldsThe pre-factor $C_{N,\beta}$ has the large-$N$ behavior (4.26), whose logarithm is $\sim\mathcal{O}(N^{2}\ln N)$ and thus strictly speaking leading with respect to $N^{2}$. However, it is just an overall constant term, and the ’dynamical’ part of the free energy is of $\sim\mathcal{O}(N^{2})$.

Here, $n^{\star}(x)$ is the minimizer of the functional (5.2) in the space of normalizable and non-negative functions $n(x)$.

We set up the minimization problem by searching for the critical pointsNote that the factor $1/2$ in front of the double integral disappears because the functional differentiation picks up two counting functions, as in the integrand we have $n(x)n(x^{\prime})$. An interesting account on functional differentiation can be found at .

for $x$ in the support of $n^{\star}(x)$.

of our Coulomb gas for $N\to\infty$? It is just given by $f=\mathcal{S}[n^{\star}(x),\kappa]\equiv\mathcal{F}_{0}[n^{\star}(x)]$ - the action evaluated at the saddle-point density.

To summarize, the main task is now to find the solution of the integral equation (5.8)

2 Disintegrate the integral equation

As a preliminary observation, note that the support of $n^{\star}(x)$ (i.e. the set of $x$-values for which $n^{\star}(x)>0$) cannot be the full real line. In the limit $x\to\infty$, the integral term

- where we used normalization of the density - which is clearly incompatible with the behavior $\sim x^{2}/2$ of the known term in the equationThis is true in general for potentials growing super-logarithmically at infinity - not just for the quadratic potential corresponding to Gaussian ensembles..

Therefore, we need to look for a solution over an interval $(a,b)$ of the real line. Indeed, a rather amusing feature of this type of integral equations - of the Carleman class - is that the support over which the solution is to be found is itself unknown, and part of the problem!

Let us now first convert the integral equation into a “simpler” one.

3 Better weak than nothing

The solution to the integral equation (5.10) can be obtained by first differentiating both sides with respect to $x$. Since $\ln|x-x^{\prime}|$ is not (strictly speaking) differentiable at $x=x^{\prime}$, we consider the derivative in the weak sense.

Let $u$ be a function in $\mathcal{L}^{1}([a,b])$. We say that $v\in\mathcal{L}^{1}([a,b])$ is a weak derivative of $u$ if

for all infinitely differentiable functions $\varphi$ with $\varphi(a)=\varphi(b)=0$. The notion of weak derivative extends the standard (strong) derivative to functions that are not differentiable, but integrable in $[a,b]$. Also, if $u$ is differentiable in the standard sense, than its weak and strong derivatives coincide - just using integration by parts.

Setting $u(x)=\int dx^{\prime}n^{\star}(x^{\prime})\ln|x-x^{\prime}|$, we can write

To solve (5.14), we invoke a theorem by Tricomi , stating that

provided that $[a,b]$ is a single (compact) support and $C$ is an arbitrary constant.

Question. Who tells me that the optimal counting function $n^{\star}(x)$ is supported on a single interval $[a,b]$? $\blacktriangleright$ There is some nice physical intuition behind this. The “thermodynamical” interpretation of the eigenvalues implies that the gas of particles is confined by a quadratic well with a single minimum (see Fig. 4.1). It is then physically reasonable to foresee that the particles will fill the single minimum of the potential. If a potential has many minima, then it is possible that $n^{\star}(x)$ “splits” into as many connected components as the number of minima of the potential. Any attempt to use (5.15) in these multiple-support cases will produce unphysical solutions.

Evaluating the principal value integral with $g(t)=t$ and imposing the normalization $\int_{a}^{b}dx\ n^{\star}(x)=1$, we get

Note that the density in (5.16) is a solution of the integral equation (5.14) between $a$ and $b$ for any choice of $a$ and $b$. How to fix the “optimal” $a$ and $b$ will be the subject of the next sections.

[Of course, do not even consider trusting us on this. You are not allowed to proceed until you have derived (5.16) yourself. Sorry.]

4 Smart tricks

Now, stare at (5.16) intensely. As promised, the function $n^{\star}(x)$ (defined for $x\in(a,b)$) indeed depends on two free parameters $a$ and $b$.

We need now to compute the intensive free energy

It will of course depend as well on the two free parameters $a$ and $b$, which arose as a Phoenix from the ashes of the integral equation (5.14).

A couple of smart tricks will make our life easier. First, we would really like to get rid of the double integral in

To do that, we multiply the saddle point equation (5.10)

by $n^{\star}(x)$ and integrate over $x$. This way we obtain

Next, we fix the Lagrange multiplier $\kappa$ by setting $x=a$ in (5.19). We obtain $\kappa=a^{2}/2-\int_{a}^{b}dx\ n^{\star}(x)\ln(x-a)$. Combining everything, we get

No more $\kappa$, and no more double integrals. Nice, uh?

5 The final touch

Inserting (5.16) into (5.21) and computing the integrals with the help of an abacusIt may be useful to first change variables $z=(x-a)/(b-a)$. The resulting integrals can then be handled by most symbolic computation programs. , we obtain

We now have our (quite ugly) intensive free energy: In the code [$\spadesuit$ integral_check.m] we provide a simple numerical confirmation that the above result is equivalent to (5.21).

All we need to do is to minimize it with respect to $a$ and $b$ - the (soft) edge points of the support of $n^{\star}(x)$.

If you do that, you will obtain the solutionThe fact that the soft edges are symmetrically located around the origin is a consequence of the symmetry of the confining potential under the exchange $x\to-x$. $a=-\sqrt{2}$ and $b=\sqrt{2}$, which imply for $n^{\star}(x)$ from (5.16) the following form

the famous Wigner’s semicircle law. Very appropriate name, given that it is not the equation of a semicircle, but rather of a semi-ellipse. The code [$\spadesuit$ Tricomi_check.m] offers a numerical verification that the semicircle indeed solves equation (5.14) for $a=-b=-\sqrt{2}$.

The primitive of the integrand is - ignoring an additive constant

6 Epilogue

What is again the interpretation of the “semicircular” $n^{\star}(x)$? It is just the equilibrium profile of a gas of many charged particles on a line, which minimizes the free energy of the gas. In the “eigenvalue” language, it represents the normalized histogram of the $N$ eigenvalues of a single (very big!) instance of the Gaussian ensemble. The property that this object also faithfully represents the spectrum averaged over many samples (i.e. $n^{\star}(x)=\langle n(x)\rangle=\rho(x)$) is called self-averaging and we will assume it to hold.

The code [$\spadesuit$ Coulomb_gas.m] provides a numerical verification of what we worked on in this Chapter and the previous one. It simulates the Coulomb gas through a simple Monte Carlo procedure, which produces the equilibrium density for long enough times. Also, a numerical check of the semicircle distribution can be performed directly, i.e. through the numerical diagonalization of random matrices, with the code [$\spadesuit$ Gaussian_finite_N_rescaled.m] (see Fig. 5.1).

Note that, at the very beginning of the derivation of (5.23), we rescaled the eigenvalues by $\sqrt{\beta N}$ (Eq. (4.2)). Therefore, in the simulations we need to perform the same rescaling of our eigenvalues by $\sqrt{\beta N}$ before comparing the histogram to the theoretical semicircle. This is in agreement with the precise statement we made in the second Question in Chapter 3, namely

As a final remark, what happens if the confining potential is not quadratic? In general, if our invariant ensemble is characterized by a joint probability density of the entries of the form

then the joint law of the eigenvalues is of the form

and the analogue of the Tricomi equation for the spectral density is

Try to solve for $n^{\star}(x)$ in the case $V(x)=x-\alpha\ln x$ ($x>0$). This will correspond to the Wishart-Laguerre ensemble of random matrices, which will be extensively discussed in Chapter 13.

Question. Do all existing random matrix ensembles have the semicircle as their average spectral density? $\blacktriangleright$ Certainly not! The spectral density is highly non-universal - i.e. it strongly depends on the ensemble you consider. This said, it is true that many ensembles share it as their spectral density for large $N$. This is the case for instance of Wigner ensembles (non-invariant), when the distribution of entries decays sufficiently fast at infinity (see ).

Question. I see that the Coulomb gas treatment is insensitive to the precise value of $\beta$. But is it possible to construct an explicit random matrix ensemble $\rho[H]$, whose eigenvalues are distributed according to a Coulomb gas with $\beta\neq 1,2,4$? $\blacktriangleright$ Yes! This has been achieved by Dumitriu and Edelman , who produced ensembles of tridiagonal matrices - hence non-invariant - with independent but not identically distributed nonzero entries, whose jpdf of eigenvalues can be nevertheless computed analytically. This jpdf turns out to be equal to the Gaussian or Wishart-Laguerre ones, albeit with a continuous Dyson index $\beta>0$ (it enters as a parameter of the distribution of the nonzero entries). These ensembles are very useful also on the numerical side: they provide a much faster way to sample GXE-distributed eigenvalues (with X=O,U,S), without having to diagonalize full Gaussian matrices!

Question. If I drop the symmetry requirements on the entries of the ensemble $(H_{ij}\neq H_{ji})$, what is the resulting analogue of the semicircle law for complex eigenvalues? $\blacktriangleright$ This is called the Girko-Ginibre (or circular) law. In essence, for any sequence of random $N\times N$ matrices whose entries are i.i.d. random variables, all with mean zero and variance equal to $1/N$, the limiting spectral density is the uniform distribution over the unit disc in the complex plane.

7 To know more…

The Gaussian ensemble for $\beta=2$. The eigenvalues can be interpreted as the positions of fermions in a harmonic trap. To understand this mapping, have a look at and references therein.

Recently, the Coulomb gas technique has been improved and modified to tackle a wealth of different problems. It all started with a beautiful calculation on the following problem: what is the probability that all the eigenvalues of a Gaussian matrix are negative? Check this paper out .

The Gaussian ensembles can also come in a variant called fixed-trace: this means that one multiplies the jpdf (4.1) by $\delta\left(\sum_{i=1}^{N}x_{i}^{2}-t\right)$, which fixes the squared trace to the value $t$ (see for details).

The normalization constant $\mathcal{Z}_{N,\beta}$ for the Gaussian ensemble can be computed for finite $N$, with simple algebraic manipulations on the so called Selberg integral

It was computed by the norwegian mathematician A. Selberg, who showed that, when it exists, it is given by

To know more about recent developments in the beautiful theory of Selberg integrals, have a look at .

Chapter 6 Time for a change

In this Chapter, we show how to compute the jpdf of eigenvalues for random matrix models - whenever possible.

Suppose we have to compute the following double integrals

with $\rho_{1}(x,y)=f(x^{2}+y^{2})$ and $\rho_{2}(x,y)=xf(x^{2}+y^{2})$. Here, $f(t)$ is a function of your choice that makes both integrals convergent.

A good strategy is to make the “polar” change of variables $\{x,y\}=\{r\cos\theta,r\sin\theta\}$ to write

where $\hat{\rho}_{1}(r,\theta)=rf(r^{2})$ and $\hat{\rho}_{2}(r,\theta)=r^{2}\cos\theta f(r^{2})$. Obviously, we had to include here the extra Jacobian factor

This is all trivial and easy. But together with the following two remarks, it is all you need to know to fully understand what happens in the RMT case, with jpdf of entries and eigenvalues all over the place.

$\hat{\rho}_{1}(r,\theta)$ (the new integrand) is nothing but $\rho_{1}(r\cos\theta,r\sin\theta)\times|J(r,\theta)|$ (the old integrand, written in terms of the new variables, times the Jacobian factor) - and similarly for $\hat{\rho}_{2}$.

The marginal $\hat{\rho}_{1}(r)=\int_{0}^{2\pi}d\theta\hat{\rho}_{1}(r,\theta)$ is easier to compute than the corresponding $\hat{\rho}_{2}(r)$. This for two reasons: i) the original $\rho_{1}(x,y)$, once expressed in the new polar variables, no longer depends on one of them $(\theta)$, and ii) also the Jacobian does not depend on $\theta$. So the integration in $\theta$ becomes trivial and gives just a constant factor $2\pi$.

2 …that is the question

Take the case of real symmetric matrices for simplicity - call them $H$ instead of $H_{s}$ from now on.

Look again at the jpdf of eigenvalues (2.15) for the GOE ensemble $(\beta=1)$

How to obtain it from the jpdf of entries in the upper triangle, $\rho[H]$

In this Chapter, we provide an answer to this outstanding question.

3 Keep your volume under control

A relatively simple calculation shows that

this defines the so-called Haar measure on the orthogonal group. The Haar measure is invariant under orthogonal conjugation, and defines a probability space on orthogonal matrices. For further information, consult .

4 For doubting Thomases…

The $\{o_{ij}\}$ are real variables. The volume we are after is

where the delta functions enforce the constraints on the columns of $O$ being orthogonal with each other, and each having unit norm.

in agreement with (6.6) for $N=2$ as it should.

5 Jpdf of eigenvalues and eigenvectors

As in section 6.1 - but this time with more variables - we are after the change of variables $H\to\{\bm{x},O\}$

On the left hand side, the jpdf of the $N(N+1)/2$ entries of $H$ in the upper triangle, including the diagonal. On the right hand side, the jpdf $\hat{\rho}$ of both eigenvalues $(N)$ and independent eigenvector components ($N(N-1)/2$, the dimension of the Stiefel manifold spanned by the orthogonal group over the reals). The number of “degrees of freedom” is OK, thanks to the mind-wrecking and highly nontrivial identity $N(N+1)/2=N+N(N-1)/2$.

Clearly, on the right hand side we had to include the Jacobian of the change of variables, which we are going to compute below. While in principle this Jacobian could depend on the full set of variables $\{\bm{x},O\}$, it turns out that it only depends on the eigenvalues $\{\bm{x}\}$, exactly as it happens for the change to polar coordinates (6.3).

In our RMT case, this Jacobian is precisely the so-called Vandermonde determinantWhy this is indeed a determinant in disguise will become clearer very shortly.,

This can be generalized to the hermitian and quaternion self-dual cases. The only difference is that the Vandermonde is then raised to the power $\beta=2,4$ respectively. We will prove this in the next Chapter.

6 Leave the eigenvalues alone

Now, stare at the right hand side of (6.12) carefully.

The joint probability density of eigenvalues and eigenvectors $\hat{\rho}(x_{1},\ldots,x_{N},O)$ is the product of two terms: the jpdf of entries - written as a function of eigenvalues and eigenvectors - times the Jacobian - which is a function of the eigenvalues alone.

Then the next question is: how can I get the jpdf of eigenvalues alone? Well, you will need to integrate out the eigenvector components $\{O\}$ in (6.12). More precisely

It is certainly possible when the original jpdf of entries, once expressed in terms of eigenvalues and eigenvector components, is itself independent of eigenvectors - in complete analogy with our previous example with $r$ and $\theta$. In this case, we would get

The prototypes of this favorable case are the rotationally invariant ensembles, see the next sectionInstead, for models with independent entries, the jpdf of entries cannot be - in general - written in terms of the eigenvalues alone. For such models, the jpdf of eigenvalues is therefore not generally known..

7 For invariant models…

We can now formulate a cute little theorem for invariant models . The proof is given below.

Note that there is no absolute value around the Vandermonde, as the eigenvalues are ordered.

Let us see how this theorem works in practice for the GOE case. We have

which needs to be compared with Eq. (2.15) for $\beta=1$ - given without proof at the time

Do the two equations (6.18) and (6.19) indeed agree, as they should? Almost.

Notice that (6.18) holds for ordered eigenvalues, while (6.19) holds for unordered eigenvalues (hence the need to include the absolute value). The two normalization constants differ indeed by a factor $N!$

8 The proof

in (6.16) come from? It is instructive to look at this derivation more closely.

Recall from (6.14) and (6.15) that (for the favorable case where one can integrate out the eigenvectors)

Chapter 7 Meet Vandermonde

The “repulsive” term between eigenvalues of invariant models $\prod_{i<j}(x_{j}-x_{i})$ can be written as a determinant, called Vandermonde in honor of the French mathematician Alexandre-Théophile Vandermonde (who never wrote it ).

The Vandermonde is clearly a completely anti-symmetric polynomial in $N$ variables: take for example $N=3$. We have $\Delta_{3}(\bm{x})=(x_{2}-x_{1})(x_{3}-x_{1})(x_{3}-x_{2})$. Now, exchange any two $x_{j}$s: for example, $x_{3}\leftrightarrow x_{2}$. We get $-\Delta_{3}(\bm{x})$ (we pick up a minus sign any time we make any exchange of two $x_{j}$s).

The Vandermonde has a quite funny property: we can understand it already on a $2\times 2$ matrix. Take

Stare at these two determinants carefully. We have just replaced the second row of the first matrix (containing first powers of $x_{1}$ and $x_{2}$) with a first degree polynomial. The result is just $3$ times the Vandermonde on the left. The $17$ has disappeared altogether! This means that you have a lot of freedom in devising a matrix whose determinant gives the Vandermonde.

More formally, the entries $x_{i}^{k}$ in the $(k+1)$th row can be replaced, up to a constant factor $a_{0}a_{1}\cdots a_{N-1}$, by a polynomial of degree $k$ of the form: $\pi_{k}(x_{i})=a_{k}x_{i}^{k}+\cdots$, where we omit terms of lower order in $x_{i}$. The important point is that these lower order terms can be absolutely anything. The result is that:

Orthogonal polynomials are an important class of polynomials $\pi_{k}(x)$ that can be especially useful to play this trick. We will discuss in Chapter 10 how this simple property can actually turn seemingly impossible calculations into feasible ones.

For instance, let us show how the Hermite and Laguerre orthogonal polynomials can be used to express the Vandermonde. For $N=3$ it is easy to see that

2 Do it yourself

We now derive the nontrivial relation (6.13) $J(H\to\{\bm{x},O\})=\Delta_{N}(\bm{x})$ for real symmetric matrices $H$. We stress that this proof does not require any assumption on the rotational invariance of the ensemble.

Pulling out a factor ${O}$ to the left and ${O}^{T}$ to the right we obtain $\delta{H}={O}(\delta\hat{H}){O}^{T}$, where

Here, $\delta{\Omega}={O}^{T}\delta{O}$ is an antisymmetric matrixObviously, you need to prove it before proceeding.. Since $\delta{H}$ and $\delta\hat{H}$ are related via an orthogonal transformation, we only have to find the Jacobian of $\delta\hat{H}\to\{\delta X,\delta{\Omega}\}$.

Noting that $\delta X$ is diagonal, we can write

This is equivalent to the following differential relations:

Don’t you see the Vandermonde trying hard to crop up here ☺?

Let us now construct the Jacobian matrix $J$ for a concrete $3\times 3$ case. The generalization to the $N\times N$ case will then appear obvious. The matrix $J$ has dimension $\frac{N(N+1)}{2}$, so it is a $6\times 6$ matrix for $N=3$. We parametrize the antisymmetric matrix $\delta\Omega$ as follows:

Swapping rows and columns, it is possible to bring this to the diagonal form, so that the determinant becomes trivial to compute. In the general $N$ case, one has:

as expected. The proof in the complex hermitian and quaternion self-dual cases is analogous and is left as an exercise.

For a nice numerical test of the Jacobian identity (7.19), we refer to , Section 3.2, while for a “back-of-the-envelope” derivation based on counting degrees of freedom, see .

We will make extensive use of the Vandermonde determinant and its properties in Chapter 10.

Chapter 8 Resolve(nt) the semicircle

If $H$ is a random matrix, then $G_{N}(z)$ is a random complex function that has poles at the locations $x_{i}$ of each eigenvalue.

The second ingredient we need is the Sokhotski-Plemelj formula

which should be interpreted as the integral relation (for a real-valued test function $\varphi(x)$ such that the integrals make sense)

For a one-liner proof, see below (around (8.6)).

Question. What is the point of introducing this identity? $\blacktriangleright$ First, stare at (8.2) carefully. You see that, on the right hand side, the imaginary part is just a delta function. So, this identity is (yet another) way of representing a delta function, as the imaginary part of a rational function (the left hand side). Knowing that the spectral density is defined in terms of a delta function $\rho(x)=\langle(1/N)\sum_{i=1}^{N}\delta(x-x_{i})\rangle$, you should be spotting an interesting connection here. More on this later.

2 Averaging

Imagine now to take the limit $N\to\infty$ of $\langle G_{N}(z)\rangle$, where we average over the distribution of the matrix $H$. This average is called resolvent, or $Green^{\prime}sfunction$, or Stieltjes transform. It is natural to assume (and can be mathematically justified) that:

the poles at $x_{i}$ merge into a continuous “cut” on the real line,

we have to “weigh” the integrand with the average density of eigenvalues $\rho(x)$ at point $x$.

The cut on the real line is therefore nothing but the support of the spectral density, and the average resolvent is defined for all complex values $z$ outside this cut (for example, outside the interval $[-\sqrt{2},\sqrt{2}]$ on the real line for the Gaussian ensemble).

If you are inclined to believe that (8.5) is very plausible (to say the least), we can now proceed smoothly.

So, if you know (or can calculate) the resolvent in the complex plane, you can from it deduce the spectral density.

All this in theory. Practice in the next section.

3 Do it yourself

We propose here a truly elementary derivation of the algebraic equation satisfied by the resolvent for the Gaussian ensemble.

Consider the partition function of the standard Gaussian ensemble, after a further rescaling $x_{i}\to x_{i}\sqrt{N}$ and ignoring prefactors

Compared to our earlier Coulomb gas treatment, we have pulled out a factor $N$ (not $N^{2}$), so that the $x_{i}$ are now of $\mathcal{O}(1)$ for large $N$. Instead of introducing a continuous counting function $n(x)$ (as we did in Chapter 4), we can directly perform the saddle point evaluation of the $N$-fold integral (8.11), obtaining for each variable $x_{i}$ the equation

Multiplying (8.13) by $\frac{1}{N(z-x_{i})}$ and summing over $i$, we get:

Adding and subtracting $z$ in the numerator, the left-hand-side $L$ becomes

As for the right-hand-side, let us define $R=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j\neq i}\frac{1}{z-x_{i}}\frac{1}{x_{i}-x_{j}}\ .$ Writing

one obtains the following self-consistency equation for $R$

Equating $L$ to $R$, we obtain as promised that the saddle-point condition (8.13) gets converted into an equation for the resolvent

This is good, but is still a differential equation for $G_{N}(z)$, while we promised an even simpler algebraic equation. It is actually easy to get rid of the differential term in (8.18) by noticing that, with $x_{i}$ of $\mathcal{O}(1)$, the resolvent as defined in (8.1) is itself of $\mathcal{O}(1)$ and therefore the term $\frac{1}{2N}G_{N}^{\prime}(z)$ is subleading for large $N$.

Taking the average, the surviving algebraic (at long last!) equation for $N\to\infty$ reads

It is instructive to solve (8.19) directly as a quadratic equation (recall that quadratic equations for complex variables admit the same solving formula as their real counterparts), yielding

From this expression, you see that i) for $|x|>\sqrt{2}$ you obtain that the density is , and ii) for $|x|<\sqrt{2}$, you need to select the $(-)$ or $(+)$ sign in front, according to whether $x>0$ or $x<0$ respectively. After choosing the right sign, you get $\rho(x)=(1/\pi)\sqrt{2-x^{2}}$ as expected.

4 Localize the resolvent

Let us now take a step back and “unpack” the definition of the resolvent in equation (8.1).

and $c_{i}^{j}$ is the $i$th component of the normalized eigenvector associated with the $j$th eigenvalue of $H$.

So, you may now ask, what should I make of these matrix elements? Well, it turns out that they contain precious information about the localization properties of the matrix ensemble they are associated with.

Simply put, the term localization refers to how “spread out” over their components the eigenvectors of a matrix are. Let us define the inverse participation ratio (IPR) of a normalized eigenvector as

Now, when the eigenvector’s components are all roughly of the same magnitude, then we must have $c_{i}^{j}\simeq 1/\sqrt{N}$, $\forall\ i$, due to normalization. Hence, we will have $\mathcal{I}_{N,j}\simeq 1/N$, and the IPR will vanish in the large $N$ limit.

If, on the other hand, the eigenvector is significantly different from zero only on a number $s$ of sites, then for those sites we will have $c_{i}^{j}\simeq 1/\sqrt{s}$, and the IPR will remain roughly equal to $1/s$ in the large $N$ limit. So, all in all, the IPR is a handy tool that tells us whether certain eigenvectors of a matrix are extended (i.e. have an extensive number of non zero components) or instead localized on a finite number of sites.

Although this may sound like a mathematical curiosity, the localization properties of matrix ensembles are related to a number of relevant features of the physical systems they describe. In particular, it is often crucial to detect the so called mobility edge, i.e. the critical eigenvalue that separates the part of the spectrum associated with extended states from the one associated with localized states. For example, it has famously been shown that the mobility edge determines the Anderson transition in electronic systems .

All in all, it should be now clear that having analytical access to the distributional properties of the IPRs corresponding to different segments of a given ensemble’s eigenvalue spectrum is a valuable thing. Luckily, this is where our diagonal elements (8.23) come to the rescue. Indeed, it has been shown in that the average value $P(x)$ of IPRs associated with states whose corresponding eigenvalues lie between $x$ and $x+dx$ can be written in the large $N$ limit as

5 To know more…

The saddle-point evaluation (8.13) based on the partition function (8.11) is clearly valid when the neglected terms in the exponent are indeed subleading $(\mathcal{O}(N))$. There are models - rotationally invariant by construction - where the Dyson index $\beta$ is allowed to scale with $N$ . These models provide explicit realizations of invariant $\beta$-ensembles, for which the resolvent equation is necessarily more involved. Ref. is also suggested for an elementary derivation of this “improved” resolvent equation in the presence of a hard wall in the spectrum.

Matrix models such as the Gaussian can be constructed introducing a fictitious time evolution (stochastic) of the entries. In this case, it is possible to show that the resolvent satisfies a partial differential equation of the Burgers type (see the beautiful paper ).

The equation for the resolvent can be given a pretty interpretation in terms of planar diagrams. Diagrammatic methods are at the heart of many beautiful results in RMT (see and ).

Chapter 9 One pager on eigenvectors

Take the GUE ensemble of $N\times N$ hermitian matrices. Any given matrix in the ensemble will have unit-norm eigenvectors having in general complex components. What is the statistics of such components?

Since eigenvalues and eigenvectors of invariant matrix models are decoupled, the only constraint on the $N$ components of an eigenvector is that its norm must be one, therefore their jpdf reads

where $C_{N}$ is a normalization constant.

It is convenient to compute the marginal distribution of a single component, say $|c_{1}|^{2}$, given by

Similarly, we can compute the jpdf of eigenvector components (this time all real numbers) of a GOE matrix.

The calculation in (9.2) is carried out by first defining an auxiliary object

such that $P_{GUE}(y)=P_{GUE}(y;1)$. Then, taking the Laplace transform with respect to $t$ to kill the delta function in (9.3)

and finally converting the 2d integrals in polar coordinates

where we have absorbed the angular constants in the overall normalization.

Inverting the Laplace transform, we obtain

where $\theta(z)$ is the Heaviside step function. Setting $t=1$ and normalizing, we obtain

Computing the average $\langle y\rangle$ in both cases

leads us to consider the scaled variable $\eta=yN$ and take the limit $N\to\infty$. This produces the scaled densities

The first of these densities is called the Porter-Thomas distribution . Note also that the Gaussian nature of the matrix ensembles has not been used anywhere in the derivation (the same densities would be obtained for any orthogonal or unitary ensemble).

The study of eigenvectors of random matrices has been recently revived due to their importance in quantum systems (see, e.g.,)

Chapter 10 Finite N𝑁N

Look back at Chapter 1, where we constructed Gaussian matrices and histogrammed their eigenvalues. For $N\to\infty$, we showed in various ways that the average spectral density converges to the semicircle law. But what happens for finite $N$? Can we compute analytically the shape of the histogram for, say, a $13\times 13$ Gaussian matrix? The answer is Yes - and not only for Gaussian matrices, but for any rotationally invariant ensemble! This is done here. We start from the case $\beta=2$, as it is much easier.

Already in Chapter 7, we mentioned that the Vandermonde determinant has some funny properties: in particular, each row in the Vandermonde matrix can be replaced by a polynomial of suitable degree, with many a priori unspecified coefficients. The freedom in choosing these polynomials is enormous. A judicious choice is the key of the celebrated orthogonal polynomial technique.

Take the jpdf of the $N$ real eigenvalues of a rotationally invariant ensemble with $\beta=2$

which is written in the ‘potential’ form (see eq. (5.30)). For example, for the Gaussian ensemble $V(x)=x^{2}/2$.

What is the goal then? To compute the average spectral density for finite $N$, i.e. the $N-1$-fold integral

where the partition function is $\mathcal{Z}_{N}=\int dx_{1}\cdots dx_{N}\prod_{i=1}^{N}e^{-V(x_{i})}|\Delta_{N}(\bm{x})|^{2}$.

Note that in (10.2) we are integrating over all variables but one. These integrals are nasty, though! The integrand does not factorize at all, so we need to find some smart trick to carry out the integration. It took a while even to the pioneers of these calculations (for instance, Gaudin and Mehta) to figure out how to proceed. The steps are as follows:

Rewrite the Vandermonde $\Delta_{N}(\bm{x})$ as a determinant of the matrix $A$, whose entries are polynomials $\pi_{k}(x)$ (to be determined), as in (7.3)

Use the general relationHereafter, inside a determinant the indices of the entries will run from $1$ to $N$.

applied to the matrix $A$ from step $1$, to write

where $\phi_{i}(x)=e^{-V(x)/2}\pi_{i}(x)$ and

which is a central object in RMT: the kernel.

Choose judiciously the (so far undetermined) polynomials $\pi_{j}(x)$. A great choice is to pick them orthonormal with respect to the weightNote that there is a factor $(1/2)$ multiplying $V(x)$ in the kernel (10.7), while there is none in the weight function of the orthonormal polynomials in (10.8). $\exp(-V(x))$

For instance, for the Gaussian (unitary) ensemble ($V(x)=x^{2}/2$) the corresponding orthonormal polynomials are

if $H_{j}(x)$ are Hermite polynomials satisfying $\int_{-\infty}^{\infty}dx\ H_{j}(x)H_{k}(x)\exp(-x^{2})=\sqrt{\pi}2^{j}j!\delta_{jk}$.

Question. What is the advantage of choosing polynomials with this “orthonormality” property? $\blacktriangleright$ Well, the reason is that the kernel $K_{N}(x,x^{\prime})$ in (10.7), if the polynomials are chosen this way, satisfies a quite amazing “reproducing” property $\int dyK_{N}(x,y)K_{N}(y,x^{\prime})=K_{N}(x,x^{\prime})\ .$ (10.10) The proof is very simple: just insert (10.7) into (10.10) and use the orthonormality relation (10.8). This property has a quite unexpected consequence, which eventually allows to carry out the multiple integrations in (10.2) in a very elegant, iterative way. Another ingredient is necessary, though, and is presented in the next Section.

2 Integrating inwards

Summarizing, we have to carry out the multiple integration in (10.2) over a jpdf, which can be written as the determinant of a kernel (see (10.6)), something like

In normal situations, this would seem a rather hopeless task. But the reproducing property of the kernel offers an unexpected way around.

First, an illuminating $2\times 2$ example, and then the full-fledged (though dry) theory. Imagine the following $2\times 2$ matrix $J_{2}(\bm{x})$, depending on the vector $\bm{x}=\{x_{1},x_{2}\}$ through a function $f(x,y)$ as follows

Suppose now that the function $f$ satisfies the ‘’reproducing” property (10.10), namely $\int f(x,y)f(y,z)\,d\mu(y)=f(x,z)$ for a certain measure $\mu(y)$. What happens to the following integral

where $q=\int d\mu(x_{2})f(x_{2},x_{2})$. We used the reproducing property to evaluate the second integral.

Maybe this short calculation is not particularly revealing, but it can be actually extended to the $N\times N$ case as follows: let $J_{N}(\bm{x})$ be an $N\times N$ matrix whose entries depend on a real vector $\bm{x}=(x_{1},x_{2},\ldots,x_{N})$ and have the form $J_{ij}=f(x_{i},x_{j})$, where $f$ is some function satisfying the “reproducing kernel” property $\int f(x,y)f(y,z)\,d\mu(y)=f(x,z)\ ,$ for some measure $d\mu(y)$. Then the following holds: