{"id":15461,"date":"2021-10-23T14:08:51","date_gmt":"2021-10-23T12:08:51","guid":{"rendered":"https:\/\/djalil.chafai.net\/blog\/?p=15461"},"modified":"2021-10-30T08:45:00","modified_gmt":"2021-10-30T06:45:00","slug":"dynamics-of-a-planar-coulomb-gas","status":"publish","type":"post","link":"https:\/\/djalil.chafai.net\/blog\/2021\/10\/23\/dynamics-of-a-planar-coulomb-gas\/","title":{"rendered":"Dynamics of a planar Coulomb gas"},"content":{"rendered":"

This short post is about a joint work with Fran\u00e7ois Bolley<\/a> and Joaqu\u00edn Fontbona<\/a> on a Dyson Ornstein-Uhlenbeck process, published a few years ago. The numerics and graphics are new.<\/p>\n

\"\"<\/a><\/p>\n

\"\"<\/a><\/p>\n

\"\"<\/a><\/p>\n

\"\"<\/a><\/p>\n

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
\n$$
\n\\mathbb{d}X^{i,N}_t
\n=\\sqrt{\\frac{2}{N}}\\mathbb{d} B^{i,N}_t
\n-2X^{i,N}_t\\,\\mathbb{d} t
\n-\\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,
\n\\quad 1\\leq i\\leq N,
\n$$
\nwhere $\\beta\\geq0$ is a parameter and $B$ is a standard Brownian motion on $(\\mathbb{R}^2)^N=\\mathbb{C}^N$.<\/p>\n

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

Barycenter (first moment) of the cloud of particles.<\/strong> 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
\n$$
\n\\mathrm{d}m_t=\\sqrt{\\frac{2}{N^2}}\\mathrm{d}W_t-2m_t\\mathrm{d}t
\n$$
\nwhere $W$ is a standard Brownian motion on $\\mathbb{R}^2=\\mathbb{C}$. In particular, for $t\\geq0$,
\n$$
\nm_t\\sim\\mathcal{N}\\Bigr(m_0\\mathrm{e}^{-2t},\\frac{1-\\mathrm{e}^{-4t}}{2N^2}I_2\\Bigr).
\n$$\n<\/p>\n

This follows from the It\u00f4 formula and the fact that the real and imaginary parts of the first moment are eigenfunctions of the Dyson Ornstein-Uhlenbeck operator.<\/p>\n

Second moment of the cloud of particles.<\/strong> 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
\n$$
\n\\mathrm{d}R_t = \\sqrt{\\frac{8}{N^2}R_t}\\mathrm{d}W_t
\n +4\\left[\\frac{1}{N}+ \\frac{\\beta}{4}\\frac{N-1}{N}-R_t \\right]\\mathrm{d}t,
\n$$
\nwhere ${(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$,
\n$$
\n\\begin{align*}
\n\\mathrm{Wasserstein}_1(\\mathrm{Law}(R_t),\\Gamma_N)
\n&\\leq \\mathrm{e}^{-4t}\\mathrm{Wasserstein}_1(\\mathrm{Law}(R_0),\\Gamma_N)\\\\
\n\\mathbb{E}[R_t\\mid R_0]
\n&= R_0\\mathrm{e}^{-4t}+\\Bigr(\\frac{1}{N}+\\frac{\\beta}{4}\\frac{N-1}{N}\\Bigr)(1-\\mathrm{e}^{-4t}).
\n\\end{align*}
\n$$
\nNote that $\\frac{\\beta}{4}$ is the second moment of the uniform distribution on the disc of radius $\\sqrt{\\frac{\\beta}{2}}$.\n<\/p>\n

This follows from It\u00f4's formula, L\u00e9vy's characterization of BM, and the fact that the second moment is, up to a constant, an eigenfunction of the Dyson Ornstein-Uhlenbeck operator.<\/p>\n

Open problem.<\/strong> Compute the spectral gap of the dynamics, or find how it depends over $N$.<\/p>\n

Julia code.<\/strong><\/p>\n

\r\nusing Plots\r\n\r\n# Define the 2D Dyson Ornstein-Uhlenbeck process\r\nBase.@kwdef mutable struct DOU2D\r\n    n::Int = 2         # number of particles\r\n    \u03b2::Float64 = 2.    # repulsion parameter  \r\n    T::Float64 = 10.   # terminal time\r\n    dt::Float64 = 1E-3 # time increment\r\n    r::Float64  = 1.   # initial conditions uniform on segment [-r,r]+3*im\r\n    m = floor(Int, T \/ dt) # number of times\r\n    x::Array{Complex{Float64},2} = zeros(Complex{Float64},m,n)\r\nend # end struct\r\n\r\nfunction compute!(X::DOU2D)\r\n    # initial condition on grid on centered circle of radius r\r\n    #Theta = range(1, X.n, length = X.n) * 2 * pi \/ X.n\r\n    #X.x[1,:] = X.r * (cos.(Theta) + im * sin.(Theta)) # unif grid on circle\r\n    # initial condition on grid on segment [-r,r] - 3*im\r\n    X.x[1,:] = range(-X.r, X.r, length = X.n) .- 3 * im\r\n    # trajectories\r\n    for i = 2:X.m\r\n     dB = (randn(1,X.n) + im * randn(1,X.n))\/sqrt(2)\r\n     X.x[i,:] = copy(X.x[i-1,:])\r\n     for j in 1:X.n # particles\r\n         X.x[i,j] += sqrt(2\/X.n) * dB[j] * sqrt(X.dt)\r\n         X.x[i,j] += -2 * X.x[i-1,j] * X.dt\r\n         for k in 1:X.n\r\n             if (j == k) continue end\r\n             X.x[i,j] += X.\u03b2 * X.dt \/ (X.n * (conj(X.x[i-1,j] - X.x[i-1,k])))\r\n         end # for \r\n     end # for\r\n    end # for\r\nend # function\r\n\r\n# process and graphics\r\ndou = DOU2D(n = 66, r = 3., T = 5., dt = 1E-4)\r\ncompute!(dou)\r\n#\r\npdou2d = plot()\r\nfor j in 1:dou.n # particles\r\n    plot!(pdou2d, real(dou.x[:,j]), imag(dou.x[:,j]), aspect_ratio =:equal, legend = false)\r\nend # for\r\nsavefig(pdou2d,\"dou2d.png\")\r\nsavefig(pdou2d,\"dou2d.svg\")\r\n#\r\nr = dou.r + 1\r\nanim = @animate for i = 1:100:dou.m\r\n    scatter(real(dou.x[i,:]), imag(dou.x[i,:]),\r\n            xlims = (-r,r), ylims = (-r,r),\r\n            aspect_ratio =:equal, legend = false)\r\nend\r\ngif(anim, \"dou2danim.gif\", fps = 10)\r\n<\/pre>\n
Further reading.<\/strong><\/p>\n