Particle Push in Magnetic Field (Boris Method)
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 \(d\mathbf{v}/dt = q/m(\mathbf{v}\times \mathbf{B})\) (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[0] = v[0] + q/m*B*v[1]*dt; v2[1] = v[1] - q/m*B*v[0]*dt; /*update position*/ x2[0] = x[0] + v2[0]*dt; x2[1] = x[1] + v2[1]*dt; /*push down*/ v[0]=v2[0]; v[1]=v2[1];
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
$$
\begin{array}{rcl}
\displaystyle\mathbf{v}^{n+1/2} &=& \mathbf{v}^{n-1/2}+\dfrac{q}{m}\left[\mathbf{E} + \mathbf{v}^{n}\times\mathbf{B}\right]\Delta t \\
\displaystyle&=& \mathbf{v}^{n-1/2}+\dfrac{q}{m}\left[\mathbf{E} + \dfrac{\mathbf{v}^{n-1/2}+\mathbf{v}^{n+1/2}}{2}\times\mathbf{B}\right]\Delta t
\end{array}
$$
This expression can be rewritten in matrix notation as
$$
\displaystyle\left[\mathbf{I}-\mathbf{R}\epsilon\right] \mathbf{v}^{n+1/2} =
\left[\mathbf{I}+\mathbf{R}\epsilon\right] \mathbf{v}^{n-1/2}+\frac{q}{m}\mathbf{E}\Delta t
$$
where I is the identity matrix, R is the unit rotation matrix given by
$$
\mathbf{R}=\dfrac{1}{B}\left(\begin{array}{ccc}
0 & B_z & -B_y\\
-B_z & 0 & B_x\\
B_y & -B_x & 0\end{array}\right)
$$
and ε is the scaling factor, \(\epsilon = \dfrac{qB}{m}\Delta t/2=\omega_c\Delta t/2\). Unfortunately, this equation is implicit, and solving it requires performing a matrix inversion. The solution is given by
$$
\displaystyle\mathbf{v}^{n+1/2}=[\mathbf{I}-\mathbf{R}\epsilon]^{-1}[\mathbf{I}+\mathbf{R}\epsilon]\mathbf{v}^{n-1/2}
+[\mathbf{I}-\mathbf{R}\epsilon]^{-1}\mathbf{E}\dfrac{q}{m}\Delta t
$$
with the matrix inverse given by Tajima as
$$
[\mathbf{I}-\mathbf{R}\epsilon]^{-1} = \dfrac{1}{\det(\mathbf{I}-\mathbf{R}\epsilon)} \left(
\begin{array}{ccc}
1+\epsilon^2b_x^2 & \epsilon b_z + \epsilon^2b_x b_y & -\epsilon b_y+\epsilon^2b_x b_z\\
-\epsilon b_z + \epsilon^2 b_xb_y & 1+\epsilon^2b_y^2 & \epsilon b_x + \epsilon^2b_yb_z\\
\epsilon b_y + \epsilon^2b_xb_z & -\epsilon b_x+\epsilon^2 b_y b_z & 1+\epsilon^2b_z^2
\end{array}\right)
$$
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, \(r_L=mv_\perp/(|q|B)=5.68\times10^{-5}\text{m}\).
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 \([\mathbf{I}-\mathbf{R}\epsilon]^-1 \simeq \mathbf{I} + \mathbf{R}\epsilon + \ldots\). By substituting and eliminating the quadratic term, we obtain $$ \begin{array}{rcl} \mathbf{v}^{n+1/2} &=&[\mathbf{I}+\mathbf{R}\epsilon][\mathbf{I}+\mathbf{R}\epsilon]\mathbf{v}+ [\mathbf{I}+\mathbf{R}\epsilon]\mathbf{E}\dfrac{q}{m}\Delta t\\ &=&\dfrac{q}{m}\mathbf{E}\dfrac{\Delta t}{2}+(\mathbf{I}+\mathbf{R}\epsilon)\left[\mathbf{v}^{n-1/2}+\dfrac{q}{m}\dfrac{\Delta t}{2}\mathbf{E}\right] \end{array} $$ 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.
Figure 3. Result obtained with the explicit method from Tajima. The result is again incorrect.
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
$$
\displaystyle \frac{\mathbf{v}^{n+1/2}-\mathbf{v}^{n-1/2}}{\Delta t} = \frac{q}{m}\left[\mathbf{E}+\frac{\mathbf{v}^{n+1/2}+\mathbf{v}^{n-1/2}}{2}\times\mathbf{B}\right]
$$
Boris noticed that we can eliminate the electric field by defining
$$
\begin{array}{rcl}
\mathbf{v}^{n-1/2} &=&\mathbf{v}^{-} – \dfrac{q\mathbf{E}}{m}\dfrac{\Delta t}{2} \qquad \qquad (1)\\
& &\\
\mathbf{v}^{n+1/2} &=& \mathbf{v}^{+} + \dfrac{q\mathbf{E}}{m}\dfrac{\Delta t}{2} \qquad\qquad (2)
\end{array}
$$
When these definitions are substituted into the original equation, we obtain pure rotation
$$
\dfrac{\mathbf{v}^{+}-\mathbf{v}^{-}}{\Delta t} = \dfrac{q}{2m}\left(\mathbf{v}^{+}+\mathbf{v}^{-}\right)\times\mathbf{B}
$$
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, \(\tan(\theta/2)=-(qB){m}\Delta t/2\). The vector form of this is \(\mathbf{t}\equiv-\hat{\mathbf{b}}\tan\theta/2 = (q\mathbf{B}/m)\Delta t/2\). The bisector vector (v prime) is then
$$
\mathbf{v}’=\mathbf{v}^{-}+\mathbf{v}^{-}\times\mathbf{t} \qquad\qquad(3)
$$
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
$$
\mathbf{v}^{+}=\mathbf{v}^{-}+\mathbf{v}’\times\mathbf{s}\qquad\qquad(4)
$$
where
$$
\mathbf{s}=\dfrac{2\mathbf{t}}{1+t^2}
$$
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.
Figure 4. Result obtained using the Boris method. Energy is conserved.
Implementation
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.
Attachments:
- ParticleIntegrator.java (source code for plasma particle integrator)
- ParticleIntegrator.zip (NetBeans workspace, zipped file)
References:
- 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
Feel free to leave a comment if you have any questions or something isn’t clear. Thank you for stopping by!
Nice article, I didn’t know about Tajima although my Boris implementation was in fact Tajima. It wasn’t too slow since the difference in needed floating point operations isn’t that serious for a 2D case.
So implicit Tajima and Boris produce the same trajectory (assuming infinite precision). Rearranging and simplyfing the Tajima expression for “v plus” gives Boris, so Boris could be derived without geometric interpretation as a rotation, but with just reducing the number of floating point operations in mind.
Correction: What I said about the relation of Tajima and Boris is for the case where no electric field is present. I am not sure for the general case.
That definitely makes sense. Thanks for pointing this out! I sort of figured that in the end the implicit and the multi-step Boris method should be the same but I didn’t really feel the motivation to prove it. I also implemented the implicit solver the way it’s written in Tajima’s book, with actual matrix inversion and matrix multiplication. I bet that bunch of terms end up cancelling out in the end. But the way it is written, this implicit method took about 2 seconds for some 10 million time steps, while the Boris method took about about 1.2. So roughly 60% slower or so (with the non-optimized versions)…
I converted this example to javascript to run in a webpage. You can see the results here.
http://dl.dropbox.com/u/5095342/PIC/pic.html
This site uses html5 so you will need a current version of your browser to see algorithm runtimes. IE9, firefox 4 or current google chrome web browser.
Thanks John! It would be great if you could add support for graphing the results. But I think this would require running the code as an applet, or alternatively server-side.
I updated the page so that it now displays a graph of the output. It uses the html5 “canvas” element to display the graph.
Nice! Thank you for doing this, John. I had no idea it was so easy to draw using Javascript (I am definitely not up to date when it comes to HTML5). I’ll try to do something similar for future posts to give visitors something to play with.
Hey there, nice article!
Quick question about your reference to the Boris method – you have it as MATT-152 whereas a quick web search found mostly references to MATT-769 – I cant seem to get any text from either report (maybe im looking in the wrong place) could you confirm which report it is?
Hi, thanks for all these nice tutorials. I have never used C++ or Java, but I am starting now, so the tutorials here are useful. One question is about vectorization – is that also a feature that can speed up things a bit. E.g. in other high-level languages I am using, to avoid loops its enough to write:
t = part.q/part.m*B*0.5*dt
instead of placing the calculation of each vector element in a loop, as you have in the available code. I use this option as much as possible, also to make calculations for many particles simultaneously. Is that possible in Java or C?
HI Elias, not that I know of, unless you start dealing with extensions such as CUDA. The reason why vector operations are so much faster in Matlab is that Matlab is an interpretive language that doesn’t precompile the code. As Matlab chugs through your code, it first translates each instruction from its syntax to the machine language before executing it. This happens even for code inside loops. So if you have a code such as
the instruction inside the loop (sum = sum + data[i]) ends up getting compiled 1000 times. This is a huge waste of computational resources. A vector statement results in only a single compilation. Java and C/C++ are both precompiled languages – C/C++ codes get translated to the machine language directly, while Java is turned into a byte code – and as such you don’t have to worry about issues like this. Unless you are actually using a vector computer, a single line syntax will not gain you much except for perhaps a slicker looking code. The way to get an actual speed up is by performing the calculation on a the graphics card (GPU) using language extensions such as CUDA. Graphics cards are in reality highly optimized vector computers that can perform operations on multiple components of an array concurrently. But that’s a topic for a future discussion…
Hi there,
I am also running similar particle codes. Do you know of how much integration time is needed for such trajectories or it is just upto the user? In other words, does it need to be some gyro periods long? Also, with regards to non-guiding center motion, such integration times could be any number…..how does one determine such numbers for these nonlinear problems (spatial dependence in B)?
Thanks
Hi Rushat, are you talking about the size of the simulation time step? If so, in a fully-kinetic simulation, where you simulate electrons and ions as particles, you will in general have two frequencies to consider: the plasma frequency and the cyclotron frequency. You need to set you simulation such that you resolve the higher of these. Also, in order for this Boris integration scheme to work, you need to take multiple steps about the orbit. I did some parametric studies as part of my dissertation work, and found (if I remember right) that the error started being somewhat acceptable with 15 steps per orbit, and reduced to below 0.1% with 75 steps per orbit. This error took into account both the displacement of the guiding center and the correct Larmor radius. In my simulation I ended up setting the time step based on the highest value of the magnetic field in the simulation, since the cyclotron frequency scales linearly with the field strength.
Hey Lubos,
Thanks for the reply. I guess I did not explain the problem correctly. The thing is I am not doing any PIC type simulation. I am just running for trajectories of particles in a given magnetic field. I am not solving it self-consistently. So, I am wondering if for my problem, what would be a sufficient ending time (normalized to gyro-frequency) for eg for protons……200*gyro-freq or 200*gyro-freq.
Thanks
I meant 200*gyro-freq or 2000*gyro-freq as the ending time??
I guess I am still confused 🙂
There is no “defined” ending time. This depends on your simulation, on what exactly you are simulating. If for instance you want to simulate one second of real time, you will need enough time steps to capture this (n_it = 1 sec / delta_t). On the other hand, PIC simulations are typically run until steady state which is defined by some parameter no longer changing. For a plume simulation, where you continuously inject particles from the source, the steady state is reached once the total particle count no longer changes. But if you preload the sim with a fixed number of particles and simply let these evolve, you will want to run the simulation until something else (velocity, temperature, mean position?) stops changing significantly.
Hi,
Nice tutorial, thanks. Do you know of a similarly well-written example of a relativistic Boris pusher? Is is as simple as using the gamma factor at the ‘old’ time step, or does something more clever have to be done?
Thanks
Hi Jason, no and sorry, I am not familiar with relativistic PIC codes. All the work I have ever done was with velocities < 1e7 m/s.
look in the Birdsall reference.
in short, yes you are correct, it is just a matter of applying gamma but it must be done in a properly centered fashion.
Charlson, do you have an example you could share?
For relativistic PIC code, where do I need to apply gamma?
Hi. I am trying to simulate particle movement in presence of a magnetic filed, But, I don’t know how to find current values. Of course, I have used the relation current=rho*velocity, But , I think this is not true. Can you help me to obviate my problem? (I am sorry for my poor English)
Hi, very helpful article you write. I download your java program and run it in NetBeans. I modified some pramaters and it worked. However I just got the larmor radius as results. How can I got the full results like positions in every time steps?
There should be a file called “trace.txt” that contains these in the directory from where you ran the program.
Here is my attempt to put it back as a matrix-vector product like in part 1, but derived by Boris method.
http://braingab.blogspot.com/2015/10/particle-pushing-with-boris-method.html
Hello,
in the Boris method, in the denominator of the s vector, the s^2 is the (sx^2+sy^2+sz^2) or the si^2 where i=x,y,z for each vector component?
The denominator is (1+t^2) if I am looking at the same equation as you. That “t” is the magnitude of the “t” vector so t^2 is tx^2+ty^2+tz^2.
Hi Lubos, I have a question. If we consider the vxB part only, is the Boris mover an alternative implementation of the Tajima’s implicit mover? They seem to solve the same original equation to push v_n-1/2 to v_n+1/2. Also, would they be stable under crazily bad dx or dt, like dt * Omega_ce ~ 1e5?
how is “s” become “2t/1+t^2”?
i have question related to boris code. code which i am writing has 3 different normalized values of velocities.now i have to convert it in ms-1.i am unable to find normalization factor with which i multiply it with to convert it in ms-1. kindly help me
in my case i am using the velocity v=(0.1,0.01,0)^T. kindky help me to convert in ms-1.
There is a typo in the equation s=2t/(1+t^2) for Boris’s algorithm. The t in the deonominator should be bolded since it is a vector. From the way it is typed, it looks like one is vec(t) and the other is t like time as in $\delta t$ you have defined earlier.
This has been a very helpful first introduction into numerical methods.
Thank you!