The Boris Algorithm: simulating relativistic particles in external EM fields
Introduction
Problem statement
The motion of charged relativistic particles subjected to the interaction with an external electromagnetic field (EM) is governed by the Lorentz force (here expressed in SI units):
\begin{equation}\label{eq:lorentz} \frac{d(\gamma\mathbf{v})}{dt} = \alpha ( \mathbf{E} + \mathbf{v}\times \mathbf{B}) \end{equation}
where \(\gamma = (1 - v^2/c^2)^{-\frac{1}{2}}\) is the Lorentz factor and \(\alpha = q/m\) with \(q\) and \(m\) the particle’s charge and mass respectively. In the non-relativistic limit, the LHS of the equation becomes:
\begin{equation}\label{eq:classical_lorentz} \frac{d\mathbf{v}}{dt} = \alpha ( \mathbf{E} + \mathbf{v}\times \mathbf{B}) \end{equation}
which is a linear first-order differential equation.
The trajectory of point particles obeying Eq. \eqref{eq:lorentz} can be analytically found for some configurations of the EM fields. However, those scenarios are often an abstraction of much more complex field configurations found in nature and technology, giving rise to results that are only qualitatively correct.
Despite being physically interesting and useful to calculate approximate (sometimes their accuracy can be so close to reality that it coincides with it for every practical purpose) quantities, the analytical solutions are not available for arbitrary EM field configurations. In addition, the motion of charged particles creates itself electrical and magnetic fields which will modify the path (and hence the creation of other EM fields) of other particles. This makes the microscopic motion of charged particles chaotic and impossible to solve analytically for every practical situation. Furthermore, additional particle-particle interactions or external fields are needed to correctly model the behaviour of real particles (such as chemical and nuclear reactions and gravity).
If we are interested only in the collective motion of a (big enough) collection of particles, we can resort to the tools of statistical mechanics (leading to plasma physics); however, if we are concerned about the motion of the single particles, we do need to resort to numerical integration of Eq. \eqref{eq:lorentz}, which is the aim of the present work.
Physical relevance of the problem
The study of both individual and collective motion of charged particles arises in several key topics of modern physics. The study of plasma has applications in astrophysics: ranging from the physics of the solar corona to the motion of individual electrons or ions with extremely high gamma factor inside the galactic magnetic field, charged particles are everywhere in the cosmos and therefore the study of their motion is of top concern as a tool to probe the evolution of the cosmos as well as the working of stellar objects.
Relatively small packets of charged particles are used in particle accelerators to understand the fundamental physics of the universe: the precision needed to create and maintain highly relativistic charged particles in a coherent packet is a formidable undertaking that needs an accurate computation of several aspects of the motion of particles in complex magnetic fields.
A technological application of plasma physics with the potential to be a long-lasting change to how humanity lives is (controlled) nuclear fusion: both experimental and computational efforts are needed to understand the conditions needed to exploit the power of nuclear fusion in a practical and economical manner, since the manipulation of high-density, high-temperature plasma of different particles immersed in strong non-uniform magnetic fields is a very complex task.
Employed method
In this work, we are not concerned with the creation of EM fields from the particles, so that we can concentrate on the study of the motion of the particles inside external EM fields. We will also ignore all other particle-particle interactions as well as any other external force.
Numerically integrating equation \eqref{eq:classical_lorentz} is simpler due to its linear nature compared to the relativistic one. However, in real-case scenarios, numerical simulations often involve particles that are at least mildly-relativistic. Therefore, in this study, we examine the solution to equation \eqref{eq:lorentz} in both its classical and ultra-relativistic limits. To achieve this, we employ the widely-used Boris algorithm, which is widely used in astrophysical simulations.