A comprehensive, educational Python toolkit covering every building block of a Particle-In-Cell (PIC) plasma simulation — from single-particle pushers to complete electrostatic and electromagnetic PIC codes, advanced diagnostics, and collisional physics, all with animated visualisations.

22 self-contained scripts · NumPy + Matplotlib only · No compilation required

• Overview — The PIC Method
• Physics Background
• Repository Structure
• Module Descriptions
  • Particle Pushers
  • Field Solvers
  • Particle–Grid Coupling
  • Complete PIC Codes
  • Plasma Physics Demos
  • Diagnostics & Analysis
  • Advanced Techniques
• Requirements
• Quick Start
• How the Pieces Fit Together
• Possible Improvements
• References
• License

The Particle-In-Cell (PIC) method is the workhorse of computational plasma physics. It represents the plasma as a collection of macro-particles and evolves them self-consistently with electromagnetic fields on a grid.

The PIC loop (every time step):

┌─────────────────────────────────────────────────────────┐
│                     PIC CYCLE                           │
│                                                         │
│  1. CHARGE / CURRENT DEPOSITION                         │
│     particles → grid  (weighting / shape functions)     │
│           ↓                                             │
│  2. FIELD SOLVE                                         │
│     Poisson (ES) or Maxwell / Yee FDTD (EM)            │
│           ↓                                             │
│  3. FIELD INTERPOLATION                                 │
│     grid → particles  (gather)                          │
│           ↓                                             │
│  4. PARTICLE PUSH                                       │
│     Boris / Vay / Higuera-Cary  (leapfrog in time)     │
│           ↓                                             │
│  5. BOUNDARY CONDITIONS                                 │
│     periodic / reflecting / absorbing                   │
│           ↓                                             │
│  6. (OPTIONAL) COLLISIONS                               │
│     Monte Carlo Collisions (MCC)                        │
│                                                         │
│  repeat ────────────────────────────────────────────→    │
└─────────────────────────────────────────────────────────┘

Each script in this repo implements one or more of these steps. The full PIC codes tie them all together.

Non-relativistic:

$$m \frac{d\mathbf{v}}{dt} = q\left(\mathbf{E} + \mathbf{v} \times \mathbf{B}\right)$$

Relativistic (using proper velocity $\mathbf{u} = \gamma \mathbf{v}$):

$$m \frac{d\mathbf{u}}{dt} = q\left(\mathbf{E} + \frac{\mathbf{u}}{\gamma} \times \mathbf{B}\right), \qquad \gamma = \sqrt{1 + \frac{|\mathbf{u}|^2}{c^2}}$$

$$\nabla^2 \varphi = -\frac{\rho}{\varepsilon_0}, \qquad \mathbf{E} = -\nabla \varphi$$

$$\frac{\partial \mathbf{B}}{\partial t} = -\nabla \times \mathbf{E}$$

$$\varepsilon_0 \frac{\partial \mathbf{E}}{\partial t} = \frac{1}{\mu_0}\nabla \times \mathbf{B} - \mathbf{J}$$

Boris_Algorithm/
│
│  ── PARTICLE PUSHERS ──────────────────────────────────
├── Non-relativistic-Boris.py    # Classical Boris pusher, B-field sweep
├── Relativistic-Boris.py        # Relativistic Boris (proper velocity u = γv)
├── Vay-Pusher.py                # Vay pusher — Boris vs Vay comparison
├── Higuera-Cary-Pusher.py       # 3-way pusher comparison (Boris/Vay/HC)
│
│  ── FIELD SOLVERS ─────────────────────────────────────
├── Poisson-Solver-1D2D.py       # 1D tridiagonal & 2D FFT Poisson solvers
├── Yee-FDTD-2D.py               # 2D TM-mode FDTD Maxwell solver
│
│  ── PARTICLE–GRID COUPLING ───────────────────────────
├── Particle-Weighting.py        # Shape functions (NGP/CIC/TSC), Esirkepov J
├── Binomial-Filter.py           # Smoothing filters & transfer functions
│
│  ── COMPLETE PIC CODES ───────────────────────────────
├── PIC-1D-TwoStream.py          # 1D electrostatic PIC (two-stream)
├── EM-PIC-1D.py                 # 1D electromagnetic PIC (laser-plasma)
├── Weibel-Instability.py        # 2D EM PIC (Weibel / filamentation)
│
│  ── PLASMA PHYSICS DEMOS ─────────────────────────────
├── Plasma-Oscillation.py        # Cold plasma wave, ω_pe measurement
├── Landau-Damping.py            # ES PIC, Landau damping rate measurement
├── Ion-Acoustic-Wave.py         # Two-species ES PIC (electrons + ions)
├── Magnetic-Mirror.py           # Converging B-field, adiabatic invariant μ
├── Debye-Sheath.py              # Plasma–wall sheath, Bohm criterion
├── Vlasov-1D1V.py               # Semi-Lagrangian Vlasov–Poisson solver
│
│  ── DIAGNOSTICS & ANALYSIS ───────────────────────────
├── Dispersion-Relation.py       # Analytical ω(k) + PIC FFT validation
├── Convergence-Analysis.py      # Δt / Δx / PPC convergence studies
├── Energy-Conservation.py       # Energy, momentum, Gauss's law checks
│
│  ── ADVANCED TECHNIQUES ──────────────────────────────
├── Monte-Carlo-Collisions.py    # PIC-MCC: elastic, ionisation, CX
├── Particle-Merging.py          # Adaptive particle merging & splitting
│
└── README.md

File: Non-relativistic-Boris.py

The classic Boris algorithm in SI units for an electron in a uniform magnetic field ($\mathbf{E} = 0$).

Boris algorithm steps:

File: Relativistic-Boris.py

Same Boris splitting, but operates on the proper velocity $\mathbf{u} = \gamma \mathbf{v}$. The rotation vector gains a $\gamma^{-}$ factor:

$$\mathbf{t} = \frac{q \mathbf{B},\Delta t}{2 m \gamma^{-}}, \qquad \gamma^{-} = \sqrt{1 + |\mathbf{u}^{-}|^2/c^2}$$

File: Vay-Pusher.py

The Vay pusher (2008) is an alternative to Boris for relativistic regimes, especially when $|\mathbf{E}| \sim c|\mathbf{B}|$. It updates the 4-velocity implicitly.

File: Higuera-Cary-Pusher.py

Three-way comparison of Boris, Vay, and Higuera-Cary pushers in crossed $\mathbf{E} \times \mathbf{B}$ fields.

File: Poisson-Solver-1D2D.py

Solves $\nabla^2\varphi = -\rho/\varepsilon_0$ and computes $\mathbf{E} = -\nabla\varphi$.

File: Yee-FDTD-2D.py

Full 2D Finite-Difference Time-Domain solver for Maxwell's equations on a staggered Yee grid (TM mode: $E_z$, $B_x$, $B_y$).

File: Particle-Weighting.py

Also includes Esirkepov current deposition (charge-conserving $\nabla\cdot\mathbf{J} + \partial\rho/\partial t = 0$).

File: Binomial-Filter.py

Digital smoothing filters for PIC current/charge densities.

Includes transfer function analysis in Fourier space.

File: PIC-1D-TwoStream.py

$$\text{CIC deposit} ;\xrightarrow{\rho}; \text{FFT Poisson} ;\xrightarrow{E}; \text{CIC gather} ;\xrightarrow{F}; \text{Leapfrog push}$$

File: EM-PIC-1D.py

$$\text{CIC current} ;\xrightarrow{J}; \text{Yee FDTD} ;\xrightarrow{E, B}; \text{CIC gather} ;\xrightarrow{F}; \text{Boris push}$$

File: Weibel-Instability.py

Two counter-streaming electron beams in 2D generate the Weibel (filamentation) instability: spontaneous growth of magnetic fields perpendicular to the beam direction.

File: Plasma-Oscillation.py

Cold plasma with sinusoidal density perturbation. Measures the plasma frequency $\omega_{pe}$ from field energy oscillation via FFT.

File: Landau-Damping.py

100 000-particle ES PIC with Maxwellian distribution + low-amplitude sinusoidal perturbation. Measures the collisionless damping rate and compares to the analytical Landau result:

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

File: Ion-Acoustic-Wave.py

Two-species (electron + ion) ES PIC with $m_i/m_e = 100$. Ion density perturbation excites an ion acoustic wave; measures IAW frequency:

$$\omega_{IAW} = \frac{k,c_s}{\sqrt{1 + k^2\lambda_D^2}}$$

File: Magnetic-Mirror.py

Single particle in a parabolic mirror field $B_z(z) = B_0(1 + z^2/L_m^2)$ with divergence-free $B_r$. Tests 5 pitch angles to demonstrate:

• Magnetic mirroring (reflection of trapped particles)
• Loss cone (particles escaping through the mirror throat)
• Conservation of the adiabatic invariant $\mu = m v_\perp^2 / (2B)$

File: Debye-Sheath.py

Plasma bounded by two conducting walls. Both electrons and ions are kinetic with absorbing wall BCs.

File: Vlasov-1D1V.py

Grid-based (not PIC!) solver for the Vlasov equation on a $(x, v)$ phase-space grid using Strang splitting and semi-Lagrangian advection.

$$\frac{\partial f}{\partial t} + v\frac{\partial f}{\partial x} + \frac{q E}{m}\frac{\partial f}{\partial v} = 0$$

File: Dispersion-Relation.py

Computes analytical dispersion relations and validates them against PIC simulation data in $(k, \omega)$ space.

PIC validation: runs a broadband 1D ES PIC, FFTs $E(x,t) \to |E(k,\omega)|^2$ and overlays the analytical curve.

File: Convergence-Analysis.py

Systematic sweep of PIC parameters on a cold plasma oscillation benchmark:

File: Energy-Conservation.py

Standalone diagnostic module monitoring:

File: Monte-Carlo-Collisions.py

Implements binary collision models for gas-discharge / low-temperature plasmas:

Uses the null-collision method for efficient sampling. Demo runs a 1D capacitively-coupled plasma (CCP) discharge with RF voltage drive.

Diagnostics: EEDF, IEDF, particle populations, cumulative collision statistics.

File: Particle-Merging.py

Adaptive macro-particle count control:

Diagnostics: PPC distribution heatmap, merge/split event counts, momentum conservation, weight histogram.

pip install numpy matplotlib scipy

Each script is self-contained — just run it:

# ── Particle Pushers ──
python Non-relativistic-Boris.py
python Relativistic-Boris.py
python Vay-Pusher.py
python Higuera-Cary-Pusher.py

# ── Field Solvers ──
python Poisson-Solver-1D2D.py
python Yee-FDTD-2D.py

# ── Particle–Grid Operations ──
python Particle-Weighting.py
python Binomial-Filter.py

# ── Full PIC Simulations ──
python PIC-1D-TwoStream.py
python EM-PIC-1D.py
python Weibel-Instability.py

# ── Plasma Physics Demos ──
python Plasma-Oscillation.py
python Landau-Damping.py
python Ion-Acoustic-Wave.py
python Magnetic-Mirror.py
python Debye-Sheath.py
python Vlasov-1D1V.py

# ── Diagnostics ──
python Dispersion-Relation.py
python Convergence-Analysis.py
python Energy-Conservation.py

# ── Advanced ──
python Monte-Carlo-Collisions.py
python Particle-Merging.py

Each script opens an animated Matplotlib window with multi-panel diagnostics (or static multi-panel plots for analysis scripts).

                         PIC SIMULATION TOOLKIT
                         ══════════════════════

  ┌────────────────────────────────────────────────────────────────┐
  │                                                                │
  │  ┌──────────────────────────────────────────────────────────┐  │
  │  │         Particle-Weighting.py  ·  Binomial-Filter.py     │  │
  │  │         (deposit ρ, J to grid · smooth)                  │  │
  │  └──────────────────────┬───────────────────────────────────┘  │
  │                         ↓                                      │
  │  ┌──────────────────────────────────────────────────────────┐  │
  │  │    Poisson-Solver-1D2D.py  (ES)   ·   Yee-FDTD-2D.py   │  │
  │  │    (solve for E)                  (solve for E, B)       │  │
  │  └──────────────────────┬───────────────────────────────────┘  │
  │                         ↓                                      │
  │  ┌──────────────────────────────────────────────────────────┐  │
  │  │         Particle-Weighting.py                            │  │
  │  │         (interpolate E, B → particle positions)          │  │
  │  └──────────────────────┬───────────────────────────────────┘  │
  │                         ↓                                      │
  │  ┌──────────────────────────────────────────────────────────┐  │
  │  │     Non-relativistic-Boris.py  ·  Relativistic-Boris.py │  │
  │  │     Vay-Pusher.py  ·  Higuera-Cary-Pusher.py            │  │
  │  │     (advance v, x)                                       │  │
  │  └──────────────────────┬───────────────────────────────────┘  │
  │                         ↓                                      │
  │  ┌──────────────────────────────────────────────────────────┐  │
  │  │     Monte-Carlo-Collisions.py  (optional collisions)     │  │
  │  │     Particle-Merging.py        (adaptive particles)      │  │
  │  └──────────────────────┬───────────────────────────────────┘  │
  │                         ↓                                      │
  │                     repeat                                     │
  │                                                                │
  ├────────────────────────────────────────────────────────────────┤
  │  COMPLETE CODES:                                               │
  │   PIC-1D-TwoStream.py      (1D ES)                            │
  │   EM-PIC-1D.py             (1D EM)                            │
  │   Weibel-Instability.py    (2D EM)                            │
  │   Debye-Sheath.py          (1D ES, bounded)                   │
  │   Ion-Acoustic-Wave.py     (1D ES, two-species)               │
  │   Monte-Carlo-Collisions.py (1D ES + MCC)                     │
  ├────────────────────────────────────────────────────────────────┤
  │  ALTERNATIVE METHODS:                                          │
  │   Vlasov-1D1V.py      (grid-based Vlasov, no particles)       │
  ├────────────────────────────────────────────────────────────────┤
  │  PHYSICS BENCHMARKS:                                           │
  │   Plasma-Oscillation.py    ω_pe measurement                   │
  │   Landau-Damping.py        Collisionless damping               │
  │   Magnetic-Mirror.py       Adiabatic invariant μ               │
  ├────────────────────────────────────────────────────────────────┤
  │  DIAGNOSTICS:                                                  │
  │   Dispersion-Relation.py   ω(k) theory vs PIC                 │
  │   Convergence-Analysis.py  Δt, Δx, PPC scaling                │
  │   Energy-Conservation.py   Conservation checks                 │
  └────────────────────────────────────────────────────────────────┘

1. CLI with argparse — Accept dt, N_particles, B, E, etc. from the command line.
2. OOP refactor — Extract BorisPusher, PoissonSolver, YeeSolver classes; share via composition.
3. Unit tests — Validate Larmor radius, energy conservation, Gauss's law, CFL stability.
4. Type hints & docstrings — NumPy-style documentation throughout.

5. Higher-order FDTD — Replace 2nd-order Yee with 4th-order spatial stencil.
6. PML boundaries — Perfectly Matched Layer instead of Mur ABC for better absorption.
7. Implicit PIC — Direct-implicit or energy-conserving implicit methods (removes $\omega_{pe} \Delta t < 2$ constraint).
8. Relativistic Weibel — Extend Weibel-Instability.py with relativistic Boris push.
9. 3D PIC — Full 3D electromagnetic PIC (computationally intensive but educational).
10. Electromagnetic dispersion — Add R-wave, L-wave, Bernstein mode dispersion curves.

11. Save animation — Export .mp4 / .gif via FFMpegWriter or PillowWriter.
12. Interactive controls — Matplotlib sliders or Plotly/Dash web app.
13. 3D phase-space — $(x, v_x, v_y)$ scatter plots for orbit topology.

14. Numba JIT — @numba.njit on inner loops for 10–100× speedup.
15. GPU acceleration — CuPy / JAX for large-scale PIC.
16. MPI decomposition — Domain decomposition for distributed-memory HPC.

1. J. P. Boris, Relativistic plasma simulation — optimization of a hybrid code, Proc. 4th Conf. Num. Sim. Plasmas (1970).
2. C. K. Birdsall & A. B. Langdon, Plasma Physics via Computer Simulation, Taylor & Francis, 2004.
3. H. Qin et al., "Why is Boris algorithm so good?", Physics of Plasmas 20, 084503 (2013). doi:10.1063/1.4818428
4. J.-L. Vay, "Simulation of beams or plasmas crossing at relativistic velocity", Physics of Plasmas 15, 056701 (2008).
5. K. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media", IEEE Trans. Antennas Propag. 14 (3), 302–307 (1966).
6. T. Zh. Esirkepov, "Exact charge conservation scheme for Particle-in-Cell simulation", Comput. Phys. Commun. 135, 144–153 (2001).
7. R. W. Hockney & J. W. Eastwood, Computer Simulation Using Particles, CRC Press, 1988.
8. A. B. Langdon, "Effects of the spatial grid in simulation plasmas", J. Comput. Phys. 6, 247–267 (1970).
9. L. D. Landau, "On the vibrations of the electronic plasma", J. Phys. (USSR) 10, 25 (1946).
10. E. S. Weibel, "Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution", Phys. Rev. Lett. 2, 83 (1959).
11. V. Vahedi & M. Surendra, "A Monte Carlo collision model for the PIC method", Comput. Phys. Commun. 87, 179–198 (1995).

This project is provided for educational purposes. Feel free to use, modify, and distribute.