The two-stream instability in a particle-in-cell code

Two beams of electrons move through each other at velocities $+v_0$ and $-v_0$, over a fixed background of ions that keeps the plasma neutral on average. A small ripple in density bunches both beams, the bunches make an electric field, and the field bunches them further. The ripple grows exponentially at a rate $\gamma$ that linear theory predicts, until electrons are trapped in the wave's potential wells and the growth stops.

A particle-in-cell (PIC) code follows this with no model beyond Newton's law and Gauss's law. Each step it spreads the particles' charge onto a grid, solves for the electric field on the grid, reads the field back at each particle, and moves the particles. Units: the plasma frequency ${\omega }_p=1$, the beam speed $v_0=1$.

For two cold beams at $\pm v_0$, each of density $n_0/2$, write $a=k{v}_0/{\omega }_p$. The dispersion relation

is the quartic ${\omega }^4-(2a^2+1){\omega }^2+a^4-a^2=0$, and its unstable root has

which is positive only for $a<1$. The cell below computes $\gamma$ from it.

The particle code

The box is one wavelength of the fastest mode, $L=2\pi /k$, with 32 grid cells and 2000 particles, 1000 per beam. The particles start evenly spaced, a quiet start, with their positions shifted by $10^{-3}\sin kx$ so that only the one mode is seeded. With random positions the instability would grow out of noise in every mode at once, and the fit below would need far more particles.

• Deposit. Each particle's charge is shared between the two nearest grid points in proportion to how close it is (cloud-in-cell). The ions add a uniform $+1$.
• Field. Gauss's law $dE/dx=\rho$, integrated along the grid by the trapezoid rule; the mean of $E$ is removed, since a periodic box has no net field.
• Push. Leapfrog: velocities live at half steps, positions at whole steps, and an electron's acceleration is $-E$.

Running it

The cell runs 250 steps to $t=25$. Every 12 steps it adds a frame to the film: a histogram of the particles in the $(x,v)$ plane. It also records the amplitude of the seeded mode at every step.

In the film the two beams are the bright lines at $v=\pm 1$. They ripple, then fold over each other into vortices: electrons caught in the wave's potential wells, going round at the bounce frequency. At that point the field stops growing.

The growth rate

While the mode is small it grows as $e^{\gamma t}$, so $\ln |E_1|$ is a straight line of slope $\gamma$. The cell fits that line by least squares between $t=8$ and $t=18$. Before $t=8$ the seed still carries the other three roots of the quartic, a decaying one and an oscillating pair, and their share falls only as $e^{-\gamma t}$. After $t=18$ the amplitude is above 0.1 and trapping has begun to slow the growth.

The fit gives $\gamma =0.35366$, $3\times 10^{-4}$ above theory. The test allows 1%, and here is why. Moving the window to $[4,12]$ or $[6,16]$ moves the fitted rate by up to 0.4%, so the choice of window alone is worth that much. Halving the time step, doubling the grid or quadrupling the particles moves it by at most 0.4% as well (posts/refs/two_stream_pic.py runs each case). A numpy implementation of the same algorithm gives the same rate to $10^{-9}$.

Saturation

The amplitude of the mode peaks at 0.568 at $t=23.8$ and then falls back: the trapped electrons now give energy back to the field as fast as they take it. From $t=0$ to the peak, the energy in the field grows by a factor of about $3\times 10^9$. The dashed line in the figure is $e^{\gamma t}$ with the theoretical rate, drawn through the fitted amplitude at $t=8$.