\documentclass[a4paper,oneside]{saip}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% defaults, keep as is, don't add anything that changes the document's defaults.
\usepackage{helvet}          % to have better matching fonts between Latex and Word (Title/Headings)
\usepackage{times}           % to have better matching fonts between Latex and Word (bulk text)
\usepackage{courierten}   % to have better matching fonts between Latex and Word (monospaced)
\usepackage[T1]{fontenc} % needed to deal with these fonts
\usepackage{amsmath,amssymb,amsfonts,cases}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% common packages needed, add or comment as needed
\usepackage{graphicx, subcaption,booktabs}  %to include figures
\usepackage{tikz}
\usetikzlibrary{positioning}
\usetikzlibrary{arrows, calc}
\usepackage{url}
\usepackage{geometry} %	need: 
\geometry{%showframe, 
top=40.0mm, bottom=27.0mm, left=25.0mm, right=25.0mm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional packages, you can add others as well for your own needs
%\usepackage{xcolor}
%\usepackage{algorithm2e}
%\usepackage[version=4]{mhchem} %for isotopes and all with \ce{^{15}_{8}O_{7}}
%\usepackage{amsmath}
%\usepackage{amssymb}
%\usepackage{stmaryrd}
%\usepackage{mathrsfs}
%\usepackage{mathabx}
%\usepackage{upgreek}
%\usepackage{comment}
%\usepackage{verbatim}
%\usepackage{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO NOT USE the hyperref package, your PDF will be included in a bigger document and the links will be broken.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}


\title{Eigenvalue determination for a Toy and Woods-Saxon potentials using an unsupervised PINN}

\author{J T Tshipi$^{1}$ and A S Cornell$^{2}$}

\affil{$^1$Department of Physical and Earth Sciences, Sol Plaatje University, Kimberley, South Africa}
\affil{$^2$University of Johannesburg, Department of Physics, PO Box 524, Auckland
Park, Johannesburg, 2006, Gauteng, South Africa}

%only put the email of the main corresponding author, not of all authors.
\email{\url{tshegofatso.tshipi@spu.ac.za}}



\begin{abstract}
Traditional numerical methods are widely used to determine eigenvalues in quantum few-body systems. However, comparatively fewer studies have explored modern machine learning approaches such as Physics-Informed Neural Networks (PINNs) for this purpose. In this work, we apply PINNs to solve the two-body radial Schrödinger equation and compute the low-lying S-wave bound states for a toy molecular potential and the Woods-Saxon potential. 
The PINN was observed to initially converge to the ground state before transitioning to higher excited states, regardless of the starting point, which is consistent with theory.
For the molecular potential, the predicted first and second bound-state energies achieved relative errors of $1.482\times10^{-2}\% $ and $6.37\times 10^{-3}\%$, respectively. In a preliminary study, the Woods-Saxon potential achieved the relative errors for the first three bound states of $2.99\times 10^{-1}\%$, $2.38\times 10^{-1}\%$ and $14.85\%$. These results are promising and demonstrate that PINNs have the potential to offer a viable alternative to traditional numerical methods for solving bound-state quantum eigenvalue problems. Especially with further optimisation of the neural network architecture.
\end{abstract}


\section{Introduction}

Bound states play an important role in our understanding as well as predicting the stable of matter from nucleons to atoms. These states occur when particle(s) are confined to a fixed region of space by an attractive force, resulting in negative energies and localized wave functions; that is, they are non-zero inside the potential and zero elsewhere. In quantum few-body problems, bound states are determined by solving the Schrödinger equation to obtain eigenvalues, corresponding to bound state energies, and eigenfunctions which describe the wave function. Determining these bound states, numerically, typically involves search algorithms, such as bisection~\cite{bisectionM} or Newton-Raphson~\cite{LanduComp} methods, to find energies alongside ordinary differential equation (ODE) solvers (like Runge–Kutta~\cite{rkmethods} or the Numerov method~\cite{numerov}) to find the corresponding wave functions.

In this proceeding we shall explore the application of the Physics-Informed Neural Network (PINN) framework to calculate low-lying S-wave bound states. This approach was successfully used in~\cite{usnn} to determine bound state eigenpairs for simple textbook potentials, i.e. the (in)finite square well and harmonic oscillator. Here we shall we investigate its efficiency by solving for bound state eigenpairs using a toy molecular potential~\cite{1996Paper,1997Paper} and the Woods-Saxon potential~\cite{woodsax}. These two potentials offer a more realistic representation of the nuclear potential since they are short-ranged, have finite width and finally their shapes at the edge changes smoothly~\cite{krane}.  

\section{Theory of Bound states}
We shall consider two spinless neutral particles that interact via a short-range central potential, the radial Schr\"{o}dinger equation describing the relative motion of the particles with reduced mass $\mu$, is
\begin{equation}
 \label{singleChannelSE}
 \left[ -\frac{\hbar^2}{2\mu}\frac{d^2}{dx^2}+\frac{\ell\left(\ell+1 \right)}{r^2}+V(r)
 \right]\psi_{\ell}(x) = E\psi_{\ell}(x). 
\end{equation}
As we are only interested in determining the S-wave bound-states, i.e. $\ell = 0$, the centrifugal barrier term will become zero, thus increasing the probability of forming a bound state. The Schr\"{o}dinger equation~(\ref{singleChannelSE})  becomes   
\begin{equation}
 \label{singleChannelSE_SWave}
 \left[ -\frac{\hbar^2}{2\mu}\frac{d^2}{dx^2}+V(r)\right]\psi(r) = E\psi(r). 
\end{equation}    
To obtain the physical solutions of~(\ref{singleChannelSE_SWave}) we will impose boundary conditions at $r = 0$ and $r \rightarrow \infty$ making this a Dirichlet problem. These kinds of problems can only have solutions at certain discrete energies, since it is not always possible to satisfy the physical boundary condition~\cite{book}. At the origin, we have $\psi(0) = 0$ this is true for all types of solutions e.g bound, resonance and scattering states. For bound states, the particle $\mu$ is trapped and thus confined to a particular region in space by an attractive force. Therefore the probability of finding the particle beyond the walls of the potential decays to zero. As such the boundary conditions we require are
\begin{equation}
    \label{BoundaryConditions}
    \psi(0) = 0 \qquad \mbox{ and } \qquad \psi(r) \xrightarrow[r \rightarrow \infty]{} 0 . 
\end{equation}
We will now use the unsupervised neural network approach proposed in~\cite{usnn} to determine the bound state energies and the corresponding wave functions, where this approach is capable of solving~(\ref{singleChannelSE_SWave}) when both $E$ and $\psi(r)$ are unknowns.

\section{Potentials} 
\subsection{Toy potential} 
The first potential we will consider is the toy molecular potential shown in Figure~\ref{fig:potentials}. The potential~\cite{1996Paper, 1997Paper} 
\begin{equation}   
    \label{toyPot}
    V (r) = 5.0\exp\left[ -0.25\left(r-3.5\right)^2\right] - 8\exp\left(r-0.2r^2\right),
\end{equation}
is primarily used as a testbed for assessing the performance of new theoretical models as it generates both bound and resonance states.
We solve~(\ref{singleChannelSE_SWave}) using atomic units~\cite{AtomicUnits} i.e. the energies are measured in MeV, distance in femtometre(fm) and $\hbar^2/(2\mu) = 1/2 \mathrm{MeV} \cdot \mathrm{fm}^{-2}$. To determine the bound state eigenvalues and the wave functions for this potential, we use the boundary conditions $\psi_{0}(0) = \psi_{0}(5) = 0$.

\subsection{Woods-Saxon potential}
The Woods-Saxon potential, see Figure~\ref{fig:potentials},  is used to model the mean field potential that the nucleons feel inside the nucleus~\cite{krane}.  
The standard form of the Wood-Saxon potential~\cite{woodsax} is given by
\begin{equation}
    \label{WoodSaxPontential}
     V(r) = - \frac{V_0}{1 + \exp\left(\frac{r-R}{d}\right)}. 
\end{equation}
Here we will use the parameters found in~\cite{Mmusi}, $R = r_0A^{\frac{1}{3}}$ with $r_0 = 1.1$fm and $d=0.6$fm. The value of the strength $V_0=31.20$ MeV. The boundary conditions  $\psi_{0}(0) = \psi_{0}(13) = 0$ will be used to determine bound states eigenpairs from equation~(\ref{singleChannelSE_SWave}).

\section{Computing bound states with PINN}
In setting up our neural network to solve to solve for the eigenvalues of the S-wave Schr\"{o}dinger equation~(\ref{singleChannelSE_SWave}) for the potentials~(\ref{toyPot}) and~(\ref{WoodSaxPontential}), we shall follow the approach detailed in \cite{usnn}, which the reader may consult for more details. In this approach the trial wave function of~(\ref{singleChannelSE_SWave}) is given by 
\begin{equation}
    \label{trialsln}
    \psi(r) = f_b + g(r)N(r,E;\theta). 
\end{equation}
The network $N(r, E; \theta)$ is a single-output neural network (NN) with two hidden layers and two inputs: $r$ and $E$. The total trainable parameters of the network are represented by $\theta$.  Figure~\ref{fig:us_ann} illustrates how both the eigenvalue $E$ and the trial solution $\psi(r)$ are determined; the eigenvalue $E$ is obtained by passing a constant vector of ones (having the same number of points as $r$) through a single linear neuron, whose weights are subsequently updated via optimisation. Thereafter, both $r$ and $E$ are fed into a feed-forward NN to produce the output $N(r, E; \theta)$. Finally, the parametric trick in~(\ref{trialsln}) is applied to estimate the wave function. The parametric function $g(r)$ plays an important role in the 
performance of the NN.  As recommended in~\cite{HMM}, it should be smooth and bounded; therefore, we adopt the form given in~\cite{usnn},
\begin{equation} 
   \label{gfunction}
   g(r) = \left(1-e^{-(r-r_{min})} \right)\left(1-e^{-(r-r_{max})} \right).
\end{equation}
Here, $r$ is the domain, and $r_{min}$ and $r_{max}$ are the end points of the domain. It can be easily verified that the trial solution~(\ref{trialsln}) satisfies the boundary conditions~(\ref{BoundaryConditions}). The bound state energies and the corresponding wave functions are determined by minimising the following composite loss function 
\begin{equation}
    \label{compositeLossFunction}
    \mathcal{L}(\theta) = \nu_{DE}\mathcal{L}_{DE} + \nu_{f}\mathcal{L}_{f}+\nu_{\lambda}\mathcal{L}_{\lambda}+\nu_{drive}\mathcal{L}_{drive},
\end{equation}
where 
\begin{equation}
    \label{compositeLossFunction_LDE}
    \mathcal{L}_{DE} = \left<\left(\left[ -\frac{\hbar^2}{2\mu}\frac{d^2}{dr^2}+V(r)\right]\psi(r) - E\psi(r) \right)^2 \right>_r  
\end{equation}
% \frac{1}{N}\sum_{i = 1}^{N}

\begin{equation}
    \label{compositeLossFunction_Lreg} 
    \mathcal{L}_{f} = \frac{1}{\left<f(r,E)^2\right>_r}, \quad  \mathcal{L}_{E} = \frac{1}{\left<E^2\right>} \quad \mbox{ and } \quad \mathcal{L}_{drive} = \left<e^{-E + c}\right>.
\end{equation}
% \subsection{Explanation of the loss terms}
$\left<\cdot \right>_r$ represents taking an average over the variable $r$ and $\left<\cdot \right>$ represents the simple sample mean. The $\mathcal{L}_{DE}$ loss ensures that the estimated eigen-pairs satisfy the Schr\"{o}dinger equation while $\mathcal{L}_{f}$ and $\mathcal{L}_{E}$ encourage the NN to learn non-trivial solutions of the wave-function solution $\psi(r)$ and the energy $E$. Finally, $\mathcal{L}_{drive}$ drives the NN to scan for larger eigenvalues. The derivatives with respect to the inputs and training weights are determined using automatic differentiation~\cite{autodiff}.
% \subsection{Training}
In the training stage,  the $\mathcal{L}_{drive}$ loss term is gradually increased by incrementing the value of $c$ in regular intervals. As mentioned earlier, this will encourage the NN to search for higher eigenpairs. At every epoch a batch of $r$ points in the interval $\left[r_{min},r_{max}\right]$ is perturbed by adding a Gaussian noise. This is done to prevent the NN from learning the solutions only at fixed intervals.

In determining the bound states eigen-pairs for the molecular potential~(\ref{toyPot}), we used a NN with 50 nodes per hidden layer and the sin activation function~\cite{usnn}. The domain was divided into $100$ equally spaced collocation points, and $10^5$ training epochs were used. We found, empirically, that setting $\nu_{DE} = \nu_{f} = \nu_{\lambda} = \nu_{drive} = 1$ in equation~(\ref{compositeLossFunction}) produced optimal results. The hyperparameters used for the Woods-Saxon potential~(\ref{WoodSaxPontential}) were the following: 380 collocation points, 300 units per hidden layer, $250\times 10^3$ training epochs. Here we set $\nu_{DE} = 10^{-7}, \nu_{f} = 10^{12}, \nu_{\lambda} = 10^4 \mbox { and }\nu_{drive} = 1$, this will ensure that all the individual loss terms~(\ref{compositeLossFunction_LDE} - \ref{compositeLossFunction_Lreg}) contribute equally to the total loss.

For both potentials, we added a Gaussian noise of $\mathcal{N}\left(\mu = 0, \sigma^2 = 0.06r_{max}\right)$ to the batch input data for every epoch and the loss function was minimized using the Adam optimizer~\cite{adamoptimizer} with the learning rate of $8\times10^3$.
% PART 2

\begin{figure*}[t]
\centering
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{saip_molecular_potential.png} 
        \caption{\label{fig:molecular_pot}Molecular potential}
    \end{subfigure}  
\hfil 
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{saip_woodsSaxon_potential.png} 
        \caption{\label{fig:ws_pot}Woods-Saxon potential}
    \end{subfigure}
\caption{\label{fig:potentials} The toy molecular potential and the Woods-Saxon potential used as input for solving the Schrödinger equation and determining the energy levels and wavefunctions of the system.}  
\end{figure*} 

\begin{figure}
\centering
\resizebox{0.65\textwidth}{!}{%
\begin{tikzpicture}[->,>=stealth', shorten >=1pt, auto, node distance=1.5cm, 
  thick, main node/.style={circle,draw,font=\sffamily\small\bfseries}] 

% Inputs

\node[main node] (x) {$r$};
\node[main node, below=1cm of x] (lambda) {$E$};
\node[main node, left = 1cm of lambda] (input1) {1};

% Hidden Layer 1
\node[main node, right=2cm of x, yshift=1cm] (h11) {$\sigma$};
\node[main node, below=0.8cm of h11] (h12) {$\sigma$};
\node[main node, below=0.8cm of h12] (h13) {$\sigma$};
\node[below=1.2cm  of h12] (xhh3) {\vdots};
\node[main node, below=0.8cm of xhh3] (h14) {$\sigma$};

% Hidden Layer 2
\node[main node, right=1.8cm of h11] (h21) {$\sigma$};
\node[main node, below=0.8cm of h21] (h22) {$\sigma$};
\node[main node, below=0.8cm of h22] (h23) {$\sigma$};
\node[below=1.2cm  of h22] (xhh4) {\vdots};
\node[main node, below=0.8cm of xhh4] (h24) {$\sigma$};

% Output
\node[main node, right=1.5cm of h22, yshift=-2em] (output) {};
\node[right=0.05mm of output] (test) {$\mathcal{N}(r,E)$};
% Parameterization
\node[draw, rectangle, minimum width=2.8cm, minimum height=0.8cm, right=2cm of output] (param) {$f_b + g(r)\mathcal{N}(r,E) $};
\node[right=0.5cm of param] (finalout) {$\psi(r,E)$};

% Connections
\draw (input1) -- (lambda);

\foreach \i in {x,lambda}{
    \foreach \j in {h11,h12,h13,h14}{
        \draw (\i) -- (\j);
    }    
}
  
\foreach \i in {h11,h12,h13,h14}{
    \foreach \j in {h21,h22,h23,h24}{
        \draw (\i) -- (\j);
    }
}


\foreach \i in {h21,h22,h23,h24}{
    \draw (\i) -- (output);
}
 
\draw (test) -- (param);
\draw (param) -- (finalout);

% Layer Labels
\node[above=1.1cm of x] {Inputs};
\node[above=0.2cm of h11,align=center] {Hidden\\Layer 1};
\node[above=0.2cm of h21,align=center] {Hidden\\Layer 2};
\node[above=2.5cm of output,align=center] {Output};
\node[above=2.3cm of param] {Parameterization};
\end{tikzpicture}    
}
\caption{The network architecture of the unsupervised NN. $\sigma$ represents non-linear activation function  }
\label{fig:us_ann}
\end{figure}

\section{Results}
\begin{figure*}[h!]
\centering
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{EigenValues.jpg}
        \caption{\label{fig:EVPotetential1}Molecular potential}
    \end{subfigure}  
\hfil 
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{WS_EigenNN.jpg}
        \caption{\label{fig:EVPotetential2}Woods-Saxon potential}
    \end{subfigure}
\caption{\label{fig:predEigenvalues}Plot of the eigenvalues' search history during training. The dashed lines are the exact eigenvalues generated by the two potentials.}    
\end{figure*} 
\begin{table}[h!]
 \begin{center}
 \begin{tabular}{@{}*4c@{}}          
\toprule[1.5pt] 
\multicolumn{1}{c}{$n$} & \multicolumn{1}{c}{$E$ (MeV)}& \multicolumn{1}{c}{Relative error(\%)}&\multicolumn{1}{c}{Ref}\\             
\bottomrule[1.5pt]    
$0$&
$ -4.571 182 814$&
  -&
 \cite{1997Paper} \\  
 &
 $-4.570504665$&
 $1.482\times 10^{-2}$&
 This work\\    
\newline 
 $1$&
$ -0.884280776$&
  -&
\cite{1997Paper}\\  
 \newline
 &
 $ -0.884224414$&
 $6.37\times 10^{-3}$&
 This work\\  
 \bottomrule[1.5pt]    
\end{tabular}
\end{center}
 \caption{The eigenvalues for the S-wave bound states for the toy potential~(\ref{toyPot})}    
 \label{table:EigenvaluesV1}    
\end{table}


\begin{table}[h!]
 \begin{center}
 \begin{tabular}{@{}*4c@{}}          
\toprule[1.5pt] 
\multicolumn{1}{c}{$n$} & \multicolumn{1}{c}{$E$ (MeV)}& \multicolumn{1}{c}{Relative error(\%)}&\multicolumn{1}{c}{Ref}\\             
\bottomrule[1.5pt]    
$0$&
$ -27.00085$&
  -&
 \cite{Mmusi} \\  
 & 
 $-26.920$&
 $2.99\times 10^{-1}$&
 This work\\    
\newline 
 $1$&
$ -16.42311$&
  -&
\cite{Mmusi}\\  
 \newline
 &
 $ -16.384$&
 $2.38\times 10^{-1}$&
 This work\\
\newline 
 $2$&
$ -3.40020$&
  -&
\cite{Mmusi}\\  
 \newline
 &
 $ -2.895$&
 $14.85$& 
 This work\\   
 \bottomrule[1.5pt]    
\end{tabular}
\end{center}
 \caption{The eigenvalues for the S-wave bound states for the potential~(\ref{WoodSaxPontential}).}    
 \label{table:EigenvaluesV2}    
\end{table}  



As shown in Figure~\ref{fig:predEigenvalues}, the NN initially converges to the lowest energy level (the ground state) before moving on to excited states. This behaviour occurs regardless of the starting point. It aligns with a well-known property of NNs and PINNs: they tend to learn low-frequency components of a solution before capturing higher-frequency ones~\cite{expertGuide}. Tables~\ref{table:EigenvaluesV1}–\ref{table:EigenvaluesV2} list the predicted eigenvalues and their corresponding absolute relative errors. Figures~\ref{fig:wavefunv1}–\ref{fig:wavefunv2} compare the predicted and theoretical wavefunctions for the two potentials. Based on the width of the $E_0$ plateau in Figure~\ref{fig:EVPotetential1}, we expect $E_0$ to have the smallest relative error. However, the results show that the network finds it easier to `learn' the first excited state. One possible explanation is that the driving term $\mathcal{L}_{drive}$ pushed the network toward higher energy levels before the ground state eigenvalue was `fully learned'.

\begin{figure*}[h!]
\centering
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{GroundState_WaveFunctions_final.jpg}
        \caption{\label{fig:GS_wave}The Ground state wave function}
    \end{subfigure}  
\hfil 
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{FirstState_WaveFunctions_final.jpg}
        \caption{\label{fig:GS_wave}The first excited state wave function}
    \end{subfigure}
\caption{\label{fig:wavefunv1}The two S-wave-functions for the toy potential~(\ref{toyPot}). The red curve is the true eigenfunction while the blue curve represents the predicted wave-function.}    
\end{figure*} 



\begin{figure*}[h]
\centering
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{WSGSWs.jpg}
        \caption{\label{fig:WSGS_wave}The Ground state wave function}
    \end{subfigure}  
\hfil 
    \begin{subfigure}[t]{0.49\textwidth}
        \centering 
        \includegraphics[width=\textwidth]{WSfirstSWs.jpg}
        \caption{\label{fig:WS2S_wave}The first excited state wave function}
    \end{subfigure}
\hfil 
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{WSsecondSWs.jpg}
        \caption{\label{fig:GS_wave1}The second excited state wave function}
    \end{subfigure}    
\caption{\label{fig:wavefunv2}The three S-wave-functions for the WoodsSaxon potential~(\ref{WoodSaxPontential}). The blue curve represents the true eigenfunction while blue is the predicted wave-function.} 
\end{figure*}

\section{Conclusion}
In this work, we applied the unsupervised PINN approach to solve for the eigenpairs of the two-body radial Schrödinger equation. The results indicate that PINNs have the potential, if properly tuned, to serve as a viable alternative to traditional numerical methods for solving bound-state quantum eigenvalue problems. The loss term weights and hyperparameters were adjusted manually. In our work, we plan to implement adaptive loss weighting using the method proposed in~\cite{bragga}, and to employ the Hyperopt package~\cite{hyperopt} for automated hyperparameter optimisation. Furthermore, since the Hamiltonian for bound states is Hermitian this implies that the wave functions are orthonormal and normalisable, we will add these properties into the loss function to make it more robust. By incorporating these approaches into the training process, it has the potential to reduce the occurrence of spurious bound-state solutions and improve prediction accuracy across different potential models.

\section*{Acknowledgments}
JTT gratefully acknowledges funding from the Absa Research Grant (Sol Plaatje 
University) for providing the equipment essential to this research. ASC is partially 
supported by the National Research Foundation of South Africa.

%Bibliography
%make sure you fill the .bib properly, e.g. edit "biblio.bib" and add all your references in there.
\bibliographystyle{IEEEtran}
\bibliography{biblio}

\end{document}