Quantum features of molecular interactions associated with time-dependent non-central potentials

The study of quantum properties of dynamical systems involving a non-central potential has attracted a considerable interest in the literature because there are plentiful applications of such systems in both physics and chemistry. In particular, quantum molecular interactions [1–7] can be described in terms of such non-central potentials. Concerning this, there are lots of reports on theoretical studies for characteristics of interacting molecules based on the Schrödinger solutions [7–21]. For instance, quantum features of molecular interactions which are subjected to ring-shaped potentials [22–25] or double ring-shaped potentials [26–31] in quantum chemistry were extensively investigated. The mechanism of interaction between deformed nuclei has been analyzed [32, 33]. Bound states of an electron in a Coulomb-type potential in the presence of an Aharonov–Bohm field have also been studied [34, 35]. The physical analysis of molecular interactions are crucial for understanding the underlying mechanism of diverse modern technology tied to materials science, such as nanotechnology, sensor efficiency, protein folding, DNA condensation, drug design, etc.

The purpose of this work is to review recent progress in the research area relevant to quantum characteristics of molecular interactions governed by generalized time-dependent non-central potentials, where non-central potentials are especially composed of Coulomb potential together with an inverse quadratic potential. We also introduce a double ring-shaped generalized non-central time-dependent potential in order to investigate such molecular interactions. To attain accurate results when we study molecular systems, it is necessary to introduce an exact Hamiltonian that yields actual time dependence of molecule behaviors. If we consider the convention that time-varying factors have usually been neglected on most studies of dynamical systems, the recent tendency [3, 36–48] for considering time dependence of physical parameters in this field may open up a new trend in the analysis of molecular interactions.

Because the potentials we treat in this work are somewhat complicated time-dependent forms, usual treatment of them in view of quantum mechanics on the basis of separation of variables method may not be available. To overcome this difficulty, other potential methods for solving quantum solutions, such as the invariant operator method [49, 50], the unitary transformation method [48], and the Nikiforov–Uvarov (NU) method [16, 19, 51], can be used.

It is well known that quantum solutions (the Schrödinger solutions) of time-dependent Hamiltonian systems are represented in terms of the eigenstates of an invariant operator of the system [49, 50]. This is the reason why the invariant operator method is useful in this time-dependent situation. In general, it is not an easy task to solve the eigenvalue equation of the invariant operator for a complicated system. Because of this difficulty, we will transform the eigenvalue equation to a simple form by a unitary operator. Then, we can easily solve the eigenvalue equation in the transformed system by taking advantage of special mathematical techniques. The NU method and the asymptotic iteration method (AIM) will be used as such mathematical manipulations in this review. Through these processes, complete solutions of the eigenvalue equation of the invariant operator can be obtained and, consequently, we can identify the corresponding eigenstates. This is the basic strategy that has been used in order to solve the Schrödinger equation of the system. Quantum properties of the system will be studied on the bases of Schrödinger solutions obtained in such a way.

Quantum problem of molecular interactions such as the Coulomb-type and harmonic interactions can be described in terms of central [36–38] or noncentral [7, 39–41] time-dependent Hamiltonians. Attractive or repulsive forces between various molecules including non-bonded atoms are responsible for a specific formation of molecular structure and its change. Hence, the study for the physical mechanism of molecular interactions plays a key role to the understanding of the structure of condensed matters and biological systems as well as chemical reactions and processes.

We consider a Hamiltonian that describes the motion of a particle subjected to a time-dependent non-central potential, which is

$$\begin{eqnarray}&&H(t)=\displaystyle \frac{1}{2\mu (t)}\left({p}^{2}+\displaystyle \frac{g(\theta ,\varphi )}{{r}^{2}}\right)-\displaystyle \frac{Z(t)}{r},\end{eqnarray} \tag{ 1 }$$

where $\mu (t)$ and Z(t) are time-dependent coefficients with the condition $Z(t)\gt 0$, and $g(\theta ,\varphi )$ is a function of θ and . In effect, $\mu (t)$ stands for effective mass that varies with time. Recently, molecular systems and other dynamical systems with time-dependent effective mass, as well as position-dependent effective mass, became a topic of research [7, 37, 40, 52–60]. In some cases, effective mass of a molecule in a system may vary through its interaction with the environment or various excitations such as energy, temperature, stress, pressure, phonon, etc. If the environment or degree of such excitations changes through time, the resultant effective mass of the molecule may also be subjected to a change in time.

The momentum operator in the spherical coordinate is represented as ${p}^{2}={p}_{r}^{2}+\tfrac{{L}^{2}}{{r}^{2}}$, where ${p}_{r}=-{\rm{i}}{\hslash }(\tfrac{\partial }{\partial r}+\tfrac{1}{r}\phantom{\rule{}{3.0ex}})$, and the total angular momentum is given by

$$\begin{eqnarray}&&{L}^{2}=-{{\hslash }}^{2}\left[\displaystyle \frac{1}{{\sin }^{2}\theta }\displaystyle \frac{{\partial }^{2}}{\partial {\varphi }^{2}}+\displaystyle \frac{1}{\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \displaystyle \frac{\partial }{\partial \theta }\right)\right].\end{eqnarray} \tag{ 2 }$$

Due to $g(\theta ,\varphi )$, the potential in the Hamiltonian, equation (1), is non-central. Such potential is useful for describing time behavior of intricate non-central molecular interactions.

Let us denote the Schrödinger state of the system as ${\rm{\Psi }}(t)$. Then, we can write the corresponding Schrödinger equation in the form

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial }{\partial t}{\rm{\Psi }}(t)=H(t){\rm{\Psi }}(t).\end{eqnarray} \tag{ 3 }$$

As mentioned in the introductory part, it is in general useful to introduce an invariant operator in order to study quantum properties of a system whose Hamiltonian is dependent on time. One can construct the invariant I(t) from the Liouville–von Neumann equation that is given by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{I}}{{\rm{d}}{t}}=\displaystyle \frac{\partial I}{\partial t}+\displaystyle \frac{1}{{\rm{i}}{\hslash }}[I,H]=0.\end{eqnarray} \tag{ 4 }$$

To find the expression of the invariant I(t) from this equation, let us put

$$\begin{eqnarray}&&I(t)=A(t)\left({p}^{2}+\displaystyle \frac{g(\theta ,\varphi )}{{r}^{2}}\right)+B(t)({{rp}}_{r}+{p}_{r}r)+C(t){r}^{2}-D(t)\displaystyle \frac{1}{r},\end{eqnarray} \tag{ 5 }$$

where ${A}({t}),{B}({t}),{C}({t}),$ and $D(t)$ are time-dependent coefficients that will be evaluated from now on. Let us suppose that $g(\theta ,\varphi )$ varies sufficiently slowly with the variation of angles θ and so that its derivatives with respect to angles can be neglected: ${\rm{d}}{g}(\theta ,\varphi )/{\rm{d}}\theta \sim {\rm{d}}{g}(\theta ,\varphi )/{\rm{d}}\varphi \sim 0$. Then, a minor mathematical procedure with this assumption after substituting equations (1) and (5) into equation (4) yields

$$\begin{eqnarray}\dot{A}(t) & = & -2{\mu }^{-1}(t)B(t)\,\,\,\,\,\dot{B}(t)=-{\mu }^{-1}(t)C(t),\,\,\,\,\,\dot{C}(t)=0,\\ \dot{D}(t) & = & -2Z(t)B(t),\,\,\,\,\,\,A(t)Z(t)={[2\mu (t)]}^{-1}D(t).\end{eqnarray} \tag{ 6 }$$

Evidently, the third relation implies that $C(t)$ is a constant: $C(t)={C}_{0}$. Then, from fundamental mathematical procedures with equation (6), other coefficients are determined in terms of C₀. The detailed evaluations for this result in

$$\begin{eqnarray}&&A(t)=\eta (t),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&B(t)={B}_{0}-{C}_{0}{\int }_{0}^{t}\displaystyle \frac{1}{\mu ({t}^{{\prime} })}{{\rm{d}}{t}}^{{\prime} },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&D(t)=\displaystyle \frac{{D}_{0}}{{A}_{0}^{1/2}}{\eta }^{1/2}(t),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\eta (t)={A}_{0}-2{B}_{0}{\int }_{0}^{t}\displaystyle \frac{1}{\mu ({t}^{{\prime} })}{{\rm{d}}{t}}^{{\prime} }+2{C}_{0}{\int }_{0}^{t}\left(\displaystyle \frac{1}{\mu ({t}^{{\prime} })}{\int }_{0}^{{t}^{{\prime} }}\displaystyle \frac{1}{\mu ({t}^{{\prime\prime} })}{{\rm{d}}{t}}^{{\prime\prime} }\right){{\rm{d}}{t}}^{{\prime} },\end{eqnarray} \tag{ 10 }$$

where ${A}_{0}=A(0)$, ${B}_{0}=B(0)$, and ${D}_{0}=D(0)$. By inserting equations (7)–(9) into equation (5), we obtain [39]

$$\begin{eqnarray}I(t) & = & \eta (t)\left({p}^{2}+\displaystyle \frac{g(\theta ,\varphi )}{{r}^{2}}\right)+\left({B}_{0}-{C}_{0}{\displaystyle \int }_{0}^{t}\displaystyle \frac{1}{\mu ({t}^{{\prime} })}{{\rm{d}}{t}}^{{\prime} }\right)({{rp}}_{r}+{p}_{r}r)\\ & & +{C}_{0}{r}^{2}-\displaystyle \frac{{D}_{0}}{{A}_{0}^{1/2}}{\eta }^{1/2}(t)\displaystyle \frac{1}{r}.\end{eqnarray} \tag{ 11 }$$

We note that we have levied some conditions in order to obtain this formula of the invariant operator as a solvable case: A₀, B₀, C₀, and D₀ are real, while the function Z(t) has the form [39]

$$\begin{eqnarray}&&Z(t)=\displaystyle \frac{{D}_{0}}{2{A}_{0}^{1/2}\mu (t)}{\eta }^{-1/2}(t).\end{eqnarray} \tag{ 12 }$$

Because the quantum wave functions of the system is represented in terms of the eigenstates of I(t), it is now necessary to evaluate its eigenvalue equation. Let us write the eigenvalue equation of the invariant I(t) as

$$\begin{eqnarray}&&I(t){\phi }_{n}(\overrightarrow{r},t)={E}_{n}{\phi }_{n}(\overrightarrow{r},t),\end{eqnarray} \tag{ 13 }$$

where E_n are the eigenvalues and ${\phi }_{n}(\overrightarrow{r},t)$ are time-dependent eigenstates. It may be not an easy task to solve this equation because I(t) given in equation (11) is a very complicated form. Hence, it is favorable to solve equation (13) after we transform it into a simple form. To do this, we consider a unitary transformation of the form

$$\begin{eqnarray}&&{\phi }_{n}^{{\prime} }(\overrightarrow{r})=U(t){\phi }_{n}(\overrightarrow{r},t),\end{eqnarray} \tag{ 14 }$$

where U(t) is a time-dependent unitary operator that is giving by

$$\begin{eqnarray}&&U(t)=\exp \left(\displaystyle \frac{{\rm{i}}{B}(t)}{2{\hslash }{A}_{0}}{r}^{2}\right)\exp \left[\displaystyle \frac{-{\rm{i}}}{2{\hslash }}\mathrm{ln}{\left(\displaystyle \frac{\eta (t)}{{A}_{0}}\right)}^{1/2}({{rp}}_{r}+{p}_{r}r)\right].\end{eqnarray} \tag{ 15 }$$

Then, we can easily verify that the eigenvalue equation reduces in terms of ${\phi }_{n}^{{\prime} }(\overrightarrow{r})$ as

$$\begin{eqnarray}&&{{UIU}}^{-1}{\phi }_{n}^{{\prime} }(\overrightarrow{r})={I}_{0}{\phi }_{n}^{{\prime} }(\overrightarrow{r})=E{\phi }_{n}^{{\prime} }(\overrightarrow{r}).\end{eqnarray} \tag{ 16 }$$

After a minor evaluation, we obtain that [39]

$$\begin{eqnarray}&&\left[{A}_{0}\left({p}^{2}+\displaystyle \frac{g(\theta ,\varphi )}{{r}^{2}}\right)+({C}_{0}{A}_{0}-{B}_{0}^{2}){r}^{2}-\displaystyle \frac{{D}_{0}}{r}\right]{\phi }_{n}^{{\prime} }(r,\theta ,\varphi )=E{\phi }_{n}^{{\prime} }(r,\theta ,\varphi ).\end{eqnarray} \tag{ 17 }$$

Note that this equation is very simple compared to the original eigenvalue equation represented in equation (13). For convenience, let us put ${\omega }_{0}={C}_{0}{A}_{0}-{B}_{0}^{2}$. Than, the solution of equation (17) varies depending on the value or sign of ${\omega }_{0}$. Because A₀, B₀, and C₀ are arbitrary real constants, we can adjust them to be ${\omega }_{0}=0$. Because of the merit of this procedure that simplifies the problem further we only consider the case ${\omega }_{0}=0$ in this work. Then, the eigenvalue equation can be rewritten as

$$\begin{eqnarray}&&-{{\hslash }}^{2}{A}_{0}\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial }{\partial r}\right)+\displaystyle \frac{1}{{r}^{2}\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \displaystyle \frac{\partial }{\partial \theta }\right)+\displaystyle \frac{1}{{r}^{2}{\sin }^{2}\theta }\displaystyle \frac{{\partial }^{2}}{\partial {\varphi }^{2}}\right]{\phi }_{{nlm}}^{{\prime} }(r,\theta ,\varphi )\\ &&\quad +\,\left(\displaystyle \frac{{A}_{0}g(\theta ,\varphi )}{{r}^{2}}-\displaystyle \frac{{D}_{0}}{r}\right){\phi }_{{nlm}}^{{\prime} }(r,\theta ,\varphi )={E}_{{nlm}}{\phi }_{{nlm}}^{{\prime} }(r,\theta ,\varphi ).\end{eqnarray} \tag{ 18 }$$

According to the invariant operator theory, the wave functions of the system are represented in terms of ${\phi }_{{nlm}}^{{\prime} }(r,\theta ,\varphi )$. Because this equation is independent of time, we can apply the separation of variables method in order to solve it. Considering this, we put

$$\begin{eqnarray}&&{\phi }_{{nlm}}^{{\prime} }(r,\theta ,\varphi )={R}_{{nl}}(r){Y}_{{lm}}(\theta ,\varphi ).\end{eqnarray} \tag{ 19 }$$

Some rearrangements after substituting this in equation (18) lead to two separated equations as [39, 40]

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}R(r)}{{{\rm{d}}{r}}^{2}}+\displaystyle \frac{2}{r}\displaystyle \frac{{\rm{d}}{R}(r)}{{\rm{d}}{r}}+\left(\displaystyle \frac{{E}_{n}}{{{\hslash }}^{2}{A}_{0}}-\displaystyle \frac{c}{{r}^{2}}+\displaystyle \frac{{D}_{0}}{{{\hslash }}^{2}{A}_{0}r}\right){R}_{{nl}}(r)=0,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}}+\cot \theta \displaystyle \frac{\partial }{\partial \theta }+\displaystyle \frac{1}{{\sin }^{2}\theta }\displaystyle \frac{{\partial }^{2}}{\partial {\varphi }^{2}}+c-\displaystyle \frac{g(\theta ,\varphi )}{{{\hslash }}^{2}}\right){Y}_{{lm}}(\theta ,\varphi )=0,\end{eqnarray} \tag{ 21 }$$

where c is a separation constant. The first equation is the radial equation, while the second is for angular parts.

In this section, we will derive the solutions for the radial part of the eigenvalue equation given in equation (20). Because this equation involves a time-dependent non-central potential, we need an alternative powerful method for treating it. Hence, we will use the method of NU [51] for this purpose. For convenience of the use of this method at a later time, we transform the variables as s = r. Then, equation (20) can be rewritten in terms of s as [40]

$$\begin{eqnarray}&&{R}^{{\prime\prime} }(s)+\displaystyle \frac{2}{s}{R}^{{\prime} }(s)+\displaystyle \frac{1}{{s}^{2}}[-{K}_{n}^{2}{s}^{2}+\varepsilon s-l(l+1)]R(s)=0,\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{nl}}}{{{\hslash }}^{2}{A}_{0}}=-{K}_{n}^{2},\,\,\,\,\,\displaystyle \frac{{D}_{0}}{{{\hslash }}^{2}{A}_{0}}=\varepsilon ,\,\,\,\,\,c=l(l+1),\end{eqnarray} \tag{ 23 }$$

with the condition ${E}_{{nl}}\lt 0$.

Before the use of the NU method for solving this equation, we will briefly mention its general methodology. The NU method was originally developed in order to reduce a second order differential equation to a hypergeometric equation through a coordinate transformation of the form $s=s(r)$. Let us consider the following general equation written in terms of s

$$\begin{eqnarray}&&{{\rm{\Psi }}}^{{\prime\prime} }(s)+\displaystyle \frac{\tilde{\tau }(s)}{\sigma (s)}{{\rm{\Psi }}}^{{\prime} }(s)+\displaystyle \frac{\tilde{\sigma }(s)}{\sigma (s)}{\rm{\Psi }}(s)=0,\end{eqnarray} \tag{ 24 }$$

where $\sigma (s)$ and $\tilde{\sigma }(s)$, and $\tilde{\tau }(s)$ are some polynomials. We note that $\tilde{\tau }(s)$ is a polynomial of s, at most, in the first order. Now we put ${\rm{\Psi }}(s)$ in the form [19]

$$\begin{eqnarray}&&{\rm{\Psi }}(s)=w(s)y(s),\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&\sigma (s){y}^{{\prime\prime} }(s)+\tau (s){y}^{{\prime} }(s)+\lambda y(s)=0,\end{eqnarray} \tag{ 26 }$$

and w(s) is represented in terms of a logarithmic derivative such that

$$\begin{eqnarray}&&\displaystyle \frac{{w}^{{\prime} }(s)}{w(s)}=\displaystyle \frac{\pi (s)}{\sigma (s)}.\end{eqnarray} \tag{ 27 }$$

Note that y(s) is expressed by a hypergeometric function and this polynomial solution is relevant to the Rodrigues relation which is [19]

$$\begin{eqnarray}&&{y}_{n}(s)=\displaystyle \frac{{B}_{n}}{\rho (s)}\displaystyle \frac{{{\rm{d}}}^{n}}{{{\rm{d}}{s}}^{n}}[{\sigma }^{n}(s)\rho (s)],\end{eqnarray} \tag{ 28 }$$

where B_n are normalization constants and $\rho (s)$ is a weight function which yields ${(\sigma \rho )}^{{\prime} }=\tau \rho$.

For further development of the NU theory, we introduce a function π and a useful parameter λ, which are [19]

$$\begin{eqnarray}&&\pi (s)=\displaystyle \frac{1}{2}[{\sigma }^{{\prime} }(s)-\tilde{\tau }(s)]\pm {\left[{\left(\displaystyle \frac{{\sigma }^{{\prime} }(s)-\tilde{\tau }(s)}{2}\right)}^{2}-\tilde{\sigma }(s)-k\sigma (s)\right]}^{1/2},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\lambda (s)=k+{\pi }^{{\prime} }(s),\end{eqnarray} \tag{ 30 }$$

where k will be derived afterwards. If we consider that the expression given inside the square root should be given in the form which is a square of a polynomial, the representation of λ should become [19]

$$\begin{eqnarray}&&\lambda (s)={\lambda }_{n}(s)=-n{\tau }^{{\prime} }(s)-\displaystyle \frac{n(n-1)}{2}{\sigma }^{{\prime\prime} }(s),\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&\tau (s)=\tilde{\tau }(s)+2\pi (s).\end{eqnarray} \tag{ 32 }$$

If we rewrite the functions as $\tilde{\tau }={\alpha }_{1}-{\alpha }_{2}s$, $\sigma =s(1-{\alpha }_{3}s)$, and $\tilde{\sigma }=-{\zeta }_{1}{s}^{2}+{\zeta }_{2}s-{\zeta }_{3}$, equation (24) becomes [42, 61].

$$\begin{eqnarray}&&\left[\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{s}}^{2}}+\displaystyle \frac{{\alpha }_{1}-{\alpha }_{2}s}{s(1-{\alpha }_{3}s)}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{s}}+\displaystyle \frac{-{\zeta }_{1}{s}^{2}+{\zeta }_{2}s-{\zeta }_{3}}{{s}^{2}{(1-{\alpha }_{3}s)}^{2}}\right]\ {\rm{\Psi }}(s)=0.\end{eqnarray} \tag{ 33 }$$

Considering this formula, we find from equation (29) that

$$\begin{eqnarray}&&\pi (s)={\pi }_{1}(s)\pm {[{\pi }_{2}(s)]}^{1/2},\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&{\pi }_{1}(s)=\displaystyle \frac{1}{2}(1-{\alpha }_{1})+\displaystyle \frac{1}{2}({\alpha }_{2}-2{\alpha }_{3})s,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}{{\pi }}_{2}({s}) & = & \left[\displaystyle \frac{1}{4}{({\alpha }_{2}-2{\alpha }_{3})}^{2}+{\zeta }_{1}-k{\alpha }_{3}\right]{s}^{2}+\left[\displaystyle \frac{1}{2}(1-{{\alpha }}_{1})({{\alpha }}_{2}-2{{\alpha }}_{3})\right.\\ & & -{{\zeta }}_{2}+{k}\phantom{\rule{}{2.25ex}}]\,{s}+\displaystyle \frac{1}{4}{{\delta }}_{1}^{2},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{\delta }_{1}={[{(1-{\alpha }_{1})}^{2}+4{\zeta }_{3}]}^{1/2}.\end{eqnarray} \tag{ 37 }$$

Here, k is given by [61]

$$\begin{eqnarray}&&k={k}_{\mathrm{1,2}}={\beta }_{1}\mp {\beta }_{2}^{1/2},\end{eqnarray} \tag{ 38 }$$

where

$$\begin{eqnarray}&&{\beta }_{1}=-\displaystyle \frac{1}{2}(1-{\alpha }_{1})({\alpha }_{2}-2{\alpha }_{3})+{\zeta }_{2}-\displaystyle \frac{1}{2}{\alpha }_{3}{\delta }_{1}^{2},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{\beta }_{2}=\displaystyle \frac{1}{4}{\alpha }_{3}^{2}{\delta }_{1}^{2}{\delta }_{2}^{2},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\delta }_{2}=2{\left\{\displaystyle \frac{1}{4}\left[{\left(\displaystyle \frac{{\alpha }_{2}}{{\alpha }_{3}}-2\right)}^{2}+{\delta }_{1}^{2}\right]+\displaystyle \frac{1}{{\alpha }_{3}}\left(\displaystyle \frac{1}{2}(1-{\alpha }_{1})({\alpha }_{2}-2{\alpha }_{3})-{\zeta }_{2}\right)+\displaystyle \frac{{\zeta }_{1}}{{\alpha }_{3}^{2}}\right\}}^{1/2}.\end{eqnarray} \tag{ 41 }$$

Let us first see for the case $k={k}_{1}={\beta }_{1}-{\beta }_{2}^{1/2}$. From the formulae given above, we can see that equation (28) becomes

$$\begin{eqnarray}&&{y}_{n}={P}_{n}^{({\delta }_{1},{\delta }_{2})}(1-2{\alpha }_{3}s),\end{eqnarray} \tag{ 42 }$$

where ${P}_{n}^{(\alpha ,\beta )}$ are Jacobi polynomials. This leads to the general solution as

$$\begin{eqnarray}&&{\rm{\Psi }}(s)=w(s)y(s),\end{eqnarray} \tag{ 43 }$$

where

$$\begin{eqnarray}&&w(s)={s}^{{\delta }_{3}}{(1-{\alpha }_{3}s)}^{{\delta }_{4}},\end{eqnarray} \tag{ 44 }$$

with

$$\begin{eqnarray}&&{\delta }_{3}=\displaystyle \frac{1}{2}(1-{\alpha }_{1}+{\delta }_{1}),\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{\delta }_{4}=\displaystyle \frac{1}{2}\left({\alpha }_{1}-\displaystyle \frac{{\alpha }_{2}}{{\alpha }_{3}}+1\right)+\displaystyle \frac{{\delta }_{2}}{2}.\end{eqnarray} \tag{ 46 }$$

In the limit ${\alpha }_{3}=0$, it is possible to use the following identities

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\alpha }_{3}\to 0}{P}_{n}^{({\delta }_{1},{\delta }_{2})}(1-2{\alpha }_{3}s)={L}_{n}^{{\delta }_{1}}(\sqrt{{\alpha }_{2}^{2}+4{\zeta }_{1}}s),\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\alpha }_{3}\to 0}{(1-{\alpha }_{3}s)}^{{\delta }_{4}}=\exp \left[\left(\displaystyle \frac{{\alpha }_{2}}{2}-\sqrt{\displaystyle \frac{{\alpha }_{2}^{2}}{4}+{\zeta }_{1}}\right)s\right],\end{eqnarray} \tag{ 48 }$$

where L_n^a are the associated Laguerre polynomials defined in [62]. Then, we obtain the solution for equation (43) in the form

$$\begin{eqnarray}&&{\rm{\Psi }}(s)=\ {s}^{{\delta }_{3}}\exp \left[\left(\displaystyle \frac{{\alpha }_{2}}{2}-\sqrt{\displaystyle \frac{{\alpha }_{2}^{2}}{4}+{\zeta }_{1}}\right)s\right]{L}_{n}^{{\delta }_{1}}\phantom{\rule{}{2.25ex}}(\sqrt{{\alpha }_{2}^{2}+4{\zeta }_{1}}s\phantom{\rule{}{2.25ex}}).\end{eqnarray} \tag{ 49 }$$

The same procedure can also be used for the evaluation of a second solution that corresponds to the second of equation (34), which is related to

$$\begin{eqnarray}&&k={k}_{2}={\beta }_{1}+{\beta }_{2}^{1/2}.\end{eqnarray} \tag{ 50 }$$

In this way, we easily have

$$\begin{eqnarray}&&{y}_{n}={P}_{n}^{(-{\delta }_{1},{\delta }_{2})}(1-2{\alpha }_{3}s),\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&w(s)={s}^{{\tilde{\delta }}_{3}}{(1-{\alpha }_{3}s)}^{{\delta }_{4}},\end{eqnarray} \tag{ 52 }$$

where

$$\begin{eqnarray}&&{\tilde{\delta }}_{3}=\displaystyle \frac{1}{2}(1-{\alpha }_{1}-{\delta }_{1}).\end{eqnarray} \tag{ 53 }$$

By comparing equation (33) with equation (22), we find

$$\begin{eqnarray}&&{\alpha }_{1}=2,\,\,\,\,\,{\zeta }_{1}={K}_{n}^{2},\,\,\,\,\,{\zeta }_{2}=\varepsilon ,\,\,\,\,\,{\zeta }_{3}=l(l+1),\end{eqnarray} \tag{ 54 }$$

while ${\alpha }_{2}={\alpha }_{3}=0$. From this, we confirm that the solutions of equation (22) can be written as

$$\begin{eqnarray}&&{R}_{{nl}}(r)={N}_{n}{r}^{\tfrac{1}{2}(\sqrt{1+4l(l+1)}-1)}{\text{}}{e}^{-{K}_{n}r}{\text{}}{L}_{n}^{\sqrt{1+4l(l+1)}}(2{K}_{n}r),\end{eqnarray} \tag{ 55 }$$

where N_n are normalization constants. Notice that this kind of solutions is well known in the literature (e.g., see [31]). These are the radial-part solutions that correspond to bound states associated with the invariant I₀. Using the normalization condition ${\int }_{0}^{\infty }{R}_{{nl}}^{* }(r){R}_{{nl}}(r){r}^{2}{\rm{d}}{r}=1$, we have

$$\begin{eqnarray}&&{N}_{n}={\left(\displaystyle \frac{{\rm{\Gamma }}(n+1)}{2{\rm{\Gamma }}(n+2\upsilon +1)}\right)}^{1/2}{\left(\displaystyle \frac{{D}_{0}}{2{{\hslash }}^{2}{A}_{0}}\right)}^{\upsilon +1}\displaystyle \frac{1}{{\left(n+\upsilon +\tfrac{1}{2}\right)}^{\upsilon +3/2}},\end{eqnarray} \tag{ 56 }$$

where $2\upsilon =\sqrt{1+4l(l+1)}$. Note that equation (55) can now be expressed as

$$\begin{eqnarray}&&{R}_{{nl}}(r)={N}_{n}{r}^{\upsilon -1/2}\exp \left(\displaystyle \frac{-{D}_{0}}{4{{\hslash }}^{2}{A}_{0}(n+\upsilon +\tfrac{1}{2})}r\right){L}_{n}^{2\upsilon }\left(\displaystyle \frac{{D}_{0}}{2{{\hslash }}^{2}{A}_{0}(n+\upsilon +\tfrac{1}{2})}r\right),\end{eqnarray} \tag{ 57 }$$

while the eigenvalues E_nl are of the form

$$\begin{eqnarray}&&{E}_{{nl}}=-\displaystyle \frac{{D}_{0}^{2}}{4{{\hslash }}^{2}{A}_{0}}{(2n+1+\sqrt{1+4l(l+1)})}^{-2},\,\,\,\,\,n=0,1,2,\cdots .\end{eqnarray} \tag{ 58 }$$

Thus, the solutions of the radial part of the eigenvalue equation have been derived. The remaining task is to evaluate the angular part solutions.

4.1. Separation of variables

In this section, we will derive the solutions of the angular part of the equation. Angular part solutions may vary depending on the kind of the given potential. Let us consider a special case that the angular function and the potential of the system are given by

$$\begin{eqnarray}&&g(\theta ,\varphi )=\displaystyle \frac{a}{{\sin }^{2}\theta }+\displaystyle \frac{b\cos \theta }{{\sin }^{2}\theta }+\displaystyle \frac{V(\varphi )}{{\sin }^{2}\theta },\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&V(\varphi )=d\,{\csc }^{2}(\nu \varphi )+e{\sec }^{2}(\nu \varphi ),\,\,\,\,\,\nu =1\,,2,\cdots ,\end{eqnarray} \tag{ 60 }$$

where a and b are positive real parameters, while d and e are real parameters. The potential (60) is similar to the ones given in [63] and [64]. According to the separation of variables method, we take the function Y in equation (19) to be

$$\begin{eqnarray}&&{Y}_{{lm}}(\theta ,\varphi )={{\rm{\Theta }}}_{{lm}}(\theta ){{\rm{\Phi }}}_{m}(\varphi ).\end{eqnarray} \tag{ 61 }$$

By substituting this equation with equation (59) into equation (21), the differential equations for angular parts can be separated, leading to

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{\rm{\Theta }}(\theta )}{{\rm{d}}{\theta }^{2}}+\cot \theta \displaystyle \frac{{\rm{d}}{\rm{\Theta }}(\theta )}{{\rm{d}}\theta }+\left(c-\displaystyle \frac{{m}^{2}}{{\sin }^{2}\theta }-\displaystyle \frac{a}{2{\hslash }{\sin }^{2}\theta }-b\displaystyle \frac{\cos \theta }{2{\hslash }{\sin }^{2}\theta }\right){\rm{\Theta }}(\theta )=0,\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{\rm{\Phi }}(\varphi )}{{\rm{d}}{\varphi }^{2}}+\left({m}^{2}-\displaystyle \frac{V(\varphi )}{{{\hslash }}^{2}}\right){\rm{\Phi }}(\varphi )=0.\end{eqnarray} \tag{ 63 }$$

These two equations will be solved separately in the next subsections.

4.2. The solutions of Θ part

In terms of the new variable $z=\cos \theta$, equation (62) can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{\rm{\Theta }}(z)}{{{\rm{d}}{z}}^{2}}-\displaystyle \frac{2z}{1-{z}^{2}}\displaystyle \frac{{\rm{d}}{\rm{\Theta }}(z)}{{\rm{d}}{z}}+\displaystyle \frac{1}{{(1-{z}^{2})}^{2}}\left(-{{cz}}^{2}-\displaystyle \frac{b}{{{\hslash }}^{2}}z+c-{m}^{2}-\displaystyle \frac{a}{{{\hslash }}^{2}}\right){\rm{\Theta }}(z)=0.\end{eqnarray} \tag{ 64 }$$

To solve this equation, we can use a method similar to the one given in [39] for the case of the Makarov potential [65]. The NU method, which we have used in the previous section, will also be used again at this stage. By comparing equation (24) with equation (64), we have

$$\begin{eqnarray}&&\tilde{\tau }=-2z,\,\,\,\,\,\sigma =1-{z}^{2},\,\,\,\,\,\tilde{\sigma }=-{{cz}}^{2}-\displaystyle \frac{b}{{{\hslash }}^{2}}z+c-{m}^{2}-\displaystyle \frac{a}{{{\hslash }}^{2}}.\end{eqnarray} \tag{ 65 }$$

This equation can be solved through a method . The formula of π is crucial in developing a theory associated with the NU method for a specific case. In this case, it is represented as

$$\begin{eqnarray}&&\pi =\pm {\left[(c-k){z}^{2}+\displaystyle \frac{b}{{{\hslash }}^{2}}z-\left(c-{m}^{2}-\displaystyle \frac{a}{{{\hslash }}^{2}}-k\right)\right]}^{1/2}.\end{eqnarray} \tag{ 66 }$$

Notice, according to the standard NU method, that the function under the square root should be a square of some polynomials. From this fact, π is given by

$$\begin{eqnarray}&&\pi =\pm {\left(\displaystyle \frac{{m}^{2}+u+a/{{\hslash }}^{2}}{2}\right)}^{1/2}z+\displaystyle \frac{b/{{\hslash }}^{2}}{{m}^{2}+u+a/{{\hslash }}^{2}},\,\,\,\,\,\mathrm{for}\,\,k=\displaystyle \frac{2c-{m}^{2}-u-a/{{\hslash }}^{2}}{2},\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&\pi =\pm {\left(\displaystyle \frac{{m}^{2}-u+a/{{\hslash }}^{2}}{2}\right)}^{1/2}z+\displaystyle \frac{b/{{\hslash }}^{2}}{{m}^{2}-u+a/{{\hslash }}^{2}},\,\,\,\,\,\mathrm{for}\,\,k=\displaystyle \frac{2c-{m}^{2}+u-a/{{\hslash }}^{2}}{2},\end{eqnarray} \tag{ 68 }$$

where $u={[{({m}^{2}+a/{{\hslash }}^{2})}^{2}-{(b/{{\hslash }}^{2})}^{2}]}^{1/2}$. Regarding $\tau =\tilde{\tau }+2\pi$, while its derivative should be negative, we can confirm that τ becomes

$$\begin{eqnarray}&&\tau =-2z\left[1+{\left(\displaystyle \frac{{m}^{2}+u+a/{{\hslash }}^{2}}{2}\right)}^{1/2}\right]-\displaystyle \frac{b}{{{\hslash }}^{2}}{\left(\displaystyle \frac{{m}^{2}+u+a/{{\hslash }}^{2}}{2}\right)}^{-1/2}.\end{eqnarray} \tag{ 69 }$$

Then, the explicit formula of λ and ${\lambda }_{{n}^{{\prime} }}$ are, respectively, written as

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{2c-{m}^{2}-a/{{\hslash }}^{2}}{2}-\displaystyle \frac{u}{2}-{\left(\displaystyle \frac{{m}^{2}+u+a/{{\hslash }}^{2}}{2}\right)}^{1/2},\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{\lambda }_{{n}^{{\prime} }}=2{n}^{{\prime} }\left[1+{\left(\displaystyle \frac{{m}^{2}+u+a/{{\hslash }}^{2}}{2}\right)}^{1/2}\right]+{n}^{{\prime} }({n}^{{\prime} }-1).\end{eqnarray} \tag{ 71 }$$

By equating equation (70) with equation (71), we have

$$\begin{eqnarray}&&(2{n}^{{\prime} }+1){\left(\displaystyle \frac{{m}^{2}+u+a/{{\hslash }}^{2}}{2}\right)}^{1/2}+\displaystyle \frac{u-{m}^{2}-a/{{\hslash }}^{2}}{2}+{n}^{{\prime} }({n}^{{\prime} }+1)=c-{m}^{2}-\displaystyle \frac{a}{{{\hslash }}^{2}}.\end{eqnarray} \tag{ 72 }$$

We now see from the relation $c=l(l+1)$ that l is written in the form $l={n}^{{\prime} }+l^{\prime}$, where

$$\begin{eqnarray}&&l^{\prime} =\displaystyle \frac{1}{\sqrt{2}}{\left\{{m}^{2}+\displaystyle \frac{a}{{{\hslash }}^{2}}+{\left[{\left({m}^{2}+\displaystyle \frac{a}{{{\hslash }}^{2}}\right)}^{2}-{\left(\displaystyle \frac{b}{{{\hslash }}^{2}}\right)}^{2}\right]}^{1/2}\right\}}^{1/2}.\end{eqnarray} \tag{ 73 }$$

The substitution of this relation into equation (58) results in

$$\begin{eqnarray}&&{E}_{{{nn}}^{{\prime} }l^{\prime} }=-\displaystyle \frac{{D}_{0}^{2}}{4{{\hslash }}^{2}{A}_{0}}{\{2n+1+{[1+4(n^{\prime} +l^{\prime} )(n^{\prime} +l^{\prime} +1)]}^{1/2}\}}^{-2}.\end{eqnarray} \tag{ 74 }$$

These are the eigenvalues of the invariant operator.

Let us now evaluate an angular part solution ${\rm{\Theta }}(\theta )$ of the Schrödinger equation from equation (64). Starting from the functions σ and π, we apply the NU method. In this case, the functions necessary when developing the NU method are [39, 40]

$$\begin{eqnarray}&&w(z)={\rho }^{1/2}(z),\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&{y}_{{n}^{{\prime} }}(z)={y}_{0,{n}^{{\prime} }}\displaystyle \frac{1}{\rho (z)}\displaystyle \frac{{{\rm{d}}}^{{n}^{{\prime} }}}{{{\rm{d}}{z}}^{{n}^{{\prime} }}}[{(1-z)}^{{n}^{{\prime} }+({m}_{1}+{m}_{2})}{(1+z)}^{{n}^{{\prime} }+({m}_{1}-{m}_{2})}],\end{eqnarray} \tag{ 76 }$$

where

$$\begin{eqnarray}&&\rho (z)={(1-z)}^{{m}_{1}+{m}_{2}}{(1+z)}^{{m}_{1}-{m}_{2}},\end{eqnarray} \tag{ 77 }$$

with ${m}_{1}={[({m}^{2}+u+a/{{\hslash }}^{2})/2]}^{1/2}$ and ${m}_{2}={[({m}^{2}-u+a/{{\hslash }}^{2})/2]}^{1/2}$. The solutions ${y}_{{n}^{{\prime} }}(z)$ are expressed in terms of the Jacobi polynomials, ${P}_{{n}^{{\prime} }}^{({m}_{1}+{m}_{2},{m}_{1}-{m}_{2})}(z)$. Substituting equations (75) and (76) with equation (77) into equations (42) and (44), the corresponding solutions are given by

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{{n}^{{\prime} }}(\theta )={{\rm{\Theta }}}_{0,{n}^{{\prime} }}{(1-\cos \theta )}^{({m}_{1}+{m}_{2})/2}{(1+\cos \theta )}^{({m}_{1}-{m}_{2})/2}{P}_{{n}^{{\prime} }}^{({m}_{1}+{m}_{2},{m}_{1}-{m}_{2})}(\cos \theta ),\end{eqnarray} \tag{ 78 }$$

where ${{\rm{\Theta }}}_{0,{n}^{{\prime} }}$ are normalization constants. We can determine ${{\rm{\Theta }}}_{0,{n}^{{\prime} }}$ from the normalization condition, $-{\int }_{-1}^{+1}{{\rm{\Theta }}}_{{n}^{{\prime} }}^{* }(z){{\rm{\Theta }}}_{{n}^{{\prime} }}(z){\rm{d}}{z}=1$:

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{0,{n}^{{\prime} }}={\left(\displaystyle \frac{(2{n}^{{\prime} }+2{m}_{1}+1){\rm{\Gamma }}({n}^{{\prime} }+1){\rm{\Gamma }}({n}^{{\prime} }+2{m}_{1}+1)}{{2}^{2{m}_{1}+1}{\rm{\Gamma }}({n}^{{\prime} }+{m}_{1}+{m}_{2}+1){\rm{\Gamma }}({n}^{{\prime} }+{m}_{1}-{m}_{2}+1)}\right)}^{1/2}.\end{eqnarray} \tag{ 79 }$$

The relation for the orthogonality of the Jacob polynomials has been used in this evaluation.

4.3. The solutions of Φ part

It is necessary to obtain another angular part solution, equation (63), for the Schrödinger equation. After some rearrangements, equation (63) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{\rm{\Phi }}(\varphi )}{{\rm{d}}{\varphi }^{2}}+\left({m}^{2}-\displaystyle \frac{{\rm{d}}\,{\csc }^{2}(\nu \varphi )+e\,{\sec }^{2}(\nu \varphi )}{{{\hslash }}^{2}}\right){\rm{\Phi }}(\varphi )=0,\,\,\,\,\,\nu =1\,,2,\cdots .\end{eqnarray} \tag{ 80 }$$

In terms of a new variable $\tilde{z}=\cos (\nu \varphi )$, this becomes

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{\rm{\Phi }}(\tilde{z})}{{\rm{d}}{\tilde{z}}^{2}}-\displaystyle \frac{\tilde{z}}{1-{\tilde{z}}^{2}}\displaystyle \frac{{\rm{d}}{\rm{\Phi }}(\tilde{z})}{{\rm{d}}\tilde{z}}+\displaystyle \frac{-{\epsilon }^{2}{\tilde{z}}^{4}+{\rm{\Omega }}{\tilde{z}}^{2}-{\rm{\Lambda }}}{{[\tilde{z}(1-{\tilde{z}}^{2})]}^{2}}{\rm{\Phi }}(\tilde{z})=0,\end{eqnarray} \tag{ 81 }$$

where

$$\begin{eqnarray}&&{\epsilon }^{2}=\displaystyle \frac{{m}^{2}}{{\nu }^{2}},\,\,\,\,\,{\rm{\Omega }}=\displaystyle \frac{{m}^{2}}{{\nu }^{2}}-\displaystyle \frac{d}{{{\hslash }}^{2}{\nu }^{2}}+\displaystyle \frac{e}{{{\hslash }}^{2}{\nu }^{2}},\,\,\,\,\,{\rm{\Lambda }}=\displaystyle \frac{e}{{{\hslash }}^{2}{\nu }^{2}}.\end{eqnarray} \tag{ 82 }$$

To solve equation (81), we use the AIM [64, 66–69]. For this, let us now see the basic concept of the AIM. The theory of AIM is originally proposed to treat homogenous linear second order differential equations of the form

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{y}_{n}(x)}{{{\rm{d}}{x}}^{2}}={\lambda }_{0}(x)\displaystyle \frac{{{\rm{d}}{y}}_{n}(x)}{{\rm{d}}{x}}+{S}_{0}(x){y}_{n}(x),\end{eqnarray} \tag{ 83 }$$

where ${S}_{0}(x)$ and ${\lambda }_{0}(x)$ are differentiable functions with the condition ${\lambda }_{0}(x)\ne 0$. By iterating the derivatives in equation (83) up to jth number where $j=0,1\,,2,\cdots$, we obtain the equations [66]

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{j+2}{y}_{n}(x)}{{{\rm{d}}{x}}^{j+2}}={\lambda }_{j}(x)\displaystyle \frac{{{\rm{d}}{y}}_{n}(x)}{{\rm{d}}{x}}+{S}_{j}(x){y}_{n}(x).\end{eqnarray} \tag{ 84 }$$

For j = 0, this recovers to equation (83). In case ${\lambda }_{0}(x)$ and ${S}_{0}(x)$ are exactly known, ${\lambda }_{j}(x)$ and S_j(x) with $j\geqslant 1$ are given by the recurrence relations [64]

$$\begin{eqnarray}&&{\lambda }_{j}(x)={\lambda }_{0}(x){\lambda }_{j-1}(x)+{S}_{j-1}(x)+\displaystyle \frac{{\rm{d}}{\lambda }_{j-1}(x)}{{\rm{d}}{x}},\end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray}&&{S}_{j}(x)={S}_{0}(x){\lambda }_{j-1}(x)+\displaystyle \frac{{{\rm{d}}{S}}_{j-1}(x)}{{\rm{d}}{x}}.\end{eqnarray} \tag{ 86 }$$

Then, from the condition for the quantization of the system, which is [64]

$$\begin{eqnarray}{{\rm{\Delta }}}_{j}(x)\equiv \left|\begin{array}{cc}{\lambda }_{j+1}(x) & {S}_{j+1}(x)\\ {\lambda }_{j}(x) & {S}_{j}(x)\end{array}\right|=0,\end{eqnarray} \tag{ 87 }$$

the solutions are derived to be [66, 67]

$$\begin{eqnarray}&&{y}_{n}(x)=\exp \left(-{\int }^{x}\varpi (\tilde{z}){\rm{d}}{z}\right)\left[{h}_{1}+{h}_{2}{\int }^{x}\exp \left({\int }^{z}[{\lambda }_{0}(t)+2\varpi (t)]{\rm{d}}{t}\right){\rm{d}}{z}\right],\end{eqnarray} \tag{ 88 }$$

where h₁ and h₂ are constants, and $\varpi (x)$ is given by

$$\begin{eqnarray}&&\varpi (x)=\displaystyle \frac{{S}_{j}(x)}{{\lambda }_{j}(x)}=\displaystyle \frac{{S}_{j-1}(x)}{{\lambda }_{j-1}(x)}\cdots =\displaystyle \frac{{S}_{0}(x)}{{\lambda }_{0}(x)}.\end{eqnarray} \tag{ 89 }$$

Notice that we have defined ${{\rm{\Delta }}}_{j}(x)$ (equation (87)) slightly differently from that of [64] without loss of generality for our convenience.

The general solutions of equation (83) are given by equation (88). The first part of equation (88) gives the polynomial solutions that are convergent and physical, whereas the second part of equation (88) gives non-physical solutions that are divergent. Although equation (88) is the general solutions of equation (83), we take the coefficient of the second part as ${h}_{2}=0$, in order to obtain square integrable solutions. Therefore, the corresponding eigenfunctions can be derived from the following wave function generators for exactly solvable potentials

$$\begin{eqnarray}&&{y}_{n}(x)={h}_{1}\exp \left(-{\int }^{x}\varpi (z){\rm{d}}{z}\right).\end{eqnarray} \tag{ 90 }$$

Now, we take the ansatz of the form [64]

$$\begin{eqnarray}&&{\rm{\Phi }}(\tilde{z})={\tilde{z}}^{\zeta }{(1-{\tilde{z}}^{2})}^{\delta }\chi (\tilde{z}),\end{eqnarray} \tag{ 91 }$$

where

$$\begin{eqnarray}&&\delta =\displaystyle \frac{1}{4}+\displaystyle \frac{1}{4}{\left(1+4\displaystyle \frac{{\rm{d}}}{{{\hslash }}^{2}{\nu }^{2}}\right)}^{1/2},\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&\zeta =\displaystyle \frac{1}{2}+\displaystyle \frac{1}{2}{\left(1+4\displaystyle \frac{e}{{{\hslash }}^{2}{\nu }^{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 93 }$$

If we insert equation (91) into equation (81), we get [64]

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}\chi (\tilde{z})}{{\rm{d}}{\tilde{z}}^{2}}={\lambda }_{0}(\tilde{z})\displaystyle \frac{{\rm{d}}\chi (\tilde{z})}{{\rm{d}}\tilde{z}}+{S}_{0}(\tilde{z})\chi (\tilde{z}),\end{eqnarray} \tag{ 94 }$$

with

$$\begin{eqnarray}&&{\lambda }_{0}(\tilde{z})=\displaystyle \frac{(2\zeta +4\delta +1){\tilde{z}}^{2}-2\zeta }{\tilde{z}(1-{\tilde{z}}^{2})},\end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&{S}_{0}(\tilde{z})=\displaystyle \frac{{(2\delta +\zeta )}^{2}-{\epsilon }^{2}}{1-{\tilde{z}}^{2}}.\end{eqnarray} \tag{ 96 }$$

Now we can calculate ${{\rm{\Delta }}}_{0}(\tilde{z})$, ${{\rm{\Delta }}}_{1}(\tilde{z})$, ${{\rm{\Delta }}}_{2}(\tilde{z})$, ⋯ defined in equation (87), regarding the quantization condition that specifies ${\epsilon }_{0}$, ${\epsilon }_{1}$, ${\epsilon }_{2}$, ⋯. The results are

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{0}(\tilde{z})={\lambda }_{1}(\tilde{z}){S}_{0}(\tilde{z})-{S}_{1}(\tilde{z}){\lambda }_{0}(\tilde{z})=0\,\,\,\,\Rightarrow \,\,\,\,{\epsilon }_{0}=\pm (2\delta +\zeta ),\end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{1}(\tilde{z})={\lambda }_{2}(\tilde{z}){S}_{1}(\tilde{z})-{S}_{2}(\tilde{z}){\lambda }_{1}(\tilde{z})=0\,\,\,\,\Rightarrow \,\,\,\,{\epsilon }_{1}=\pm (2\delta +\zeta +2),\end{eqnarray} \tag{ 98 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{2}(\tilde{z})={\lambda }_{3}(\tilde{z}){S}_{2}(\tilde{z})-{S}_{3}(\tilde{z}){\lambda }_{2}(\tilde{z})=0\,\,\,\,\Rightarrow \,\,\,\,{\epsilon }_{2}=\pm (2\delta +\zeta +4),\end{eqnarray} \tag{ 99 }$$

and so on. In this way, we can have the general expression of eigenvalues as

$$\begin{eqnarray}&&{\epsilon }_{{n}^{{\prime\prime} }}=\pm (2\delta +\zeta +{2n}^{{\prime\prime} }),\,\,\,\,\,{n}^{{\prime\prime} }=0,1,2,\cdots .\end{eqnarray} \tag{ 100 }$$

Then, we have

$$\begin{eqnarray}&&{{m}}^{2}={{\nu }}^{2}{{\epsilon }}_{{n}^{{\prime\prime} }}^{2}={{\nu }}^{2}{(2\delta +\zeta +{2n}^{{\prime\prime} })}^{2}.\end{eqnarray} \tag{ 101 }$$

By substituting equation (73) with equation (101) for m into equation (74), we can easily identify the complete eigenvalues ${E}_{{nn}^{\prime} n^{\prime\prime} }$.

In terms of ${\tilde{z}}^{2}$, equation (94) can be transformed to a Gauss hypergeometric differential equation given in equation (15.10.1) of [70]. From this, we identify the solution of equation (94) to be [64]

$$\begin{eqnarray}&&\chi (\tilde{z}){=}_{2}{F}_{1}\left(-{n}^{{\prime\prime} },2\delta +\zeta +{n}^{{\prime\prime} }+\frac{1}{2};\zeta +\displaystyle \frac{1}{2};{\tilde{z}}^{2}\right).\end{eqnarray} \tag{ 102 }$$

Regarding this, the angular wave functions, i.e., the solutions of equation (80), are expressed in the form [64]

$$\begin{eqnarray}&&{\rm{\Phi }}({\varphi })={N}_{n^{\prime\prime} }{(\cos {\nu \varphi })}^{{\zeta }}{(\sin {\nu \varphi })}_{}^{2{\delta }}{}_{2}{{F}}_{1}\left(-{n}^{{\prime\prime} },2\delta +\zeta +{n}^{{\prime\prime} }+\displaystyle \frac{1}{2};\zeta +\displaystyle \frac{1}{2};{\cos }^{2}(\nu \varphi )\right),\end{eqnarray} \tag{ 103 }$$

where ${{ \mathcal N }}_{{n}^{{\prime\prime} }}$ are the normalization constants of the angular wave functions ${\rm{\Phi }}(\varphi )$. The constants can be calculated from the condition ${\int }_{0}^{2\pi }{| {\rm{\Phi }}(\varphi )| }^{2}{\rm{d}}\varphi =1.$ From the use of this condition with the integral relation including the hypergeometric function [70, 71]:

$$\begin{eqnarray}&&{\int }_{0}^{1}{z}^{\gamma -1}{(1-z)}^{s-\gamma }{[{}_{2}{F}_{1}(-n,n+s;\gamma ;z)]}^{2}{\rm{d}}{z}=\displaystyle \frac{n!}{(s+2n)}\displaystyle \frac{{[{\rm{\Gamma }}(\gamma )]}^{2}{\rm{\Gamma }}(n+s-\gamma +1)}{{\rm{\Gamma }}(n+s){\rm{\Gamma }}(s+\gamma )},\end{eqnarray} \tag{ 104 }$$

we obtain [64]

$$\begin{eqnarray}&&{{ \mathcal N }}_{{n}^{{\prime\prime} }}={\left(\displaystyle \frac{\left(2\delta +\zeta +2{n}^{{\prime\prime} }+\tfrac{1}{2}){\rm{\Gamma }}(2\delta +\zeta +{n}^{{\prime\prime} }+\tfrac{1}{2}){\rm{\Gamma }}(\delta +{n}^{{\prime\prime} }+\tfrac{1}{2}\right)}{{n}^{{\prime\prime} }!{\left[{\rm{\Gamma }}(\zeta +\tfrac{1}{2})\right]}^{2}{\rm{\Gamma }}(2\delta +{n}^{{\prime\prime} }+1)}\right)}^{1/2}.\end{eqnarray} \tag{ 105 }$$

Thus, the complete solutions for ${\rm{\Phi }}(\varphi )$ part are found and they have been represented in terms of the function ₂F₁. Because the relation between the hypergeometric function ₂F₁ and the general Heun function is known in the literature [72], the solutions ${\rm{\Phi }}(\varphi )$ (or ${\rm{\Phi }}(\tilde{z})$ in equation (81)) can also be represented in terms of the Heun function.

The full eigenstates of the invariant operator I₀ that belongs to the transformed system are given by multiplying the three individual solutions for separated parts derived in the previous sections. From equation (57), equation (78) with equation (79), and equation (103) with equation (105), we confirm that the eigenstates for the noncentral system characterized by equation (59) with equation (60) are represented in the form

$$\begin{eqnarray}{\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}^{{\prime} }(r,\theta ,\varphi ) & = & {N}_{n}{{\rm{\Theta }}}_{0,{n}^{{\prime} }}{{ \mathcal N }}_{{n}^{{\prime\prime} }}{r}^{\upsilon -1/2}\exp \left(\displaystyle \frac{-{D}_{0}}{4{{\hslash }}^{2}{A}_{0}(n+\upsilon +\tfrac{1}{2})}r\right){L}_{n}^{2\upsilon }\left(\displaystyle \frac{{D}_{0}}{2{{\hslash }}^{2}{A}_{0}(n+\upsilon +\tfrac{1}{2})}r\right)\\ & & \times \,{(1-\cos \theta )}^{({m}_{1}+{m}_{2})/2}{(1+\cos \theta )}^{({m}_{1}-{m}_{2})/2}{P}_{{n}^{{\prime} }}^{({m}_{1}+{m}_{2},{m}_{1}-{m}_{2})}(\cos \theta )\\ & & \times \,{(\cos \nu \varphi )}^{\zeta }{(\sin \nu \varphi )}^{2\delta }{}_{2}{{F}}_{1}(-{n}^{{\prime\prime} },2\delta +\zeta +{n}^{{\prime\prime} }+\,1/2;\zeta +1/2;{\cos }^{2}\nu \varphi ).\end{eqnarray} \tag{ 106 }$$

Now the complete normalized eigenfunctions of the invariant I that belongs to the original system, i.e., the untransformed system, are given from the relation

$$\begin{eqnarray}&&{\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(r,\theta ,\varphi )={U}^{-1}{\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}^{{\prime} }(r,\theta ,\varphi ).\end{eqnarray} \tag{ 107 }$$

The evaluation using equation (15) leads to

$$\begin{eqnarray}{\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(r,\theta ,\varphi ) & = & {N}_{n}{{\rm{\Theta }}}_{0,{n}^{{\prime} }}{{ \mathcal N }}_{{n}^{{\prime\prime} }}{\left(\displaystyle \frac{{A}_{0}D(t)}{{D}_{0}\eta (t)}\right)}^{\upsilon +1/2}{r}^{\upsilon -1/2}\exp \left(\displaystyle \frac{-{iB}(t)}{2{\hslash }\eta (t)}{r}^{2}\right)\\ & & \times \,\exp \left(\displaystyle \frac{-D(t)}{4\eta (t){{\hslash }}^{2}(n+\upsilon +\tfrac{1}{2})}r\right){L}_{n}^{2\upsilon }\left(\displaystyle \frac{D(t)}{2\eta (t){{\hslash }}^{2}(n+\upsilon +\tfrac{1}{2})}r\right)\\ & & \times \,{(1-\cos \theta )}^{({m}_{1}+{m}_{2})/2}{(1+\cos \theta )}^{({m}_{1}-{m}_{2})/2}{P}_{{n}^{{\prime} }}^{({m}_{1}+{m}_{2},{m}_{1}-{m}_{2})}(\cos \theta )\\ & & \times \,{(\cos \nu \varphi )}^{\zeta }{(\sin \nu \varphi )}_{}^{2\delta }{}_{2}{{F}}_{1}(-{n}^{{\prime\prime} },2\delta +\zeta +{n}^{{\prime\prime} }+1/2;\zeta +1/2;{\cos }^{2}\nu \varphi ).\end{eqnarray} \tag{ 108 }$$

In the next section, this formula will be used when we represent the wave functions of the system.

We know from the Lewis–Riesenfeld theory [49, 50] that the wave functions are represented in terms of the eigenststes ${\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(r,\theta ,\varphi )$ of the invariant operator, which are given in equation (108). Aside from these, the remaining factors that constitute the wave functions are phases ${\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)$ of the quantum wave. Hence, we now turn our attention to the derivation of ${\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)$. These phases satisfy the relation

$$\begin{eqnarray}&&{\hslash }\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)=\langle {\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)| \left({\rm{i}}{\hslash }\displaystyle \frac{\partial }{\partial t}-H\right)| {\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)\rangle .\end{eqnarray} \tag{ 109 }$$

By using the unitary relation given in equation (14), this equation can also be expressed in terms of transformed eigenstates such that

$$\begin{eqnarray}&&{\hslash }\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)=\langle {\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}^{{\prime} }(r,\theta ,\varphi )| \left(-\displaystyle \frac{{I}_{0}}{2\mu (t)\eta (t)}\right)| {\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}^{{\prime} }(r,\theta ,\varphi )\rangle .\end{eqnarray} \tag{ 110 }$$

Then, a minor evaluation leads to obtain the phases in the form

$$\begin{eqnarray}&&{\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)=-{E}_{{nn}^{\prime} n^{\prime\prime} }{\int }_{0}^{t}\displaystyle \frac{1}{2\mu ({t}^{{\prime} })\eta ({t}^{{\prime} })}{{\rm{d}}{t}}^{{\prime} },\end{eqnarray} \tag{ 111 }$$

where ${E}_{{nn}^{\prime} n^{\prime\prime} }$ are given by equation (74) with equations (73) and (101).

The exact solutions of the original Schrödinger equation, equation (3), associated with the Hamiltonian, equation (1), can now be written as

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{}_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}}(r,\theta ,\varphi ,t)={{\rm{e}}}^{i{\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)}{\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(r,\theta ,\varphi ),\end{eqnarray} \tag{ 112 }$$

where the full expressions of ${\phi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(r,\theta ,\varphi )$ and ${\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)$ are given in equations (108) and (111), respectively. These are the core results concerning quantum wave expression of the system. The wave functions, equation (112), satisfy the Schrödinger equation, although these solutions are somewhat complicated. These solutions allow us to study the various quantum characteristics of the system, such as the temporal evolution of canonical variables and other observables, spectrum of the quantum energy, the uncertainty relation, and the density of probability.

Let us see for a particular case that $\mu (t)$ is given by

$$\begin{eqnarray}&&\mu (t)={\mu }_{0}/[1+f\sin ({\omega }_{m}t)],\end{eqnarray} \tag{ 113 }$$

where ${\mu }_{0}$ and f are real constants, ${\omega }_{m}$ is a frequency. We assume that $f\ll 1$. Then Z(t) in equation (12) becomes

$$\begin{eqnarray}&&Z(t)=\displaystyle \frac{{D}_{0}}{2{A}_{0}^{1/2}\mu (t)}{\left[{A}_{0}-2{B}_{0}\displaystyle \frac{f+{\omega }_{m}t-f\cos ({\omega }_{m}t)}{{\mu }_{0}{\omega }_{m}}+\displaystyle \frac{{C}_{0}}{{\mu }_{0}^{2}{\omega }_{m}^{2}}{[f+{\omega }_{m}t-f\cos ({\omega }_{m}t)]}^{2}\right]}^{-1/2}.\end{eqnarray} \tag{ 114 }$$

From figures 1 and 2, we can confirm that the probability densities in this case somewhat disperse with time. The scale of $| {{\rm{\Psi }}}_{n,n^{\prime} n^{\prime\prime} }(r,\theta ,\varphi ,t){| }^{2}$ at a certain position r, appeared in these figures, is very small, even though we consider that the upper part in the graphs has been truncated near the origin. There are two main sets of reasons for this. At first, because the scale of r-axis in these figures is very large, we can conclude that the probability density is distributed through a relatively broad range in r. This may lower the scale of the probability density at a certain r. The second reason is that the graphs in figures are not plotted for entire angle ranges but for specific values of θ and . Because the angle part of the wave function is normalized through the entire solid angle, the probability density for a particular angle would be (very) lower than unity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Quantum features of molecular interactions governed by a time-dependent non-central potential that is composed of the Coulomb potential and the inverse quadratic potential have been analyzed. In the Hamiltonian of the system, given in equation (1), $g(\theta ,\varphi )$ guarantees the potential as being non-central. The consideration of time dependence of parameters are crucial in order to gain high precision results in analyzing molecular interactions for diverse quantum characteristics of actual molecular systems.

As you can see, the dependence of coefficients on time makes the Hamiltonian of the system a complicated form. Alternative powerful techniques for managing the system are necessary in such situation in order to solve Schrödinger solutions. Among various mathematical methods, we have used the invariant operator method together with the unitary transformation method, the NU method, and the AIM.

Using the Lewis–Riesenfeld theory, we have obtained the invariant of the system. The invariant operator has been transformed to a simple form by a unitary transformation. Quantum solutions in the transformed system were easily obtained, because the invariant operator in the transformed system is a time-independent simple one. The NU method and the AIM are used for solving eigenvalue equation of the invariant operator in the transformed system. The double ring-shaped generalized non-central time-dependent potential is considered as a particular case (see equation (59) with equation (60)). From the inverse transformation of quantum solutions in the transformed system, that are easily identified, we have obtained the complete quantum solutions in the original system (i.e., in the un-transformed system). The full wave functions represented in equation (112) allow us to know useful information concerning the quantum characteristics of the system, such as the time evolution of canonical variables and other observables, spectrum of the quantum energy, the uncertainty relation, and the density of probability.

For a long time, special attention has been paid to evaluating wave functions and their measurements because the exact knowledge of them is the starting point for understanding the quantum behavior of a given dynamical system (for example, see [43, 73], and references therein). Among various kinds of knowledge that can be attained from the wave functions, quantum phases ${\xi }_{{{nn}}^{{\prime} }{n}^{{\prime\prime} }}(t)$ of the wave packets given in equation (111) are important. Note that the interference of the system is determined by the quantum phase. Wave functions are sensitive to phase changes, especially when the characteristic dimension of a system is sufficiently small. The observation of genuine quantum phenomena, with interference patterns that are strictly determined by a wave equation, can be realized for particular molecular systems such as biomolecules [74]. These are accompanied by macroscopic phenomena such as quantum superpositions of energy eigenstates, decoherences, entanglements, coherent energy transports in molecules (or in organic molecules), and other novel effects which are mainly nonclassical. Topological quantum interference, as a fascinating mesoscopic effect, enables us to analyze thermal processes, quantum coherences, resonant tunnelings, and parity effects of tunnel splittings in magnetic molecular clusters [75]. In particular, interferences associated with the geometric phase is a promising research topic associated with a topological phenomenon that can be applied to phase-shifting techniques and the tunneling resonances [76–78].

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant No.: NRF-2016R1D1A1A09919503).