\documentclass[../main.tex]{subfiles}
\graphicspath{{\subfix{../images/}}}

\begin{document}
\section{Methods}\label{sec:methods}
\subsection{Theoretical background}\label{sec:theoretical_background}
In this experiment we study a particle with a charge $q$. To do so, we need to consider the forces acting on that particle. These forces are given by the Penning trap's electric field and magnetic field, an illustration can be found in fig. \ref{fig:penning_trap}. 
\begin{figure}
    \centering
    \includegraphics[width=0.6\linewidth]{images/penning_trap.pdf}
    \caption{The basic configuration of a hyperbolic Penning trap \cite{penning:fig}. Electric potential $V_{0}$ is applied to the electrodes in end caps (top and bottom), and in the ring (middle). The magnetic field is denoted as $\mathbf{B}$, the distance from the center to the end caps are denoted as $z_{0}$, and the radial distance as $r_{0}$.}
    \label{fig:penning_trap}
\end{figure}
The electric field $\mathbf{E}$ in our model is related to the electric potential $V$ through 
\begin{equation}\label{eq:electric_field_potential}
    \mathbf{E} = - \nabla V,
\end{equation} %
%
where we define the electric potential as
\begin{equation}\label{eq:electric_potential}
    V(x, y, z) = \frac{V_{0}}{2 d^{2}} (2z^{2} - x^{2} - y^{2}).
\end{equation} %
%
The characteristic dimension $d = \sqrt{z_{0}^{2} + r_{0}^{2} / 2}$ determine the scale of the region between the electrodes. The magnetic field is homogeneous and defined as
\begin{equation}\label{eq:magnetic_field}
    \mathbf{B} = B_{0} \hat{e}_{z} = (0, 0, B),
\end{equation} %
%
where $B_{0} > 0$ determines the strength of the field. The electric potential $V_{0}$ is applied to the electrodes, where the end caps are positively charged and the ring is negatively charged. The particle is confined by the electric field in the z-direction, however, it is not confined in the radial direction (xy-plane). The magnetic field is necessary to ensure the particle is fully confined in the Penning trap, and will force the particle to move in a circular orbit. 

First, we consider Newton's second law \eqref{eq:newton_second}, to determine the position of the particle. 
\begin{equation}\label{eq:newton_second}
    m \ddot{\mathbf{r}} = \sum_{i} \mathbf{F}_{i}
\end{equation} % 
%
In addition, we introduce the Lorentz force \eqref{eq:lorentz_force}, which describes the force acting on the particle.
\begin{equation}\label{eq:lorentz_force}
    \mathbf{F} = q \mathbf{E} + q \mathbf{v} \times \mathbf{B},
\end{equation} %
%
where $\ddot{\mathbf{r}} = d^{2} \mathbf{r} / dt^{2}$. We combine eq. \eqref{eq:newton_second} and eq. \eqref{eq:lorentz_force} in
\begin{equation}\label{eq:newton_lorentz}
    m \ddot{\mathbf{r}} = (q \mathbf{E} + q \mathbf{v} \times \mathbf{B}),
\end{equation} %
%
and derive the differential equations in appendix \ref{sec:eq_motion}. We can rewrite and simplify these equations as
\begin{align}
    \label{eq:motion_x}
    \ddot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x &= 0, \\
    \label{eq:motion_y}
    \ddot{y} - \omega_{0} \dot{x} - \frac{1}{2} \omega_{z}^{2} y &= 0, \\
    \label{eq:motion_z}
    \ddot{z} + \omega_{z}^{2} z &= 0,
\end{align} %
%
where
\begin{equation*}
    \omega_{0} \equiv \frac{q B_{0}}{m}, \quad \omega_{z}^{2} \equiv \frac{2 q V_{0}}{m d^{2}}.
\end{equation*} %
%
We find the general solution for eq. \eqref{eq:motion_z}
\begin{equation}\label{eq:eq_general}
    z(t) = c_{1} e^{i \omega_{z} t} + c_{2} e^{-i \omega_{z} t},
\end{equation} %
derived in appendix \ref{sec:eq_general}. Continuing, we will use a Calcium ion with a single positive charge. That is, we assume the charge of the particle is $q > 0$.

Since eq. \eqref{eq:motion_x} and eq. \eqref{eq:motion_y} are coupled, we want to rewrite them as a single differential equation. We derive this in appendix \ref{sec:eq_complex}, and the resulting equation is given by
\begin{equation}\label{eq:single_differential}
    \ddot{f} + i \omega_{0} \dot{f} - \frac{1}{2} \omega_{z}^{2} f = 0. 
\end{equation} %
%
Eq. \eqref{eq:single_differential} has a general solution given by
\begin{equation}\label{eq:general_solution}
    f(t) = A_{+}e^{-i(\omega_{+} t + \phi_{+})} + A_{-}e^{-i(\omega_{-} t + \phi_{-})}.
\end{equation} %
%
The amplitude $A_{+}$ and $A_{-}$ are positive, the phases $\phi_{+}$ and $\phi_{-}$ are constant, and the angular rate is given by
\begin{equation}\label{eq:angular_rate}
    \omega_{\pm} = \frac{\omega_{0} \pm \sqrt{\omega_{0}^{2} - 2 \omega_{z}^{2}}}{2}.
\end{equation} 
We find the physical coordinates at a given time $t$ using 
\begin{equation*}
    x(t) = \text{Re}f(t), \quad y(t) = \text{Im}f(t),
\end{equation*} %
and eq. \eqref{eq:general_solution}. We can rearrange the right hand side of the expression derived in \ref{sec:eq_coord}, to find the physical coordinates
\begin{align}\label{eq:physical_coord}
    x(t) &= A_{+} \cos(\omega_{+} t + \phi_{+}) + A_{-} \cos(\omega_{-} t + \phi_{-}) \\
    y(t) &= - A_{+} \sin(\omega_{+} t + \phi_{+}) - A_{-} \sin(\omega_{-} t + \phi_{-})
\end{align} %
% Problem 4
However, to obtain a bound on the particle's movement in the radial direction, we have to ensure that 
\begin{align*}
    |f(t)| &= \sqrt{(x(t))^{2} + (y(t))^{2}}. \\
\end{align*} %
Which means we have to put constraint on the values of $\omega_{0}$ and $\omega_{z}^{2}$ in order to find a bounded solution $|f(t)| < \infty$. We derive an expression for this in appendix \ref{sec:bounded_solution}. Since we have assumed the particle charge $q > 0$, and the mass $m > 0$, the constraints are
\begin{equation}
    \frac{q}{m} > \frac{4 V_{0}}{d^{2} B_{0}^{2}}.
\end{equation} %
% 
When $t \rightarrow \infty$, the upper and lower limits are
\begin{align}
    \label{eq:upper_b}
    R_{+} &= A_{+} + A_{-}, \\
    \label{eq:lower_b}
    R_{-} &= |A_{+} - A_{-}|, 
\end{align}
which we derive in appendix \ref{sec:upper_lower_bound}.

In this experiment, we will use the initial conditions for the particle given in table \ref{tab:initial_particle_cond}. From eq. \eqref{eq:motion_z} and the initial conditions, we find the specific solution
\begin{equation*}
    z(t) = z_{0} \cos (\omega_{z} t)
\end{equation*}
For the particle movement in the radial direction, we find the specific solution using eq. \eqref{eq:general_solution} with 
\begin{align*}
    A_{+} &= \frac{v_{0} + \omega_{-} x{0}}{\omega_{-} - \omega_{+}}, & A_{-} &= - \frac{v_{0} + \omega_{+} x{0}}{\omega_{-} - \omega_{+}}, \\
    \phi_{+} &= 0, & \phi_{-} &= 0.
\end{align*}

We will also consider the numerical simulation of multiple particles, confined in a Penning trap. The particles will experience a repelling force from each other, giving the electric field 
\begin{equation}\label{eq:electric_field_interaction}
    \mathbf{E} = k_{e} \sum_{j = 1}^{n} q_j \frac{\mathbf{r} - \mathbf{r}_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}}.
\end{equation} %
The Coulomb constant $k_{e}$ value can be found in table \ref{tab:constants}. 

For multiple particles we have to modify the equations of motion, by adding a term for the force a given particle experience at a given point. When we scale eq. \eqref{eq:electric_field_interaction} by charge and mass, we get a new set of equations of motion
\begin{align}
    \label{eq:coulomb_motion_x}
    \ddot{x}_{i} - \omega_{0,i} \dot{y}_{i} - \frac{1}{2} \omega_{z,i}^{2} x_{i} - k_{e} \frac{q_{i}}{m_{i}}\sum_{j \neq i} q_j \frac{x_{i} - x_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}} &= 0, \\
    \label{eq:coulomb_motion_y}
    \ddot{y}_{i} - \omega_{0,i} \dot{x}_{i} - \frac{1}{2} \omega_{z,i}^{2} y_{i} - k_{e} \frac{q_{i}}{m_{i}}\sum_{j \neq i} q_j \frac{y_{i} - y_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}} &= 0, \\
    \label{eq:coulomb_motion_z}
    \ddot{z}_{i} + \omega_{z,i}^{2} z_{i} - k_{e} \frac{q_{i}}{m_{i}}\sum_{j \neq i} q_j \frac{z_{i} - z_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}} &= 0,
\end{align} %
where $i$ and $j$ denote the particle indices. When we include a time-dependence to the applied potential we make a replacement of the initial electric potential
\begin{equation}\label{eq:pertubation}
    V_{0} \rightarrow V_{0} (1 + f \cos (\omega_{V} t)),
\end{equation}
where $f$ denotes the amplitude and $\omega_{V}$ the angular rate.

\subsection{Algorithms and implementation}\label{sec:algo_implementation}
When we consider a multi-particle system, we adapt the notation in our derived equations. The force acting on a particle in the Penning trap, which experience particle interaction, can be written as
\begin{equation}\label{eq:lorentz_force_interaction}
    \mathbf{F}_{i}(t, \mathbf{v}_{i}, \mathbf{r}_{i}) = q_{i} \mathbf{E}(t, \mathbf{r}_{i}) + q_{i} \mathbf{v}_{i} \cross \mathbf{B} - \mathbf{E}_{p}(t, \mathbf{r}_{i}),
\end{equation}
where $i$ and $j$ still denotes the particle indices, and $\mathbf{E}_{p}$ denotes the force exerted on particle $i$ by particle $j$.
\begin{equation}\label{eq:electric_force_interaction}
    \mathbf{E}_{p}(t, \mathbf{r}_{i}) = k_{e} q_{i} \sum_{j \neq i} 
    q_{j} \frac{\mathbf{r}_{i} - \mathbf{r}_{j}}{| \mathbf{r}_{i} - \mathbf{r}_{j} |^{3}}.
\end{equation} %
In addition, we adapt the notation of \eqref{eq:newton_second}, and define the first derivative of the particle position as
\begin{align}
    \label{eq:r_derivative}
    \frac{d \mathbf{r}_{i}}{dt} = \dot{\mathbf{r}}_{i} &= \mathbf{v}_{i}, \\
    \label{eq:v_derivative}
    \frac{d \mathbf{v}_{i}}{dt} = \dot{\mathbf{v}}_{i} &= \frac{\mathbf{F}_{i} (t, \mathbf{v}_{i}, \mathbf{r}_{i}}{m_{i}},
\end{align}

We first implemented the forward Euler method, using the expression for a coupled system given in appendix \ref{sec:algo_euler}. We define the forward Euler algorithm in \ref{algo:forward_euler}. 
\begin{figure}[H]
    \begin{algorithm}[H]
    \caption{Forward Euler method}
    \label{algo:forward_euler}
        \begin{algorithmic}
            \Procedure{Evolve forward Euler}{$particles, dt$}
            \State $N \leftarrow \text{Number of particles in } particles$
            \State $a \leftarrow \text{Calculate } \frac{\vb{F_i}}{m_i} \text{ for each particle in } particles$
            \For{ $i = 1, 2, \ldots , N$ }
            \State $particles_i.\vb{r} \leftarrow particles_i.\vb{r} + dt \cdot particles_i.\vb{v}$
            \State $particles_i.\vb{v} \leftarrow particles_i.\vb{v} + dt \cdot a_i$
            \EndFor
            \EndProcedure
        \end{algorithmic}
    \end{algorithm}
\end{figure}

We also implemented the 4th order Runge-Kutta (RK4) method, using the expression given in appendix \ref{sec:algo_rk4}. We define the RK4 algorithm in \ref{algo:rk4}. $\mathbf{F}$ does not take any arguments, however, the total force acting on the particle is calculated using the value of position and velocity within $particles$. 
\begin{figure}
    \begin{algorithm}[H]
    \caption{RK4 method}
    \label{algo:rk4}
        \begin{algorithmic}
            \Procedure{Evolve RK4}{$particles, dt$}
            \State $N \leftarrow \text{Number of particles inside the Penning trap}$
            \State $orig\_p \leftarrow \text{Copy of particles}$ 
            \State $tmp\_p \leftarrow \text{Array of particles of size }N$
            \State $\vb{k}_{\vb{r}} \leftarrow \text{2D array of vectors of size } 4 \cross N$
            \State $\vb{k}_{\vb{v}} \leftarrow \text{2D array of vectors of size } 4 \cross N$

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},1,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},1,i} \leftarrow \frac{\vb{F}_i}{m_i}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{r},1,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{v},1,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},2,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},2,i} \leftarrow \frac{\vb{F}_i}{m_i}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{r},2,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{v},2,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},3,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},3,i} \leftarrow \frac{\vb{F}_i}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + dt \cdot \vb{k}_{\vb{r},3,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + dt \cdot \vb{k}_{\vb{v},3,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},4,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},4,i} \leftarrow \frac{\vb{F}}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + \frac{dt}{6} 
                \cdot \left( \vb{k}_{\vb{r},1,i} + 2 \vb{k}_{\vb{r},2,i} 
                + 2 \vb{k}_{\vb{r},3,i} + \vb{k}_{\vb{r},4,i} \right)$

            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + \frac{dt}{6} 
                \cdot \left( \vb{k}_{\vb{v},1,i} + 2 \vb{k}_{\vb{v},2,i} 
                + 2 \vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} \right)$
            
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Final update}
            \EndProcedure
        \end{algorithmic}
    \end{algorithm}
\end{figure}
For a system of particles, the positional value of one particle affects the rest of the particles. We calculated the next step for all particles, updated their values in a copy, and updated the positional value based on the original values. We simulated the movement of particles confined in a Penning trap. All simulations used the initial conditions for particle 1 and 2 given in table \ref{tab:initial_particle_cond}. 
 
\subsection{Testing and error analysis}
We implemented a test suite to validate our results, details can be found in appendix \ref{sec:testing}. To analyze our results, we compared the relative error of the analytical solution with the result from the implemented forward Euler method and RK4 method. We simulated four times using time steps given in table \ref{tab:time_steps}, and estimated the error convergence rate $r_{\text{err}}$ for both methods using 
\begin{equation}
    r_{\text{err}} = \frac{1}{3} \sum_{k=2}^{4} \frac{\log(\Delta_{\text{max},k} / \Delta_{\text{max}, k-1})}{\log(h_{k} / h_{k-1})}. 
\end{equation}
The maximum error of a simulation $k$ with a step size $h_{k}$ is given by
\begin{equation*}
    \Delta_{\text{max},k} = \max_{i} |\mathbf{r}_{i, \text{exact}} - \mathbf{r}_{i}|,
\end{equation*}
where $\mathbf{r}_{i, \text{exact}}$ is the analytical solution, and $\mathbf{r}_{i}$ is the computed result.

\subsection{Tools}
The numerical methods are implemented in C++, and parallelized using \verb|OpenMP| \cite{openmp:2018}. In addition, we used the profiler \verb|Scalasca| \cite{scalasca} and \verb|Score-P| \cite{scorep} event tracking, during development. We used the Python library \verb|matplotlib| \cite{hunter:2007:matplotlib} to produce all the plots, and \verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the figures. 
\end{document}
