Displaced Harmonic Oscillator V ∼ min [(x + d)², (x − d)²] as a Benchmark Double-Well Quantum Model

Abstract

For the displaced harmonic double-well oscillator, the existence of exact polynomial bound states at certain displacements $d$ is revealed. The N-plets of these quasi-exactly solvable (QES) states are constructed in closed form. For non-QES states, the Schrödinger equation can still be considered “non-polynomially exactly solvable” (NES) because the exact left and right parts of the wave function (proportional to confluent hypergeometric function) just have to be matched in the origin.

1. Introduction

In the traditional realistic applications of quantum mechanics one picks up a tentative one-dimensional potential $V\left(x\right)$ in a way inspired, say, by the principle of correspondence. Subsequently, one solves the corresponding “user-friendly” (i.e., ordinary differential) bound-state Schrödinger equation

$$-\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}{\psi }_n\left(x\right)+V\left(x\right){\psi }_n\left(x\right)=E_n{\psi }_n\left(x\right),{\psi }_n\left(x\right)\in L^2\left(\mathbb{R}\right),n=0,1,\dots$$

(1)

numerically, arriving at the acceptance or rejection of the model in dependence on the agreement or disagreement between the predictions and the measurements [1].

The task becomes simplified when the needs of the quantum phenomenology are met by a function $V\left(x\right)$ which is analytic. Then it makes sense to restrict attention to the potentials which are exactly solvable (ES, [2]) or at least partially exactly alias “quasi-exactly” solvable (QES, [3]). Indeed, the “redundant” analyticity requirement may prove useful, e.g., when one treats the elementary Schrödinger Equation (1) as a mere methodical laboratory needed, say, during the study of supersymmetric systems [4].

Whenever a suitable ES or QES model is not available, people feel motivated to assign the solvability status also to some non-analytic models $V\left(x\right)$. Typically, these models are of the piecewise constant square-well type [5]. In applications, their decisive formal merit is that their bound-state wave functions ${\psi }_n\left(x\right)$ can be defined, using the matching method, in terms of elementary trigonometric functions.

In between the analytic and non-analytic extremes one finds a few parity-symmetric models $V\left(x\right)=V(-x)$ in which the analyticity is only broken in the origin. Under this assumption it was possible to assign the QES solvability status to the spatially symmetrized quartic anharmonic oscillator [6] or to the symmetrized sextic anharmonic oscillator [7].

In the non-polynomial exact solvability (NES) context the non-QES exact spectra and bound states of the parity-symmetric models proved obtainable, for a few parity-symmetric potentials, via the $x=0$ matching of their exact special-function left and right parts ${\psi }_n^{(\pm )}\left(x\right)$ of the wave functions. This enabled one to assign the NES status to the symmetrized exponential potential $V\left(x\right)\sim exp\left|x\right|$ (with ${\psi }_n^{(\pm )}\left(x\right)$ expressed in terms of Bessel functions [8,9]), to its asymptotically decreasing alternative [10] and/or to the symmetrized Morse potential $V\left(x\right)$ (i.e., to the sum $V\left(x\right)$ of the two exponentials for which ${\psi }_n^{(\pm )}\left(x\right)$ appeared proportional to a confluent hypergeometric function [11,12]).

In our present paper we intend to complement these results by showing that in some of the parity-symmetric quantum models and, in particular, for

$$V\left(x\right)=V_{\left[d\right]}^{}\left(x\right)=\left\{\begin{array}{cc}{(x-d)}^2,\hfill & x>0,\hfill \\ {(x+d)}^2,\hfill & x<0\hfill \end{array}\right.$$

(2)

in Equation (1), the QES and NES regimes can coexist. One might also add that disproportionately, the ES quantum systems with the single-well harmonic-oscillator potential $V\left(x\right)=V^{\left(HO\right)}\left(x\right)=x^2$ is discussed in virtually any textbook on quantum mechanics (cf., e.g., [13]). In contrast, the most straightforward twinned, d-displaced generalization (2) of the harmonic-oscillator model with $d\in \mathbb{R}\backslash \left\{0\right\}$ is never mentioned. In what follows, we intend to fill the gap.

1.1. Motivation in Physics

In the context of physics the above-mentioned NES constructions were motivated by the needs of molecular physics so that the non-analyticity of $V\left(x\right)$ in the origin appeared inessential. Nevertheless, even in this area of physics the analyticity of the benchmark models $V\left(x\right)$ re-emerged as essential, mainly due to the related possibility of a qualitative analysis of quantum dynamics using the mathematical tools of the classical theory of catastrophes [14].

The simplest forms of the theory (viz., the so called fold and cusp catastrophes) proved mostly needed [15]. Let us just remind the readers briefly that in this language the equilibrium of a classical dynamical system is identified with a minimum of a suitable “Lyapunov” potential $V\left(x\right)$ [16]. This enables one to classify all of the qualitatively different patterns of the bifurcation of the equilibria in terms of the parameter-dependence of the minima.

By far the most interesting and useful example of such a bifurcation pattern is called cusp catastrophe. In this case one conventionally chooses

$$V\left(x\right)=V^{\left(cusp\right)}(x,a,b)=x^4+a{x}^2+bx.$$

(3)

It is then easy to verify that such a Lyapunov function possesses a single minimum at $a>0$ but two minima (i.e., a double-well shape) at $a<0$ and at a not too large b. It is also worth noting that in the a-b plane the boundary of transition between the two regimes is a spiked, cusp-shaped curve [17]. Thus, one can use the formalism for a qualitative description of dynamics, say, of the Josephson junction between the two Bose–Einstein condensates [18].

After one turns attention to the genuine quantum forms of dynamics, the work with the classical a-b plane geometry of model (3) ceases to be satisfactory. One can expect that the classical bifurcation phenomena become smeared out by the quantization. In particular, in the double-well case, due to the tunneling, one would only have to move far away from the fine-tuned cusp singularity in order to be able to distinguish, experimentally, between the low-lying states which remain localized, in general, either near the left minimum $R_{-}$ or near the right minimum $R_{+}$ of $V^{\left(cusp\right)}(x,a,b)$.

In the latter, deep-double-well dynamical regime the tunneling becomes suppressed since the growth of the distance $R_{+}$–$R_{-}$ between the minima becomes accompanied by their increasing separation mediated by a growing and thickening barrier [19]. The approximate calculation of the bound states may be then facilitated by the use of a fairly efficient local harmonic-oscillator approximation,

$$V^{\left(cusp\right)}(R_{\pm }+\xi ,a,b)=V^{\left(cusp\right)}(R_{\pm },a,b)+{\omega }_{\pm }^2{\xi }^2+\mathcal{O}\left({\xi }^3\right).$$

More or less equivalently, one may replace $V^{\left(cusp\right)}(x,a,b)$ by a global double-well harmonic oscillator,

$$V(x,a,b)=min[{\omega }_{-}^2{(x-R_{-})}^2+D_{-},{\omega }_{+}^2{(x-R_{+})}^2+D_{+}].$$

(4)

Whenever the difference $|{\omega }_{-}-{\omega }_{+}|$ or $|V^{\left(cusp\right)}(R_{-},a,b)-V^{\left(cusp\right)}(R_{+},a,b)|$ becomes large, the low-lying states remain localized just near one of the minima. The realization of a genuine quantum effect of bifurcation may only be expected to occur when ${\omega }_{-}^2\approx {\omega }_{+}^2$ and $D_{-}\approx D_{+}$. Thus, we arrive at the spatially symmetric toy-model potential

$$V^{}(x,d,\omega )={\omega }^{2}min[{(x+d)}^2,{(x-d)}^2]$$

(5)

in which we may scale out $\omega \to 1$ and obtain the ultimate double-well version of the one-parametric potential (2) of our present interest.

It is worth adding that our asymptotically quadratic one-dimensional potential (2), viz., $V_{}(x,d,1)=x^2-2d\left|x\right|+d^2$ with $d\ne 0$ represents in fact just a single-well interaction at a negative displacement $d<0$, and it only acquires the more interesting double-well shape at a positive $d>0$ (see its sample in Figure 1 where $d=2$).

1.2. Motivation in Mathematics

At the not too large displacements d, the elementary one-dimensional potential (2) can be perceived as an example of a smooth perturbation of an ubiquitous harmonic oscillator. Still, the non-analyticity of $V\left(x\right)$ in the origin looks “suspicious”. For this reason, in spite of the relative boundedness of the perturbation, the model seems to never be recalled in the context of perturbation theory. It sounds almost like a paradox in that in the textbooks on quantum mechanics the role of the most popular illustrative example of the perturbation of harmonic oscillator is almost invariably played by the (nicely analytic) quartic harmonicity ∼$x^4$ yielding the vanishing radius of convergence of the resulting Rayleigh-Schrödinger perturbation series [20].

In our recent comment on quartic anharmonicities [6] we pointed out that the not too rational (perhaps, purely emotionally motivated) insistence of the major part of the physics community on the strict analyticity of the one-dimensional phenomenological potentials $V\left(x\right)$ did also cause problems in the monograph [3]. Indeed, the author of this very nice review of the so called quasi-exactly solvable (QES) Schrödinger equations (offering also an extensive list of further relevant references) did not imagine that besides the best known sextic-polynomial QES potential $V\left(x\right)$, the QES status can be also assigned, in a way explained in [6], to its quartic-polynomial predecessor, provided only that we admit its non-analyticity at $x=0$.

Recently we returned to the QES problem and imagined that even the quartic polynomials $V\left(x\right)=A\left|x\right|+B{x}^2+C{\left|x\right|}^3+x^4$ which were studied and assigned the QES status in the latter reference need not still represent the “first nontrivial” case. Our subsequent study of the problem resulted in the conclusion that the role of the simplest nontrivial QES quantum model can in fact be played by the spiked harmonic oscillator (2).

A concise constructive demonstration of this assertion is to be presented in the following section.

2. Polynomial Solvability

Our model (2) is important as it offers a qualitative description of the phenomena connected with the changes of the topology of the forces as characterized by the so called pitchfork bifurcation of equilibria in classical mechanics [16]. In quantum context Formula (2) interpolates between the exactly solvable limit $d\to +\infty$ (yielding the doubly degenerate spectrum of energies $E_n^{\pm }[\omega ,\infty ]=(2n+1)\omega$, $n=0,1,\dots$, where the superscript marks the parity of the related wave functions) and another exactly solvable limit $d\to 0$ providing the same, but non-degenerate, spectrum of energies $E_n[\omega ,0]=(2n+1)\omega$, $n=0,1,\dots$ with parities ${(-1)}^n$. In between the wave functions have to be derived as satisfying the physical boundary conditions in infinity as well as in the origin.

2.1. The Matching of States $\psi \left(x\right)$ in the Origin

The phenomenological appeal of the genuine quantum bifurcations often forces us to turn attention to approximate methods. Pars pro toto let us mention our recent studies [19,21] of the phenomenon of tunneling in which an analytic (in fact, polynomial) potential $V\left(x\right)$ possessing an N-plet of the local minima $V\left(R_j\right)$, $j=1,2,\dots ,N$ has been replaced, in the low-lying spectrum regime, by a locally adapted harmonic-oscillator well as sampled, at $N=2$, by Equation (4). Beyond the $N=2$ harmonic-oscillator approximation as well as beyond $N=2$, the construction of the wave-functions using the matching of the logarithmic derivatives at the $(N-1)$-plet of non-analyticities becomes purely numerical.

In this sense our present $N=2$ toy model is exceptional. At any (i.e., not necessarily just at bound-state) energy E its exact and correct wave functions ${\psi }^{(\pm )}(x,E)$ (which vanish in infinity and could be called, therefore, the Jost solutions) are proportional to certain special functions which can be considered known (see their explicit form given a bit later below). Thus, one only has to match the logarithmic derivatives of the resulting left branch ${\psi }^{(-)}(x,E)$ and right branch ${\psi }^{(+)}(x,E)$ of $\psi \left(x\right)$ in the origin,

$$\underset{x\to 0}{lim}{\psi }^{(+)}(x,E)=\underset{x\to 0}{lim}{\psi }^{(-)}(x,E),\underset{x\to 0}{lim}{\left({\psi }^{(+)}\right)}^{\prime }(x,E)=\underset{x\to 0}{lim}{\left({\psi }^{(-)}\right)}^{\prime }(x,E).$$

(6)

As long as the ultimate, properly matched wave function has to be of an even or odd parity, we may write $\psi (x,E)=\psi (x,E,p)={(-1)}^p\psi (-x,E,p)$ with $p=0$ or $p=1$, respectively. Thus, under certain non-singularity conditions in the origin we may introduce an arbitrary regular ad hoc normalization

$$\underset{x\to 0}{lim}{\psi }^{(+)}(x,E,0)={\mathcal{N}}_0\ne 0,\underset{x\to 0}{lim}{\left({\psi }_n^{(+)}\right)}^{\prime }(x,E,1)={\mathcal{N}}_1\ne 0$$

(7)

and arrive at the two alternative “shooting method” matching rules

$$\underset{x\to 0}{lim}{\left({\psi }^{(+)}\right)}^{\prime }(x,E,0)=0,\underset{x\to 0}{lim}{\psi }^{(+)}(x,E,1)=0.$$

(8)

In this spirit also the above-cited samples of a few other NES constructions of bound states can be perceived as the implementations of the same matching rules alias transcendental Equation (8). In all of these cases, these equations determine the whole spectrum of the bound state energies $E=E_n$ in an implicit manner.

2.2. The Simplest Polynomial Solution

Our interest in one-dimensional potential (2) was born when we noticed that it is quasi-exactly solvable (QES). By definition [3] this means that it admits the existence of exact elementary wave functions at certain parameters and energies. One can easily verify, by insertion, that at a special shift $d=d^{\left(QES\right)}=-1$ the nodeless function

$${\psi }_0^{\left(QES\right)}\left(x\right)=\left\{\begin{array}{c}\left(1+x\right)e^{-\left(x^2/2+x\right)},x>0,\hfill \\ \left(1-x\right)e^{-\left(x^2/2-x\right)},x<0\hfill \end{array}\right.$$

(9)

belongs to $L^2\left(\mathbb{R}\right)$, obeys the standard matching conditions (i.e., has a continuous logarithmic derivative in the origin) and satisfies our Schrödinger Equations (1) + (2) at $E=E_0^{\left(QES\right)}=3$ exactly. Thus, this expression represents the even-parity ground-state wave function of the model.

A non-constructive proof of the QES property of our model is elementary. In the single-well scenario we choose the special value of $d=d^{\left(QES\right)}=-1$ and inserted the nodeless function

$${\psi }_0^{\left(QES\right)}\left(x\right)=\left\{\begin{array}{c}\left(1+x\right)e^{-\left(x^2/2+x\right)},x>0,\hfill \\ \left(1-x\right)e^{-\left(x^2/2-x\right)},x<0\hfill \end{array}\right.$$

(10)

(which belongs to $L^2\left(\mathbb{R}\right)$) in Schrödinger equation

$$-\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\psi \left(x\right)+V\left(x\right)\psi \left(x\right)=E\psi \left(x\right).$$

(11)

Clearly, function (10) obeys the standard matching conditions in the origin (i.e., it has a continuous logarithmic derivative there) and it solves our bound-state problem (2) + (11) at $E=E_0^{\left(QES\right)}=3$. Thus, QES expression (10) represents the exact even-parity ground-state wave function. In Figure 2 our QES choice of parameters has been used for illustration purposes. The similar elementary QES feature also occurs after one moves to the double-well dynamical regime.

Let us add that in the standard literature the QES property has only been assigned to the 1D potentials which were manifestly analytic along the whole real line [22,23,24,25,26,27,28]. The present, non-analytic generalization of the QES property has only been proposed in Ref. [6]. The non-analytic QES construction has been preformed there for the phenomenologically important class of quartic anharmonic oscillators. Our present real QES wave function (10) may be also recalled as providing a closed-form illustration of the fact that besides the potential, also our wave functions are non-analytic at $x=0$ since

$$\frac{d^3}{d{x}^3}{\psi }_0^{\left(QES\right)}\left(0^{+}\right)=2\ne \frac{d^3}{d{x}^3}{\psi }_0^{\left(QES\right)}\left(0^{-}\right)=-2.$$

(12)

We are persuaded that the descriptive merits of the similar generalized QES constructions overweight their not too essential feature of non-analyticity at $x=0$. In this sense the existing literature would certainly deserve to be completed.

2.3. Systematic Approach

After one inserts potential (2) in Schrödinger Equation (11) it is possible to assume the existence of the even-parity and/or odd-parity QES bound-state solutions in the most general form of the standard normalizable polynomial ansatz such that, on the positive half-axis,

$${\psi }^{\left(QES\right)}\left(x\right)=\psi \left(x\right)=e^{-x^2/2+dx}\times \sum _{k=0}^N{a}_k{x}^k,x\ge 0,a_N\ne 0.$$

(13)

We may accept the normalization convention

$$\left\{\begin{array}{cc}\psi \left(0\right)=1,& {\psi }^{\prime }\left(0\right)=0,\\ \psi \left(0\right)=0,& {\psi }^{\prime }\left(0\right)=1,\end{array}\right.\mathrm{i}.\mathrm{e}.,\left\{\begin{array}{cccc}{a}_0=1,& a_1=-d,& & \mathrm{parity}=\mathrm{even},\\ a_0=0,& a_1=1,& a_2=-d,& \mathrm{parity}=\mathrm{odd}.\end{array}\right.$$

(14)

With $a_{N+1}=a_{N+2}=\dots =0$ and after appropriate insertions and elementary algebra the QES solvability of the model appears equivalent to the validity of the set of linear recurrences

$$(E-1-2n)a_n+2d(n+1)a_{n+1}+(n+1)(n+2)a_{n+2}=0,n=0,1,\dots ,N.$$

(15)

From the last item with $n=N$ we get the QES value of energy $E=2N+1$ so that we are left with an N-plet of relations between the $N-1$ unknown coefficients $a_2$, $a_3$, …, $a_N$ and an unknown, QES-compatible value (or rather a multiplet of values) of the shift $d=d^{\left(QES\right)}$.

Starting from the choice of $N=0$ and $E=1$ we find that the odd solution cannot exist while the well known even-parity solution only exists at $d=0$. In the next step with $N=1$ and $E=3$ the single constraint $2a_0+2d{a}_1=0$ only admits the pure $d=0$ harmonic oscillator in the odd-parity case while the first nontrivial QES solutions with $d=\pm 1$ emerges in the even-parity scenario.

Once we skip the harmonic-oscillator solutions and demand that $d\ne 0$ we may set $N=2$ and $E=5$ in Equation (15). Then we have to satisfy the set of two relations

$$4a_0+2d{a}_1+2a_2=0,2a_1+4d{a}_2=0.$$

(16)

This leads to the two easy odd-parity solutions with $d=d_{\pm }=\pm 1/\sqrt{2}$ and $a_2=-d$. They become accompanied by the equally easy even-parity solutions with $a_2=1/2$ and $d=\pm \sqrt{5/2}$.

Next we choose $N=3$ and find that the odd-parity construction degenerates to the triplet of relations

$$2+2d{\mathit{a}}_2+3{\mathit{a}}_3=0,d+{\mathit{a}}_2=0,{\mathit{a}}_2+3d{\mathit{a}}_3=0$$

with solution $a_3=1/3$, $d=\pm \sqrt{3/2}$ and $a_2=-d$. In parallel, the even-parity conditions

$$-2d+2d{\mathit{a}}_2+3{\mathit{a}}_3=0,3-d^2+{\mathit{a}}_2=0,{\mathit{a}}_2+3d{\mathit{a}}_3=0$$

lead to the slightly less trivial quadruplet of eligible QES shifts

$$d_{\pm ,\pm }=\pm \sqrt{\frac{9\pm \sqrt{57}}{4}}.$$

(17)

Each one of them defines the related two coefficients

$$a_2=\frac{2d^2}{2d^2-1},a_3=\frac{2d}{6d^2-3}.$$

(18)

2.4. General Case

At any N we may study the QES sets of coupled nonlinear algebraic equations using the computer-assisted elimination technique of Gröbner bases. In the even-parity case this algorithm generates certain polynomials $P^{\left(N\right)}\left(d\right)$ which specify the QES-compatible shifts d as their zeros (see

Table 1. QES even-parity values of shifts d are defined, implicitly, as zeros of polynomials $P^{\left(N\right)}\left(d\right)$. A sample of the Gröbnerian elimination is added.

N	Polynomial ${\mathit{P}}^{\left(\mathit{N}\right)}\left(\mathit{d}\right)$	Elimination of ${\mathit{a}}_{\mathit{N}}$
2	$-5+2d^2$	$-1+2a_2=0$
3	$3-9d^2+2d^4$	$-8d+2d^3+3a_3=0$
4	$27-28d^2+4d^4$	$11-2d^2+12a_4=0$
5	$-15+75d^2-40d^4+4d^6$	$29d-19d^3+2d^5+15a_5=0$
⋮	…	…

). The odd-parity QES solutions appear determined, similarly, via the values of $a=a_2=-d$ which coincide with the real zeros of other polynomials $Q^{\left(N\right)}\left(a\right)$ (see

Table 2. QES values $a=a_2=-d$ determined as zeros of polynomials $Q^{\left(N\right)}\left(a\right)$ (odd-parity case).

N	Polynomial ${\mathit{Q}}^{\left(\mathit{N}\right)}\left(\mathit{a}\right)$	Elimination of ${\mathit{a}}_{\mathit{N}}$
2	$-1+2a^2$	$-a+a_2=0$
3	$-3+2a^2$	$-1+3a_3=0$
4	$3-12a^2+4a^4$	$5a-2a^3+6a_4=0$
5	$15-20a^2+4a^4$	$7-2a^2+30a_5=0$
⋮	…	…

).

The most compact representation of the above set of results may be obtained when one imagines that the even- and/or odd-parity boundary condition (14) at $x=0$ may be reinterpreted as the respective initial, additional (i.e., $n=-1$) condition for recurrences (15). This enables us to solve these recurrences in terms of compact formulae.

Theorem 1.

In the even-parity QES case the general closed-form solution of recurrences (15) reads

$$a_k=a_k^{\left(N\right)}=\frac{{(-1)}^k}{k!(k-1)!}det{\mathcal{P}}^{(N,k)}\left(d\right),k=(1,2,)3,4,\dots .$$

(19)

Its d-dependence is encoded in the tridiagonal k by k matrices

$${\mathcal{P}}^{(N,k)}\left(d\right)=\left(\begin{array}{cccccc}d\hfill & 1\hfill & 0\hfill & 0\hfill & \dots \hfill & 0\hfill \\ 2N\hfill & 2d\hfill & 2\hfill & 0\hfill & \ddots \hfill & ⋮\hfill \\ 0\hfill & 2N-2\hfill & 4d\hfill & 6\hfill & \ddots \hfill & 0\hfill \\ 0\hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & 0\hfill \\ ⋮\hfill & \ddots \hfill & 0\hfill & 2(N-k+3)\hfill & 2(k-2)d\hfill & (k-1)(k-2)\hfill \\ 0\hfill & \dots \hfill & 0\hfill & 0\hfill & 2(N-k+2)\hfill & 2(k-1)d\hfill \end{array}\right).$$

(20)

These formulae must be complemented by the constraint $a_{N+1}^{\left(N\right)}=0$, i.e., by the polynomial algebraic equation

$$det{\mathcal{P}}^{(N,N+1)}\left(d\right)=0$$

(21)

which represents an implicit definition of all of the admissible values of the specific, QES-compatible shift parameters $d=d^{\left(QES\right)}$.

Proof. 

One may proceed by mathematical induction. It is only necessary to keep in mind that in matrix (20) the first row reflecting the initial condition looks (and is) anomalous. That is why the form and validity of the $k=1$ solution $a_1=-d$ and also of the $k=2$ solution $a_2=d^2-N$ are to be verified separately. □

Theorem 2.

In the odd-parity QES case the closed-form solution of recurrences (15) reads

$$a_{k+1}=a_{k+1}^{\left(N\right)}=\frac{{(-1)}^k}{(k+1)!k!}det{\mathcal{Q}}^{(N,k)}\left(d\right),k=1,2,\dots$$

(22)

with the d-dependence encoded in the tridiagonal k by k matrices

$${\mathcal{Q}}^{(N,k)}\left(d\right)=\left(\begin{array}{ccccc}2d\hfill & 2\hfill & 0\hfill & \dots \hfill & 0\hfill \\ 2N-2\hfill & 4d\hfill & 6\hfill & \ddots \hfill & ⋮\hfill \\ 0\hfill & 2N-4\hfill & \ddots \hfill & \ddots \hfill & 0\hfill \\ ⋮\hfill & \ddots \hfill & \ddots \hfill & 2(k-1)d\hfill & k(k-1)\hfill \\ 0\hfill & \dots \hfill & 0\hfill & 2(N-k+1)\hfill & 2kd\hfill \end{array}\right).$$

(23)

This must be complemented by the algebraic polynomial-equation definition $a_{N+1}^{\left(N\right)}=0$, i.e.,

$$det{\mathcal{Q}}^{(N,N)}\left(d\right)=0$$

(24)

yielding all of the admissible values of the specific, QES-compatible shift parameters $d=d^{\left(QES\right)}$.

Proof. 

Once one keeps in mind that in the odd-parity cases we have $a_2=-d$ at all N and k, the proof remains elementary and proceeds, along similar lines as above, by mathematical induction. □

3. Special-Function Solvability

In the broader context of quantum physics, the amount of attention which is currently paid to polynomially solvable (traditionally called “exactly solvable”, ES) Schrödinger equations seems disproportionate. Especially in comparison with the fact that the attention virtually disappears in the NES cases, i.e., immediately behind the latter ES area. The paradox is striking because the sophisticated commercial software (such as MAPLE or MATHEMATICA, etc.) as well as the freely available NIST review [29] enable the physicists to work with elementary functions as easily as with the special functions.

Within the NES approach one can decide to treat the terminating and non-terminating confluent hypergeometric series on equal footing. As an immediate consequence of such an amendment of the paradigm one obtains an extension of the class of the user-friendly interactions. In our present paper we accept the idea and intend to describe one of its most natural implementations in the one-dimensional quantum-dynamical scenario.

In this context one can feel inspired by the classical Thom’s concepts of stability and of its loss which are closely related to the algebraic geometry of the minima of certain analytic Lyuapunov’s functions $V\left(x\right)$ [30]. For this reason the classical dynamics near a cusp-related bifurcation (i.e., near an “abrupt” transition between one and two stable minima of $V\left(x\right)$) is traditionally described via a representative analytic, anharmonic-oscillator (AHO) polynomial (3). It is widely believed that the polynomial Lyapunov’s functions $V\left(x\right)$ are best suited for the study of the Thom’s classical catastrophes and that even the constructions of their quantum analogues should not break the correspondence with these polynomials. Thus, even though our present model $V_{\left[d\right]}^{}\left(x\right)$ could have also pretended to play the role of a cusp-related Lyapunov function (especially in the simplest nontrivial bifurcation scenario with $b=0$ in AHO model), it seemed disqualified by its non-analyticity in the origin.

3.1. Single-Well Dynamical Regime with $D<0$

Very recently we imagined that any quantization of a classical catastrophe must replace the usual schematic time-dependent trajectories $x=x\left(t\right)$ in a classical space of states by certain selfadjoint operators $\widehat{x}\left(t\right)$ in a suitable infinite-dimensional physical Hilbert space ${\mathcal{H}}^{\left(P\right)}$ [31,32,33]. Thus, after quantization the concepts such as bifurcation (etc.) must necessarily be transferred from the domain of algebraic geometry to the area of functional analysis and spectral theory. Naturally, this implies that the non-analyticity of $V_{\left[d\right]}^{}\left(x\right)$ becomes inessential while its solvability may start playing the role of a significant merit in comparison with a purely numerical nature of AHO Schrödinger equations [34,35,36,37]. In other words, the twinned harmonic oscillator (2) with a variable real displacement d and with two minima if $d>0$ seems suitable and acceptable as a benchmark quantum model with catastrophic features.

Besides the serendipitous and rather formal QES property as discussed above, a turn of attention to the fully general twinned harmonic-oscillator model (2) is recommended from a phenomenological perspective. Another independent reason would be again mathematical. After an elementary change of variables one reveals that Equation (1) coincides with the two branches of the confluent hypergeometric equation which are matched together in the origin. In this sense our model (1) + (2) resembles square wells and remains solvable via the matching of logarithmic derivatives of $\psi \left(x\right)$ at $x=0$. In the literature the conceptual appeal of the complete and, in principle, arbitrarily precise solvability based on the use of confluent hypergeometric special functions may be found advocated by multiple authors [1]. In what follows we will defend and accept such an NES point of view, demonstrating its applicability to the one-parametric family of potentials (2) in both of their one-centered and two-centered realizations.

Once we recall the left-right symmetry (i.e., parity-symmetry, $\mathcal{P}$-symmetry) of our Hamiltonian, we may temporarily restrict our attention to the differential Equation (1) on one of the half-lines (i.e., say, with $x\in (0,\infty )$). Then we may change variables as follows,

$$x-d=\sqrt{z}\in (-d,\infty ),\psi \left(x\right)=exp(-z/2)w\left(z\right).$$

Whenever we choose a negative $d=-{\mu }^2<0$, this converts Equation (1) into Kummer’s equation

$$z\frac{d^{2}w\left(z\right)}{{dz}^2}+(b-z)\frac{dw\left(z\right)}{dz}-aw\left(z\right)=0,z\in ({\mu }^4,\infty )$$

(25)

(cf. [29], §13.2(i)) with parameters $a=(1-E)/4$ and $b=1/2$. As long as $z\ge {\mu }^2>0$ we may ignore the branch-point singularity of Equation (25) at $z=0$ and we may write down the general solution as a superposition of two independent confluent hypergeometric series, viz.,

$$w_{\left(even\right)}\left(z\right)=M(a,b,z)={}_1{F}_1(a,b,z)=$$

$$=\sum _{s=0}^{\infty }\frac{{\left(a\right)}_s}{{\left(b\right)}_{s}s!}{z}^s=1+\frac{a}{b}z+\frac{a(a+1)}{b(b+1)2!}{z}^2+\cdots$$

(26)

and

$$w_{\left(odd\right)}\left(z\right)=z^{1-b}M(a-b+1,2-b,z).$$

(27)

Subsequently, we have to satisfy the asymptotic physical boundary condition. This forces us to combine the two independent solutions $w_{\left(even\right)}\left(z\right)$ and $w_{\left(odd\right)}\left(z\right)$ of Equation (25) in the following unique manner,

$$w_{\left(physical\right)}\left(z\right)=U\left(a,b,z\right)=\frac{\Gamma \left(1-b\right)}{\Gamma \left(a-b+1\right)}M\left(a,b,z\right)+$$

$$+\frac{\Gamma \left(b-1\right)}{\Gamma \left(a\right)}{z}^{1-b}M\left(a-b+1,2-b,z\right).$$

(28)

The choice of this function with the power-law asymptotical behavior $U\left(a,b,z\right)\sim z^{-a}$ valid in the complex plane of $z=x+{\mu }^2$ with the left cut from $z=0$ to $-\infty$ yields the correct, asymptotically vanishing wave functions at positive $x\ge 0$ in closed form,

$$\psi \left(x\right)=\psi (x,E)=$$

$$=exp(-{(x+{\mu }^2)}^2/2)U\left((1-E)/4,1/2,{(x+{\mu }^2)}^2\right).$$

(29)

The formula remains manifestly dependent on the value of energy E. After the matching of the left and right branches of our analytic solution in the origin, the spectrum of bound states will be obtained as split in the even-parity and the odd-parity states. Their energies $E_n^{(\pm )}$ will coincide with the roots of the respective transcendental secular equations, viz.,

$$\psi (0,E^{(+)})=1{\psi }^{\prime }(0,E^{(+)})=0$$

(30)

or

$$\psi (0,E^{(-)})=0{\psi }^{\prime }(0,E^{(-)})=1.$$

(31)

Here the prime marks the derivative with respect to the first argument.

The numerical performance of the determination of the energies based on the latter equations is entirely routine and reflects our expectations that with the growth of the negative shift ${\mu }^2=-d$ the narrowing of the potential well will push the energy levels upwards. Thus, the standard HO ground state with the energy equal to one will be pushed up to $E_0=E_0(-{\mu }^2)$ with $E_0(-0.25)\approx 1.3349$ and $E_0(-1)=3$ (this is quasi-exact state), etc. The similar accelerated growth moves the second excited state from initial $E_2\left(0\right)=5$ more quickly to $E_2(-1)\in (8.658,8.659)$, etc. The bracketing property of the latter type of results is highly welcome in computations of course.

3.2. Double-Well Dynamical Regime with $D>0$

In the double-well regime with positive $d={\nu }^2>0$ the presence of the square-root branch-point of Equation (25) at $x={\nu }^2$, i.e., inside the interval of physical coordinates would make the direct use of the asymptotically correct closed-form solutions (28) slightly uncomfortable. The key new technicality is that one must keep the choice of the proper branch of the coordinate function $x=x\left(z\right)=d\pm \sqrt{z}$ under explicit control.

For this reason we recommend the use of an alternative bound-state construction in which one recalls the $\mathcal{P}$-symmetry of our Schrödinger Equation (1) and employs a variant of the so called shooting method. In the language of mathematics this means that one simply separates the spatially even wave functions (for which we set $\psi \left(0\right)=1$ and ${\psi }^{\prime }\left(0\right)=0$ as an initial boundary condition) from the spatially odd solutions (for which we have $\psi \left(0\right)=0$ and ${\psi }^{\prime }\left(0\right)=1$)—this is to be compared with Equations (30) and (31), respectively. In each of these two subcases one arrives again at the Kummer’s Equation (25) which is now defined in an extended and partially doubly covered interval of $z\in (0,\infty )$ and $z\in (0,{\nu }^4)$.

In contrast to preceding subsection we will now solve the initial-value alias Cauchy problem knowing that we know its general solution in confluent-hypergeometric form

$$\psi \left(x\right)=\psi (x,E)=C_1{w}_{\left(even\right)}\left[z\left(x\right)\right]+C_2{w}_{\left(odd\right)}\left[z\left(x\right)\right].$$

(32)

In the next stage of development we fix the ratio of the coefficients $C_j$ via the above-mentioned initial conditions imposed upon the both-sided limits of $\psi \left(x\right)$ and ${\psi }^{\prime }\left(x\right)$ at $x=0$ (i.e., at $z={\nu }^4$ when taken below its square-root branch). The subsequent necessary control of the choice of the branch of function (27) (i.e., of the analytic-function component $w_{\left(odd\right)}\left[z\left(x\right)\right]$ of our solution (32) at $z=0$) is not difficult because we have $b=1/2$ (i.e., a square-root branch) in definition (27).

In the last step we may finally recall the standard oscillation theorems [38] which

• relate the quantum number $n=0,1,\dots$ of a bound-state energy $E_n$ to the number $Z_n=n$ of the finite nodal zeros $x_j\in (-\infty ,\infty )$ of ${\psi }_n\left(x\right)$;
• enable us to perturb $E_n\to E_n^{(+)}=E_n+{\delta }^2$, localize the related maximal real (and, say, positive) finite nodal zero $x_{\left(max\right)}^{(+)}>0$ and find the one-to-one correspondence between the exact-bound-state-energy limits of $\delta \to 0$ and of $x_{\left(max\right)}^{(+)}\to \infty$.

The latter correspondence enables us to localize, numerically, the exact energy, with arbitrary preselected precision, via a trial-and-error choice of $E_n^{(\pm )}$, to be iteratively amended by a tentative, smaller and smaller decrease (whenever the maximal real node $x_{\left(max\right)}^{(+)}$ keeps growing) or increase (whenever the maximal node happens to leave the real line and become complex).

4. The Displacement-Dependence of the Spectrum

The numerical results of the above-outlined algorithm are sampled in Figure 3 and by

Table 3. The shift-dependence of the ground-state energy.

Shift	Energy Estimates
$d$	Lower	Upper
0	1	1
0.25	0.768972	0.768974
0.50	0.635528	0.635530
0.75	0.590300	0.590301
1.00	0.618910	0.618920
1.50	0.801493	0.801494
2.00	0.951410	0.951420
∞	1	1

,

Table 4. The shift-dependence of the 1st-excited-state energy.

Shift	Energy Estimates
$d$	Lower	Upper
0	3	3
0.25	2.483910	2.483920
0.50	2.060760	2.060770
0.75	1.724710	1.724720
1.00	1.468460	1.468470
1.50	1.157479	1.157480
2.00	1.035760	1.035770
∞	1	1

and

Table 5. Shift-dependence of the 2nd-excited-state energy (at $d=1$ the state is quasi-exact).

Shift	Energy Estimates
$d$	Lower	Upper
0	5	5
0.25	4.34600	4.34700
0.50	3.79410	3.79420
0.75	3.34470	3.34471
1	3	3
1.50	2.64860	2.64870
2.00	2.73500	2.73510
∞	3	3

. The key advantage of such a special-function-based version of the well known shooting algorithm is that in our model the trial-and-error functions $\psi \left(x\right)$ are known in the exact (albeit rather clumsy, computer-stored) special-function confluent-hypergeometric form (which is closed and analytic but, unfortunately, still too long for a printed display). This enables us to keep the numerical errors under active control and to construct, systematically, the pairs of the lower and upper estimates of the bound-state energies.

In Table 3 sampling the ground-state results, the non-monotonicity of the $d-$dependence of the energy is slightly counterintuitive but still easy to understand. It is obvious that the growth of the positive shift $d>0$ leads, at the beginning, to the broadening of the well without a significant change of its size near the origin. This implies the possibility of an approximation of the potential by a broader effective well, i.e., in a solvable approximation, by $V^{\left(eff\right)}={\omega }^2{x}^2$ with $\omega <1$ yielding the smaller ground state energy. Later, i.e., at the larger shifts d the central repulsive spike starts growing more quickly, separating the two HO wells and making them only weakly coupled. This must lead, approximately, to the double degeneracy of the low-lying spectrum and, in particular, to the unavoidable return of the ground state energy to its initial value at the sufficiently large ds.

Table 4 illustrates that as long as the wave function of the first excited state vanishes in the origin, it is much less sensitive to the influence of the central spike. According to our tests, the energy just decreases with the growth of d and, in contrast to the behavior of the ground state energy, it never drops below $E=1$ at the large shifts $d\gg 1$. Marginally we may add that the test is facilitated by the observation that at the relevant spectral bound $E=1$, the wave-function solutions with the correct odd-parity behavior in the origin may be further simplified and expressed in terms of error functions. This simplification of the general solution also facilitates the demonstration of the complex-value nature of the asymptotic nodal zero at any $E<1$.

Table 5 shows that the d-dependence of the second excited level (i.e., its decrease followed by the increase) is similar to its ground-state predecessor. Figure 3 complements these results by a typical sample of the first three wave functions at $d=3/2$.

5. Conclusions

In conclusion, let us re-emphasize that in the current literature there exists a large number of papers and studies devoted to the experiments and realistic models of quantum physical systems described by a simplified, one-dimensional (i.e., ordinary differential) Schrödinger equation (i.e., often, by our Equation (1) which is defined on the whole real line of the coordinate). Still, in contrast to the dominance of the scene by the papers describing the bound states lying in a single-well potential, one finds perceivably less papers dealing with the potentials having a double-well shape. Naturally, one of the reasons is purely formal. One often finds that many methods of solving Schrödinger equations (and, among them, the perturbation-expansion methods) happen to fail when applied to the double-well systems.

This leads to a remarkable conflict between the prevailing, purely mathematically motivated preference of the single-well models in various methodical considerations, and the frequently encountered lack of some sufficiently efficient construction recipes, especially in the context of the double-well scenarios which are much more interesting, in the phenomenological context, mainly due to the emerging manifestations of the influence of the tunneling effect.

In such a context we sought, in our present paper, for a nontrivial support of the existing tendencies towards an extension of the class of the “exactly solvable” 1D potentials $V\left(x\right)$ from their completely analytic subset to its non-analytic generalizations. We are persuaded that such a tendency is natural and well motivated. We also believe that in the nearest future our approach might find some immediate and successful applications, e.g., during some of the popular and, usually, partially non-numerical descriptions of the physical systems (as studied, e.g., in [39]) in which, in the language of mathematics, the double-well structure of the potential $V\left(x\right)$ is necessary, but in which the phenomenon of the avoided level crossing is expected to occur while being not so easily detected by the mere finite-precision experimental measurements of the spectra.

Last but not least let us mention that among the physical systems to which our present methodical innovation could also appear, in the nearest future, to be immediately applicable are also, potentially, all of the various more-dimensional realistic multi-well models as encountered, e.g., in quantum chemistry [40] and/or in the quantum-computation-motivated mesoscopic physics of the systems of quantum dots [41]. Indeed, under the widespread interaction-analyticity constraint, just the approximate solutions of the partial differential Schrödinger equations are usually available in these systems.

In the methodical context, a strengthening and/or independent complement of our present arguments may be also found in Ref. [11] where we sampled the methodical as well as practical gains of the approach by introducing a spiked, centrally symmetrized 1-D potential of the Morse-oscillator type. In spite of its non-analyticity in the origin the model was still shown to exhibit several features of the more conventional complete exact solvability.

For the similar and generic 1-D interactions characterized by the various incomplete forms of solvability (such as QES) the situation appears mathematically more complicated because the presence of the non-analyticities at some coordinates may restrict, severely, the practical as well as methodical applicability of the incomplete constructions in quantum physics. One can feel encouraged by the existence of parallels between the one-dimensional and three-dimensional central interactions. The phenomenological appeal of the latter models is usually perceived as independent of the presence or absence of a non-analyticity in the origin. Pars pro toto, let us just mention the singular nature of the so-called shape-invariant realizations of supersymmetry in quantum mechanics [2].

There are no reasons for an asymmetric treatment of solvability (and, in particular, of the QES models) in one and three dimensions. As long as in 3-D models the singularities in the origin are currently tolerated, a transfer of the same “freedom of non-analyticity” to 1-D systems is desirable. In our present paper we showed that such a change of paradigm is also productive (new models may be constructed) and satisfactory (in the QES setting, the classification of models may be formulated as starting from the most elementary quadratic polynomials).

After a removal of the QES constraints, the next-to-elementary, NES nature of our present bound-state problem (1) + (2) can be perceived, on a purely methodical level, as an inspiration of possible piecewise-analytic shifts, matchings and generalizations of many other traditional phenomenological interactions. Even if we restrict attention just to the present example (2) offering a tunable transition between a single- and double-minimum shape, we would like to emphasize that it possesses several specific advantages.

First of all, the exact analytic form of the related NES wave functions offers an undeniable amendment, e.g., in comparison with the conventional polynomial options as sampled by the famous Mexican-hat (MH) central potential

$$V^{\left(MH\right)}\left(\right|\overrightarrow{r}\left|\right)=-a|\overrightarrow{r}{|}^2+{\left|\overrightarrow{r}\right|}^4.$$

Thus, one of the domains of possible applicability of our model (2) might be sought in the advanced field theories of self-interacting systems possessing the so called Nambu-Goldstone bosons in which the use of benchmark $V^{\left(MH\right)}\left(\right|\overrightarrow{r}\left|\right)$ remains popular in spite of its purely numerical nature (cf., e.g., picture Nr. 7 in [42]). Not insisting too much on the analyticity of $V^{\left(MH\right)}\left(\right|\overrightarrow{r}\left|\right)$, the author of loc. cit. emphasized that with $a>0$, the MH dynamical ansatz leads to the solutions which “have lower symmetry than the Lagrangian, and contain mass zero bosons”.

We believe that in comparison with $V^{\left(MH\right)}\left(\right|\overrightarrow{r}\left|\right)$, our present model $V_{\left[d\right]}^{}\left(x\right)$ could prove competitive, in spite of its non-analyticity, not only in the above-mentioned area but also in various more straightforward applications ranging from the elementary quantum mechanics of atoms and molecules up to various topics in condensed matter physics. In all of these areas the possible phenomenological use of our NES quantum model (1) + (2) might complement the purely numerical approaches which are needed when the currently more popular polynomial potentials are used.