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.