Illustration of Dyson Brownian Motion

The Dyson Brownian motion is defined as $$\begin{equation}\label{eq:20151124a} X_{t + dt} = X_t + H \sqrt{dt} \end{equation}$$ with $H$ a random matrix from the Gaussian Unitary Ensemble of $n \times n$ Hermitian matrices. It is then well-known that the dynamics of the eigenvalues $\lambda_i(t)$ of $X_t$ is described by the process $$\begin{equation}\label{eq:20151124b} d\lambda_i = \sum_{ j \neq i} \frac{1}{\lambda_i - \lambda_j} dt + dB_i \end{equation}$$ where $B_1, \ldots, B_n$ are independent Brownian processes. In this post I illustrate the process $\eqref{eq:20151124b}$ numerically.
I illustrate $\eqref{eq:20151124b}$ in the case $n = 5$. The Mathematica code I use is
Clear[RandomMatrixGUE] 
(* This code calculates a random matrix from the GUE *) 
RandomMatrixGUE[NN_] := Module[{A, B, H}, 
A = RandomVariate[NormalDistribution[], {NN, NN}]; 
B = RandomVariate[NormalDistribution[], {NN, NN}]; 
H = 1/2 (A + Transpose[A]) + I / 2 (B - Transpose[B]); 
H]

(* Simulates a path of the Dyson Brownian motion 
in: 
The dimension of the matrices is size x size. 
T is the maximum time to simulate. 
NSteps is in how many steps you simulate. 
out: a list of matrices X_t *)
SimulateOnePath[Size_, T_, NSteps_] := Module[{dt = T/NSteps , path}, 
path = Table[RandomMatrixGUE[Size] Sqrt[dt], {NSteps}]; 
path = Accumulate[path]; 
path = Prepend[path, ConstantArray[0, {Size, Size}]]; 
path]

(* plot all eigenvalues of the path X_t *) 
(* the code contains an ugly double loop *) 
PlotEigenValues[X_, dt_] := Module[{t = 0 , i, j, lambdapath, eig, NSteps = Length[X], n = Length[X[[1]]]}, 
lambdapath = ConstantArray[{}, n]; 
For[i = 1, i <= NSteps, i++, 
eig = Eigenvalues[X[[i]]]; 
eig = Sort[eig]; 
For[ j = 1, j <= n, j++, lambdapath[[j]] = Append[lambdapath[[j]], {t, eig[[j]]}]]; 
t = t + dt]; 
ListPlot[lambdapath, Joined -> True]] 

T = 1; 
NSteps = 1000; 
X = SimulateOnePath[5, T, NSteps]; 
PlotEigenValues[X, T/NSteps] 

This produces the following graph.

The eigenvalues $\lambda_1, \lambda_2, \ldots, \lambda_5$ as function of time

One can see in the above graph that the paths for $\lambda_1, \lambda_2, \ldots, \lambda_5$ never cross. This is because the drift $\frac{1}{\lambda_i - \lambda_j}$ becomes very large if the eigenvalues $\lambda_i$ and $\lambda_j$ become close. In the random matrix literature, the phenomenon is called eigenvalue repulsion. Some physicists say that the eigenvalues behave as fermions. In contrast, independent Brownian motions $B_i$ often cross, see graph below.

5 paths of the Brownian motion

Further reading

1. A proof of $\eqref{eq:20151124b}$ follows quite easily if one uses 1st and 2nd order perturbation formulas from quantum mechanics on the expansion $\lambda_i ( X_t + H \sqrt{dt} )$. Details can for example be found in this blog post by Terence Tao.
2. In the next post I perform a similar simulation for the Ornstein-Uhlenbeck process $$\begin{equation*} dX_t = -\alpha X_t dt + H \sqrt{dt} \end{equation*}$$ There the repulsive effect is even more striking because the process tries to move $X_t$ back to $0$ if $\alpha > 0$ .