\subsection{Discretization of the Schrödinger equation}
All numerical methods for solving differential equations rely on the discretization of continuous variables. In this project, we will discretize both space and time. We will begin with the two-dimensional time-dependent Schrödinger equation:
\begin{align}
    i\hbar \frac{\partial}{\partial t} \psi(x,y,t) &= -\frac{\hbar^2}{2m} \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) \psi(x,y,t) \nonumber\\&+ V(x,y) \psi(x,y,t). 
\end{align}
We define a grid with spacing \(h\) in both the \(x\) and \(y\) directions to discretize space. The grid points are given by \(x_i = i h\) and \(y_j = j h\), where \(i\) and \(j\) are integers. The wave function at these grid points is denoted as \(\psi_{i,j}(t) = \psi(x_i, y_j, t)\) . We also discretize time with a time step of size \(\Delta t\), such that \(t_n = n \Delta t\), where \(n\) is an integer. The wave function at time \(t_n\) is denoted as \(\psi_{i,j}^n = \psi_{i,j}(t_n)\).
For simplicity, we set both \(\hbar = \qty{1}{\per \s}\) and \(m = 1\). This yields the dimensionless form of the Schrödinger equation:
\begin{align}
     i \frac{\partial}{\partial t} u(x,y,t) &= - \frac{\partial^2}{\partial x^2} u(x,y,t) - \frac{\partial^2}{\partial y^2} u(x,y,t) \nonumber\\&+ v(x,y) u(x,y,t). 
\end{align}
In this form, we substitute \(u\) for \(\psi\) and \(v\) for \(V\) to eliminate constants. We can then reobtain the physical probability density, $P = \psi^\ast \psi$, by appropriately normalizing the wave function:
\begin{equation}
    P(x,y,t) = \frac{u^\ast(x,y,t) u(x,y,t)}{\int \int u^\ast(x,y,t) u(x,y,t) \, dx \, dy}.
\end{equation}
In practice, the integral in the denominator is equal to the sum of all grid points multiplied by the area element \(h^2\).

Using the Crank-Nicolson scheme, we can approximate the time evolution of the wave function. The Crank-Nicolson method is a second-order accurate, unconditionally stable implicit finite difference method in both space and time. The update rule for the wave function using the Crank-Nicolson scheme is given by the following formula \cite{liuCrankNicolsonMethod2022}:
\begin{equation}
    \left( I + \frac{i \Delta t}{2} H \right) u^{n+1} = \left( I - \frac{i \Delta t}{2} H \right) u^n,
\end{equation}
where \(I\) is the identity matrix, \(H\) is the discretized Hamiltonian operator, and \(u^n\) and \(u^{n+1}\) are the wave functions at time steps \(n\) and \(n+1\), respectively. The Hamiltonian operator \(H\) can be discretized using finite difference approximations for the second derivatives:
\begin{align}
    H_{i,j} &= -\frac{1}{h^2} \left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j} \right) \nonumber\\&+ v_{i,j} u_{i,j}. 
\end{align}
By defining the constant as \(r \equiv \frac{i \Delta t}{2 h^2}\), we can rewrite the update rule as follows:
\begin{align}
    &u_{i,j}^{n+1} + r \left( u_{i+1,j}^{n+1} + u_{i-1,j}^{n+1} + u_{i,j+1}^{n+1} + u_{i,j-1}^{n+1} - 4u_{i,j}^{n+1} \right) \nonumber\\&= u_{i,j}^n - r \left( u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n \right)\nonumber\\& - i \Delta t v_{i,j} \frac{u_{i,j}^{n+1} + u_{i,j}^n}{2}. 
\end{align}
We can rearrange this equation to form a linear system that can be solved for \(u^{n+1}\) at each time step. We define the index mapping \(k = i + j M\), where \(M\) is the number of grid points in one dimension and \(M = \frac{1}{h} + 1\). This allows us to rewrite the two-dimensional grid as a one-dimensional array \cite{anderskvellestadProject5FYS3150}. The resulting linear system can be expressed in matrix form as follows:
\begin{equation}
    A u^{n+1} = B u^n,
\end{equation}
where \(A\) and \(B\) are matrices that depend on the discretization parameters and the potential \(v\). Only the internal points of the problem are contained in \(A,\,B,\,u^n\). Boundary conditions are applied separately. For this project, we will use Dirichlet boundary conditions and set the wave function to zero at the grid boundaries.
Matrix  \(A\) is tridiagonal with additional diagonals due to the two-dimensional nature of the problem. The matrix \(B\) has a similar structure. To solve for \(u^{n+1}\), first calculate $b^n = B u^n$, then solve the linear system \(A u^{n+1} = b^n\) using an appropriate numerical method, such as LU decomposition or iterative solvers.

\subsection{Initial conditions and potential}
To simulate the time evolution of the wave function, we must specify the initial conditions and the potential, \(v(x,y)\). A common choice for the initial wave function is a Gaussian wave packet, defined as follows \cite{anderskvellestadProject5FYS3150}:
\begin{align}
    u(x,y,0) &= \exp\left( -\frac{(x-x_C)^2}{2\sigma_x^2} - \frac{(y-y_C)^2}{2\sigma_y^2}\right) \nonumber\\&\cdot \exp \left( i(p_x x + p_y y) \right),
\end{align}
where \((x_C, y_C)\) is the center of the wave packet and where  \(\sigma_x, \sigma_y)\) are the widths of the packet in the \(x\) and \(y\) directions, respectively, and where the initial momenta in the respective directions are represented by \((p_x, p_y)\). The potential, \(v(x,y)\), can be chosen based on the physical system being modeled. For example, a double-slit-like potential can be used. To implement slits in the potential grid, we use the algorithm outlined in \cref{alg:add_slits}. This algorithm modifies the potential grid, \(V\), by adding walls of a specified thickness and position while leaving apertures for the slits. An initialization function finalizes the potential setup by calculating the matrices \(A\) and \(B\) based on the modified potential.


\begin{algorithm}[H]
\caption{Add Slits to Potential Grid}
\label{alg:add_slits}
\begin{algorithmic}[1]
\Require $h$, $M$, $V$, $V_{\text{wall}}$, $ \text{wall\_thickness}$, $\text{wall\_position}$, $\text{slit\_aperture}$, $\text{slit\_separation}$, $\text{num\_slits}$
\State $x_{\text{start}} \gets \left\lfloor \dfrac{\text{wall\_position} - 0.5\,\text{wall\_thickness}}{h} \right\rfloor$
\State $x_{\text{end}}   \gets \left\lfloor \dfrac{\text{wall\_position} + 0.5\,\text{wall\_thickness}}{h} \right\rfloor$
\State $y_c \gets M/2$
\State $\text{slit\_half} \gets \left\lfloor \dfrac{0.5\,\text{slit\_aperture}}{h} \right\rfloor$
\State $\text{pitch} \gets \text{slit\_aperture} + \text{slit\_separation}$
\State $\text{pitch}_{\text{idx}} \gets \text{pitch}/h$

\For{$i = x_{\text{start}}$ \textbf{to} $x_{\text{end}}$}
    \For{$j = 0$ \textbf{to} $M-1$}
        \State $\text{in\_slit} \gets \text{false}$
        \For{$s = 0$ \textbf{to} $\text{num\_slits}-1$}
            \State $\Delta \gets \left(s - \dfrac{\text{num\_slits}-1}{2}\right)\,\text{pitch}_{\text{idx}}$
            \State $c_s \gets \left\lfloor y_c + \Delta + 0.5 \right\rfloor$ \Comment{rounded center index}
            \If{$j \in [\,c_s - \text{slit\_half}, \; c_s + \text{slit\_half}\,]$}
                \State $\text{in\_slit} \gets \text{true}$
                \State \textbf{break}
            \EndIf
        \EndFor
        \If{\textbf{not} $\text{in\_slit}$}
            \State $V[i,j] \gets V_{\text{wall}}$
        \EndIf
    \EndFor
\EndFor

\State \Call{initialize\_potential}{}

\end{algorithmic}
\end{algorithm}

\subsection{Implementation details}
To optimize performance, the numerical methods were implemented in C++, compiled with optimization flags, and compiled using the \texttt{g++} compiler. Version \texttt{g++ (GCC) 15.2.1 20251111 (Red Hat 15.2.1-4)} was used for compilation. The Armadillo library  \cite{sandersonArmadilloEfficientFramework2025,sandersonPracticalSparseMatrices2019} was used for all linear algebra operations because it provides efficient implementations of matrix and vector operations. Matrices \(A\) and \(B\) are constructed as sparse \cite{sandersonPracticalSparseMatrices2019}, complex matrices to save memory and improve computational efficiency. We solved the linear system \(A u^{n+1} = b^n\) using the \texttt{spsolve} function from the Armadillo library. This function is optimized for sparse matrices and uses an implementation provided by SuperLU \cite{demmelSuperLUUsersGuide1999}.

To allow for easy experimentation with different parameters, the program reads a configuration file in TOML format during runtime. The configuration file specifies parameters in four main categories:
Discretization parameters: \(h, \Delta t\), and total simulation time \(T\); Initial wave packet parameters (center position, width, and momentum); Potential parameters (wall height, wall thickness, slit aperture, slit separation, and number of slits); Output parameters (output file name). The program uses Armadillo's built-in functions to parse the wave function results at each time step and write them to a file in \texttt{arma\_ascii} format.

All analyses and visualizations of the results are performed in Python using the NumPy \cite{harrisArrayProgrammingNumPy2020} and Matplotlib \cite{hunterMatplotlib2DGraphics2007} libraries. The probability density is computed from the wave function data, and various plots are generated to illustrate the time evolution of the wave packet and the effects of the potential. A self-developed parser for the Armadillo ASCII format reads the output files generated by the C++ program, as described in \cref{alg:pythonparser}.

\begin{algorithm}[H]
\caption{Load Armadillo ASCII Complex Cube}\label{alg:pythonparser}
\begin{algorithmic}[1]

\Procedure{LoadCube}{path, slices\_first}

    \State txt $\gets$ read file at \textit{path}
    \State lines $\gets$ split(txt)
    \State header $\gets$ first two lines
    \State $(n_{\text{rows}}, n_{\text{cols}}, n_{\text{slices}}) \gets$ parse(header[1])

    \State rest $\gets$ remaining lines

    \For{each line in rest}
        \State remove all ``('' and ``)''
        \State apply regex 
        \[
           \texttt{([+-]?\d[\d.eE+-]*)\ ,\ ([+-]?\d[\d.eE+-]*)}
        \]
        which captures two floating-point values \(r\) and \(i\), possibly with
        signs or scientific notation
        \State replace with 
        \[
           \texttt{\$1+\$2j}
        \]
        yielding Python-style complex numbers \(r + ij\)
        \State append transformed line to fixed
    \EndFor

    \State arr$_{2D}$ $\gets$ \Call{numpy.loadtxt}{fixed}
    \State arr$_{3D}$ $\gets$ reshape(arr$_{2D}$, $(n_{\text{slices}}, n_{\text{cols}}, n_{\text{rows}})$)

    \If{slices\_first}
        \State \Return arr$_{3D}$
    \Else
        \State \Return transpose(arr$_{3D}$)
    \EndIf

\EndProcedure

\end{algorithmic}
\end{algorithm}

\subsubsection*{Use of Artificial Intelligence}
In this project, artificial intelligence in the form of large language models (LLMs) was employed to assist with various stages of the research and development process. The LLMs generated initial drafts of code snippets that the researcher reviewed and refined for correctness and efficiency. GitHub Copilot, an LLM, was used for this task. LLMs were also used to edit the report's language by providing grammar and style suggestions. For this task, the DeepL Write service by DeepL GmbH was used.

