Notebooks · posts/trapped-electron-bounce-period

The bounce period of an electron trapped in a plasma wave

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

In the frame of a large-amplitude plasma wave, an electron near the wave's phase velocity moves like a pendulum. Released 0.1 radian from the bottom of a trough, it takes $6.287115$ time units to bounce back, which is $4 K (sin^{2} 0.05)$ from the complete elliptic integral, and $2 π (1 + θ_{0}^{2} /16)$ to six digits. Electrons with speed below $2 ω_{B} / k$ at the trough stay in it; faster ones pass over the crests.

T_{B} = \frac{4}{ω _{B}} K (sin^{2} \frac{θ _{0}}{2}), ω_{B} = k \frac{e ϕ _{0}}{m}

Landau damping is the linear theory of a wave losing energy to the electrons that move at nearly its phase velocity. It assumes the wave is small enough that those electrons go straight through it. O'Neil (Phys. Fluids 8, 2255,

asked what happens when the wave is large: the resonant electrons are

caught in its potential troughs and bounce back and forth there, and the damping stops after about one bounce. The time scale of that nonlinear stage is the bounce period, and this post computes it.

In the frame moving with the wave, the potential is static, $ϕ (x) = ϕ_{0} cos k x$ , and an electron obeys $m x ¨ = - e k ϕ_{0} sin k x$ . That is the pendulum equation. In units with $k = m = e = 1$ it reads $x ¨ = - ϕ_{0} sin x$ , and small oscillations about the bottom of a trough have the bounce frequency $ω_{B} = ϕ_{0}$ .

The program integrates the equation with velocity Verlet and measures the period. The references are a closed form, the complete elliptic integral $K$ , which the program computes with the arithmetic-geometric mean, and posts/refs/trapped-electron-bounce-period.py, which gets $K$ from scipy's ellipk and the period from scipy's solve_ivp without either of the program's methods.

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 35–77

PI = 3.141592653589793
# Verlet time step; its phase error per period is about (omega_B dt)^2 / 24 of a period
DT = 0.002
# steps allowed before a run gives up looking for a turning point
NMAX = 100000

#| One velocity-Verlet run of x'' = -phi0 sin(x) from rest at x0.
#| The electron swings to the far side and comes back to rest there after
#| half a period; the time of that turning point, found by interpolating
#| where v crosses zero, doubled, is the bounce period.
def bounce_period(x0, phi0, dt) {
    x = x0
    v = 0.0
    t = 0.0
    a = -phi0 * sin(x)
    for i to NMAX {
        vprev = v
        x = x + v * dt + 0.5 * a * dt * dt
        anew = -phi0 * sin(x)
        v = v + 0.5 * (a + anew) * dt
        a = anew
        t = t + dt
        if vprev < 0.0 and v >= 0.0 {
            frac = -vprev / (v - vprev)
            return 2.0 * (t - dt + frac * dt)
        }
    }
    return -1.0
}

#| The complete elliptic integral of the first kind with parameter m = k^2,
#| K(m) = pi / (2 M(1, sqrt(1 - m))), where M is the arithmetic-geometric
#| mean; each step of the mean doubles its correct digits.
def ellip_k(m) {
    a = 1.0
    b = sqrt(1.0 - m)
    while abs(a - b) > 1.0e-15 * a {
        an = 0.5 * (a + b)
        b = sqrt(a * b)
        a = an
    }
    return PI / (2.0 * a)
}

PI = 3.141592653589793
# Verlet time step; its phase error per period is about (omega_B dt)^2 / 24 of a period
DT = 0.002
# steps allowed before a run gives up looking for a turning point
NMAX = 100000

#| One velocity-Verlet run of x'' = -phi0 sin(x) from rest at x0.
#| The electron swings to the far side and comes back to rest there after
#| half a period; the time of that turning point, found by interpolating
#| where v crosses zero, doubled, is the bounce period.
def bounce_period(x0, phi0, dt) {
    x = x0
    v = 0.0
    t = 0.0
    a = -phi0 * sin(x)
    for i to NMAX {
        vprev = v
        x = x + v * dt + 0.5 * a * dt * dt
        anew = -phi0 * sin(x)
        v = v + 0.5 * (a + anew) * dt
        a = anew
        t = t + dt
        if vprev < 0.0 and v >= 0.0 {
            frac = -vprev / (v - vprev)
            return 2.0 * (t - dt + frac * dt)
        }
    }
    return -1.0
}

#| The complete elliptic integral of the first kind with parameter m = k^2,
#| K(m) = pi / (2 M(1, sqrt(1 - m))), where M is the arithmetic-geometric
#| mean; each step of the mean doubles its correct digits.
def ellip_k(m) {
    a = 1.0
    b = sqrt(1.0 - m)
    while abs(a - b) > 1.0e-15 * a {
        an = 0.5 * (a + b)
        b = sqrt(a * b)
        a = an
    }
    return PI / (2.0 * a)
}

The period at $θ_{0} = 0.1$

A pendulum released from rest at angle $θ_{0}$ has the period

T = \frac{4}{ω _{B}} K (m), m = sin^{2} \frac{θ _{0}}{2},

where $K (m) = \int_{0}^{π /2} d φ / 1 - m sin^{2} φ$ . For a small release angle the series $K (m) = \frac{π}{2} (1 + m /4 + 9 m^{2} /64 + \dots)$ gives $T = \frac{2 π}{ω _{B}} (1 + θ_{0}^{2} /16 + 11 θ_{0}^{4} /3072 + \dots)$ : the period grows with amplitude because the restoring force $sin x$ is weaker than $x$ .

At $θ_{0} = 0.1$ and $ϕ_{0} = 1$ the elliptic integral gives $T = 6.2871145$ , and $2 π (1 + θ_{0}^{2} /16) = 6.2871123$ is short by $2.25 \times 1 0^{- 6}$ , the size of the next term. The Verlet run lands within $1 0^{- 6}$ of the elliptic value: its phase error is $(ω_{B} Δ t)^{2} /24 = 1.7 \times 1 0^{- 7}$ of a period, and the test allows twice that.

[2]lines 100–130

TH0 = 0.1
phi0 = 1.0
omega_B = sqrt(phi0)
T_num = bounce_period(TH0, phi0, DT)
T_ellip = 4.0 / omega_B * ellip_k(sin(TH0 / 2.0) * sin(TH0 / 2.0))
T_series = 2.0 * PI / omega_B * (1.0 + TH0 * TH0 / 16.0)
show T_num
show T_ellip
show T_series
# Verlet's phase error per period, as a fraction of it: (omega_B dt)^2 / 24
verlet_phase = omega_B * omega_B * DT * DT / 24.0
show verlet_phase

test "the arithmetic-geometric mean gives K(1/2) = 1.854074677301372, scipy's ellipk(0.5)" {
    assert_near(ellip_k(0.5), 1.854074677301372, 1.0e-14)
}
test "4 K(sin^2(0.05)) = 6.287114549310, scipy's value" {
    # 0.05 is half the release angle
    assert_near(TH0 / 2.0, 0.05, 1.0e-15)
    assert_near(T_ellip, 6.287114549310, 1.0e-11)
}
test "the Verlet phase error (omega_B dt)^2 / 24 is 1.7e-7 of a period" {
    assert_near(verlet_phase, 0.002 * 0.002 / 24.0, 1.0e-20)
}
test "the Verlet period is within twice the Verlet phase error, 2 (dt^2/24) T, of 4K" {
    assert_lt(abs(T_num - T_ellip), 2.0 * DT * DT / 24.0 * T_ellip)
}
test "2 pi (1 + theta0^2/16) = 6.287112298 is below 4K by 2.25e-6, the size of the theta0^4 term" {
    assert_near(T_series, 6.287112297997, 1.0e-11)
    assert_near(T_ellip - T_series, 2.0 * PI * 11.0 * pow(TH0, 4.0) / 3072.0, 2.0e-9)
}

TH0 = 0.1
phi0 = 1.0
omega_B = sqrt(phi0)
T_num = bounce_period(TH0, phi0, DT)
T_ellip = 4.0 / omega_B * ellip_k(sin(TH0 / 2.0) * sin(TH0 / 2.0))
T_series = 2.0 * PI / omega_B * (1.0 + TH0 * TH0 / 16.0)
show T_num
show T_ellip
show T_series
# Verlet's phase error per period, as a fraction of it: (omega_B dt)^2 / 24
verlet_phase = omega_B * omega_B * DT * DT / 24.0
show verlet_phase

test "the arithmetic-geometric mean gives K(1/2) = 1.854074677301372, scipy's ellipk(0.5)" {
    assert_near(ellip_k(0.5), 1.854074677301372, 1.0e-14)
}
test "4 K(sin^2(0.05)) = 6.287114549310, scipy's value" {
    # 0.05 is half the release angle
    assert_near(TH0 / 2.0, 0.05, 1.0e-15)
    assert_near(T_ellip, 6.287114549310, 1.0e-11)
}
test "the Verlet phase error (omega_B dt)^2 / 24 is 1.7e-7 of a period" {
    assert_near(verlet_phase, 0.002 * 0.002 / 24.0, 1.0e-20)
}
test "the Verlet period is within twice the Verlet phase error, 2 (dt^2/24) T, of 4K" {
    assert_lt(abs(T_num - T_ellip), 2.0 * DT * DT / 24.0 * T_ellip)
}
test "2 pi (1 + theta0^2/16) = 6.287112298 is below 4K by 2.25e-6, the size of the theta0^4 term" {
    assert_near(T_series, 6.287112297997, 1.0e-11)
    assert_near(T_ellip - T_series, 2.0 * PI * 11.0 * pow(TH0, 4.0) / 3072.0, 2.0e-9)
}

A deeper trough bounces faster

The only frequency in the problem is $ω_{B} = ϕ_{0}$ , so quadrupling the wave amplitude halves the bounce period at the same release angle. The program runs $ϕ_{0} = 4$ with the time step halved, which keeps $ω_{B} Δ t$ the same, and the period comes out half of the $ϕ_{0} = 1$ period to rounding.

[3]lines 142–147

T_strong = bounce_period(TH0, 4.0, DT / 2.0)
show T_strong
test "at phi0 = 4 the period is half, 3.1435573, as solve_ivp finds" {
    assert_near(T_strong, 0.5 * T_num, 1.0e-12)
    assert_near(T_strong, 3.143557274657, 1.0e-6)
}

At large amplitude the period diverges at the separatrix

An electron released at rest near the crest, $θ_{0} \to π$ , takes longer and longer to start moving, and $K (m)$ grows like $\frac{1}{2} ln \frac{16}{1 - m}$ as $m \to 1$ . The first figure puts the Verlet periods beside $4 K$ at six release angles; the program checks each against the elliptic integral, and the reference script checks $4 K$ against solve_ivp at the same angles.

[4]lines 160–189

AMPS = [0.5, 1.0, 2.0, 2.5, 3.0, 3.1]
T_amp = zeros(len(AMPS))
T_k = zeros(len(AMPS))
for j to len(AMPS) {
    T_amp[j] = bounce_period(AMPS[j], phi0, DT)
    T_k[j] = 4.0 * ellip_k(sin(AMPS[j] / 2.0) * sin(AMPS[j] / 2.0))
}
show T_amp
show T_k

# the curve 4K(sin^2(theta0/2)) on a fine grid of release angles, stopping short of pi
N_CURVE = 150
th_c = linspace(0.02, 3.13, N_CURVE)
T_c = zeros(N_CURVE)
for j to N_CURVE {
    T_c[j] = 4.0 * ellip_k(sin(th_c[j] / 2.0) * sin(th_c[j] / 2.0))
}
figure("Bounce period against release angle, phi0 = 1", "release angle theta0 (radians)", "bounce period $\\omega_B T$")
plot(th_c, T_c, "$4K(\\sin^2(\\theta_0/2))$")
scatter(AMPS, T_amp, "velocity Verlet, dt = 0.002")
hline("2 pi, the small-amplitude period", 2.0 * PI)

test "at each release angle the Verlet period matches 4K to 2e-6 relative" {
    for j to len(AMPS) {
        assert_lt(abs(T_amp[j] / T_k[j] - 1.0), 2.0e-6)
    }
}
test "4K at theta0 = 3.1 is 21.0392237, as solve_ivp finds" {
    assert_near(T_k[len(AMPS) - 1], 21.0392236680, 1.0e-8)
}

AMPS = [0.5, 1.0, 2.0, 2.5, 3.0, 3.1]
T_amp = zeros(len(AMPS))
T_k = zeros(len(AMPS))
for j to len(AMPS) {
    T_amp[j] = bounce_period(AMPS[j], phi0, DT)
    T_k[j] = 4.0 * ellip_k(sin(AMPS[j] / 2.0) * sin(AMPS[j] / 2.0))
}
show T_amp
show T_k

# the curve 4K(sin^2(theta0/2)) on a fine grid of release angles, stopping short of pi
N_CURVE = 150
th_c = linspace(0.02, 3.13, N_CURVE)
T_c = zeros(N_CURVE)
for j to N_CURVE {
    T_c[j] = 4.0 * ellip_k(sin(th_c[j] / 2.0) * sin(th_c[j] / 2.0))
}
figure("Bounce period against release angle, phi0 = 1", "release angle theta0 (radians)", "bounce period $\\omega_B T$")
plot(th_c, T_c, "$4K(\\sin^2(\\theta_0/2))$")
scatter(AMPS, T_amp, "velocity Verlet, dt = 0.002")
hline("2 pi, the small-amplitude period", 2.0 * PI)

test "at each release angle the Verlet period matches 4K to 2e-6 relative" {
    for j to len(AMPS) {
        assert_lt(abs(T_amp[j] / T_k[j] - 1.0), 2.0e-6)
    }
}
test "4K at theta0 = 3.1 is 21.0392237, as solve_ivp finds" {
    assert_near(T_k[len(AMPS) - 1], 21.0392236680, 1.0e-8)
}

Trapped and passing electrons

The energy $E = v^{2} /2 - ϕ_{0} cos x$ is conserved. An electron at the bottom of a trough ( $x = 0$ ) with speed $v_{0}$ reaches the crest ( $x = \pm π$ , where the potential energy is $+ ϕ_{0}$ ) only if $v_{0}^{2} /2 - ϕ_{0} > ϕ_{0}$ , that is $v_{0} > 2 ϕ_{0} = 2 ω_{B}$ . The curve $E = ϕ_{0}$ , $v = \pm 2 ω_{B} cos (x /2)$ , is the separatrix: inside it electrons are trapped and bounce, outside they stream past the crests. The trapped band is $4 ω_{B} / k$ wide in velocity, and it is the electrons in this band that O'Neil's theory follows.

The program finds the boundary without using the energy: it launches electrons from the trough, marks an electron as passing when it gets past $x = π$ , and bisects on $v_{0}$ . The boundary comes out at $2 ω_{B}$ to the bisection's tolerance.

[5]lines 209–295

# how long a launched electron is followed: several bounce periods
T_FOLLOW = 40.0
DT_ORBIT = 0.005

#| 1 if an electron launched from the trough at speed v0 gets past the
#| crest at x = pi within T_FOLLOW, 0 if it stays trapped.
def passes(v0, phi0) {
    x = 0.0
    v = v0
    a = 0.0
    nstep = round(T_FOLLOW / DT_ORBIT)
    for i to nstep {
        x = x + v * DT_ORBIT + 0.5 * a * DT_ORBIT * DT_ORBIT
        anew = -phi0 * sin(x)
        v = v + 0.5 * (a + anew) * DT_ORBIT
        a = anew
        if x > PI {
            return 1.0
        }
    }
    return 0.0
}

# bisection on the launch speed: lo is trapped, hi passes
N_BISECT = 30
v_lo = 1.5
v_hi = 2.5
for i to N_BISECT {
    v_mid = 0.5 * (v_lo + v_hi)
    if passes(v_mid, phi0) > 0.5 {
        v_hi = v_mid
    } else {
        v_lo = v_mid
    }
}
v_trap = 0.5 * (v_lo + v_hi)
show v_trap

test "the trapped-passing boundary at the trough is 2 omega_B = 2, to 1e-4" {
    assert_near(v_trap, 2.0 * omega_B, 1.0e-4)
}

# phase-space orbits: six trapped ones inside the separatrix and two passing ones above it
V_LAUNCH = [0.4, 0.8, 1.2, 1.6, 1.9, 1.99, 2.1, 2.5]
# 3200 steps of 0.005 = 16 time units, longer than the slowest trapped orbit here (v0 = 1.99, period 14.8)
N_ORBIT = 3200
# keep every 25th point, to hold the figure data down
EVERY = 25
figure("Phase space of electrons in the wave frame, phi0 = 1", "x (wave phase, radians)", "v", -3.3, 3.3, -3.0, 3.0)
for k to len(V_LAUNCH) {
    ox = zeros(N_ORBIT // EVERY)
    ov = zeros(N_ORBIT // EVERY)
    x = 0.0
    v = V_LAUNCH[k]
    a = 0.0
    for i to N_ORBIT {
        if i % EVERY == 0 {
            ox[i // EVERY] = x
            ov[i // EVERY] = v
        }
        x = x + v * DT_ORBIT + 0.5 * a * DT_ORBIT * DT_ORBIT
        anew = -phi0 * sin(x)
        v = v + 0.5 * (a + anew) * DT_ORBIT
        a = anew
        # fold passing orbits back into one wavelength
        if x > PI {
            x = x - 2.0 * PI
        }
    }
    scatter(ox, ov, "orbits from the trough")
}
N_SEP = 201
xs = linspace(-PI, PI, N_SEP)
# the separatrix as one closed curve: the upper half left to right, then the lower half back
sep_x = zeros(2 * N_SEP)
sep_v = zeros(2 * N_SEP)
for i to N_SEP {
    sep_x[i] = xs[i]
    sep_v[i] = 2.0 * omega_B * cos(xs[i] / 2.0)
    sep_x[2 * N_SEP - 1 - i] = xs[i]
    sep_v[2 * N_SEP - 1 - i] = -sep_v[i]
}
plot(sep_x, sep_v, "separatrix $v = \\pm 2\\omega_B\\cos(x/2)$")
style("separatrix $v = \\pm 2\\omega_B\\cos(x/2)$", "dash")

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

# how long a launched electron is followed: several bounce periods
T_FOLLOW = 40.0
DT_ORBIT = 0.005

#| 1 if an electron launched from the trough at speed v0 gets past the
#| crest at x = pi within T_FOLLOW, 0 if it stays trapped.
def passes(v0, phi0) {
    x = 0.0
    v = v0
    a = 0.0
    nstep = round(T_FOLLOW / DT_ORBIT)
    for i to nstep {
        x = x + v * DT_ORBIT + 0.5 * a * DT_ORBIT * DT_ORBIT
        anew = -phi0 * sin(x)
        v = v + 0.5 * (a + anew) * DT_ORBIT
        a = anew
        if x > PI {
            return 1.0
        }
    }
    return 0.0
}

# bisection on the launch speed: lo is trapped, hi passes
N_BISECT = 30
v_lo = 1.5
v_hi = 2.5
for i to N_BISECT {
    v_mid = 0.5 * (v_lo + v_hi)
    if passes(v_mid, phi0) > 0.5 {
        v_hi = v_mid
    } else {
        v_lo = v_mid
    }
}
v_trap = 0.5 * (v_lo + v_hi)
show v_trap

test "the trapped-passing boundary at the trough is 2 omega_B = 2, to 1e-4" {
    assert_near(v_trap, 2.0 * omega_B, 1.0e-4)
}

# phase-space orbits: six trapped ones inside the separatrix and two passing ones above it
V_LAUNCH = [0.4, 0.8, 1.2, 1.6, 1.9, 1.99, 2.1, 2.5]
# 3200 steps of 0.005 = 16 time units, longer than the slowest trapped orbit here (v0 = 1.99, period 14.8)
N_ORBIT = 3200
# keep every 25th point, to hold the figure data down
EVERY = 25
figure("Phase space of electrons in the wave frame, phi0 = 1", "x (wave phase, radians)", "v", -3.3, 3.3, -3.0, 3.0)
for k to len(V_LAUNCH) {
    ox = zeros(N_ORBIT // EVERY)
    ov = zeros(N_ORBIT // EVERY)
    x = 0.0
    v = V_LAUNCH[k]
    a = 0.0
    for i to N_ORBIT {
        if i % EVERY == 0 {
            ox[i // EVERY] = x
            ov[i // EVERY] = v
        }
        x = x + v * DT_ORBIT + 0.5 * a * DT_ORBIT * DT_ORBIT
        anew = -phi0 * sin(x)
        v = v + 0.5 * (a + anew) * DT_ORBIT
        a = anew
        # fold passing orbits back into one wavelength
        if x > PI {
            x = x - 2.0 * PI
        }
    }
    scatter(ox, ov, "orbits from the trough")
}
N_SEP = 201
xs = linspace(-PI, PI, N_SEP)
# the separatrix as one closed curve: the upper half left to right, then the lower half back
sep_x = zeros(2 * N_SEP)
sep_v = zeros(2 * N_SEP)
for i to N_SEP {
    sep_x[i] = xs[i]
    sep_v[i] = 2.0 * omega_B * cos(xs[i] / 2.0)
    sep_x[2 * N_SEP - 1 - i] = xs[i]
    sep_v[2 * N_SEP - 1 - i] = -sep_v[i]
}
plot(sep_x, sep_v, "separatrix $v = \\pm 2\\omega_B\\cos(x/2)$")
style("separatrix $v = \\pm 2\\omega_B\\cos(x/2)$", "dash")

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

The period at θ0​=0.1

A deeper trough bounces faster

At large amplitude the period diverges at the separatrix

Trapped and passing electrons

The period at $θ_{0} = 0.1$