Notebooks · posts/two-stream-instability-in-a-particle-code

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

Adam Zettel · September 25, 2026 · written with an agent

Two cold electron beams passing through each other are unstable. A particle-in-cell code with 2000 particles measures the growth rate of the fastest mode as 0.35366 $ω_{p}$ ; cold-beam linear theory gives $1/ (22) = 0.35355$ . The phase-space film shows the beams rolling up into vortices of trapped electrons, which is where the growth stops.

γ_{m a x} = \frac{ω _{p}}{2 2} at k v_{0} = 3/8 ω_{p}

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 $γ$ 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 $ω_{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} / ω_{p}$ . The dispersion relation

1 = \frac{1}{2} [\frac{1}{( ω - a ) ^{2}} + \frac{1}{( ω + a ) ^{2}}]

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

γ^{2} = \frac{1}{2} [8 a^{2} + 1 - (2 a^{2} + 1)],

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

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 37–59

PI = 3.141592653589793
# the fastest-growing wavenumber of two equal cold beams: k v0 = sqrt(3/8) omega_p
KMODE = sqrt(3.0 / 8.0)

#| The cold-beam growth rate γ at a = k·v0/ωp, from the quartic above;
#| zero for a ≥ 1.
def cold_growth(a) {
    g2 = 0.5 * (sqrt(8.0 * a * a + 1.0) - (2.0 * a * a + 1.0))
    if g2 < 0.0 {
        return 0.0
    }
    return sqrt(g2)
}
gamma_th = cold_growth(KMODE)
show gamma_th

test "the cold-beam growth rate at k v0 = sqrt(3/8) is 1 / (2 sqrt 2) = 0.353553" {
    assert_near(gamma_th, 1.0 / (2.0 * sqrt(2.0)), 1.0e-15)
}
test "the root of the quartic at a = 0.48 grows at 0.335693, and at a = 1.2 it is stable (numpy.roots)" {
    assert_near(cold_growth(0.48), 0.3356926515941623, 1.0e-13)
    assert_near(cold_growth(1.2), 0.0, 0.0)
}

The particle code

The box is one wavelength of the fastest mode, $L = 2 π / k$ , with 32 grid cells and 2000 particles, 1000 per beam. The particles start evenly spaced, a quiet start, with their positions shifted by $1 0^{- 3} sin k x$ 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 $d E / d x = ρ$ , 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$ .

[2]lines 78–176

NPART = 2000
NGRID = 32
DT = 0.1
NSTEP = 250
# so the run ends at t = 25
T_END = NSTEP * DT
SEED = 1.0e-3
LBOX = 2.0 * PI / KMODE
DX = LBOX / NGRID

#| Wrap a position into the periodic box [0, LBOX).
def wrap(x) {
    y = x
    if y < 0.0 {
        y = y + LBOX
    }
    if y >= LBOX {
        y = y - LBOX
    }
    return y
}

#| Charge density on the grid: +1 from the ions, minus each electron shared
#| between its two nearest grid points (cloud-in-cell).
def deposit(xp: Arr, rho: Arr) {
    ng = len(rho)
    for j to ng {
        rho[j] = 1.0
    }
    w = ng / len(xp)
    for i to len(xp) {
        s = xp[i] / DX
        j0 = floor(s)
        f = s - j0
        j1 = j0 + 1.0
        if j1 >= ng {
            j1 = j1 - ng
        }
        rho[j0] = rho[j0] - w * (1.0 - f)
        rho[j1] = rho[j1] - w * f
    }
}

#| E from dE/dx = ρ by the trapezoid rule, with its mean removed.
def solve_field(rho: Arr, ef: Arr) {
    ng = len(rho)
    ef[0] = 0.0
    for j to ng - 1.0 {
        ef[j + 1] = ef[j] + 0.5 * (rho[j] + rho[j + 1]) * DX
    }
    mean = sum(ef) / ng
    for j to ng {
        ef[j] = ef[j] - mean
    }
}

#| The field at position x, interpolated from its two nearest grid points.
def gather(x, ef: Arr) {
    ng = len(ef)
    s = x / DX
    j0 = floor(s)
    f = s - j0
    j1 = j0 + 1.0
    if j1 >= ng {
        j1 = j1 - ng
    }
    return ef[j0] * (1.0 - f) + ef[j1] * f
}

#| The amplitude of the box's fundamental in E: (2/ng)·|∑ⱼ ef[j]·exp(−i·KMODE·j·DX)|.
def mode_amplitude(ef: Arr) {
    ng = len(ef)
    c = 0.0
    s = 0.0
    for j to ng {
        c = c + ef[j] * cos(KMODE * j * DX)
        s = s + ef[j] * sin(KMODE * j * DX)
    }
    return 2.0 / ng * sqrt(c * c + s * s)
}

xp = zeros(NPART)
vp = zeros(NPART)
HALF = NPART / 2.0
for i to HALF {
    x0 = (i + 0.5) * LBOX / HALF
    xp[i] = wrap(x0 + SEED * sin(KMODE * x0))
    vp[i] = 1.0
    xp[HALF + i] = wrap(x0 - SEED * sin(KMODE * x0))
    vp[HALF + i] = -1.0
}
rho = zeros(NGRID)
ef = zeros(NGRID)
deposit(xp, rho)
solve_field(rho, ef)
show sum(rho) / NGRID
test "the box is neutral: the mean charge density is zero" {
    assert_near(sum(rho) / NGRID, 0.0, 1.0e-12)
}

NPART = 2000
NGRID = 32
DT = 0.1
NSTEP = 250
# so the run ends at t = 25
T_END = NSTEP * DT
SEED = 1.0e-3
LBOX = 2.0 * PI / KMODE
DX = LBOX / NGRID

#| Wrap a position into the periodic box [0, LBOX).
def wrap(x) {
    y = x
    if y < 0.0 {
        y = y + LBOX
    }
    if y >= LBOX {
        y = y - LBOX
    }
    return y
}

#| Charge density on the grid: +1 from the ions, minus each electron shared
#| between its two nearest grid points (cloud-in-cell).
def deposit(xp: Arr, rho: Arr) {
    ng = len(rho)
    for j to ng {
        rho[j] = 1.0
    }
    w = ng / len(xp)
    for i to len(xp) {
        s = xp[i] / DX
        j0 = floor(s)
        f = s - j0
        j1 = j0 + 1.0
        if j1 >= ng {
            j1 = j1 - ng
        }
        rho[j0] = rho[j0] - w * (1.0 - f)
        rho[j1] = rho[j1] - w * f
    }
}

#| E from dE/dx = ρ by the trapezoid rule, with its mean removed.
def solve_field(rho: Arr, ef: Arr) {
    ng = len(rho)
    ef[0] = 0.0
    for j to ng - 1.0 {
        ef[j + 1] = ef[j] + 0.5 * (rho[j] + rho[j + 1]) * DX
    }
    mean = sum(ef) / ng
    for j to ng {
        ef[j] = ef[j] - mean
    }
}

#| The field at position x, interpolated from its two nearest grid points.
def gather(x, ef: Arr) {
    ng = len(ef)
    s = x / DX
    j0 = floor(s)
    f = s - j0
    j1 = j0 + 1.0
    if j1 >= ng {
        j1 = j1 - ng
    }
    return ef[j0] * (1.0 - f) + ef[j1] * f
}

#| The amplitude of the box's fundamental in E: (2/ng)·|∑ⱼ ef[j]·exp(−i·KMODE·j·DX)|.
def mode_amplitude(ef: Arr) {
    ng = len(ef)
    c = 0.0
    s = 0.0
    for j to ng {
        c = c + ef[j] * cos(KMODE * j * DX)
        s = s + ef[j] * sin(KMODE * j * DX)
    }
    return 2.0 / ng * sqrt(c * c + s * s)
}

xp = zeros(NPART)
vp = zeros(NPART)
HALF = NPART / 2.0
for i to HALF {
    x0 = (i + 0.5) * LBOX / HALF
    xp[i] = wrap(x0 + SEED * sin(KMODE * x0))
    vp[i] = 1.0
    xp[HALF + i] = wrap(x0 - SEED * sin(KMODE * x0))
    vp[HALF + i] = -1.0
}
rho = zeros(NGRID)
ef = zeros(NGRID)
deposit(xp, rho)
solve_field(rho, ef)
show sum(rho) / NGRID
test "the box is neutral: the mean charge density is zero" {
    assert_near(sum(rho) / NGRID, 0.0, 1.0e-12)
}

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.

[3]lines 190–236

# the film: a phase-space histogram of 48 x-bins by 32 v-bins over |v| < VMAX
NXB = 48
NVB = 32
VMAX = 2.0
EVERY = 12
hist = zeros(NVB, NXB)
#| Count the particles in each (v, x) cell of the histogram; row 0 is v = +VMAX.
def film(xp: Arr, vp: Arr, h: Arr) {
    for r to NVB {
        for c to NXB {
            h[r, c] = 0.0
        }
    }
    for i to len(xp) {
        c = floor(xp[i] / LBOX * NXB)
        r = floor((VMAX - vp[i]) / (2.0 * VMAX) * NVB)
        if r >= 0.0 {
            if r < NVB {
                h[r, c] = h[r, c] + 1.0
            }
        }
    }
}

times = zeros(NSTEP)
amp = zeros(NSTEP)
energy = zeros(NSTEP)
# the leapfrog stagger: take the velocities back half a step
for i to NPART {
    vp[i] = vp[i] + 0.5 * DT * gather(xp[i], ef)
}
figure("Phase space of the two beams (press play)", "x", "v", 0.0, LBOX, -VMAX, VMAX)
for n to NSTEP {
    deposit(xp, rho)
    solve_field(rho, ef)
    times[n] = n * DT
    amp[n] = mode_amplitude(ef)
    energy[n] = 0.5 * DX * sum(ef * ef)
    if mod(n, EVERY) == 0.0 {
        film(xp, vp, hist)
        plot2d(hist, "electrons per cell")
    }
    for i to NPART {
        vp[i] = vp[i] - DT * gather(xp[i], ef)
        xp[i] = wrap(xp[i] + DT * vp[i])
    }
}

# the film: a phase-space histogram of 48 x-bins by 32 v-bins over |v| < VMAX
NXB = 48
NVB = 32
VMAX = 2.0
EVERY = 12
hist = zeros(NVB, NXB)
#| Count the particles in each (v, x) cell of the histogram; row 0 is v = +VMAX.
def film(xp: Arr, vp: Arr, h: Arr) {
    for r to NVB {
        for c to NXB {
            h[r, c] = 0.0
        }
    }
    for i to len(xp) {
        c = floor(xp[i] / LBOX * NXB)
        r = floor((VMAX - vp[i]) / (2.0 * VMAX) * NVB)
        if r >= 0.0 {
            if r < NVB {
                h[r, c] = h[r, c] + 1.0
            }
        }
    }
}

times = zeros(NSTEP)
amp = zeros(NSTEP)
energy = zeros(NSTEP)
# the leapfrog stagger: take the velocities back half a step
for i to NPART {
    vp[i] = vp[i] + 0.5 * DT * gather(xp[i], ef)
}
figure("Phase space of the two beams (press play)", "x", "v", 0.0, LBOX, -VMAX, VMAX)
for n to NSTEP {
    deposit(xp, rho)
    solve_field(rho, ef)
    times[n] = n * DT
    amp[n] = mode_amplitude(ef)
    energy[n] = 0.5 * DX * sum(ef * ef)
    if mod(n, EVERY) == 0.0 {
        film(xp, vp, hist)
        plot2d(hist, "electrons per cell")
    }
    for i to NPART {
        vp[i] = vp[i] - DT * gather(xp[i], ef)
        xp[i] = wrap(xp[i] + DT * vp[i])
    }
}

The growth rate

While the mode is small it grows as $e^{γ t}$ , so $ln ∣ E_{1} ∣$ is a straight line of slope $γ$ . 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^{- γ t}$ . After $t = 18$ the amplitude is above 0.1 and trapping has begun to slow the growth.

The fit gives $γ = 0.35366$ , $3 \times 1 0^{- 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 $1 0^{- 9}$ .

[4]lines 255–295

#| The least-squares slope of ln(y) against t over the samples with t0 ≤ t ≤ t1.
def log_slope(t: Arr, y: Arr, t0, t1) {
    n = 0.0
    st = 0.0
    sy = 0.0
    stt = 0.0
    sty = 0.0
    for i to len(t) {
        if t[i] >= t0 - 1.0e-9 {
            if t[i] <= t1 + 1.0e-9 {
                ly = log(y[i])
                n = n + 1.0
                st = st + t[i]
                sy = sy + ly
                stt = stt + t[i] * t[i]
                sty = sty + t[i] * ly
            }
        }
    }
    return (n * sty - st * sy) / (n * stt - st * st)
}
T_FIT0 = 8.0
T_FIT1 = 18.0
gamma_fit = log_slope(times, amp, T_FIT0, T_FIT1)
gamma_early = log_slope(times, amp, 4.0, 12.0)
gamma_mid = log_slope(times, amp, 6.0, 16.0)
show gamma_fit
show gamma_fit / gamma_th - 1.0
show gamma_early
show gamma_mid

test "the fitted growth rate is the cold-beam rate 0.353553 to within 1%" {
    assert_near(gamma_fit / gamma_th, 1.0, 0.01)
}
test "the same algorithm in numpy (posts/refs/two_stream_pic.py) fits 0.353661097" {
    assert_near(gamma_fit, 0.353661097275761, 1.0e-9)
}
test "fits over [4, 12] and [6, 16] stay within 0.4% of theory" {
    assert_near(gamma_early / gamma_th, 1.0, 0.004)
    assert_near(gamma_mid / gamma_th, 1.0, 0.004)
}

#| The least-squares slope of ln(y) against t over the samples with t0 ≤ t ≤ t1.
def log_slope(t: Arr, y: Arr, t0, t1) {
    n = 0.0
    st = 0.0
    sy = 0.0
    stt = 0.0
    sty = 0.0
    for i to len(t) {
        if t[i] >= t0 - 1.0e-9 {
            if t[i] <= t1 + 1.0e-9 {
                ly = log(y[i])
                n = n + 1.0
                st = st + t[i]
                sy = sy + ly
                stt = stt + t[i] * t[i]
                sty = sty + t[i] * ly
            }
        }
    }
    return (n * sty - st * sy) / (n * stt - st * st)
}
T_FIT0 = 8.0
T_FIT1 = 18.0
gamma_fit = log_slope(times, amp, T_FIT0, T_FIT1)
gamma_early = log_slope(times, amp, 4.0, 12.0)
gamma_mid = log_slope(times, amp, 6.0, 16.0)
show gamma_fit
show gamma_fit / gamma_th - 1.0
show gamma_early
show gamma_mid

test "the fitted growth rate is the cold-beam rate 0.353553 to within 1%" {
    assert_near(gamma_fit / gamma_th, 1.0, 0.01)
}
test "the same algorithm in numpy (posts/refs/two_stream_pic.py) fits 0.353661097" {
    assert_near(gamma_fit, 0.353661097275761, 1.0e-9)
}
test "fits over [4, 12] and [6, 16] stay within 0.4% of theory" {
    assert_near(gamma_early / gamma_th, 1.0, 0.004)
    assert_near(gamma_mid / gamma_th, 1.0, 0.004)
}

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 1 0^{9}$ . The dashed line in the figure is $e^{γ t}$ with the theoretical rate, drawn through the fitted amplitude at $t = 8$ .

[5]lines 307–350

peak = 0.0
t_peak = 0.0
for n to NSTEP {
    if amp[n] > peak {
        peak = amp[n]
        t_peak = times[n]
    }
}
show peak
show t_peak
e_peak = 0.0
for n to NSTEP {
    e_peak = max(e_peak, energy[n])
}
show e_peak / energy[0]
test "the field energy grows by a factor of 2.9857e9 from t = 0 to its peak (numpy run)" {
    assert_near(e_peak / energy[0] / 2985666941.723498, 1.0, 1.0e-6)
}

test "the mode peaks at amplitude 0.56818 at t = 23.8, as in the numpy run" {
    assert_near(peak, 0.5681768094146153, 1.0e-6)
    assert_near(t_peak, 23.8, 1.0e-9)
}
test "after the peak the amplitude falls: at t = 24.9 it is 0.4754" {
    assert_near(times[NSTEP - 1], 24.9, 1.0e-9)
    assert_near(amp[NSTEP - 1], 0.47542542181643355, 1.0e-5)
}

# the theory line through the run's amplitude at t = 8, the start of the fit
I_FIT = round(T_FIT0 / DT)
line = zeros(NSTEP)
for n to NSTEP {
    line[n] = amp[I_FIT] * exp(gamma_th * (times[n] - times[I_FIT]))
}
figure("Amplitude of the seeded mode", "$t$ ($1/\\omega_p$)", "$|E_1|$", 0.0, 25.0, 1.0e-4, 1.0)
yscale("log")
plot(times, amp, "PIC, 2000 particles")
plot(times, line, "linear theory, $e^{\\gamma t}$")
style("linear theory, $e^{\\gamma t}$", "dash")
vline("fit from", T_FIT0)
vline("fit to", T_FIT1)

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

peak = 0.0
t_peak = 0.0
for n to NSTEP {
    if amp[n] > peak {
        peak = amp[n]
        t_peak = times[n]
    }
}
show peak
show t_peak
e_peak = 0.0
for n to NSTEP {
    e_peak = max(e_peak, energy[n])
}
show e_peak / energy[0]
test "the field energy grows by a factor of 2.9857e9 from t = 0 to its peak (numpy run)" {
    assert_near(e_peak / energy[0] / 2985666941.723498, 1.0, 1.0e-6)
}

test "the mode peaks at amplitude 0.56818 at t = 23.8, as in the numpy run" {
    assert_near(peak, 0.5681768094146153, 1.0e-6)
    assert_near(t_peak, 23.8, 1.0e-9)
}
test "after the peak the amplitude falls: at t = 24.9 it is 0.4754" {
    assert_near(times[NSTEP - 1], 24.9, 1.0e-9)
    assert_near(amp[NSTEP - 1], 0.47542542181643355, 1.0e-5)
}

# the theory line through the run's amplitude at t = 8, the start of the fit
I_FIT = round(T_FIT0 / DT)
line = zeros(NSTEP)
for n to NSTEP {
    line[n] = amp[I_FIT] * exp(gamma_th * (times[n] - times[I_FIT]))
}
figure("Amplitude of the seeded mode", "$t$ ($1/\\omega_p$)", "$|E_1|$", 0.0, 25.0, 1.0e-4, 1.0)
yscale("log")
plot(times, amp, "PIC, 2000 particles")
plot(times, line, "linear theory, $e^{\\gamma t}$")
style("linear theory, $e^{\\gamma t}$", "dash")
vline("fit from", T_FIT0)
vline("fit to", T_FIT1)

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors