Part 8, Chapter 3

PIC Simulations

Particle-in-cell computational methods for plasma kinetics

3.1 The PIC Algorithm

The particle-in-cell (PIC) method is the workhorse of computational plasma physics. Rather than discretizing the full distribution function on a velocity-space grid (which is prohibitively expensive in 3D-3V), PIC represents the plasma by a collection of computational "superparticles," each representing many physical particles. The algorithm cycles through four steps each timestep:

Push particles: Advance particle positions and velocities using the equations of motion in the local fields.
Deposit charge/current: Interpolate particle charges and currents onto the spatial grid.
Solve fields: Solve Maxwell's equations (or Poisson's equation in the electrostatic limit) on the grid.
Interpolate fields: Interpolate the grid fields back to particle positions for the next push.

The equations of motion for each particle are:

$$\frac{d\mathbf{x}_p}{dt} = \mathbf{v}_p, \qquad m_p\frac{d\mathbf{v}_p}{dt} = q_p\left(\mathbf{E}(\mathbf{x}_p) + \mathbf{v}_p \times \mathbf{B}(\mathbf{x}_p)\right)$$

3.2 The Boris Pusher

The Boris algorithm is the standard method for advancing particle velocities in electromagnetic fields. It is second-order accurate, time-reversible, and exactly preserves the particle speed during the magnetic rotation step (volume-preserving in phase space).

The algorithm splits the velocity update into three substeps. First, a half electric-field kick:

$$\mathbf{v}^- = \mathbf{v}^{n-1/2} + \frac{q\Delta t}{2m}\mathbf{E}^n$$

Then a magnetic rotation. Define the rotation vectors:

$$\mathbf{t} = \frac{q\mathbf{B}^n \Delta t}{2m}, \qquad \mathbf{s} = \frac{2\mathbf{t}}{1 + |\mathbf{t}|^2}$$

$$\mathbf{v}' = \mathbf{v}^- + \mathbf{v}^- \times \mathbf{t}, \qquad \mathbf{v}^+ = \mathbf{v}^- + \mathbf{v}' \times \mathbf{s}$$

Finally, a second half electric-field kick:

$$\mathbf{v}^{n+1/2} = \mathbf{v}^+ + \frac{q\Delta t}{2m}\mathbf{E}^n$$

The key property is that $|\mathbf{v}^+| = |\mathbf{v}^-|$ exactly, so the magnetic field does no work, as it should physically. This prevents secular energy drift due to the magnetic push.

3.3 Stability and Resolution Requirements

Electrostatic PIC simulations have stringent resolution requirements. The grid spacing must resolve the Debye length:

$$\Delta x < \lambda_D = \sqrt{\frac{\epsilon_0 T_e}{n_e e^2}}$$

Failure to resolve the Debye length leads to numerical (grid) heating: the alias of short-wavelength modes folds energy to long wavelengths, artificially heating the plasma. The timestep must satisfy the CFL-like condition:

$$\omega_{pe}\,\Delta t < 2$$

For electromagnetic PIC, additional constraints arise from the Courant condition for light-wave propagation: $c\,\Delta t < \Delta x$.

Finite-size particles: To reduce discrete particle noise, each superparticle is given a spatial extent (cloud or spline shape) rather than being a point charge. This smooths the short-range Coulomb interaction, suppressing unphysical close encounters while faithfully representing the collective long-range behavior. Common choices include first-order (linear, or CIC) and higher-order B-spline shapes.

3.4 Numerical Heating and Noise

PIC simulations suffer from discrete particle noise that scales as $1/\sqrt{N_{ppc}}$ where$N_{ppc}$ is the number of particles per cell. The electric field fluctuation level is:

$$\frac{\epsilon_0 \langle E^2 \rangle}{n_e T_e} \sim \frac{1}{N_{ppc} \cdot (n_e \lambda_D^d)}$$

where $d$ is the spatial dimension. This noise acts as a collisional drag and diffusion, giving an effective collision frequency:

$$\nu_\text{num} \sim \frac{\omega_{pe}}{N_{ppc} \cdot (n_e \lambda_D^d)}$$

For the simulation to faithfully reproduce collisionless physics, we need the numerical collision rate to be much smaller than any physical rate of interest over the simulation duration.

Interactive Simulation

Minimal 1D Electrostatic PIC: Two-Stream Instability

Python

script.py123 lines

import numpy as np

# 1D Electrostatic PIC simulation of two-stream instability
# Two counter-streaming electron beams in a neutralizing ion background

# Parameters
L = 2 * np.pi       # System length
Ng = 64              # Number of grid cells
Np = 4000            # Number of particles
dt = 0.1             # Timestep (in units of 1/omega_pe)
Nt = 200             # Number of timesteps
dx = L / Ng

# Initialize two counter-streaming beams
v_beam = 3.0         # Beam velocity
v_th = 0.5           # Thermal spread

x = np.random.uniform(0, L, Np)
v = np.zeros(Np)
# First beam (half particles moving right)
v[:Np//2] = v_beam + v_th * np.random.randn(Np//2)
# Second beam (half particles moving left)
v[Np//2:] = -v_beam + v_th * np.random.randn(Np//2)

# Add small perturbation to seed instability
x += 0.01 * np.sin(2 * np.pi * x / L)
x = x % L

# Charge deposition (CIC - Cloud-In-Cell)
def deposit(x_p, Ng, L):
    dx = L / Ng
    rho = np.zeros(Ng)
    xn = x_p / dx
    j = np.floor(xn).astype(int) % Ng
    w = xn - np.floor(xn)
    for p in range(len(x_p)):
        rho[j[p]] += (1.0 - w[p])
        rho[(j[p]+1) % Ng] += w[p]
    rho = rho * (-1.0) * Ng / len(x_p) + 1.0  # Normalize, add ion background
    return rho

# Solve Poisson: d^2 phi / dx^2 = -rho/eps0 (normalized)
def solve_poisson(rho, dx):
    Ng = len(rho)
    rho_k = np.fft.fft(rho)
    k = 2 * np.pi * np.fft.fftfreq(Ng, dx)
    phi_k = np.zeros(Ng, dtype=complex)
    phi_k[1:] = rho_k[1:] / (k[1:]**2)
    phi = np.real(np.fft.ifft(phi_k))
    # E = -dphi/dx (spectral)
    E_k = -1j * k * np.fft.fft(phi)
    E = np.real(np.fft.ifft(E_k))
    return E

# Field interpolation (CIC)
def interpolate(E_grid, x_p, Ng, L):
    dx = L / Ng
    xn = x_p / dx
    j = np.floor(xn).astype(int) % Ng
    w = xn - np.floor(xn)
    E_p = (1.0 - w) * E_grid[j] + w * E_grid[(j+1) % Ng]
    return E_p

# Track energy and phase space snapshots
KE_history = []
FE_history = []
snapshots = [0, Nt//4, Nt//2, 3*Nt//4, Nt-1]
phase_data = {}

for step in range(Nt):
    # Deposit charge
    rho = deposit(x, Ng, L)

# Solve for E field
    E_grid = solve_poisson(rho, dx)

# Interpolate E to particles
    E_p = interpolate(E_grid, x, Ng, L)

# Push particles (leapfrog)
    v += E_p * dt
    x += v * dt
    x = x % L

# Compute energies
    KE = 0.5 * np.sum(v**2) / Np
    FE = 0.5 * np.sum(E_grid**2) * dx
    KE_history.append(KE)
    FE_history.append(FE)

if step in snapshots:
        # Sample particles for phase space plot
        idx = np.random.choice(Np, min(500, Np), replace=False)
        phase_data[step] = (x[idx].copy(), v[idx].copy())

KE_arr = np.array(KE_history)
FE_arr = np.array(FE_history)
time = np.arange(Nt) * dt

print("=== 1D Electrostatic PIC: Two-Stream Instability ===")
print(f"Particles: {Np}, Grid cells: {Ng}, Timesteps: {Nt}")
print(f"Beam velocity: +/- {v_beam:.1f} v_th, Thermal spread: {v_th:.1f}")
print(f"dx/lambda_D ~ {dx/v_th:.2f}")
print(f"omega_pe * dt = {dt:.2f}")
print(f"")
print(f"{'Time':>8} {'KE':>12} {'Field E':>12} {'Total E':>12}")
print("-" * 48)
for i in range(0, Nt, Nt//10):
    print(f"{time[i]:8.1f} {KE_arr[i]:12.4f} {FE_arr[i]:12.6f} {KE_arr[i]+FE_arr[i]:12.4f}")

# Growth rate estimate from field energy
mid = Nt // 4
if FE_arr[mid] > FE_arr[0] and FE_arr[0] > 0:
    gamma_est = 0.5 * np.log(FE_arr[mid] / max(FE_arr[0], 1e-20)) / (time[mid])
    print(f"\nEstimated growth rate: gamma ~ {gamma_est:.3f} omega_pe")
    print(f"Theoretical: gamma ~ {np.sqrt(3)/2 * (1.0)**(1/3):.3f} omega_pe (for cold beam)")

print(f"\nPhase space snapshots taken at t = {[s*dt for s in snapshots]}")
for step in snapshots:
    xp, vp = phase_data[step]
    v_spread = np.std(vp)
    print(f"  t={step*dt:6.1f}: v_spread={v_spread:.3f}, min_v={np.min(vp):.2f}, max_v={np.max(vp):.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Previous Next →