Simulation of Charged Particle Orbits in
Fusion Plasmas
Bachelor's thesis in Engineering Physics

MATHIAS HOPPE
AYLWIN IANTCHENKO
INGRID STRANDBERG

Department of Applied Physics
CHALMERS UNIVERSITY OF TECHNOLOGY
Göteborg, Sweden 2015





Bachelor’s thesis in Engineering Physics

Simulation of Charged Particle Orbits
in Fusion Plasmas

Mathias Hoppe
Aylwin Iantchenko
Ingrid Strandberg

Department of Applied Physics
CHALMERS UNIVERSITY OF TECHNOLOGY

Göteborg, Sweden 2015



SIMULATION OF CHARGED PARTICLE ORBITS IN FUSION PLASMAS

Mathias Hoppe
Aylwin Iantchenko
Ingrid Strandberg

© Mathias Hoppe, Aylwin Iantchenko, Ingrid Strandberg, 2015

Bachelor’s thesis TIFX02-15-41
Department of Applied Physics
Chalmers University of Technology
SE-412 96 Göteborg
Sweden
Telephone: +46 (0)31-772 1000

Cover:
A 3D model of the International Thermonuclear Experimental Reactor (ITER) with
charged particle orbits, each of duration 15µs, plotted in different colors. The orbits
were calculated using the tools developed for this project.

Chalmers Reproservice
Göteborg, Sweden 2015



SIMULATION OF CHARGED PARTICLE ORBITS IN FUSION PLASMAS

Mathias Hoppe
Aylwin Iantchenko
Ingrid Strandberg
Department of Applied Physics
Chalmers University of Technology

Abstract
When designing a fusion device, knowledge of the particle motion inside the fusion
plasma is crucial. The charged plasma particles are confined inside the device using
a strong magnetic field, which influences particle motion. Particle trajectories can
therefore be obtained by numerically solving the equations of motion for a charged
particle in the confining magnetic field.

With the simulation tool developed as a part of this project, charged particle
orbits are studied. Especially, the properties of the so called banana and passing
orbit topologies are studied and the observed results explained using theoretical
models. We find expressions that approximately describe the width of the banana
and passing orbits and the location of the banana orbit’s mirror points.

The orbit dependence on mass, charge and energy is investigated and an expres-
sion for the particle’s deviation from a field line is derived. Also, the cause for
banana orbits forming is studied and their occurrence is shown to depend on how
the particle’s velocity vector is directed. Finally, the two computational methods
used, where either the particle or its guiding-center is followed, are compared with
respect to both energy conservation and computational time. The guiding-center
approach is shown to greatly reduce computational cost.

Sammandrag
När en fusionsanordning ska designas krävs kunskap om hur partiklarna som utgör
fusionsplasmat rör sig. De laddade plasmapartiklarna hålls instängda i fusionsanord-
ningen med hjälp av starka magnetfält som påverkar partiklarnas rörelse. Partikel-
banorna kan därför beräknas genom att ställa upp och numeriskt lösa rörelsekva-
tionerna för laddade partiklar i magnetfält.

Partikelbanor studeras med hjälp av det simuleringsverktyg som speciellt utveck-
lats för detta projekt. I synnerhet studeras egenskaperna hos så kallade banan-
och övergångsbanor. Med hjälp av teoretiska modeller förklaras de gjorda observa-
tionerna och uttryck för bland annat banornas bredd samt läget för banan-banans
spegelpunkter tas fram.

Banans beroende av massa, laddning och energi undersöks och ett uttryck för
partikelns avvikelse från en fältlinje härleds. Även orsaken till att banan-banor upp-
står studeras och deras uppkomst visar sig bero på hur partikelns hastighetsvektor
är riktad. Slutligen jämförs de två beräkningsmetoderna som används, där antingen
partikeln eller dess bancentra följs, med avseende på energikonservering och beräkn-
ingshastighet. Att följa bancentrat visar sig vara beräkningsmässigt effektivare.

i



Acknowledgements
We would first like to thank our supervisors: Eero Hirvijoki, for tirelessly and with
dedication teaching us the theory of the subject of this thesis, as well as pushing us to
do our very best, and István Pusztai, for constantly giving insightful and enlightening
comments, even when commitments outside of Gothenburg and Sweden kept him
busy. Thank you both for always being at hand to give answers whenever we had
questions.

We would also like to thank Professor Tünde Fülöp, for being so kind and making
us feel very welcome into the electromagnetic field theory research group.

The Authors, Göteborg, May 2015

ii



Contents

1 Introduction 1
1.1 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Analytical theory of particle motion 4
2.1 Particle Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1 Particle motion in a homogeneous magnetic field . . . . . . . . . 4
2.1.2 Particle drifts in inhomogeneous magnetic fields . . . . . . . . . 6

2.2 Guiding-center motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Derivation of the guiding-center Lagrangian . . . . . . . . . . . . 10
2.2.2 Derivation of the guiding-center equations of motion . . . . . . . 13

3 Numerical methods 16
3.1 Magnetic field data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Time integration method . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Domain check . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Program workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4 Simulation of particle orbits 24
4.1 Observed orbit topologies . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.1.1 Orbit width . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.2 Deviation from field lines for banana orbits . . . . . . . . . . . . 30
4.1.3 Banana mirror points . . . . . . . . . . . . . . . . . . . . . . . . 33

4.2 Explanation of orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.2.1 Effect of mass, charge and energy . . . . . . . . . . . . . . . . . 37
4.2.2 The mirroring effect . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2.3 Mirror point dependence of ξ . . . . . . . . . . . . . . . . . . . . 41

4.3 Comparison between simulation methods . . . . . . . . . . . . . . . . . 43
4.3.1 Energy Conservation . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.2 Computational time . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.3.3 Particle trajectories . . . . . . . . . . . . . . . . . . . . . . . . . 46

5 Summary and Discussion 50

References 52

A Derivation of the charged particle Lagrangian 54

B Gauge invariance of the Lagrangian 56

C Magnetic moment: an adiabatic invariant 57

D Running the program 59
D.1 Downloading and compiling . . . . . . . . . . . . . . . . . . . . . . . . 59
D.2 Command line arguments . . . . . . . . . . . . . . . . . . . . . . . . . 59
D.3 Parameters from file . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
D.4 Output data format . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

iii



“You do not follow the path of a hummingbird by looking at every flap of its wings”

Alain J. Brizard (paraphrased)

iv



1 Introduction
Together with increasing electricity demands and climate change, the need for an
environmentally friendly substitute to fossil fuel is necessary. One of the most
promising long term candidates is nuclear fusion, which generates energy by fusing
lighter nuclei into heavier ones. The energy released during fusion comes from the
difference in the binding energies of the initial and final states. The goal of controlled
fusion research is to extract this energy.

In order to induce these reactions, very high temperatures are needed. These
processes occur naturally in stars, including our Sun, in which temperatures are of
the order 10 million K [1]. In order to produce nuclear fusion on Earth however,
temperatures of the order 100 million K are required [2]. This is because we use a
different fusion reaction than the dominant processes in the Sun. Additionally, there
is a tremendously high pressure in the core of the Sun which cannot be produced in
fusion devices. As fuel used for fusion devices must be heated to extreme temper-
atures, a substantial problem becomes apparent: there are no materials capable of
withstanding such high temperatures.

A star is held together by its own gravitational field. On Earth we provide
confinement by using a property of the fusion fuel which arises at high temperatures.
As the temperature rises in any material, atoms start to move faster. This movement
gives rise to the different states of matter, such as solid, liquid and gaseous. At high
enough temperatures atoms start colliding and cause electrons to be released, thus
creating an ionized gas state called a plasma state. The property of being ionized
makes the plasma susceptible to the effects of electric and magnetic fields.

Applying a straight magnetic field in a plasma will make the constituent charged
particles move in a helical orbit around the magnetic field lines. This straight mag-
netic field is sufficient for containing charged particle motion that is perpendicular
to the field lines, but if a particle has a velocity component parallel with the field
lines, it will drift towards the front and back walls of the device. Having charged
particles repeatedly collide with the device walls would cause severe damage, as well
as lead to an unacceptably high rate of energy loss.

A better approach is instead to bend the magnetic field into a torus shape, so that
particles are lead in a circular orbit inside the device. This is the approach taken
today in fusion research and it is utilized in the most successful types of fusion
devices: tokamaks [3] and stellarators [4]. There are several operational differences
between these two kinds of devices, but the idea of a bent torus-shaped magnetic
field is used in both. In this project, magnetic equilibrium data from the tokamak
International Thermonuclear Experimental Reactor (ITER) will be used. A model
of ITER can be seen in Fig. 1.

1



1.1 Purpose

Figure 1: The outline of the inner walls of the International Thermonuclear Experimental
Reactor.

The confinement of the fusion plasma comes with additional difficulties to those
just mentioned. The hot plasma contains high energy free electrons and ions, so one
of the difficulties is to trap the high energy particles in such a way that they deposit
their energy inside the core of the plasma. This is of paramount importance in order
to maintain the high temperature, and also to prevent highly energetic particles
from hitting the device walls, and thus causing damage to them. Accordingly, the
magnetic geometry has to be meticulously optimized for this to be done successfully.
Because of this, knowledge of the particle motion inside the plasma is necessary when
designing a fusion device. To achieve this, the equations of motion for a charged
particle in a magnetic field need to be solved. This has been the main focus of the
project.

1.1 Purpose
In this project, the trajectories of charged particles inside a tokamak fusion device
will be simulated and studied. Simulations will be done using two different methods,
except for when performance restrictions prevent the use of either of the methods.
Observations made in the simulations will then be explained using theoretical models,
and the two methods of simulation will be compared in terms of performance and
agreement. Simulations will be performed by numerically solving the equations of
motion applicable to the system, and in order to do this a specialized simulation
tool must be developed.

1.2 Limitations
The goal is to follow collisionless, single particle orbits within a given static, inhomo-
geneous magnetic field. The collisionless description means neglecting the possibility
of collision between particles. To take collisions into account stochastic methods are
necessary, which is outside the scope of this project. We restrict the problem to de-
terministic (time-reversible) processes. For high energy particles collisionless theory
is often well justified, since the collision frequency is proportional to T−3/2, where T
is the kinetic energy of the particle [5]. The frequency of collisions is typically much
lower than the frequency of particle transit in the device.

2



1.3 Method

Time-varying fields will not be considered. This is because in real tokamaks
the confining magnetic field varies slowly enough to be considered constant for our
application. There are fluctuating fields as well due to turbulence and instabilities,
but that is also outside the scope of the thesis.

In the plasma, each charged particle contributes to the electric and magnetic
fields. If there are many particles they can significantly affect the fields. We will
not consider this “feedback”, which essentially means that only a small number
of particles will be taken into account. When particle motion in real tokamak
geometry is studied, we will focus on the effects of the magnetic field, excluding
electric fields. Furthermore, particles with high enough energy that they require
a relativistic description will not be covered. The exclusion of electric fields and
relativistic particles is mainly due to that including these phenomena would make the
calculations more complex, yet it would not add essentially to the understanding.

1.3 Method
In this project both analytical mechanics and numerical methods will be used to
simulate and visualize the motion of charged particles in the inhomogeneous mag-
netic field of a tokamak device. Given the magnetic vector field, the initial position
and velocity of the charged particle as well as its mass and charge, the particle’s
trajectory inside the reactor will be traced.

First, in Section 2, the theory for charged particle motion in electromagnetic fields
will be reviewed. We will begin in Section 2.1 by solving the equations of motion
in a simple, straight, homogeneous magnetic field. This gives the well-known result
that the particle gyrates around a field line as it moves along it. After this, we
solve for an inhomogeneous magnetic field, discovering particle drifts. Since the
gyrations are not of interest on a larger scale, we wish to decouple the gyro-motion
from the non-trivial guiding-center motion. Using methods of analytical mechanics
in Section 2.2, a Lagrangian for the guiding-center and its following equations of
motion are derived.

Because the resulting equations of motion cannot be solved analytically, a nu-
merical ordinary differential equation (ODE) solver was written using the C pro-
gramming language. In Section 3 the magnetic field data is described, the included
algorithms are explained, and the overall program workflow is presented.

In Section 4, the simulation results found using the tool developed are presented.
In Section 4.1, the trajectories for different kinds of particles are shown in figures.
An analysis of the observed particle motion follows in Section 4.2. Finally, we
compare the results obtained from calculations of the particle motion versus the
guiding-center motion. The performance of the numerical solver in the two different
cases is also compared. It becomes clear why the guiding-center approach is needed:
it reduces the computational cost tremendously.

3



2 Analytical theory of particle motion
Since the plasma consists of a large number of charged particles such as electrons,
hydrogen ions and alpha particles, each of these particles will be affected by a force
from the electric and magnetic field. This force is the well-known Lorentz force and
it takes the form

F = q (E + v ×B) , (1)
where q is the particle charge, E the electric field, v the particle velocity and B the
magnetic field.

Since the goal of this project is to study the motion of various charged particles
within the fusion plasma, we make several assumptions that should not affect the
qualitative conclusions drawn from the later analysis. First we assume that there
are no other forces but the Lorentz force acting on the particles. This is valid
since gravitational effects are clearly negligible due to the ratio mg/qvB (which is
a rough estimate of gravitational effects compared to magnetic effects) being of the
order 10−10 at most, for an alpha particle.

To simplify the treatment of the problem, we neglect electric fields, though these
can be added in a straightforward way if desired. The total force acting on a plasma
particle is then

F = qv ×B.

Using Newton’s second law we can now write the equations of motion for a charged
particle in a fusion plasma as {

v̇ = q
m
v ×B,

ṙ = v.
(2)

These equations will be solved analytically for special cases in the following two sub-
sections. First, in Section 2.1.1, they are solved in a static, homogeneous magnetic
field after which the resulting expressions for the particle motion is discussed. Then,
in Section 2.1.2, the particle motion is studied in a more general, inhomogeneous
magnetic field.

To decrease computational cost in numerical calculations, it is of interest to de-
couple the gyro-motion of the particle from its guiding-center motion. This is done
using the guiding-center transformation, and in Section 2.2 the Lagrangian and its
associated equations of motion for the guiding-center will be derived.

2.1 Particle Motion
We want to gain an understanding of the particle motion in strongly magnetized
plasmas. The magnetic fields present in fusion devices are quite complicated and
require numerical methods for tracing the charged particles, but starting out with
a simple magnetic field will give us a qualitative understanding of what happens to
the particle. We develop this step by step. The simplest possible set-up is a straight
homogeneous field, which will be examined first.

2.1.1 Particle motion in a homogeneous magnetic field
For a particle in a straight static magnetic field, B = Bx̂, where B is the magnetic
field strength (which is assumed constant for now), the equations of motion (2) can

4



2.1 Particle Motion

be written as

v̇x = 0, (3)

v̇y =
q

m
vzB, (4)

v̇z = − q

m
vyB, (5)

ṙ = v. (6)

With the initial conditions r0 = (x0, y0, z0) and v0 = (vx0, vy0, vz0) the first of these
equations is easily solved to give

x(t) = vx0t+ x0.

To solve for y and z, we begin by differentiating both (4) and (5) with respect to
time in order to get

v̈y =
q

m
v̇zB,

v̈z = − q

m
v̇yB.

Substituting the expressions for v̇y and v̇z into these equations, defining Ω = |q|B/m
and reordering the terms, gives us two second order differential equations of the
familiar forms

v̈y + Ω2vy = 0,

v̈z + Ω2vz = 0.

These equations have periodic solutions that can be written as

vy = Re
[
CeiΩt

]
, (7)

vz = Re
[
DeiΩt

]
, (8)

where C and D are complex integration constants. Using Eq. (4) we get the relation
between C and D as C = i sgn(q)D, where sgn(q) denotes the sign of the charge q.
This essentially means that vy and vz are offset from each other by a phase angle of
90◦, and allows us to determine |C| by adding the squares of vy and vz. If we define
the perpendicular velocity to be v⊥ ≡ (v2y + v2z)

1/2, we get(
vy(t)

)2
+
(
vz(t)

)2
= |C|2(cos2 Ωt+ sin2Ωt) = v2⊥ ≡

√
v2y0 + v2z0,

where the last identity holds since C is a constant. Note that this means v⊥ will
remain constant. Since vx is also constant, the kinetic energy is constant: the
magnetic field does no work on the particle.

Because C is a complex constant, the above is not enough to fully determine its
value. We can however write C in polar form, offset by some phase angle ϕ, as
C = v⊥e

iϕ. Dividing vy(0) = vy0 by vz(0) = vz0 then gives

ϕ = − sgn(q) arctan
(
vz0
vy0

)
.

5



2.1 Particle Motion

Finally, by integrating (7) and (8) in time and using the initial conditions y(0) = y0
and z(0) = z0 we arrive at the expressions determining the particle position in the
y-z plane,

y(t) = ygc − Re
[
i
v⊥
Ω
ei(Ωt+ϕ)

]
,

z(t) = zgc + sgn(q)Re
[
v⊥
Ω
ei(Ωt+ϕ)

]
,

(9)

where ygc = y0 + (v⊥/Ω)Re(eiϕ) and zgc = z0 − (v⊥/Ω) sgn(q)Re(eiϕ). These equa-
tions carry some very useful information. The first thing we notice is that the motion
of the particle in the magnetic field is in the shape of a helix around the magnetic
field line, illustrated in Fig. 2. As we will see, the qualitative behavior remains
similar even for more complicated magnetic fields. It also motivates the important
guiding-center transformation which will be discussed later.

Figure 2: The motion of an electron in a straight magnetic field. A magnetic field line
has been visualized for reference.

From the expressions above we also notice the importance of the three quantities
ygc, zgc and v⊥/Ω ≡ ρ. The first two denote the position of the center of the
gyromotion made by the particle. This point, (ygc, zgc), is often referred to as the
guiding-center or gyro-center. The third quantity is the gyration radius, commonly
referred to as the Larmor radius and denoted ρ, which determines the size of the
helical orbit in the plane perpendicular to the field line. With these quantities it is
useful to define

x = X + ρ, (10)
where x is the particle position, X is the guiding-center and ρ is the gyro-radius
vector. This will be used in following sections, starting with 2.1.2.

2.1.2 Particle drifts in inhomogeneous magnetic fields
From the theory of straight magnetic fields above, we can easily conclude that a
particle with even a very small velocity component parallel to the magnetic field
would cause the particle to follow the field line until some physical boundary stops
the particle. Therefore, in order to contain the particle for a sufficiently long time,
a device of considerable length would be required. This is obviously not practical,

6



2.1 Particle Motion

and so instead the device is bent into a torus (see Fig. 1). This bending of the
device however, also requires the magnetic field to be curved along with it, making
the forces acting on the particle different. When the magnetic field is curved, its
strength should have a spatial variation due to its divergence-free nature, which will
cause yet another force acting on the particle.

In this project, knowledge about particle drifts will be needed later on when
particle orbits are analyzed, though we will mostly need some basic results of this
theory. Because of this, the theory derived in this section is only derived briefly,
with several details left out. A more exhaustive analysis of particle drifts can be
found for example in the book by Chen (1984) [6].

Constant force
In order to simplify later calculations, let us first study what happens to the particle
when we apply a general force F = Fxx̂+Fyŷ+Fzẑ. Keeping the straight magnetic
field B = Bx̂ for now, we add the force F to the equations of motion (2),

v̇x =
1

m
Fx,

v̇y =
q

m
vzB +

1

m
Fy,

v̇z = − q

m
vyB +

1

m
Fz.

By requiring F to be constant, we can solve these equations rather easily in the
same way we solved (2) earlier. This yields the solutions

vx =
Fx

m
t+ vx0,

vy = Re
[
v⊥e

iΩt+ϕ
]
+

Fz

qB
,

vz = sgn(q)Re
[
iv⊥e

iΩt+ϕ
]
− Fy

qB
.

The motion in the x-direction (parallel to the magnetic field lines) has now become
accelerating, just as expected from Newtons second law. This component will not
be of great concern to us, as it is along the magnetic field line, and so will not
cause the particle to deviate from the field lines. However, we also found that in
the y- and z-directions a velocity perpendicular to the field lines is induced. It is
interesting to note that this velocity causes the particle’s guiding-center to gradually
move across the magnetic field, or in other words, causes it to drift. The velocity
vgc = (1/qB)(Fzŷ − Fyẑ) is therefore referred to as a drift velocity.

It is possible to obtain a general formula for the drift velocity, vgc, which holds
also for general forces F and magnetic fields B = B(x, y, z, t)b̂. To obtain this
formula, we start from the vector form of the equations of motion,

mv̇ = qv ×B + F .

Since the mv̇ term only gives the gyro-motion, which we already know about, we can
omit it. We form the vector product of the resulting equation with B and rewrite

7



2.1 Particle Motion

the expression a bit,

q(v ×B)×B + F ×B = 0 ⇐⇒
q(B2v − v||B

2b̂) = F ×B,
(11)

where v|| = v · b̂. We assume now that the velocity v can be decomposed into one
part parallel with the magnetic field, v|| = v||b̂, and one part perpendicular to the
magnetic field, vF . The left-hand side of Eq. (11) then simply becomes qvF , and
we can write an expression for the drift velocity,

vF =
F ×B

qB2
. (12)

As mentioned above, this formula is general and applies to any force F and magnetic
field B.

∇B drift
Let us study the magnetic field B = B(y, z)x̂, which only varies in the ŷ and ẑ
directions, so that ∇B ⊥ B. It is not possible to analytically find an exact solution
for B of a general form. However, since the Larmor radius is much smaller than
the typical length scale over which B varies, B can be expanded in a Taylor series
around the particle orbit’s guiding-center, to get an approximate expression for the
drift velocity. Denote the magnetic field scale length LB = B/|∇B|. The condition
is then ρ ≪ LB, which is well satisfied in fusion devices for thermal particles. Since
we are only interested in the drift velocity here, the velocity component in the x
direction will be ignored.

The expansion for our magnetic field around the guiding-center is

B =

[
B(ygc, zgc) + (y − ygc)

∂B

∂y
+ (z − zgc)

∂B

∂z
+O(ρ2)

]
x̂.

Using Equations (4) and (5), and the expansion above to first order, we get the new
equations

mv̇y = Fy = qvz

[
B(ygc) + (y − ygc)

∂B

∂y
+ (z − zgc)

∂B

∂z

]
,

mv̇z = Fz = −qvy

[
B(ygc) + (y − ygc)

∂B

∂y
+ (z − zgc)

∂B

∂z

]
.

Since we are only interested in the drifts of the guiding-center, not the actual gy-
rations, we average the force over a gyration period. From Eq. (9) we easily find
the gyration period time to be τ = 2π/Ω. Denoting the average force as ⟨F ⟩, and
substituting the expressions for vy and vz from (7) and (8) to get an approximate

8



2.1 Particle Motion

expression, we find

⟨Fy⟩ =
1

τ

∫ τ

0

−|q|v⊥ sin(Ωt+ ϕ)

[
B(ygc) + ρ sin(Ωt+ ϕ)

∂B

∂y
+

+ ρ sgn(q) cos(Ωt+ ϕ)
∂B

∂z

]
dt = −|q|v⊥ρ

2

∂B

∂y
,

⟨Fz⟩ =
1

τ

∫ τ

0

−qv⊥ cos(Ωt+ ϕ)

[
B(ygc) + ρ sin(Ωt+ ϕ)

∂B

∂y
+

+ ρ sgn(q) cos(Ωt+ ϕ)
∂B

∂z

]
dt = −|q|v⊥ρ

2

∂B

∂z
.

We can now set the force F∇B = (0, ⟨Fy⟩, ⟨Fz⟩) and use Eq. (12) to find an expression
for the drift velocity. Noting that the choice for the orientation of B was entirely
arbitrary, as long as the gradient is perpendicular to it, we can find the general
expression for the drift velocity due to a perpendicular gradient

v∇B =
mv2⊥
2qB2

(b̂×∇B). (13)

This equation predicts that the guiding-center drift velocity of the particle will be
in a direction perpendicular to both b̂ and ∇B.

A physical picture for the ∇B drift is shown in Fig. 3, and follows from the fact
that the local radius of curvature of the gyro-orbit is smaller on the side of the orbit
with a larger magnetic field, and correspondingly, the radius is larger on the side
with the smaller magnetic field. If the trajectories are calculated and plotted for
such orbits, a net drift perpendicular to both b̂ and ∇B can be seen [5].

Figure 3: Ion ∇B drift motion. The gradient is directed upwards, and the dot indicates
that the magnetic field points outwards from the page. The resulting guiding-center drift
vgc is leftward, perpendicular to both ∇B and B.

Curvature drift
Apart from the ∇B drift, drift motion due to the curved geometry of the field will
also arise. This drift motion comes from the centripetal force acting on the particle

9



2.2 Guiding-center motion

which can be written as
F c =

mv2∥
Rc

r̂ =
mv2∥Rc

R2
c

, (14)

where Rc is the local curvature radius of the field. We also assume the field strength
B to be locally constant. Inserting (14) into (12) gives us the curvature drift velocity

vc =
F ×B

qB2
=

mv2∥
qB2

Rc ×B

R2
c

.

The radius of curvature can be written as Rc/R
2
c = −(b̂ · ∇)b̂. Define b̂ · ∇ = ∇∥,

the gradient along b̂. The curvature drift velocity can now be expressed as

vc =
mv2∥
qB

b̂×∇∥b̂. (15)

In this form it is evident that the curvature drift is caused by parallel gradients in
b̂ (with a gradient scale length L ≫ ρ).

2.2 Guiding-center motion
In the previous sections the particle motion has been described, showing the small
gyration together with the overall guiding-center motion. However, in some sit-
uations it is desirable to focus only on the average displacement of the particle,
neglecting the gyration. It is thus of interest to find how to separate these two
different components of the motion, and only solve for the guiding-center motion.
Because we are interested in how the particle moves on a much larger time scale
than the gyrofrequency, and a much larger length scale than the Larmor radius, a
Lagrangian for the guiding-center can be developed, from which the new equations
of motion can be found.

2.2.1 Derivation of the guiding-center Lagrangian
The regular Lagrangian for a charged particle is

L = T − U =
1

2
mv2 − qϕ+ qv ·A. (16)

For a complete derivation of this Lagrangian, see Appendix A. The idea is now to
derive a Lagrangian for which the gyro-motion is separated from the guiding-center
motion. This will be done on the basis of slowly varying fields and a small Larmor
radius. For a more thorough derivation we refer too the papers by Cary & Brizard
(2009) [7] and Littlejohn (1983) [8]. The derivation follows the steps done in both
of these papers, with some smaller modifications.

The Lagrangian we have right now is the same as that shown in Eq. (16), except
we do not keep the electric potential term. It can be rewritten as

L = [qA+mv] · ẋ− 1

2
mv2, (17)

where v = |v|. We want to find coordinates that transform the Lagrangian into the
desired form. As introduced in Eq. (10), let the particle position be

x = X + ρ,

10



2.2 Guiding-center motion

where X is the guiding-center position and ρ the Larmor vector. This has the total
time derivative

ẋ = Ẋ + ρ̇.

Let us now insert the new coordinates into our Lagrangian (17). Also, expand A(x)
in a Taylor series around the guiding-center position X.

L = q[A+ ρ · ∇A+O(ρ2)] · (Ẋ + ρ̇) +mv · (Ẋ + ρ̇)
1

2
mv2 =

= qA · Ẋ + qA · ρ̇+ qρ · ∇A · Ẋ + qρ · ∇A · ρ̇+mv · (Ẋ + ρ̇)
1

2
mv2 +O(ρ2).

(18)

Note that the vector potential A is now evaluated at the guiding-center position,
that is A = A(X). Moving onwards, we will use the identities

∇A · Ẋ = Ẋ ×B + Ẋ · ∇A,

A · ρ̇ =
d
dt (A · ρ)− Ẋ · ∇A · ρ,

ρ · ∇A · ρ̇ =
1

2
(ρ× ρ̇ ·B) +

1

2

d
dt (ρ · ∇A · ρ̇)− 1

2
ρ · (Ẋ · ∇∇A) · ρ︸ ︷︷ ︸

O(ρ2)

.

(19)

When substituting Equations (19) into the Lagrangian (18), one term will be can-
celed. To continue, we will use the gauge invariance of the Lagrangian, namely that
the gauge transformation

L → L+
dF

dt
,

will not affect the equations of motion. A proof of this statement is given in Ap-
pendix B. This means that the two total time derivatives from (19) can be removed.
We are then left with

L = qA · Ẋ + qρ · (Ẋ ×B) +
q

2
(ρ× ρ̇ ·B) +mv · (Ẋ + ρ̇)− 1

2
mv2 +O(ρ2). (20)

We continue by decomposing the velocity v into

v = v∥b̂+ v⊥ĉ, (21)

where the unit vector b̂ is in the direction of the magnetic field lines. The perpendic-
ular unit vector ĉ is rotating with the angular velocity ζ̇ = Ω. The vector directions
are illustrated in Fig. 4. There is also a unit vector â, â = b̂× ĉ in the gyro-vector
direction. The gyro-angle is denoted by ζ. Expressions for â and ĉ in terms of the
gyro-angle ζ and a pair of fixed-frame, orthogonal unit vectors ê1 and ê2 satisfying
b̂ = ê1 × ê2 is shown in Eq. (22).

11



2.2 Guiding-center motion

ê1
ζ

b̂

â

ĉ

ê2

Figure 4: Coordinate unit vectors. A fixed frame is composed by ê1 and ê2. The vector
b̂ is pointing in the direction of the magnetic field lines. The vector ĉ is rotating with the
angular velocity ζ̇ = Ω. The vector â = b̂× ĉ is pointing in the gyrovector direction.

â = cos(ζ)ê1 − sin(ζ)ê2,

ĉ = − sin(ζ)ê1 − cos(ζ)ê2.
(22)

We also note

ρ =
v⊥
Ω
â,

ρ̇ =
d

dt

(
v⊥
Ω
â

)
=

v⊥
Ω

d

dt
(â) =

v⊥
Ω
ζ̇ĉ.

(23)

Insert (21) and (23) into the Lagrangian (20). Using the triple product rule and
perpendicularity between vectors, after a few algebraic manipulations, we find the
following Lagrangian:

L =
[
qA+mv∥b̂

]
· Ẋ +

m2v2⊥
2qB

ζ̇ − 1

2
mv2.

Furthermore, let J = m2v2⊥/2qB and note that mv2/2 = mv2∥/2 + µB, where µ =

mv2⊥/2B is the magnetic moment. The final Lagrangian becomes

Lgc(X, v∥, µ, ζ; t) =
[
qA+mv∥b̂

]
· Ẋ + Jζ̇ − m

2
v2∥ − µB. (24)

Notice that we have carried out a transformation from the Lagrangian depending
on the phase space coordinates (x,v; t) to the guiding-center Lagrangian depending
on the new coordinates (X, v∥, µ, ζ; t).

As we can see, the guiding-center Lagrangian does not contain the coordinate ζ,
since it is a cyclic coordinate. From here it follows that ∂Lgc/∂ζ̇ is a constant of
motion [9]. Evaluating this derivative gives

∂Lgc

∂ζ̇
= J =

m

q
µ, (25)

12



2.2 Guiding-center motion

thus the magnetic moment can be regarded as a constant of motion for the guiding-
center. Keep in mind that µ is really an adiabatic invariant, so it is only a constant
of motion for the particle to first order approximation. A more in-depth look on this
together with another derivation of the invariance of µ can be found in Appendix C.

The kinetic energy of the guiding-center motion is conserved in this case. This is
due to the lack of an electric field. An energy gain is otherwise expected from work
done by the electric field on the guiding-center. Including a time dependence of the
magnetic field also adds to the energy. More details and calculations of this can be
found in Northrop (1963) [10].

2.2.2 Derivation of the guiding-center equations of motion
The equations of motion are obtained from the Lagrangian (24) by applying the
Euler-Lagrange equations for each of the phase-space coordinates X, v∥, µ and ζ.
Starting with ζ we get that

∂Lgc

∂ζ̇
=

m

q
µ,

∂Lgc

∂ζ
= 0,

and thus the equation of motion associated with ζ becomes

d
dt

( ∂Lgc

∂ζ̇

)
− ∂Lgc

∂ζ
= 0 =⇒ µ̇ = 0,

which indicates that µ is constant as stated earlier.
Moving on to the next parameter, µ, we get that

∂Lgc

∂µ̇
= 0,

∂Lgc

∂µ
=

m

q
ζ̇ −B,

so the equation of motion associated with µ is

ζ̇ =
qB

m
= Ω.

Next, the derivatives of Lgc with respect to v∥ and v̇∥ are

∂Lgc

∂v̇∥
= 0,

∂Lgc

∂v∥
= mb̂ · Ẋ −mv∥,

and from the Euler-Lagrange equations we find

v∥ = b̂ · Ẋ,

which states that v∥ is the parallel velocity as was mentioned at the beginning.
We continue now with the derivation of the equation of motion for X. We begin

with the first partial derivative which gives

∂Lgc

∂X
= q∇(A · Ẋ) +mv∥∇(b̂ · Ẋ)− µ∇B =

= Ẋ ×
[
qB +mv∥(∇× b̂)

]
+ Ẋ · ∇

(
qA+mv∥b̂

)
− µ∇B,

13



2.2 Guiding-center motion

where we used the vector identity

∇
(
A · Ẋ

)
=
(
A · ∇

)
Ẋ +

(
Ẋ · ∇

)
A+ Ẋ ×

(
∇×A

)
+A×

(
∇× Ẋ

)
,

and the fact that the operator ∇ acting on Ẋ is zero because Ẋ is not explicitly
position dependent. The next derivative is

∂Lgc

∂Ẋ
= qA+mv∥b̂,

thus

d
dt

( ∂Lgc

∂Ẋ

)
= q∇A · Ẋ +mv̇∥b̂+mv∥∇b̂ · Ẋ.

Using the Euler-Lagrange equations

d
dt

∂Lgc

∂Ẋ
=

∂Lgc

∂X
,

we get the equations of motion

mv̇∥b̂ = Ẋ ×
[
qB +mv∥(∇× b̂)

]
− µ∇B.

We have now arrived at the following equations of motion:

µ̇ = 0, ζ̇ = Ω, v∥ = b̂ · Ẋ,

mv̇∥b̂ = Ẋ ×
[
qB +mv∥(∇× b̂)

]
− µ∇B.

(26)

Let us now manipulate the last one of these. By defining the effective magnetic field

B∗ = ∇×
(
A+

mv∥
q

b̂

)
= B +∇×

mv∥
q

b̂, (27)

the last equation of Eq. (26) can be rewritten as

mv̇∥b̂ = −µ∇B + qẊ ×B∗. (28)

It is possible to solve Eq. (28) in terms of Ẋ by doing a cross product of each term
with b̂. Using v∥ = b̂ · Ẋ we get

Ẋ = v∥
B∗

B∗
||
− µ

q
∇B × b̂

B∗
||
, (29)

where B∗
|| = b̂ ·B∗. Now, using the identity

(b̂×B∗)× b̂ = (b̂ · b̂)B∗ − (B∗ · b̂)b̂,

we can write
B∗ = (B∗ · b̂)︸ ︷︷ ︸

B∗
||

b̂+ (b̂×B∗)× b̂.

14



2.2 Guiding-center motion

In the second term B∗, insert Eq. (27). This gives

B∗ = B∗
||b̂+

mv∥
q

[
b̂×

(
∇× b̂

)]
× b̂.

Inserting the result from the following identity

∇(b̂ · b̂)︸ ︷︷ ︸
= 0

= 2(b̂ · ∇)b̂+ 2b̂× (∇× b̂) =⇒ b̂× (∇× b̂) = (b̂ · ∇)b̂,

gives
B∗ = B∗

||b̂−
mv∥
q

b̂ · ∇︸ ︷︷ ︸
∇∥

b̂× b̂.

Inserting this and µ = mv2⊥/2B into Eq. (29) we finally arrive at

Ẋ = v∥b̂+
mb̂

qB∗
||
×

(
v2⊥
2B

∇B + v2∥∇∥b̂

)
. (30)

In this equation for the guiding-center, the first term, v∥b̂, gives the motion along the
field lines while the second term gives the guiding-center drifts and can be compared
with the drift equations (13) and (15).

In a similar way the Eq. (28) can be solved in terms of v̇∥ by doing a scalar
multiplication with B∗. This results in the expression

v̇∥ = −µB∗

mB∗
||
· ∇B. (31)

The guiding-center equations of motion are thus given by Eq. (30) and (31).

15



3 Numerical methods
Since the purpose of this project is to simulate charged particle motion in a fusion
plasma, knowledge about various numerical methods is necessary. In this section we
will look at the theory and parts constituting the simulation program developed.

As mentioned earlier, equations of motion that include a complicated magnetic
field have to be solved numerically. The magnetic field from ITER, to be used in
the following simulations, certainly falls into this category. The magnetic field data
used in this project was given in matrix form with values separated by double spaces,
and had been calculated with a magnetic equilibrium solver that solves the Grad-
Shafranov equations [11], which describe a magnetohydrodynamic equilibrium in a
toroidally symmetric system. The magnetic field is accordingly toroidally symmetric.
A second data file was also provided, containing information about the wall shape of
the device, which is also toroidally symmetric. A brief introduction to the contents
of the data is given in Section 3.1.

The numerical integration of the equations of motion is done by an ODE solver.
The solver implemented here uses a Runge-Kutta method, which will be described
in Section 3.2. In Section 3.3, we proceed to present the method used to determine
whether the simulated particle collides with the device walls. Finally, in Section 3.4,
an overview of the simulation tool developed as part of this project is given.

All plots were made using the programming language Python and the library
matplotlib [12].

3.1 Magnetic field data
In order to make physically relevant simulations, magnetic field data from ITER
was used in the simulation program. The data file contained three matrices, the
field in cylindrical coordinates: BR(R, z), Bϕ(R, z) and Bz(R, z). A contour plot of
the projection of the magnetic field lines onto the R− z plane, with the device wall
superimposed over the field, is shown in Fig. 5a. The R−z plane is the poloidal plane.
Because of the toroidal symmetry of the device and the magnetic field, the azimuthal
coordinate has been neglected in Fig. 5, as the most interesting information comes
from looking at the field in the poloidal plane. Concentric, closed field lines fill up
most of the device. At the edges, mostly the left side and the bottom, the field
lines are open. Note especially the position of the magnetic field axis in 5a, which is
approximately at R = 6.7m, z = 0.5m. The entire device wall, plotted in 3D, can
be seen in Fig. 1. The wall data consisted of a number of (R, z)-points indicating
the wall contour.

16



3.2 Time integration method

3 4 5 6 7 8 9

R (m)

−6

−4

−2

0

2

4

6

z
(m

)

(a) Magnetic field lines

4 5 6 7 8

R (m)

−4

−2

0

2

4

z
(m

)

0

1

2

3

4

5

6

7

8

9

10

F
ie
ld

st
re
n
gt
h
(T

)

(b) Magnetic field strength

Figure 5: Two different views of the ITER magnetic field. In (a), the magnetic field lines
in the poloidal plane, with directional arrows, have been plotted, while in (b) the magnetic
field strength B = |B| is shown, also in the poloidal plane. In both plots, the ITER wall
contour has been superimposed to show the boundaries of the device. Note especially how
the field strength decreases with R.

Another important quantity of the magnetic field is the magnetic field strength
B = |B|. The field strength plays an important role for the guiding-center equations
of motion, and has been plotted in a color map in Fig. 5b.

When analyzing the system analytically, one can use the fact that the mag-
netic field of ITER behaves approximately like that for a toroid wound in current-
carrying wire.1 Using Ampere’s law, we find an approximation for the magnetic
field strength in a point with cylindrical coordinates (R, z) to be (see for example
Cheng (2014) [13])

B(R) ≈ B0R0

R
, (32)

where B0 is the magnetic field strength at the arbitrary point R0. This expression
predicts that the field strength will be inversely proportional to the radial distance
from the axis of symmetry, of the device, i.e. the middle of the toroid. By looking at
Fig. 5b we can rather easily convince ourselves of the validity of this approximation,
as the field seems to decrease monotonically with R.

3.2 Time integration method
The basic idea for any routine for solving an ordinary differential equation (ODE)
with initial value boundary conditions (called an initial value problem, IVP) is to
start at the initial values, and then take steps in the direction dictated by the func-

1In reality, there are additional features to the magnetic field, but these are generally small
compared to the toroidal magnetic field and can, for our purposes, be neglected.

17



3.2 Time integration method

tion corresponding to the derivative. This is the function f in Eq. (33). For small
steps, a good approximation to the underlying differential equation is achieved [14].

Consider an IVP over the time interval [t(0), t(f)],ż = f(t, z),

z
(
t(0)
)
= z(0).

(33)

To obtain a numerical approximation of the solution z, the interval [t(0), t(f)] is
divided into N equal subintervals. Mesh points t(j) = t(0) + jh, j = 1, . . . , N are
selected, where

h =
t(f) − t(0)

N
,

is the step size. One of the simplest numerical procedures to solve an IVP is the
Euler forward method, which is a first-order integration scheme. Its formula is

z(n+1) = z(n) + hf(t(n), z(n)).

The Euler forward method is not recommended for practical use due to its poor
accuracy compared to other methods with an equivalent step size, it is also not
very stable. However, information from several Euler-type steps can be combined
to obtain a higher order method. This is the basis for Runge-Kutta methods.

Runge-Kutta numerical methods are a family of one-step methods for solving first
order ODEs, using function values in multiple stages within one step. The Euler
forward method is equivalent to a first stage Runge-Kutta method. One of the most
powerful and most used Runge-Kutta methods is the four-stage (RK4), which is
of order four: O(h4). This is often used in conjunction with an adaptive step size
algorithm. A common variant is the Runge-Kutta-Fehlberg method (RKF45) that
has an error estimate of order five, O(h5). That is the one we will use. A more
in-depth look at this method will now follow. For further information we refer the
reader to Numerical Recipes in C [15].

An explicit Runge-Kutta (RK) method of the mth stage is given by

z(n+1) = z(n) + h
m∑
i=1

amiki,

with

k1 = f
(
t(n),z(n)

)
,

k2 = f
(
t(n) + c2h,z

(n) + ha21k1

)
,

...

km = f

t(n) + cmh, z
(n) + h

m−1∑
i=1

amiki

 ,

for coefficients aij, bk and cj. These are arranged in a so called Butcher tableau, as
shown in Table 1. The specific coefficients used for this method are the Fehlberg
parameters that can be seen in Table 2. The parameters are acquired by solving

18



3.2 Time integration method

Table 1: A generic Butcher Tableau.

0

c2 a21
c3 a31 a32
... ... ... . . .
cm as1 am2 · · · am,m−1 am,m

b1 b2 · · · bm−1 bm
b̂1 b̂2 · · · b̂m−1 b̂m

Table 2: Fehlberg parameters for the RKF45 method.

1/4 1/4 0 0 0 0 0

3/8 3/32 9/32 0 0 0 0

12/13 1932/2197 −7000/2197 7296/2197 0 0 0

1 439/216 −8 3680/513 −845/4104 0 0

1/2 −8/27 2 −3544/2565 1859/4104 −11/40 0

25/216 0 1408/2565 2197/4104 −1/5 0
16/135 0 6656/12825 28561/56430 −9/50 2/55

a system of algebraic equations that comes from comparing coefficients of Taylor
series expansions of ż = f(t, z) and the k-vectors.

Accuracy is improved with decreasing step size h, but this comes with the draw-
back of longer computation time. To optimize the step size, an adaptive step size
control algorithm is implemented. At each step, two different approximations of the
solution are made and compared, one of order four (denoted as z(n+1)) and one of
order five (denoted as ẑ(n+1)). The estimated error for each vector component

εj =
∣∣∣z(n+1)

j − ẑ
(n+1)
j

∣∣∣ ,
is compared to a desired accuracy ε0. If any εj > ε0, the step size is decreased, and
if any εj < ε0, the step size is increased. The reason for calculating separate errors
for each component is that their values can differ in orders of magnitude, as in our
application. So it is necessary to have separate measures for the error to ensure that
the adaptive step size algorithm functions properly.

Here a fourth-order Runge-Kutta method with five stages is used along with a
fifth-order method with six stages. For six stages there are six k-vectors:

k1 = f
(
t(n), z(n)

)
,

k2 = f
(
t(n) + c2h,z

(n) + ha21k1

)
,

...

k6 = f

t(n) + c6h, z
(n) + h

5∑
j=1

a6jkj

 .

19



3.3 Domain check

The embedded fourth-order formula is

z(n+1) = z(n) + h

5∑
i=1

biki,

and a better value is determined using the fifth-order method

ẑ(n+1) = ẑ(n) + h
6∑

i=1

b̂iki.

Using the coefficients from the Butcher tableau, the error estimate for each com-
ponent εj can be calculated as

εj =
∣∣∣z(n+1)

j − ẑ
(n+1)
j

∣∣∣ =
∣∣∣∣∣∣

6∑
i=1

(bi − b̂i)(ki)j

∣∣∣∣∣∣ .
The optimal step size change is given by

hopt =

βh
(
ε0
ε

)1/5
, εj ≥ ε0

βh
(
ε0
ε

)1/4
, εj < ε0

(34)

where β ≃ 1 is a “safety factor”, because our error estimate is not exact. Usually
β = 0.8 or β = 0.9. This is used to further ensure that a small enough step is taken.

The correctness of our implementation of this solver was tested by solving the
Lotka-Volterra equations [16], also known as the predator-prey equations. The so-
lution to these equations is well known.

3.3 Domain check
In order to produce realistic simulations of particle trajectories, it is important to
remember that there exists a constraint on the particle position; the particle has to
be inside the device at all times. Thus an algorithm to check whether the particle
has collided with the device wall, needs to be constructed. One of the simplest ways
to do this is by checking whether certain line segments intersect [17].

Assume that the path the particle takes to travel from one point u0 = (x0, y0) to
u1 = (x1, y1) can be represented by a parametrization

u(t) = u0 + t(u1 − u0) , t ∈ [0, 1].

Now, suppose that the wall contour is represented by the coordinates (x̃i, ỹi)
N
i=0,

where N is the total number of points. Assume that a parametrization can be done
of the line connecting the point ni: v0 = (x̃i, ỹi) and the closest neighbor (in a certain
direction) v1 = (x̃i+1, ỹi+1). The parametrization for this line is then of the form

v(s) = v0 + s(v1 − v0) , s ∈ [0, 1].

If the particle path is intersecting this part of the contour, there exist a point p ∈ Ω
such that u(t) = p = v(s) equivalent to u0 + t(u1 − u0) = v0 + s(v1 − v0). Here Ω
is the domain of the device. By using each coordinate representation in x’s and y’s,
the equation can be rewritten in matrix form as

20



3.3 Domain check

(
x1 − x0 x̃i − x̃i+1

y1 − y0 ỹi − ỹi+1

)
︸ ︷︷ ︸

A

(
t
s

)
︸︷︷︸

x

=

(
x0 − x̃i

y0 − ỹi

)
︸ ︷︷ ︸

b

. (35)

This matrix-equation has a unique solution if the determinant of A is nonzero. Thus

det(A) =
∣∣∣∣∣x1 − x0 x̃i − x̃i+1

y1 − y0 ỹi − ỹi+1

∣∣∣∣∣ = (x1 − x0)(ỹi − ỹi+1)− (x̃i − x̃i+1)(y1 − y0) ̸= 0.

The solution is then given by x = A−1b where A−1 is

A−1 =
1

det(A)

(
ỹi − ỹi+1 x̃i+1 − x̃i

y0 − y1 x1 − x0

)
.

By doing the matrix multiplication the resulting equations become(
t
s

)
=

1

det(A)

(
ỹi − ỹi+1 x̃i+1 − x̃i

y0 − y1 x1 − x0

)(
x1 − x0 x̃i − x̃i+1

y1 − y0 ỹi − ỹi+1

)
=

=
1

det(A)

(
(ỹi − ỹi+1)(x1 − x0) + (x̃i+1 − x̃i)(y1 − y0)
(y0 − y1)(x̃i − x̃i+1) + (x1 − x0)(ỹi − ỹi+1)

)
.

Whenever there exists a solution to (35) such that

{
t = (ỹi − ỹi+1)(x1 − x0) + (x̃i+1 − x̃i)(y1 − y0)
s = (y0 − y1)(x̃i − x̃i+1) + (x1 − x0)(ỹi − ỹi+1)

t, s ∈ [0, 1],

for 0 ≤ i ≤ N , then the particle path intersects the device contour, which indicates
that the particle has collided with the device.

21



3.4 Program workflow

3.4 Program workflow
In the previous two subsections, the algorithms we have written were presented.
Here, the workflow of the complete program will be described and illustrated by a
flow chart. First, two things need to be mentioned.

To be able to get magnetic field data for every single point in the device the
discrete data has to be interpolated. A two-dimensional interpolation library called
Interp2d [18] was used for this purpose. This library also contains functions for
differentiation, which were used to calculate derivatives of the magnetic field for the
guiding-center method.

Second, the equations of motion were derived for a Cartesian coordinate system.
This was done to avoid the more complicated expressions for derivatives of fields
in curvilinear coordinates. However, by virtue of the toroidal symmetry of the
tokamak, cylindrical coordinates are preferred. The magnetic field data was given
in cylindrical coordinates, so coordinate transformations had to be done. This step
is not shown explicitly in the flow chart.

The program works as follows:

• Read input data containing the simulation time tend, particle mass, charge,
initial position and velocity. The magnetic and domain data files are also
given as inputs, as well as a flag indicating whether you want to solve for
particle or guiding-center motion.

• Read and store this data, passing the field data to the interpolator.

• Check if the initial position is inside the device. If no, stop.

• If yes, choose which problem to solve and store the respective initial values.

• Start the time integrator, run until it reaches tend.

• Check if memory is allocated for solution data, if not, allocate.

• Take one step with the integrator by calculating all variables for the time t+h,
where h is the time step length.

• Check if the error is small enough. If not, calculate the next position and
velocity for a smaller time step.

• If yes, check if the new position is inside the allowed domain.

• If yes, repeat until tend is reached. Then write the solution data to file.

• If no, exit with a message saying the particle has hit the wall and write the
solution data acquired so far to file.

22



3.4 Program workflow

Start

Input

Read domain and magnetic field data
Initialize interpolator

Initial position
inside domain?Stop no Guiding-center?

yes

Convert particle initial values
to guiding-center initial values

yes

Save initial values
to solution vector

no

t > tend?

Print solution
to csv-file

yes

Memory
allocated?

no
Run RKF45

yes

Allocate memory for
solution vector

no

ε < ε0?

Decrease step size
no

New position
inside domain?

yes

no

yes

Stop

All of our code can be found at https://github.com/eerosdisciples. The
QR-code leads to this website. For more thorough instructions on how to run the
program, please consult Appendix D.

23

https://github.com/eerosdisciples


4 Simulation of particle orbits
As part of this project, a computer program for calculating particle orbits was devel-
oped, using the theory presented in the previous sections. In this section we will now
use that tool to study the orbits of different charged particles that may be found in
fusion plasmas. The purpose of this section is to understand how the properties of
different particles affect their motion within the tokamak, and also to explain the
reasons for the peculiarities of their orbits.

In order to make physically relevant simulations, we will need to know certain
properties of fusion plasmas. The simulation tool used requires particle parameters
such as mass, charge and initial velocity to be given. For this reason we should
find the typical values of these parameters in fusion plasmas. Since, by definition,
a plasma is ionized, it will contain free electrons (mass m ≈ 5.49 · 10−4 u, charge
q = −e) and ions. In ITER, a plasma consisting of the two hydrogen isotopes
deuterium (2H or D+) and tritium (3H or T+) will be used, leading us to believe
that a large portion of the ions within the plasma will have mass m ≈ 2 - 3 u and
charge q = e. In a typical fusion plasma of temperature 1 − 10 keV [19], these
particles will have energies in the region 1− 10 keV. However, it turns out that the
characteristics of the particle orbits for D+ and T+ become visible (compared to the
scale of the device) first at higher energies, and so most particles in this section will
be given energies in the MeV scale. The use of these high energies are not entirely
unphysical though, as such high energy ions can actually originate from neutral
beams used for plasma heating.

Other common particles within fusion plasmas are the helium ions called alpha
particles (mass m ≈ 4 u, charge q = 2e), which are produced in the nuclear reactions.
They are generated at the energy of 3.5MeV. Because they are fusion products they
will have much more kinetic energy than the average plasma temperature.

Apart from the particles just described, certain impurities will also be present. In
particular, we can expect to find carbon ions, C6+ (mass m ≈ 12 u, charge q = 6e),
and tungsten ions, W56+ (mass m ≈ 184 u, charge q = 56e), which originate from
the components facing the plasma and from the wall of the reactor. Since these
particles do not originate from nuclear reactions, they will have energies in the same
order as the plasma temperature, 1− 10 keV.

Except for the electron, all particles described so far have positive charges. In
order to illustrate an interesting property of the orbit, we will however need a particle
with a negative charge. The electron could in principle be used, but due to its small
mass the effect would require us to look very closely at the particle orbit in order
to see it. Because of this, we also introduce a particle called protide. This particle
is a hydrogen isotope, consisting of one proton and two electrons, thus having mass
m ≈ 1 u and charge q = −e. In a real fusion plasma, this particle will be extremely
rare and short-lived, if at all present. The use of the particle in this project is rather
motivated by its mass and charge, which happens to give interesting results.

4.1 Observed orbit topologies
Knowing what particles we may expect to find in a plasma, and that the magnetic
field will lock them to an orbit, a natural question would be what these orbits will
look like. From the theory derived in Section 2, we would expect the orbits to look a
bit different for different particles, depending on their mass, charge, initial position

24



4.1 Observed orbit topologies

and initial velocity. Using the simulation tool developed in this project, we have been
able to simulate both the particle and guiding-center orbits for different particles,
and in the present section our findings from these simulations will be presented. A
deeper analysis of these observations will be postponed until Section 4.2.

A simple test of the program in which we place an alpha particle at R = 8.0287m
and z = 0.2538m, and give it a kinetic energy of E = 3.5MeV, yields the orbits
shown in Fig. 6, depending on how its initial velocity vector is directed. The simula-
tion program calculates either the particle’s coordinates (the black orbits in Fig. 6)
or the corresponding guiding-center coordinates (the red orbits in Fig. 6) in three
dimensions (hence the three-dimensional plots in the upper part of the figure), but
as can be seen from the lower parts of the figure, the orbits may advantageously be
plotted in the two-dimensional R-z (poloidal) plane, which turns out to be much
easier to work in. The reason for this is simply that the magnetic field used is
symmetric in the azimuthal angular coordinate (assuming cylindrical or similar co-
ordinates are used), and so our three-dimensional geometry can be reduced to a
simpler two-dimensional one. It can be worth noting that on top of the banana
motion, the orbits have a precession around the tokamak. This information is lost
in the 2D figure. This may be important, because as the trapped electrons do their
precession around the tokamak, they can get in resonance with drift waves and can
generate instabilities (so called trapped electron modes). Nonetheless, this is not
encompassed by this project and we proceed to the analysis of only 2D figures.

Fig. 6 shows the two most common orbits found in tokamaks, labeled banana and
passing orbits (Fig. 6a/6c and 6b/6d, respectively). Particles following a banana or-
bit are commonly referred to as trapped particles, while particles following a passing
orbit are simply referred to as passing, or circulating, particles. These are the only
two orbit topologies we will encounter, and we will see that they are closely linked.
As mentioned above, the only difference between the orbits is how the particle’s
initial velocity is directed. A measure of this is the quantity ξ, that can be defined
as

ξ = v∥/v, (36)
which will become a powerful tool in the analysis of the transitions between banana
and passing orbits. We will only be interested in the value of ξ at the beginning of
the orbit, and so will denote it by ξ0 in order to distinguish it from the continuously
changing ξ.

As can be seen in Fig. 6 (with some effort), the particle orbits form bent helices
that follow the corresponding guiding-center orbits. A magnification of the particle’s
orbit is also shown in Fig. 7, where the helix shape is more prominent. The magnifi-
cation reveals an interesting feature of the particle’s orbit, which we shall return to
later, namely the fact that the guiding-center seems to slow down at the “banana
tips” (called mirror points), i.e. the particle’s velocity parallel to the magnetic field
lines decreases.

4.1.1 Orbit width
The two simplest parameters that can be varied in the program are the mass and
charge of the particle to be simulated. This section will be dedicated to the study
of how these two parameters affect the particle’s orbit. In order to have physically
relevant results, as far as possible, some the most common kinds of particles found

25



4.1 Observed orbit topologies

(a) (b)

4 5 6 7 8 9

R (m)

−4

−2

0

2

4

z
(m

)

(c)

4 5 6 7 8 9

R (m)

−6

−4

−2

0

2

4

6

z
(m

)

(d)

Figure 6: Example orbits for an alpha particle with energy 3.5MeV. All four plots show
the orbit that the particle follows (black), as well as the orbit that the guiding-center
follows (red), both of which are calculated using the tool developed as part of this project.
The orbits presented in these figures are classified as either banana orbits ((a) and (c)),
or passing orbits ((b) and (d)). As can be seen, 2D plots are usually easier to work with,
and also give a hint about where the orbits get their names from.

26



4.1 Observed orbit topologies

Figure 7: Magnification of the banana orbit for an alpha particle with kinetic energy
3.5MeV. Here, the helix shape of the particle’s orbit becomes apparent, motivating the
guiding-center transformation. As can be seen, the velocity parallel to the magnetic field
appears to decrease as the particle approaches the mirror points of the banana.

in fusion plasmas will be used in simulations. These are alpha particles, deuterium
(D+), tritium (T+), carbon (C6+) and tungsten (W56+). Note that since these are
part of the plasma, they are all ionized and thus possess positive charges.

In Fig. 8, four different particles (D+, T+, C6+, W56+) have been simulated
using the guiding-center method for the same energies E = 1MeV, initial position
(R = 8m, z = 0m) and same ξ0 = 0.5. The result is the four banana orbits shown,
which all have different widths. From this figure we can immediately draw one
important conclusion, namely that the banana width increases with an increased
mass. This can be seen by comparing the orbits for D+ (mass m ≈ 2 u) and T+

(mass m ≈ 3 u), of which the T+ orbit is somewhat wider. Knowing this, we can
then also conclude that the width decreases for an increase in charge magnitude,
which is easily seen by looking at the orbit for W56+. The banana width must
in other words be proportional to some ratio between mass and charge. A simple
calculation reveals that this ratio cannot simply be m/q, and so we should instead
expect it to be of the form mα/qβ, for α, β > 0. This will be analyzed in more detail
later in Section 4.2.

27



4.1 Observed orbit topologies

5.5 6.0 6.5 7.0 7.5

R (m)

2.0

2.2

2.4

2.6

2.8

3.0

3.2

3.4

z
(m

)

W 56+

C
6+

D+

T+

Figure 8: Four different guiding-center orbits for particles with energy 1MeV, but with
different mass/charge ratios. As can be seen, increased mass appears to result in a wider
banana orbit, while increased charge appears to result in the opposite.

Fig. 8 shows several interesting properties of banana orbits, but one thing that
cannot be seen is the effect of the charge sign. In order to study this we must turn to
the somewhat unrealistic choice of the protide particle, consisting of one proton and
two electrons, thus having mass m ≈ 1 u and charge q = −e.2 Passing this particle,
and a regular proton, to the guiding-center simulation program, giving them the
same energies E− = E+ = 1MeV and initial positions, yields the orbits shown in
Fig. 9 for ξ0 = 0.5 and ξ0 = 0.7, respectively. Arrows have been introduced in the
plots to show in which directions the particles’ guiding-centers follow the orbits.

As can be seen in Fig. 9a, depending on the sign of the particle’s charge, the
guiding-center banana orbit will turn either inwards (positive charge) or outwards
(negative charge) at the mirror points. For the passing orbits in Fig. 9b we instead
see how the sign of the charge causes the orbit to extend further towards the center
of the fusion device. This gives a hint for what causes the banana to have a width,
as it seems something is pushing the particle in different directions. We will see
later that this is due to particle drifts, as discussed in Section 2.1.2. On a side note,
if the sign of the initial value of the velocity is changed, you find orbits similar to
when the sign of the particle charge is changed.

2The reason we are not simply using an electron is that its small mass causes the width of
its banana orbit to make the effects we want to study virtually invisible, unless we zoom in very
closely. Zooming in would however instead cause us to miss the effect of the protons charge sign

28



4.1 Observed orbit topologies

4 5 6 7 8 9

R (m)

−6

−4

−2

0

2

4

6
z
(m

)

proton

protide

(a) Banana orbits

4 5 6 7 8 9

R (m)

−6

−4

−2

0

2

4

6

z
(m

)

proton

protide

(b) Passing orbits

Figure 9: Examples of banana and passing orbits of an alpha particle and a protide,
plotted in the same figures for reference. In both figures the particles’ initial positions
were (R, z) = (8.0287, 0.2538)m and their energies were E = 3.5MeV. Note how the
banana width of the protide is much thinner than that of the alpha particle. Also note
that directions are opposite in the banana orbit, due to the difference of the particle
charges.

Importance of energy for banana width
So far we have simulated several different particles with the same energies. It turns
out however, that also the energy of the particle affects the banana orbit width.
In Fig. 10, three different alpha particles (mass m ≈ 4 u, charge q = 2e) were
simulated with energies E1 = 3.5MeV, E2 = 1MeV and E3 = 350 keV. The results
strongly remind of those in Fig. 9, and so it seems that whatever causes the banana
width is not just affected by the particle’s mass and charge, but also by it’s speed
(since E = mv2/2). We see in Fig. 10 how the banana width increases by an
increase in energy, and so it would be reasonable to think that the banana width is
proportional to some positive power of the energy. This will later be shown to be
true in Section 4.2.

29



4.1 Observed orbit topologies

6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4

R (m)

2.0

2.2

2.4

2.6

2.8

3.0

3.2

3.4
z
(m

)
350 keV

1 MeV

3.5 MeV

Figure 10: The guiding-center orbits for three different alpha particles with kinetic
energies 3.5MeV (red), 1MeV (green) and 350 keV (blue). The value of ξ0 was ξ0 = 0.5.
As can be seen, the banana width decreases as the particle’s kinetic energy is decreased.

4.1.2 Deviation from field lines for banana orbits
Earlier in Section 2.1.2 we derived the different drifts that exist in a magnetized
plasma. In this section we will continue with the observation of particle orbits, by
studying how the different drifts affect the motion in terms of deviation from field
lines, for a banana orbit.

With known coordinates of a certain field line, we define the deviation as the
perpendicular distance from this field line to the trajectory of the particle. The
coordinates of a certain field line can be obtained by removing all but the parallel
motion in the direction of the magnetic field, however the particle is still experiencing
the same mirror force, and will hence have a reflection for the same azimuthal angle.

In the simulations, all particles had energies of 10 keV as expected in a thermal
plasma, apart from the alpha particle which had an energy of 3.5MeV. The devi-
ation was studied for one initial position of about (R, z) = (6.8, 0)m and a value
of ξ0 = −0.18. The result can be found in Figures 11 and 12. Due to the large
number of points that needed to be calculated, only the guiding-center method was
used in Fig. 12. Note that since each particle will have a different velocity, a plot of
the deviation as a function of time would result in displaced curves. Therefore, to
assist the comparison, we parametrize the particle paths along the field line by the

30



4.1 Observed orbit topologies

index i. This way, we can more easily compare the deviations of the orbits from the
magnetic field line.

At a first look we notice the rather unexpected local minimum which occurs
between the points of maximum deviation for each particle. We also see that the
minimum becomes less clear and finally disappears for lower values of the maximum
deviation. When the particles have been reflected at the mirror point, they move
backwards, but now on a different field line. Thus the result we see in Fig. 11 and 12
show that the field lines are not circularly symmetric around the center of the device.
The longer the jump to a different field line, the greater the change of the field line
curvature.

The next thing we recognize is that the deviation, and hence the drifts, are very
different for particles with different mass, charge and energy. For the T+ and D+

particles, we have the same charge, but a difference in mass of about one atomic
mass unit. For these, we see a slight change in the distance from the field line,
where the maximum deviation is 18% lower for the deuterium particle compared
with tritium. We can conclude that when the mass is increased, so is the deviation
from the field line.

However for the tungsten particle, having the highest mass, the deviation is quite
low. The reason for this must be a higher value of charge, and thus a higher value
of charge decreases the deviation.

Next, if we study the mass-to-charge ratio, we notice that while this value is
the same for deuterium, carbon and alpha particles, the particle trajectories differ
substantially. When looking at a fixed value of energy, we see that both the charge
and the mass are increased by a factor of six for the carbon particle when compared
to deuterium, while the value of the maximum deviation decreases by about 57%.
This means that the deviation value cannot be proportional to m/q, but either charge
or mass, or both, must be raised to some power, as concluded in Section 4.1.1. The
same observations hold for the tungsten particle, where the maximum deviation has
decreased by about 85.3%, when compared to tritium, while the mass-to-charge
ratio has increased by about 10%.

As a third observation, in Figures 11 and 12 we notice how the energy affects the
deviation from the field lines. For the alpha particle with an energy of 3.5MeV in
Fig. 11, the deviation is ten times greater than for the other particles with lower
energy in Fig. 12. Note the different scales in the figures. This suggests that the
particle’s deviation from the field lines increases with the energy.

In Fig. 11, both the guiding-center and the particle motion method has been used
to simulate the alpha particle. Here we see that guiding-center seems to follow the
center of the particle gyration almost perfectly, and the small differences that occur
are probably due to error from the numerical calculations of the deviation from the
field lines. This topic, comparing the guiding-center method and the particle motion,
will be studied in detail in Section 4.3.

31



4.1 Observed orbit topologies

0 20 40 60 80 100

i

0.00

0.05

0.10

0.15

0.20

0.25

0.30

D
e
v
ia

ti
o
n

(m
)

α-GC

α-PM

Figure 11: The deviation from the field lines for the alpha particle, during one banana
orbit, simulated with both methods. The energy used for the simulations was 3.5MeV.
The index i is used to parametrize the particle paths along the field line.

0 20 40 60 80 100

i

0.000

0.005

0.010

0.015

0.020

0.025

0.030

0.035

D
e
v
ia

ti
o
n

(m
)

C6+

T+

D+

W56+

e

Figure 12: The deviation from the field lines for D+, T+, C6+, W56+ and e− particles
during one banana orbit. The energies used for the simulations was 10 keV. The index i
is used to parametrize the particle paths along the field line.

32



4.1 Observed orbit topologies

4.1.3 Banana mirror points
Despite varying mass, charge and energy for particles, the banana mirror points
remained constant in all figures above. This is also clearly seen in Figures 8 and 10.
It follows that the type of orbit (banana or passing) is not affected by any of these
parameters, but rather by some other parameter. It was mentioned briefly in the
introduction to Section 4.1 that it is actually the direction of the particle’s initial
velocity that determines whether it will follow a banana or passing orbit. It was
also mentioned that the quantity ξ0 = v∥/v can be used as a measure for how
the particle’s initial velocity is directed. We will now study how variations of this
parameter affects both the particle and its guiding-center orbits.

In Fig. 13, alpha particles with three different values of ξ0 have been simulated
at constant energy E = 3.5MeV. The left part of the figure shows particle orbits,
while the right part shows guiding-center orbits. The particle orbits in Fig. 13a
seem to agree fairly well with the guiding-center orbits in Fig. 13b, possibly with
the exception of the orbits for ξ0 = 0.6. Here, the particle’s orbit comes closer to the
point we call the orbit’s transition point, i.e. the point separating banana orbits from
passing orbits. A suggested reason of why this deviation occurs is that the guiding-
center approximation might have a small error compared with the exact particle
motion. We remind that first order guiding-center theory is used here, for better
accuracy higher order theory could be used. For larger values of ξ0 this error seems
to increase, which affects the mirror point coordinates and results in the observed
deviation. This will be discussed in detail later in Section 4.3.

As can be seen in Fig. 13, a banana orbit will be successively stretched for in-
creasing values of ξ0, until its tips touch and the banana orbit turns into a passing
orbit. It seems that the banana mirror point depends on ξ0, and to understand
why, we should look at how the particle/guiding-center energy behaves. First, let
us split the total energy E into two terms, E = E∥ + E⊥, where E∥ = mv2∥/2 is the
“parallel energy” and E⊥ = mv2⊥/2 = µB the “perpendicular energy”. These two
terms tell us how much energy is stored in the motion along the magnetic field lines
and perpendicular to the magnetic field lines, respectively. We already know from
energy conservation that the particle’s total energy must remain constant through
the entire orbit, thus we would not get much interesting information from just look-
ing at E. When we split the energy into parallel and perpendicular terms however,
as has been done in Fig. 14a and Fig. 14b, we observe how these two terms actually
vary over the course of the orbit and constantly balance each other.

There are several points of interest in Fig. 14, but possibly the most significant
one is the fact that the parallel energy terms for the banana orbits drop to zero in
the mirror points. Since E = E∥ + E⊥, the perpendicular energy term E⊥ must
balance out the parallel term there, and thus accumulates all of the 3.5MeV energy
in the mirror points. We can also conclude from this information that, since the
particle mass m is constant, the parallel velocity v∥ must also drop to zero at the
mirror points. Likewise, the particle’s perpendicular speed v⊥ must increase to its
maximum at the mirror points. By looking at the orbits, and keeping in mind
that the orbit is approximately locked to the magnetic field lines, this allows us
to conclude that the particle (or rather its guiding-center) must turn back in the
direction it came. In other words, the parallel speed v∥ = v · b̂ must change sign.

For the orbits resulting from ξ0 = 0.7 on the other hand, the parallel energy
term never goes to zero. We see in Fig. 14 how, for both the particle and guiding-

33



4.1 Observed orbit topologies

4 5 6 7 8

R (m)

−4

−2

0

2

4

z
(m

)

(a) Particle motion method

4 5 6 7 8

R (m)

−4

−2

0

2

4

z
(m

)

(b) Guiding-center method

Figure 13: Three different particle orbits for an alpha particle simulated using both the
particle motion method (a) and the guiding-center method (b). For each orbit, the particle
was placed in (R, z) = (8, 0)m and given a kinetic energy of 3.5MeV. The values of ξ0 were
ξ0 = 0.5 (red), ξ0 = 0.6 (green) and ξ0 = 0.7 (blue). As can be seen, the orbits simulated
using the particle motion method comes closer to the transition point between the banana
and passing orbit. This is expected, and due to properties of the guiding-center method.

center orbits, the E⊥ term never attains the full energy of 3.5MeV. The particle
will therefore never have to stop, in contrast to the other orbits, but can continue
moving along the magnetic field lines, thus forming a passing orbit.

Another thing worth noting is the quite large temporal offset between the energy
components for the particle orbit resulting from ξ0 = 0.6 and the energy components
for the corresponding guiding-center orbit. As we will see later, this quite large offset
appears for the same reasons as the corresponding particle orbit in Fig. 13a was
elongated, compared to the guiding-center orbit resulting from ξ0 = 0.6 in Fig. 13b.
This offset in the energy however gives us a hint that just a small change in mirror
point position will have significant effects on the orbit time (i.e. the time it takes
for the particle to return to its initial position), when the mirror point is close to
the transition point.

34



4.1 Observed orbit topologies

Total Energy µB mv2||/2
......

0 10 20 30 40 50 60 70 80

Time (µs)

0

1

2

3

4
E
n
er
gy

(M
eV

)

ξ0 = 0.5

ξ0 = 0.6

ξ0 = 0.7

(a) Particle motion method

0 10 20 30 40 50 60 70 80

Time (µs)

0

1

2

3

4

E
n
er
gy

(M
eV

)

ξ0 = 0.5

ξ0 = 0.6

ξ0 = 0.7

(b) Guiding-center method

Figure 14: The total energy (dotted lines), parallel energy (E∥, dashed lines), and per-
pendicular energy (E⊥, solid lines) components for the orbits of Fig. 13. The value of ξ0
was varied between ξ0 = 0.5, ξ0 = 0.6 and ξ0 = 0.7. It is particularly interesting to note
how the E∥ term goes to zero for the two banana orbits (red and green), while for the
passing orbit E∥ > 0.

The successive transition into a passing orbit observed in Fig. 13 for increasing
ξ0 leads to the conclusion that there must be some value of ξ0, let’s call it ξc, that
is such that for all ξ0 < ξc a banana orbit is formed, while for all ξ0 > ξc a passing

35



4.1 Observed orbit topologies

orbit is formed. In Fig. 15 this value has been plotted as a function of the radial
distance to the center of the fusion device. When generating the figure, an electron
with energy E ≈ 1.25 keV was placed in z ≈ 0.5m, which is approximately the z-
coordinate of the magnetic axis as shown in Fig. 5a. Had the particle been perfectly
aligned with the magnetic axis, a value of ξc = 0 would have been expected for
all radii on the “inner” side of the magnetic axis. The guiding-center method was
used for generation, as the particle motion method would require prohibitively long
computation time.

For each radial position there are two values of ξc that determine the kind of
orbit topology, since if ξc is one transition value, so is −ξc. Because of this, the
dependence is symmetric around the ξc = 0 line. We also see that the border of
the region forms a parabola, suggesting |ξc| ∝ (R − RM)1/2, where RM is the R
coordinate of the magnetic axis.

6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0

R (m)

−0.6

−0.4

−0.2

0.0

0.2

0.4

0.6

ξ c BANANA

PASSING

PASSING

Figure 15: The ratio ξc = v∥/v plotted as a function of R with z aligned with the magnetic
axis. The particle used for simulations was an electron with energy E ≈ 1.25 keV. When
the particle’s ξ0 is within the enclosed region in the figure, a banana orbit will be formed.
Otherwise, it will follow a passing orbit.

36



4.2 Explanation of orbits

4.2 Explanation of orbits
In the previous section we observed some of the fundamental properties of the two
most common orbit topologies found in tokamaks. Most of the observations made
so far have yet to be explained; this is the purpose of this section. First we will
try to use the theory derived in Sections 2.1.2 and 2.2.2 to explain the observed
effects on the width of the banana orbits. Next, the banana width will be related to
the orbit’s deviation from magnetic field lines in order to explain the observations of
Sections 4.1.1 and 4.1.2. Finally, the banana orbit’s mirroring effect will be examined
and an explanation for the observations will be given.

4.2.1 Effect of mass, charge and energy
After having seen that the banana width varies with both mass, charge and energy
in Section 4.1.1, we would now like to understand what causes the banana width.
The answer comes from the drifts derived in Sections 2.1.2 and 2.2.2. Note that the
derivations made in those sections are merely two kinds of expressions for the same
drifts.

Since drifts are perpendicular to the magnetic field lines, they will force the
particle to move across field lines instead of just following them perfectly. The
guiding-center velocity, as derived in Section 2.2.2, was found to be

Ẋ = v∥b̂+
mb̂

qB∗
||
×

(
v2⊥
2B

∇B + v2∥∇∥b̂

)
. (37)

In this expression we have a parallel velocity term, v∥b̂, and a drift term. The charge
sign dependence on the banana as observed in Fig. 9 then easily gets its explanation,
as the sign of the charge determines the sign of the drift term. Depending on the
charge sign, the guiding-center drift will be either to the left or right from the field
line. Thus, for the positive charge sign of the proton, we saw how the guiding-center
drifted towards the left side of the field line. For the protide on the other hand, the
drift velocity would have been directed in the opposite direction, and thus it drifts
off to the right of the field line.

Though the effect of the charge sign is easily seen from Eq. (37), the effect of
mass, charge magnitude and energy is not as straight forward. In order to estimate
the drifts we are now going to approximate the distance that the particle deviates
from the field lines.

Since our fusion device is toroidally symmetric, Noether’s theorem states that
there is some corresponding conserved quantity. This happens to be the canonical
angular momentum pϕ, given by

pϕ(R, z) = mRvϕ + qRAϕ, (38)

where R is the distance from the axis of symmetry of the device to the particle, vϕ
the azimuthal component of the velocity and Aϕ the azimuthal component of the
magnetic vector potential. The poloidal magnetic flux between the magnetic axis
and a given magnetic surface (also called flux surface) is 2πΨ, where the magnitude
of Ψ is

Ψ(R, z) = −RAϕ(R, z).

37



4.2 Explanation of orbits

ρp
Ψ

Ψ∗

Figure 16: The definition of the poloidal Larmor radius, ρp, as the distance between
an orbit and a magnetic surface, in terms of points on a given magnetic surface, and the
particle’s guiding-center orbit.

The canonical angular momentum (38) can then be written as

pϕ(R, z) = mRvϕ − qΨ. (39)

By definition, Ψ is constant on a magnetic surface in an axially symmetric magnetic
field, so the particle must drift across magnetic surfaces if vϕ ≈ v∥ changes (which
it does, as we have seen in Section 4.1.3). In other words, the variation of Ψ and v∥
must balance each other out.

We can introduce a quantity Ψ∗ defined as

Ψ∗ = −pϕ
q
.

Combining this with Eq. (39) gives us

Ψ∗ = Ψ− mRvϕ
q

. (40)

This quantity Ψ∗, often referred to as the “drift surface”, has a dimension of magnetic
flux. For a given pϕ, the magnetic surface for which the poloidal magnetic flux Ψ is
equal to Ψ∗ is the location of a particle with vϕ ≈ v∥ = 0. That is, at the mirror
points (banana tips). This means that for a trapped particle the drift surface Ψ∗
is the magnetic surface that goes through the banana tips because there vϕ ≈ v∥
vanishes. For a passing orbit v∥ never vanishes, so the orbit never passes the drift
surface, as illustrated in Fig. 16. This orbit is either always outside or inside the
drift surface, depending on the sign of q.

As we established above in conjunction with Eq. (39), Ψ must also vary if v∥
varies. Here Ψ does not label a specific flux surface, instead Ψ(t) denotes the value
of Ψ along the particle trajectory.

We introduce a vector ρp that is locally normal to the flux surfaces, and its length
ρp is the distance between the particle location Ψ and the drift surface Ψ∗, as is also

38



4.2 Explanation of orbits

shown in Fig. 16. The magnitude of ρp is

ρp = ρp ·
∇Ψ

|∇Ψ|
, (41)

and by construction
ρp · ∇Ψ = Ψ−Ψ∗.

The vector ∇Ψ is naturally normal to a flux surface. Using Eq. (40) we find

ρp · ∇Ψ =
mRvϕ

q
,

which combined with Eq. (41) gives

ρp =
mRvϕ
q|∇Ψ|

.

We can use that the poloidal magnetic field can be expressed as Bp = ∇ϕ × ∇Ψ,
with |∇ϕ| = 1/R to find that |∇Ψ| = RBp. Thus

ρp =
mvϕ
qBp

≈
mv∥
qBp

.

This expression holds many similarities to that of the regular Larmor radius, ρ, as
defined in Section 2.1, but this expression contains the poloidal field strength Bp
and parallel velocity v∥ unlike ρ. For this reason, we call ρp the poloidal Larmor
radius. The poloidal Larmor radius is the typical distance a trapped particle moves
away from the flux surface because of drifts.

We can estimate the orbit width δ from the variation of ρp along the orbit due
to the variation of v∥, denoted by ∆v∥:

δ ≈
m∆v∥
qBp

.

With this expression for the orbit width, we are able to estimate the effect of
mass, charge magnitude and energy. Since ∆v∥ =

√
2∆E∥/m we can write the orbit

width as
δ =

√
m

q

√
2∆E∥

Bp
, (42)

where the ratio
√
m/q appears. It was something like this that was predicted in

Sections 4.1.1 and 4.1.2, and when numbers are put in it also appears to agree with
the simulations of Section 4.1.2 quite well. From Eq. (42) we calculate that the
guiding-center orbit widths for D+ and T+ should be fairly similar, and differing
only by a small amount. The widths for C6+ and especially W56+ on the other hand,
should be much smaller. The width of the C6+ orbit compared to the orbit of D+

should differ by about 59%, which is in good agreement with the simulation that
gave the value as 57%. For W56+ and T+, we estimate that the former’s orbit should
deviate by 86% less than the latter’s (simulated value was 85.3%).

The observations made in Fig. 10 are also fairly well explained by Eq. (42). As
we can see, the orbit width should be proportional to (∆E∥)

1/2. This seems likely,
since we can at least conclude that the width should be proportional to some power
of E between 0 and 1.

39



4.2 Explanation of orbits

4.2.2 The mirroring effect
In Section 4.1.3, we saw how the parallel and perpendicular components of the
guiding-center energy vary, which determines whether a banana or a passing orbit
is followed. We also saw that in a banana orbit, the parallel velocity v∥ goes to zero
and is forced to change sign at the mirror points, forcing the particle to turn back.
In this section we will elaborate on the reasons for the particle to follow a certain
orbit, and relate the observations made in 4.1.3 to known mathematical theory.

Let us first recall the expression for the particle’s energy,

E = E∥ + E⊥ =
mv2∥
2

+ µB. (43)

In Section 2.2.1, and especially Eq. (25), we learned that the magnetic moment
µ = mv2⊥/2B is an adiabatic invariant, which means that E⊥ should vary due to
changes in the magnetic field strength. Looking at the magnetic field strength B
in Fig. 5b, we see that as the particle moves in the direction of decreasing R, the
magnetic field strength B increases, thus increasing E⊥ along with it. Since the total
energy must remain constant, E∥ is therefore forced to decrease, which is consistent
with the observations of Section 4.1.3. When the energy components vary, so must
the speed components v∥ and v⊥. At a point with a sufficiently high magnetic field
strength, the perpendicular energy term would completely dominate and force E∥
to go to zero, bringing v∥ to zero with it. At this point, which becomes the banana
mirror point, it would not be possible for the particle to continue along its path into
regions of even greater B, and so it is forced to go in some other direction. Since
the orbit is approximately locked to a magnetic field line (only approximately since
drifts will cause it to move across field lines) there is only one other possible direction
in which the particle can go, i.e. the same way where it came from. Moving in that
direction, combined with the drift effects described in Section 4.2.1, causes the orbit
to look like a banana in the poloidal plane.

For certain values of ξ0, the adiabatic invariant µ will be sufficiently small for E⊥
never to dominate completely. In this case, v∥ never goes to zero and is therefore
never forced to change sign. The particle will just continue to follow a magnetic
field line and form a passing orbit. This was the case for the orbit with ξ0 = 0.7 in
Fig. 14, where we could see how the E⊥ term remained below 3.5MeV at all times
for this initial condition.

So how would we go about making the particle go further along the field line,
effectively making the banana longer? The key is in balancing the two energy terms
E∥ and E⊥ appropriately. Assuming the initial position of the particle to remain
unchanged, the magnetic field would obviously not change if we change any other
particle parameter. We can therefore make the parallel energy E∥ “last longer” by
increasing its initial value, thus allowing the particle to go further before losing all
its parallel velocity. This is done simply by increasing the particle’s initial velocity
in the direction parallel to the magnetic field. Note however that increasing E⊥ at
the same time, by an equal amount, would cause no change to the banana length.

Now, let E0 denote the total energy of the system (which will remain constant in
time). We may then rewrite Eq. (43) as

E∥ = E0 − µB.

40



4.2 Explanation of orbits

For a banana orbit, as we have seen, the E∥ term will vanish in the mirror points
of the banana. If we let Bmax be the maximum strength of the magnetic field along
the magnetic field line that the particle follows, we find that the particle will follow
a banana orbit only if

E0 − µBmax < 0.

Rewritten, this can be stated as requiring that µ > E0/Bmax, for the particle to
follow a banana orbit. But if we let Bmin be the minimum magnetic field strength
the particle experiences along the orbit, and note that µ can be expressed in terms
of v⊥ and B in this point, as µ = mv2⊥/2Bmin, we can use the fact that µ is constant
and that v∥ = 0 (implying v⊥ = v) in the mirror points to substitute E0 for mv2/2
and find the condition for the particle to follow a banana orbit∣∣∣∣v⊥v

∣∣∣∣ >√Bmin

Bmax
.

Now, using the fact that v2⊥ = v2−v2∥ and the definition of ξ in (36), we may express
a condition for the initial value of ξ,

|ξ0| <
√
1− Bmin

Bmax
, (44)

which will cause the particle to follow a banana orbit.

4.2.3 Mirror point dependence of ξ

As mentioned in Section 3.1, a simple approximation for the magnetic field strength
is

B(R) ≈ B0R0

R
, (45)

where B0 is the field strength at the arbitrary point R0. Using (45) we can rewrite (44)
to find an expression for the location of the mirror points of the banana orbit. We
already concluded that the magnetic field strength attains its maximum value, Bmax,
in the banana mirror points. If we call the radial position of the mirror point Rmax,
Eq. (44) becomes an equality in this point. it now expresses the value of ξ0 that gives
a mirror point in R = Rmax. Here we can use the magnetic field approximation (45)
to approximate both Bmin and Bmax, and substituting these into (44) and solving
for Rmax, we find

Rmax =
(
1− ξ20

)
Rmin. (46)

For a particle starting with its z coordinate aligned with the magnetic field axis,
Rmin becomes the initial R of the particle. In Fig. 17, the mirror points for an alpha
particle starting in (R, z) = (7.9, 0.5)m (almost aligned with the magnetic axis),
plotted as a function of ξ0 is shown. Simulation was done using the guiding-center
method. The figure also includes the theoretical curve predicted from Eq. (46) with
Rmin = 7.9m, and as can be seen the simulated curve agrees quite well. The small
deviation is most likely due to that the orbit is not initiated with z exactly aligned
with the magnetic axis.

This dependence is interesting as it gives us a very simple (yet apparently quite
good) method for predicting whether a particle will follow a banana or a passing
orbit, just from knowing its initial position (given z is aligned with the magnetic

41



4.2 Explanation of orbits

axis) and ξ0. In addition to what has already been said, it can be noted that for
certain values of ξ0 and Rmin, the predicted Rmax will never be reached, as it lies
too close to the center of the device. This is the case for all passing orbits, as can
be seen from Eq. (44).

−0.4 −0.2 0.0

ξ0

6.0

6.5

7.0

R
m
a
x

(m
)

α Theory

Figure 17: The mirror point as a function of initial ξ for an alpha particle compared
with the theoretical position given by Rmax = R0(1 − ξ2). The coordinates (R0, Z0) for
the initial position used in the simulation were (7.9, 0.5)m.

In Eq. (44) we could, instead of just studying any mirror point, study the transi-
tion point ξc, which sets the limit for a banana and a passing orbit. The simulation
results of its value for an electron, starting with z aligned with the magnetic axis,
was shown in Fig. 15. Here, we can start from Eq. (44) and note that for a ba-
nana orbit reaching just to the transition point, Rc, the maximum magnetic field
strength is Bmax ≈ B0R0/Rc. Substituting this approximation, along with a similar
approximation for Bmin, we find

|ξc| =
√
1− Bmin

Bmax
=

√
1− Rc

Rmin
. (47)

This gives the value of ξ0 that separates the banana orbits from the passing orbits.
The transition values for ξ0(R), as given by (47), were shown in Fig. 15 and the
estimated Rmin dependence for ξc ∝ (Rmin − RM)1/2 given there appears to be
approximately true. In Eq. (47), we also find support for the statement that ξc should
go to zero as Rmin approaches alignment with the magnetic axis in the z direction.

42



4.3 Comparison between simulation methods

Since magnetic field lines are similar to ellipses converging towards the same point,
being perfectly aligned with the magnetic axis would allow the particle to start from
Rmin = RM. As Rmin goes to RM, so would Rc. When Rmin = RM finally, we would
also have Rc = RM and so ξc would be zero.

Another interesting conclusion can be drawn from Eq. (47) regarding the mirror
point location, as predicted by Eq. (46). By definition we will have a passing orbit
for all values of ξ0 > ξc, but since the orbit is locked to a magnetic field line, it should
still pass the transition point (neglecting drifts). The “mirror point” (which is rather
a furthest point for the passing orbit) as predicted by Eq. (46) is no longer valid, and
must therefore only apply to banana orbits. For passing orbits, the “mirror point” is
instead approximately equal to the transition point Rc. Due to the increased effect
of particle drifts however, the orbit will in reality pass through a point R < Rc.

4.3 Comparison between simulation methods
Previously we have introduced the guiding-center formalism, separating the small
scale gyro-motion from the average motion of the particle, in which our main interest
lies. In this way we have been able to simulate the motion of the particle throughout
the device without needing to calculate the gyration at each point. However, up to
this point we have not yet discussed the quality of this transformation regarding the
physics and the resulting trajectory of the particle. In this section, we will compare
the two equations with each other in terms of energy conservation, orbit topology
and also computational time.

4.3.1 Energy Conservation
The first topic to be studied in this section will be the energy conservation as a
function of time for both methods. The aim is to observe the difference between the
methods concerning energy conservation for a fixed choice of initial parameters.

The energy conserving property of solutions for particle motion and guiding-center
motion is compared in Fig. 18. The trajectories were simulated over a time period
of 0.8ms, with a fixed error tolerance ε0 = 10−7. It is clear from the plot that for
the particle motion a numerical error accumulates, causing the energy to increase.
The energy change in the guiding-center approach is not appreciable during the
simulation time.

The simulations used to produce Fig. 18 bring other numerical advantages of
the guiding-center motion into light. The particle simulation took 63.3 s to execute,
producing 1 980 815 data points. Meanwhile, the guiding-center simulation execution
time was 1.1 s and produced 4 474 data points. The amount of steps taken to produce
the sizable number of data points for the particle motion is what causes the energy
level to deviate, since the RKF45 method is not inherently energy conserving and
thus causes the error to accumulate. The fact that the RKF45 method is the cause
of this inaccuracy makes it possible to regulate the conservation of energy depending
on method and initial parameters used. This will be utilized in the following section.

43



4.3 Comparison between simulation methods

0 100 200 300 400 500 600 700 800

Time (µs)

3.2

3.3

3.4

3.5

3.6

3.7

3.8

E
n

er
gy

(M
eV

)

(a) Energy for particle motion. The energy is not conserved, it accumulates.

0 100 200 300 400 500 600 700 800

Time (µs)

3.2

3.3

3.4

3.5

3.6

3.7

3.8

E
n

er
gy

(M
eV

)

(b) Energy for guiding-center motion. The energy is seemingly conserved.

Figure 18: Energy for alpha particles plotted over time for the numerical solutions of
the particle motion and guiding-center motion. All parameters were the same for both
simulations.

4.3.2 Computational time
In this section we will study the program execution time for the two methods, as we
increase the time for which we want to follow the particle. In this way, this section
will give another motivation of the guiding-center method, as well as motivating
the use of particles with higher energies, compared with the values expected in a
thermal plasma.

If a comparison of the computational time between the methods is to be done
correctly, it is important that both of the equations result in the same energy con-
servation. The amount of which the initial and final energy values deviate, will be
regarded as a calibration parameter, putting both of the algorithms on equal ground.

As presented in Section 3.2, the optimal step size hopt is calculated with re-
spect to the fixed error tolerance ε0, see Eq. (34). The smaller the step size, the
greater the accuracy in the result compared with the theoretical prediction. Hence,
a smaller value of ε0 will lead to a better conservation of energy. Since the num-
ber of steps taken differ greatly between the implementation of the guiding-center
equations and the equations of particle motion, it follows that the latter is far more
dependent on the step size and thus the error tolerance ε0. In the comparison that

44



4.3 Comparison between simulation methods

follows, we will fix the order of energy conservation to about six digits, so that
Einitial/Efinal = 1± 0.5 · 10−6.

For a start, let us define the simulation length tend as the time for which the
orbit should be traced. Likewise, we define the wall-clock time as the actual time it
takes to run the program for a certain simulation. The wall-clock time taken for a
particular method to calculate the needed data points, as a function of simulation
length, is given by Fig. 19. There we see that the wall-clock time for the guiding-
center method is much less, and increases slower, than for the corresponding particle
motion algorithm.

0 100 200 300 400 500

Simulation length (µs)

0

5

10

15

20

25

30

W
al
l-
cl
o
ck

ti
m
e
(s
)

GCM

Particle motion

Figure 19: Comparison between the particle motion method and the guiding-center
motion method of how much wall-clock time calculations require to complete. On the x
axis is the simulation length, i.e. tend, and on the y axis is the wall-clock time it took
for the program to carry out the calculations. A linear curve fit has been made to the
data which gives an inclination of kPM ≈ 53 300 for particle motion and kGCM ≈ 1 100 for
guiding-center motion. That means the guiding-center method implemented is almost 50
times faster than the particle motion method.

One thing we notice when running simulations of particles with different energies,
is that the wall-clock time increases with energy. The simulation time it takes for
a low energy particle to travel a certain distance is higher than the simulation time
it would take for a particle with higher energy to cover the same distance. As the
simulation time increases, so does the number of data points calculated for this
particular particle. Since there is a limit to how many data po