# The Hill determinant method revisited

J. Killingbeck, G. Jolicard, A. Grosjean

► To cite this version:

J. Killingbeck, G. Jolicard, A. Grosjean. The Hill determinant method revisited. Journal of Mathematical Chemistry, 2019, 57 (1), pp.107-118. 10.1007/s10910-018-0939-9. hal-02122919

HAL Id: hal-02122919

https://hal.science/hal-02122919

Submitted on 31 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Hill determinant method revisited

J. P. Killingbeck<fnref n="1" />,<fnref n="2" /> · G. Jolicard<fnref n="2" /> · A. Grosjean<fnref n="2" />

Received: 22 March 2018 / Accepted: 21 June 2018 / Published online: 1 August 2018
© The Author(s) 2018

## Abstract

A Hill method is shown to give accurate energies and expectation values for several double and triple well sextic potential problems, thus contradicting previous criticisms of the method.

**Keywords** Hill determinant · Energy level · Expectation value

## 1 Introduction

Besides the works which have used the Hill determinant method to obtain numerical energy levels for the Schrodinger equation there also exist many critical works which have examined the validity and range of applicability of the method [1–10]. Several of these works have claimed to expose serious defects of the method. The main defects imputed to the method are:

(1) that it does not work for polynomial potentials which include an $x^6$ term

(2) that it cannot give accurate information about the wavefunction

(3) that it can converge to the wrong energy

In the present work we effectively disarm these criticisms by developing a Hill method which gives accurate energies and expectation values for several complicated problems which involve a sextic potential. Section 1 describes some modern modifications to the Hill determinant method and summarizes the main criticisms which have been made of the method. Section 2 sets out the mathematical algorithm which are used in the calculations of the present work and Sects. 3 and 4 present a range of numerical results, often for the problems for which the method has been claimed to fail. The problems treated in Sect. 3 are mainly double well ones. Section 4 treats some subtle phenomena which can arise for triple well problems and which clearly demon-

<footnote>J. P. Killingbeck
j.p.killingbeck@hull.ac.uk</footnote>

<footnote n="1">¹ Department of Engineering, Robert Blackburn Building, University of Hull, Hull HU6 7RX, UK</footnote>

<footnote n="2">² Observatoire de Besançon - Institut UTINAM CNRS UMR 6213, Université de Franche-Comté, 41 Bis Avenue de L'Observatoire, BP1615, 25010 Besançon Cedex, France</footnote>
strate the effectiveness of the Hill method of the present work. A combination of the Hill method with hypervirial perturbation theory (HVPT) is shown to be particularly effective in identifying particular states. Section 5 gives a brief conclusion.

## 2 The criticisms of the hill determinant method

We shall treat in turn the three points listed in Sect. 1.

### 2.1 The case of a sextic potential

The standard Hill determinant approach to the Schrodinger equation

$$-\alpha D^2\psi + \Sigma V(M)x^M\psi = E\psi \quad (1)$$

is to write the wavefunction as a product of an exponential term and a series

$$\psi(x) = \exp(-F(x)) \Sigma A(N)x^N \quad (2)$$

Inspection of the works which claim that the method will not work correctly for a sextic potential term [7–10] shows that they deliberately present a very specialized version of the method, with the function F(x) chosen specifically in the form ax² + bx⁴ and (more importantly) with the parameters chosen to give the wavefunction exactly the correct asymptotic form to fit to the sextic term in the potential. Inserting this choice into the Schrodinger equation (1) then produces a recurrence relation for the A coefficients which can be represented as a matrix eigenvalue equation for a tridiagonal matrix, the Hill determinant for this case. It is by studying the zeros of this determinant that some authors claim to find deficiencies in the Hill determinant approach. However (unsurprisingly) with a different form of F(x) the resulting modified Hill determinant is not restricted to tridiagonal form and gives satisfactory numerical results even for sextic potentials [11–13].

### 2.2 Lack of information about the wavefunction

In fact many Hill determinant calculations in the literature have not attempted to use an exact asymptotic exponential factor in the wavefunction and have not treated the problem of finding the A coefficients as a banded matrix eigenproblem. They have made the simple universal choice F(x) = βx²/2 in the exponential factor, varying β to optimize the results. They have also adopted a point of view more appropriate to the construction of a power series solution to the Schrodinger equation, thus avoiding an explicit study of the associated Hill determinant (which would in any case not have the tridiagonal form used by critical authors). In early calculations the series appearing in Eq. (2) was actually summed for each trial energy by some authors [14–16]. However, it was discovered that summing the series is not necessary. Having derived the recurrence relation satisfied by the A(N) it suffices to use the criterion A(N, E) = 0

<figure><img src="image_3.png" /></figure>
for a sufficiently high N in order to find an eigenvalue E, thus studying the individual coefficients rather than the full sum of the series. This result was initially observed empirically [2, 17, 18] but Killingbeck explained it theoretically [19, 20] and also demonstrated it computationally [21], producing a “Hill-series method” by incorporating the Hill determinant method in a simple power series method. The Dirichlet boundary condition $\psi(x)=0$ can then be imposed at any chosen finite distance $x = X_0$, with the Hill determinant result being obtained when $X_0$ tends to infinity. With this Hill-series method it becomes possible to locate accurately the nodes of the wavefunction for any given state [22] and also to calculate expectation values directly without quadrature, so that in this generalized form the Hill determinant approach CAN give useful information about the wavefunction.

## 2.3 Convergence to the “wrong” energy

By varying the parameter $\beta$ it has been found possible to find energy levels for two different potentials (differing only in the sign of one potential coefficient) within the same calculation [23–25]. An apparently reasonable choice of $\beta$ can thus give what appears to be a wrong eigenvalue (since it belongs to the other potential); we suspect that this might be the origin of the claim about convergence to the “wrong” energy. This feature was discovered in the 1980s but this early work has perhaps been overlooked in some recent studies of the number of low lying states of a sextic potential which have a wavefunction representable by a finite series times an exponential factor [26]. For a quartic potential such solutions had earlier been found in [27] to be possible for a set of specific $V_4$ values in the potential $x^2 + V_4x^4$. Leach [28] cited some early work on this problem of quasi-exact solubility for the sextic potential but a modern detailed treatment of the problem was given by Handy and coworkers using their highly effective new development of the Hill determinant method, the orthogonal polynomial projection quantization method (OPPQ) [29, 30]. In the present work we shall use only the particular case of infinite $X_0$ in our Hill-series method [21] and thus treat only the standard Hill determinant approach which has been subjected to criticism in the literature. Because we do not construct the Hill determinant associated with our recurrence relation we shall usually refer to the technique simply as the Hill method. Section 7 gives more details about this feature.

## 3 A summary of the recurrence relation approach to the theory

In the calculations of the present work we study the particular case of the bound states arising from a Schrodinger equation with an even parity polynomial potential

$$ -\alpha D^2 \psi + (\Sigma V(M)x^{2M}) \psi = E \psi \quad (3) $$

and take the wavefunction in the form

$$ \psi(x) = \exp(-\beta x^2/2) \Sigma A(N)x^{2N+P} \quad (4) $$
Inserting this postulated wavefunction into (3) gives the recurrence relation for the A(N);

$$ (2N + P + 1)(2N + P + 2)\alpha A(N + 1) = [(4N + 2P + 1)\alpha\beta - E]A(N) - \alpha\beta^2 A(N - 1) + \sum V(M)A(N - M) \quad (5) $$

The integer P is the parity indicator (0 for even parity, 1 for odd parity). To initiate the recurrence relation we set $A(0) = 1$ and use a trial E value, seeking to adjust E so as to render A(N) zero for some sufficiently high N value. It is clear from Eq. (5) that the A(N) will decrease rapidly in magnitude as N increases, so that A(N) for a high N might be beyond computer underflow. However, this can be avoided by introducing the two new coefficients AE(N) and AV(N), which are the derivatives of A(N) with respect to E and to V(1), respectively. If we write Eq. (5) for brevity in the condensed form

$$ C(N + 1)A(N + 1) = C(N, E)A(N) - \alpha \beta^2 A(N - 1) + \sum V(M)A(N - M) \quad (6) $$

then inspection of (5) shows that the AE(N) and AV(N) satisfy the equations

$$ C(N + 1)AE(N + 1) = C(N, E)AE(N) - \alpha \beta^2 AE(N - 1) + \sum V(M)AE(N - M) - A(N) \quad (7) $$

$$ C(N + 1)AV(N + 1) = C(N, E)AV(N) - \alpha \beta^2 AV(N - 1) + \sum V(M)AV(N - M) + A(N - 1) \quad (8) $$

The initial values of AE(0) and AV(0) are both zero, along with the value $A(0) = 1$. Since all three recurrence relations use the same coefficients they can be used simultaneously. At a high N value we correct the trial E value by the Newton's rule formula

$$ E \rightarrow E + DE \quad \text{with} \quad DE = -A(N)/AE(N) \quad (9) $$

and find the expectation value $X2$ of $x^2$ from the Eq. [30]

$$ X2 = -AV(N)/AE(N) \quad (10) $$

Since all the calculated quantities depend only on the ratios of coefficients, we can avoid any overflow or underflow problems by dividing all the current coefficients by $A(N)$ at each step, thus leaving unchanged the ratios which are needed in Eqs. (9) and (10). There is a further drastic simplification which can be made and which is relevant to the comments of several critical authors who compare the N value needed to obtain accurate eigenvalues by the Hill method with the dimension of the matrix needed by other approaches (we note again that this criticism is adopting the approach via the actual Hill determinant). The storage array dimension needed for our approach is typically 5, whatever the value of N. Suppose, for example, that we have only non-zero coefficients $V(1), V(2)$ and $V(3)$ in the potential (sufficient to give a sextic

<figure><img src="image_4.png" /></figure>
potential). Then it is clear from (5) that we calculate A(N + 1) from the previous terms A(N) down to A(N−3),so that at each step only 5 terms are involved. We introduce numerical labels which omit the N, so that A(N + 1) down to A(N−3) are now denoted by A(1) down to A(−3),with a similar change for the AE and AV coefficients. We can then use only 5 storage spaces for each of the A, AE and AV coefficients, simply moving along the numerical values at each step while using the scaling factor 1/A(1) simultaneously for all the coefficients. At a chosen high N value we then use the values of A(1),AE(1) and AV(1) in Eqs. (9) and (10). To ensure that the initial E value is near to an energy eigenvalue in our calculations we have used two different approaches. Within the Hill method itself a preliminary scan over a selected range of E values can be made, picking out the E values at which F(E) = A(1) changes sign for some high N value and so providing a starting list of approximate eigenvalues to be used in the full calculation. Alternatively, for the case of a system involving several potential wells a preliminary HVPT calculation can often be performed, similarly giving a list of approximate energy eigenvalues. We combine these two approaches in our specimen calculations.

**4 Numerical results for double well systems**

The program which implements the algorithm described in the previous section is
very simple and permits the calculation both of the energy levels of a polynomial
potential and of the expectation value of any even power of x which appears in the
potential. To find ⟨x<sup>4</sup⟩, for example, we use A(N−2) instead of A(N−1) as the last
term in Eq. (8). Except when a different α value is stated we used the value α = 1 in
all the reported calculations. To check the method we carried out a test calculation for
the potential x<sup>2</sup> + x<sup>4</sup> and found that the virial relation E = 2⟨x<sup>2</sup⟩ + 3⟨x<sup>4</sup>⟩ was obeyed
accurately by the three calculated quantities. We then tested the method on examples
from several previously published works on polynomial potentials, improving the
accuracy of most of the published results; Chhajlany and Malnev [27] gave a useful
set of exact reference eigenvalues for the sextic perturbed oscillator. In particular we
applied the method described here to the sextic potential examples used in some of the
papers cited earlier which criticized the Hill determinant approach and which proposed
“improved” versions of it. In all cases our method extended to double precision the
results, often of lower accuracy, presented as illustrations of those improved methods.
The results given for several sextic double well potentials using the efficient orthogonal
polynomial projection method in [29] can also be reproduced using our version of
the standard Hill method. An analysis of the results of the Rayleigh–Ritz matrix
calculations of [31] for a double well problem with α = 1/2 and the potential V(x) =
−17.5√2x<sup>2</sup> + x<sup>6</sup> revealed several interesting results not noted in [31], which used only
an even parity basis and calculated the ground state energy and the ⟨x<sup>2</sup⟩ and ⟨x<sup>6</sup⟩ values
(but not ⟨x<sup>4</sup⟩)). The ground state energy which was obtained by a matrix calculation
in [31] was E = −40.52625282343567, with ⟨x<sup>2</sup⟩> = 2.72643167389269 These results
agree with ours but we also calculated ⟨x<sup>4</sup⟩> for this state to be 7.836655839176 and
the energy of the lowest odd parity state to be −40.52625279631701. This reveals
the small energy level splitting. It also shows the localisation of the ground state
<table><caption>**Table 1** E, $\langle x^2 \rangle$ and $\langle x^4 \rangle$ for the potential $-45x^2+x^6$, using $N=2000$ and $\beta=10$</caption><thead><tr><th>E</th><th>&lang;x<sup>2</sup>&rang;</th><th>&lang;x<sup>4</sup>&rang;</th></tr></thead><tbody><tr><td>-103.028012221359</td><td>3.720987936532</td><td>14.4177167326914</td></tr><tr><td>-103.028012221276</td><td>3.720987936585</td><td>14.4177167329559</td></tr><tr><td>-77.9932968326146</td><td>3.400925261909</td><td>13.223438309890</td></tr><tr><td>-77.9932967741217</td><td>3.400925297326</td><td>13.223438483054</td></tr><tr><td>-54.9153731701885</td><td>3.052164697119</td><td>11.964549031550</td></tr><tr><td>-54.9153571664213</td><td>3.052173792676</td><td>11.964592148268</td></tr><tr><td>-34.0411628329423</td><td>2.660163153695</td><td>10.580868079421</td></tr><tr><td>-34.0390572685100</td><td>2.661253846336</td><td>10.585882698980</td></tr></tbody></table>

wavefunction in the potential wells, since $\langle x^2 \rangle^2$ is only five percent less than $\langle x^4 \rangle$. We note that the use of the virial theorem would have permitted the calculation of $\langle x^6 \rangle$ directly from E and $\langle x^2 \rangle$ in [31].

As a new and more challenging test of our Hill method we have applied it some deep double well potential of the form $V(x)=x^6-Vx^2$. Table 1 shows some specimen energies and expectation values for the case $V=45$. Since the method treats states of each parity separately it can give the very small level splittings between even and odd partner states and show how it increases with increasing state number. The results also show that the expectation value $\langle x^2 \rangle$ is slightly larger for the odd partner state, since the odd wavefunction has to be zero at the origin and so is pushed out a little more towards the well bottom. For the ground state the $\langle x^2 \rangle$ value and the small percentage difference between $\langle x^2 \rangle^2$ and $\langle x^4 \rangle$ shows that the wavefunction is concentrated near the bottom of the potential wells. For this double well problems with small splittings it was not difficult to find the approximate trial energy to use to locate each state, since HVPT with origin at the bottom of the right hand well gives both the state number and the average of the even and odd level energies for each pair, with an accuracy of between five and ten decimal digits. In [32] the authors asserted that the Hill determinant method does not converge for the cases $V=8$ and $V=4$. The Hill method of this work gives double precision energies for the first ten states for both these $V$ values with the simple fixed choice of $\beta=5$ and dimension $N=200$. For the triple well problems treated in the next section there can exist two states of the same parity with a very small energy difference, so that special procedures are needed to find the two energies separately.

## 5 Two triple well problems

We have pointed out that critics of the Hill determinant method have often specifically used a function of the form $F(x)=ax^2+bx^4$ in the exponential factor of Eq. (4). Although we do not use this factor in our postulated wavefunction we can use it for the different purpose of constructing an interesting test potential for our method. By
making the special choice $F(x) = x^2/2 + Bx^4/4$ we find that the wavefunction $\psi = \exp(-F(x))$ satisfies the Schrodinger equation

$$ -D^2\psi + [ (1-3B)x^2 + 2Bx^4 + B^2x^6 ] \psi = \psi \quad (11) $$

with eigenvalue 1, regardless of the value of B. This invariance of the eigenvalue has been pointed out previously [33, 34]. The Hill method of this work does indeed give the lowest eigenvalue 1 for several trial positive B values. Inspection of Eq. (11) shows that if we take the terms involving B as a perturbation of the $x^2$ potential, then the energy perturbation series in powers of B should sum to 1 plus zero. It is only possible to calculate a few terms of the perturbation series to sufficiently high precision but they do indeed show this zero energy shift. The first order terms in B can be seen to be zero by using the two oscillator groundstate expectation values $\langle x^2 \rangle = 1/2$ and $\langle x^4 \rangle = 3/4$, so that $2B \langle x^4 \rangle - 3B \langle x^2 \rangle = 0$. Supersymmetry theory [35] predicts that by changing the coefficient of $x^2$ in Eq. (11) from $1-3B$ to $1+3B$ (with B a positive number) we should obtain a partner potential for which the nth energy eigenvalue is 2 less than the $(n+1)$th eigenvalue of the potential in Eq. (11). The Hill method of this work confirms this prediction to double precision accuracy. It is clear that if B is negative then the energy perturbation series will still give the perturbed groundstate energy to be 1 for the potential in Eq. (11). However, for a negative B our postulated wavefunction $\exp[-F(x)]$ diverges at infinity and so cannot give a normalised wavefunction, even though the potential in Eq. (11) becomes a triple well potential which actually will have properly normalised bound state wavefunctions. The accurate Hill method of this work shows that there is a very small extra term which makes the energy slightly greater than 1 for the groundstate and which represents a non-analytic correction term to be added to the sum of the perturbation series. Table 2 shows how this term varies with the value of B. By a numerical analysis of the energy shift $\epsilon = E-1$ we find that its variation with B is well described by the non-analytic function $A \exp(-0.5/B)$, with a roughly optimized A value of 0.82. Thus halving B roughly squares the value of $E-1$, as can be seen from the table. The existence of the small non-analytic energy term was predicted theoretically in [36], which ascribed the existence of this term to a very small tunnelling effect between the outer and inner potential wells. Although the energy acquires an extra non-analytic term for negative B values our results showed that $\langle x^2 \rangle$ for the ground state varies smoothly across $B=0$ and has the perturbation series with the three leading terms $\langle x^2 \rangle = 0.5-0.75B+3B^2$. This is possible because, although the potential develops an outer pair of extra potential wells when B is negative, the groundstate wavefunction remains concentrated in the inner well, which varies smoothly as B varies through zero. The triple well potential which arises for negative B value displays yet more subtle properties when we look at the excited states. To give a particularly clear example of this behaviour we take the value $B=-0.01$. For this value the non-analytic groundstate energy term is so small that it is invisible to a double precision calculation of the energy but the triple well structure has an outer well so deep that it can contain several bound states. In such a situation we have a degenerate triple well potential

$$ V_1(x) = 1.03x^2 - 0.02x^4 + 0.0001x^6 \quad (12) $$
<table><caption>**Table 2** The ground state energy shift $E-1$ and $\langle x^2 \rangle$ for the potential $(1-3B)x^2+2Bx^4+B^2x^6$ with $B<0$, using $\beta=5$ and $N=1200$</caption><thead><tr><th>B</th><th>E-1</th><th>⟨x<sup>2</sup>⟩</th></tr></thead><tbody><tr><td>-0.02</td><td>1.231 (−11)</td><td>0.516379751224</td></tr><tr><td>-0.03</td><td>5.0781423 (−8)</td><td>0.525888883736</td></tr><tr><td>-0.04</td><td>3.2464302 (−6)</td><td>0.536672394571</td></tr><tr><td>-0.05</td><td>3.9178546 (−5)</td><td>0.548886614039</td></tr><tr><td>-0.06</td><td>2.0539406 (−4)</td><td>0.561941132068</td></tr><tr><td>-0.07</td><td>6.6875921 (−4)</td><td>0.574578581725</td></tr><tr><td>-0.08</td><td>1.6174884 (−3)</td><td>0.585516196259</td></tr><tr><td>-0.09</td><td>3.2105262 (−3)</td><td>0.593909260841</td></tr><tr><td>-0.10</td><td>5.5522813 (−3)</td><td>0.599435651280</td></tr></tbody></table>

which gives a severe test for the Hill method. Quoting to only 2 decimal places, the symmetric potential $V_1(x)$ has maxima of 15.84 at $|x|=5.90$ and outer well minima of 2.98 at $|x|=9.92$. To obtain a preliminary view of the bound states which exist for this triple well system we have done two HVPT calculations, one based on $x=0$ and the other based on the bottom of the right hand outer well. Table 3 includes the results of the two HVPT calculations and shows that for this degenerate triple well potential some of the energies are indeed the same for both inner and outer wells. A large number of digits in an HVPT energy indicates a highly localized wavefunction for the energy level concerned. We might thus anticipate some kind of resonance effect for these levels, leading to a splitting when the full space region is studied. A Hill method scan with a very small energy increment was then carried out and revealed the existence of close triplets of energy levels associated with the triple energies arising in the HVPT calculation. There are also some low energy single levels which arise naturally from a small perturbation of the central well oscillator potential $1.03x^2$. The results of the scan are shown in Table 3 along with the HVPT results which gave the information about the narrow regions to be scanned. To deal with very close energy values the Newton’s method equation of Eq. (9) is provided with a shift limiting procedure. A small maximum allowed shift DEM is chosen and then the normally calculated shift DE is divided by $1+|R|$, where $R=DE/DEM$. This preserves the sign of the calculated DE but limits its magnitude to be less than DEM. This is necessary because of a defect in Newton’s method; if by chance the current E estimate is near a maximum of $A(N, E)$ then the calculated DE will be very large and the calculation might jump to a distant root of the function. The scan is started just above and then just below the closely spaced pair in order to find first one level of the pair and then the other. The results of the scan are shown in Table 4 and reveal that some levels have associated wavefunctions which are strongly concentrated in the outer wells, since the intervening potential maximum inhibits tunnelling effects for low energy states. The subtle effects arising for this system are clearly shown by the Hill method of the present work.

A simple matrix diagonalisation model can be used to describe our results for this triple well system. Denoting the localised wavefunctions of a given energy $A$ in the left (L), centre(C) and right(R) wells by $|L\rangle$, $|C\rangle$ and $|R\rangle$ we can set up a $3 \times 3$ Hamiltonian

<figure><img src="image_5.png" /></figure>
<table><caption>**Table 3** Results of the HVPT calculations and the Hill method scan process for the triple well system with B = −0.01</caption><thead><tr><th>HVPT</th><th>HILL SCAN (even, odd)</th></tr></thead><tbody><tr><td>1.000000000000</td><td></td></tr><tr><td>2.969026973019</td><td></td></tr><tr><td>[4.8739185603]</td><td>4.8739185(29e, 61o, 92e)</td></tr><tr><td>6.70977209</td><td></td></tr><tr><td>[8.4703909]</td><td>8.470(379e, 391o, 403e)</td></tr><tr><td>10.14759</td><td></td></tr><tr><td>[11.7298]</td><td>11.7(284e, 298o, 314e)</td></tr></tbody></table>

β and N are as for Table 1. The energies common to the x = 0 and x = 9.9229 HVPT calculations are shown in a square bracket. The other results are for the x = 0 calculation

<table><caption>**Table 4** Hill method energy levels and ⟨x²⟩ values for the triple well system with B = −0.01</caption><thead><tr><th>E</th><th>⟨x²⟩</th></tr></thead><tbody><tr><td>1.000000000000</td><td>0.50782027158</td></tr><tr><td>2.969026973019</td><td>1.557438867</td></tr><tr><td>4.873918528791</td><td>49.464</td></tr><tr><td>4.873918560318</td><td>96.246610034</td></tr><tr><td>4.873918591846</td><td>49.464</td></tr><tr><td>6.709772091912</td><td>3.9015277685</td></tr><tr><td>8.470378758659</td><td>48.2534456</td></tr><tr><td>8.470390902062</td><td>91.2659381126</td></tr><tr><td>8.470403083054</td><td>48.253914</td></tr><tr><td>10.147591892047</td><td>6.7454954006</td></tr><tr><td>11.728311284376</td><td>46.70524882</td></tr><tr><td>11.729823470610</td><td>84.967393928</td></tr><tr><td>11.731342268074</td><td>46.75438894</td></tr></tbody></table>

matrix with diagonal matrix elements all equal to A and with small off diagonal elements ⟨ClHIL⟩ and ⟨ClHIR⟩ both equal to a small number B representing tunnelling effects between inner and outer wells. ⟨LlHlR⟩ and all the overlap matrix elements can be taken to be effectively zero. It is then easy to find three energy eigenfunctions; (1,0,−1)ᵀ has the energy A, (1, √2, 1)ᵀ has the energy A+B√2 and (1, −√2, 1)ᵀ has the energy A−B√2. The Hill method results show this evenly spaced structure and could even yield a value for the tunnelling parameter B in the simple matrix model. We note that [37] gives a useful set of numerical results for several triple well potentials which have a sextic potential term. We checked that our Hill method can not only give the small energy level splittings for the potentials in [37] but can also give extra expectation values results which indicate the localisation of the wavefunctions.
## 6 Conclusion

The efficient version of the Hill determinant method presented here permits calculation of both energies and expectation values and can deal with the sextic potential terms which various critics have said that it cannot handle. It uses a simple exponential factor $exp(-βx^2/2)$ in the wavefunction, rather than the more complicated form which the critics have imputed to it in their comments. Critics have also demanded that a very small dimension of Hill determinant (typically 20) should be required to prove that the method gives converged results. We have pointed out in this work that there is no difficulty in using a very large apparent dimension in order to ensure double precision convergence, since the actual array dimension used is typically 5. We note that even with $N=2000$ each calculation takes about 3 s on a laptop computer. The numerical results presented here for subtle triple well problems should have served to disarm much of the criticism of the Hill determinant method which has appeared in the literature and also to make more widely known a useful extension of the method. The approach which we adopt here shows that by imposing the simple condition $A(N, E)=0$ on the coefficients of the power series we can obtain accurate information about energies and expectation values (and thus implicitly about the wavefunction). We have chosen the value of $β$ empirically but some authors have attempted to find approximate formulae to give a suitable $β$ value for given potential parameters [38, 39]. Our main purpose has been to show that, despite previous criticisms, the Hill method works well for complicated bound state calculations involving sextic potential terms. In the present work we have tested the Hill method by exploring several interesting phenomena associated with double and triple well bound states in one dimension. However, several extensions of the Hill method described here are possible. By setting $P=L+1$ and subtracting a term $L(L+1)α$ from the coefficient of $A(N+1)$ in Eq. (5) we can deal with states of angular momentum $L$. By using complex $A$ and $AE$ coefficients and a complex $β$ in the three recurrence relations of the present Hill method we have elsewhere shown that it is possible to calculate the complex energies of resonant states for a potential such as $x^2-λx^4$ [40].

## 7 Note added in revision (may 2018)

The referees have asked the authors to explain further how the present method overcomes the problems raised by the critical works [1–10]. As we have already indicated, the opponents of the Hill determinant method set up what in rhetoric would be called a “straw man” target FIRST by giving the wavefunction an exponential factor containing $x^2$ and $x^4$ terms which are specifically chosen to give the exact asymptotic form to fit to the $x^6$ term in a sextic potential and SECOND by then explicitly evaluating the associated Hill determinant numerically and searching for its zeros to find energy levels. The essential central claim made by the critical works is that the calculation does not produce a correctly square integrable wavefunction associated with each energy level. We thank one of the referees for stressing that this defect of the method has ONLY been established for the very special “exact” exponential factor case described above. Even within the Hill determinant calculation the problem can be cured by using

<figure><img src="image_6.png" /></figure>
a different exponential factor. Thus when some authors refer to “the well-known con-
vergence difficulties of THE Hill determinant method” they are falsely extending to
ALL possible exponential factors a phenomenon which has ONLY been shown to be
relevant for ONE particular choice of parameters.

However, the method of the present work uses a quite different approach which DOES NOT EVEN USE any of the procedures criticized in [1–10] To calculate expectation values directly the method does NOT need to use a quadrature or any procedure requiring the explicit presence of a normalised wavefunction; instead it uses an implicit differentiation process which essentially defines an expectation value as a perturbation response coefficient. Because the present author had already previously generalised the method to such an extent that it can be used to calculate energy levels and expectation values for the case of Dirichlet boundary conditions imposed at ANY finite x value, it is clear that within that generalised method the boundary value of infinite x is the appropriate one to use in order to simulate the calculations made in the various works which have previously dealt with the simple Hill determinant method, which DOES NOT have that flexibility in its boundary conditions. That infinite boundary condition method is the one described and applied throughout this work. Since the method given here does NOT evaluate the Hill determinant it would be inappropriate to call it a Hill determinant method and so it has been called the Hill method, because it DOES use the typical Hill wavefunction form involving an exponential factor multiplied by a series. In choosing the parameter $\beta$ in the simple exponential factor $\exp(-\beta x^2/2)$ in the wavefunction it is, of course, possible to undertake a sequence of calculations in which $\beta$ is varied until the dimension $N$ needed is made as small as possible. Clearly this optimised $\beta$ would vary as the potential varies, thus requiring a sequence of optimisations. In actual calculations it is found that even a simple value such as $\beta=5$ or 10 works for many different problems if $N$ is made sufficiently large; using $N=2000$ is very much more than sufficient to give speedy and accurate results in normal double precision (with an average calculation time of 3 s). In some trial calculations we found that the rules suggested in [38, 39] can often reduce the required $N$ value to several hundred. However, to avoid an extra preliminary calculation we find that use of a sufficiently high $N$ is a much more simple approach, since it renders the choice of $\beta$ much less critical. As explained in Sect. 2, the storage arrays used in the calculations of this work are actually ONLY of dimension 5, with the 5 elements being overwritten as the calculation proceeds up to the chosen $N$.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

References

1. S.N. Biswas, K. Datta, R.P. Saxena, P.K. Srivastava, V.S. Varma, Phys. Rev. D **4**, 3617 (1971)

2. A. Hautot, M. Nicolas, J. Phys. A Math. Gen. **16**, 2953 (1983)

3. A. Hautot, Phys. Rev. D **33**, 437 (1986)

4. R.N. Chaudhuri, M. Mondal, Phys. Rev. A **40**, 6080 (1989)

5. F.M. Fernandez, Phys. Rev. A **44**, 3336 (1991)
6. M. Znojil, Phys. Lett. A 155, 83 (1991)

7. R.N. Chaudhuri, Phys. Rev. D 31, 2687 (1985)

8. M. Tater, J. Phys. A Math. Gen. 20, 2483 (1987)

9. M. Tater, A.V. Turbiner, J. Phys. A Math. Gen. 26, 697 (1983)

10. M. Znojil, Phys. Rev. D 34, 1224 (1986)

11. S.N. Biswas, K. Datta, R.P. Saxena, P.K. Srivastava, V.S. Varma, J. Math. Phys. 14, 1190 (1973)

12. K. Banerjee, Proc. R. Soc. A 364, 265 (1978)

13. M. Bansal, S. Srivastava, Vishwamittar, Phys. Rev. A 44, 8012 (1991)

14. D. Secrest, K. Cashion, J.O. Hirschfelder, J. Chem. Phys. 37, 830 (1962)

15. J. Killingbeck, Phys. Lett. A 84, 95 (1981)

16. A. Palma, G. Campoy, Phys. Lett. A 121, 221 (1987)

17. C.A. Ginsburg, Phys. Rev. Lett. 48, 839 (1982)

18. C.J. Tymczak, G.S. Japaridze, C.R. Handy, X.Q. Wang, Phys. Rev. Lett. 80, 3673 (1998)

19. J. Killingbeck, J. Phys. A Math. Gen. 18, L1025 (1985)

20. J. Killingbeck, J. Phys. A Math. Gen. 20, 1285 (1987)

21. J. Killingbeck, J. Phys. A Math. Gen. 21, 111 (1988)

22. J.P. Killingbeck, J. Phys. A Math. Theor. 43, 085302 (2009)

23. J. Killingbeck, J. Phys. A Math. Gen. 19, 2903 (1986)

24. R.N. Chaudhuri, J. Phys. A Math. Gen. 21, 567 (1988)

25. M. Znojil, Phys. Lett. A 116, 207 (1986)

26. M.S. Child, S.H. Dong, X.G. Wang, J. Phys. A Math. Gen. 33, 5653 (2000)

27. S.C. Chhajlany, V.N. Malnev, J. Phys. A Math. Gen. 23, 3711 (1990)

28. P.G.L. Leach, Physica 17D, 331 (1985)

29. C.R. Handy, D. Vrinceanu, J. Phys. A Math. Theor. 46, 135202 (2013)

30. C.R. Handy, D. Vrinceanu, R. Gupta, J. Phys. A Math. Theor. 47, 295203 (2014)

31. P. Koscik, A. Okopinska, J. Phys. A Math. Theor. 40, 10851 (2007)

32. C.R. Handy, D. Vrinceanu, C. Marth, H.A. Brooke, Mathematics 3, 104 (2015)

33. I.W. Herbst, B. Simon, Phys. Lett. 78B, 304 (1978)

34. P. Amore, F.M. Fernandez, Ann. Phys. 385, 1 (2017)

35. J. Killingbeck, Phys. Lett. A 115, 253 (1986)

36. R.J. Damburg, R.K. Propin, J. Phys. A Math. Gen. 15, 3481 (1982)

37. N. Aquino, J. Garza, G. Campoy, A. Vela, Rev. Mex. de. Fis. 57, 46 (2011)

38. A. Hautot, A. Magnus, J. Comp. Appl. Math. 5, 3 (1979)

39. A.N. Drozdov, J. Phys. A Math. Gen. 28, 445 (1995)

40. J.P. Killingbeck, J. Phys. A Math. Theor. 40, 9017 (2007)