using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations
gr();
We study the dynamics of many electric charges in a potential well. We restrict our attention to 1D: picture an infinitely long wire with carges on it. We will see that as the number of charges becomes large, we can determine the limiting distribution using an additive Riemann–Hilbert problem, but one which the interval its posed on is solved for.
Consider a point charge in a well V(x)=x2/2, initially located at λ0. The dynamics of the point charge are governed by dλdt=−V′(λ)=−λ
V = x -> x^2/2
Vp = x -> x
λ⁰ = 2.3 # initial location
prob = ODEProblem((λ,_,t) -> -Vp(λ), λ⁰, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);
@gif for t=0.0:0.05:7.0
plot(V, range(-5, stop=5, length=100); label="potential", title="t = $t")
scatter!([λ(t)] ,[0.0]; label="charge")
end
In the limit the charge reaches an equilibrium: it no longer varies in time. I.e., it reaches a point where dλdt=0, which is equivalent to solving 0=−V′(λ)=−λ
Suppose there are now two charges, λ1 and λ2. The effect on the first charge λ1 is to repulse away from λ2 via via: dλ1dt=1λ1−λ2
λ⁰ = [1.2,2.3] # initial location
prob = ODEProblem((λ,_,t) -> [1/(λ[1] - λ[2]), 1/(λ[2] - λ[1])], λ⁰, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);
@gif for t=0.0:0.05:7.0
scatter(λ(t) ,zeros(2); label="charges", xlims=(-5,5), title="t = $t")
end
Adding in a potential well and we get an equilbrium again: dλ1dt=1λ1−λ2−V′(λ1)dλ2dt=1λ2−λ1−V′(λ2)
λ₀ = [1.2,2.3] # initial location
prob = ODEProblem((λ,_,t) -> [1/(λ[1] - λ[2]) - Vp(λ[1]), 1/(λ[2] - λ[1]) - Vp(λ[2])], λ₀, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);
@gif for t=0.0:0.05:7.0
plot(V, range(-5, stop=5, length=100); label="potential", title="t = $t")
scatter!(λ(t) ,zeros(2); label="charges", xlims=(-5,5), title="t = $t")
end
The limiting distribution is given by 0=1λ1−λ2−V′(λ1)0=1λ2−λ1−V′(λ2)
Each charge repulses every other charge, so we end up needing to sum over them all: dλkdt=N∑j=1j≠k1λk−λj−V′(λk)
N = 100
λ⁰ = randn(N) # initial location
prob = ODEProblem((λ,_,t) -> [sum(1 ./ (λ[k] .- λ[[1:k-1;k+1:end]])) - Vp(λ[k]) for k=1:N], λ⁰, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);
@gif for t=0.0:0.05:7.0
plot(V, range(-5, stop=5, length=100); label="potential", title="t = $t")
scatter!(λ(t) ,zeros(N); label="charges", xlims=(-5,5), title="t = $t")
end
As the number of charges becomes large, they spread off to infinity. In the case of V(x)=x2, we can renormalize by dividing by N so they stay bounded: μk=λk√N
@gif for t=0.0:0.05:7.0
scatter(λ(t)/sqrt(N) ,zeros(N); label="charges", xlims=(-2,2), title="t = $t")
end
This begs questions: why does it balance out at ±√2? Why does it have a nice histogram precisely like √2−x2π:
histogram(λ(10.0)/sqrt(N); nbins=20, normalize=true, label="histogram of charges")
plot!(x -> sqrt(2-x^2)/(π), range(eps()-sqrt(2.0); stop=sqrt(2)-eps(), length=100), label="semicircle")
Plugging in λk=√Nμk, we get a dynamical system for μk: dμkdt=1NN∑j=1j≠k1μk−μj−μk
It is convenient to represent the point charges by Dirac delta functions: wN(x)=1NN∑k=1δμk(x)
Formally (see a more detailed explanation below), wN(x) tends to a continuous limit as N→∞, which we have guessed from the histogram to be w(x)=√2−x2π for −√2<x<√2. We expect this limit to satisfy the same equation as wn, that is Hw(x)=−xπ
Indeed:
x = Fun(-sqrt(2) .. sqrt(2))
w = sqrt(2-x^2)/π
norm(hilbert(w) + x/π)
UndefVarError: norm not defined Stacktrace: [1] top-level scope at In[16]:3
Why is it [−√2,√2]?
We thus want to choose the interval [a,b] so that there exists a w(x) satisfying
As we saw last lecture, there exists a bounded solution to H[−b,b]u=−x/π, namely u(x)=√b2−x2π. The choice b=√2 ensures that ∫b−bu(x)dx=1, hence u(x)=w(x).
This is beyond the scope of the course, but the convergence of wN(x) to w(x) is known as weak-* convergence. A simple version of this is that ∫dcwN(x)dx→∫dcw(x)dx
a = -0.1; b= 0.3;
w = Fun(x -> sqrt(2-x^2)/π, a .. b)
sum(w) # integral of w(x) between a and b
0.1790059168895313
length(filter(λ -> a ≤ λ ≤ b, λ(10.0)/sqrt(N)))/N # integral of w_N(x) between a and b
0.18
Another varient of describing weak-* convergence is that the Cauchy transforms converge, that is, for z on any contour γ surrounding a and b (now the support of w) we have ∫bawN(x)x−zdx→∫baw(x)x−zdx
x = Fun(-sqrt(2) .. sqrt(2))
w = sqrt(2-x^2)/π
z = 0.1+0.2im
cauchy(w, z)
0.1949409042727309 + 0.013681505291622225im
(sum(1 ./((λ(10.0)/sqrt(N)) .- z))/N)/(2π*im)
0.19542123392652283 + 0.013724448552242337im