We evaluate the accuracy of the resulting solutions by assessing the conservation of probability in the system with and without a double slit, i.e., with \(V_\text{wall} = 0\)  and in the presence of such a potential. We simulate the systems for a total time of \qty{0.008}{\s} in steps of \(\Delta t = \qty{2.5e-5}{\s}\). Spatial discretization is carried out using a step size of \(h = \num{0.005}\). In the setup with a non-zero potential wall, the initial wave packet is centered at \(x=\num{0.25}\) and \(y=\num{0.5}\), with widths of \(\sigma_x=\num{0.05}\) and \(\sigma_y=\num{0.10}\), respectively.  Its momenta are \(p_x=\num{200.0}\) and \(p_y=\num{0.0}\), resulting in strongly directed motion along \(x\). 
A double-slit barrier with a thickness of \(\num{0.02}\) is placed at \(x=\num{0.5}\).  
Each slit has an aperture of \(\num{0.05}\) and is separated by \(\num{0.05}\). The wall potential is effectively infinite at \(V_{\text{wall}}=\num{1.0e10}\). In the case of a constant potential, the wave packet has widths of \(\sigma = \num{0.05}\) in both \(x\) and \(y\). Otherwise, the setup is identical. The resulting wave functions are evaluated for loss of probability.
\begin{equation}
    \delta = 1 - \iint u^\ast(x,y,t) u(x,y,t) \mathrm{d}x \mathrm{d}y,
\end{equation}
where \(\iint u^\ast(x,y,0) u(x,y,0) \mathrm{d}x \mathrm{d}y = 1\) was fixed using normalization. In a perfect physical system, this quantity is conserved; however, due to limited numerical accuracy, especially in the derivatives and solving of the linear system, it is not conserved. Therefore, it is a good indication of the algorithm's potential accuracy. \Cref{fig:prob_conservation} shows the resulting analysis. The maximum deviation is about \num{2e-15}, demonstrating performance that is only one order of magnitude worse than floating-point precision. Notably, the accuracy of the system without a double slit is lower than the corresponding values with the perturbation. Additionally, the simulation's dependence on potential is evident, precluding any general attestation of performance for arbitrary potentials. Nevertheless, it can be assumed that the simulation will produce very accurate results for many potentials, attesting to the numerical accuracy of the results.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{include/probability_conservation.pdf}
    \caption{Loss of total probability $\delta$ in dependence on the simulation duration for two systems of different potential.}
    \label{fig:prob_conservation}
\end{figure}

In addition to the quantitative evaluation of performance, we conduct a qualitative analysis of the system's evolution. To correctly assess the evolution of the wave packet, we perform a double slit experiment on a system for a total duration of \qty{0.002}{\s}. Three time points are sampled to demonstrate the wave packet's interference with itself behind the double slit. The results of this analysis are presented in \cref{fig:t_0_evolution,fig:t_1_evolution,fig:t_2_evolution} in Appendix \cref{app:evolv}.
We observe the expected qualitative behavior with respect to interference and the correct diffusion of the wave function. Therefore, we can conclude that the simulation framework behaves correctly. To evaluate the evolution of the probability distribution, we produced animations of the entire runtime for each simulation. These animations can be accessed via the repository of this paper, though they are not essential to understanding the paper.

Lastly, we demonstrate the flexibility of our approach to simulating the numerical evolution of such a system by simulating the probability of detecting a particle on a screen inserted into the system at a specific time. We computed the probability density function (PDF) for a screen position of \(y = \num{0.8}\) at \(T = \qty{0.002}{\s}\) for three systems: a single slit, a double slit, and a triple slit. The detection PDF is given by the slice
\begin{equation}
    P(x) = \frac{u^\ast(x,\num{0.8},\qty{0.002}{\s}) u(x,\num{0.8},\qty{0.002}{\s})}{\int \mathrm{d}x u^\ast(x,\num{0.8},\qty{0.002}{\s}) u(x,\num{0.8},\qty{0.002}{\s})},
\end{equation}
where the integral is evaluated as a sum over the total time slice. The resulting PDF is displayed in \cref{fig:screen}. Once again, the interference of the wave packet with itself is clearly visible, leading to the expected result. More importantly, this shows how such a numerical simulation can be used to compute arbitrary metrics because the entire numerical solution of the wave function is saved to a file that can be repurposed in multiple ways.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{include/detector_screens.pdf}
    \caption{PDF of detecting the wave packet at a certain position on a screen inserted at $y=\num{0.8}$ and $T = \qty{0.002}{\s}$.}
    \label{fig:screen}
\end{figure}