text:
# 1 Introduction

## 1 Introduction

In the past few years supersymmetric quantum mechanics based on shape invariance and intertwining concepts has manifested a significant progress -. Its technique started to influence not only the traditional branches of physics such as atomic, nuclear and high energy physics, which originally stimulated its emergence , but also the classical areas of mathematical physics and the theory of differential equations. Recently, in - the detailed investigation of factorization technique has been performed for one specific form of second order differential equations (SODE) with polynomial coefficients, admitting polynomial solutions based on known Rodrigue’s formula . We choose similar initial arguments to construct explicitly a wide class of QM potentials for 1D Schrödinger equations admitting (after separation of asymptotic behaviour of the wave function) polynomial solutions.

Our consideration and analysis in sections 2 and 3 shares a common subject with the known Natanzon papers , but differs from them by being more general because of the fact that we work with Papperitz rather than hypergeometric equation with arbitrary positions of the singular points (in complex domain). Moreover, in addition to transformations described in , where the spectral parameters are preserved, we include into consideration the cases when we change the roles of spectral parameters in original polynomial system and in the Schrödinger equation (generalizing in this way the known consideration of the Coulomb case, see the discussion in the text).

One more thing we would like to mention is that Turbiner’s approach  to generalized Bocher problem is very near to ours for exactly solvable case (though in  the possibility of the transformation to appropriate Sturm-Liouville problem is merely mentioned rather than investigated in full details). But for quasi-solvable cases our approach differs from those in  because we do not use the factorization ideas and investigate symmetry properties but concentrate on the connection of the appropriate polynomials with the corresponding family of Schrödinger equations. It allows us to explicitly define the additional relations for the operator to be zero-grading (in terms of ) and gives the way of regular construction of quasi-solvable potentials (all inside some definite family) with an arbitrary chosen number of algebraically constructed eigenstates.

The last remark we should make is that we do not intend to perform the comparison of the proposed approach with all other known methods of the construction of exactly solvable quantum potentials as this should be definitely the topic of a review paper rather than original research paper. But, in our opinion, the proposed approach will definitely influence the reviewing of the results obtained in other ways (especially for quasi-solvable cases) and further establishment of its inherent relations with other methods.

The paper is organized as follows. In the Section 2 we start with the special case leading to polynomial solutions, namely the polynomial family introduced in -. Section 3 is devoted to the explicit construction of the Schrödinger equation corresponding to that polynomial family, and to the presentation of some classification scheme as well as the discussion of its relation to the known solvable cases. In Section 3 we represent the generalization of the proposed scheme in some directions, namely the application to the partial differential equations and to higher order polynomial coefficients, and we demonstrate one non-trivial irreducible example and make some concluding remarks on the applicability of the considered approach.

## 2 One Construction of Polynomial Solutions for SODE

Practically all solvable 1D quantum problems correspond to those which can be transformed to the equation of hypergeometric type, and that gives eigenvalues of bound states by the requirement of finite hypergeometric series, thus being a polynomial of a given order . We start the consideration from the case leading precisely to that known situation, though from a little bit another point of view. Recently - it was shown that the eigenvalue problem for the operator of the form

$`\widehat{}={\displaystyle \frac{1}{W(x)}}{\displaystyle \frac{d}{dx}}\left(A(x)W(x){\displaystyle \frac{d}{dx}}\right)`$ (1)
leads to polynomial solutions with special requirements for the functions $`A,W`$. Namely, if we choose $`A(x)`$ as a polynomial of at most a second order, let us define,

$$A(x)=a_0+a_1x+a_2x^2,$$
(2)
and $`W(x)`$ as a non-negative function such that $`\frac{1}{W(x)}\frac{d}{dx}\left(A(x)W(x)\right)`$ is at most a first order polynomial

$$B(x)=b_0+b_1x,$$
(3)
then we can construct an orthogonal polynomial family being a solution for the eigenvalue problem, namely for the operator $`\widehat{}`$,

$`{\displaystyle \frac{1}{W(x)}}{\displaystyle \frac{d}{dx}}\left(A(x)W(x){\displaystyle \frac{d\mathrm{\Phi }_n(x)}{dx}}\right)+\gamma _n\mathrm{\Phi }_n(x)=0.`$ (4)
The polynomials given by the classical Rodrigue’s formula

$`\mathrm{\Phi }_n(x)={\displaystyle \frac{a_n}{W(x)}}\left({\displaystyle \frac{d}{dx}}\right)^n\left(A^n(x)W(x)\right)`$ (5)
are orthogonal with respect to the weight function W(x) on the interval $`(a,b)`$, chosen such that the following conditions hold

$$A(a)W(a)=A(b)W(b)=0.$$
(6)
This can be shown by classical consideration as exposed e.g. in . In this way we can choose the interval as inside the roots of the polynomial $`A(x)`$, if the latter one possesses the real roots, or the roots of $`W(x)`$, including infinity points, in the case of the functions tending to zero at infinity, or some combination of both alternatives. The eigenvalues $`\gamma _n`$ turn out to be given as -

$`\gamma _n=n\left({\displaystyle \frac{(A(x)W(x))^{}}{W(x)}}\right)^{}{\displaystyle \frac{n(n1)}{2}}A^{\prime \prime }(x),`$ (7)
which in case (2) and (3) is equal to expression (11).

At that point, having enough information, we can write and solve the equation for $`W(x)`$, namely

$`{\displaystyle \frac{d}{dx}}\left(A(x)W(x)\right)B(x)W(x)=0`$ (8)
As the equation is linear ODE of the first order its solution has the form, explicitly

$`W(x)={\displaystyle \frac{C}{A(x)}}\mathrm{exp}\left\{{\displaystyle \frac{B(x)dx}{A(x)}}\right\}={\displaystyle \frac{C\mathrm{exp}\left\{\frac{b_0+b_1x}{a_0+a_1x+a_2x^2}𝑑x\right\}}{a_0+a_1x+a_2x^2}},`$ (9)
and, of course, the integral we wrote is easy to calculate

$`{\displaystyle \frac{b_0+b_1x}{a_0+a_1x+a_2x^2}𝑑x}={\displaystyle \frac{b_1\mathrm{log}(a_0+a_1x+a_2x^2)}{2a_2}}`$ (10)
$`{\displaystyle \frac{\left(2a_2b_0+a_1b_1\right)\mathrm{arctan}(\frac{a_1+2a_2x}{\sqrt{a_{1}^{}{}_{}{}^{2}+4a_0a_2}})}{a_2\sqrt{a_{1}^{}{}_{}{}^{2}+4a_0a_2}}},`$
but for our purposes it is convenient to use it in the form we represent $`W(x)`$ in equation (9).

It is evident from the direct substitution that in the case (2) and (3) the eigenvalues $`\gamma _n`$ are given explicitly by

$`\gamma _n=n(b_1+a_2(n1)),`$ (11)
while the equation for the polynomials becomes

$`(a_0+a_1x+a_2x^2)\mathrm{\Phi }_n^{\prime \prime }(x)+(b_0+b_1x)\mathrm{\Phi }_n^{}(x)`$ (12)
$`n(b_1+a_2(n1))\mathrm{\Phi }_n(x)=0`$
The equation (12) is SODE with 3 singular points at the roots of $`A(x)`$ (we denote them as $`x_1,x_2`$) and infinity. It is easy to check that if the roots of $`A(x)`$ are different the equation is of hypergeometric type (generally speaking Papperitz equation ), whilst it is confluent hypergeometric when the roots are coincident.

The question naturally appearing is what sort of quantum mechanical problems could be associated with the polynomial family we described. To answer it, there is a natural way, namely we can try to implement adjusted pair of coordinate transformation and similarity transformation for the equation (12) in such a way, as to obtain the constant coefficient at second derivative and zero coefficient at first derivative. The resulting equation will be of Schrödinger type. Let us perform this program.

## 3 Transformation to the Schrödinger equation

Let us forget for a while that $`A(x)`$ and $`B(x)`$ are polynomials. We have the equation

$`A_2(x)y^{\prime \prime }(x)+A_1(x)y^{}(x)+\epsilon y(x)=0`$ (13)
with arbitrary coefficient functions $`A_1(x),A_2(x)`$. First, we make the variable change

$`x=F(u),`$
$`{\displaystyle \frac{d}{dx}}={\displaystyle \frac{1}{F^{}(u)}}{\displaystyle \frac{d}{du}},`$
$`{\displaystyle \frac{d^2}{dx^2}}={\displaystyle \frac{1}{F^2(u)}}{\displaystyle \frac{d^2}{du^2}}{\displaystyle \frac{F^{\prime \prime }(u)}{F^3(u)}}{\displaystyle \frac{d}{du}},`$ (14)
and choose the transformation in a form allowing to introduce some yet undefined but prescribed function of new coordinate $`\omega (u)`$ (which we could define later for the sake of the most convenient choice)

$$\omega ^2(u)[F^{}(u)]^2=A_2(x)$$
(15)
to get

$`\omega ^2(u)y^{\prime \prime }(u)+y^{}(u)\omega (u)\left({\displaystyle \frac{2A_1(F(u))A_2^{}(F(u))}{\sqrt{A_2(F(u))}}}+\omega ^{}(u)\right)`$ (16)
$`+\gamma y(u)=0`$
where prime means the differentiation with respect to the function’s argument. Now we implement the similarity transformation $`Y(u)=\mathrm{exp}(\chi (u))y(u)`$ and choose the function $`\chi (u)`$ in such a way as to kill the term with the first derivative, so that we must have

$`\chi {}_{}{}^{}(u)={\displaystyle \frac{1}{2\omega (u)}}({\displaystyle \frac{A_{2}^{}{}_{}{}^{}(F(u))}{2\sqrt{A_2(F(u))}}}{\displaystyle \frac{A_1(F(u))}{\sqrt{A_2(F(u))}}}\omega ^{}(u)).`$ (17)
Then, the equation is transformed into

$`\omega ^2(u)Y^{\prime \prime }(u)+Y(u)(\gamma {\displaystyle \frac{A_1(F(u))^2}{4A_2(F(u))}}+{\displaystyle \frac{A_1(F(u))A_{2}^{}{}_{}{}^{}(F(u))}{2A_2(F(u))}}`$ (18)
$`{\displaystyle \frac{3A_{2}^{}{}_{}{}^{}(F(u))^2}{16A_2(F(u))}}+{\displaystyle \frac{\omega ^{}(u)^2}{4}}+{\displaystyle \frac{A_{2}^{}{}_{}{}^{\prime \prime }(F(u))}{4}}{\displaystyle \frac{A_{1}^{}{}_{}{}^{}(F(u))}{2}}{\displaystyle \frac{\omega (u)\omega ^{\prime \prime }(u)}{2}})`$
which can be considered as Schrödinger type equation if we manage to identify and to separate some free constant parameter in it at $`Y(u)`$ (playing the role of ”energy”), after the division of both terms by $`\omega ^2(u)`$. As we see in the trivial choice of $`\omega (u)=1`$ we simply get the Schrödinger equation with the energy $`\gamma `$ and the potential $`V(u)`$ given by

$`V(u)=+{\displaystyle \frac{4A_1(F(u))^28A_1(F(u))A_2^{}(F(u))3A_2^{}(F(u))^2}{16A_2(F(u))}}`$ (19)
$`{\displaystyle \frac{A_1^{}(F(u))}{2}}{\displaystyle \frac{A_2^{\prime \prime }(F(u))}{4}}.`$
One known case when we have to choose $`\omega (x)`$ not equal to unity is the Coulomb potential as we will see below. Now we can use the fact that both $`A_2(x),A_1(x)`$ are polynomials, choosing them in accordance with (2,3) as $`A_2(x)=a_2x^2+a_1x+a_0`$ and $`A_1(x)=b_1x+b_0`$ to obtain the equation

$`Y^{\prime \prime }(u)+Y(u)[{\displaystyle \frac{a_2b_1+2\gamma }{2\omega ^2(u)}}+{\displaystyle \frac{\omega ^{}(u)^2}{4\omega ^2(u)}}{\displaystyle \frac{\omega ^{\prime \prime }(u)}{2\omega (u)}}`$
$`{\displaystyle \frac{3\left(a_1+2a_2F(u)\right)^2}{16\omega ^2(u)\left(a_0+a_1F(u)+a_2F(u)^2\right)}}`$ (20)
$`+{\displaystyle \frac{\left(a_1+2a_2F(u)\right)\left(b_0+b_1F(u)\right)}{2\omega ^2(u)\left(a_0+a_1F(u)+a_2F(u)^2\right)}}`$
$`{\displaystyle \frac{\left(b_0+b_1F(u)\right)^2}{4\omega ^2(u)\left(a_0+a_1F(u)+a_2F(u)^2\right)}}]=0`$
Now, we start to analyze first systematically the simplest case putting function $`\omega (u)`$ to be unity. Then, for the potential of the Schrödinger equation we have

$`V(u)={\displaystyle \frac{a_2}{2}}+{\displaystyle \frac{b_1}{2}}+{\displaystyle \frac{3\left(a_1+2a_2F(u)\right)^2}{16\left(a_0+a_1F(u)+a_2F(u)^2\right)}}+`$ (21)
$`{\displaystyle \frac{\left(a_1+2a_2F(u)\right)\left(b_0+b_1F(u)\right)}{2\left(a_0+a_1F(u)+a_2F(u)^2\right)}}+{\displaystyle \frac{\left(b_0+b_1F(u)\right)^2}{4\left(a_0+a_1F(u)+a_2F(u)^2\right)}}.`$
We can explicitly find the dependence $`x=F(u)`$ by solving the equation (15). Then, taking the inverse function, we get for $`x`$

$`x=F(u)=\{\begin{array}{cc}\frac{a_1\mathrm{sinh}(\frac{\sqrt{a_2}u}{2})^2+\sqrt{a_0a_2}\mathrm{sinh}(\sqrt{a_2}u)}{a_2},& a_20,𝒟0\\ \frac{a_1+(a_1+2a_2)e^{\sqrt{a_2}u}}{2a_2},& a_20,𝒟=0\\ \sqrt{a_0}u+\frac{a_1u^2}{4},& a_2=0,a_10\\ \frac{1}{\sqrt{a_0}}u,& a_2=0,a_1=0,a_00\end{array}`$
where $`𝒟=a_1^24a_0a_2`$ is the discriminant for $`A_2(x)`$.

The explicit expression for the quantum potential can be obtained after the substitution of (3) into expression (21) but for the general case the resulting formula becomes too complicated. Specifying the values of the parameters $`a_i,b_i`$ it is possible to show that almost all known solvable cases in quantum mechanics except the Coulomb potential, are inside the potential family we constructed. The latter one will be analyzed later when we try to construct new potentials by choosing a nontrivial $`\omega (u)`$ function. So let us first investigate the case $`\omega (u)=1`$.

We shall classify the different cases by the roots of the polynomial $`A_2(x)`$. We have two topologically different cases when $`A_2(x)`$ is not degenerate, namely when $`A_2(x)`$ has two different roots and the discriminant $`𝒟0`$, and the case of one root with degeneracy 2, $`𝒟=0`$. As one can easily see in the space of the parameters $`a_i`$ the first case fills the region inside and outside the conical surface for which the equation is $`𝒟=0`$. So we call the appropriate cases regular and irregular (for $`𝒟=0`$), respectively. We shall refer to them sometimes as Jacobi and Morse cases, based on the name of the appropriate polynomials (for the first one) and solutions (second one).

The additional cases appear in (3) as a result of the degeneracy of the polynomial $`A_2(x)`$, it could be of the first order (we call this case Laguerre’s case; $`a_2=0`$) and of the zero order (further referred to as Hermite case; $`a_2=a_1=0`$).

The first remark we would like to make is that, our polynomial family has five parameters whereas spectra depend on top power coefficients of $`A_2(x),A_1(x)`$ only, namely $`a_2,b_1`$ (see equation (11)), and so we have the evident freedom of choosing some parameters without loosing characteristic features of the problem. Obviously, we can change parameters $`b_0,a_1`$ simply by the trivial change of the origin in $`x`$ variable and the scale on it. Then, if we choose definite values of $`a_2,b_1`$ (one of them could be considered as a scale for energy and could be chosen e.g. as unity) we obtain a two parametric family of polynomials and one parametric family with full isospectrality property. The variation of the parameter $`a_0`$ then will lead to different non-trivial cases we mentioned.

So, we start with the regular case $`𝒟0`$. Then the change of variables $`xu`$ is given by the top line in (3), and the orthogonal polynomials have the weight function given by

$`W(x)={\displaystyle \frac{\mathrm{exp}\left\{\frac{\left(2a_2b_0a_1b_1\right)}{a_2𝒟}\mathrm{arctan}(\frac{a_1+2a_2x}{𝒟})\right\}}{(a_0+a_1x+a_2x^2)^{1\frac{b_1}{2a_2}}}}`$ (22)
As this family of polynomials has no commonly used name we will refer to it as to generalized Jacobi polynomials, ordinary Jacobi case corresponds to $`a_0=1,a_1=0,a_2=1`$,$`b_0+b_1=2p,b_0b_1=2q`$ corresponding to symmetrically chosen real roots at $`\pm 1`$ and the interval of orthogonality $`[1,1]`$.

The quantum potential has the following general form ($`z=\sqrt{a_2}u`$)

$`V(u)={\displaystyle \frac{A+B\mathrm{sinh}z+C\mathrm{cosh}z+D\mathrm{sinh}2z+E\mathrm{cosh}2z}{(2\sqrt{a_0a_2}\mathrm{cosh}z+a_1\mathrm{sinh}z))^2}}`$ (23)
with the coefficients $`A,B,C,D,E`$ expressed through the original ones as

$`A=`$ $`a_{2}^{}{}_{}{}^{2}\left(5a_{1}^{}{}_{}{}^{2}20a_0a_2+8b_{0}^{}{}_{}{}^{2}\right)+\left(3a_{1}^{}{}_{}{}^{2}4a_0a_2\right)b_{1}^{}{}_{}{}^{2}`$
$`+\mathrm{\hspace{0.33em}2}a_2\left(3a_{1}^{}{}_{}{}^{2}+12a_0a_24a_1b_0\right)b_1,`$
$`B=`$ $`8\sqrt{a_0a_2}\left(b_12a_2\right)\left(2a_2b_0a_1b_1\right),`$
$`C=`$ $`4a_1\left(2a_2b_1\right)\left(2a_2b_0a_1b_1\right),`$ (24)
$`D=`$ $`4a_1\sqrt{a_0a_2}\left(a_2b_1\right)^2,`$
$`E=`$ $`\left(a_{1}^{}{}_{}{}^{2}+4a_0a_2\right)\left(a_2b_1\right)^2.`$
It is straightforward to see that the potential family (23) includes Pöschl-Teller potentials (both ordinary and modified), Scarf-like potentials, Rosen-Morse and Manning-Rosen potentials  at appropriate choice of the parameters.

In singular case $`D=0`$ it is more convenient to introduce other parameters rather than $`a_i`$, namely $`a_0=\alpha ^2,a_2=\beta ^2,a_1=2\alpha \beta `$, automatically satisfying the degeneracy condition, and then the weight function becomes

$`W(x)=\left(\alpha +\beta x\right)^{2+b_1}\mathrm{exp}\left\{{\displaystyle \frac{b_1\alpha b_0\beta }{\beta ^2(\alpha +\beta x)}}\right\},`$ (25)
and the potential reads through newly introduced coefficient $`A,B,C`$

$`V(u)=A+Be^{\beta u}+Ce^{2\beta u},`$
$`A={\displaystyle \frac{\left(b_1\beta ^2\right)^2}{4\beta ^2}},`$
$`B={\displaystyle \frac{\left(b_1\alpha b_0\beta \right)\left(b_12\beta ^2\right)}{2\alpha \beta ^2}}`$ (26)
$`C={\displaystyle \frac{\left(b_1\alpha b_0\beta \right)^2}{4\alpha ^2\beta ^2}}.`$
This evidently corresponds to the Morse class of potentials .

The case when $`A_2(x)`$ becomes the first order polynomial ($`a_2=0`$), gives for the weight the following formula

$`W(x)=e^{\frac{b_1x}{a_1}}\left(a_0+a_1x\right)^{1+\frac{b_0}{a_1}\frac{a_0b_1}{a_1^2}},`$ (27)
and for the potential

$`V(u)=A+Bu^2+Cu^2,`$ (28)
with the new parameters $`A,B,C`$ (not to be confused with those obtained above) expressed through the old ones as

$`\begin{array}{cc}A& =\frac{b_1(a_1b_0a_0b_1)}{2a_1^2},\\ B& =\frac{(a_1^22a_1b_0+2a_0b_1)(3a_1^22a_1b_0+2a_0b_1)}{4a_1^4},\\ C& =\frac{b_1^2}{16}.\end{array}`$
The resulting potentials, as we see, are the combination of harmonic oscillator (HO) potential plus centrifugal potential $`B/u^2`$. And for the sake of completeness it is worthwhile to mention, that the case when $`A_2(x)`$ is a constant ($`a_2=a_1=0`$) corresponds to ordinary HO case with the oscillator position shifted by $`b_0/b_1`$.

Before going further, let us construct the explicit representation for the wave functions of the appropriate Schrödinger equation and let us discuss the bound states within this approach. As we made two subsequent transformations to obtain the Schrödinger equation, the solution in terms of polynomials has the form

$`Y(u)=\mathrm{exp}\left\{\chi (u)\right\}P_n(F(u))`$ (29)
where $`\chi (u)`$ is given by (17) and $`F(u)`$ is given by (15). The energy corresponding to this eigenfunction turns out to be the sum of $`\gamma _n`$ ( see equation (11) ) and some constant factor depending upon the parameters $`a_i,b_i`$ and leading to the shift of the energy’s origin. As the family of orthogonal polynomials has infinity and countable number of members, though the class of the potentials includes not only those which grow indefinitely at infinity (supporting bound states only), but also such with a finite number of levels (and finite ionization energy), we have to understand what is the condition for a bound state in the system. Indeed, this is very simple in the discussed case, the function $`W(u)`$ gives the asymptotic behaviour of the ground state wave function and as the point transformation $`x=F(u)`$ could be non-trivial, the resulting high order polynomials $`P_n(F(u))`$ can have growing behaviour at infinity which might more than compensate that of $`W(u)`$ and thus makes $`Y(u)`$’s norm infinite. Therefore, the condition for the bound states is simply the condition of a finite norm of $`Y(u)`$. Of course, an interesting question appears whether the non-normalizable solutions of polynomial type correspond to some physically significant features of the system, e.g. to quasi-bound states (long-lived localized states) embedded in the continuum, but we will not discuss it here.

Now, we can consider other possible choices of the function $`\omega (u)`$. This is stimulated by the known sequence of transformations for the Coulomb problem , where the first step is the change of scale in a way to put ”energy” parameter into potential function with subsequent transformation of the original Schrödinger equation into hypergeometric equation.

As we can see from the equation (3) the problem is pure algebraic and there are several ways to try to obtain the free parameter which could be interpreted as energy. The first one is to choose $`\omega (u)`$ to satisfy the equations either $`\omega ^{}(u)/\omega (u)=k`$ or $`\omega ^{\prime \prime }(u)/\omega (u)=k`$, that will produce constant factor due to fractions including the appropriate ratios in equation (3). The second one takes place when $`\omega ^2(u)=A_2(F(u))^k`$, $`k>0`$, that could lead to free coefficient in potential due to cancellation of some denominators in (3). We shall not pursue these cases any further, but treat the most important case instead. Namely, if the order of polynomial $`A_2(x)`$ is less than two, then new possible cases also appear, as we will see, for $`k=1`$, which turns out to be precisely the Coulomb case that we discuss below. The last and more special case is realized when $`A_2(x)`$ has different real roots and the coefficients $`b_0,b_1`$ are chosen in such a way as to construct common divisor (of the first order) for both numerators and denominators in potential. In this case choosing $`\omega (u)`$ in the form of $`\omega (u)=(F(u)x_1)`$ we also obtain free coefficient in potential. We will not consider all the above mentioned cases in detail here but restrict ourselves to one specific choice of $`A_2(x)`$, stimulated by the Coulomb problem.

Let us assume that $`A_2(x)=x`$ and we will choose the $`\omega (u)`$ in the form $`\omega ^2(u)=F(u)^k`$. Then, if we choose $`k=1`$, it is easy to see from the equation (15) that the point canonical transformation turns out to be identity, and we get the standard Coulomb case

$`Y^{\prime \prime }(u)+\left[{\displaystyle \frac{b_1^2}{4}}+{\displaystyle \frac{b_0}{2u^2}}{\displaystyle \frac{b_0^2}{4u^2}}{\displaystyle \frac{b_0b_1}{2u}}+{\displaystyle \frac{\gamma }{u}}\right]Y(u)=0.`$ (30)
The case $`k=1`$ is also a special case here, so we have the following expression for the variable change $`x=e^u`$ and the Schrödinger equation takes the form

$`Y^{\prime \prime }(u)+[{\displaystyle \frac{(b_01)^2}{4}}{\displaystyle \frac{\left(b_0b_12\gamma \right)}{2}}e^u{\displaystyle \frac{b_1^2}{4}}e^{2u},]Y(u)=0`$ (31)
which is the case of the quantum Morse potential.

For $`k\pm 1`$ we obtain after integration of the equation (15) the following expression for the point transformation,

$`x=\left[{\displaystyle \frac{(k1)(uC_1)}{2}}\right]^{\frac{2}{1k}}.`$ (32)
The substitution of the last expression leads to a fairly complicated form of the Schrödinger equation for $`Y(u)`$, but as one can show, there are no more cases except those we mentioned, where it is possible to obtain a free parameter in the role of ”energy”.

The successful implementation of the construction of solvable Schrödinger potentials, as we already have seen, was due to the evident existence of polynomial solution of the equation (12). It is possible to find the generalization in more complicated cases, which is the topic of the following section.

## 4 Some generalizations of the approach

As we saw, the main feature of the system considered was that the second order differential operator $`\widehat{}`$ (see equations (1) and (12)) preserves the linear subspace $`_n`$ of the polynomials of order $`n`$ for all $`n`$. It was due to the special adjustment of the orders of polynomial coefficients with the order of appropriate derivatives. This idea, of course, can be implemented not only for a special case of SODE and polynomial coefficients up to the second order, but in much more general case of linear PDEs. Indeed, we can construct the following general form of the $`n`$-th order linear differential operator $`\widehat{}`$ with the same property, so that it preserves the space when acting in the space of the polynomials of $`m`$ variables $`\stackrel{}{x}=\{x_1,\mathrm{}x_m\}`$. The general form of the appropriate linear PDEs reads

$`\widehat{}Y(\stackrel{}{x})={\displaystyle \underset{j=0}{\overset{n}{}}}P_{j+N}(\stackrel{}{x})_{\kappa _j}Y(\stackrel{}{x})=\mathrm{\hspace{0.17em}0},`$ (33)
where $`N`$ is some non-positive integer number (when $`N<0`$ we have degenerate cases, see analogous discussion on SODE in the previous section). We introduce multi-index $`\kappa _j`$ as $`j`$-th order partial derivative over arbitrary combinations of variables in a standard way, by

$`_{\kappa _j}={\displaystyle \frac{^{i_1}}{x_1^{i_1}}}\mathrm{}{\displaystyle \frac{^{i_m}}{x_m^{i_m}}},j=i_1+\mathrm{}i_m.`$ (34)
We also say that the weight of $`\kappa _j`$ equals $`j`$ and write it as $`\mathrm{\#}\kappa _j=j`$. It is worthwhile to point out that we even do not need to demand the commutativity of the derivatives, so that the same consideration could be applied for non-commutative case (quantum groups and quantum algebras see, e.g. ). Moreover, we can consider the operators which do not preserve such finite dimensional spaces, but map one into another (with higher dimension). The latter is just the permission for $`N`$ to be positive.

In the case of standard partial derivatives, when we are looking for the solution in terms of $`k`$-th order polynomials in the ring $`F[x_1,x_2,..,x_m]`$, every term in the sum in (33) maps the argument into the space spanned by the monomials $`x^{\kappa _{N+k}}=x_1^{i_1}x_2^{i_2}\mathrm{}x_{m}^{}{}_{}{}^{i_m},\{i_1+\mathrm{}i_m=N+k\}`$. The latter space is a finite dimensional vector space and a direct sum of the spaces of symmetric homogeneous polynomials corresponding to different permutations of indices for the monomials written above. We denote the space spanned by the definite monomials of order $`k`$ as $`{}_{}{}^{(k)}𝒯_{[i_1\mathrm{}i_m]}^{}`$). Then we can write down the expression for the dimension for the image space for operator action

$`dim_k^N={\displaystyle \underset{j=0}{\overset{k+N}{}}}{\displaystyle \underset{\stackrel{i_1,\mathrm{},i_m=0}{i_1+\mathrm{}+i_m=m}}{\overset{m}{}}}dim^{(j)}𝒯_{[i_1\mathrm{}i_m]}.`$ (35)
Now, the construction of the k-th order polynomial solution

$`Y_k(\stackrel{}{x})={\displaystyle \underset{j=0}{\overset{k}{}}}C_{\kappa _j}x^{\kappa _j}`$ (36)
leads simply to the linear algebraic problem for non-trivial solution for the coefficients $`C_{\kappa _j}`$.

At this point two different cases are possible. The first one is realized if $`N0`$, that is the maximal order of polynomial coefficient is less or equal to the PDE’s order. In this case we can always satisfy the system of equations because the number of linear homogeneous equations for $`C_{\kappa _j}`$ is precisely equal to $`dim_k^0`$. Then, the non-triviality condition is the condition of zero determinant for corresponding matrix obtained from equating all coefficients at monomials of type $`x^{\kappa _i}`$ to zero, and this gives us spectral parameter for the polynomial family, namely the quantization condition imposed on the coefficient $`P_0(x)`$. Then, considering the problem over the field of complex numbers, we can always construct a polynomial family in this case for some quantized value of the coefficient $`P_0(x)`$. In the contrary, when we have the condition $`N>0`$, we are still obliged to fulfill $`dim_k^N`$ conditions but only for $`dim_k^0`$ coefficients $`C_{\kappa _j}`$. The resulting system becomes overcomplete which simply means that we can construct some separate polynomial solutions of the equation (33) for only a few levels, maybe even one, that is for definite choice of $`n`$ and, additionally for special values of some of the coefficients in the coefficient polynomials.

It is very interesting to mention here that in the 1-D case considered in the previous section for SODE, the appropriate matrix turns out to be upper-three-diagonal, with additional relationship between elements, so that its determinant for $`j`$-th order system has the following form (preserving notions for coefficient of $`A_2(x),A_1(x)`$)

$`\left(\begin{array}{ccccc}& \gamma & b_0\mathrm{\hspace{0.33em}\hspace{0.33em}\; \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\; 2}a_0& \mathrm{}& 0\\ & 0& \gamma +b_1\mathrm{\hspace{0.33em}\hspace{0.33em}\; 2}b_0& \mathrm{}& 0\\ & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ & 0& \mathrm{}& (j1)b_0+(j1)(j2)a_1& j(j1)a_0\\ & 0& \mathrm{}& \gamma +(j1)b_1+(j1)(j2)a_2& jb_0+j(j1)a_1\\ & 0& \mathrm{}& 0& \gamma +jb_1+j(j1)a_2\end{array}\right)`$
that leads to one degenerate eigenvalue (rather than $`j+1`$) for a given $`j`$-th order polynomial.

What realizations of the scheme described above could be successfully used for the construction of the solutions for Schrödinger equation except those which we demonstrated in Section 2? There are three evident but not easy ways. The first and the easiest one is the consideration of 1D problems with higher order polynomials, and we intend to demonstrate this in one example, putting away its full description as a subject of separate publication. We consider the construction of the polynomial solutions starting from the third order polynomials to give the representation of the problems emerging there. Let us have the equation of the form

$`x^3y^{\prime \prime }(x)+\alpha (x^21)y(x)+(\beta x+\gamma )y(x)=0.`$ (37)
Then, in the same manner as we did in Section 2 we map the equation to the Schrödinger equation via two subsequent point transformations $`x=4/u^2`$, and the gauge transformation

$`\chi (u)=C_1+{\displaystyle \frac{\alpha u^4}{64}}+{\displaystyle \frac{32\alpha }{2}}\mathrm{log}u.`$ (38)
The resulting equation has the form

$`Y^{\prime \prime }(u)+\left(\gamma V(u)\right)Y(u)=0,`$ (39)
$`V(u)={\displaystyle \frac{(3\alpha \alpha ^2)u^2}{8}}+{\displaystyle \frac{u^6\alpha ^2}{256}}+{\displaystyle \frac{4\alpha ^28\alpha 16\beta +3}{4u^2}}.`$ (40)
The solution will be given by the formula (29), but now we have to construct the polynomials in a non-trivial way, because the Rodrigue’s formula is no more applicable. So we start to look directly for the polynomial solution of the equation (37). Let us define $`n`$-th order polynomial as $`p_n(x)=_{i=0}^nc_ix^i`$. When we substitute this anzatz into equation (37) and equate to zero the coefficients at all orders of independent variable $`x`$, we obtain $`n+1`$ equations for $`n`$ coefficients $`c_i,i=0\mathrm{}n1,\gamma `$ (the last coefficient $`c_n`$ should be chosen in order to satisfy the standardization condition). Therefore the system is overcomplete and for a non-trivial solution we must specify some additional coefficient in a unique way or to impose one additional condition on some parameters of the system. In our case we have the only choice to add $`\beta `$ to the list of coefficients to be found. Let us find then the solution for example for $`n=1,2`$ explicitly starting from the case $`n=1`$.

Solving the set of equations for $`c_0,\beta ,\gamma `$, two solutions $`c_0^\pm =\pm \mathrm{\hspace{0.17em}1},\beta =\alpha ,\gamma ^\pm =\alpha `$ can be constructed that giving for $`Y_1(u)`$

$`Y_1^\pm (u)=e^{\frac{\alpha u^4}{64}}\left(1+{\displaystyle \frac{4}{u^2}}\right)u^{\frac{32\alpha }{2}},(\gamma =\alpha ).`$ (41)
We have to impose the finiteness condition for the norm of the solution, so that in our case $`\alpha <1/4`$. The last requirement follows from simultaneous demand on proper behaviour at infinity that leads to $`\alpha <0`$ and integrability at $`u=0`$. But as one can see, in the interval $`1/2<\alpha 1/4`$ the potential becomes repulsive as $`u0`$. The last property means that we have to impose the boundary condition $`Y(u)|_{u=0}=0`$ that leads to the restriction $`\alpha 1/2`$ on the admittable region for the variation of $`\alpha `$.

Then, both solutions will be, as one can see, bound states of the system and the solution for larger $`\gamma `$ will have zero not only at $`u=0`$ but also at $`u=2`$, which evidently corresponds to the first excited state, whereas the first one corresponds to the ground state.

In a similar way, for $`n=2`$ we obtain the condition $`\beta =2\left(1+\alpha \right)`$ and the following three solutions for $`c_0,c_1,\gamma `$

$`c_0^{(1)}={\displaystyle \frac{\alpha }{\alpha +1}},c_1^{(1)}=0\gamma ^{(1)}=0,`$ (42)
$`c_0^{(\pm )}={\displaystyle \frac{\alpha }{\alpha +2}},c_1^{(\pm )}=\pm {\displaystyle \frac{\sqrt{2\alpha \left(2\alpha +3\right)}}{\alpha +2}},\gamma ^{(\pm )}=\pm 2\sqrt{2\alpha \left(2\alpha +3\right)}.`$
Then, the appropriate eigenfunctions are given by

$`Y^{(2)}(u)=e^{\frac{\alpha u^4}{64}}u^{\frac{3}{2}\alpha }\left({\displaystyle \frac{16}{u^4}}{\displaystyle \frac{\alpha }{\alpha +1}}\right),`$ (43)
$`Y^{(\pm )}(u)=e^{\frac{\alpha u^4}{64}}u^{\frac{3}{2}\alpha }\left({\displaystyle \frac{16}{u^4}}+{\displaystyle \frac{\alpha }{\alpha +2}}\pm {\displaystyle \frac{4\sqrt{2\alpha \left(2\alpha +3\right)}}{(\alpha +2)u^2}}\right).`$
The similar consideration as we performed for $`n=1`$ shows that admittable region for the parameter $`\alpha `$ now is given by $`\alpha 5/2`$. Then the constructed eigenstates represent the ground and the first two excited states for the potential

$`V(u)={\displaystyle \frac{\alpha ^2u^6}{256}}{\displaystyle \frac{\alpha \left(\alpha 3\right)u^2}{8}}+{\displaystyle \frac{4\alpha ^2+24\alpha +35}{4u^2}}.`$ (44)
The remarkable feature of the example we considered was that we construct some eigenstates corresponding to a given potential using polynomial anzatz of a given order. The solution for those cases corresponding to the ground state and to the excited states was given by the number of the roots of the

polynomials of a given order. This is in contradistinction with the standard solvable situation where polynomial’s order is equal to the quantum number because for the classical polynomial family, the $`n`$-th order polynomial has precisely $`n`$ roots on an interval where the family is defined, unlike in the general case, where the real polynomial can have complex (nonreal) roots.

¿From the last we can make some useful conclusions. As we see, if the order of polynomial coefficient functions is different from the order of the equation, the polynomial solutions in general do not exist except for the special values of the parameters. In the example that has been considered, we have to put one additional condition on the parameter to construct non-trivial solution. Nevertheless, after appropriate restriction of the region for some parameters ($`\alpha `$ in the considered example) and by fixing the value for some other (expressing $`\beta `$ through the $`\alpha `$ in the discussed case) we were able to construct eigenfunctions of some low lying eigenstates

Now, we return to the discussion of other possible generalizations of the proposed approach.

The second admissible choice is to consider, in the same way as we did, the Schrödinger equations solvable in momentum representation. As it is easy to understand, physically interesting potentials of polynomial type with the order greater than two will correspond to higher order differential equations, so that we can e.g. ask which Schrödinger equations can be constructed, based on the polynomial solutions for the equation like this $`A_4(x)y^{\prime \prime \prime \prime }(x)+A_3(x)y^{\prime \prime \prime }+\mathrm{}=0`$? We stop the discussion of this possibility at this point leaving it for the future publications too.

The last and the most complicated case corresponds to the consideration of genuine higher dimensional problems and appropriate PDEs. The most difficult obstacle to be overcome there is the necessity to perform transformations from original equation for polynomials to the equation of the Schrödinger type. Although the theory of characteristics is applicable in this case, the resulting equation, at first glance, could be hardly interpreted in terms of the Schrödinger type. As for the latter one we must demand the existence of pure constant term included into coefficient function at zero-derivative (see the discussion after equation (18)). Nevertheless, this direction is of great importance for deeper insight into integrability and solvability problems in quantum mechanics.

## 5 Discussion

In summary, we demonstrated that one can construct explicit formulae for the family of the orthogonal polynomials depending on five parameters, and thus we can associate with them the family of isospectral potentials (isospectrality with respect to 3 free parameters) which include almost all known quantum mechanically exactly solvable potentials. Some generalization of the approach to higher dimensional equations (PDEs) as well as to the higher order ODEs has been proposed.

We may conclude with a little bit speculative thought. If one would put forward the requirement for the bound states of a quantum system (in analytical case) as a demand on polynomial type of the reduction of the wave function (that seems to be reasonable and evidently fulfilled for all up-to-date known 1D solvable cases), then the immediate conclusion follows that the proposed approach (with its generalizations discussed) includes all analytically solvable cases.

## Acknowledgements

The authors are grateful to the referees for drawing our attention to the papers of Natanzon and Turbiner and for stimulating the improvement of the manuscript. This research has been supported by the Ministry of Science and Technology of the Republic of Slovenia. GK acknowledges also the postdoctoral research grant of the Slovenian Science Foundation.