Project-3/latex/sections/methods.tex

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

\begin{document}
\section{Theory and methods}
% Problem 1
When we study the Penning traps effect on a particle with a charge $q$, we need to consider the forces acting on the particle. We can use Newton's second law \eqref{eq:newton_second} to determine the position of a particle as a function of time. In addition, we introduce the Lorentz force \eqref{eq:lorentz_force}, which describes the force on the particle. The position can be described by %
%
\begin{equation}\label{eq:newton_lorentz}
    m \ddot{\mathbf{r}} = (q \mathbf{E} + q \mathbf{v} \times \mathbf{B}).
\end{equation}
%
Using eq. \eqref{eq:newton_lorentz} we derive the differential equations in sec. \ref{sec:eq_motion}, and can rewrite them 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} %
%
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 sec. \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$.

% Problem 2
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 sec. \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}

% Problem 3
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 rate is given by
\begin{equation*}
    \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 derived expression in \ref{sec:eq_coord}, and 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*} %
When $t \rightarrow \infty$, the upper and lower limits are
\begin{align}
    R_{+} &= A_{+} + A_{-} \\
    R_{-} &= |A_{+} - A_{-}|, 
\end{align}
derived in sec. \ref{sec:bounded_solution}.

% Problem 3
To obtain the bounded solution, we need to consider the expression
\begin{equation*}
    \alpha = (\omega_{+} - \omega_{-}) t +( \phi_{+} -  \phi_{-}).
\end{equation*}
When $t \rightarrow \infty$ the constant phases will not affect the expression, we obtain a bounded solution if we consider the rate such that $|f(t)| < \infty$. That is, we need $\omega_{0}^{2} - 2 \omega_{z}^{2}$ in eq. \eqref{eq:angular_rate} to be real.
\begin{align*}
    \omega_{0}^{2} - 2 \omega_{z}^{2} = 0, \quad \omega_{0} > 2 \omega_{z}
\end{align*}
\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 solution 
Physical properties given by newtons second law \eqref{eq:newton_second}

The particle moves and its position can be determined using newton. where the electric field





\begin{table}[H]
    \centering
        \begin{tabular}[c]{lll}
            Constant & 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}

\subsection{Dealing with a multi--particle system}

For a multi-particles system, we need to modify $\vb{F}$ to account for the 
force of other particles in the system acting upon each other. To do that, we
add another term to $\vb{F}$
\begin{equation}
    \vb{F}_i(t, \vb{v}_i, \vb{r}_i) = q_i \vb{E}(t, \vb{r}_i) + q_i \vb{v}_i \cross \vb{B} - \vb{E}_p(t, \vb{r}_i),
\end{equation}
where $i$ and $j$ are particle indices and
\begin{equation}
    \vb{E}_p(t, \vb{r}_i) = q_i k_e \sum_{j \neq i} 
    q_j \frac{\vb{r_i} - \vb{r_j}}{\left| \vb{r_i} - \vb{r_j} \right|^3}.
    \label{eq:}
\end{equation}

Newton's second law for a particle $i$ is then
\begin{equation}
    \frac{d^2\vb{r}_i}{dt^2} = \frac{\vb{F}_i\left(t, \frac{d\vb{r}_i}{dt}, \vb{r_i}\right)}{m_i},
    \label{eq:newtonlaw2}
\end{equation}
We can then rewrite the second order ODE from equation~\ref{eq:newtonlaw2} 
into a set of coupled first order ODEs. 
We now rewrite Newton's second law of motion as
\begin{equation}
    \begin{split}
        \frac{d\vb{r}_i}{dt} &= \vb{v}_i \\
        \frac{d\vb{v}_i}{dt} &= \frac{\vb{F}_i(t, \vb{v}_i, \vb{r}_i)}{m_i}.
    \end{split}
    \label{eq:coupled}
\end{equation}

\subsection{Forward Euler}

For a particle $i$, the forward Euler method for a coupled system is 
expressed as
\begin{equation}
    \begin{split}
        \vb{r}_{i,j+1} &= \vb{r}_{i,j} + h \cdot \frac{d\vb{r}_{i,j}}{dt} = \vb{r}_{i,j} + h \cdot \vb{v}_{i,j} \\
        \vb{v}_{i,j+1} &= \vb{v}_{i,j} + h \cdot \frac{\vb{v}_{i,j}}{dt} = \vb{v}_{i,j} + h \cdot \frac{\vb{F}\left( t_{j}, \vb{v}_{i,j}, \vb{r}_{i,j} \right)}{m},
    \end{split}
\end{equation}
for particle $i$ where $j$ is the current time step of the particle, 
$m$ is the mass of the particle, and $h$ is the step length.

When dealing with a multi-particle system, we need to ensure that we do not
update the position of any particles until every particle has calculated their
next step. An easy way of doing this is to create a copy of all the particles,
then update the copy, and when all the particles have calculated their next 
step, simply replace the particles with the copies. 
Algorithm~\ref{algo:forwardeuler} provides an overview on how that can be achieved. % Make this better

\begin{figure}[H]
    \begin{algorithm}[H]
    \caption{Forward Euler method}
    \label{algo:forwardeuler}
        \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}

\subsection{4th order Runge-Kutta}

For a particle $i$, we can express the 4th order Runge-Kutta (RK4) method as
\begin{equation}
    \begin{split}
        \vb{v}_{i,j+1} &= \vb{v}_{i,j} + \frac{h}{6} \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) \\
        \vb{r}_{i,j+1} &= \vb{r}_{i,j} + \frac{h}{6} \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),
    \end{split}
\end{equation}
where
\begin{equation}
    \begin{split}
        \vb{k}_{\vb{v},1,i} &= \frac{\vb{F}_i(t_j, \vb{v}_{i,j}, 
            \vb{r}_{i,j})}{m} \\
        \vb{k}_{\vb{r},1,i} &= \vb{v}_{i,j} \\
        \vb{k}_{\vb{v},2,i} &= \frac{\vb{F}_i(t_j+\frac{h}{2}, \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2}, \vb{r}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{r},1,i}}{2})}{m} \\
        \vb{k}_{\vb{r},2,i} &= \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2} \\
        \vb{k}_{\vb{v},3,i} &= \frac{\vb{F}_i(t_j+\frac{h}{2}, \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},2,i}}{2}, \vb{r}_{i,j} 
            + h \frac{\cdot \vb{k}_{\vb{r},2,i}}{2})}{m} \\
        \vb{k}_{\vb{r},3,i} &= \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},2,i}}{2} \\
        \vb{k}_{\vb{v},4,i} &= \frac{\vb{F}_i(t_j+h, \vb{v}_{i,j} 
            + h \cdot \vb{k}_{\vb{v},3,i}, \vb{r}_{i,j} 
            + h \cdot \vb{k}_{\vb{r},3,i})}{m} \\
        \vb{k}_{\vb{r},4,i} &= \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2}.
    \end{split}
\end{equation}

In order to find each $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$, 
we need to first compute all $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$ 
for all particles, then we can update the particles in order to compute 
$\vb{k}_{\vb{r},i+1}$ and $\vb{k}_{\vb{v},i+1}$. In order for the algorithm
to work, we need to save a copy of each particle before starting so that we
can update the particles correctly for each step.

This approach would require 8 loops to be able to complete the calculation since
we cannot update the particles until after all $\vb{k}$ values have been 
computed, however if we create a temporary array that holds particles, we can
put the updated particles in there, and then use that array in the next loop,
and would reduce the required amount of loops down to 4.

\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} + \vb{k}_{\vb{r},2,i} 
                + \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} + \vb{k}_{\vb{v},2,i} 
                + \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}
    $\vb{F}$ in the algorithm does not take any arguments as it uses the 
    velocities and positions of the particles inside the array $particles$ to
    calculate the total force acting on particle $i$.
\end{figure}


\subsection{Testing the simulation}

\subsection{Relative error and error convergance rate}

\subsection{Tools}
The numerical methods implemented in C++, are parallelized using OpenMP \cite{openmp:2018}. We used the Python library matplotlib \cite{hunter:2007:matplotlib} to produce all the plots, and seaborn \cite{waskom:2021:seaborn} to set the theme in the figures.

\end{document}