\documentclass[a4paper,oneside]{saip}
\usepackage{comment}
\usepackage{float}
\usepackage{subcaption}
\usepackage{xcolor}
\usepackage{graphicx}
\newcommand{\p}{\partial}
\renewcommand{\thesubfigure}{Fig \thefigure(\alph{subfigure}):}
\captionsetup[subfigure]{labelformat=simple}
\begin{document}


\title{
Solution of the Two Dimensional Schr\"odinger Equation using 16/8 Daubechies Wavelet Scaling Functions as Basis Set
}

\author{Moritz Braun \,\orcidlink{0000-0001-8710-7561}$^{1}$,
   Obiageli  L. Ezenwachukwu\,\orcidlink{0000-0001-5044-8565}$^{1}$}

\affil{$^{1}$\,Department of Physics, University of South Africa , PO Box 392, 003 UNISA, South Africa}

\email{\url{moritz.braun@gmail.com}}



\begin{abstract}
In this contribution we demonstrate the suitability of the 16/8 Daubechies wavelets scaling function for solving the two dimensional Schr\"odinger equation using the example of the isotropic two dimensional harmonic oscillator.  We show that calculations with this basis set converge quickly both for the 
eigenvalues as well as the wave functions. A convergence order of approximately 9  is established.
\end{abstract}


\section{Introduction}\label{sec:intro}
Symlet8 or sym8 wavelets are  specialized wavelets within wavelet theory valued for their near-symmetry, smoothness, and high-order vanishing moments\cite{Graps1995}.
The sym8 wavelet is applied within the Discrete Wavelet Transform (DWT) to split signals into detailed and approximation components.
The improved symmetry of sym8 is preferred over db8 in image processing, edge detection, and analysis where signal phase information is important\cite{Amid2023}.\\
In this contribution we make use of the associated scaling functions modified to satisfy periodic boundary bonditions as 1D basis set to solve the two dimensional Schr\"odinger equation making use of higher order  Gauss Legendre integration within the Galerkin approach employing 2D basis functions  constructed as tensor product of the one-dimensional ones.
\subsection{Definition of the Basis Set}
We use two dimensional basis  functions based on shifted and scaled scaling sym8 wavelet  functions  
\begin{equation}	
f_i(x)=\phi\left( x/h+x_{\rm max} /h -i\right)/\sqrt{h}
\end{equation}
with  $h=x_{\rm max}/N$, where $x_{\rm max}$ is the distance from the origin to he boundary.\\
We however modify the functions by giving them periodic boundary conditions via
\begin{equation}
	g_i(x)=\sum_{k=-K_{\rm max},} ^ {K_{\rm max}-1}
f_i(x+kL)\, 
\end{equation}
with $L=2x_{\rm max}$.  For $x_{\rm max}$ large enough very this will have a negligible  effect. It can be shown, 
that these basic functions  are orthonormal.
The  two dimensional basis functions are given by
	\begin{equation}
	F_{ij}(x,y)=g_{i}(x) g_{j}(y)
     \end{equation}
\subsection{Hamiltonian}
The Hamiltonian for the  two-dimensional Schr\"odinger equation is given by 
	\begin{equation}
		H=-\frac{\p^2}{\p x^2}
		  -\frac{\p^2}{\p y^2} 
		+V(x,y)
	\end{equation}
and the  matrix element is obtained as  
\begin{equation}
	h_{ij,kl}=   \int\!\!\!  \int dx dy F_{ij}
	\left[
	-\frac{\p^2}{\p x^2}
		  -\frac{\p^2}{\p y^2} 
		+V(x,y)
		\right]
	F_{kl}
\end{equation}
Because of the product structure of 2D basis functions this becomes
\begin{equation}
	h_{ij,kl}=q_{ik}\delta_{jl}+\delta_{ik} q_{jl}+v_{ij,kl}
\end{equation}
with 
\begin{equation}
v_{ij,kl}=\int\!\!\!  \int dx dy F_{ij}V(x,y)F_{kl}
\end{equation}
and
\begin{equation}
	q_{ij}=\int dx g_i'(x) g_j'(x)
\end{equation}
If $V$ is the sum of a function of $x$ and a function of $y$ the potential also factorizes into a sum of two  tensor
  products. However this assumption is not made for the  calculation.
\section{Details of Calculation}
The double indices like $ij$ will be transformed into
one-dimensional indices before solving  the eigenvalue probem.\\
The calculations were done with a Python code of 130 lines using the modules pywt\cite{Lee2019}, scipy\cite{jones2001} and  numpy\cite{HarrisCR2020}.\\ The integrals above are evaluated via repeated Gauss Legendre integration of order 11 in 1D  and 2D. The tensor products were evaluated via the Kronecker numpy.kron function in Python. The eigenvalue problem was solved using the module eigh from scipy.linalg. The Python source code of 130 lines is available on request from moritz.braun@gmail.com.
\section{Results}
We consider the two-dimensional harmonic oscillator
	$V(x,y)=x^2+y^2$ with the eigenvalues $2\nu$ 
	for $\nu = 1,2,\cdots$ and degeneracy  
	$\nu$ for level $\nu$. 
	\\
In the table below the resulting 10 lowest eigenvalues are given
for $x_{\rm max}=8$ and for $N=12,16$ and 20 from left to right:
\medskip 
\begin{table}[H]
\centering % the table itself is centred
\setlength{\tabcolsep}{0.5em} % for the horizontal padding
\renewcommand{\arraystretch}{1.2}% for the vertical padding
\begin{center}
%	\caption{Results of 10 Lowest Eigenvalues}
\begin{tabular}{l c c c}
\hline
\hline
$\nu$ &        &           &    \\
\hline
1 & 2.00029488 &  2.00001603 &  2.00000158\\
2 & 4.00191520 &  4.00011385 &  4.00001107\\
3 & 4.00191520 &  4.00011385 &  4.00001107\\
4 & 6.00353437 &  6.00021149 &  6.00002052\\
5 & 6.01034282 &  6.00070568 &  6.00006927\\
6 & 6.01034396 &  6.00070585 &  6.00006930\\
7 & 8.01196025 &  8.00080302 &  8.00007865\\
8 & 8.01196025 &  8.00080302 &  8.00007865\\
9 & 8.03803395 &  8.00307480 &  8.00031353\\
10 & 8.03803395 &  8.00307480 &  8.00031353\\
\hline
\end{tabular}
\end{center}
\caption{Lowest ten eigenvalues obtained for N=12,16 and 20 from left to right}\end{table}

\vspace{5mm}
\noindent
The convergence is quite fast.\\
We define $dE_1(N)=E_1(N)-2$. A least square fit to
    \begin{equation}
      \log dE_1(N) = a +b \log N
      \end{equation}
      for $N=8,10,12,14,16,18$ and $20$
      results in the plot below:
      \begin{figure}[H]
      \centering
    \includegraphics[width=0.75\textwidth]{gnuplotFitdE1.png}
        \caption{Linear Least Square Fit of 
    	$a+b \log N$ to $\log dE_1(N)$ }   
   \end{figure}
\subsection{ Plots of wave functions}
Below  in Fig. 2 a) - 2	 h) heat maps for  the wave functions of the lowest 8 states of the 2D harmonic oscillator are shown.
%%%%%%%
\begin{figure}[H]
\centering
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave0.png}
    \caption{Heat map of ground state wave function}
    \label{fig:wf0}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave1.png}
    \caption{Heat map of first excited state}
    \label{fig:wf1}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave2.png}
    \caption{Heat map of second excited state}
    \label{fig:wf2}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave3.png}
    \caption{Heat map of third excited state}
    \label{fig:wf3}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave4.png}
    \caption{Heat map of fourth excited state}
    \label{fig:wf4}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave5.png}
    \caption{Heat map of fifth excited state}
    \label{fig:wf5}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave6.png}
    \caption{Heat map of sixth excited state}
    \label{fig:wf6}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=\linewidth]{wave7.png}
    \caption{Heat mat of seventh excited state}
    \label{fig:wf7}
\end{subfigure}
\end{figure}
\noindent
Below a heatmap of the difference between the numerical and anylytical ground state wave function for $N=20$ is shown:\\
\begin{figure}[H]
\includegraphics[width=0.7\textwidth]{dwf_N=20.png}
\caption{Heat map of difference between numerical and analytical ground state wave function}
\end{figure}
\section{Conclusions}
The two dimensional isotropic  harmonic oscillator as a test case has been   solved with wavelet scaling code. The convergence is quite rapid. In future we intend
to use this two dimensional  code for more interesting  potentials.
\section*{Acknowledgements}
This work was supported by the University of South Africa.  
\bibliographystyle{IEEEtran}
\bibliography{PhDRef}
\end{document}
