Miklos Rontoab and
Eli Pollak*a
aChemical and Biological Physics Department, Weizmann Institute of Science, 76100 Rehovot, Israel. E-mail: eli.pollak@weizmann.ac.il
bSchool of Chemistry, University of Leeds, Leeds, LS2 9JT, UK
First published on 18th September 2020
The accurate determination of tunneling splittings in chemistry and physics is an ongoing challenge. However, the widely used variational methods only provide upper bounds for the energy levels, and thus do not give bounds on the gap between them. Here, we show how the self-consistent lower bound theory developed previously can be applied to provide upper and lower bounds for tunneling splitting between symmetric and antisymmetric doublets in a symmetric double-well potential. The tight bounds are due to the very high accuracy of the lower bounds obtained for the energy levels, using the self-consistent lower bound theory. The accuracy of the lower bounds is comparable to that of the Ritz upper bounds. This is the first time that any theory gave upper and lower bounds to tunneling splittings.
The theoretical study of tunneling splittings in molecular systems has a long history. Widely used semiclassical approaches are based on the Wenzel–Kramers–Brillouin (WKB)4 and instanton methods.5 The WKB approach has been applied successfully to several model potentials, such as one- and multidimensional symmetric double-wells,16–18 one-dimensional asymmetric potentials,19 ring-puckering vibrations,15 and in combination with vibrational perturbation theory.20 A path-integral molecular dynamics method was used for double-well potentials in ref. 21. A perturbative instanton approach was used for non-rigid molecules in ref. 10. The tunneling splittings of malonaldehyde and water clusters were studied by a path-integral formulation22 and a ring-polymer implementation of instanton theory.23,24 These two systems were also studied by a wide range of approaches using accurate potential energy surfaces based on semiclassical,25 WKB,26 diffusion Monte Carlo (DMC),13,27,28 and multiconfigurational time-dependent Hartree (MCTDH) methods.29 Ab initio “on-the-fly” techniques combined with semiclassical trajectories were also employed for malonaldehyde30 and for the deep tunneling splittings in ammonia.31
Another class of theoretical methods used for the study of tunneling splittings are variational strategies, which are based on the diagonalization of the Hamiltonian matrix. For example, a calculation of tunneling splitting in phosphine was performed by exploiting symmetries in ref. 12. Malonaldehyde was studied using the Lanczos construct.14
It is then not surprising that bounding tunneling splittings is an old challenge. Simon32 worked out the leading asymptotics for an infinitely large one-dimensional potential. Kirsch and Simon33 derived lower bounds on eigenvalue splittings. Hemmen and Wreszinski34 derived an upper bound for the tunneling rate of a large quantum spin. Abramovich35 gave a lower bound comparison theorem. A workshop on the problem was held in 2006.36 More recently, Yu and Yang37 provided lower bounds for the gap between two adjacent eigenvalues. Semiclassical theories and estimates abound.38,39 However, even to date, a general theoretical framework for providing computational upper and lower bounds on tunneling splittings is lacking; this is the topic of this paper.
The main problem is that upper bounds to eigenvalues, provided by the Ritz–MacDonald variational procedure are not sufficient. It is necessary to obtain accurate lower bounds and it is the combination of the upper and lower bounds on the eigenvalues that gives upper and lower bounds on the gap between the eigenvalues. In a recent series of publications,40–43 we have developed a highly accurate Self-Consistent Lower Bound Theory (SCLBT), which can provide lower bounds to eigenvalues that are of similar and sometimes even have greater accuracy than those of the upper bounds obtained with the Ritz–MacDonald variational method.44,45 The theory generalizes and improves upon Temple's lower bound expression.46 The SCLBT is based on the application of the Lanczos algorithm47 for the tridiagonalization of matrices. Through the Lanczos transformation, explicit relations can be obtained between the overlaps of the exact and approximate wavefunctions, and these are then used extensively to derive the SCLBT expressions.
Thus far, the theory, described in detail in ref. 42 and 43, has been applied successfully to some nontrivial lattice models. The aim of this paper is to further explore the theory and apply it to the challenging problem of finding tight bounds to eigenvalue differences such as in tunneling splitting of eigenvalue energies. For many years, lower bounds were not sufficiently accurate.48–63 Especially when considering tunneling doublets, their accuracy was too poor, the gap between the lower bounds and the true energies was typically much larger than the actual energy differences between the levels. Upper and lower bounds for level differences are only meaningful if the accuracy of the upper and lower bounds is comparable and tighter than the energy splitting between the doublet states. In this paper, we demonstrate that the SCLBT is sufficiently accurate to determine tunneling splittings in a double-well potential. The lower bounds are as accurate as the Ritz–MacDonald upper bounds, allowing us to accurately bound the tunneling splitting energies.
The working expressions of the SCLBT are reviewed briefly in Section 2. The initial input required to apply the SCLBT is a not very accurate set of lower bounds for the eigenvalues of interest. Previously, we have used the Weinstein lower bound expression48 for this purpose. However, the conditions of validity of the Weinstein bound are not trivial (see Appendix A in ref. 43) and it is not always obvious how to determine without further information if the Weinstein lower bounds are truly valid. In addition, the SCLBT limits the highest state to be considered by the requirement that the corresponding Ritz eigenvalue (λj) needs to be less than the exact eigenvalue of the nearest state εj+1. The application of the method is predicated on an objective way for determining if and when these conditions hold. In previous implementations this was not a limitation; however, in the present case of a double-well potential we found that for states above the barrier the Weinstein lower bound could not be applied, as there was no way to determine when it was valid. To overcome this problem, we suggest here the application of a semiclassical estimate of the energy levels in question. The semiclassical action quantization condition (n + 1/2)h provides a reasonable estimate of the location of energy levels, especially those above the potential barrier separating between the wells. This information is then used to obtain the required initial lower bound to the energy levels as well as to determine the highest state for which the theory can be applied. This strategy is robust, since the SCLBT is not very sensitive to the accuracy of this initial input.
The application of the SCLBT to a quartic double-well potential supporting two tunneling doublets is described in Section 3. Using 11 Lanczos basis states, the lowest tunneling doublet was found to have a splitting 0.00199 ≤ ΔE ≤ 0.00224 (a.u.) while the second doublet was bounded as 0.123 ≤ ΔE ≤ 0.154. The numerically exact tunneling splittings are 0.00210 and 0.13562 for the two doublets, respectively.64 The tightness of the lower and upper bounds allows us to improve upon the estimate of the tunneling splitting value itself by taking the average of the upper and lower bounds (0.00212 and 0.13832 for the two doublets), proving highly accurate means. The paper ends with a summary of our conclusions.
Ĥ|φj〉 = εj|φj〉, j = 1, 2, …, | (2.1) |
(2.2) |
Hj,k = 〈Ψj|Ĥ|Ψk〉. | (2.3) |
(2.4) |
ĤL = LĤL | (2.5) |
ĤL|Φ(L)j〉 = λ(L)j|Φ(L)j〉, j = 1, …, L, | (2.6) |
(σ(L)j)2 = 〈Φ(L)j|Ĥ2 − (λ(L)j)2|Φ(L)j〉. | (2.7) |
The derivation of the SCLBT is described in detail in ref. 43; here, we give the final working expressions only. A key element of the theory is that the set of basis vectors |Ψj〉 is constructed using the Lanczos methodology. That is, for a chosen initial basis vector |Ψ1〉 the next vector is defined as
(2.8) |
(2.9) |
The set |Ψj〉, j = 1, …, L is by construction an orthonormal set and the resulting tridiagonal Hamiltonian is diagonalized to give the eigenvalues and eigenfunctions as defined in eqn (2.6).
As discussed in the introduction, a central restriction of the lower bound theory is that the highest state to be considered, denoted as L* has the property
(2.10) |
ε(L)j,W = λ(L)j − σ(L)j ≤ εj. | (2.11) |
We assume an L-dimensional projected space, and the highest eigenvalue for which the restriction of eqn (2.10) is assured is denoted by L*. Then the lower bound expression is
(2.12) |
(2.13) |
(2.14) |
(2.15) |
(2.16) |
(2.17) |
In these equations, notation ε−,j stands for a lower bound for the j-th state, λ(L)j,min represents a minimal upper bound for the j-th state, δkj is the Kronecker delta function, and θ(x) is the unit step function (defined as θ(x) = 0 if x ≤ 0 and unity otherwise). The implementation of the SCLBT is described in detail in ref. 42 and 43 and discussed again in detail in Section 3, where we present its application for obtaining upper and lower bounds for tunneling splittings.
If the maximal state is L* = 2, then fj,max = 0 (see eqn (2.15)) and the lower bound expression reduces to the simpler form of
(2.18) |
The first expression on the right-hand side is the generalization of Temple's lower bound for the ground state derived in ref. 40 and 41. The last expression is Temple's lower bound for the ground state.
(3.1) |
The lowest eigenvalues, reported in ref. 64, together with additional data on some higher lying states are listed in Table 1.
Symmetry | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
Sym | −6.64272703 | −2.45118605 | 0.41561275 | 3.87734184 | 8.76868217 | 14.4892323 |
Antisym | −6.64062823 | −2.3155705 | 1.67849653 | 6.19186637 | 11.53966844 | 17.60097423 |
Sym – scl | — | — | 0.10084667 | 3.83155489 | 8.75593471 | 14.47969607 |
Antisym – scl | — | — | 1.75893168 | 6.18035294 | 11.52914990 | 17.59218750 |
The initial wavefunctions chosen for the symmetric and antisymmetric states are normalized double Gaussians
(3.2) |
(3.3) |
(σ(L)j)2 = HL+1,L2〈ΨL|Φ(L)j〉2. | (3.4) |
Using these results, we compared the Weinstein lower bound with the Ritz upper bound as a function of the dimensionality. The results for the ground-state doublet are shown in Fig. 1. As expected, the upper bound is a decreasing function of the dimensionality of the projected space, while the Weinstein lower bounds are increasing. The Weinstein lower bounds are low quality as compared to the upper bounds; however, they are sufficiently accurate to use as inputs to the SCLBT.
The same is true for the first excited state doublet, shown in Fig. 2. Here, initially, the symmetric Weinstein lower bound decreases with the dimensionality, which is due to the fast decrease of the Ritz eigenvalues. In this region, we cannot know without additional information if the Weinstein lower bound is valid, but the fast decrease in the Ritz eigenvalue indicates that it might be too high. However, from L = 4 onwards, the Ritz eigenvalue plot flattens out and the Weinstein lower bound increases; thus, it can be assumed that the Weinstein lower bound is correct. Although for dimensions L = 2 and 3 there is no objective way to determine if the lower bound is valid, the additional information from L ≥ 5 shows that it is actually valid in these cases.
However, when considering the second excited state and higher energies difficulties arise. The Weinstein estimate for the lower bound of the second excited state (j = 3) is shown in Fig. 3. Two problems are highlighted by this figure. The first is that it is impossible to determine from the figure where the Weinstein lower bound is valid. In parts of the region it is greater than the correct second excited state energy indicated by the dotted black line. The semiclassical estimates for the third excited state energy are also shown. The second aspect is that, as highlighted in the previous section, the condition for the validity of the SCLBT is that the Ritz eigenvalue needs to be less than the adjacent higher energy. As can be seen in Fig. 3, this is only satisfied at L ≥ 10. This is the reason why we had to use L = 11 to obtain the accurate lower bounds reported below.
The next step is to increase the basis set to L = 3. Again, as a first step, the Weinstein lower bound for the first two energies obtained from the L = 3 calculation can be used, then the iteration can be performed using the improved lower bound for the ground state; this is continued until convergence is achieved. In this way, one obtains the SCLBT lower bounds and the Temple lower bounds as a function of the dimensionality of the Lanczos basis set used. These, together with the upper bounds are reported in Table 2 for the symmetric and antisymmetric states.
Type | L = 2 | L = 3 | L = 4 | L = 5 | L = 6 | L = 7 | L = 8 | L = 9 | L = 10 | L = 11 |
---|---|---|---|---|---|---|---|---|---|---|
ub, s | −6.5280 | −6.6209 | −6.6376 | −6.6413 | −6.6421 | −6.64238 | −6.64251 | −6.64259 | −6.64264 | −6.642663 |
ub, a | −5.9441 | −6.6292 | −6.6380 | −6.6396 | −6.64009 | −6.64031 | −6.64044 | −6.64052 | −6.640565 | −6.640588 |
lb, s, ST | −6.8336 | −6.7760 | −6.6876 | −6.6545 | −6.6464 | −6.64425 | −6.64355 | −6.64324 | −6.64306 | −6.642945 |
lb, s, T | −6.9110 | −6.8633 | −6.7193 | −6.6612 | −6.6487 | −6.64584 | −6.64490 | −6.64434 | −6.64387 | −6.643536 |
lb, a, ST | −6.9107 | −6.7232 | −6.6578 | −6.6455 | −6.6426 | −6.64170 | −6.64133 | −6.64110 | −6.64093 | −6.640822 |
lb, a, T | −7.0180 | −6.7552 | −6.6651 | −6.6483 | −6.64454 | −6.64350 | −6.64287 | −6.64217 | −6.64160 | −6.641257 |
The gaps between the upper and lower bounds and the exact eigenvalue are defined as λj − εj for the upper bound and εj − ε−,j for the lower bounds. The ratio of the gaps of the lower and upper bounds using the SCLBT is plotted in Fig. 4. The typical value of the gap ratio is between 3 and 5. This is not low, but still sufficiently low to provide non-trivial bounds on the tunneling splitting as shown in Fig. 5. For L = 11 the bounds on the tunneling splitting are 0.00184 ≤ ΔE ≤ 0.00236. In comparison, using the Temple lower bounds gives 0.00141 ≤ ΔE ≤ 0.00295. The numerically exact tunneling splitting is 0.00210.
The strategy we employed was to use a semiclassical estimate for the energies, based on the action integral equal to (n + 1/2)h. This might be an inaccurate estimate for the tunneling doublets, but for the doublets this is not a problem, as we have shown that the Weinstein lower bound is appropriate. For the states above the barrier separating the two wells the two turning points for E ≥ 0 can be found directly as and the action integral is
(3.5) |
This same strategy was used to identify when the SCLBT condition expressed in eqn (2.10) was satisfied. For j = 3 the symmetric eigenvalue with L = 10 is 3.606246, which is 0.23 less than the semiclassical eigenvalue. In the antisymmetric case, the eigenvalue at L = 10 is 4.205842 while the semiclassical energy is 6.180353; thus, eqn (2.10) is satisfied and from L = 10 we can use the SCLBT with L* = 3. For L = 11 the fourth symmetric eigenvalue is 8.266892, which is less than the semiclassical energy of 8.755935 while for the antisymmetric manifold the fourth eigenvalue is 7.19180, which is much less than the semiclassical energy of 11.529150.
As a consequence, for L = 10 and L = 11 we can now obtain lower bounds also for the second tunneling doublet. First, we consider the case of L* = 3. This can be used for both the L = 10 and L = 11 dimensional Lanczos space. These results are given in the first two rows of Table 3. There is a noticeable improvement of the lower bounds as compared with those derived at L* = 2 and as the number of Lanczos vectors are increased from L = 10 to L = 11. The related gap ratios are presented in Table 4, and it is obvious that increasing L* to 3 significantly increases the quality of the lower bounds. Instead of the gap ratios in the order of ∼3–5 at L* = 2, they are in the order of ∼1–2 in this case. Obviously, the upper and lower bounds on the tunneling splittings are also improved accordingly.
L | ε−,1,s | ε−,1,a | ε−,2,s | ε−,2,a | ΔE1 | ΔE2 |
---|---|---|---|---|---|---|
10 | −6.642877 | −6.640737 | −2.469336 | −2.323896 | 0.001898 ≤ ΔE ≤ 0.002312 | — |
11 | −6.642826 | −6.640699 | −2.464559 | −2.322290 | 0.001963 ≤ ΔE ≤ 0.002238 | 0.119644 ≤ ΔE ≤ 0.153684 |
11 | −6.642836 | −6.640674 | −2.470646 | −2.318977 | 0.001989 ≤ ΔE ≤ 0.002249 | 0.122957 ≤ ΔE ≤ 0.159771 |
L | L* | R1,s | R1,a | R2,s | R2,a |
---|---|---|---|---|---|
10 | 3 | 1.628740 | 1.704331 | 1.638289 | 1.316617 |
11 | 3 | 1.534887 | 1.753226 | 1.445370 | 1.430994 |
11 | 4 | 1.696747 | 1.124009 | 2.103232 | 0.725466 |
For L* = 4, the results (last rows in Tables 3 and 4) become ambiguous. While a significant improvement can be seen for the antisymmetric lower bound energies, the symmetric case does not improve. The reason for this is that λ(11)4,s = 8.266892 is very close to the lower bound for the fifth state (8.768682). In the antisymmetric case, λ(11)4,a = 7.191803, which is substantially lower than the lower bound of 11.539668. According to eqn (2.13), when the Lanczos eigenvalue is too close to the next highest energy, the lower bound is not very accurate. This is also the reason why the L* = 4 case cannot be used to obtain lower bounds to the j = 3 state in each symmetry manifold. The resulting lower bounds are much too low.
Finally, given the high quality of the lower bounds, one can use the average of the upper and lower bounds to obtain a significantly improved estimate of the energy. For example, (λ(11)1,s + ε11,3−,1,s)/2 = −6.642745 and this can be compared with the numerically exact value of 6.642727. The mean value is nearly an order of magnitude more accurate than the estimate obtained from either the upper or lower bound separately. The same applies for the other states.
However, the SCLBT is more complex than the upper bound theory, where only a basis set and the diagonalization of the resulting Hamiltonian matrix is required. The larger the basis set, the more accurate the eigenvalues are. For the lower bounds this is not sufficient. Some initial input is required. In this paper we showed how the Weinstein lower bounds and, if necessary, semiclassical energy estimates can be used. In addition, the SCLBT limits the energy levels for which lower bounds can be obtained according to the condition in eqn (2.10). We do note though that if the inequality does not hold, then the upper bound estimate is rather poor as well. In this sense, it can be expected that for all levels where the upper bound is relatively accurate, in the sense that its distance from the exact eigenvalue is much less than the gap between the adjacent states; the same applies for the lower bounds obtained with the SCLBT.
The system studied here is one-dimensional. A next step is to show that the SCLBT also works well for multidimensional tunneling systems, where the density of states is much higher and thus, calculations may be more challenging. The most restricting condition for application of the lower bound theory is expressed in eqn (2.10) which demands that the Ritz upper bound for the L-th state needs to be lower than the true energy of the (L + 1)-th state. Especially when the density of states is high, this condition can be challenging. We do note that typically this situation is not characteristic of the ground and lowest excited state tunneling doublets in molecules since the tunneling splittings are typically much smaller than the lowest harmonic frequency characterizing the molecule. However, further studies are needed to ascertain the extent of the challenge.
In addition, in the multidimensional case, the required basis set has to be much larger, the semiclassical estimate of levels in multidimensional systems is not always straightforward and other strategies might be needed for the initial lower bound estimates. However, we do believe that the present results are a significant step forward in the development of the lower bound theory and obtaining upper and lower bounds to energy differences.
This journal is © The Royal Society of Chemistry 2020 |