Exactly solvable Schr\"odinger equation with double-well potential for hydrogen bond

TL;DR

This paper presents an exactly solvable double-well Schr"odinger equation model for hydrogen bonds, using confluent Heun functions, and demonstrates its effectiveness by modeling hydrogen bonds in KHCO3 with good experimental agreement.

Contribution

The authors develop a new exactly solvable double-well potential model for hydrogen bonds using confluent Heun functions, facilitating analytical solutions and practical applications.

Findings

Analytical solutions for hydrogen bond potentials are obtained.

Model shows good agreement with experimental data for KHCO3.

The potential is suitable for both ordinary and low-barrier hydrogen bonds.

Abstract

We construct a double-well potential for which the Schr\"odinger equation can be exactly solved via reducing to the confluent Heun's one. Thus the wave function is expressed via the confluent Heun's function. The latter is tabulated in {\sl {Maple}} so that the obtained solution is easily treated. The potential is infinite at the boundaries of the final interval that makes it to be highly suitable for modeling hydrogen bonds (both ordinary and low-barrier ones). We exemplify theoretical results by detailed treating the hydrogen bond in $KHCO_3$ and show their good agreement with literature experimental data.

Equations

$$\frac{{\hslash }^2}{2m}\frac{d^2\psi (x)}{d{x}^2}+\left[E-U(x)\right]\psi (x)=0$$

$$ϵ=\frac{8m{L}^{2}E}{{\hslash }^2{\pi }^2}y=\frac{\pi x}{2L}$$

$${\psi }_{yy}^{\prime \prime }(y)+\left[ϵ-U(y)\right]\psi (y)=0$$

$$U(y)=h{\tan }^{2}y+c\frac{\sin y}{{\cos }^{2}y}-b{\sin }^{2}y+a\sin y$$

$$z(z-1){\phi }_{zz}^{\prime \prime }(z)+\left[\alpha z^2+(\beta +\gamma +2-\alpha )z-\left(\beta +1\right)\right]{\phi }_z^{\prime }(z)+$$

$$\left[(\mu +\nu )z-\mu \right]\phi (z)=0$$

$$\delta =\mu +\nu -\alpha \frac{\beta +\gamma +2}{2}$$

$$\eta =\frac{\alpha \left(\beta +1\right)}{2}-\mu -\frac{\beta +\gamma +\beta \gamma }{2}$$

$$\phi (z{)}^{general}=A\mathrm{HeunC}\left(\alpha ,\beta ,\gamma ,\delta ,\eta ;z\right)+B{z}^{-\beta }\mathrm{HeunC}\left(\alpha ,-\beta ,\gamma ,\delta ,\eta ;z\right)$$

$$z=\frac{1+s}{2}$$

$$\left(1-s^2\right){\phi }_{ss}^{\prime \prime }(s)-\left[\frac{\alpha }{2}{s}^2+(\beta +\gamma +2)s+\left(\gamma -\beta -\frac{\alpha }{2}\right)\right]{\phi }_s^{\prime }(s)-$$

$$\left(\frac{\nu -\mu }{2}+\frac{\nu +\mu }{2}s\right)\phi (s)=0$$

$$s=\sin y$$

$${\phi }_{yy}^{\prime \prime }(y)-\left[\frac{\alpha }{2}{\sin }^{2}y+(\beta +\gamma +2)\sin y+\left(\gamma -\beta -\frac{\alpha }{2}\right)\right]{\phi }_y^{\prime }(y)\frac{1}{\cos y}-$$

$$\left(\frac{\nu -\mu }{2}+\frac{\nu +\mu }{2}\sin y\right)\phi (y)=0$$

$$\psi (y)=\phi (y)\exp \left(\frac{\alpha }{4}\sin y\right){\left(\cos y\right)}^{\frac{\beta +\gamma +1}{2}}{\left[\tan \left(\frac{\pi }{4}+\frac{y}{2}\right)\right]}^{\frac{\beta -\gamma }{2}}$$

$$\gamma =-\left(\beta +1\right)$$

$$h={\left(\beta +\frac{1}{2}\right)}^2$$

$$b=\frac{{\alpha }^2}{16}$$

$$c=\beta +\frac{1}{2}$$

$$a=\frac{\nu +\mu }{2}-\frac{\alpha }{4}$$

$$ϵ=\frac{\mu -\nu }{2}-\frac{1}{4}{\left(2\beta +1+\frac{\alpha }{2}\right)}^2$$

$$\alpha =-4\sqrt{b}$$

$$\beta =-\frac{1}{2}+\sqrt{h}$$

$$U(y)=h{\tan }^{2}y+\frac{\sqrt{h}\sin y}{{\cos }^{2}y}-b{\sin }^{2}y+a\sin y$$

$$\psi(y)=\exp\left(-\sqrt{b}\sin y\right)\left[\tan\left(\frac{\pi}{4}+\frac{y}{2}\right)\right]^{\sqrt{h}}{\rm HeunC}\Biggl{(}-4\sqrt{b},-\frac{1}{2}+\sqrt{h},$$

$$-\left(\frac{1}{2}+\sqrt{h}\right),2a,\frac{3}{8}-a-b-\frac{h}{2}-\epsilon;\frac{1+\sin y}{2}\Biggr{)}$$

$${\rm HeunC}\Biggl{(}-4\sqrt{b},-\frac{1}{2}+\sqrt{h},-\left(\frac{1}{2}+\sqrt{h}\right),2a,\frac{3}{8}-a-b-\frac{h}{2}-\epsilon;1\Biggr{)}=0$$

$${\psi }_i(y{)}^{normalized}=\frac{1}{Z_i}{\psi }_i(y)$$

$$Z_0=56;Z_1=-3338;Z_2=-41;Z_3=48$$

Full text

††journal: Chemical Physics Letters

,

1 Introduction

The search for a convenient analytically tractable double-well potential for the Schrödinger equation (SE) is a problem with long history (see e.g., [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13] and refs. therein). SE with a double-well potential is used in many well-known applications in physics and chemistry beginning from the inversion of the ammonia molecule and ending by heterostructures, Bose-Einstein condensates and superconducting circuits (see refs. in [6], [7], [8], [9]). Besides it is a basic model for the description of ordinary hydrogen bonds (HB) and so-called low-barrier hydrogen bonds (LBHB) [14], [15], [16], [17]. There are heated debates for more than twenty years that the latter may play an important role in biology (in particular in enzyme catalysis) [18], [19]. The characteristic feature of LBHB is that a proton occupies a central position in the double well and the zero-point vibrational energy is close to the barrier top. Thus the commonly used semiclassical approximation (WKB) is principally inapplicable to the description of LBHB. The researches have to resort to numerical solutions of SE that greatly restricts the possibilities for scanning of the parameter space and makes the analysis to be cumbersome and inconvenient for perception. For this reason a relevant exact solution of SE would be highly useful for the analysis of LBHB nature.

The potentials from [1], [2], [3], [4], [5] do not allow obtaining the natural eigenfunctions for the Schrödinger’s operator. For this reason the authors expand the solution of SE over some suitable full sets of functions. The coefficients in the resulting series are determined from the appropriate recurrence relations. The latter makes the obtained solutions to be rather inconvenient for applications. The potentials from [6], [7], [8], [9] are beneficially distinguished in this regard. For them the natural eigenfunctions for the Schrödinger’s operator can be expressed via the confluent Heun’s function (CHF) [20], [21]. Thus the obtained solution of SE can be expanded in a series over natural eigenfunctions of Schrödinger’s operator for the given potential. As a result the coefficients in the series become independent of one another. A lot of potentials for SE are shown to be solvable with the help of CHF [22]. The CHF is a known and by now well described special function which is a solution of the confluent Heun’s equation [20], [21], [23], [24]. In particular CHF is tabulated in Maple. The latter makes its usage to be a routine procedure. This fact renders the obtained solution of SE very convenient for applications. Recently the exact solution of the Smoluchowski equation for reorientational motion in Maier-Saupe double well potential was obtained via CHF [25], [26]. Here we apply similar approach to SE with a practically important type of double-well potential suitable for the description of hydrogen bonds.

Most of the potentials investigated with the help of CHF ([6], [7], [8], [12], [13]) belong to the so-called quasi-exactly solvable ones for which only part of the spectrum can be computed exactly. The quantization is obtained from the requirement for the wave functions to contain the factors that are terminated polynomials. For the case of CHF it is reduced to a couple of the so-called polynomiality conditions [23], [24]. One of them defines the quantization of the energy spectrum and the other sets the values of the parameters for the potential at which it is possible. Thus only at some peculiar values of the barrier height the exact energy levels and the corresponding wave functions can be computed. The potential investigated in [9] is beneficially distinguished in the regard that it is exactly solvable. It enables one to determine the bound state energies for an arbitrary set of potential parameters. The boundary conditions are obtained from a very ingenious analysis of the wavefunction behavior at infinity along the line of [27].

Unfortunately the potentials studied in [6], [7], [8], [9] are finite at the boundaries of infinite interval that makes them inapplicable to the investigation of HB. In the latter case the potential for the proton is naturally restricted by the location of heavy atoms (e.g., $O-H\cdot\cdot\cdot\cdot\cdot O$ or $N-H\cdot\cdot\cdot\cdot\cdot O$) and can be both symmetric or asymmetric even in at first sight symmetric cases like $O-H\cdot\cdot\cdot\cdot\cdot O$ (see e.g. [28]). For LBHB the separation of heavy atoms is shorter and the bond is stronger (e.g., $N-H\bullet\bullet\bullet O$) [18]. The location of the heavy atoms makes the interval to be finite and the potential to be infinite at the boundaries. For this reason in the present paper we construct an exactly solvable potential for which the whole spectrum can be computed like for that from [9] but which is infinite at the boundaries of the finite interval. The quantization results from the boundary conditions rather than from the requirement of terminated polynomials. The potential is infinite at the boundaries and thus is highly suitable for modeling HB and LBHB. We exemplify our solution by detailed calculation of the energy levels and corresponding proton wave functions for HB in $KHCO_{3}$ and show good agreement of theoretical results with literature experimental data.

The paper is organized as follows. In Sec. 2 the problem under study is formulated. In Sec. 3 the solution of SE is presented and the eigenvalues are found. In Sec. 4 the results for $KHCO_{3}$ are discussed and the conclusions are summarized. In Appendix some methodical comments are added.

2 Schrödinger equation and the potential

We consider the one-dimensional SE

[TABLE]

where $U(x)$ is a double-well potential that is infinite at the boundaries of the finite interval $x=\pm L$. We introduce the dimensionless energy $\epsilon$ and the dimensionless distance $y$

[TABLE]

so that $-\pi/2\leq y\leq\pi/2$. The dimensionless SE takes the form

[TABLE]

For the potential we choose the model form

[TABLE]

Here $h$ and $c$ determine both the barrier width and the asymmetry, $b$ is the barrier height parameter and $a$ is the asymmetry parameter (along with $c$). An example of this potential (with particular choice $c=\sqrt{h}$) is depicted in Fig.1.

The variation of $a$ practically does not touch upon the variation of the barrier width while the variation of $c$ alters the latter considerably.

3 Reduction of Schrödinger equation to the confluent Heun’s one

CHF $\varphi(z)\equiv{\rm HeunC}\left(\alpha,\beta,\gamma,\delta,\eta;z\right)$ (in the Maple notation) is a solution of the confluent Heun’s equation

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

It should be mentioned that there is the second independent solution of (5) but it is superfluous for the problem under consideration and is discarded further. Mathematically it means that in the general solution of (5)

[TABLE]

we set the constant $B$ to be zero $B=0$. Physically it stems from the fact that the second solution is inconsistent with the requirement that for the infinite at the boundaries potential the wave function must be zero there.

We introduce a new variable $s$ by the relationship

[TABLE]

The equation takes the form

[TABLE]

[TABLE]

We introduce a new variable $y$ by the relationship

[TABLE]

The equation takes the form

[TABLE]

[TABLE]

Finally we introduce a new function $\psi(y)$ by the relationship

[TABLE]

We set

[TABLE]

and denote

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From (14) and (15) we express the parameters $\beta$ and $\alpha$ via those of the potential $h$ and $b$ but various choice of the signs is possible. The requirement of the agreement of the obtained theoretical results with literature experimental ones (see next Sec.) leads us to the correct choice

[TABLE]

[TABLE]

The other choice of the sign for $\beta$ is inconsistent with the requirement that for the infinite at the boundaries potential the wave function must be zero there. The other choice of the sign for $\alpha$ leads to an nonphysical solution that has no bound states. As a result we obtain SE (3) with the potential

[TABLE]

whose solution is the wave function

[TABLE]

[TABLE]

The potential (21) is infinite at the boundaries. Thus the wave function must be zero there. This requirement at $y=-\pi/2$ is automatically satisfied due to $\tan$ in (22). The requirement at $y=\pi/2$ can be satisfied only by a specific constraint imposed on ${\rm HeunC}$ and thus yields the boundary condition for determining the energy levels. As a result we obtain the equation for eigenvalues by setting $y=\pi/2$

[TABLE]

Its solutions form the spectrum of eigenvalues $\epsilon_{n}$ where n=0,1,2,...\ for the energy $\epsilon$.

It should be noted that in a few cases the solution of the confluent Heun’s equation can be reduced to more simple functions (see [29], [30], [31], [32] and refs. therein). However such cases impose severe constraints on the parameters of the potential. For instance the case 1. from [32] in our notation requires $a=\sqrt{b}$. The resulting potential with such asymmetry is very far from what we need for the description of hydrogen bonds. Other cases from [32] relate the parameters of the potential with the energy $\epsilon$ so that exact solutions can be obtained only for some specific energy levels. In our opinion a realistic asymmetric double-well potential can not be warped to satisfy above mentioned constraints.

4 Results and discussion

Fig.1 shows that the parameters of the potential (21) can be chosen to provide precise description of the energy levels structure for a specific hydrogen bond. In Fig.1 the energy levels for the hydrogen bond in $KHCO_{3}$ crystal (experimental data are taken from [28]) are presented. The energy levels for the proton stretching mode along an ordinary hydrogen bond usually form a pair of doublets within the wells. The $O-H\cdot\cdot\cdot\cdot\cdot O$ in $KHCO_{3}$ crystal presents a typical example [28]. The ground-state splitting $\Delta E_{01}=E_{1}-E_{0}=216\ {\rm cm^{-1}}$ was determined from large scattering cross-section of protons. The transition frequencies for the upper states $\Delta E_{02}=E_{2}-E_{0}=2475\ {\rm cm^{-1}}$ and $\Delta E_{03}=E_{3}-E_{0}=2820\ {\rm cm^{-1}}$ were determined with help of the infrared and Raman spectroscopies. The distance between the minima of the potential $\approx 0.6\ \AA$ is known from the crystal structure. These experimental values are obtained from our dimensionless ones if we take $h=35$, $b=180$, $a=-8$ and set $L=0.68\ \AA$ (with taking into account that for a proton $m=1\ {\rm amu}$). The latter means that the distance between the minima is $\approx 0.65\ \AA$ that is in good agreement with the above mentioned estimate. The normalized wave functions are obtained from (22) as follows

[TABLE]

where

[TABLE]

They are depicted in Fig.2. The functions $\psi_{0}(y)^{normalized}$ and $\psi_{1}(y)^{normalized}$ correspond to the split ground state. The ground state splitting is rather large ($216\ {\rm cm^{-1}}$ in dimensional units). It results from the asymmetry of the potential. For the symmetric potential it would be smaller and the tunneling frequency would be $18\ {\rm cm^{-1}}$ [28]. Taking into account that $k_{B}T\approx 200\ {\rm cm^{-1}}$ at $T=300^{\circ}K$ one can conclude the functions describe the states with the population of the same order of magnitude (up to the factor $e^{-1}$) at room temperature. The highly delocalized functions $\psi_{2}(y)^{normalized}$ and $\psi_{3}(y)^{normalized}$ describe excited states. It should be stressed that the latter are rather close to the barrier top (see Fig.1). Thus the functions $\psi_{2}(y)^{normalized}$ and $\psi_{3}(y)^{normalized}$ can not be described in semiclassical approximation (WKB). The exact analytic representation of these wave functions is the merit of the present approach.

The potential in Fig.1 is slightly asymmetric (the left well is deeper than the right one). The wave functions are largely (but not entirely) localized in the left and right wells respectively. It is in accordance with the qualitative fact that a small fraction of the wave functions is delocalized in the other well [28]. The asymmetry of the potential leads to tiny delocalization of the probability density for the corresponding wave functions and is known to play a key role in spectroscopy and proton-transfer dynamics [28]. We obtain excellent agreement of our wave functions with those calculated numerically in [28] for a suitable double-well potential. However it should be stressed that although the numerical analysis of SE with a double-well potential can yield identical results with those of the analytic treatment the latter has a substantial advantage. The parameter space of a double-well potential is extremely vast and to find their necessary combination to describe experimental data for the energy levels structure of a specific hydrogen bond with the help of the numerical solution of SE seems to be a daunting task. In our case the solution of the obtained equation (23) with the help of Maple software package facilitates the problem enormously. The latter strategy provides an efficient tool for scanning of the parameter space of the model. Another advantage of our analytic solution over the numerical one is the flexibility in the implementation of subsequent analysis (calculation of transition rates for proton-transfer dynamics, their corrections due to tunneling effect, etc.). We do not touch upon this matter in the present article and leave it for further work. Nevertheless it is obvious that the availability of the exact analytic expression for the wave function enables one much wider possibilities in treating the system under investigation than purely numerical solution of SE.

One can conclude that the suggested double-well potential for the Schrödinger equation is exactly solvable. The analytic expressions for the wave functions are obtained via the confluent Heun’s one. The equation for the eigenvalues (i.e., energy levels) can be easily solved numerically with the help of Maple software package. Besides the potential is infinite at the boundaries of the final interval that makes it to be highly suitable for modeling hydrogen bonds (both ordinary and low-barrier ones). The parameter space for the potential is rich enough to provide precise quantitative description of relevant experimental data on hydrogen bonds.

5 Appendix

Although Maple treats (23) very efficiently it proves to be too time consuming in evaluating integrals with CHF, e.g. those needed for the normalization of the wave functions. The author of the present article can not explain the inefficiency of $int(...,numeric)$ in this case and refer the reader to [23] where some drawbacks of ${\rm HeunC}$ realization in Maple are expertly discussed. In contrast with the help of the series representation of ${\rm HeunC}$ (with large but finite number of terms) one can obtain a very efficient tool for evaluating the corresponding integrals. The series expansion of ${\rm HeunC}$ [20], [21], [24] in our case takes the form

[TABLE]

[TABLE]

The coefficient $v_{k}$ is determined from the three-term recurrence relation

[TABLE]

with the initial conditions $v_{-1}=0$, $v_{0}=1$. Here

[TABLE]

[TABLE]

[TABLE]

We find that the visual identity of the series representation of ${\rm HeunC}$ with that given by Maple can be attained if we take some large but finite number $N$ for terms in the series instead of the infinity in (26). $N$ can be small for $\psi_{0}(y)$ but increases with the number $i$ in $\psi_{i}(y)$ and for $\psi_{4}(y)$ it should be $N>230$. Having made sure that both representations of ${\rm HeunC}$ are visually identical we (making use of the series representation of ${\rm HeunC}$) perform the numerical calculations of the required integrals by standard methods. If one avoids the enclosed cycles in calling the three-term recurrence relation then the calculation proceeds within seconds on any PC of even low power.

Acknowledgements. The author is grateful to Dr. Yu.F. Zuev for helpful discussions. The work was supported by the grant from RFBR N 15-44-02230.