Particles Simulation in EM Fields
In many fields of physics, it is important to study the dynamics of single particles, but there are few cases in which an analytical solution is possible.
For example, if I consider a charged particle in motion interacting with the electromagnetic field, the equation that describes its dynamics is the following:
\begin{equation} \label{eq:1} \frac{d\mathbf{v}}{dt}=\frac{q}{m}(\mathbf{E}+\mathbf{v}\times\mathbf{B}) \end{equation}
The analytical solution is possible in a few cases, in particular with simple $ \mathbf{E} $ and $ \mathbf{B} $ configurations. It is therefore essential to resort to numerical methods, in particular in this case numerically integrating the equation of motion.
The first step towards the numerical solution is the rewriting of SI units in a dimensionless form, to avoid any numerical cancellations due to the different scales of the quantities present.
\begin{equation}\mathbf{\mathbf{\tilde{v}}}=\frac{\mathbf{v}}{v_0} \;\;\;\;\tilde{t}=\frac{t}{t_0} \;\; \; \tilde{\mathbf{E}}=\frac{\mathbf{E}}{E_0} \; \;\;\tilde{\mathbf{B}}=\frac{\mathbf{B}}{B_0}\end{equation}
dove $v_0,\;t_0,\; E_0,\; B_0$ they are constants with the same dimensions as the initial variable to obtain pure dimensionless numbers.
I can rewrite $\mathbf{v},\;t,\;\mathbf{E},\;\mathbf{B}$ as
\begin{equation} \mathbf{v}=\tilde{\mathbf{v}}\;v_0 \; \; \; t=\tilde{t}\;t_0 \;\; \; \mathbf{E}=\tilde{\mathbf{E}}\;E_0 \;\; \;\mathbf{B}=\tilde{\mathbf{B}}\;B_0 \end{equation} I substitute in the equation \eqref{eq:1}:
\begin{equation} \frac{v_0d\mathbf{\tilde{v}}}{t_0d\tilde{t}}=\frac{q}{m}(\mathbf{\tilde{E}}E_0+\mathbf{\mathbf{\tilde{v}}}v_0\times\mathbf{\tilde{B}}B_0) \quad \quad \end{equation}
\begin{equation} \frac{d\mathbf{\tilde{v}}}{d\tilde{t}} =\frac{qt_0}{mv_0}(\mathbf{\tilde{E}}E_0+\mathbf{\tilde{v}}v_0\times\mathbf{\tilde{B}}B_0) \end{equation}
\begin{equation}\frac{d\mathbf{\tilde{v}}}{d\tilde{t}}=\frac{qt_0 E_0}{mv_0}\mathbf{E}+\frac{qt_0}{mv_0}v_0 B_0 v \times \mathbf{B} \end{equation}
\begin{equation}\frac{qt_0 E_0}{mv_0}=1 \;\;\;\;\; t_0=\frac{mv_0}{q E_0}\end{equation}
I substitute $t_0$ in the coefficient before $\mathbf{v} \times \mathbf{B}:$
\begin{equation}\frac{qB_0}{m } \frac{m v_0}{q E_0}=\frac{B_0}{E_0}{v_0}=1 \;\;\;\;\; E_0=B_0v_0\end{equation}
In this way I get the dimensionless form of the Lorentz equation (omitting the tildes):
\begin{equation}\frac{d\mathbf{v}}{dt}=\mathbf{E}+\mathbf{v} \times \mathbf{B} \end{equation}
It doesn’t look much different than the first version but now the quantities are dimensionless.
With this non-dimensionalization procedure, I obtain three reference scales of the quantities I deal with using numerical methods.
\begin{equation}t_0 = \frac{m v_0} {q E_0}\; [s] \qquad v_0 = \frac{E_0}{B_0} \; \left[\frac{m}{s}\right] \qquad l_0 = v_0\;t_0 \; [m]\end{equation}
Choice of the numerical integrator
Regarding the choice of the integration method of the differential equation, I compared three different integrators: Runge-Kutta 2nd order, Runge-Kutta 4th order and as the system is Hamiltonian a symplectic integrator that preserves the total system energy.
The symplectic integrators Position-Verlet and Velocity-Verlet often used to integrate the equations of motion are not suitable as one of the hypotheses of these methods is that the force is independent of the velocity.
I, therefore, used the Boris Pusher a symplectic integrator suitable for integrating the equation of motion of a particle subject to the Lorentz force. \
To demonstrate the advantages of the Boris integrator I will analyze a simple case that is the motion of a charged particle immersed in a constant magnetic field.
One particle and magnetic field B
I consider two cases one in which the velocity of the particle is parallel to B and in which the velocity is perpendicular to B.
Parallel Velocity
\begin{equation}\mathbf{v}_0=(0,0,v_z) \qquad v_z=1\end{equation} \begin{equation}\mathbf{B}=(0,0,B_z) \qquad B_z=1\end{equation}
\[\begin{cases} \frac{dv_x}{dt}= 0 \\ \frac{dv_y}{dt}= 0 \\ \frac{dv_z}{dt}= 0 \end{cases}\]So the $ \mathbf {v} $ will be constant, equal to the initial $ v_0 $ and the particle will move uniformly in the direction $ z $.
Perpendicular velocity
\begin{equation}\mathbf{v}_0=(v_x,0,0) \qquad v_x=1 \end{equation}
\begin{equation}\mathbf{B}=(0,0,B_z) \qquad B_z=1\end{equation}
\begin{equation}
\begin{cases}
\frac{dv_x}{dt}= 0 \
\frac{dv_y}{dt}= - Bv_x\
\frac{dv_z}{dt}=0
\end{cases}
\end{equation}
I equate the Lorentz force with the centripetal force:\begin{equation}F_L=F_c \quad vB=\frac{v^2}{R} \quad B=\frac{v}{R} \quad R=\frac{v}{B}=1\end{equation}
Now I calculate the rotation period by substituting R:\begin{equation}T=\frac{2\pi}{\omega}=\frac{2\pi R}{v} \qquad T=\frac{2\pi}{v}\frac{v}{B} = \frac{2\pi}{B} = 2\pi\end{equation}
The numerical results coincide with the expected analytical results if the units of measurement are reported correctly.
Choosing as arbitrary values $E = 0.01 \; \frac{V}{m}$ e $B = 10^{-8} \; T $ and therefore considering an electron as a particle $q = 1.6 \times 10^{-19} \;C $ e $m = 9 \times 10^{-31} \; kg$ we obtain the scales with respect to which we have dimensioned the problem:
\begin{equation}t_0 = \frac{m v_0} {q E_0} = 5.625 \times 10^{-4} \; s \qquad l_0 = \frac{E_0t_0}{B_0} = 562.5 \; m \qquad v_0 = \frac{l_0}{t_0} = 10^6 \; \frac{m}{s}\end{equation}
For example we wanted to convert from dimensionless quantity to SI unit:
- the radius of the circular orbit: $R = 1 \times l_0 = 562.5 \; m $
- the period of the circular orbit: $T = 2\pi \times t_0 = 3.53 \; s $
- the kinetic energy of the circular orbit: $ E = \frac{1}{2}m \times v_0^2 = 4.5 \times 10^{-19} \;J$
Comparison between Rk2, Rk4, and Boris Pusher
I evaluate the different behavior of the three integrators considering how the energy of the particle varies during integration. The physical system considered conserves the initial energy so I expect that the numerical simulation will also preserve it.\ \begin{equation}\Delta E=\frac{|E-E_0|}{E_0} = \frac{|\mathbf{v}^2 - \mathbf{v}_0^2|}{\mathbf{v}_0^2}\end{equation} I evaluate the relative error on the energy considering only the kinetic energy of the particle as the magnetic field B does not supply energy as it is always perpendicular to the direction of displacement, therefore it affects only the direction of the velocity.
The error for RK2 and RK4 reaches relevant values and would continue to grow, while Boris’ algorithm maintains a constant negligible error between a $ 10 ^ {- 16} $ and $ 10 ^ {- 14} $. We verify that the use of a symplectic supplement allows the conservation of the initial energy.
ExB Drift
Analytical Solution
Now let’s consider the motion of a moving particle in the presence of an electric field $ \mathbf {E} $ in addition to the magnetic field $\mathbf{B}$.
\begin{equation}\mathbf{B}=(0,0,B_z) \end{equation}
\begin{equation}\mathbf{E}=(0,E_y,0) \end{equation}
I can write the particle $ \mathbf {v} $ as the sum of a parallel component and a perpendicular to $\mathbf{B}$:
$\mathbf{v} = \mathbf{v_{\parallel}} + \mathbf{v_{\perp}}$.
In turn $\mathbf{v_{\parallel}} = \mathbf{v_d} + \mathbf{v_g}(t)$
assuming that the drift speed is constant over time.
The idea is to lead us back to a simpler motion by making a change of reference system.
In fact, I can rewrite \eqref {eq: 1} as:
\begin{equation}
\frac{d\mathbf{v_g}}{dt} = \mathbf{E} + (\mathbf{v_d} + \mathbf{v_g}) \times \mathbf{B}\end{equation}
\begin{equation}
\frac{d\mathbf{v_g}}{dt} = \left(\mathbf{E} + \mathbf{v_d} \times \mathbf{B}\right) + \mathbf{v_g} \times \mathbf{B}
\end{equation}
I can cancel the term with the electric field by imposing a certain $\mathbf{v_d}$:
\begin{equation}
\mathbf{E} + \mathbf{v_d} \times \mathbf{B} = 0
\end{equation}
I apply the external product to both members and using a vector identity, I get:
\begin{equation}
\mathbf{v_d} = \frac{\mathbf{E} \times \mathbf{B}}{|\mathbf{B}|^2}
\end{equation}
I moved to a reference system whose $ \mathbf {E} $ field is null, so I find a simple gyration motion:
\begin{equation}
\frac{d\mathbf{v_g}}{dt} = \mathbf{v_g} \times \mathbf{B}
\end{equation}
If I consider $\mathbf{B}=(0,0,B_z)$ e $\mathbf{E}=(0,E_y,0)$ :
\begin{equation}
\mathbf{v_d} = \frac{\mathbf{E} \times \mathbf{B}}{|\mathbf{B}|^2} = \frac{E_y}{B_z} \vec{i}
\end{equation}
I can rewrite the equations of motion as:
\begin{equation}
\begin{cases}
\frac{d(v_x - v_d)}{dt}=v_y B_z \
\frac{dv_y}{dt}=-(v_x - v_d)B_z + E_y \
\frac{dv_z}{dt}=0
\end{cases}
\end{equation}
The solution will be:
\begin{equation}
\begin{cases}
v_x = \frac{E_y}{B_z} (1 + \cos(t))
v_y = - \frac{E_y}{B_z}\sin(t)
v_z = 0
\end{cases}
\end{equation}
The trajectory of the particle is a cycloid along the x-direction. The motion is unexpected having an electric field acting along $ y $ we would expect the particle to accelerate in this direction.
In reality, when the speed of the particle increases, the force resulting from the interaction with $ \mathbf {B} $ also increases, always normal to the trajectory.
The charge stopped at the origin undergoes an acceleration along y due to $ \mathbf {E} $. As soon as the charge begins to move the magnetic force acts in the direction normal to the trajectory.
This is until the end of the period when the resulting force is zero, so the charge stops, and the motion is repeated periodically.
Integrator Convergence
In the previous case of gyration, we studied the energy conservation properties of the three integrators, now we would like to study their convergence properties. To do this I will compare the analytical solution found in the previous section with the numerical solution. This will be repeated for different timesteps $\Delta t$.
As values for $\mathbf{E}$, $\mathbf{B}$ e $\mathbf{v_0}$ of the particles we will use: \begin{equation}\mathbf{B}=(0,0,1) \end{equation} \begin{equation}\mathbf{E}=(0,\frac{1}{2},0) \end{equation} \begin{equation}\mathbf{v_0} = (1,0,0)\end{equation}
Let’s verify that RK2 and Boris Pusher are integrators of the $ 2 ^ or $ order: a factor of 10 for $ \Delta t $ results in a reduction of 100 in the error. RK4 is of the $ 4^o $ order so a factor 10 of the $ \Delta t $ corresponds to a reduction of $ 10 ^4 $ for the error.
Numerical Solution
As a supplement we use Boris Pusher, as despite a lower precision than RK4 we get a sufficient accuracy: $10^{-4}$ con timestep $\Delta t = 10^{-2}$.
Come valori per $\mathbf{E}$, $\mathbf{B}$ e $\mathbf{v_0}$ della particella useremo:
\begin{equation}\mathbf{B}=(0,0,1) \end{equation}
\begin{equation}\mathbf{E}=(0,\frac{1}{2},0) \end{equation}
\begin{equation}\mathbf{v_0} = (1,0,0)\end{equation}
It can be seen that the particle does not have a net gain in kinetic energy. We can explain this thanks to the derivation of the previous section: if we move in a reference system that moves with $ \mathbf {v_d} $ we cancel the field $ \mathbf {E} $, so the particle interacts only with $ \mathbf {B} $. $ \mathbf {B} $ does not do work on the particle, in fact, it does not modify the velocity modulus but only the direction. To be more precise, this change of reference system can only be done if the relativistic invariant $ E ^ 2 - B ^ 2 $ is conserved, so $ | \mathbf {E} \ | <| \mathbf {B} | $. If $ | \mathbf {E} | > | \mathbf {B} | $ the invariant is no longer conserved and I cannot make this change of reference system.
Particles near an X point
I place $ N = 10 ^ 4 $ particles in a square $ L \times L $ with uniformly distributed position and velocity. In particular, the velocity of each will be $ | v | = 0.1 $ and a random direction obtained with a random angle between $ 0 $ and $ 2 \pi $.
The particles are immersed in an electric and magnetic field:
\begin{equation}\vec{E}=\left(0,0,\frac{1}{2} \right)\end{equation}
\begin{equation}\vec{B}=\left(\frac{y}{L},\frac{x}{L},0\right) \end{equation}
I write the equations of motion:
\begin{equation}
\begin{cases}
\frac{dv_x}{dt}= 0 \
\frac{dv_y}{dt}= 0 \
\frac{dv_z}{dt}= \frac{1}{2} + \frac{x}{L}v_x - \frac{y}{L}v_y
\end{cases}
\end{equation}
The particles perform a peculiar motion, in fact the particles located in the central area of the square are accelerated more towards $ z $ giving effect to a sort of jet of particles.
Let’s check for which particles in our domain it is possible to change the reference system:
\begin{equation} |\mathbf{E}| < |\mathbf{B}|\end{equation}
\begin{equation}\frac{1}{2} < \sqrt{\frac{x^2}{L^2} + \frac{y^2}{L^2}} \end{equation}
\begin{equation} x^2 + y^2 > \frac{L^2}{4}\end{equation}
so for particles outside the circumference with radius $ R = \frac {L} {2} $ it is possible to cancel the field $ \mathbf {E} $ and they will not be accelerated. While the internal particles will accelerate towards $ z $.
Conclusions
We have considered three different examples of particles in electric and magnetic fields, in particular for the last case the numerical simulation has allowed us to understand the non-trivial dynamics of the system.
Code
If you’re interested take a look at the code implemented on my Github.
References
I based my work on two different courses taken at in my master:
- Numerical algorithms for physics
- Plasma and fluid matter physics
both held by prof. A. Mignone
