Press "Enter" to skip to content

Dynamics of a planar Coulomb gas

This short post is about a joint work with François Bolley and Joaquín Fontbona on a Dyson Ornstein-Uhlenbeck process, published a few years ago. The numerics and graphics are new.

Above, two numerical experiments with $N=66$ particles and initial conditions on a uniform grid on a circle, and on a segment. Data is displayed as moving particles as well as trajectories. The one started from a segment is run ten times slower than the one started from a circle. These interacting particles on $\mathbb{R}^2=\mathbb{C}$ follow the Dyson Ornstein-Uhlenbeck dynamics
=\sqrt{\frac{2}{N}}\mathbb{d} B^{i,N}_t
-2X^{i,N}_t\,\mathbb{d} t
-\frac{\beta}{N}\sum_{j\neq i}\frac{X^{j,N}_t-X^{i,N}_t}{|X^{i,N}_t-X^{j,N}_t|^2}\mathbb{d}t,
\quad 1\leq i\leq N,
where $\beta\geq0$ is a parameter and $B$ is a standard Brownian motion on $(\mathbb{R}^2)^N=\mathbb{C}^N$.

The long time behavior of the particle system is
where $X^N_\infty$ follows the $\beta$ Ginibre planar Coulomb gas
\prod_{1\leq i < j\leq N}|x_i-x_j|^\beta \mathrm{d}x_1\cdots\mathrm{d}x_N, $$ while the high dimensional limit at equilibrium is the uniform law on a disc $$ \frac{1}{N}\sum_{i=1}^N\delta_{X^{i,N}_\infty} \overset{\mathrm{weak}}{\underset{N\to\infty}{\longrightarrow}} \frac{\mathbf{1}_{z\in\mathbb{C}:|z|\leq\sqrt{\frac{\beta}{2}}}}{\pi}. $$

Barycenter (first moment) of the cloud of particles. The stochastic process ${(m_t)}_{t\geq0}={(\frac{1}{N}\sum_{i=1}^nX^{i,N}_t)}_{t\geq0}$ is an Orsntein-Uhlenbeck process solving
where $W$ is a standard Brownian motion on $\mathbb{R}^2=\mathbb{C}$. In particular, for $t\geq0$,

This follows from the Itô formula and the fact that the real and imaginary parts of the first moment are eigenfunctions of the Dyson Ornstein-Uhlenbeck operator.

Second moment of the cloud of particles. The process ${(R_t)}_{t\geq0}={(\frac{1}{N}\sum_{i=1}^n|X^{i,N}_t|^2)}_{t\geq0}$ is a Cox-Ingersoll-Ross process ${(R_t)}_{t\geq0}$ solving the stochastic differential equation
\mathrm{d}R_t = \sqrt{\frac{8}{N^2}R_t}\mathrm{d}W_t
+4\left[\frac{1}{N}+ \frac{\beta}{4}\frac{N-1}{N}-R_t \right]\mathrm{d}t,
where ${(W_t)}_{t\geq0}$ is a real standard Brownian motion. In particular, its invariant distribution is the Gamma law $\Gamma_N=\mathrm{Gamma}(a,\lambda)$ with shape parameter $a=N+\frac{\beta}{4}N(N-1)$ and scale parameter $\lambda=N^2$. Moreover, for any initial condition $R_0$ and $t \geq 0$,
&\leq \mathrm{e}^{-4t}\mathrm{Wasserstein}_1(\mathrm{Law}(R_0),\Gamma_N)\\
\mathbb{E}[R_t\mid R_0]
&= R_0\mathrm{e}^{-4t}+\Bigr(\frac{1}{N}+\frac{\beta}{4}\frac{N-1}{N}\Bigr)(1-\mathrm{e}^{-4t}).
Note that $\frac{\beta}{4}$ is the second moment of the uniform distribution on the disc of radius $\sqrt{\frac{\beta}{2}}$.

This follows from Itô’s formula, Lévy’s characterization of BM, and the fact that the second moment is, up to a constant, an eigenfunction of the Dyson Ornstein-Uhlenbeck operator.

Open problem. Compute the spectral gap of the dynamics, or find how it depends over $N$.

Julia code.

using Plots

# Define the 2D Dyson Ornstein-Uhlenbeck process
Base.@kwdef mutable struct DOU2D
    n::Int = 2         # number of particles
    β::Float64 = 2.    # repulsion parameter  
    T::Float64 = 10.   # terminal time
    dt::Float64 = 1E-3 # time increment
    r::Float64  = 1.   # initial conditions uniform on segment [-r,r]+3*im
    m = floor(Int, T / dt) # number of times
    x::Array{Complex{Float64},2} = zeros(Complex{Float64},m,n)
end # end struct

function compute!(X::DOU2D)
    # initial condition on grid on centered circle of radius r
    #Theta = range(1, X.n, length = X.n) * 2 * pi / X.n
    #X.x[1,:] = X.r * (cos.(Theta) + im * sin.(Theta)) # unif grid on circle
    # initial condition on grid on segment [-r,r] - 3*im
    X.x[1,:] = range(-X.r, X.r, length = X.n) .- 3 * im
    # trajectories
    for i = 2:X.m
     dB = (randn(1,X.n) + im * randn(1,X.n))/sqrt(2)
     X.x[i,:] = copy(X.x[i-1,:])
     for j in 1:X.n # particles
         X.x[i,j] += sqrt(2/X.n) * dB[j] * sqrt(X.dt)
         X.x[i,j] += -2 * X.x[i-1,j] * X.dt
         for k in 1:X.n
             if (j == k) continue end
             X.x[i,j] += X.β * X.dt / (X.n * (conj(X.x[i-1,j] - X.x[i-1,k])))
         end # for 
     end # for
    end # for
end # function

# process and graphics
dou = DOU2D(n = 66, r = 3., T = 5., dt = 1E-4)
pdou2d = plot()
for j in 1:dou.n # particles
    plot!(pdou2d, real(dou.x[:,j]), imag(dou.x[:,j]), aspect_ratio =:equal, legend = false)
end # for
r = dou.r + 1
anim = @animate for i = 1:100:dou.m
    scatter(real(dou.x[i,:]), imag(dou.x[i,:]),
            xlims = (-r,r), ylims = (-r,r),
            aspect_ratio =:equal, legend = false)
gif(anim, "dou2danim.gif", fps = 10)

Further reading.

    Leave a Reply

    Your email address will not be published.

    This site uses Akismet to reduce spam. Learn how your comment data is processed.

    Syntax · Style · .