In the previous article, we discussed how to integrate charged particle velocity in the presence of electric field. We now include the magnetic component. As you already know, magnetic field causes charged particles to rotate about the field line. This component is thus often called v×B (or v cross B) rotation.
Simple implementation for E=0, Forward Difference
First let’s consider a case without any electric field. The governing equation for velocity is then (here we are using bold typeface for vector quantities). A simple implementation of the integrator for 2D case with a magnetic field in the k direction appears to be:
/*grab magnetic field at current position*/ B=EvalB(x); /*get new velocity at n+1*/ v2 = v + q/m*B*v*dt; v2 = v - q/m*B*v*dt; /*update position*/ x2 = x + v2*dt; x2 = x + v2*dt; /*push down*/ v=v2; v=v2;
Here velocity is assumed to exist at half times in the spirit of the leapfrog scheme (i.e. v = v[n-1/2] and v2 = v[n+1/2]). If you implement this method (as done in UpdateVelocityForward in the attached code), you’ll see right away that it doesn’t work. As shown in Figure 1 below, the result is quite similar to what happened earlier in the case of forward difference integrator for the B=0 case. The particle keeps gaining energy, and instead of completing closed circles around the guiding center, it is continuously spiraling away.
Figure 1. Incorrect result obtained with the forward integration of vxB rotation. The particle is gaining numerical energy, as shown by its orbit spiraling away. The analytical result is a closed circular orbit at the Larmor radius, which is shown by the solid blue line.
Lorentz Integrator 1: Tajima’s Implicit Method
Getting the right solution requires taking approach similar to what was done previously for the electrostatic case. Instead of updating the velocity from time “n-1/2″ to “n+1/2″ using the velocity at “n-1/2″, we should use the average velocity at time “n”. This modifies our Lorentz force integrator to
This expression can be rewritten in matrix notation as
where I is the identity matrix, R is the unit rotation matrix given by
and ε is the scaling factor, . Unfortunately, this equation is implicit, and solving it requires performing a matrix inversion. The solution is given by
with the matrix inverse given by Tajima as
This method is implemented in the attached Java integrator as UpdateVelocityTajimaImplicit. As you can see in Figure 2 below, it indeed works. Not only is the energy conserved, but the computed Larmor radius is also right. We used 0.01T for the magnetic strength, the particle is an electron, and has initial tangential velocity of 100,000m/s. The Larmor radius, .
Figure 2. Result obtained with the implicit method from Tajima. This time, the particle trajectory is integrated correctly, and the energy is conserved.
Unfortunately, as you can see, this method is rather complicated, and involves a significant amount of calculation. At millions of particles per simulation, these calculations quickly add up into slow code performance…
Lorentz Integrator 2: Tajima’s Explicit Method
Tajima introduced a method that can be used if a small enough time step is selected such that ε<<1. In that case . By substituting and eliminating the quadratic term, we obtain
But as shown in Figure 3, this method is also incorrect for time steps practical to kinetic plasma simulations. This can be clearly seen by looking at the equation. In the absence of electric field, this integrator is identical to the original forward difference.
Lorentz Integrator 3: Boris Method
So are we stuck using the expensive implicit method? No, not quite. In 1970, Boris described an elegant alternative, which is now commonly known as the Boris Method. Boris method is the de facto standard for particle pushing in plasma simulation codes. Again, we are solving
Boris noticed that we can eliminate the electric field by defining
When these definitions are substituted into the original equation, we obtain pure rotation
Boris next utilized some basic geometry (see Figure 4-4a in Birdsall, p. 62) to derive the expression for performing the rotation. The first step is to find the vector bisecting the angle formed by the pre- and the (to be yet computed) post-rotation velocity. The angle through which the velocity will rotate in the given time step is, from geometry, . The vector form of this is . The bisector vector (v prime) is then
This “v prime” vector is perpendicular to both the magnetic field (the vector “t”) and the vector from “v minus” to “v plus”, the post-rotation velocity we are looking for. This connecting vector is again obtained from geometry as the cross product of “v prime” and a new vector “s”. This vector “s” is just a version of the rotation vector “t” scaled to satisfy the requirement that magnitude of velocity remains constant in the rotation. Mathematically speaking
To implement the Boris method, first obtain “v minus” by adding half acceleration to the initial velocity, per Equation 1. Then perform the full rotation according to Equations 3 and 4. Finally, add another half acceleration, as given by Equation 2.
All four methods described here are implemented in the attached Java program, ParticleIntegrator.java (Netbeans Workspace, zip). Java is actually a great language for scientific computing. Although in the past Java used to be much slower than C/C++, this is no longer case. The syntax is very similar to C++, but Java comes bundled with a large standard library which makes it very easy to do things such as add multithreading to your code. All this is possible in C/C++, but requires downloading additional libraries, compiling them for your architecture, and making sure they play well with your makefiles/workspace. None of this is necessary in Java, as it already comes preinstalled.
The program implements a particle push integrator than advances a single particle in 3D electromagnetic fields for a specified number of time steps. The fields at the particle velocity are evaluated by calling functions EvalE and EvalB. Currently these return constant values, but feel free to experiment to make the fields variable. A function called UpdateVelocity acts as a wrapper that calls the appropriate implementation. The Boris method is coded up in the attached Java program in the function UpdateVelocityBoris. To switch between the different methods, simply uncomment the appropriate call on line 72 in function UpdateVelocity. To switch from an electron to a heavier ion, simply modify the charge and mass in the Particle class definition on line 84.
- ParticleIntegrator.java (source code for plasma particle integrator)
- ParticleIntegrator.zip (NetBeans workspace, zipped file)
- Boris, J.P., The acceleration calculation from a scalar potential, Plasma Physics Laboratory, Princeton Univeristy, MATT-152, March 1970
- Boris, J.P., Relativistic plasma simulation-optimization of a hybrid code, Proceeding of Fourth Conference on Numerical Simulations of Plasmas, November 1970
- Birdsall, C.K., and Langdon, A.B., “Plasma Physics Via Computer Simulations“, Institute of Physics Publishing, Bristol and Philadelphia, 1991
- Tajima, T., “Computational Plasma Physics“, Westview Press, 2004