Part 8, Chapter 1

Nonlinear Dynamics

Solitons, BGK modes, and wave breaking in plasmas

1.1 Solitons and the KdV Equation

Ion acoustic waves in a plasma can steepen nonlinearly (the crests travel faster than the troughs) while simultaneously being dispersed. When these two effects balance exactly, the result is a stable, localized wave packet that propagates without changing shape: a soliton.

Starting from the ion fluid equations with Boltzmann electrons and expanding in small amplitude $\epsilon$ using the reductive perturbation method (stretched coordinates $\xi = \epsilon^{1/2}(x - c_s t)$,$\tau = \epsilon^{3/2}t$), one obtains the Korteweg-de Vries (KdV) equation:

$$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + \frac{\partial^3 u}{\partial x^3} = 0$$

The nonlinear term $u\,\partial_x u$ causes wave steepening, while the dispersive term $\partial_x^3 u$ spreads the wave. The exact one-soliton solution is:

$$u(x,t) = 3V\,\text{sech}^2\!\left(\frac{\sqrt{V}}{2}(x - Vt - x_0)\right)$$

Key properties: (1) The amplitude is proportional to the velocity $V$ -- taller solitons move faster. (2) The width scales as $1/\sqrt{V}$ -- taller solitons are narrower. (3) Two solitons pass through each other and emerge unchanged (elastic collision). (4) The KdV equation is integrable and has infinitely many conserved quantities.

1.2 BGK Modes

Bernstein-Greene-Kruskal (BGK) modes are exact nonlinear solutions of the Vlasov-Poisson system. Unlike linear waves that eventually break when their amplitude grows, BGK modes are self-consistent stationary (in the wave frame) structures.

In the wave frame moving at velocity $v_\phi$, the distribution function depends only on the energy:

$$\mathcal{E} = \frac{1}{2}m(v - v_\phi)^2 + q\phi(x)$$

Particles with $\mathcal{E} > q\phi_{max}$ are "passing" (they stream over the potential hills), while particles with $\mathcal{E} < q\phi_{max}$ are "trapped" -- they bounce back and forth in the potential wells. The trapped particle distribution is free to be specified (it is not determined by the passing particle distribution), giving the BGK construction its flexibility:

$$f(\mathcal{E}) = \begin{cases} f_\text{passing}(\mathcal{E}) & \mathcal{E} > q\phi_{max} \\ f_\text{trapped}(\mathcal{E}) & \mathcal{E} < q\phi_{max} \end{cases}$$

One then solves Poisson's equation self-consistently:

$$\frac{d^2\phi}{dx^2} = -\frac{e}{\epsilon_0}\left(n_i(\phi) - n_e(\phi)\right)$$

where the densities are obtained by integrating the distribution over velocity. BGK modes manifest as phase-space vortices (also called electron or ion holes), which are commonly observed in simulations and space plasmas.

1.3 Wave Breaking and Phase-Space Vortices

When a plasma wave grows to large amplitude, the fluid description breaks down because the electron oscillation velocity approaches the phase velocity of the wave. At this point, particles become trapped by the wave potential and the wave "breaks." The wave breaking limit for a cold plasma wave is:

$$E_{WB} = \frac{m_e \omega_{pe} v_\phi}{e}$$

For a thermal plasma, the Coffey limit gives a lower threshold. Wave breaking leads to efficient particle trapping and acceleration, which is the basis of laser wakefield acceleration, where electrons surf the broken plasma wave and gain energy.

1.4 Modulational Instability

A uniform wave train can be unstable to modulations of its envelope. For Langmuir waves, the envelope evolution is governed by the nonlinear Schrodinger equation:

$$i\frac{\partial A}{\partial t} + P\frac{\partial^2 A}{\partial x^2} + Q|A|^2 A = 0$$

where $P$ is the group velocity dispersion and$Q$ is the nonlinear frequency shift. Modulational instability occurs when $PQ > 0$, with growth rate:

$$\gamma = |K|\sqrt{2Q|A_0|^2 P - P^2 K^2}$$

where $K$ is the modulation wavenumber and$A_0$ is the pump amplitude. The instability leads to the formation of localized wave packets -- envelope solitons -- and ultimately to Langmuir collapse in higher dimensions, where wave energy focuses to ever-smaller scales.

Interactive Simulation

KdV Soliton Propagation and Phase-Space Vortex

Python

script.py124 lines

import numpy as np

# ==============================
# Part 1: KdV Soliton Propagation
# ==============================
# Solve: du/dt + u*du/dx + d^3u/dx^3 = 0
# using pseudo-spectral method (FFT for spatial derivatives, RK4 in time)

N = 256
L = 40.0
dx = L / N
x = np.linspace(0, L, N, endpoint=False)

# Initial condition: two solitons with different amplitudes
V1 = 2.0    # Fast soliton
V2 = 0.5    # Slow soliton
x01 = 8.0
x02 = 20.0

def soliton(x_arr, V, x0):
    arg = np.sqrt(V) / 2 * (x_arr - x0)
    return 3 * V / np.cosh(arg)**2

u0 = soliton(x, V1, x01) + soliton(x, V2, x02)

# Wavenumbers for spectral derivatives
k = 2 * np.pi * np.fft.fftfreq(N, dx)

def rhs(u):
    u_hat = np.fft.fft(u)
    # du/dx (spectral)
    du_dx = np.real(np.fft.ifft(1j * k * u_hat))
    # d^3u/dx^3 (spectral)
    d3u_dx3 = np.real(np.fft.ifft(-1j * k**3 * u_hat))
    return -u * du_dx - d3u_dx3

# RK4 time integration
dt = 0.005
Nt = 1000
T_final = Nt * dt

u = u0.copy()
snapshots = {}
snap_times = [0, Nt//4, Nt//2, 3*Nt//4, Nt-1]
for step in range(Nt):
    if step in snap_times:
        snapshots[step * dt] = u.copy()

k1 = rhs(u)
    k2 = rhs(u + 0.5*dt*k1)
    k3 = rhs(u + 0.5*dt*k2)
    k4 = rhs(u + dt*k3)
    u = u + dt/6 * (k1 + 2*k2 + 2*k3 + k4)

snapshots[Nt * dt] = u.copy()

# Conservation check
mass_initial = np.sum(u0) * dx
mass_final = np.sum(u) * dx
energy_initial = np.sum(u0**2) * dx
energy_final = np.sum(u**2) * dx

print("=== KdV Soliton Collision ===")
print(f"Grid: N={N}, L={L:.0f}, dt={dt:.4f}, T={T_final:.1f}")
print(f"Soliton 1: amplitude={3*V1:.1f}, velocity={V1:.1f}")
print(f"Soliton 2: amplitude={3*V2:.1f}, velocity={V2:.1f}")
print(f"")
print(f"Conservation check:")
print(f"  Mass:   initial={mass_initial:.6f}, final={mass_final:.6f}, error={abs(mass_final-mass_initial)/abs(mass_initial):.2e}")
print(f"  Energy: initial={energy_initial:.6f}, final={energy_final:.6f}, error={abs(energy_final-energy_initial)/abs(energy_initial):.2e}")
print()

print(f"{'Time':>8} {'Max u':>10} {'Position of max':>16}")
print("-" * 38)
for t in sorted(snapshots.keys()):
    u_snap = snapshots[t]
    idx_max = np.argmax(u_snap)
    print(f"{t:8.2f} {u_snap[idx_max]:10.3f} {x[idx_max]:16.2f}")

# ==============================
# Part 2: Phase-Space Vortex (Electron Hole)
# ==============================
print()
print("=== Phase-Space Vortex (BGK Electron Hole) ===")
print()

# Construct a BGK-like phase space portrait
# Potential well: phi(x) = phi0 * sech^2(x/w)
phi0 = 1.0
w = 2.0
x_ps = np.linspace(-10, 10, 100)
phi_x = phi0 / np.cosh(x_ps / w)**2

# Passing particles: total energy E > phi0
# Trapped particles: 0 < E < phi0
# In wave frame, E = 0.5*v^2 - phi(x) for electrons (q=-e, but we use normalized units)

print(f"Potential: phi(x) = {phi0:.1f} * sech^2(x/{w:.1f})")
print(f"Trapping width in velocity: delta_v = 2*sqrt(phi0) = {2*np.sqrt(phi0):.3f}")
print(f"Bounce frequency: omega_b ~ sqrt(phi0)/w = {np.sqrt(phi0)/w:.3f}")
print()

# Phase space structure
v_max = 4.0
Nv = 60
v_arr = np.linspace(-v_max, v_max, Nv)

print("Phase space density (passing + trapped regions):")
print(f"{'x':>8} {'v_trap_min':>12} {'v_trap_max':>12} {'phi(x)':>10}")
print("-" * 46)
for i in range(0, len(x_ps), 10):
    xi = x_ps[i]
    phi_val = phi_x[i]
    if phi_val > 0.01:
        v_trap = np.sqrt(2 * phi_val)
        print(f"{xi:8.2f} {-v_trap:12.3f} {v_trap:12.3f} {phi_val:10.4f}")
    else:
        print(f"{xi:8.2f} {'---':>12} {'---':>12} {phi_val:10.4f}")

print()
print("The separatrix in phase space forms a cat-eye (vortex) pattern.")
print("Trapped particles execute closed orbits inside the separatrix.")
print("Passing particles stream over the potential hump.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Previous Next →