1 Introduction

One of the main objectives of quantum mechanics is to determine the exact analytical solutions of Schrödinger equation with a given potential for the quantum mechanical system. This is due to the fact that the complete wave function involves implicitly all the required information to completely define the observable quantities of the quantum system under consideration; furthermore, a thorough analysis of the diatomic molecular spectroscopy has been traced back to the bound state energy spectra of the atomic and molecular systems under investigation.

The KF potential, is defined as follows [1, 2]:

$$V\left(r\right)= \frac{-A}{r}+ \frac{B}{{r}^{2}}, A=2{D}_{o}{\text{r}}_{\text{e}}, B={D}_{o}{{r}_{e}}^{2},$$
(1)

where \({D}_{o}\) is the bond interaction dissociation energy between two atoms in the DM separated by an equilibrium internuclear distance \({r}_{e}\). This potential is a superposition of the well-known attractive Coulomb potential at large distances and a repulsive centrifugal potential barrier at short distances. Their combination generates an effective potential well, which serves as a potential model to investigate the molecular structure, the molecular energy spectra, the chemical interactions, and the internuclear vibration [3, 4] of the diatomic molecules. This potential belongs to the class of the diatomic molecular potentials [5] for its ability to specify the range of attraction and the strength of the repulsive interactions.

Different methods have been used to study diatomic energy spectra using diatomic molecular potentials: the NU method [6,7,8,9,10,11,12], the exact method [13], the shifted 1/N expansion [14], the exact quantization rule method [15], the SUSY approach [16], the tridiagonal J-matrix representation [17, 18], and algebraic methods [19, 20].

Z.Yalçin et al. [21] have obtained the eigen-energies of the excited 1 and 3 S states of the He atom with the molecular KF potential by using the hyper-spherical harmonics method. Furthermore, they have demonstrated that the KF potential well is sharper and deeper than the pure Coulombic potential well. Also, E. Ikhdair et al. [22] have determined the exact solutions of the D-dimensional Schrödinger equation with the KF potential by means of the standard method. Moreover, Shi-hai Dong and Gus-hua Sun [23] have obtained the exact analytical solutions of the radial Schrödinger equation with the KF potential in D dimensions. In addition, they have studied the dependence of the energy eigenvalues on the spatial dimension D. The KF potential has been extensively investigated by means of the factorization method [24], the asymptotic iteration method [25], among other analytical methods.

The aim of the present work is to take a different approach and solve Schrödinger equation analytically with the KF potential by using the NU method and to obtain the bound state energy levels for the homonuclear and heteronuclear diatomic molecules N2, H2, O2, NO, HCl, CH and LiH.

This work is organized as follows: The NU method is introduced briefly in Sect. 2. In Sect. 3, the exact eigenvalues are obtained and their corresponding radial eigenfunctions are analytically formulated by solving the radial Schrödinger equation within the context of the NU method for KF potential. In Sect. 4, we performed numerical calculations of the rovibrational energy levels for the above-mentioned DMs. Also, a graphical representation and discussion of the effective potential well are given. Finally, our work is concluded in Sect. 5 with a summary and concluding remarks.

2 The Nikiforov-Uvarov Method

The NU method [26, 27] is dedicated to solving the hypergeometric type equations of the form:

$${\psi }^{{\prime }{\prime }}\left(x\right)+ \frac{\tilde{\tau }\left(x\right)}{\sigma \left(x\right)}{\psi }^{{\prime }}\left(x\right)+ \frac{\tilde{\sigma }\left(x\right)}{{\sigma }^{2}\left(x\right)}\psi \left(x\right)=0,$$
(2)

where \(\tilde{\tau }\left(x\right)\) is a polynomial of degree at most one, \(\sigma \left(x\right)\) and \(\tilde{\sigma }\left(x\right)\) are polynomials of degree at most two, \(\psi \left(x\right)\) is a hypergeometric-type function, via orthogonal polynomials. Furthermore, the function \(\psi \left(x\right)\) takes the form:

$$\psi \left(x\right)= {\upvarphi }\left(x\right)\text{Y}\left(x\right),$$
(3)

Equation (2) can be reformulated to a hypergeometric type of equation of the form:

$$\sigma \left(x\right){Y}^{{\prime }{\prime }}\left(x\right)+ \tau \left(x\right){Y}^{{\prime }}\left(x\right)+ \lambda Y\left(x\right)=0,$$
(4)

to obtain its solutions. We introduce the following function of a maximum first degree:

$$\pi \left(x\right)= \frac{{\sigma }^{{\prime }}\left(x\right)- \tilde{\tau }\left(x\right)}{2} \pm \sqrt{{\left(\frac{{\sigma }^{{\prime }}\left(x\right)- \tilde{\tau }\left(x\right)}{2}\right)}^{2}-\tilde{\sigma }\left(x\right)+k\sigma \left(x\right)}.$$
(5)

The maximum first degree of \(\pi \left(x\right)\) imposes a restriction on the expression under the radical, its discriminant must be equal to zero. Equating the discriminant to zero yields a quadratic equation for \(k\), which can be solved algebraically to get the possible values of \(k\) and the corresponding values of \(\pi \left(x\right)\). It is worth noting that, according to the method, we must select the values of \(\pi \left(x\right)\) for which \({\tau }^{{\prime }}\left(x\right)<0\).

The hypergeometric type of function \(Y\left(x\right)\) in Eq. (3) is obtained by Rodrigus formula:

$${Y\left(x\right)=Y}_{n}\left(x\right)= \frac{{B}_{n}}{\rho \left(x\right)}\frac{{d}^{n}}{d{x}^{n}}\left({\sigma }^{n}\left(x\right)\rho \left(x\right)\right),$$
(6)

where \({B}_{n}\) is the normalization constant, \(\rho \left(x\right)\) is the weight function satisfying the conditions:

$${\nu }^{{\prime }}\left(x\right)= \frac{\tau \left(x\right)}{\sigma \left(x\right)}\nu \left(x\right), \nu \left(x\right)= \sigma \left(x\right)\rho \left(x\right)$$
(7)
$$\tau \left(x\right)= \tilde{\tau }\left(x\right)+2\pi \left(x\right).$$
(8)

According to the NU method, the parameter \(\lambda\) is defined by:

$$\lambda =k+ {\pi }^{{\prime }}\left(x\right).$$
(9)

Accordingly, the eigenvalues, Eq. (4), can be obtained from the following equation:

$$\lambda = {\lambda }_{n}= -n{\tau }^{{\prime }}\left(x\right)- \frac{n\left(n-1\right)}{2}{\sigma }^{{\prime }{\prime }}\left(x\right) n=\text{0,1},\text{2,3},\dots.$$
(10)

It is worth noting that \(\lambda\) and \({\lambda }_{n}\) are obtained from a particular result of \({Y\left(x\right)=Y}_{n}\left(x\right).\) The function \(\varphi \left(x\right)\), in Eq. (3), is obtained from the following logarithmic derivative:

$$\frac{{\varphi }^{{\prime }}\left(x\right)}{\varphi \left(x\right)}=\frac{\pi \left(x\right)}{\sigma \left(x\right)}.$$
(11)

3 Bound state eigenvalues and their corresponding Radial Eigenfunctions for the KF potential by means of NU Method

The radial equation for KF potential is of the form:

$$\frac{{d}^{2}{R}_{nl}\left(r\right)}{d{r}^{2}}+ \frac{2}{r}\frac{d{R}_{nl}\left(r\right)}{dr}+ \frac{2\mu }{{\hslash }^{2}{r}^{2}}\left({E}_{nl}{r}^{2}+Ar-B- \frac{{\hslash }^{2}l\left(l+1\right)}{2\mu }\right){R}_{nl}\left(r\right)=0,$$
(12)

where \(\mu\) is the reduced mass, \(l\) is the orbital angular momentum quantum number, and \(n\) is the vibrational quantum number. Introducing the following parameters:

$$\alpha =\frac{2\mu }{{\hslash }^{2}}, \beta =l(l+1), {E}_{nl}= -\omega , \omega >0$$
(13)

the radial equation can be expressed as

$$\frac{{d}^{2}{R}_{nl}\left(r\right)}{d{r}^{2}}+ \frac{2}{r}\frac{d{R}_{nl}\left(r\right)}{dr}+\frac{1}{{r}^{2}}\left(-\alpha \omega {r}^{2}+\alpha Ar-\left(\alpha B+\beta \right)\right){R}_{nl}\left(r\right)=0.$$
(14)

Comparing Eq. (14) with Eq. (2), we obtain the following:

$$\tilde{\tau }\left(r\right)=2, \sigma \left(r\right)=r, \tilde{\sigma }\left(r\right)= -\alpha \omega {r}^{2}+ \alpha Ar-(\alpha B+\beta )$$
(15)

Hence, by means of Eq. (5) we can write \(\pi \left(r\right)\) as follows:

$$\pi \left(r\right)= -\frac{1}{2}\pm \frac{1}{2}\sqrt{4\alpha \omega {r}^{2}+4\left(k-\alpha A\right)r+4\left(\alpha B+\beta \right)+1}.$$
(16)

To determine the constant \(k\), the discriminant must be equal to zero, and as a result we obtain

$${k}_{\pm }= \alpha A \pm \sqrt{\alpha \omega \left(4\left(\alpha B+\beta \right)+1\right)}.$$
(17)

Inserting Eq. (17) into Eq. (16), \(\pi \left(r\right)\) can be rewritten as

$$\pi \left(r\right)=\left\{\begin{array}{c}-\frac{1}{2}\pm \frac{1}{2}\left[\sqrt{4\alpha \omega }r+\sqrt{4\left(\alpha B+\beta \right)+1}\right] ,for {k}_{+}\\ -\frac{1}{2}\pm \frac{1}{2}\left[\sqrt{4\alpha \omega }r-\sqrt{4\left(\alpha B+\beta \right)+1}\right] ,for {k}_{-}\end{array}\right\}$$
(18)

Accordingly, we select a particular expression for \(\pi \left(r\right)\) which meets the condition \({\tau }^{{\prime }}\left(r\right)<0\) as

$$\pi \left(r\right)= -\frac{1}{2}\left(1-\gamma +\sqrt{4\alpha \omega }r\right),\gamma =\sqrt{4\left(\alpha B+\beta \right)+1}.$$
(19)

From Eq. (8), \(\tau \left(r\right)\) is given by:

$$\tau \left(r\right)=\left(\gamma +1\right)-\sqrt{4\alpha \omega }r.$$
(20)

From Eq. (9) and Eq. (10), the parameters \(\lambda\) and \({\lambda }_{n}\) are given by:

$$\lambda = \alpha A-(\gamma +1)\sqrt{\alpha \omega },$$
(21)
$${\lambda }_{n}=2n\sqrt{\alpha \omega }.$$
(22)

Equating the two equations above, we obtain

$$\omega =\frac{\alpha {A}^{2}}{{\left(\left(2n+1\right)+\gamma \right)}^{2}}.$$
(23)

Substituting for the values of the parameters \(\alpha\), \(\gamma\), and\(\omega\) we get the bound state eigenvalues in the form:

$${E}_{nl}= \frac{-2\mu {A}^{2}}{{\hslash }^{2}{\left(2n+1+\sqrt{{\left(2l+1\right)}^{2}+\frac{8\mu B}{{\hslash }^{2}}}\right)}^{2}}, n=\text{0,1},2,\dots ,l=\text{0,1},2,\dots ,n.$$
(24)

Equation (24) is identical to the energy eigenvalue equation provided by Eq. (28) of Ref. [28], if we take the substitution \(N=3\). Hence, we finally obtain the energy eigenvalue equation:

$${E}_{nl}= -\frac{8\mu }{{\hslash }^{2}}{\left(\frac{{D}_{o}{r}_{e}}{\left(2n+1\right)+\sqrt{{\left(2l+1\right)}^{2}+\frac{8\mu {D}_{o}{{r}_{e}}^{2}}{{\hslash }^{2}}}}\right)}^{2}.$$
(25)

To determine the normalized radial eigenfunctions \({R}_{nl}\) in terms of the orthogonal associated Laguerre polynomials, consider Eq. (11) from which we can recast the logarithmic derivative as follows

$$ln \varphi \left(r\right)= \int \frac{\pi \left(r\right)}{\sigma \left(r\right)}dr$$
(26)

Inserting the expressions for \(\pi \left(r\right)\) and \(\sigma \left(r\right)\) into the above equation we yield

$$\varphi \left(r\right)= {N}_{nl}{r}^{\left(\frac{\gamma -1}{2}\right)}{e}^{-\sqrt{\alpha \omega }r}.$$
(27)

Applying the same logarithmic derivative transformation to Eq. (7) yields

$$ln\left(\nu \left(r\right)\right)=\int \frac{\tau \left(r\right)}{\sigma \left(r\right)}dr,$$
(28)

from which we get

$$\rho \left(r\right)= c{r}^{\gamma }{e}^{-\delta r} , \delta =\sqrt{4\alpha \omega }.$$
(29)

By means of Eq. (6), \({Y}_{n}\left(r\right)\) is given by

$${Y}_{n}\left(r\right)= {B}_{n}{r}^{-\gamma }{e}^{\delta r}\frac{{d}^{n}}{d{r}^{n}}\left({r}^{n+\gamma }{e}^{-\delta r}\right),$$
(30)

which can be compared to the well-known Rodrigues formula for the orthogonal associated Laguerre polynomials [29]

$${L}_{n}^{k}\left(\xi \right)=\frac{{\xi }^{-k}{e}^{\xi }}{n!}\frac{{d}^{n}}{d{\xi }^{n}}\left({\xi }^{n+k}{e}^{-\xi }\right),$$
(31)

by taking \({B}_{n}= \frac{1}{n!}\), \(k=\gamma\), \(\xi =\delta r\). Hence, we obtain:

$${Y}_{n}\left(r\right)={L}_{n}^{\gamma }\left(\delta r\right).$$
(32)

Finally, from Eqs. (27) and (32) the normalized radial eigenfunctions can be expressed in terms of the orthogonal Laguerre polynomials as follows

$${R}_{nl}\left(r\right)= {N}_{nl}{r}^{\left(\frac{\gamma -1}{2}\right)}{e}^{-\frac{\delta }{2}r}{L}_{n}^{\gamma }\left(\delta r\right).$$
(33)

With the aid of the normalization condition

$$\underset{0}{\overset{\infty }{\int }}{r}^{2}{\left({R}_{nl}\left(r\right)\right)}^{2}dr=1,$$
(34)

and the result of Ref. [29]

$${\int }_{0}^{\infty }{t}^{\alpha +\beta }{e}^{-t}{\left({L}_{n}^{\alpha }\left(t\right)\right)}^{2}dt= {\sum }_{r=0}^{n}{\left(\begin{array}{c}\beta \\ n-r\end{array}\right)}^{2}\frac{\varGamma (\alpha +\beta +r+1)}{r!},$$
(35)

the normalization constant \({N}_{nl}\) takes the form:

$${N}_{nl}= {\delta }^{\frac{\gamma }{2}+1}\sqrt{\frac{n!}{(2n+\gamma +1)\varGamma (n+\gamma +1)}}.$$
(36)

Substituting Eq. (36) into Eq. (33), we finally get

$${R}_{nl}\left(\rho \right)= \sqrt{\frac{{\delta }^{3}\varGamma (n+1)}{(2n+\gamma +1)\varGamma (n+\gamma +1)}}{\rho }^{\left(\frac{\gamma -1}{2}\right)}{e}^{-\frac{\rho }{2}}{L}_{n}^{\gamma }\left(\rho \right), \rho =\delta r.$$
(37)

4 Results and discussion

The energy levels for the diatomic molecules HCl, CH, N2, H2, O2, NO, and LiH are computed numerically and tabulated in Tables 1, 2, 3, 4 and 5, by means of spectroscopic parameters presented in Table 6. The eigenvalues obtained in Tables 1, 2, 3 and 4 are in excellent agreement with the values computed previously by means of alternative methods [25], [30,31,32], [34]. To the best of our knowledge, our work is the first to quantify the energy levels of N2, H2, and CH DMs in the context of the KF interaction. With the aim of providing a comparative benchmark against which future studies can be compared, we have only included our data for these DMs in Table 6.

Table 1 Spectroscopic parameters for different diatomic molecules in the ground state
Table 2 Energy eigenvalues (in eV) of KF molecular potential for different values of n and \(l\)for ground state \({\text{O}}_{2}\) diatomic molecule, where \(\hslash c=\) 1973.29 eV Å
Table 3 Energy eigenvalues (in eV) of KF molecular potential for different values of n and \(l\) for ground state LiH diatomic molecule, where \(\hslash c=\) 1973.29 eV Å
Table 4 Energy eigenvalues (in eV) of KF molecular potential for different values of n and \(l\) for ground state HCL diatomic molecule, where \(\hslash c=\) 1973.29 eV Å
Table 5 Energy eigenvalues (in eV) of KF molecular potential for different values of n and \(l\) for ground state NO diatomic molecule, where \(\hslash c=\) 1973.29 eV Å
Table 6 Energy eigenvalues (in eV) of KF molecular potential for different values of n and \(l\) for ground states \({\text{N}}_{2}\), \({\text{H}}_{2}\), and CH diatomic molecules, where \(\hslash c=\) 1973.29 eV Å
Fig. 1
figure 1

Kratzer-Feus-type potential for different diatomic molecules.

In Fig. 1, the KF potential is plotted for several DMs; its graphical representation reveals the nature of the chemical bond and the molecular behavior at r = \({r}_{e}.\) As shown in Fig. 1, the potential becomes infinite as the internuclear separation approaches zero due to the internuclear repulsion. On the other hand, as the molecule decomposes, the internuclear separation tends to infinity, and consequently, the potential vanishes. According to Tables 1, 2, 3, 4 and 5, and Fig. 1, the eigenvalues of N2 are greater than those of the other diatomic molecules, since its potential depth is sharper and larger than the other ones. Therefore, we detect an increase in the energy eigenvalues as the molecular potential depth increases.

5 Conclusion

In the present paper, we solved Schrödinger equation for KF potential to get energy eigenvalues and the corresponding normalized total wavefunctions. Numerical results are represented for some well-known diatomic molecules. The findings of this work demonstrate the remarkable precision with which the KF potential models the diatomic molecular structures and quantifies their bound state energy levels. Our results can be extended into applications in future research studies.