BOSTON UNIVERSITY
Department of Astronomy
Center for Space Physics

Computational Plasma Physics

Particle-in-cell, hybrid, fluid, and electromagnetic FDTD codes running on thousands of CPU cores — combined with rigorous analytic theory — are the engine of our research.

Our Simulation Philosophy

Plasma physics problems in space are often too complex for purely analytic approaches, and too rich in kinetic physics for simple fluid models. Our group occupies a productive niche: we develop and deploy codes that capture kinetic, fluid, and wave physics at the scales where the interesting phenomena occur and measurements are made. We use four main classes of simulation, described below, along with strong analytic theory to enable us to understand complex physical processes.

Our codes are run on national high-performance computing (HPC) facilities. Graduate students in the group typically spend their first year becoming comfortable with the codes and HPC environment, and spend subsequent years extending the codes and using them to solve new problems. Many of our simulation codes were built from scratch within the group, giving students genuine expertise in both plasma physics and scientific software.

10⁴+
CPU cores per large run
10³–10¹⁰
Particles in our largest simulations
2D & 3D
Full three-dimensional codes
30+
Years of PIC code development

The scale of these runs is enormous. Our recent high-latitude electrojet simulation (Oppenheim, Dimant, Koontaweepunya, Green & Evans, GRL 2026) followed the plasma on a grid of 32,768 × 512 × 512 — roughly 8.6 billion cells — carrying more than 40 billion computational particles (about 8.7 billion electrons and 34 billion ions). At 0.15 m resolution the box stretches nearly 5 km along the geomagnetic field yet spans only ~77 m across it, so a single run resolves everything from centimeter-scale kinetic waves to kilometer-long structures at once. Pushing all of that through tens of milliseconds of plasma evolution — many thousands of time steps — keeps thousands of CPU cores busy for days on national supercomputers, and produces datasets measured in terabytes.

Ion-density turbulence in a vertical plane spanning the electrojet. Ion-density turbulence in horizontal planes at a sequence of altitudes.

Output from that simulation: ion density in a vertical plane across the electrojet (top) and in horizontal planes at a sequence of altitudes (bottom), as Farley–Buneman turbulence grows from noise into fully developed waves. Click either movie to play. From Oppenheim et al., GRL 2026.

Our Four Simulation Approaches

■ Particle-in-Cell (PIC)

PIC codes represent plasma electrons and ions as individual computational particles, moving under self-consistent electromagnetic fields. At each time step: (1) particle charges and currents are deposited onto a spatial grid; (2) Poisson’s equation (electrostatic) or Maxwell’s equations (electromagnetic) are solved on the grid; (3) particles are pushed with the Lorentz force. This first-principles approach captures wave–particle interactions, resonant heating, Landau damping, and nonlinear saturation mechanisms invisible to fluid models. Our 2D and 3D electrostatic PIC codes, developed in-house over three decades, have produced some of the highest-resolution kinetic simulations of ionospheric turbulence ever achieved. We also deploy electromagnetic PIC for chromospheric and auroral problems.

■ Hybrid Simulations

In many problems, ions can be treated as fluid particles while electrons require a kinetic treatment — or vice versa. Our hybrid codes exploit this separation: typically, ions are represented as macroparticles while electrons are treated as a massless, charge-neutralizing fluid (or vice versa). This approach dramatically reduces computation cost relative to full PIC while retaining the kinetic physics that matters most. Hybrid simulations in our group have been particularly valuable for studying coupled Farley–Buneman / gradient-drift instabilities in the equatorial E region (Young, Oppenheim & Dimant 2017), where the separation of scales between ions and electrons can be exploited.

■ Multi-Fluid Simulations

For large-scale problems and multi-species plasmas — such as the partially ionized solar chromosphere containing electrons, protons, neutral hydrogen, helium, and heavier species — we use multi-fluid codes that track each species with its own fluid equations. These codes sacrifice single-particle kinetics in exchange for the ability to simulate large spatial domains and long time scales at reasonable computational cost. Our multi-fluid chromospheric codes (Evans et al. 2023, 2025) were the first to simulate the thermal Farley–Buneman instability in the solar chromosphere over physically realistic domain sizes, capturing the turbulent heating rate and its dependence on solar conditions.

■ Electromagnetic FDTD

Finite-difference time-domain (FDTD) methods solve Maxwell’s equations directly on a spatial grid, advancing electric and magnetic fields forward in time. Our group applies electromagnetic FDTD codes to model how radio waves propagate, scatter, and refract through structured ionospheric plasma — including the turbulent irregularities produced by the instabilities we simulate with PIC. The open-source radio wave propagation code released by Green, Longley, Oppenheim & Young (2025, Frontiers in Astronomy and Space Sciences) implements this approach.

Relative perturbed plasma density from a 2-D PIC/hybrid simulation of coupled Farley-Buneman/gradient-drift turbulence in the equatorial E region.
PIC/hybrid simulation of coupled Farley–Buneman / gradient-drift turbulence. Time evolution of the relative perturbed density (δn/n₀) in a 2-D PIC/hybrid simulation of the equatorial E region: ions are evolved kinetically via a particle-in-cell (PIC) method while electrons are treated as an inertialess fluid. Gradient-drift waves grow along the background density gradient, and Farley–Buneman turbulence develops in the central density trough once the local electric field exceeds threshold. Click to play. From Young, Oppenheim & Dimant (2017).

National Supercomputing Resources

Our group holds competitive allocations on national HPC systems through NSF ACCESS (formerly XSEDE) and uses BU’s on-campus Shared Computing Cluster (SCC) for development and moderate-scale production runs.

NSF ACCESS — Stampede3, Anvil

Large allocation for production-scale PIC and multi-fluid simulations. Our kinetic electrojet simulations routinely use thousands of cores for multi-day runs.

BU Shared Computing Cluster (SCC)

On-campus cluster for development, parameter surveys, testing, and course computation. Enables rapid iteration before committing to national facility runs.

GPU Computing

We are porting our innermost PIC loop kernels to CUDA/HIP for GPU-accelerated nodes, leveraging the high arithmetic throughput of modern accelerated architectures.

Some of Our Techniques

Selected Simulation Publications

Computational Plasma Physics

This work involves writing codes that run on thousands of cores and studying plasma physics from first principles, spanning both solar and space physics.

Research in the group runs through the BU Astronomy PhD program: bu.edu/astronomy/graduate.

← Back to Home