Electrostatic Plasma Simulation: Electrons

APAM1601y
Columbia University

Introduction

References:
	Berkeley Plasma Simulation Group
	Prof. Charles "Ned" Birdsall
	Plasma Physics Via Computer Simulation, by Birdsall and Langdon, (Institute of Physics, 1991)

This notebook calculates a (very small) "particle-in-cell" (PIC) simulation of a one-dimensional plasma. Electrons and ions interact collectively by way of electrostatic forces. This is a very demanding computation. PIC simulations are usually performed on large supercomputers and require many hours (or days) of computation. Current research using these types of plasma codes are 3D, electromagnetic, and defined with complicated geometries.

Equations

The equations of motion for the ith particle are
	(∂_t) (Subscript[x, i]) = (Subscript[v, i])
	(∂_t) (Subscript[v, i]) = ((Subscript[q, i])/(Subscript[m, i])) E((Subscript[x, i]), t)

The equation for the electric field
	E(x, t) = – (∂_x) Φ(x, t)
	(∂_(x,x))Φ = – ρ(x, t) /(Subscript[ϵ, 0])
where ρ is the total charge density, including both the electrons and the ions. 
But, we can speed the calculations if we assume the ions are so heavy that they do not move. This restricts our problems to considering electron dynamics only.

Physics

A plasma is characterized by a length scale called the "Debye length", a time scale by the (inverse) "plasma frequency", and a velocity scale associated with the "thermal velocity".

In MKS units,
	(Subscript[ω, p]) = (Sqrt[((Subscript[q, e])^2)n/(Subscript[ϵ, 0])(Subscript[m, e])]), proportional to square-root of the density
	(Subscript[v, the]) = (Sqrt[(Subscript[kT, e])/(Subscript[m, e])]), proportional to the square-root of the temperature
	(Subscript[λ, D])= (Subscript[v, the])/(Subscript[ω, p]), the ratio of the thermal speed to the plasma frequency

In order for us to properly simulate a plasma, we need a "box" much larger than (Subscript[λ, D])and a simulation that lasts much longer than ((Subscript[ω, p])^(-1)). We also require the numerical time step, Δt, to be much shorter than ((Subscript[ω, p])^(-1)).

Numerics

For simplicity, we'll make the following assignments:
	• Δx = Δt = 1 (the spatial and temporal grid points are integers)
	• ((Subscript[q, e])/(Subscript[m, e])) = -1 for electrons
	• ((Subscript[q, i])/(Subscript[m, i])) = 1/25 for ions. (The ions are actually much heavier!)

With Δx = 1, the total length of our "box" of particles must be L >> 1.

With Δt = 1 and with (Subscript[q, e])/(Subscript[m, e]) = 1, then Δt (Subscript[ω, p])= (Sqrt[(Subscript[q, e])n/(Subscript[ϵ, 0])])<< 1. If we have 2000 electrons, then n = (10^3)/L ~ 1. Thus, we need to set the parameter (Subscript[q, e])/(Subscript[ϵ, 0]) << 1.

With (Subscript[v, the]) = (Subscript[λ, D])/(Subscript[ω, p]), then (Subscript[v, the]) ~ 1.

Basic Definitions: Heavy stationary ions

SetDirectory[NotebookDirectory[]]

boxLength = 500; (* Length of plasma and number of grid points *)

numElectrons = 4000; (* also number of ions *)

qOverϵ = 4.0*^-4; (* (Subscript[q, e])/(Subscript[ϵ, 0]) << 1 *)

electronTemp = 0.25; (* (Subscript[kT, e]); (Subscript[v, the]) = (Sqrt[(Subscript[kT, e])]) *)

Print["Average Density = ", nAvg = numElectrons/boxLength];
Print["Electron Thermal Velocity = ", vthe = Sqrt[electronTemp]];
Print["Plasma Frequency = ", ωp = Sqrt[(numElectrons/boxLength) qOverϵ]];
Print["Debye Length = ", λD = vthe/ωp];
Print["Number of particles per cell = ", nAvg];
Print["Number of Debye lengths in box = ", boxLength/λD];
Print["Number of time steps per plasma oscillation = ",1/ωp];

Electrons

Each particle is represented by a coordinate, {(Subscript[x, i]), (Subscript[v, i])}.

electrons=Table[{RandomReal[{0,boxLength}], RandomVariate[NormalDistribution[0.0,vthe]]},{numElectrons}];

Dimensions[electrons]

ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Electron Phase Space"]

Plot[PDF[NormalDistribution[0,vthe],x],{x,-2,2},Filling->Axis];
Histogram[electrons[[All,2]],Automatic,"PDF"];
pltElectDistribution = Show[%,%%, AxesLabel -> {"v","F(v)"}]

The electron charge density is evaluated on the grid points. We use "triangular weighting" to map the electron's position (in-between the grid points) to the charge density at each point.

rho0 = Table[0.0, {boxLength+1}]; (* including 0...boxLength *)
rhoT = rho0;

addρ[x_Real] := Module[{loc},
loc = Floor[x];
rhoT[[loc+1]] += (loc+1-x);
rhoT[[loc+2]] += (x - loc);]

getρ[e_List] := Module[{},
rhoT = rho0;
addρ /@ First[Transpose[e]];
rhoT[[1]] += rhoT[[boxLength+1]];
Drop[rhoT,-1]]

rhoElectron = getρ[electrons]; // Timing

Plus @@ rhoElectron (* How many electrons? *)

pr1=ListPlot[rhoElectron]

pr2=ListPlot[MovingAverage[rhoElectron,Floor[λD]],Joined->True,PlotStyle->{Red,Thickness[0.01]}]

Show[pr1,pr2,PlotLabel->"Electron Density",AxesLabel->{"x","ρ"},PlotRange->{0,2 nAvg}]

Notice large fluctuations in the electron density exists. Why? Answer: the amount of "noise" in our simulation depends upon the number of particles. On average, we have only a few particles for cell. The (random) statistical fluctuations scale as ±(Sqrt[n]), where n is the average number of particles per cell. For n = 4, we get ± 50% noise!

Ions

We take the ion density to be a constant and equal to the average electron density.

ionDensity = numElectrons/(boxLength) //N

Finding the Electric Field

pr1=ListPlot[ionDensity-rhoElectron]

pr2=ListPlot[MovingAverage[ionDensity-rhoElectron,Floor[λD]],Joined->True,PlotStyle->{Green,Thickness[0.01]}]

Show[pr1,pr2,PlotLabel->"Charge Density",AxesLabel->{"x","ρ"}]

kFourier = Table[2 Sin[2.0 π n/(2 boxLength)], {n, 0, boxLength-1}]; (* This represents the "effective" wavenumber used to invert Poisson's equation with the FFT. *)

kFourier2 = kFourier^2;

ListPlot[(1/(Drop[kFourier2,1])),PlotRange->{0,25},PlotLabel->"1/(k^2)",AxesLabel->{"k",""},Joined->True]

For a given electron and ion density, this computes the electric field using an FFT to invert Poisson's equation...

getE[rhoE_] := Module[{Φ},
Φ = {0.0} ~Join~ (Drop[Fourier[ionDensity - rhoE],1]/Drop[kFourier2,1]);
Φ = Re[InverseFourier[qOverϵ Φ]];
-(RotateLeft[Φ] - RotateRight[Φ])/2]

eField = getE[rhoElectron]; // Timing

ListPlot[eField,PlotLabel->"Electric Field",AxesLabel->{"x","E"}]

Of course, with a "noisy" charge density, we also get a "noisy" electric field!

Electron and Ion Step

In this section, we use the "Euler-Cromer" method to advance the position and velocity of a particle for a given electric field.

stepE[eField_List][{x_,v_}] := Module[{loc,locp, force, vNew},
loc = Floor[x];
locp = Mod[loc+1,boxLength]+1;
force = eField[[loc+1]] (loc+1 - x) + eField[[locp]] (x - loc);
(* advance velocity *)
vNew = v - force;
{Mod[x + vNew,boxLength], vNew}];

(* 1st Time is slower *) (ParallelMap[stepE[eField][#]& , electrons];) // Timing

(electrons1 = ParallelMap[stepE[eField][#]& , electrons];) // Timing

View of the change in the particles' phase-space positions...

ListPlot[electrons1-electrons,Frame->True]

Many Steps

Now, we can repeat these calculations for many times steps. Be warned: as we increase the number of particles, this becomes a very long computation...

numSteps = 60;

(* How many plasma periods? *) (numSteps  5)ωp / (2 π)

pInd1 = 0;ProgressIndicator[Dynamic[pInd1],{1,numSteps}]
Timing[Do[
Do[rhoElectron=getρ[electrons];eField=getE[rhoElectron];electrons=ParallelMap[stepE[eField][#1]&,electrons];,{5}];
pInd1 = i;
ePlot1[i] = GraphicsRow[{ListPlot[eField,PlotLabel->"Electric Field",AxesLabel->{"x","E"},PlotRange->{-0.01,0.01},Joined -> True],ListPlot[electrons,AxesLabel->{"x","v"},PlotRange->{{1,boxLength},{-1.5,1.5}},PlotLabel->"Electron Phase Space"]}],{i, numSteps}];]

AnimationVideo[Show[ePlot1[i]], {i, 1, numSteps, 1}]
Export["Efield-e1.mov",%,"QuickTime"]

pr1=ListPlot[rhoElectron];
pr2=ListPlot[MovingAverage[rhoElectron,Floor[λD]],Joined->True,PlotStyle->{Red,Thickness[0.01]}];
Show[pr1,pr2,PlotLabel->"Electron Density",AxesLabel->{"x","ρ"}]

The mean-squared average velocity of the electrons and ions...

(Plus @@ (Sqrt[#^2]& /@ Last[Transpose[electrons]]))/numElectrons

ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Electron Phase Space"]

Electron Beam and Electrons

What happens with a different initial condition?

In this section, we examine counter-drifting electrons.

beamFraction = 0.2;

electrons=Table[{RandomReal[{0,boxLength}], RandomVariate[NormalDistribution[0.0,vthe]]},{numElectrons-Floor[beamFraction numElectrons]}] ~Join~Table[{RandomReal[{0,boxLength}],3 vthe  (1+RandomReal[{-0.1,0.1}])},{Floor[beamFraction numElectrons]}];

Plot[PDF[NormalDistribution[0,vthe],x],{x,-2,2},Filling->Axis];
Histogram[electrons[[All,2]],Automatic,"PDF"];
Show[%,%%, AxesLabel -> {"v","F(v)"}]

ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Election Phase Space"]

Many Steps

numSteps = 150;

(* How many plasma periods? *) (numSteps  5)ωp/(2 π)

eFieldSave = {};

pInd2 = 0;ProgressIndicator[Dynamic[pInd2],{1,numSteps}]
Timing[Do[
Do[rhoElectron=getρ[electrons];
eField=getE[rhoElectron];
electrons=ParallelMap[stepE[eField][#1]&,electrons];
,{5}];
eFieldSave = Append[eFieldSave, eField];
pInd2 = i;ePlot2[i] = GraphicsRow[{ListPlot[eField,PlotLabel->"Electric Field",AxesLabel->{"x","E"},PlotRange->{-0.02,0.02}, Joined -> True],ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Electron Phase Space",PlotRange->{-1.5,2.0}]}],{i,numSteps}];]

AnimationVideo[Show[ePlot2[i]], {i, 1, numSteps, 1}]
Export["Efield-e2.mov",%,"QuickTime"]

ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Electron Phase Space"]

ListPlot[eField,PlotLabel->"Electric Field",AxesLabel->{"x","E"}, Joined -> True]

Dimensions[eFieldSave]

ListDensityPlot[Transpose[eFieldSave],AspectRatio -> 1/GoldenRatio, FrameLabel -> {"time", "Length"}]

Dimensions[Fourier[eFieldSave]]
Dimensions[Abs[Fourier[eFieldSave][[1;;numSteps/2,1;;boxLength/2]]]]
Max[Abs[Fourier[eFieldSave][[1;;numSteps/2,1;;boxLength/2]]]]

ListDensityPlot[Abs[Fourier[eFieldSave][[1;;Floor[numSteps/4],1;;Floor[boxLength/10]]]],AspectRatio -> 1/GoldenRatio, FrameLabel -> {"k", "ω"}, PlotRange -> {0,0.2}]

pr1=ListPlot[rhoElectron];
pr2=ListPlot[MovingAverage[rhoElectron,Floor[λD]],Joined->True,PlotStyle->{RGBColor[1,0,0],Thickness[0.01]}];
Show[pr1,pr2,PlotLabel->"Electron Density",AxesLabel->{"x","ρ"}]

Counter-Drifting Electrons

What happens with a different initial condition?

In this section, we examine counter-drifting electrons.

electrons=Table[{RandomReal[{0,boxLength}],(RandomInteger[{0,1}]-0.5) vthe 4 (1+RandomReal[{-0.1,0.1}])},{numElectrons}];

Plot[PDF[NormalDistribution[0,vthe],x],{x,-2,2},Filling->Axis];
Histogram[electrons[[All,2]],Automatic,"PDF"];
pltInitialElectronDist = Show[%,%%, AxesLabel -> {"v","F(v)"}]

ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Election Phase Space"]

Many Steps

numSteps = 150;

(* How many plasma periods? *) (numSteps  5)ωp / (2 π)

eFieldSave = {};

pInd3 = 0;ProgressIndicator[Dynamic[pInd3],{1,numSteps}]
Timing[Do[
Do[rhoElectron=getρ[electrons];
eField=getE[rhoElectron];electrons=ParallelMap[stepE[eField][#1]&,electrons];
,{5}];
eFieldSave = Append[eFieldSave, eField];
pInd3 = i;ePlot3[i] = GraphicsRow[{ListPlot[eField,PlotLabel->"Electric Field",AxesLabel->{"x","E"},PlotRange->{-0.04,0.04},Joined -> True],ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Electron Phase Space",PlotRange->{-2.0,2.0}]}],{i,numSteps}];]

AnimationVideo[Show[ePlot3[i]], {i, 1, numSteps, 1}]
Export["Efield-e3.mov",%,"QuickTime"]

ListPlot[electrons,AxesLabel->{"x","v"},PlotLabel->"Electron Phase Space"]

Plot[PDF[NormalDistribution[0,vthe],x],{x,-2,2},Filling->Axis];
Histogram[electrons[[All,2]],Automatic,"PDF", ChartStyle -> Green];
Show[%,pltInitialElectronDist, AxesLabel -> {"v","F(v)"}]

ListPlot[eField,PlotLabel->"Electric Field",AxesLabel->{"x","E"}]

ListDensityPlot[Transpose[eFieldSave],AspectRatio -> 1/GoldenRatio, FrameLabel -> {"time", "Length"}]

Max[Abs[Fourier[eFieldSave][[1;;numSteps/2,1;;boxLength/2]]]]

ListDensityPlot[Abs[Fourier[eFieldSave][[1;;Floor[numSteps/4],1;;Floor[boxLength/10]]]],AspectRatio -> 1/GoldenRatio, FrameLabel -> {"k", "ω"}, PlotRange -> {0,0.5}]

pr1=ListPlot[rhoElectron];
pr2=ListPlot[MovingAverage[rhoElectron,Floor[λD]],Joined->True,PlotStyle->{RGBColor[1,0,0],Thickness[0.01]}];
Show[pr1,pr2,PlotLabel->"Electron Density",AxesLabel->{"x","ρ"}]

Summary

A one-dimensional "particle-in-cell" simulation of a plasma was defined using Mathematica. PIC simulations represent a direct model of a plasma allowing computation of dynamical and statistical properties of plasmas. We examined three particularly easy examples determined by the initial conditions of the distributions of ions and electrons. (Many more examples are possible–even with this simple model!) In one example, we looked at "noise" in a stable plasma. In the second example, we looked at the evolution of counter-streaming electrons that excite the so-called "two-stream" instability. The initial kinetic energy of the electrons excites growing electrostatic oscillations. Particle energy is converted into electric-field energy. In the third example, we looked at a "clump" of plasma and observed "ambipolar" electron confinement and the plasma Debye length.