Let {n} be a large natural number, and let {M_n} be a matrix drawn from the Gaussian Unitary Ensemble (GUE), by which we mean that {M_n} is a Hermitian matrix whose upper triangular entries are iid complex gaussians with mean zero and variance one, and whose diagonal entries are iid real gaussians with mean zero and variance one (and independent of the upper triangular entries). The eigenvalues {\lambda_1(M_n) \leq \ldots \leq \lambda_n(M_n)} are then real and almost surely distinct, and can be viewed as a random point process {\Sigma^{(n)} := \{\lambda_1(M_n),\ldots,\lambda_n(M_n)\}} on the real line. One can then form the {k}-point correlation functions {\rho_k^{(n)}: {\bf R}^k \rightarrow {\bf R}^+} for every {k \geq 0}, which can be defined by duality by requiring

\displaystyle \mathop{\bf E} \sum_{i_1,\ldots,i_k \hbox{ distinct}} F( \lambda_{i_1}(M_n),\ldots,\lambda_{i_k}(M_n))

\displaystyle = \int_{{\bf R}^k} \rho_k^{(n)}(x_1,\ldots,x_k) F(x_1,\ldots,x_k)\ dx_1 \ldots dx_k

for any test function {F: {\bf R}^k \rightarrow {\bf R}^+}. For GUE, which is a continuous matrix ensemble, one can also define {\rho_k^{(n)}(x_1,\ldots,x_k)} for distinct {x_1<\ldots<x_k} as the unique quantity such that the probability that there is an eigenvalue in each of the intervals {[x_1,x_1+\epsilon],\ldots,[x_k,x_k+\epsilon]} is {(\rho_k^{(n)}(x_1,\ldots,x_k)+o(1))\epsilon^k} in the limit {\epsilon\rightarrow 0}.

As is well known, the GUE process is a determinantal point process, which means that {k}-point correlation functions can be explicitly computed as

\displaystyle \rho^{(n)}_k(x_1,\ldots,x_k) = \det( K^{(n)}(x_i,x_j) )_{1 \leq i,j \leq k}

for some kernel {K^{(n)}: {\bf R} \times {\bf R} \rightarrow {\bf C}}; explicitly, one has

\displaystyle K^{(n)}(x,y) := \sum_{k=0}^{n-1} P_k(x) e^{-x^2/4}P_k(y) e^{-y^2/4}

where {P_0, P_1,\ldots} are the (normalised) Hermite polynomials; see this previous blog post for details.

Using the asymptotics of Hermite polynomials (which then give asymptotics for the kernel {K^{(n)}}), one can take a limit of a (suitably rescaled) sequence of GUE processes to obtain the Dyson sine process, which is a determinantal point process {\Sigma} on the real line with correlation functions

\displaystyle \rho_k(x_1,\ldots,x_k) = \det( K(x_i,x_j) )_{1 \leq i,j \leq k} \ \ \ \ \ (1)

where {K} is the Dyson sine kernel

\displaystyle K(x,y) := \frac{\sin(\pi(x-y))}{\pi(x-y)}. \ \ \ \ \ (2)

A bit more precisely, for any fixed bulk energy {-2 < u < 2}, the renormalised point processes {\rho_{sc}(u) \sqrt{n} ( \Sigma^{(n)} - \sqrt{n} u )} converge in distribution in the vague topology to {\Sigma} as {n \rightarrow \infty}, where {\rho_{sc}(u) := \frac{1}{2\pi} (4-u^2)^{1/2}_+} is the semi-circular law density.

On the other hand, an important feature of the GUE process {\Sigma^{(n)} = \{\lambda_1,\ldots,\lambda_n\}} is its stationarity (modulo rescaling) under Dyson Brownian motion

\displaystyle d\lambda_i = dB_i + \sum_{j \neq i} \frac{dt}{\lambda_i-\lambda_j}

which describes the stochastic evolution of eigenvalues of a Hermitian matrix under independent Brownian motion of its entries, and is discussed in this previous blog post. To cut a long story short, this stationarity tells us that the self-similar {n}-point correlation function

\displaystyle \rho^{(n)}_n(t,x) := t^{-n/2} \rho^{(n)}_n(x/\sqrt{t})

obeys the Dyson heat equation

\displaystyle \partial_t \rho^{(n)}_n = \frac{1}{2} \sum_{i=1}^n \partial_{x_i}^2 \rho^{(n)}_n - \sum_{1 \leq i,j \leq n;i\neq j} \partial_{x_i} \frac{\rho^{(n)}_n}{x_i-x_j}

(see Exercise 11 of the previously mentioned blog post). Note that {\rho^{(n)}_n} vanishes to second order whenever two of the {x_i} coincide, so there is no singularity on the right-hand side. Setting {t=1} and using self-similarity, we can rewrite this equation in time-independent form as

\displaystyle -\frac{1}{2} \sum_{i=1}^n \partial_i (x_i \rho^{(n)}_n) = \frac{1}{2} \sum_{i=1}^n \partial_{x_i}^2 \rho^{(n)}_n - \sum_{1 \leq i,j \leq n;i\neq j} \partial_{x_i} \frac{\rho^{(n)}_n}{x_i-x_j}.

One can then integrate out all but {k} of these variables (after carefully justifying convergence) to obtain a system of equations for the {k}-point correlation functions {\rho^{(n)}_k}:

\displaystyle -\frac{1}{2} \sum_{i=1}^k \partial_i (x_i \rho^{(n)}_k) = \frac{1}{2} \sum_{i=1}^k \partial_{x_i}^2 \rho^{(n)}_k - \sum_{1 \leq i,j \leq k;i\neq j} \partial_{x_i} \frac{\rho^{(n)}_k}{x_i-x_j}

\displaystyle - \sum_{i=1}^k \partial_{x_i} \int_{\bf R} \frac{\rho^{(n)}_{k+1}(x_1,\ldots,x_{k+1})}{x_i-x_{k+1}}\ dx_{k+1},

where the integral is interpreted in the principal value case. This system is an example of a BBGKY hierarchy.

If one carefully rescales and takes limits (say at the energy level {u=0}, for simplicity), the left-hand side turns out to rescale to be a lower order term, and one ends up with a hierarchy for the Dyson sine process:

\displaystyle 0 = \frac{1}{2} \sum_{i=1}^k \partial_{x_i}^2 \rho_k - \sum_{1 \leq i,j \leq k;i\neq j} \partial_{x_i} \frac{\rho_k}{x_i-x_j} \ \ \ \ \ (3)

\displaystyle - \sum_{i=1}^k \partial_{x_i} \int_{\bf R} \frac{\rho_{k+1}(x_1,\ldots,x_{k+1})}{x_i-x_{k+1}}\ dx_{k+1}.

Informally, these equations show that the Dyson sine process {\Sigma = \{ \lambda_i: i \in {\bf Z} \}} is stationary with respect to the infinite Dyson Brownian motion

\displaystyle d\lambda_i = dB_i + \sum_{j \neq i} \frac{dt}{\lambda_i-\lambda_j}

where {dB_i} are independent Brownian increments, and the sum is interpreted in a suitable principal value sense.

I recently set myself the exercise of deriving the identity (3) directly from the definition (1) of the Dyson sine process, without reference to GUE. This turns out to not be too difficult when done the right way (namely, by modifying the proof of Gaudin’s lemma), although it did take me an entire day of work before I realised this, and I could not find it in the literature (though I suspect that many people in the field have privately performed this exercise in the past). In any case, I am recording the computation here, largely because I really don’t want to have to do it again, but perhaps it will also be of interest to some readers.

The basic tool here is cofactor expansion, which we write as follows. Given an {n \times n} matrix {A_n = (a_{i,j})_{1 \leq i,j \leq n}}, we can expand the determinant {\det(A_n)} as

\displaystyle \det(A_n) = a_{n,n} \det(A_{n-1}) - \sum_{i=1}^{n-1} a_{i,n} \det( R_i( A_{n-1}, X ) ) \ \ \ \ \ (4)

where {X := (a_{n,1},\ldots,a_{n,n-1})} is the bottom left row of {A_n}, {A_{n-1} := (a_{i,j})_{1 \leq i,j \leq n-1}} is the top left minor of {A_{n-1}}, and {R_i(A_{n-1},X)} is the {n-1 \times n-1} matrix formed by replacing the {i^{th}} row of {A_{n-1}} by {X}.

As an example of how this cofactor expansion is used in determinantal point processes, let us recall Gaudin’s lemma, which we state slightly informally:

Lemma 1 (Gaudin’s lemma) Let {K: {\bf R} \times {\bf R} \rightarrow {\bf R}} be a (sufficiently nice) function obeying the idempotent relation

\displaystyle K(x,z) = \int_{\bf R} K(x,y) K(y,z)\ dy \ \ \ \ \ (5)

and the trace formula

\displaystyle \int_{\bf R} K(x,x)\ dx = n. \ \ \ \ \ (6)

For each {k}, let {\rho_k} be the correlation function

\displaystyle \rho_k(x_1,\ldots,x_k) := \det( K(x_i,x_j) )_{1 \leq i,j \leq k}.

Then we have

\displaystyle \int_{\bf R} \rho_{k+1}(x_1,\ldots,x_k,x_{k+1})\ dx_{k+1} = (n-k) \rho_k(x_1,\ldots,x_k).

Proof: From (4) we have

\displaystyle \rho_{k+1}(x_1,\ldots,x_k,x_{k+1}) = K(x_{k+1},x_{k+1}) \rho_k(x_1,\ldots,x_k) \ \ \ \ \ (7)

\displaystyle - \sum_{i=1}^{k} K(x_i,x_{k+1}) \det( R_i( A_k, X(x_{k+1}) ) )

where {A_k := (K(x_i,x_j))_{1 \leq i,j \leq k}} and {X(x_{k+1}) := (K(x_{k+1},x_1),\ldots,K(x_{k+1},x_k))}. Integrating using (6), we conclude that

\displaystyle \int_{\bf R} \rho_{k+1}(x_1,\ldots,x_k,x_{k+1})\ dx_{k+1} = n \rho_k(x_1,\ldots,x_k)

\displaystyle - \sum_{i=1}^{k} \int_{\bf R} K(x_i,x_{k+1}) \det( R_i( A_k, X(x_{k+1}) )\ dx_{k+1} ).

By the multilinearity of determinant, we can write

\displaystyle \int_{\bf R} K(x_i,x_{k+1}) \det( R_i( A_k, X(x_{k+1}) ) )\ dx_{k+1}

\displaystyle = \det( R_i(A_k, \int_{\bf R} K(x_i,x_{k+1}) X(x_{k+1})\ dx_{k+1} ) ).

But from (5) we have

\displaystyle \int_{\bf R} K(x_i,x_{k+1}) X(x_{k+1})\ dx_{k+1} = (K(x_i,x_1),\ldots,K(x_i,x_k))

and so

\displaystyle R_i(A_k, \int_{\bf R} K(x_i,x_{k+1}) X(x_{k+1})\ dx_{k+1} ) = A_k,

and the claim follows. \Box

To establish (3), we will need some identities that are specific to the sine kernel (2) that are reminiscent of (5) and (6). The analogue of (6) is very easy:

\displaystyle K(x,x) = 1. \ \ \ \ \ (8)

The sine kernel actually does obey (5), but we will need the following more general identity:

Lemma 2 For any real numbers {x,y,z}, one has

\displaystyle \int_{\bf R} K(x,t) K(y,t) \frac{dt}{z-t} = \frac{K(x,y) - \cos(\pi(z-y)) K(x,z)}{z-y} \ \ \ \ \ (9)

(using the principal value interpretation of the integral) when {y,z} are distinct. In particular, by l’Hopital’s rule one has

\displaystyle \int_{\bf R} K(x,t) K(y,t) \frac{dt}{y-t} = - \partial_y K(x,y) \ \ \ \ \ (10)

Note that if one multiplies by {z} in (9) and then sends {z \rightarrow \infty}, one recovers (5).

Proof: This can be verified by a routine contour integration, but we will give a more Fourier analytic proof instead. We will omit some (routine) details concerning justifying the various interchanges of integrals we will be using here. By a limiting argument we may take {x,y,z} to be distinct. Observe that

\displaystyle K(x,y) = \int_{-1/2}^{1/2} e^{2\pi i \alpha(x-y)}\ d\alpha \ \ \ \ \ (11)

and so the left-hand side of (9) can be rewritten as

\displaystyle \int_{-1/2}^{1/2} \int_{-1/2}^{1/2} e^{2\pi i (\alpha x + \beta y)} \int_{\bf R} \frac{e^{-2\pi i (\alpha+\beta) t}}{z-t}\ dt d\alpha d\beta.

From the basic identity {\int_{\bf R} \frac{\sin x}{x}\ dx = \frac{\pi}{2}} we have

\displaystyle \int_{\bf R} \frac{e^{-2\pi i \theta t}}{z-t}\ dt = \hbox{sgn}(\theta) \pi i e^{-2\pi i \theta z}

for any {\theta \in {\bf R}}, so the preceding expression simplifies to

\displaystyle \pi i \int_{-1/2}^{1/2} \int_{-1/2}^{1/2} e^{2\pi i (\alpha (x-z) + \beta (y-z))} \hbox{sgn}(\alpha+\beta)\ d\alpha d\beta.

We split this as the difference of

\displaystyle \pi i \int_{\bf R} \int_{-1/2}^{1/2} e^{2\pi i (\alpha (x-z) + \beta (y-z))} \hbox{sgn}(\alpha+\beta)\ d\alpha d\beta. \ \ \ \ \ (12)

and

\displaystyle \pi i \int_{|\beta| > 1/2} \int_{-1/2}^{1/2} e^{2\pi i (\alpha (x-z) + \beta (y-z))} \hbox{sgn}(\alpha+\beta)\ d\alpha d\beta. \ \ \ \ \ (13)

For (12), we make the change of variables {\gamma := \alpha+\beta} and rewrite this as

\displaystyle \pi i \int_{-1/2}^{1/2} e^{2\pi i \alpha (x-y)} \int_{\bf R} e^{2\pi i \gamma (y-z)} \hbox{sgn}(\gamma)\ d\gamma d\alpha.

We have (in the distributional sense at least) that

\displaystyle \pi i \int_{\bf R} e^{2\pi i \gamma (y-z)} \hbox{sgn}(\gamma)\ d\gamma = \frac{1}{z-y}

and so by (9) this term simplifies to {\frac{K(x,y)}{z-y}}.

Now we turn to (13). Here we observe that in the region of integration we have {\hbox{sgn}(\alpha+\beta)=\hbox{sgn}(\beta)}, and so this expression can be factored as

\displaystyle \pi i (\int_{|\beta| > 1/2} e^{2\pi i \beta(y-z)} \hbox{sgn}(\beta)\ d\beta) (\int_{-1/2}^{1/2} e^{2\pi i \alpha(x-z)}\ d\alpha).

The final factor here is {K(x,z)} thanks to (9), while a brief computation shows (in the distributional sense, at least) that

\displaystyle \pi i (\int_{|\beta| > 1/2} e^{2\pi i \beta(y-z)} \hbox{sgn}(\beta)\ d\beta) = \frac{\cos(\pi(z-y))}{z-y}

and the claim follows.

Now we turn to the verification of (3). It will suffice to establish the identity

\displaystyle \int_{\bf R} \frac{\rho_{k+1}(x_1,\ldots,x_{k+1})}{x_i-x_{k+1}}\ dx_{k+1} = \frac{1}{2} \partial_{x_i} \rho_k(x_1,\ldots,x_k) \ \ \ \ \ (14)

\displaystyle - \sum_{1 \leq j \leq k;i\neq j} \frac{\rho_k(x_1,\ldots,x_k)}{x_i-x_j}

for each {i=1,\ldots,k}. By a continuity argument we may take the {x_1,\ldots,x_k} to be distinct.

Fix {i}. By (7), the left-hand side of (14) can be written as the difference of

\displaystyle \rho_k(x_1,\ldots,x_k) \int_{\bf R} \frac{K(x_{k+1},x_{k+1})}{x_i-x_{k+1}}\ dx_{k+1}

and

\displaystyle \sum_{j=1}^{k}\det( R_j( A_k, \int_{\bf R} K(x_j,x_{k+1}) X(x_{k+1}) \frac{dx_{k+1}}{x_i-x_{k+1}}) ).

The first expression vanishes thanks to (8), so we turn to the second. From Lemma 2 we have

\displaystyle \int_{\bf R} K(x_j,x_{k+1}) X(x_{k+1}) \frac{dx_{k+1}}{x_i-x_{k+1}} = \frac{1}{x_i-x_j} ( K(x_j,x_a) )_{a=1}^k

\displaystyle - \frac{\cos(\pi(x_i-x_j))}{x_i-x_j} ( K(x_i,x_a) )_{a=1}^k

when {i \neq j}, and

\displaystyle \int_{\bf R} K(x_j,x_{k+1}) X(x_{k+1}) \frac{dx_{k+1}}{x_i-x_{k+1}}) = -( \partial_i K(x_i,x_a) )_{a=1}^k

when {i=j}. In the former case, we see from row operations that

\displaystyle \det( R_j( A_k, \int_{\bf R} K(x_j,x_{k+1}) X(x_{k+1}) \frac{dx_{k+1}}{x_i-x_{k+1}}) ) = \frac{1}{x_i-x_j} \det(A_k)

while from the Leibniz rule and symmetry (and (8)) we have

\displaystyle \partial_{x_i} \det(A_k) = 2 \det( R_i( A_k, ( \partial_i K(x_i,x_a) )_{a=1}^k ) );

putting all of these facts together, we obtain (14) as required.

Remark 1 Presumably, a similar exercise can be performed for the Airy process that governs the asymptotics of extreme eigenvalues of GUE (which include the famous Tracy-Widom law), but I did not attempt this.