Example: Particle Gyrating in Uniform Magnetic Field

This is a simple demonstration of using a BorisIntegrator to calculate the trajectory of a gyrating particle.

NOTE: Eventually, there will be a higher-level interface so one doesn’t have to deal with what specific integrator is used for tracing.

As always, let’s start by importing some libraries. matplotlib and numpy are pretty universal. scipy.constants is for physical constants. Finally ggcmpy.tracing is where the tracing functionality currently lives.

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np
from scipy import constants

import ggcmpy.tracing

Setting up a Uniform Magnetic Field

ggcmpy.tracing provides UniformField for this purpose, which is initialized with a constant \(\vec B_0\) and/or \(\vec E_0\). The UniformField just returns those magnetic and electric field independent of position.

B_0 = 1e-8  # [T]
fields = ggcmpy.tracing.UniformField(B_0=np.array([0.0, 0.0, B_0]))

Setting up the Boris Integrator

The BorisIntegrator needs to know charge and mass of the particle (an electron here), as well as the electromagnetic fields – we use the uniform fields from above.

The integration starts at initial position \(\vec x_0\) and initial velocity \(\vec v_0\). We also need to specify the end time \(t_{max}\) (chosen to be one gyro period), and for the time being the time step \(dt\) as well.

q = -constants.e
m = constants.m_e
x0 = np.array([0.0, 0.0, 0.0])  # [m]
v0 = np.array([0.0, 100.0, 0.0])  # [m/s]
om_ce = np.abs(q) * B_0 / m  # [rad/s]
t_max = 2 * np.pi / om_ce  # one gyroperiod # [s]

boris = ggcmpy.tracing.BorisIntegrator(fields, q, m)
df = boris.integrate(x0, v0, t_max)
# This limits the time step to 1% of the gyroperiod, ie., at least 100 steps per gyration
df2 = boris.integrate(x0, v0, t_max, gyro_max=0.01)

The integrate() function returns a Pandas DataFrame, which is essentially a table, where each row is the particle data at a given time step. The columns are time, particle position x, y, z, and particle velocity vx, vy, vz.

df

	time	x	y	z	vx	vy	vz
0	0.000000	0.000000	0.000000	0.0	0.000000	100.000000	0.0
1	0.000357	-0.010215	0.032515	0.0	-57.187658	82.033968	0.0
2	0.000714	-0.037189	0.053346	0.0	-93.826609	34.591437	0.0
3	0.001072	-0.071230	0.055009	0.0	-96.751722	-25.280512	0.0
4	0.001429	-0.100106	0.036907	0.0	-64.911944	-76.068650	0.0
5	0.001786	-0.113442	0.005542	0.0	-9.747964	-99.523752	0.0
6	0.002143	-0.106445	-0.027813	0.0	48.918661	-87.217914	0.0
7	0.002501	-0.081630	-0.051175	0.0	90.007801	-43.572879	0.0
8	0.002858	-0.047913	-0.056149	0.0	98.755279	15.728791	0.0
9	0.003215	-0.017410	-0.040947	0.0	72.017946	69.378782	0.0
10	0.003572	-0.001081	-0.011032	0.0	19.403078	98.099544	0.0

Plotting

Pandas has built-in plotting functionality – here is the particle position in the \(x-y\) plane, which as expected is a circle.

The columns are stored in numpy arrays, so one can of course just plot them by hand with matplotlib or whatever, too.

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
df.plot(x="x", y="y", c="time", kind="scatter", ax=axs[0], title="Default time step")
df2.plot(
    x="x", y="y", c="time", kind="scatter", ax=axs[1], title="Smaller gyro time step"
)
axs[0].set_aspect("equal")
axs[1].set_aspect("equal")

Energy Conservation

The Boris pusher makes sure that the magnetic Lorentz force just changes the direction of the particle, but keeps the speed exactly (to machine precision) the same. So in this case, where there is no electric field, kinetic energy of the particle should be exactly conserved – and we can check that.

The plot below indeed shows constant kinetic energy.

df["E"] = 0.5 * constants.m_e * np.linalg.norm(df[["vx", "vy", "vz"]].values, axis=1)
df.plot(x="time", y="E", title="Energy of the particle over time");