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 \) .