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

\begin{document}
\section{Derivation of equations}\label{sec:derivations}
% Problem 1
\subsection{Equations of motion}\label{sec:eq_motion} %
First, we need to define the velocity of the particle
\begin{align*}
    \mathbf{v} \equiv \frac{d \mathbf{r}}{dt} &= \bigg( \frac{dx}{dt}, \frac{dy}{dt}, \frac{dz}{dt} \bigg).
\end{align*} %
%
We can rewrite the velocity as $\dot{r} = (\dot{x}, \dot{y}, \dot{z})$, and find the cross product
\begin{align*}
    q \mathbf{v} \cross \mathbf{B} &= 
    q \begin{vmatrix}
        \hat{e}_{x} & \hat{e}_{y} & \hat{e}_{z} \\
        \dot{x} & \dot{y} & \dot{z} \\
        0 & 0 & B_{0}
    \end{vmatrix}
    = q \big( B_{0} \dot{y}, -B_{0} \dot{x}, 0 \big).
\end{align*} %
%
We are considering an ideal Penning traps, where we define the electric potential as
\begin{align*}
    V(x, y, z) &= \frac{V_{0}}{2 d^{2}}(2z^{2} - x^{2} - y^{2}).
\end{align*} %
%
The relationship between the electric field $\mathbf{E}$ and the electric potential of the field is given by
\begin{align*}
    \mathbf{E} &= - \nabla V \\
    &= - \bigg( \frac{dV}{dx}, \frac{dV}{dy} \frac{dV}{dz} \bigg) \\
    &= \frac{V_{0}}{d^{2}} \big( x, y, -2z \big).
\end{align*} %
%
We can now express the Lorentz force as
\begin{align*}
    \mathbf{F} &= q \mathbf{E} + q \mathbf{v} \cross \mathbf{B} \\
    &= \frac{q V_{0}}{d^{2}} \big( x, y, -2z \big) + \big(q B_{0} \dot{y}, -q B_{0} \dot{x}, 0 \big),
\end{align*} %
%
and insert it into Newtons equation \eqref{eq:newton_second}. We get
\begin{align*}
    \ddot{\mathbf{r}} &= \bigg( \frac{q V_{0}}{m d^{2}} x, \frac{q V_{0}}{m d^{2}} y, -\frac{2 q V_{0}}{m d^{2}} z \bigg) + \bigg(\frac{q B_{0}}{m} \dot{y}, -\frac{q B_{0}}{m} \dot{x}, 0 \bigg),
\end{align*} %
%
which can be written as
\begin{align*}
    \ddot{x} &= \frac{q V_{0}}{m d^{2}} x + \frac{q B_{0}}{m} \dot{y}, \\
    \ddot{y} &= \frac{q V_{0}}{m d^{2}} y - \frac{q B_{0}}{m} \dot{x}, \\
    \ddot{z} &= -\frac{2 q V_{0}}{m d^{2}} z.
\end{align*} %
%
If we define
\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*} %
%
the equations of motion can be written as
\begin{align*}
    \ddot{x} &= \frac{1}{2} \omega_{z}^{2} x + \omega_{0} \dot{y}, \\
    \ddot{y} &= \frac{1}{2} \omega_{z}^{2} y - \omega_{0} \dot{x}, \\
    \ddot{z} &= -\omega_{z}^{2} z. \\
\end{align*} %
%
\subsection{General solution}\label{sec:eq_general}
We consider the characteristic equation of a second order differential equation \cite{lindstrom:2016:ch10:5},  
\begin{align*}
    r^{2} + \omega_{z}^{2} &= 0 \\
    r &= \pm \sqrt{- \omega_{z}^{2}}.
\end{align*} %
The characteristic equation has two complex roots
\begin{equation*}
     r_{1} = - i \omega_{z}, \quad r_{2} = i \omega_{z},
\end{equation*}
which give us solutions in the general form
\begin{equation*}
    z = c_{1} e^{i \omega_{z} t} + c_{2} e^{-i \omega_{z} t}.
\end{equation*}
In addition, for a complex number $z = a + ib$, we can define $e^{z} \equiv e^{a} (\cos{b} + i \sin{b})$ \cite{lindstrom:2016:ch3}. We can rewrite the general solution as
\begin{align*}
    c_{1} e^{i \omega_{z} t} + c_{2} e^{-i \omega_{z} t} &= c_{1} (\cos{\omega_{z} t} + i \sin{\omega_{z} t}) \\
    & \quad + c_{2} (\cos{\omega_{z} t} - i \sin{\omega_{z} t} \\
    &= E \cos{\omega_{z} t} + i F \sin{\omega_{z} t}.
\end{align*} %
%
\subsection{Complex function}\label{sec:eq_complex} %
In sec. \ref{sec:eq_motion} we found the differential equations for $\ddot{x}$ and $\ddot{y}$. To derive a single differential equation, we introduce the complex function $f(t) = x(t) + iy(t)$, which gives us
\begin{align*}
    0 &= \Big(\ddot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x \Big) + i \Big(\ddot{y} + \omega_{0} \dot{x} - \frac{1}{2} \omega_{z}^{2} y \Big) \\
    &= \ddot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x + i\ddot{y} + i\omega_{0} \dot{x} - i \frac{1}{2} \omega_{z}^{2} y \\
    &= \ddot{x} + i\ddot{y} + i\omega_{0} \dot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x - i \frac{1}{2} \omega_{z}^{2} y. \\ 
\end{align*}
Using the definition $i = \sqrt{-1}$, we can rewrite 
\begin{equation*}
    i \omega_{0} \dot{x} + (-1) \omega_{0} \dot{y} = i \omega_{0} \dot{x} + i^{2} \omega_{0} \dot{y}.
\end{equation*} 
This gives us a single differential equation
\begin{align*}
    0 &= \ddot{x} + i\ddot{y} + i\omega_{0} (\dot{x} + i \dot{y}) - \frac{1}{2} \omega_{z}^{2} x - i \frac{1}{2} \omega_{z}^{2} y \\
    &= \ddot{f} + i \omega_{0} \dot{f} - \frac{1}{2} \omega_{z}^{2} f.
\end{align*}
%
\subsection{Physical coordinates}\label{sec:eq_coord}
We can rewrite eq. \eqref{eq:general_solution}, as
\begin{align*}
    f(t) &= A_{+}e^{-i(\omega_{+} t + \phi_{+})} + A_{-}e^{-i(\omega_{-} t + \phi_{-})} \\
    &= A_{+}(\cos{(\omega_{+} t + \phi_{+})} - i \sin{(\omega_{+} t + \phi_{+})}) \\
    % \numberthis \label{eq:general_solution_trig}
    & \quad + A_{-}(\cos{(\omega_{-} t + \phi_{-})} - i \sin{(\omega_{-} t + \phi_{-})}).
\end{align*} %
%
\subsection{Upper and lower bounds}\label{sec:upper_lower_bound}
To obtain the upper and lower bounds of the particle's distance from the origin, we first find an expression for the second norm defined as $|f(t)| = \sqrt{(x(t))^{2} + (y(t))^{2}}$. 
\begin{align*}
    (x(t))^{2} &= \big( A_{+}\cos(\omega_{+} t + \phi_{+}) + A_{-}\cos(\omega_{-} t + \phi_{-}) \big)^{2} \\
    &= A_{+}^{2} \cos^{2}(\omega_{+} t + \phi_{+}) \\
    & \quad + 2 A_{+}A_{-} \cos(\omega_{+} t + \phi_{+})\cos(\omega_{-} t + \phi_{-}) \\
    & \quad + A_{-}^{2}\cos^{2}(\omega_{-} t + \phi_{-}), \\
\end{align*} %
\begin{align*}
    (y(t))^{2} &= \big( - A_{+} \sin(\omega_{+} t + \phi_{+}) - A_{-} \sin(\omega_{-} t + \phi_{-}) \big)^{2} \\
    &= A_{+}^{2} \sin^{2}(\omega_{+} t + \phi_{+}) \\
    & \quad + 2 A_{+}A_{-} \sin(\omega_{+} t + \phi_{+})\sin(\omega_{-} t + \phi_{-}) \\
    & \quad + A_{-}^{2} \sin^{2}(\omega_{-} t + \phi_{-}).
\end{align*} %
We insert these expressions, and find 
\begin{align*}
    |f(t)| &= \sqrt{(x(t))^{2} + (y(t))^{2}} \\
    &= \sqrt{A_{+}^{2} + 2 A_{+} A_{-} \cos^{2}(\alpha) + A_{-}^{2}},
\end{align*} %
where $\alpha = (\omega_{+} - \omega_{-}) t +( \phi_{+} -  \phi_{-})$. If we set $\alpha = 0$ we get $\cos(0) = 1$, and obtain the upper bound
\begin{align*}
    R_{+} &= \sqrt{A_{+}^{2} + 2 A_{+} A_{-} + A_{-}^{2}} \\
    &= \sqrt{(A_{+} + A_{-})^{2}} \\
    &= A_{+} + A_{-}.
\end{align*}
If $\alpha = \pi$ we get $\cos(\pi) = -1$, and find the lower bound
\begin{align*}
    R_{-} &= \sqrt{A_{+}^{2} - 2 A_{+} A_{-} + A_{-}^{2}} \\
    &= \sqrt{(A_{+} - A_{-})^{2}} \\
    &= |A_{+} - A_{-}|.
\end{align*} %
%
\subsection{Bounded solution}\label{sec:bounded_solution}
To find a bounded solution, we need to consider the angular rate in eq. \eqref{eq:angular_rate}. Specifically the square-root expression
\begin{align*}
    \sqrt{\omega_{0}^{2} - 2 \omega_{z}^{2}}.
\end{align*}
A bounded solution can be found when this expression is greater than zero. Meaning
\begin{align*}
    \omega_{0}^{2} &> 2 \omega_{z}^{2}.
\end{align*} %
%
We can now use the definition of $\omega_{0}$ and $\omega_{z}^{2}$ to find an expression related to the Penning trap parameters, and the particle properties, to get
\begin{align*}
    \Big( \frac{q B_{0}}{m} \Big)^{2} &> 2 \frac{2 q V_{0}}{m d^{2}} \\
    \frac{q^{2} B_{0}^{2}}{m^{2}} &> \frac{4 q V_{0}}{m d^{2}}.
\end{align*} %

\section{Values used in simulation}
\subsection{Specific analytical solution}\label{sec:spec_analytical}
To compare our implementation of the Penning trap and particle, we have to specify initial conditions of the system to compare with the analytical solution. The initial conditions of a particle with a single positive charge, can be found in table \ref{tab:initial_particle_cond}. 
\begin{table}[H]
    \centering
    \begin{tabular}[c]{ll}
        $\mathbf{r}(0)$ & $\dot{\mathbf{r}}(0)$ \\
        \hline
        $x(0)$ = $x_{0}$ & $\dot{x}(0)$ = $0$ \\
        $y(0)$ = $0$ & $\dot{y}(0)$ = $v_{0}$ \\
        $z(0)$ = $z_{0}$ & $\dot{z}(0)$ = $0$ \\
        \hline
    \end{tabular}
    \caption{Initial values of a particle with a single charge $q$, and mass $m$, confined in a Penning trap.}
    \label{tab:initial_particle_cond}
\end{table}

\begin{table}[H]
    \centering
        \begin{tabular}[c]{lll}
            Property & Value \\
            \hline
            $B_{0}$ & $1.00 T$ = $9.65 \cross 10^{1} \frac{u}{(\mu s) e}$ \\
            $V_{0}$ & $25.0 mV$ = $2.41 \cross 10^{6} \frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
            $d$ & $500 \ \mu m$ = \\
            \hline
        \end{tabular}
    \caption{Default configuration of the Penning trap, where the value of T and V can be found in table \ref{tab:constants}.}
    \label{tab:penning_config}
\end{table}

\begin{table}[H]
    \centering
        \begin{tabular}[c]{lll}
            Property & Value \\
            \hline
            $q$ & $1.00 T$ & $9.65 \cross 10^{1} \frac{u}{(\mu s) e}$ \\
            $m$ & $25.0 mV$ & $2.41 \cross 10^{6} \frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
            \hline
        \end{tabular}
    \caption{Default configuration of the Penning trap, where the value of T and V can be found in table \ref{tab:constants}.}
    \label{tab:particle_config}
\end{table}

\begin{table}[H]
    \centering
    \begin{tabular}[c]{lll}
        $n_{k}$ & Time steps & Step size \\
        \hline
        $n_{1}$ & $4000$ & $0.0125$ \\
        $n_{2}$ & $8000$ & $0.00625$ \\
        $n_{3}$ & $16000$ & $0.003125$ \\
        $n_{4}$ & $32000$ & $0.0015625$ \\
        \hline
    \end{tabular}
    \caption{Number of steps used in a given simulation $k$, where the step size corresponds to $h_{1} = 50 / n_{k} \mu s$.}
    \label{tab:time_steps}
\end{table}

\subsection{Numbers and units}\label{sec:numbers_units}
\begin{table}[H]
    \centering
    \begin{tabular}[c]{lll}
        Constant & Value & Unit \\
        \hline
        $k_{e}$ (Coulomb) & $1.38935333 \cross 10^{5}$ & $\frac{u (\mu m)^{3}}{(\mu s)^{2} e^{2}}$ \\
        T (Tesla) & $9.64852558 \cross 10^{1}$ & $\frac{u}{(\mu s) e}$ \\
        V (Volt) & $9.64852558 \cross 10^{7}$ & $\frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
        \hline
    \end{tabular}
    \caption{Value of the Coulomb constant ($k_{e}$), and the SI units for magnetic field strength ($T$) and electric potential ($V$). The base units are given by length in micrometre ($\mu m$), time in microseconds ($\mu s$), mass in ($u$), and charge in elementary charge ($e$).}
    \label{tab:constants}
\end{table}

\section{Algorithm implementation and testing}\label{sec:algo}
\subsection{Forward Euler}\label{sec:algo_euler}
For a particle $i$, at time step $j$, the forward Euler method for a coupled system can be expressed as
\begin{align*}
    \mathbf{r}_{i,j+1} &= \mathbf{r}_{i,j} + h \frac{d \mathbf{r}_{i,j}}{dt} = \mathbf{r}_{i,j} + h \mathbf{v}_{i,j} \\
    \mathbf{v}_{i,j+1} &= \mathbf{v}_{i,j} + h \frac{\mathbf{v}_{i,j}}{dt} = \mathbf{v}_{i,j} + h \frac{\mathbf{F}(t_{j},\mathbf{v}_{i,j}, \mathbf{r}_{i,j} )}{m_{i}},
\end{align*}
$m_{i}$ is the mass of the particle, and $h$ is the step length.

\subsection{4th order Runge-Kutta}\label{sec:algo_rk4}
For a particle $i$, at time step $j$, the 4th order Runge-Kutta method for a coupled system can be expressed as
\begin{align*}
    \mathbf{v}_{i,j+1} &= \mathbf{v}_{i,j} + \frac{h}{6} (\mathbf{k}_{\mathbf{v},1,i} + 2 \mathbf{k}_{\mathbf{v},2,i} + 2 \mathbf{k}_{\mathbf{v},3,i} + \mathbf{k}_{\mathbf{v},4,i} ), \\
    \mathbf{r}_{i,j+1} &= \mathbf{r}_{i,j} + \frac{h}{6} (\mathbf{k}_{\mathbf{r},1,i} + 2 \mathbf{k}_{\mathbf{r},2,i} + 2\mathbf{k}_{\mathbf{r},3,i} + \mathbf{k}_{\mathbf{r},4,i}),
\end{align*} 
where
\begin{align*}
    \mathbf{k}_{\mathbf{r},1,i} &= \mathbf{v}_{i,j} \\
    \mathbf{k}_{\mathbf{v},1,i} &= \frac{\mathbf{F}_{i}(t_{j}, \mathbf{v}_{i,j}, \mathbf{r}_{i,j})}{m_{i}}, \\
    \mathbf{k}_{\mathbf{r},2,i} &= \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},1,i}}{2}, \\
    \mathbf{k}_{\mathbf{v},2,i} &= \frac{\mathbf{F}_{i}(t_{j}+\frac{h}{2}, \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},1,i}}{2}, \mathbf{r}_{i,j} + h \frac{\mathbf{k}_{\mathbf{r},1,i}}{2})}{m_{i}}, \\
    \mathbf{k}_{\mathbf{r},3,i} &= \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},2,i}}{2}, \\
    \mathbf{k}_{\mathbf{v},3,i} &= \frac{\mathbf{F}_{i}(t_{j}+\frac{h}{2}, \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},2,i}}{2}, \mathbf{r}_{i,j} + h \frac{\mathbf{k}_{\mathbf{r},2,i}}{2})}{m_{i}}, \\
    \mathbf{k}_{\mathbf{r},4,i} &= \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},1,i}}{2}, \\
    \mathbf{k}_{\mathbf{v},4,i} &= \frac{\mathbf{F}_{i}(t_{j}+h, \mathbf{v}_{i,j} + h \mathbf{k}_{\mathbf{v},3,i}, \mathbf{r}_{i,j} + h \mathbf{k}_{\mathbf{r},3,i})}{m_{i}}.
\end{align*}
In order to find each $\mathbf{k}_{\mathbf{r},i}$ and $\mathbf{k}_{\mathbf{v},i}$, we need to first compute all $\mathbf{k}_{\mathbf{r},i}$ and $\mathbf{k}_{\mathbf{v},i}$ for all particles, then update the particle values in order to compute $\mathbf{k}_{\mathbf{r},i+1}$ and $\mathbf{k}_{\mathbf{v},i+1}$. 

\subsection{Test of implementation}\label{sec:testing}
We implemented a test suite, to validate the implementation during code development. In addition, we have implemented functionality to get informative output while testing the code. Further instructions with code can be found in the github repo \cite{github:repo}.

\end{document}