\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}