Computational Methods in Plasma Physics 2019

These are notes, lecture slides and code for use in PPPL’s Graduate Summer School, Computational Methods in Plasma Physics mini-course. All material is copyrighted by Ammar Hakim and released under the Creative Commons CC BY License.


Vast majority of plasma physics is contained in the Vlasov-Maxwell equations that describes the evolution of a particle distribution \(f_s(t,\mathbf{x},\mathbf{v})\) function in 6D phase-space. The particles move in electromagnetic fields that come from two sources: external coils and electrodes and the fields generated by the motion of the particles themselves. This makes the problem highly nonlinear: the fields tell the particles how to move, and the particles motion itself modifies the fields.

\[\frac{\partial f_s}{\partial t} + \nabla_\mathbf{x} \cdot (\mathbf{v}f_s) + \nabla_\mathbf{v} \cdot (\mathbf{a}_s f_s) = \left( \frac{\partial f_s}{\partial t} \right)_c\]

where \(\mathbf{a}_s = q_s/m_s(\mathbf{E}+\mathbf{v}\times\mathbf{B})\) is the acceleration due to Lorentz force on the particle and the right-hand side describes the effect of collisions. Of course, the electromagnetic fields are determined from Maxwell equations

\[\begin{split}\epsilon_0\mu_0 \frac{\partial \mathbf{E}}{\partial t} - \nabla_\mathbf{x} \times \mathbf{B} &= -\mu_0 \sum_s q_s \int_{-\infty}^{\infty} v f_s d\mathbf{v}^3 \\ \frac{\partial \mathbf{B}}{\partial t} + \nabla_\mathbf{x} \times \mathbf{E} &= 0 \\ \nabla_\mathbf{x}\cdot\mathbf{E} &= \frac{1}{\epsilon_0}\sum_s q_s \int_{-\infty}^{\infty} f_s d\mathbf{v}^3 \\ \nabla_\mathbf{x}\cdot\mathbf{B} &= 0.\end{split}\]

Suitable modifications are required to account for relativistic effects.

Having written the basic equations it may appear that the field of plasma physics is complete. However, this is far from the case. The Vlasov-Maxwell system is a formidable system of coupled, nonlinear equations and describes vast physics that spans an enormous range of temporal and spatial scales. Everything from electron oscillations to slow, resistive evolution of objects in near equilibrium is contained within them. Hence, over the decades many approximations to the Vlasov-Maxwell equations have been developed that are often more tractable, allowing one to obtain reasonable results in specific physical situations. Of course, the frontier in computational and theoretical plasma physics remains the complete kinetic understanding of plasma from the full Vlasov-Maxwell equations, and much of the recent research focus has been on solving them (or asymptotic approximations in strong magnetic fields) directly.

In this mini-course we will look at the key modern computational techniques to handle various problems in plasma physics. This field remains an active area of research, and in the recent years has undergone a renaissance of sorts. My goal will be give you a flavor of the methods and point out key literature in this field. In particular we will study:

  • Near first-principles simulations methods in which the Vlasov-Maxwell equation is solved directly. We will study two methods: the finite-difference-time-domain particle-in-cell method (FDTD PIC) (Lectures 1 and 2)
  • Solving fluid equations that are obtained from taking moments of the Vlasov-Maxwell equations and making a closure to approximate the moments not evolved by the fluid equations. Here we will study the modern approach based on the theory of hyperbolic PDEs and Riemann solvers. We will also look at the special requirements of solvers that are required to study fusion problems (Lecture 3).
  • Directly discretizing the Vlasov-Maxwell equations as a PDE in 6D. This is an emerging area of active research and has applications to study of turbulence in fusion machines and also exploring fundamental plasma physics in phase-space. (Lecture 4)

Through these topics we will look at the general applied mathematics literature for solvers for ordinary differential equations (ODEs) and partial differential equations (PDEs), with emphasis of those techniques that are useful in plasma physics.

Github repo for lectures and code

You can get all the tex files for the slides, source for this website and the code I ran in class by cloning the Github repo for this class.

To run the code you need to install the Gkeyll software. Please see the Gkeyll website for instructions on how to get the code and the postprocessing tools needed. In brief:

  • Install the Python (mini)conda package manager

  • Install the Gkeyll code using:

    conda install -c gkyl gkeyll
  • Install the postprocessing tool as:

    conda install -c gkyl postgkyl

Lecture 1: Introduction and Specialized Ordinary Differential Equations Solvers

PDF of Lecture 1 slides


See files code/lec1/sho-fwd-euler.lua and code/lec1/sho-mid-point.lua for forward Euler and mid-point schemes for harmonic oscillator problem. You can run the code from the directory it is kept as:

gkyl sho-mid-point.lua

To plot the trajectory do:

pgkyl -f sho-mid-point_ptclData.bp traj -e 90 -a 0.0 --fix-aspect

You can get help for options for the traj command by doing:

pgkyl traj --help

To plot the exact trajectory and the computed trajectory do:

pgkyl -f sho-mid-point_ptclData.bp -f sho-mid-point_exactData.bp traj -e 90 -a 0.0 --fix-aspect

You can save the animate to an mp4 file by passing the –save option to the traj command. This requires that you have the ffmpeg package installed.


Read above introductory summary for overall background. It is important to understand how to derive conservation and other properties of the Vlasov-Maxwell system. Much of modern computational plasma physics is focused on inventing schemes that preserve at least some of these properties of the continuous system. See any plasma text book or the following excerpt from [Juno2019] for proofs.

The concept of phase-space volume preserving and symplectic schemes can be more easily understood by looking at the example of a simple harmonic oscillator

\[\frac{d^2z}{dt^2} = -\omega^2 z\]

where \(\omega\) is the oscillation frequency.

To fully understand the physics behind these concepts one needs to understand the Lagrangian and Hamiltonian formulation of mechanics. For example, see text book of Goldstein or first volume of Landau and Lifshitz, Mechanics. An overview of Hamiltonian mechanics using noncanonical coordinates as applied to single particle motion is given in Section II of [CaryBrizard2009].

A good description of various ODE solvers and their properties is given in Chapter 2 of [DurranBook]. Also see for formulas of the Strong-Stability preserving RK methods and their stability regions.

Several ODE schemes have been designed to handle stiff sources and in particular, diffusion terms arising from discretization of diffusion equations. See [Abdulle2013] and also [Meyer2013] for description of these schemes. In particular, the scheme by Meyer at al is to be preferred to it superior stability properties.

The ODE solvers described above are low order, that is second or third order. Some recent work attempts to construct very high order schemes (10-15th order!) that essentially makes the issues of conservation and other numerical errors mostly moot. For example, see [ReinSpiegel] for a 15th order scheme for use in gravitational N-body simulations. Such very high-order schemes have not found use in plasma-physics yet, mainly as the Maxwell solvers used in PIC codes are mostly second-order anyway. However, it is possible that these very high-order methods are useful in orbit codes.

Lecture 2: The Boris algorithm and FDTD and FV schemes for Maxwell equations

PDF of Lecture 2 slides. Solution to the problem of finding \(\mathbf{A}\) if \(\mathbf{A} = \mathbf{R} + \mathbf{A}\times\mathbf{B}\) is here.


See files code/lec2/lorentz-boris.lua

You can play with this file to do various static or time-dependent electromagnetic fields. For example, motion in a constant magnetic field, in a field with a gradient, and in a driven system. See field specification in this write up for both non-resonant and resonant drivers.


Particle-in-cell methods are based on pushing macro-particles. These represent the motion of characteristics in phase-space, along which the distribution function is conserved. The macro-particle equations-of-motion are

\[\begin{split}\frac{d\mathbf{x}}{dt} &= \mathbf{v} \\ \frac{d\mathbf{v}}{dt} &= \frac{q}{m}(\mathbf{E} + \mathbf{v}\times\mathbf{B})\end{split}\]

The most widely used method to solve this system of ODEs is the Boris algorithm. See this excerpt from Birdsall and Langdon book for details on how to implement this efficiently.

The Boris algorithm is surprisingly good: it is a second-order, time-centered method that conserves phase-space volume. However, the error in phase-velocity (that is there is an error in time-period of orbits) accumulates linearly, as we saw for the harmonic oscillator. See [Qin2013] for proofs that the Boris algorithm is not symplectic but conserves phase-space volume.

The relativistic Boris algorithm does not compute the correct \(\mathbf{E}\times\mathbf{B}\) velocity. This can be corrected for and still maintain the volume-preserving property and was done in [HigueraCary2017].

The Yee-cell preserves the underlying geometric structure of Maxwell equations, and ensures that the divergence relations are maintained in the case of vacuum fields. In a plasma, however, current deposition needs to be done carefully to ensure current continuity is satisfied. See [Esirkepov2001], for example.

For extension of standard FDTD method to complex geometries, see, for example [Nieter2009] and other references. Recent research has focused on developing finite-element based PIC codes (that maintain geometric structure of Maxwell equations), but these are usually very expensive to run and very complex to develop.

Sometimes finite-volume schemes are also used to solve Maxwell equations. These may have some advantages and disadvantages compared to standard FDTD schemes. For example, FV usually do not conserve energy and find it hard to satisfy divergence relations. For a comparison of FV and FDTD methods see this page.

A comprehensive review of structure preserving algorithms for use in plasma physics is provided by [Morrison2017]. It has numerous references to the literature and should be consulted to develop a detailed understanding of such schemes.

Lecture 3: The cutting-edge of computational plasma physics research

I showed slides from Eric Shi’s thesis defense. Please see his thesis for full details on the gyrokinetic algorithm, validation with LAPD and also first results from simulating the NSTX scrap-off-layer plasma.

The Big Picture

Originally, I had intended to spend time talking about fluid solvers in this lecture. However, based on the extensive questions in the last lecture I decided the time would be better spent in looking at a few big-picture issues. We will return to numerics in Lecture 4.

In particular:

  • What are the cutting-edge research questions in computational plasma physics?
  • What is the relationship between modern numerical methods and experiments and observations? (That is, why care about this stuff in the first place? Can simulations predict rather than postdict?)
  • How to incorporate “real-world” effects into simulations? (For example, boundary conditions, atomic physics, etc)

One can look at computational physics in two ways: as an end in itself, and as a tool for applications. Both of these are important!

As an end in itself:

  • The first sits between applied mathematics and theoretical physics. The goal is to design efficient numerical methods to solve equations from theoretical physics.
  • The goal here is the numerical method itself: what are its properties? Does it faithfully represent the underlying physics? Does it run efficiently on modern computers? Research into modern numerical methods (including structure preserving methods) fall into this category.
  • Usually, besides the fun of solving complex equations (and writing code), the goal is to gain deeper understanding of underlying physics. Some theoretical questions can only be answered with computer simulations.
  • This is a perfectly legitimate research area even if no connection to experiments is made, but only satisfies the curiosity of the researchers and helps one gain a better understanding of the physics.

As a tool for applications:

  • The second is to look at the computational physics as providing tools to understand/design experiments or observations.
  • Note that a large number of routine calculations are needed to build modern experiments (heat-transfer, structural analysis, basic fluid mechanics, equilibrium and stability calculations, etc). Such routine calculations are no longer cutting edge research topics.

At the intersection of cutting-edge computational physics and modern plasma physics is a set of Billion Dollar Questions. (In general, one should not put currency values to such things).

These Billion Dollar Questions need huge investments in experimental and observational programs as well as the very latest in computational physics research.

Space Physics Examples: Parker Solar Probe

  • Paker Solar Probe. “The primary science goals for the mission are to trace how energy and heat move through the solar corona and to explore what accelerates the solar wind as well as solar energetic particles.”
  • The Probe will collect detailed measurements of electric and magnetic fields as well as detailed distribution functions of particles.
  • The solar wind plasma is nearly collisionless. It is likely that a proper understanding of kinetic physics (at the level of the Vlasov-Maxwell equations) will be needed to fully understand the physical processes.
  • Cutting-edge simulations will be critcial to this. Serious research into numerics of Vlasov-Maxwell needs to be done and very large simulations need to be run.

Many other missions are active and planned: BepiColombo to Mercury; Juno to Jupiter.

  • Much of the deep understanding of plasma processes in solar system planets (magnetospheres, ionosphere) can only be gained from detailed modeling: global kinetic modelling is likely impossible. How to incorporate some kinetic effects into fluid models?

Fusion Physics Examples: Building a working thermonuclear fusion reactor.

  • The Iter project aims to build the world’s largest tokamak, a “magnetic bottle” to contain super-hot plasma and heat it to ignition temperatures.

There are other major fusion efforts around the world:

There are major unsolved problems in the basic physics of fusion machines. Most of these can only be answered by large-scale computing and much of the numerical tools have not yet been fully developed.

The Scientific Discovery through Advanced Computing program in fusion has large projects that address the very serious Billion Dollar Question: will controlled fusion be eventually possible?

  • The numerics research here is focussed on gyrokinetic and even full kinetic understanding of fundamental turbulence and transport processes in the tokamak. These equations are very difficult to solve!
  • Disruptions are dangerous processes that can “kill” certain fusion machine: large-scale MHD simulations are needed. Significant new research is being done in new numerical methods and application of existing MHD codes to such problems.
  • Runaway electrons (relativistic high-energy electron beams) can drill holes in fusion machines. See SCREAM project and special PPCF issue.
  • Very serious! Will need huge kinetic calculations. Also, the formulation of self-consistent coupling betwen the runaway electrons and MHD is not complete. See review by [Boozer2015].

These are only selection of problems I am directly familiar with. I hope it gives you a flavor and understanding why computational plasma physics is such a serious and important field!

Lecture 4: Discontinuous Galerkin Schemes

PDF of Lecture 4 slides.


The discontinuous Galerkin method is a type of finite-element scheme in which the solution across cell boundaries can be discontinuous. The method was originally invented for elliptic equations by Nitsche in 1971 (paper is in German). However, the key paper on application to a hyperbolic PDE was written by Reed and Hill in 1973. The latter paper has more than 2100 citations. A later key paper on extension of the method to nonlinear systems of equations was by Cockburn and Shu (JCP 41, 199-224 1998).

The key ingredient of DG is Galerkin minimization on a finite-dimensional discontinuous space. This space is usually piecewise continuous polynomials, but could be other types of functions too.

I introduced the idea of weak-equality, which is central to a proper understanding of the DG scheme. This weak-equality concept, and idea of recovery built from it, can be used to construct schemes for both advection and diffusion equation that is higher-order accurate than standard DG in smooth regions. Some form of limited recovery is needed when the solution is not smooth.

One can solve the full Vlasov-Maxwell equations using DG. See [Juno2019] for example. The scheme described there conserves energy exactly, but does not conserve momentum. This is a typical feature of solvers for Vlasov equations: one can either conserve momentum or energy but not both. The construction of a momentum and energy conserving scheme remains an ongoing research problem.


[Juno2019](1, 2) Juno, J., Hakim, A., TenBarge, J., Shi, E., Dorland, W. (2018). “Discontinuous Galerkin algorithms for fully kinetic plasmas”, Journal of Computational Physics, 353, 110–147.
[CaryBrizard2009]Cary, J. R., & Brizard, A. J. (2009). “Hamiltonian theory of guiding-center motion”. Reviews of Modern Physics, 81 (2), 693–738.
[DurranBook]Dale E. Durran, “Numerical Methods for Fluid Dynamics”, Springer. Second Edition.
[Abdulle2013]Abdulle, A., & Vilmart, G. (2013). “PIROCK: A swiss-knife partitioned implicit–explicit orthogonal Runge–Kutta Chebyshev integrator for stiff diffusion–advection–reaction problems with or without noise”. Journal of Computational Physics, 242 (C), 869–888.
[Meyer2013]Meyer, C. D., Balsara, D. S., & Aslam, T. D. (2014). “A stabilized Runge–Kutta–Legendre method for explicit super-time-stepping of parabolic and mixed equations”. Journal of Computational Physics, 257 (PA), 594–626.
[ReinSpiegel]Rein, H., & Spiegel, D. S. (2014). ias15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits. Monthly Notices of the Royal Astronomical Society, 446(2), 1424–1437.
[Qin2013]Qin, H., Zhang, S., Xiao, J., Liu, J., Sun, Y., & Tang, W. M. (2013). “Why is Boris algorithm so good?” Physics of Plasmas, 20 (8), 084503–5.
[HigueraCary2017]Higuera, A. V., & Cary, J. R. (2017). “Structure-preserving second-order integration of relativistic charged particle trajectories in electromagnetic fields”. Physics of Plasmas, 24 (5), 052104–7.
[Esirkepov2001]Esirkepov, T. Z. (2001). “Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor”, Computer Physics Communications, 135, 144–153.
[Nieter2009]Nieter, C., Cary, J. R., Werner, G. R., Smithe, D. N., & Stoltz, P. H. (2009). “Application of Dey–Mittra conformal boundary algorithm to 3D electromagnetic modeling”. Journal of Computational Physics, 228 (21), 7902–7916.
[Morrison2017]Morrison, P. J. (2017). Structure and structure-preserving algorithms for plasma physics. Physics of Plasmas, 24 (5), 055502–21.
[Boozer2015]Boozer, A. H. (2015). “Theory of runaway electrons in ITER: Equations, important parameters, and implications for mitigation”. Physics of Plasmas, 22 (3), 032504–18.