Particle Pushing with Boris Method Matrix-Vector

This is going to a fairly technical post. I'm using the Boris method of charged particle pushing for my current project and just thought I'd share my notes. I first saw the Boris method described here : https://www.particleincell.com/2011/vxb-rotation/. Which provides a good explanation of why it was developed and compares it with other methods. 

However, the last time I did this I desired a simple form of computation based on a matrix/vector multiplication and addition. Something like

\[v^{n+1/2} = R v^{n-1/2} + A \]

Where R is the result of rotation in a magnetic field (not the R mentioned earlier in the linked article), and A is due to electric acceleration and interaction between electric and magnetic accelerations. And \(v^{n}\) is the velocity of the particle being stepped in time at the nth time marker. The position of the particle would be stepped in a leap-frog method as described in those documents.

I don't know if it has been put in this form before, it could have so I'm not trying to take credit for anything. I'm just following my notes. I start with the final equations listed in the link.

\[v^{n+1/2} = v^{n-1/2} + \frac{q \delta t}{m} \left (E + \frac{v^{n-1/2} + v^{n+1/2}}{2} \times  B\right ) \]

\[v^{n-1/2} = v^{-} - \frac{q \delta t}{2m} E \] 

\[v^{n+1/2} = v^{+} + \frac{q \delta t}{2m} E \]

\[v^{+} = v^{-} + \frac{q \delta t}{2m} \left (v^{+} + v^{-} \right ) \times B \]

They give a series of steps to velocity calculation. But I wanted a single multiply and add. I preferred renaming \(t = \phi \) for the next two steps. I think of \(\phi\) as like a angular vector of rotation due to the magnetic field.

\[\phi \equiv h B \]

\[h \equiv \frac{q \delta t}{2m} \]

\[v' = v^{-} + v^{-} \times \phi \]

\[v^{+} = v^{-} + v' \times \phi \frac{2}{1+ \phi^2}\]

Ok, just plugging step 3 into step 4.

\[ v^{+} = v^{-} + \left (v^{-} + v^{-} \times \phi \right ) \times \phi \frac{2}{1+ \phi^2} \]

Then start expanding the terms.

\[ v^{+} = v^{-}  + \frac{2}{1+ \phi^2} \left ( v^{-} \times \phi + (v^{-} \times \phi ) \times \phi \right )\]

\[ v^{+} = v^{-}  + \frac{2}{1+ \phi^2} \left ( v^{-} \times \phi + \phi \times (\phi \times v^{-}) \right )\]

Then use the vector identity \(a \times (b \times c) = (a \cdot c) b - (a \cdot b) c \), and \(\phi \cdot \phi = \phi^2 \)

\[ v^{+} = v^{-}  + \frac{2}{1+ \phi^2} \left ( v^{-} \times \phi + \phi (\phi \cdot v^{-} ) - \phi^2 v^{-} \right )\]

Now, pull out the \(v^{-}\). This will look really strange because of the 'dangling' vector operations. However, those can be converted into a matrix, which when a vector is multiplied through will have the desired result. But I will wait on computing that matrix.

\[ v^{+} = v^{-} \left (1  + \frac{2}{1+ \phi^2} \left ( (\times \phi)  + \phi (\phi \cdot  ) - \phi^2  \right ) \right ) \]

But this has to be put in terms of the actual current and next velocity. Basically, you can plug into the equations above to get it in terms of \(v^{n-1/2}\) and \(v^{n+1/2}\). Also re-arranged the terms a bit.

\[ v^{n+1/2} = \left ( v^{n-1/2} + \frac{q \delta t}{2m} E \right ) \left (\left (1 - \frac{2 \phi^2 }{1+ \phi^2} \right )  + \frac{2}{1+ \phi^2} \left ( (\times \phi )  + \phi (\phi \cdot  )  \right ) \right )  + \frac{q \delta t}{2m} E\]

Now expand the first term

\[ v^{n+1/2} = v^{n-1/2} \left (\left (1 - \frac{2 \phi^2 }{1+ \phi^2} \right )  + \frac{2}{1+ \phi^2} \left ( (\times \phi )  + \phi (\phi \cdot  )  \right ) \right )  + \frac{q \delta t}{2m} E \left (\left (1 - \frac{2 \phi^2 }{1+ \phi^2} \right )  + \frac{2}{1+ \phi^2} \left ( (\times \phi )  + \phi (\phi \cdot  )  \right ) \right ) + \frac{q \delta t}{2m} E\]

Ok, I won't torture you any further with this. I'm going to skip to the end and just tell you what the matrices are in terms of the fields. 

\[R = \left ( 1 - \frac{2 h^2 B^2}{1 + h^2 B^2} \right ) I + \frac{2}{1 + h^2B^2}( h (\times B) + h^2 (BB)) \]

\[ I = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\  0 & 0 & 1 \\ \end{bmatrix}\]

\[ (\times B) = \begin{bmatrix}0 & B_3 & -B_2 \\ -B_3 & 0 & B_1 \\  B_2 & -B_1 & 0 \\ \end{bmatrix}\]

\[ (BB) = \begin{bmatrix} B_1 B_1 & B_1 B_2 & B_1 B_3 \\  B_2 B_1 & B_2 B_2 & B_2 B_3 \\  B_3 B_1 & B_3 B_2 & B_3 B_3 \\ \end{bmatrix}\]

\[A = h \left ( 2 - \frac{2 h^2 B^2}{1 + h^2 B^2} \right ) E + \frac{2}{1 + h^2B^2}\left ( h^2 (E \times B) + h^3 (E \cdot B) B \right ) \]

The impressiveness of this method can maybe be seen now being 3rd order in the delta time (the h constant), even taking into account for \(E\times B\) drift and energy being gained or lost to the fields through \(E\cdot B\). While at the same time conserving energy in the pure rotation part.

The idea is R and A would be pre-calculated once per step (assuming it is time-dependent) as a function of position given the fields E and B, with a total of 12 parameters. Then, when particles are pushed the matrix at position \(x^{n}\) is looked up and a multiply and add is done to update the velocity, which is a total of 9 multiply and 9 add operations per particle.