Comparison of an improved self-consistent lower bound theory with Lehmann’s method for low-lying eigenvalues

Ritz eigenvalues only provide upper bounds for the energy levels, while obtaining lower bounds requires at least the calculation of the variances associated with these eigenvalues. The well-known Weinstein and Temple lower bounds based on the eigenvalues and variances converge very slowly and their quality is considerably worse than that of the Ritz upper bounds. Lehmann presented a method that in principle optimizes Temple’s lower bounds with significantly improved results. We have recently formulated a Self-Consistent Lower Bound Theory (SCLBT), which improves upon Temple’s results. In this paper, we further improve the SCLBT and compare its quality with Lehmann’s theory. The Lánczos algorithm for constructing the Hamiltonian matrix simplifies Lehmann’s theory and is essential for the SCLBT method. Using two lattice Hamiltonians, we compared the improved SCLBT (iSCLBT) with its previous implementation as well as with Lehmann’s lower bound theory. The novel iSCLBT exhibits a significant improvement over the previous version. Both Lehmann’s theory and the SCLBT variants provide significantly better lower bounds than those obtained from Weinstein’s and Temple’s methods. Compared to each other, the Lehmann and iSCLBT theories exhibit similar performance in terms of the quality and convergence of the lower bounds. By increasing the number of states included in the calculations, the lower bounds are tighter and their quality becomes comparable with that of the Ritz upper bounds. Both methods are suitable for providing lower bounds for low-lying excited states as well. Compared to Lehmann’s theory, one of the advantages of the iSCLBT method is that it does not necessarily require the Weinstein lower bound for its initial input, but Ritz eigenvalue estimates can also be used. Especially owing to this property the iSCLBT method sometimes exhibits improved convergence compared to that of Lehmann’s lower bounds


Theory
General considerations. In this section, lower bound theories are considered for the energy levels of quantum mechanical systems and for simplicity, real valued functions are assumed. Let us consider a Hermitian Hamiltonian operator Ĥ with energy eigenstates |ϕ j � . The corresponding true energy levels ε j can be determined from the time-independent Schrödinger equation where the energy eigenvalues are in ascending order, with the ground state denoted by j = 1 (instead of j = 0 , typically accepted in the literature). The exact representation of the Hamiltonian operator using an orthonormal basis set {|� j � = 1, 2, . . .} can be written as with However, in numerical calculations the full Hamiltonian matrix H cannot be used: it needs to be represented in a finite basis set, with L states. Let us define a projection operator onto this L-dimensional subspace as and the projection onto the complementary space as Q L such that the combination is the identity operator in the full Hilbert space as The discretized Hamiltonian projected onto the finite basis can be written as and we assume that it can be diagonalized as where | L,j � are normalized eigenfunctions and L,j are real eigenvalues. According to the Ritz-Macdonald variational principle, any L,j eigenvalue gives an upper bound to the exact eigenvalue ε j as For each L,j eigenvalue, the corresponding variance σ 2 L,j is defined as (1) H|ϕ j � = ε j |ϕ j �, j = 1, 2, . . . , H jk = � j |Ĥ| k �.  showing that when using a Lánczos basis set, the variances can be obtained from the matrix elements of the Hamiltonian, without the explicit need to calculate elements of the Hamiltonian squared. Naturally, these are implicit in the Lánczos basis vectors.

Weinstein and Temple lower bounds.
To discuss lower bound theories on an equal formal footing, let us introduce the square of the overlap of the jth eigenfunction in the projected space with the exact kth eigenfunction as Lower bounds can be obtained by using a Cauchy-Schwartz inequality in the form of where Q is a projection operator. Inserting into the Cauchy-Schwartz inequality and rearranging gives and with the assumption that a L,jj ≥ 1/2 , the Weinstein lower bound can be obtained as As discussed in Ref. 51, this assumption is further restricted by the accepted condition for the validity of the Weinstein lower bound 5 : To obtain the Temple lower bound formula, let us first define a "residual energy" ¯ L,j for each eigenstate in the projected space such that the Ritz eigenvalue can be rewritten as Using the overlap defined in Eq. (13), the residual energy is given by the relation where ε k is the exact eigenvalue and δ jk is a Kronecker delta. Thus, from Eq. (19) the diagonal elements of the overlap can be written as The Temple lower bound expression can be directly obtained by inserting this identity into Eq. (16), giving (9) σ 2 L,j = �� L,j | Ĥ 2 − 2 L,jÎ |� L,j �.
where the previously unknown overlap a L,jj is replaced by the residual energy ¯ L,j . Although the residual energy is unknown as well, its lower bound estimates can be obtained directly. For example, from the definition in Eq. (20) for j = 1 it can be readily found that ¯ L,j ≥ ε 2 for the ground state. A natural choice is then to substitute the Weinstein lower bound for the corresponding state; thus, a practical method is obtained for the calculation of lower bounds by Temple's formula.
SCLBT. By substituting Q =Q L -the projection operator for the complementary space as defined in Eq. (5)-into the Cauchy-Schwartz inequality in Eq. (14) and using the residual energy, an improved lower bound inequality can be obtained: This bound is better than the Temple one, because the denominator in the second term is always greater than one. Although the overlaps a L,jk are unknown, they can be determined by considering Eq. (12) obtained from the Lánczos construct. Using the identity the overlaps can be rewritten in terms of eigenvalues, true energies, and variances as Inserting this expression into Eq. (23) gives a lower bound as 50,51 where T L,j is defined as with ε L,j as the lower bound to the jth eigenvalue based on the L-dimensional space. This lower bound expression is a significant improvement over Temple's lower bound. The lower bound provided by these expressions can be further improved by defining tighter bounds to the residual energy ¯ L,j and to the overlap a L,kk as discussed below.
iSCLBT. In order to improve the SCLBT, a better estimate to the residual energy ¯ L,j needs to be found. Let us first define an energy η L,L * with the condition where L * is the highest energy level, for which Eq. (28) can be satisfied also for the associated Ritz eigenvalue L,L * . The value of L * is a parameter that can be varied according to the quality and validity of the Weinstein lower bound. Using these notations, the residual energy for the jth state can be rewritten as which, considering that the second term with the infinite sum in Eq. (29) is positive, can be rewritten as For the working equations of the iSCLBT, estimates for η L,L * and a L,kk need to be obtained. For the Lehmann expression, discussed in detail in the next section, it can be assumed that L * is the highest value for which the simple Weinstein lower bound condition given in Eq. (18) is satisfied, which implies that However, for the iSCLBT η L,L * = L,L * can be used instead of using the Weinstein lower bound for ε L * . Typically, there is a range l min ≤ L ≤ l max of L values for which L * is the highest state for which the Weinstein lower bound is valid 6 . According to the Ritz theorem, L,L * is a decreasing function of L. However, the highest possible value, min L,L * , is required, which still satisfies Eq. (18): this is the worst Ritz estimate for ε L * , that is still lower than ε L * +1 . Thus, the validity condition for the iSCLBT for L * is weaker than the Weinstein lower bound condition, which enables to obtain lower bounds to excited states for lower L values compared to Lehmann's theory and provides improved lower bound estimates.
For the SCLBT 50,51 , the upper bound for a L,kk is given as which, using the Lánczos construct, can be rewritten as However, a tighter bound can be obtained by considering the Cauchy-Schwartz inequality (Eq. (14)) using the projection operator Q L without introducing the residual energy such that where T L,k is given in Eq. (27). Rearranging this expression gives a new bound for a L,kk as The difference between this expression and Eq. (35) is the appearance of unity in the denominator, which provides a tighter upper bound to a L,kk . When the Weinstein lower bound is valid, L,k−1 ≤ ε k ≤ L,k ; thus, the maximal a L,kk as a function of ε k monotonically increases from 0 at ε k = L,k−1 to 1 at ε k = L,k . Using this expression and considering the above mentioned conditions for η L,L * , Eq. (32) can be maximized as where the exact eigenvalues have been replaced by their lower bounds ε L,k . The lower bound equation (31) can be rewritten as The Lehmann method requires the matrices of Ĥ 2 and Ĥ in the L-space, and the lower bound eigenvalues can be obtained by the diagonalization of Eq. (41). According to the condition in Eq. (42), a non-trivial lower bound needs to be estimated for the state ε L+1 , which can typically be a Weinstein lower bound. Nevertheless, when a Lánczos basis is used, the full Ĥ 2 matrix in the projected space is not required, but only the variances σ 2 L,j associated with the respective Ritz eigenvalues. This can be shown by multiplying Eq. (41) by � L,k | , which gives such that The variances can be obtained from Eq. (12), that is, all the information is in the matrix elements of the Lánczos representation of the Hamiltonian only.

Results
For the comparison of the lower bounds calculated using the iSCLBT method with those of the SCLBT and Lehmann methods, the same two lattice systems were used as in Ref. 50 where i denotes the lattice sites, c † i,σ and c i,σ are creation and annihilation operators for the electron in site i with spin σ =↑, ↓ , respectively, n i,σ = c † i,σ c i,σ is the number operator for the spin state σ on site i, ǫ is the "on-site" energy, which is set to zero in this case, t is the hopping energy between the nearest neighbors, and U > 0 is the Coulomb repulsion experienced by two electrons occupying the same site. The other model was the Heisenberg system 54 , describing a set of spin-1/2 particles on a lattice, given by the Hamiltonian where s i and s j are the spin-1/2 operators on the lattice sites i and j, respectively, i, j denotes the nearest-neighbor pairs only, and J is the coupling (or "exchange") constant weighting the "exchange term" s i s j . The detailed analysis of these systems is outside of the scope of this study: they were chosen because they are well suited for diagonalization using the Lánczos algorithm. The data used in this study was obtained using the H diagonalization software 55 and available in Ref. 50. In these calculations, the lattices were always considered periodic and the simulation cell was limited to a finite number N of sites with periodic boundary conditions. The diagonal ( α ) and off-diagonal ( β ) Lánczos coefficients of a Heisenberg square lattice with a unit cell of 5 × 6 and of a Hubbard square lattice model at half-filling on a unit cell of 4 × 4 were used. A tridiagonal matrix was created from the α and β coefficients and then, diagonalized using the double-precision Lapack dsyev subroutine 56 for real, symmetric matrices. The L,j energy eigenvalues become unstable, that is, they started to increase and fluctuate, at L > 118 for the Hubbard and at L > 84 for the Heisenberg model. The Weinstein lower bounds were calculated using Eq. (17), and their ranges of validity 5 were calculated according to Eq. (18). As the exact eigenvalues ε j were not known, the lowest stable eigenvalues M,j were used, at M = 117 for the Hubbard and M = 83 for the Heisenberg system. However, these values are sufficient, as only an estimation of the lowest L is required, from where the Weinstein lower bound can be considered valid. As L * = 2 is the lowest reference level for both the iSCLBT and Lehmann calculations, it is sufficient to determine the ranges of validity from L,2 . The graphs of the eigenvalues as a function of dimensionality L were intersected by the line from the validity condition in Eq. (18), than the next highest integer L provided the lowest limit of validity. Figures 1 and S1 show the ranges of validity of the Weinstein lower bounds for the Hubbard and Heisenberg models, respectively. These ranges were used to determine the range of validity of the SCLBT and Lehmann calculations. As can be seen in Figures 1 and S1 at L,5 there are only a few remaining valid states; thus, the highest achievable level is L * = 5 for both the Hubbard and the Heisenberg systems.
Given the Weinstein lower bounds and so the corresponding residual energies, we first compared the performance of the iSCLBT described in the previous section with that of the previous SCLBT implementation given in Refs. 50,51. The quality of the lower bounds from these methods was compared using their gap ratios, defined as the ratio of the distances of the lower bound ε L,j and the eigenvalue L,j to the true energy level ε j as where M,j is the lowest stable eigenvalue.
For the ground-state residual energy, ¯ L,j , the first-excited state Weinstein lower bound ε We 2 was used for the SCLBT method based on Eq. (26). In the iSCLBT method, the ε L,j bounds are calculated iteratively, by substituting the previously calculated lower bound back to the expression, until convergence is achieved. The residual energy ¯ L,j was estimated from the Weinstein lower bound ε We L,j . The SCLBT 51 improves Eq. (26) by self-consistently considering several states up to L * ; as shown in Ref. 52 increasing L * results in tighter lower bounds. For the iSCLBT method, Eqs. (27) and (38)-(40) were solved simultaneously and iteratively. The min L,L * value in Eqs. (38) and (40) was determined based on the validity of the Weinstein lower bound: this was the highest eigenvalue where the Weinstein lower bound was still valid. In the first step of the iteration, the ε L,j values in Eq. (27) were estimated from the Weinstein lower bound ε We L,j and the subsequent iteration steps use the previously calculated ε L,j values. Although any other estimated value could be used, such as the result ε L−1,j from the previous state, www.nature.com/scientificreports/ the calculation converged to the same value, regardless of the initial input; thus, the method does not explicitly rely on the Weinstein lower bound. The lowest L ini value where proper convergence was achieved was slightly higher than the L value from where the Weinstein lower bound was valid. A comparison of the SCLBT and iSCLBT results at L * = 2 , L * = 3 , L * = 4 , and L * = 5 is shown in Fig. 2. The noticeable variation of the functions is due to the fluctuation in the variances obtained from the Lánczos construct. As mentioned earlier, the SCLBT calculations start converging slightly later than the point where the Weinstein lower bound becomes valid. For example, the Weinstein lower bound for the third state ε We 3 becomes valid at L ≥ 62 , as shown in Fig. 1, while the iSCLBT at L * = 4 starts converging at L = 69 , as shown in Fig. 2. When in Eq. (39) the ratio 4f max L,j /A max L,j < −1 , the lower bound becomes complex and thus, invalid. When the  www.nature.com/scientificreports/ condition 4f max L,j /A max L,j > −1 is satisfied, an appropriate convergence can be obtained for the lower bounds within a few iterations. At L * = 2 the iSCLBT simultaneously calculates lower bounds for the two lowest lying energy levels. It can be seen that the iSCLBT results are a significant improvement over the lower bounds from the SCLBT method 50 even at L * = 2 . With the increase in L * , the iSCLBT results improve further and converge: the lower bound gaps for L * = 4 (gray line) and L * = 5 (green line) are very close to each other. The Heisenberg results for the ground-state lower bounds are significantly better than the Hubbard results, because the variance of the ground state Heisenberg eigenvalue is approximately three orders of magnitude smaller than that of the Hubbard system. While results from SCLBT 50 are significantly better in the Heisenberg case than those of the Hubbard model case, the iSCLBT results are more accurate for both models, because the variances for the excited states are comparable.
The lower bound gap ratios calculated by the iSCLBT and the Lehmann theory were compared as well, using both lattice systems. Both methods can provide lower bounds not only for the ground state, but also for higher excited states. However, the number of treatable excited states is limited by the Lánczos construct, which provides reliable eigenvalues for the lowest lying energy levels only. While the iSCLBT does not depend on the Weinstein lower bound explicitly, the Lehmann theory relies on it. Theoretically, the Weinstein lower bound is a monotonically increasing function of L where it satisfies the condition of validity in Eq. (18). However, in practical calculations, due to the fluctuations of the variances with the increasing basis set size, the Weinstein lower bound function is not monotonic so that one uses the largest Weinstein lower bound. Specifically, if for basis set dimension L + 1 the Weinstein lower bound is lower than for L, one uses the lower bound for L, etc., similar to the strategy employed in the SCLBT 50 . For the Lehmann theory, the corrected Weinstein lower bounds corresponding to L * were used for the Lehmann pole ρ in Eq. (51). This equation provides L * roots, which were calculated using the Mathematica software (Wolfram Research). The iSCLBT and Lehmann results were considered in the same ranges. As can be seen in Fig. 1, L,5 is the highest eigenvalue where the variance and the Weinstein lower bound calculated from it are valid. Thus, the highest achievable state is L * = 5 , that is, the ground state and the four lowest lying excited states. Figure 3 shows a comparison of the ground-state lower bound gap ratios for the Hubbard model calculated by the iSCLBT and Lehmann methods as a function of dimension L for different L * values, up to L * = 5 . The fluctuation of the gap ratios for both methods is due to the fluctuation in the variances. As indicated by the blue line, the iSCLBT gap ratios at L * = 2 are significantly better than the corresponding Lehmann results. The Lehmann gap ratio at L * = 2 is slightly better than the results of the SCLBT indicated by the violet line shown in the left panel of Fig. 2 and the shape of the two graphs is similar. This is because both of these methods are based on the Weinstein lower bounds. As indicated by the red line in Fig. 3, the lower bounds improve significantly at L * = 3 for both the iSCLBT and Lehmann calculations; in this case, the iSCLBT results are still superior to those of the Lehmann method. However, while the iSCLBT can provide good quality lower bounds at lower L values, the Lehmann method requires a larger basis set size to achieve similar performance. The Lehmann results are comparable to those of the iSCLBT above L ≃ 80 in this case. At L * = 4 and L * = 5 both methods converge: the Lehmann lower bounds are slightly better than those calculated by the iSCLBT. The converged gap ratios at L * = 5 are 1.880 for the iSCLBT and 1.640 for the Lehmann method. The ground-state results indicate that the increase in L * , the highest state to be considered, provides significant improvement in the lower bounds with a noticeable convergence for both lower bound methods, whose accuracy becomes similar.
A comparison of the first excited state lower bound gap ratios for the Hubbard model obtained from the iSCLBT and the Lehmann theory is shown in Fig. 4. In the left panel, the blue line indicates lower bound gap ratios from the iSCLBT at L * = 2 ; however, while the iSCLBT method can calculate lower bounds for both states (ground and first-excited), the quality of the roots provided by the Lehmann theory for the first excited state at  Figure 6 shows results for the third and fourth excited states for the Hubbard model. The quality of the Lehmann results was very low for the third excited state (left panel) even at L * = 4 . The quality of the corresponding Weinstein lower bound directly determines the quality of the Lehmann bound, as it explicitly depends on it. In this case the Weinstein bound of the fourth excited state, ε We L,4 has a local minimum at L = 90 and even when corrected as discussed above, it provides poor quality results. However, the iSCLBT can still provide good quality lower bounds regardless of the quality of the Weinstein lower bound, as shown in the left panel of Fig. 6. The convergence behavior of the iSCLBT result is similar to that of the previous cases, with a final gap ratio of 3.545. As can be seen in Fig. 1, the highest energy level for which lower bounds could be obtained for the Hubbard model was the fourth excited state. The right panel of Fig. 6 shows lower bound gap ratios for both methods at L * = 5 . Although the Weinstein lower bound ε We L,5 has a local minimum as well, for  As can be seen in Figs. 3, 4, 5 and 6, the increase in L * always results in an improvement of the lower bound for both methods. While the iSCLBT can provide lower bounds for all energy levels at a given L * reference level, the quality of the lower bounds calculated by the Lehmann method at the lowest L * value is extremely low; furthermore, the lower bounds are sharply increasing with the increase in dimensionality L. Generally, the Lehmann method can only provide L * − 1 appropriate roots for a calculation using L * states. The convergence behavior of the iSCLBT and Lehmann methods is similar in all cases. The Ritz gaps together with the final iSCLBT and Lehmann gap ratios for the Hubbard model are summarized in Table 1. The final values were calculated at L = 106, where the lower bound gaps were still numerically reliable. It can be seen that as one goes up the ladder of states all gaps and gap ratios increase. As mentioned before, the quality of the third excited state gap ratio calculated by the Lehmann method is very low.
The same calculations were performed for the Heisenberg model; the results can be found in the Supplementary Information. As shown in Fig. S1, the Lánczos construct provides less eigenvalues for the Heisenberg model; thus, the number of treatable energy levels is lower. The Heisenberg results further confirm the convergence behavior demonstrated in the Hubbard system. The ground-state results shown in Fig. S2 are better than those of the Hubbard model because, as discussed before, the Heisenberg eigenvalues have smaller variances for the ground state. Similar to the Hubbard case, the iSCLBT is noticeably better at L * = 2 and L * = 3 . The converged gap values are 1.625 for the iSCLBT and 1.250 for the Lehmann method. For the first excited state lower bound shown in Fig. S3, the iSCLBT at L * = 2 provides better gap ratios at low L values. This does not indicate that the lower bounds are better in this region; it rather shows that the Ritz gaps are still wide at low L values. From L = 70 the two methods exhibit similar convergence characteristics as that of the ground-state case: at L * = 3 the iSCLBT is better than the Lehmann result, while at higher L * values the Lehmann method is slightly better. The final lower bound gap values are 2.148 for the iSCLBT and 1.716 for the Lehmann method. The second and third excited state lower bounds can be seen in Figs. S4 and S5, respectively. The behavior and convergence of the iSCLBT and Lehmann results for the second excited state are similar to those of the Hubbard model. The converged lower bound gap ratios are 2.363 and 1.991 for the iSCLBT and Lehmann methods, respectively. In the Heisenberg model, the Lehmann method can provide reasonable lower bounds for the third excited state; however, these values are monotonically increasing with dimensionality L. As there is only one remaining value at L * = 5 , these results cannot be considered fully converged. The Ritz and lower bound gaps and final gap ratios   Table 2. It should be noted that for the fourth excited state only one point could be used, as the Lánczos construct provided only a few numerically reliable eigenvalues.

Discussion
In this paper, a further improvement of the existing SCLBT method was presented, termed iSCLBT. The lower bounds for the Ritz eigenvalues obtained from the iSCLBT were compared with those from the SCLBT and from the Lehmann methods. The framework provided by the Lánczos construct enabled the comparison of these methods on an equal formal footing. For the comparison, two model systems, a 5 × 6 Heisenberg and a 4 × 4 Hubbard lattice Hamiltonian were used. The eigenvalues and variances were determined by the Lánczos algorithm by using the H diagonalization software. By defining tighter bounds for the residual energy and for the diagonal elements of the overlap, the iSCLBT was improved over its previous implementation. First, the iSCLBT was compared with the SCLBT method for the ground-state energy. The improved theory provided significantly better lower bounds even at L * = 2 . Then, the lower bounds for the low-lying energy levels were calculated using the iSCLBT and the Lehmann method up to the fourth excited state. The effect of the highest considered level L * on the quality of the lower bounds was analyzed. Based on the analysis of the results, the following conclusions can be drawn: • The definition of tighter bounds for the residual energy and for the diagonal elements of the overlap further improved the SCLBT. • Both the iSCLBT and the Lehmann method can provide lower bounds that are significantly better than the Weinstein or Temple bounds. • The quality of lower bounds improves with the increase in L * , the highest considered level.
• The iSCLBT and Lehmann methods exhibit similar performance and convergence behavior. This is not unexpected, considering that in Ref. 42, using a finite Hamiltonian construct and assuming a Lánczos basis set, we could find conditions under which the two theories would be identical. However, in general, and under the conditions of the present theory which differs from the one presented in Ref. 42, the Lehmann theory and iSCLBT are formally different. • Compared to the iSCLBT, the Lehmann method requires a larger basis set and, as the Lehmann pole is estimated from the Weinstein lower bound, the quality of the Lehmann bounds is strongly affected by the quality of the Weinstein lower bound. • The numerical implementation of iSCLBT is simpler than that of the Lehmann method and does not require the Weinstein lower bounds.
Both the iSCLBT and Lehmann methods are suitable to provide high quality lower bounds for the low-lying energy levels for the studied lattice systems. The number of treatable Ritz eigenvalues is determined by the size of the Lánczos basis set.