Magnetic fields in Physibounce

One of the things that has frustrated me about Physibounce is that when the magnetic field is on, the radius of the cyclotron motion of the particle isn't constant. Here I will write about why this problem occurred and what can be done to fix it.

In Physibounce, the magnetic fields are constant in time (other than turning on or off). When the magnetic field is turned on, the magnetic field can either point into or out of the screen. Let me call this the \hat{z} direction. The equations of motions for the particles in Physibounce are

 \frac{d\vec{v}}{dt} = \frac{\vec{F}}{m}+\omega_c \vec{v}\times\hat{z}

 \frac{d\vec{x}}{dt} = \vec{v}.

where \omega_c = \frac{eB}{m}, with B is the magnitude of the magnetic field, e is the charge of the particle, and m the mass. F is all forces acting on the particle that are not the magnetic field. Before now, Physibounce used an semi-implicit Euler method. Let us see why the semi-implicit Euler method is bad. For simplicity, take the case when there are no forces other than the magnetic field applied. The update equation for semi-implicit Euler is

 \vec{v}(t+\Delta t) = \vec{v}(t) +\omega_c\Delta t(v_y \hat{x}-v_x\hat{y})

 \vec{x}(t+\Delta t) = \vec{x}(t) + \vec{v}(t+\Delta t)\Delta t.

Assuming I measure distances from the center of the circular motion of the particle, ideally the radius of the orbit does not change. Since the magnetic force is causing the circular motion, it must be that m|\vec{v}(t)|^2/|\vec{x}(t)| = e |\vec{v}(t)| B, or equivalently, |\vec{v}(t)| = \omega_c |\vec{x}(t)|. Thus, enforcing that the radius of orbit does not change is equivalent to enforcing that kinetic energy is conserved as \frac{1}{2} m |\vec{v}(t)|^2 = \frac{1}{2} m \omega_c^2|\vec{x}(t)|^2 (remember that magnetic fields do no work). Consider the velocity squared at the updated timestep

 |\vec{v}(t+\Delta t)|^2 = |\vec{v}(t)|^2 + \omega_c^2\Delta t^2 |\vec{v}(t)|^2.

Thus, every timestep, a mistake of \frac{1}{2} m|\vec{v}(t)|^2\Delta t^2\omega_c^2 is made to the kinetic energy. Since this is proportional to the current kinetic energy, this means that the error gets worse and worse as time goes on, which is why semi-implicit Euler is bad for modeling magnetic forces. This is also why in Physibounce the radius of the orbit increases for a particle in a magnetic field.

One solution to the problem would be to use polar coordinates. Then the semi-implicit Euler update would be

 r(t+\Delta t) = r(t)

 \theta(t+\Delta t) = \theta(t) + \omega_c \Delta t.

which looks quite simple. In this form, it is obvious that the radius of rotation is conserved. However, if I wanted to use this method, I would have to write all of the other forces of Physibounce in polar coordinates. This would be quite a bit of work to work out, so I did not choose this method.

The solution to this problem while still using Cartesian coordinates is to use a leap-frog-like scheme, specifically one called the Boris pusher. The idea of the Boris scheme is to push the velocity by a half timestep using only the force \frac{F}{m}, then make a rotation due to the magnetic field, and then push the velocity another half timestep using the same force as before [1]. As before, the position is updated with the updated velocity. I calculate the following quantities

 \vec{v}^- = \vec{v}^{n} + \frac{\vec{F}^n}{m}\frac{\Delta t}{2}

 \vec{v}^+ = \vec{v}^{-} + \frac{\omega_c \Delta t}{1+\omega_c^2\Delta t^2/4}(v^{-}_y\hat{x}-v^-_x\hat{y})-\vec{v}^-\frac{\omega_c^2 \Delta t^2}{2(1+\omega_c^2\Delta t^2/4)}

Then I update

 \vec{v}^{n+1} = \vec{v}^+ + \frac{\vec{F}^n}{m}\frac{\Delta t}{2}.

I will again consider the case when the magnetic field is the only force, but the more general case is considered in [1]. In this case, the update has the nice property that

 |\vec{v}(t+\Delta t)|^2 = |\vec{v}(t)|^2.

Thus, kinetic energy is conserved, and hence the radius of orbit of the particle is constant. This creates a more realistic simulation and does not make the computation of all the other forces in Physibounce any more complicated. There is a minimal amount of code required to incorporate it into a simulation already using a semi-implicit Euler method. These pluses make it ideal for solving the problem in Physibounce.

The Boris pusher is widely used in plasma physics simulations as it is quite fast (in terms of number of flops required) and represents the dynamics without too much error. It is used in N-body and PIC simulations.

I have updated Physibounce (first update in about a year!) to reflect this.

1. Qin, H. et al., 2013. Why is Boris algorithm so good? Phys. of Plas. 20-8.

Random numbers from a distribution and the Lambert W function

When I was making Physibounce, one of the things that took me a while to figure out is something that is unfortunately barely noticeable in gameplay. When the speed of light becomes small, I wanted the speed distribution of particles to be a 2D Maxwell-Juttner distribution rather than the non-relativistic Maxwell speed distribution. The probability density function for particles allowed to be close to the speed of light is:

f(v) = \frac{\gamma(v)^4 m^2 c^2 v}{k_B^2 T^2(1+\frac{mc^2}{k_B T})} \exp\left(-\frac{mc^2}{k_B T}\left(\sqrt{1+\gamma(v)^2\frac{v^2}{c^2}}-1\right)\right).

Where v is the velocity of the particle, \gamma(v) = (1-\frac{v^2}{c^2})^{-1/2} is the Lorentz factor, k_B is the Boltzmann constant, m is the mass of the mass of the particle, T is the temperature of the system, and c is the speed of light. This distribution will look like a Maxwell distribution in the c\to\infty limit, and is normalized so that \int_0^c f(v) dv = 1.

The simplest way to generate random numbers following a distribution is the acceptance-rejection method. This is what I do in Physibounce when the speed of light is not changed. However, I was not able to find an expression for a suitable (and efficient) bound (or maximum value) for the distribution given above. Because of this, I decided to use inverse sampling instead. To do this, I first find the cumulative distribution function as

 F = \int_0^v f(v') dv'.

Amazingly, Mathematica was able to do this analytically. By the definition, it is clear F takes on values in the interval [0,1]. Thus, if I find a relation that expressed v in terms of F (equivalent to inverting the function), I have an expression that takes uniform values on the interval [0,1] and maps them to a velocity in the interval [0,c] that follows the desired distribution. Mathematica was also able to do the inversion, but the result was in terms of the Lambert W function.

The Lambert W function is the function that satisfies

 z = W(z) e^{W(z)}.

So it seems reasonable something like this would come up when inverting the cumulative distribution function. An important thing about this function is that it is double valued! There are two branches. When I asked Mathematica to solve the problem, it chose a branch and the branch it chose turned out to be the wrong one. I wondered for a while why I was not getting reasonable velocity values when I plugged in random numbers, and after a bit I realized I needed to choose the other branch.

The next problem I ran into was that I did know of a way to actually evaluate the Lambert W function in actionscript. So I turned to the asymptotic expansions. It turned out I was in a regime where the argument of the Lambert W tended to be large. I took the first few terms in the expansion about \infty for the Lambert W. I verified that over the parameter ranges I could expect in Physibounce that the difference between the Lambert W and the expansion was less than 1%. Since it's pretty hard to look at the speeds of a bunch of objects and determine the distribution of their speeds, I thought this was more than good enough, though I did include a check to make sure no particles were created with speed greater than c. In the end the code ended up only being a few lines, but quite a bit of thought went into each of the lines!

See the ActionScript code here.