Charged Particle Trajectories in Electric and Magnetic Fields
In the wake of the first week of the course about plasma physics, we will use some simple integrations to plot different trajectories of charged particles in magnetic and electric fields.
Charged particles in constant magnetic fields¶
The equation of motion for a charged particle in a magnetic field is as follows:
$$ \frac{d \vec{v}}{ dt } = \frac{q}{m} (\vec{v} \times \vec{B}) $$
We choose to put the particle in a field that is written
$$ \vec{B} = B \vec{e_x} $$
We thus expect the particle to rotate in the $(y, z)$ plane while moving along the $x$ axis.
Let's see how we can implement this using the integrators found in the scientific python stack.
# imports
import numpy as np
from scipy.integrate import ode
To use ode, we need to define the derivatives of velocity, which is computed using the Lorentz force and of position, which is just velocity:
def newton(t, Y, q, m, B):
"""Computes the derivative of the state vector y according to the equation of motion:
Y is the state vector (x, y, z, u, v, w) === (position, velocity).
returns dY/dt.
"""
x, y, z = Y[0], Y[1], Y[2]
u, v, w = Y[3], Y[4], Y[5]
alpha = q / m * B
return np.array([u, v, w, 0, alpha * w, -alpha * v])
Now that we have defined these functions, we can intialize the integrators and their methods.
r = ode(newton).set_integrator('dopri5')
Finally, we need to define initial conditions:
t0 = 0
x0 = np.array([0, 0, 0])
v0 = np.array([1, 1, 0])
initial_conditions = np.concatenate((x0, v0))
Let's now set the conditions on our integrators and solve the problem using time stepping. We assume the charged particle has unit mass and unit charge.
r.set_initial_value(initial_conditions, t0).set_f_params(1.0, 1.0, 1.0)
positions = []
t1 = 50
dt = 0.05
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
positions.append(r.y[:3]) # keeping only position, not velocity
positions = np.array(positions)
Let's check the integration results. We expect a circle in the $(y, z)$ plane:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(positions[:, 1], positions[:, 2])
[<matplotlib.lines.Line2D at 0x790cba8>]
This looks as expected! Now, let's do some 3D:
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot3D(positions[:, 0], positions[:, 1], positions[:, 2])
B1 = np.array([x0[0], x0[1], -1])
B2 = np.array([60, 0, 0])
B_axis = np.vstack((B1, B1 + B2))
ax.plot3D(B_axis[:, 0],
B_axis[:, 1],
B_axis[:, 2])
plt.xlabel('x')
plt.ylabel('y')
ax.set_zlabel('z')
ax.text3D((B1 + B2)[0], (B1 + B2)[1], (B1 + B2)[2], "B field")
<mpl_toolkits.mplot3d.art3d.Text3D at 0x7f93dd8>
Here we get the classical result: in a magnetic field, particles loop around the magnetic field axis.
We can do an animation of this type of movement using the code below:
from tempfile import NamedTemporaryFile
import base64
from matplotlib import animation
VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
f = NamedTemporaryFile(suffix='.mp4', delete=False)
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
f.flush()
video = open(f.name, "rb").read()
f.close()
anim._encoded_video = base64.b64encode(video).decode('utf-8')
return VIDEO_TAG.format(anim._encoded_video)
from IPython.display import HTML
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_html(anim))
FRAMES = 50
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
def init():
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# animation function. This is called sequentially
def animate(i):
current_index = int(positions.shape[0] / FRAMES * i)
ax.cla()
ax.plot3D(positions[:current_index, 0],
positions[:current_index, 1],
positions[:current_index, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=FRAMES, interval=100)
# call our new function to display the animation
display_animation(anim)
What happens when the charge is negative?¶
Another interesting animation to do is when the charge is negative.
r.set_initial_value(initial_conditions, t0).set_f_params(-1.0, 1.0, 1.0)
positions = []
t1 = 50
dt = 0.05
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
positions.append(r.y[:3]) # keeping only position, not velocity
positions = np.array(positions)
FRAMES = 50
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
def init():
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# animation function. This is called sequentially
def animate(i):
current_index = int(positions.shape[0] / FRAMES * i)
ax.cla()
ax.plot3D(positions[:current_index, 0],
positions[:current_index, 1],
positions[:current_index, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=FRAMES, interval=100)
# call our new function to display the animation
display_animation(anim)
We see that the particle turns in the other direction, compared to the previous case!
What happens when the field is stronger?¶
Another visualization can be done with the same exact initial conditions but with a stronger field.
r.set_initial_value(initial_conditions, t0).set_f_params(1.0, 1.0, 10.0)
positions2 = []
t1 = 50
dt = 0.05
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
positions2.append(r.y[:3]) # keeping only position, not velocity
positions2 = np.array(positions2)
Let's now visualize the previous trajectory, as well as the new one together:
FRAMES = 50
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
def init():
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# animation function. This is called sequentially
def animate(i):
current_index = int(positions.shape[0] / FRAMES * i)
ax.cla()
ax.plot3D(positions[:current_index, 0],
positions[:current_index, 1],
positions[:current_index, 2])
ax.plot3D(positions2[:current_index, 0],
positions2[:current_index, 1],
positions2[:current_index, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=FRAMES, interval=100)
# call our new function to display the animation
display_animation(anim)
When the field is stronger, the radius of the oscillation, called the Larmor radius, decreases. Also, the frequency of oscillation increases.
The formulas for these quantities are:
$$ \rho = \frac{m v_\perp}{|q| B} $$
$$ \omega_L = \frac{|q| B}{m} $$
Charged particles in constant electric and magnetic fields¶
Another phenomenon of interest is when charged particles are subject to a constant electric field and a constant magnetic field. The expected behaviour is that the electric field will introduce a drift, while the magnetic field will just make the particles loop around the field lines.
Here, we will consider that we have:
$$ \vec{B} = B \vec{e_x} $$
$$ \vec{E} = E \vec{e_y} $$
We need to redefine our integrators from the previous sections:
def newton(t, Y, q, m, B, E):
"""Computes the derivative of the state vector y according to the equation of motion:
Y is the state vector (x, y, z, u, v, w) === (position, velocity).
returns dY/dt.
"""
x, y, z = Y[0], Y[1], Y[2]
u, v, w = Y[3], Y[4], Y[5]
alpha = q / m
return np.array([u, v, w, 0, alpha * B* w + E, -alpha * B * v])
All other parameters stay the same.
r = ode(newton).set_integrator('dopri5')
r.set_initial_value(initial_conditions, t0).set_f_params(1.0, 1.0, 1.0, 10.)
positions = []
t1 = 50
dt = 0.05
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
positions.append(r.y[:3])
positions = np.array(positions)
FRAMES = 50
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
def init():
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# animation function. This is called sequentially
def animate(i):
current_index = int(positions.shape[0] / FRAMES * i)
ax.cla()
ax.plot3D(positions[:current_index, 0],
positions[:current_index, 1],
positions[:current_index, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=FRAMES, interval=100)
# call our new function to display the animation
display_animation(anim)
The expected drift is given by the cross product of $\vec{E}$ and $\vec{B}$. In our case this is $\vec{e_y} \times \vec{e_x}$ so as observed above, the drift goes in the negative z direction. It should be noted that the drift does not depend on the charge. This is due to the fact that when averaging the $v_\perp$ value over one gyration of the particle, one gets the following equation:
$$ < m \frac{d \vec{v_\perp}}{dt}> = \vec{0} = q \left( \vec{E_\perp} + <\vec{v_\perp}> \times \vec{B}\right) $$
In the above equation, $q$ gets simplified, so it has no influence on the average drift of the charged particle.
If instead of E there is an external force, such as gravity, the drift is again in the direction of the force cross the magnetic field. Another important case of drift is when the field lines of B are curved with nonconstant radius: in this case, the particles looping around the field lines of B will feel centrifugal force and will thus escape from the field line.
It should be noted that except for the $E \times B$ drift, all other drifts depend on the charge of the particle.
Final simulation: non-constant B and E field with many particles¶
To finish off this notebook, we will proceed to simulate many particles at the same time in a constant B field, but with a non-constant periodic E field.
First, we write a function to integrate the trajector of a single particle:
def e_of_x(x):
return 10 * np.sign(np.sin(2 * np.pi * x / 25))
def compute_trajectory(m, q):
r = ode(newton).set_integrator('dopri5')
r.set_initial_value(initial_conditions, t0).set_f_params(m, q, 1.0, 10.)
positions = []
t1 = 200
dt = 0.05
while r.successful() and r.t < t1:
r.set_f_params(m, q, 1.0, e_of_x(r.y[0]))
r.integrate(r.t+dt)
positions.append(r.y[:3])
return np.array(positions)
positions = []
for m, q in zip([1, 0.1, 1, 0.1], [1, 1, -1, -1]):
positions.append(compute_trajectory(m, q))
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for position in positions:
ax.plot3D(position[:, 0], position[:, 1], position[:, 2])
plt.xlabel('x')
plt.ylabel('y')
ax.set_zlabel('z')
<matplotlib.text.Text at 0x80a2278>
Let's do an animation with these particles:
FRAMES = 100
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
def init():
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# animation function. This is called sequentially
def animate(i):
current_index = int(positions[0].shape[0] / FRAMES * i)
ax.cla()
for position in positions:
ax.plot3D(position[:current_index, 0],
position[:current_index, 1],
position[:current_index, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=FRAMES, interval=100)
# call our new function to display the animation
display_animation(anim)
The movement is quite difficult to understand. Given the previous sections, can you guess why the 4 particles behave the way they do?
Conclusions¶
This exploration in magnetic and electric fields has enabled us to visualize some of the important properties that charged particles exhibit when in a magnetic and electric field. Hopefully, this has increased your physical intuition about these phenomena.
