\section{Introduction}

\noindent It is well-known that the exact analytical solution of the
hyperradial Schr\"{o}dinger equation, in any arbitrary spatial dimension $%
(D\geq 2),$ for its bound state energy levels, is fundamental to
understanding of molecular spectrum of a diatomic molecules in a
nonrelativistic quantum mechanics, since the resulting wave function
contains all the necessary information to describe a quantum system fully.
There are only a few potentials for which the radial Schr\"{o}dinger
equation can be solved explicitly for all $n$ and all $l$ quantum numbers$.$
So far, some useful methods have been developed, such as supersymmetry
(SUSY) [1], the Nikiforov-Uvarov method (NU) [2], the Pekris approximation
[3], asymptotic iteration method (AIM) [4], variational [5], the hypervirial
perturbation [6], the shifted $1/N$ expansion (SE) and the modified shifted $%
1/N$ expansion (MSE) [7], exact quantization rule (EQR) [8], perturbative
formalism [9], polynomial solution [10] and wave function ansatz method [11]
to solve the radial Schr\"{o}dinger equation exactly or quasi-exactly for $%
l\neq 0$ within a given potential.

The kratzer and modified Kratzer-type potentials [12] that we study in the
present work are important molecular potentials describe the interaction
between two atoms and have attracted a great of interest for some decades in
the history of quantum chemistry. Until now, these potentials have been used
extensively to describe the molecular structure and interactions [13]. In
this work, we attempt to apply the recently proposed EQR method [8], as an
alternative approach, to obtain analytically ($l\neq 0$) exact energy levels
and the corresponding wave functions of the $D$-dimensional $(D\geq 2)$
Schr\"{o}dinger equation with the Kratzer-type potentials. Further, we
calculate numerically these bound state energy levels for $LiH,$ $CH,$ $HCl,$
$CO,$ $NO,$ $O_{2},$ $N_{2}$ and $I_{2}$ diatomic molecules for different
values of $n$ and $l$ quantum numbers. Further, we give numerical tests
using energy calculations for interdimensional degeneracy, i.e., $%
(n,l,D)\rightarrow (n,l\pm 1,D\mp 2)$ corresponding to the confined $D=2-4$
diimensional Kratzer's potential.

This paper is organized as follows. In Section \ref{EQR} the EQR method is
reviewed. In Section \ref{SOTK}, the $D$-dimensional ($D\geq 2$)
Schr\"{o}dinger equation is solved by this method with $l\neq 0$ quantum
numbers to obtain the energy eigenvalues and the corresponding hyperradial
wave functions of the Kratzer-type potential. In Section \ref{ID} the
interdimensional degeneracy is introduced. Numerical calculations of the
energy levels and test of interdimensional degeneracy for many typical
diatomic molecules $CO,$ $NO,$ $CH,$ $HCl,$ $LiH,$ $N_{2},$ $O_{2},$ and $%
I_{2}$ are performed for various values of $n$ and $l$ quantum numbers to
compare with those obtained by other methods in Section \ref{NA}. Finally,
some concluding remarks are given in Section \ref{CR}.

\section{Exact Quantization Rule Method}

\label{EQR}A brief outline to the EQR method is presented with an extension
to the $D$-dimensional space ($D\geq 2$). The details can be found in Ref.
[8]. The EQR has recently been proposed to solve the one-dimensional ($1D$)
Schr\"{o}dinger equation:

\begin{equation}
\frac{d^{2}}{dx^{2}}\psi (x)+k(x)^{2}\psi (x)=0,\text{ \ }k(x)=\frac{1}{%
\hbar }\sqrt{2\mu \left[ E-V(x)\right] },
\end{equation}
where $\mu $ represents the reduced mass of the two interacting particles, $%
k(x)$ is the momentum and the potential $V(x)$ is a piecewise continuous
real function of $x.$ Equation (1) can be reduced into the Riccati equation:

\begin{equation}
\frac{d}{dx}\phi (x)+\phi (x)^{2}+k(x)^{2}=0,
\end{equation}
where $\phi (x)=\psi (x)^{-1}d\psi (x)/dx$ is the logarithmic derivative of
wave function $\psi (x).$ Due to the Sturm-Liouville theorem, the $\phi (x)$
decreases monotonically with respect to $x$ between two turning points,
where $E\geq V(x).$ Specifically, as $x$ increases across a node of the wave
function $\psi (x),$ where $E\geq V(x),$ $\phi (x)$ decreases to $-\infty ,$
jumps to $+\infty ,$ and then decreases again. By carefully studying the $1D$
Schr\"{o}dinger equation, Ma and Xu [8] obtained the EQR as

\begin{equation}
\int\limits_{x_{a}}^{x_{b}}k(x)dx=N\pi +\int\limits_{x_{a}}^{x_{b}}\phi (x)
\left[ \frac{dk(x)}{dx}\right] \left[ \frac{d\phi (x)}{dx}\right] ^{-1}dx,
\end{equation}
where $x_{a}$ and $x_{b}$ are two turning points determined by letting $%
E=V(x),$ $N$ is the number of nodes of \ $\phi (x)$ in the region $E\geq
V(x) $ and it is larger by one than the number of nodes of wave function $%
\psi (x),$ i.e. $N=n+1.$ Further, Ma and Xu [8] have generalized the EQR to $%
3D$ radial Schr\"{o}dinger equation with a spherically symmetric potential.
Here, we extend the $3D$-generalization of the EQR given by Ref. [8] into
the $D$-dimensional ($D\geq 2$) Schr\"{o}dinger equation. It is well-known
that in the $D$-dimensional hyperradial Schr\"{o}dinger equation, for any
arbitrary spherically symmetric potential, can be simply written as [14]
\[
\left[ \frac{d^{2}}{dr^{2}}+\frac{\left( D-1\right) }{r}\frac{d}{dr}-\frac{%
l(l+D-2)}{r^{2}}+\frac{2\mu }{\hbar ^{2}}\left( E-V(r)\right) \right] \psi
_{n,l,m}(r,\Omega _{D})=0,
\]
\begin{equation}
\psi _{n,l,m}(r,\Omega _{D})=\psi _{n,l}(r)Y_{l}^{m}(\Omega _{D}),
\end{equation}
where $l$ denotes the orbital angular momentum quantum number, $\mu $ is the
reduced mass, $\psi _{n,l}(r)$ is the hyperradial part of the wave function
and $Y_{l}^{m}(\Omega _{D})$ is the hyperspherical harmonics. Furthermore,
to remove the first derivative from the above equation, we define a new
radial wave function $R(r)$ by means of equation
\begin{equation}
\psi _{n,l}(r)=r^{-(D-1)/2}R(r),
\end{equation}
which will give the radial wave function $R(r)$ satisfying the wave equation
\begin{equation}
\left\{ \frac{d^{2}}{dr^{2}}-\frac{\left[ 4l\left( l+D-2\right) +\left(
D-1\right) \left( D-3\right) \right] }{4r^{2}}+\frac{2\mu }{\hbar ^{2}}%
\left( E-V(r)\right) \right\} R(r)=0.
\end{equation}
Further, the hyperradial equation (6) can be written into a simple analogy
of the $2D$ and $3D$ radial Schr\"{o}dinger equation given in Ref. [8] as
[15]
\begin{equation}
\frac{d^{2}}{dr^{2}}R(r)=-\frac{2\mu }{\hbar ^{2}}\left[ E-V_{\text{eff}}(r)%
\right] R(r),\text{ \ }V_{\text{eff}}(r)=\frac{\left( \eta ^{2}-1/4\right)
\hbar ^{2}}{2\mu r^{2}}+V(r),
\end{equation}
with the parameter
\begin{equation}
\eta =l+\frac{1}{2}\left( D-2\right) ,
\end{equation}
which is a linear combination of the spatial dimensions $D$ and the angular
momentum quantum number $l.$ Thus, the EQR can be generalized to the $D$%
-dimensional ($D\geq 2$) hyperradial Schr\"{o}dinger equation, with the new
effective potential form given by (7) as

\begin{equation}
\int\limits_{r_{a}}^{r_{b}}K(r)dr=N\pi +Q_{c},\text{ }K(r)=\frac{2\mu }{%
\hbar ^{2}}\left[ E-V_{\text{eff}}(r)\right] ,
\end{equation}
where
\begin{equation}
Q_{c}=\int\limits_{r_{a}}^{r_{b}}\phi (r)\left[ \frac{dk(r)}{dr}\right] %
\left[ \frac{d\phi (r)}{dr}\right] ^{-1}dr,
\end{equation}
is called the quantum correction, $N\pi $ is the contribution from the nodes
of wave function and $K(r)$ being the momentum between the two turning
points $r_{a}$ and $r_{b}.$ It has also been found that this quantum
correction is independent of the number of nodes of wave function for
solvable quantum system. Thus, the quantum correction in Eq. (10) can be
replaced by the ground-state calculations:
\begin{equation}
Q_{c}=\int\limits_{r_{a}}^{r_{b}}\phi _{0}(r)\left[ \frac{dk_{0}(r)}{dr}%
\right] \left[ \frac{d\phi _{0}(r)}{dr}\right] ^{-1}dr=\pi q.
\end{equation}

\section{Solutions of the Kratzer and Modified Kratzer-Type Potentials}

\label{SOTK}Considerable interest has recently been shown in the
Kratzer-Fues [12,16-19] as a model to describe internuclear vibration of a
diatomic molecule [20-23]. This potential can be expressed in the form
\begin{equation}
V(r)=-D_{e}\left( \frac{2r_{e}}{r}-\frac{r_{e}^{2}}{r^{2}}\right) ,
\end{equation}
where $D_{e}$ is the dissociation energy between two atoms in a solid and $%
r_{e}$ is the equilibrium internuclear length. This type of the potential in
(12) has been studied through the smooth transformation [16], the algebraic
approach [17], and the AIM [18]. The standard Kratzer's potential has been
recently modified by adding a $D_{e}$ term to the potential [24]. A new type
of this potential is called the modified Kratzer potential:

\begin{equation}
V(r)=D_{e}\left( \frac{r-r_{e}}{r}\right) ^{2},
\end{equation}
which is shifted in amount of $D_{e}.$ The above potential has been studied
in $D$-dimensions through the polynomial solutions [25] and also the
selection of a suitable wave function ansatz [15]. The modified Kratzer
potential has also been studied before in the $3D$ through the NU method
[24]. This potential can be simply taken as
\begin{equation}
V(r)=\frac{A}{r^{2}}-\frac{B}{r}+C,
\end{equation}
where $A=D_{e}r_{e}^{2},$ $B=2D_{e}r_{e}$ and $C=D_{e}$ [15,25,26]$.$ For
this spherically symmetric modified Kratzer potential, the effective
potential in Eq. (7) can be arranged to

\begin{equation}
V_{\text{eff}}(r)=\frac{2\mu A+\frac{1}{4}\left[ (M-2)^{2}-1\right] \hbar
^{2}}{2\mu r^{2}}-\frac{B}{r}+C,\text{ \ }M=D+2l.
\end{equation}
Further, by introducing the notation

\begin{equation}
\Lambda (\Lambda +1)\hbar ^{2}=2\mu A+\frac{1}{4}\left[ (M-2)^{2}-1\right]
\hbar ^{2},\text{ \ }
\end{equation}
from which we obtain a physically acceptable solution

\begin{equation}
\Lambda =\frac{-1+\sqrt{(M-2)^{2}+8\mu A/\hbar ^{2}}}{2},
\end{equation}
and finally the effective potential becomes

\begin{equation}
V_{\text{eff}}(r)=\frac{\Lambda (\Lambda +1)\hbar ^{2}}{2\mu r^{2}}-\frac{B}{%
r}+C.
\end{equation}
On the other hand, by taking $\phi (r)=R(r)^{-1}dR(r)/dr,$ the Riccati
equation reads

\begin{equation}
\frac{d}{dr}\phi (r)+\phi (r)^{2}=-\frac{2\mu }{\hbar ^{2}}\left[ E-C-\frac{%
\Lambda (\Lambda +1)\hbar ^{2}}{2\mu r^{2}}+\frac{B}{r}\right] .
\end{equation}
Now, we apply the EQR [8,26,27] given by (9) to solve the above Riccati
equation. At first, we must determine the turning points $r_{a}$ and $r_{b}$
from solving the relation $E=V_{\text{eff}}(r)$ with $r_{a}<$ $r_{b}.$ Thus,
we find

\begin{equation}
r_{a}=\frac{1}{-2(E-C)}\left[ B-\sqrt{B^{2}+\frac{2\Lambda (\Lambda
+1)(E-C)\hbar ^{2}}{\mu }}\right] ,
\end{equation}

\begin{equation}
r_{b}=\frac{1}{-2(E-C)}\left[ B+\sqrt{B^{2}+\frac{2\Lambda (\Lambda
+1)(E-C)\hbar ^{2}}{\mu }}\right] .
\end{equation}
The momentum between two turning points $r_{a}$ and $r_{b}$ can be
calculated by $k(r)=\sqrt{2\mu \left[ E_{0}-V_{\text{eff}}(r)\right] }/\hbar
$ as

\begin{equation}
k(r)=\frac{1}{\hbar r}\sqrt{-2\mu (E-C)(r-r_{a})(r_{b}-r)}.
\end{equation}
Since the phase angle $\phi (r)$ decreases as $r$ increases, then the
solution with one node and no pole only has a unique solution $\phi
_{0}(r)=\alpha _{1}r^{-1}+\alpha _{2}$ with $\alpha _{1}>0$ due to the
monotonic property$.$ Substituting this solution into the Riccati equation
(19), we find the ground-state solution as

\begin{equation}
\phi _{0}(r)=\frac{\Lambda +1}{r}-\frac{\mu B}{(\Lambda +1)\hbar ^{2}},\text{
with }E_{0}=C-\frac{\mu B^{2}}{2(\Lambda +1)^{2}\hbar ^{2}}.
\end{equation}
Evidently $\phi _{0}(r)$ is negative when $r\rightarrow \infty $, so that
the solution satisfies the physically admissible boundary condition. In
order to obtain the energy levels of all the bound-states, making use of the
integrals given in Appendix A, we are able to calculate the quantum
correction from Eq. (11) as

\[
Q_{c}=\frac{\sqrt{-2\mu (E_{0}-C)}}{2(\Lambda +1)^{2}\hbar ^{3}}%
\int\limits_{r_{a}}^{r_{b}}\frac{\left[ (\Lambda +1)^{2}\hbar ^{2}-\mu Br%
\right] }{r\sqrt{(r-r_{a})(r_{b}-r)}}\left[ (r_{a}+r_{b})r-2r_{a}r_{b}\right]
dr
\]

\begin{equation}
=\left[ \Lambda -\sqrt{\Lambda (\Lambda +1)}\right] \pi ,
\end{equation}
where the following formulas

\begin{equation}
r_{a}+r_{b}=-\frac{B}{(E_{0}-C)},\text{ \ }r_{a}r_{b}=-\frac{\Lambda
(\Lambda +1)\hbar ^{2}}{2\mu (E_{0}-C)},\text{ }\sqrt{-2\mu (E_{0}-C)}=\frac{%
\mu B}{(\Lambda +1)\hbar },
\end{equation}
have been used in finding Eq. (24). On the other hand, for any $n$ state, we
calculate the left hand side of Eq. (9) for $V_{\text{eff}}(r)$ given in Eq.
(18) using Eq. (22) as
\begin{equation}
\int\limits_{r_{a}}^{r_{b}}K(r)dr=\left[ \frac{\mu B}{\sqrt{-2\mu (E_{n}-C)}%
\hbar }-\sqrt{\Lambda (\Lambda +1)}\right] \pi .
\end{equation}
Further, the substitution of Eqs (24) and (26) into Eq. (9), enables us to
obtain the energy levels of all bound states in any arbitrary dimension $%
D\geq 2$ as

\[
E_{nl}^{(D)}=C-\left( \frac{2\mu }{\hbar ^{2}}B^{2}\right) \left[ 2n+1+\sqrt{%
(D+2l-2)^{2}+\frac{8\mu }{\hbar ^{2}}A}\right] ^{-2}
\]
\begin{equation}
=C-\frac{\mu }{2\hbar ^{2}}\left( \frac{B}{\widetilde{n}}\right) ^{2},\text{
}\widetilde{n}=N+\Lambda ,\text{ \ }n,l=0,1,2,\cdots
\end{equation}
where $N=n+1$ was used with $n$ is the number of nodes in the wave functions
$R(r).$

Therefore, for the $D$-dimensional Kratzer's potential ($C=0$), the energy
spectra become

\begin{equation}
E_{nl}^{(D)}=-\left( \frac{2\mu }{\hbar ^{2}}D_{e}^{2}r_{e}^{2}\right) \left[
n+\frac{1}{2}+\sqrt{\left( l+\frac{D}{2}-1\right) ^{2}+\frac{2\mu }{\hbar
^{2}}D_{e}r_{e}^{2}}\right] ^{-2},
\end{equation}
and also for the modified Kratzer-type ($C=D_{e}$) potential read

\begin{equation}
E_{nl}^{(D)}=D_{e}-\left( \frac{2\mu }{\hbar ^{2}}D_{e}^{2}r_{e}^{2}\right) %
\left[ n+\frac{1}{2}+\sqrt{\left( l+\frac{D}{2}-1\right) ^{2}+\frac{2\mu }{%
\hbar ^{2}}D_{e}r_{e}^{2}}\right] ^{-2},\text{ }E_{nl}^{(D)}<D_{e}.
\end{equation}
In what follows, we want to obtain the corresponding normalized hyperwave
function. At first, we define

\begin{equation}
\nu (\nu +D-2)=l(l+D-2)+\frac{2\mu A}{\hbar ^{2}},
\end{equation}

\begin{equation}
\rho =2\kappa r,\text{ }\tau =\frac{\mu B}{\hbar ^{2}\kappa },\text{ }\kappa
=\sqrt{-2\mu (E-C)/\hbar ^{2}},
\end{equation}
where for bound state case, states with $E<C$, we have real numbers for the
parameters $\kappa $ and $\tau .$ Thus, the hyperradial equation (4)
transforms to a standard form given by

\begin{equation}
\psi _{n,l}^{\prime \prime }(\rho )+\frac{\left( D-1\right) }{\rho }\psi
_{n,l}^{\prime }(\rho )+\left[ -\frac{1}{4}+\frac{\tau }{\rho }-\frac{\nu
(\nu +D-2)}{\rho ^{2}}\right] \psi _{n,l}(\rho )=0,
\end{equation}

\begin{equation}
\nu =-\frac{(D-2)}{2}+\frac{1}{2}\sqrt{\left( 2l+D-2\right) ^{2}+\frac{8\mu A%
}{\hbar ^{2}}},
\end{equation}
whose solutions are given as [19,25]

\begin{equation}
\psi _{n,l}(\rho )=N(n,l)\rho ^{\nu }e^{-\rho /2}{}_{1}F_{1}(-\tau +\nu +%
\frac{D-1}{2},2\nu +D-1;\rho )
\end{equation}
from which we can obtain the quantum condition as $-\tau +\nu +\frac{D-1}{2}%
=-n,$ $n=0,1,2,\cdots $ $[28,29].$ Furthermore, we can renormalize the bound
state solutions in (34) since these solutions are finite and convergent for
both $\rho =0$ $(r=0)$ and $\rho \rightarrow \infty .$

On the other hand, using the following relationship between the associated
Laguerre function and the hypergeometric function:

\begin{equation}
L_{n}^{\alpha }(z)=\frac{\Gamma (\alpha +n+1)}{n!\Gamma (a+1)}%
{}_{1}F_{1}(-n,\alpha +1;z),
\end{equation}
we may further rewrite the wave functions as

\begin{equation}
\psi _{n,l}(\rho )={\cal N}(n,l)\rho ^{\nu }e^{-\rho /2}{}L_{n}^{(2\nu
+D-2)}(\rho ),
\end{equation}
where ${\cal N}(n,l)$ is the normalizing factor to be determined by the
renormalization condition

\begin{equation}
\int\limits_{0}^{\infty }\left| \psi _{n,l}(\rho )\right| ^{2}r^{(D-1)}dr=1.
\end{equation}
To find the normalizing factor, let us review the generalized Coulomb-like
integral $J_{n,\alpha }^{(\beta )}(z)$ for noninteger $\alpha $ and $\beta $
derived in [30] as

\[
J_{n,\alpha }^{(\beta )}(z)=\int\limits_{0}^{\infty }e^{-z}z^{\alpha +\beta
} \left[ L_{n}^{(\alpha )}(z)\right] ^{2}dz=\frac{\Gamma (\alpha +n+1)}{%
\Gamma (n+1)}\sum\limits_{j=0}^{n}(-1)^{j}\frac{\Gamma (n-j-\beta )}{\Gamma
(-j-\beta )}
\]

\begin{equation}
\times \frac{\Gamma (\alpha +\beta +j+1)}{\Gamma (\alpha +j+1)}\frac{1}{%
\Gamma (j+1)\Gamma (n-j+1)},\text{ }{\cal R}\text{e}(\alpha +\beta +1)>0.
\end{equation}
In the present case, we have $\beta =1.$ This $J_{n,\alpha }^{(1)}$ has only
two nonzero contributions (for $j=n-1$ and $n$) in the sum (38) because of
the gamma functions of negative integers. The result $J_{n,\alpha }^{(1)}=%
\frac{(2n+\alpha +1)\Gamma (\alpha +n+1)}{\Gamma (n+1)}$ enables us to write
for our case

\begin{equation}
J_{n,2\nu +D-2}^{(1)}=\frac{(2\nu +2n+D-1)(2\nu +n+D-2)!}{n!},
\end{equation}
which gives the normalization factor

\[
{\cal N}(n,l)=\left[ \frac{n!(2\kappa )^{D}}{(2\nu +2n+D-1)(2\nu +n+D-2)!}%
\right] ^{1/2},
\]

\begin{equation}
\kappa =\frac{\mu B}{\hbar ^{2}(n+\Lambda +1)},\text{ }\Lambda =\nu +\frac{%
(D-3)}{2}.
\end{equation}
Therefore, we can finally obtain the re-normalized hyper wave function for
the modified Kratzer potential as

\begin{equation}
\psi _{n,l,m}(r,\Omega _{D})=\left[ \frac{n!}{(\nu _{1}+2n+1)(\nu _{1}+n)!}%
\right] ^{1/2}(2\kappa )^{\nu _{2}}r^{\nu }e^{-\kappa r}{}L_{n}^{(\nu
_{1})}(2\kappa r)Y_{l}^{m}(\Omega _{D}),
\end{equation}
with

\begin{equation}
2\kappa =\frac{8\mu D_{e}r_{e}}{\hbar ^{2}(2n+\nu _{1}+1)},\text{ }\nu _{1}=%
\sqrt{\left( 2l+D-2\right) ^{2}+\frac{8\mu }{\hbar ^{2}}D_{e}r_{e}}\text{ ,}
\end{equation}
\begin{equation}
\nu _{2}=\frac{\nu _{1}+2}{2},\text{ }\nu =\frac{\nu _{1}-(D-2)}{2}.
\end{equation}
which is found to be consistent with our previous findings in Ref. [25]
using the polynomial solution.

\section{Interdimensional Degeneracy}

\label{ID}From Eq. (27), it can be seen that two interdimensional states are
degenerate whenever [31]

\begin{equation}
(n,l,D)\rightarrow (n,l\pm 1,D\mp 2).
\end{equation}
Thus, a knowledge of $E_{nl}^{(D)}$ for $D=2$ and $D=3$ provides the
information necessary to find $E_{nl}^{(D)}$ for other higher dimensions.

For example, $E_{0,4}^{(2)}=E_{0,3}^{(4)}=E_{0,2}^{(6)}=E_{0,1}^{(8)}.$ This
is the same transformational invariance described for bound states of free
atoms and molecules [32,33] and demonstrates the existence of
interdimensional degeneracies among states of the confined Kratzer
potential. \

\section{Numerical Applications}

\label{NA}In this work, we have given numerical calculations for some
typical diatomic molecules to compare with those obtained by other numerical
exact methods given by the AIM [18] and the NU method [24]. In this regard,
we have calculated the bound state energy levels with Eq. (28) using the
Kratzer's potential for $CO,$ $NO,$ $HCl,$ $LiH,$ $O_{2}$ and $I_{2}$
diatomic molecules for selected $n,$ $l$ quantum numbers with parameter
values given in Tables 1. These parameters values are taken from Refs.
[34-37]. Therefore, we give the explicit values of Kratzer's energy levels
for different values of $n$ and $l$ in Tables 2-3. We further compare these
results with those obtained by the AIM [18]. In tables 4-5, we have also
tested the interdimensional degeneracy ($D=2-4$) using the energy
calculations for $CO,$ $NO,$ $HCl,$ $LiH,$ $O_{2},$ and $I_{2}$ diatomic
molecules. On the other hand, for the recently proposed modified Kratzer
potential, we have used the parameters for $CO,$ $NO,$ $CH$ and $N_{2}$
diatomic molecules given in Table 6 to compare with the NU method [24]. For
this case, in Tables 7-8, we give the energy levels for $CO,$ $NO,$ $CH$ and
$N_{2}$ diatomic molecules for various $n,$ $l$ quantum numbers by using Eq.
(29).

\section{Concluding Remarks}

\label{CR}In this study, we have used a different approach to find the $%
l\neq 0$ analytical and numerical solutions to the $D$-dimensional
hyperradial Schr\"{o}dinger equation with Kratzer-type potentials for
various diatomic molecules within the framework of the exact quantization
rule. For such potentials, the problem is simply reduced to a Coulomb
potential plus an inverse quadratic potential term. The exact eigensolutions
for this particular case have been obtained, in a similar way as the
Hydrogenic solutions [30]. Further, we have calculated the exact bound state
energy eigenvalues and the corresponding normalized hyper wave functions for
various diatomic molecules for any $l$ angular momentum quantum number bound
by an exactly solvable Kratzer-type potential. The advantage of the present
method is that it enables one to find the energy eigenvalues directly in a
simple way. The method presented in this work can be applied to find the
energy eigenvalues and the corresponding eigenfunctions of the
Schr\"{o}dinger equation within a given potential for various diatomic
molecules.

It is obvious that Eqs. (27) and (36) (with Eqs. (31) and (40)) reduce to
the well-known eigenvalues and eigenfunctions for the $3D$ Kratzer-Fues
potential when $D=3$ [38,39]. Further, they reduce to the Coulomb-like
solutions in $D$-dimensions if $A=C=0$ [40]. Similarly, for $B=2Z>0$ and $%
A=\lambda ,$ $Z$ and $\lambda $ are real constants [41].

\acknowledgments This research was partially supported by the Scientific and
Technical Research Council of Turkey. S. M. Ikhdair wishes to dedicate this
work to his family for their love and assistance.\newpage

