Computational Methods in Plasma Physics 2021

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.

Note

If you attended my computational physics lecture at the SULI 2022 summer program and are looking for the slides, please click here. For more detailed slides and notes, see below. Also see my notes and slides from a previous class I taught.

Background

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: (i) external coils and electrodes, and (ii) 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 particle 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.

For a detailed course on plasma kinetic theory please see the web-page of the Kinetic Theory Course taught by Dellar, Schekochihin, and Fouvry. In particular, see Schekochihin’s notes. The statistical mechanics of self-gravitating system is rather different from plasmas as gravitational “charge” (i.e. mass) comes only in a single “sign” (i.e. positive). See Fouvry’s notes for this important case of Vlasov-Poisson equations.

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 using the finite-difference-time-domain particle-in-cell method (FDTD PIC).
  • 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.
  • 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.

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

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.

Lecture 1: Introduction and the Big-Picture

PDF of Lecture 1 slides

Summary

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 Big Picture

To understand why computational methods have become central in plasma physics we need to look at a few big-picture issues.

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.

Fusion Physics Examples: Building a working thermonuclear fusion reactor.

Nuclear reactors (both fission and fusion) are the only power-source (besides fossil-fules) that can supply uninterupped and cheap power. No working or prototype fusion power plant exists as of today.

  • 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 focused 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 between the runaway electrons and MHD is not complete. See review by [Boozer2015].

Space Physics Examples

  • 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 critical 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 modeling is likely impossible. How to incorporate some kinetic effects into fluid models?

Extreme Plasma Astrophysics

An emerging field of plasma physics is the study of the plasma environment around compact astrophysical objects (neutron stars, black-holes). See for example: http://kyleparfrey.org/

This is an extremely challenging field: need to incorporate General Relativistic effects into plasma equations. Coupling between gravitational effects and plasma leads to extremely energetic events. releasing huge amounts of energy.

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 2: Special ODE Integrators, FDTD Scheme

PDF of Lecture 2 slides

Designing ODE solvers

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.

Particles in an electromagnetic field, FDTD methods

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: FDTD Scheme for Maxwell equations, hyperbolic equations

PDF of Lecture 3 slides

Solving Maxwell equations with FDTD scheme

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 essence, the electric field is a vector quantity (associated with lines) while the magnetic field is a bi-vector quantity (associated with surfaces). Hence, the most natural representation on a discrete grid utilizes this geometric fact to build a consistent scheme.

To couple the plasma to the field currents and charges need to be “deposited” on the grid in a careful manner. Current continuity needs to be 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.

The literature on geometric and Lagrangian and Hamiltonian methods is difficult for most plasma physicist to understand. The classic Dover textbook by Lovelock and Rund “Tensors, Differential Forms, and Variational Principles”.

Lecture 4: Hyperbolic equations, Finite-volume schemes

PDF of Lecture 4 slides

Hyperbolic Partial Differential Equations

Hyperbolic equations describe a broad class of physical problems and are essentially characterized by finite propagation speed of disturbances. Examples of hyperbolic equations include Maxwell equations, Euler equations for ideal fluids and ideal MHD equations.

To solve hyperbolic equations one needs to use special methods, in particular when shocks and other nonlinear phenomena are present (rarefaction waves, contact discontinuities, compression waves). These methods go by the name of “shock capturing schemes” and were originally developed in the aerospace community to solve for transonic and supersonic flows on airplanes and re-entry vehicles. They are widely used in astrophysics, but not so much in studying fusion plasmas.

Petr Lax, a mathematician from NYU’s Courant Institute won the 2005 Abel Prize (sort of “Noble Prize for Mathematics”) in part for his contribution to the theory and computations of hyperbolic PDEs.

A leading text-book on finite-volume methods to solve hyperbolic equations is by Randy LeVeque [LeVeque2002]. His earlier “Green Book” [LeVeque1992] is a very good introduction to the mathematics of hyperbolic equations and somewhat better than the 2002 “Red Book”.

We will focus on finite-volume and discontinuous Galerkin schemes for partial differential equations (PDEs), specifically fluid mechanics (Euler equations) and plasma physics (MHD equations, multi-fluid equations and the Vlasov-Maxwell system). We will look at schemes suited to shock dominated flows rather than problems on resistive time-scales. Vast majority of laboratory, space and astrophysical problems have complex interactions of shocks, turbulence and magnetic fields and the schemes in these lectures will help you solve equation to study such phenomena.

As reading for the lectures please read Chapters 1, 2, 3 and 5 in [LeVeque1992]. Some notes and references to material not covered in class:

  • van Dyke’s “Album of Fluid Motion” is an excellent source of beautiful pictures of fluid flow. See this link for a PDF of an older version of the book.
  • See eigensystem of Euler equations listed here and Maxwell equations here.
  • For ideal MHD equations the eigensystem is very complex. A listing is in [RyuJones1995] but it may be a good idea to rederive this and cross-check.

References

[Juno2019]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. http://doi.org/10.1016/j.jcp.2017.10.009
[CaryBrizard2009]Cary, J. R., & Brizard, A. J. (2009). “Hamiltonian theory of guiding-center motion”. Reviews of Modern Physics, 81 (2), 693–738. http://doi.org/10.1103/RevModPhys.81.693
[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. http://doi.org/10.1016/j.jcp.2013.02.009
[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. http://doi.org/10.1016/j.jcp.2013.08.021
[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. http://doi.org/10.1093/mnras/stu2164
[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. http://doi.org/10.1063/1.4818428
[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. http://doi.org/10.1063/1.4979989
[Esirkepov2001](1, 2) 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](1, 2) 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. http://doi.org/10.1016/j.jcp.2009.07.025
[Morrison2017](1, 2) Morrison, P. J. (2017). Structure and structure-preserving algorithms for plasma physics. Physics of Plasmas, 24 (5), 055502–21. http://doi.org/10.1063/1.4982054
[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. http://doi.org/10.1063/1.4913582
[LeVeque1992](1, 2) R.J. LeVeque. “Numerical methods for conservation laws”. Birkhäuser, 1992.
[LeVeque2002]R.J. LeVeque. “Finite volume methods for hyperbolic problems”. Cambridge University Press, 2002.