Plasma Physics Programs

Interactive simulations for computational plasma physics — particle motion, wave phenomena, and kinetic instabilities

These Python programs run directly in your browser using Pyodide (WebAssembly Python). The first run downloads the Python environment (~15MB). Click "Run" to execute!

Charged Particle Motion

Boris pusher algorithm for charged particle gyration and E×B drift in uniform fields

import numpy as np
import matplotlib
matplotlib.use('AGG')
import matplotlib.pyplot as plt

# Charged Particle Motion in Electromagnetic Fields
# Boris pusher algorithm for accurate long-term integration

# Physical parameters
q = 1.602e-19    # electron charge (C)
m = 9.109e-31    # electron mass (kg)
B0 = 0.1         # magnetic field strength (T), along z
E0 = 1000.0      # electric field strength (V/m), along x

# Initial conditions
v0 = 1e6         # initial speed (m/s)
v = np.array([v0, 0.0, 0.0])  # initial velocity along x
r = np.array([0.0, 0.0, 0.0]) # initial position

# Derived quantities
omega_c = abs(q) * B0 / m     # cyclotron frequency
r_L = m * v0 / (abs(q) * B0)  # Larmor radius
T_c = 2 * np.pi / omega_c     # cyclotron period

print(f"Charged Particle Motion in E x B Fields")
print(f"=" * 45)
print(f"Cyclotron frequency: omega_c = {omega_c:.3e} rad/s")
print(f"Larmor radius: r_L = {r_L*1000:.4f} mm")
print(f"Cyclotron period: T_c = {T_c:.3e} s")
print(f"E x B drift velocity: v_d = {E0/B0:.1f} m/s")

# Time integration using Boris pusher
dt = T_c / 100   # 100 steps per cyclotron period
N_steps = 800     # ~8 cyclotron periods
t = np.zeros(N_steps)
pos = np.zeros((N_steps, 3))
vel = np.zeros((N_steps, 3))

pos[0] = r
vel[0] = v

B = np.array([0.0, 0.0, B0])

for n in range(N_steps - 1):
    t[n+1] = t[n] + dt
    E = np.array([E0, 0.0, 0.0])

# Boris pusher
    qdt_2m = q * dt / (2.0 * m)
    v_minus = vel[n] + qdt_2m * E
    t_vec = qdt_2m * B
    s_vec = 2.0 * t_vec / (1.0 + np.dot(t_vec, t_vec))
    v_prime = v_minus + np.cross(v_minus, t_vec)
    v_plus = v_minus + np.cross(v_prime, s_vec)
    vel[n+1] = v_plus + qdt_2m * E
    pos[n+1] = pos[n] + vel[n+1] * dt

# Plotting
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0f172a')
for ax in axes.flat:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

# Trajectory (x-y plane)
ax = axes[0, 0]
ax.plot(pos[:, 0]*1e3, pos[:, 1]*1e3, color='#f59e0b', linewidth=1.5, alpha=0.8)
ax.plot(pos[0, 0]*1e3, pos[0, 1]*1e3, 'o', color='#22d3ee', markersize=8, label='Start')
ax.set_title('Particle Trajectory (x-y plane)', color='#e2e8f0', fontsize=13, fontweight='bold')
ax.set_xlabel('x (mm)', color='#94a3b8')
ax.set_ylabel('y (mm)', color='#94a3b8')
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
ax.set_aspect('equal')

# Velocity components
ax = axes[0, 1]
ax.plot(t*1e9, vel[:, 0]/1e6, color='#f87171', linewidth=1.5, label='v_x')
ax.plot(t*1e9, vel[:, 1]/1e6, color='#22d3ee', linewidth=1.5, label='v_y')
ax.set_title('Velocity Components', color='#e2e8f0', fontsize=13, fontweight='bold')
ax.set_xlabel('Time (ns)', color='#94a3b8')
ax.set_ylabel('Velocity (10^6 m/s)', color='#94a3b8')
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Position vs time
ax = axes[1, 0]
ax.plot(t*1e9, pos[:, 0]*1e3, color='#f87171', linewidth=1.5, label='x(t)')
ax.plot(t*1e9, pos[:, 1]*1e3, color='#22d3ee', linewidth=1.5, label='y(t)')
ax.set_title('Position vs Time', color='#e2e8f0', fontsize=13, fontweight='bold')
ax.set_xlabel('Time (ns)', color='#94a3b8')
ax.set_ylabel('Position (mm)', color='#94a3b8')
ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Kinetic energy
KE = 0.5 * m * np.sum(vel**2, axis=1) / q  # in eV
ax = axes[1, 1]
ax.plot(t*1e9, KE, color='#a78bfa', linewidth=2)
ax.set_title('Kinetic Energy', color='#e2e8f0', fontsize=13, fontweight='bold')
ax.set_xlabel('Time (ns)', color='#94a3b8')
ax.set_ylabel('Energy (eV)', color='#94a3b8')

plt.suptitle('Charged Particle in E x B Fields (Boris Pusher)', color='#e2e8f0', fontsize=15, fontweight='bold', y=1.02)
plt.tight_layout()

Click Run to execute the Python code

First run will download Python environment (~15MB)

Key Equations of Plasma Physics

Fundamental Plasma Parameters

Debye Length:

$$\lambda_D = \sqrt{\frac{\varepsilon_0 k_B T_e}{n_e e^2}}$$

Plasma Frequency:

$$\omega_{pe} = \sqrt{\frac{n_e e^2}{\varepsilon_0 m_e}}$$

Cyclotron Frequency:

$$\omega_c = \frac{|q|B}{m}$$

Larmor Radius:

$$r_L = \frac{m v_\perp}{|q| B}$$

Vlasov Equation

$$\frac{\partial f}{\partial t} + \vec{v} \cdot \nabla f + \frac{q}{m}(\vec{E} + \vec{v} \times \vec{B}) \cdot \frac{\partial f}{\partial \vec{v}} = 0$$

Collisionless kinetic equation for the distribution function $ f(\vec{x}, \vec{v}, t) $

Langmuir Wave Dispersion

$$\omega^2 = \omega_{pe}^2 + 3k^2 v_{th}^2$$

Landau damping rate:

$$\gamma_L \approx -\sqrt{\frac{\pi}{8}} \frac{\omega_{pe}}{(k\lambda_D)^3} \exp\!\left(-\frac{1}{2k^2\lambda_D^2}\right)$$

Magnetohydrodynamics

MHD Momentum:

$$\rho\frac{d\vec{v}}{dt} = -\nabla p + \vec{J} \times \vec{B}$$

Alfven Speed:

$$v_A = \frac{B}{\sqrt{\mu_0 \rho}}$$

Magnetic Pressure:

$$p_B = \frac{B^2}{2\mu_0}$$

Plasma Beta:

$$\beta = \frac{2\mu_0 n k_B T}{B^2}$$

Guiding Center Drifts

E x B Drift

$$\vec{v}_E = \frac{\vec{E} \times \vec{B}}{B^2}$$

Grad-B Drift

$$\vec{v}_{\nabla B} = \frac{m v_\perp^2}{2qB^3} \vec{B} \times \nabla B$$

Curvature Drift

$$\vec{v}_R = \frac{m v_\parallel^2}{qB^2} \frac{\vec{R}_c \times \vec{B}}{R_c^2}$$