Initialization

In the following, we will discuss how to initialise the various models from a common set of initial data, namely

• the guiding center position $X$ in physical coordinates,
• the particle energy $W$ in eV,
• the particle mass $m$ in kg,
• the particle charge number,
• the gyro angle $\theta \in [0, 2\pi]$,
• the pitch angle $\alpha \in [0, \pi / 2]$.

For charged particles we need to convert the energy into a velocity vector. For the Pauli particle we need to compute the coordinates of the guiding centre $X$ of the particle and construct a velocity vector that only holds the component of the velocity that is parallel to the magnetic field, while the perpendicular component goes into the magnetic moment $\mu$. For the guiding center, we need the same coordinate transformation as for the Pauli particle, but we only need the absolute value of the velocity in direction of the magnetic field $u$ and the magnetic moment $\mu$.

Make sure to read about the Normalization before proceeding.

Velocities

Set $e = e' \hat{e}$ with $e' = 1.602 176 634 \cdot 10^{-19}$ and $\hat{e} = \mathrm{C}$ and similarly $m = m' \hat{m}$ with $\hat{m} = \mathrm{kg}$. Some common values for the mass are

The absolute value of the velocity is obtained from the particle energy by $\vert v \vert^2 = 2 W / m$. In normalised units, we have $W = W' \hat{W}$ and $v = v' \hat{v}$ where $\hat{W} = e \mathrm{V}$ and $\hat{v} = \hat{l} \hat{\omega}_c = \hat{l} \hat{B} e / m$.

Let us compute the normalise kinetic energy,

\[\frac{\vert v' \vert^2}{2} = \frac{W' \hat{W}}{m \hat{v}^2} = \frac{W' e \mathrm{V} m^2}{m \hat{l}^2 \hat{B}^2 e^2} .\]

Usually, we have $\hat{B} = T = \mathrm{kg} / \mathrm{C} \mathrm{s}$. Recall that $\mathrm{V} = \mathrm{kg} \mathrm{m}^2 / \mathrm{C} \mathrm{s}^2$ and let us set $\hat{l} = l_0 \mathrm{m}$, then

\[\frac{\vert v' \vert^2}{2} = \frac{m' W'}{e' l_0^2} .\]

Thus for the absolute value of the velocity we have

\[\vert v' \vert = \frac{1}{l_0} \sqrt{ 2 W' \, \frac{m'}{e'} } .\]

The pitch angle $\alpha$ determines the distribution of the kinetic energy into the perpendicular and parallel components by

\[\begin{aligned} W_{\perp}' &= W' \sin \alpha , \qquad W_{\parallel}' &= W' \, (1 - \sin \alpha) , \end{aligned}\]

and accordingly

\[\begin{aligned} \vert v_{\perp}' \vert &= \frac{1}{l_0} \sqrt{ 2 W_{\perp}' \, \frac{m'}{e'} } , \qquad \vert v_{\parallel}' \vert &= \frac{1}{l_0} \sqrt{ 2 W_{\parallel}' \, \frac{m'}{e'} } . \end{aligned}\]

The ElectromagneticFields.jl package provides three functions a⃗(t,x), b⃗(t,x), c⃗(t,x) that can be used to construct the velocity vector $v$. The function b⃗ returns the unit vector of the magnetic field, thus

\[v_{\parallel}' = \vert v' \vert \, (1 - \sin \alpha) \, b .\]

The functions a⃗ and c⃗ return two unit vectors that span the plane perpendicular to the magnetic field, thus

\[v_{\perp}' = \vert v' \vert \, \sin \alpha ( - a \, \sin \theta - c \, \cos \theta ) ,\]

where $\theta$ is the gyro angle. The magnetic moment is computed as

\[\mu = \frac{ \vert v_{\perp}' \vert^2 }{2 \vert B' \vert} = \frac{ \vert v' \vert^2 \sin \alpha }{2 \vert B' \vert} .\]

Gyro-radius Vector and Position

In order to compute the particle position, we need to construct the gyro radius vector $\rho$, which is given by

\[\rho = \frac{b \times v}{\omega_c} = \frac{b \times \hat{v} v'}{\hat{\omega_c} \omega_c'} = l_{0} \, \frac{b \times v'}{\vert B' \vert} .\]

Note that

\[\omega_c = \frac{e \vert B \vert}{m} = \frac{e \vert \hat{B} B' \vert}{m} = \frac{e \vert \hat{B} \vert}{m} \vert B' \vert = \hat{\omega}_c \omega_c' ,\]

as

\[\hat{\omega}_c = \frac{e \vert \hat{B} \vert}{m} ,\]

so that

\[\omega_c' = \vert B' \vert .\]

The normalized gyro radius vector reads

\[\rho' = \frac{\rho}{l_0} = \frac{b \times v'}{\vert B' \vert} ,\]

so that the normalised particle position is

\[x' = X' + \rho' .\]

The gyro phase $\theta$ is defined as the angle, measured in the clockwise sense, between $a$ and $\rho$, so that

\[\frac{b \times v'}{\vert v' \vert \, \sin \alpha} = a \, \cos \theta - c \, \sin \theta .\]

Guiding Center Coordinates

Instead of the particle energy $W$ and pitch angle $\alpha$, we can also compute all initial conditions from the guiding center initial data, i.e., parallel velocity $u = \vert v_{\parallel} \vert$, magnetic moment $\mu$ and gyro angle $\theta$.

Start by computing the absolute value of the perpendicular velocity by

\[\vert v_\perp \vert = \sqrt{ 2 \mu \vert B \vert } ,\]

and the total velocity by

\[\vert v \vert = \sqrt{ \vert v_{\parallel} \vert^2 + \vert v_\perp \vert^2 } .\]

With that, we can compute the perpendicular and total energy,

\[\begin{aligned} W_{\perp}' &= \frac{e l₀^2 \vert v_\perp \vert^2}{2m} , \qquad W' &= \frac{e l₀^2 \vert v \vert^2}{2m} . \end{aligned}\]

The pitch angle is obtained as

\[\alpha = \arcsin \frac{W_{\perp}'}{W'} .\]

With that, the rest of the quantities can be computed in the same way as above.

Initial Conditions

ChargedParticleDynamics.jl provides the InitialConditions module for the computation of the above quantities. Three functions are provided that return the initial conditions for the different models. For charged particles the initial conditions are $(x', v')$, for the Pauli particle we have $(X', v_{\parallel}', \mu)$ and for the guiding center we use $(X', u', \mu)$ with $u' = \vert v_{\parallel}' \vert$.

Example

As an example, let us consider a deuteron in an ITER-like analytical equilibrium (obtained from ElectromagneticFields.SolovevITER). The guiding center position is $[7, 0, 0]$, the energy is $1 \, \mathrm{MeV}$, and the pitch angle is $\pi / 2$.

Below we plot the guiding center position (in red) and the particle position (in blue) in the poloidal plane for various gyro angles $\in [0, 2\pi]$.