Lower Bounds for Coulombic Systems

As of the writing of this paper, lower bounds are not a staple of quantum chemistry computations and for good reason. All previous attempts at applying lower bound theory to Coulombic systems led to lower bounds whose quality was inferior to the Ritz upper bounds so that their added value was minimal. Even our recent improvements upon Temple’s lower bound theory were limited to Lanczos basis sets and these are not available to atoms and molecules due to the Coulomb singularity. In the present paper, we overcome these problems by deriving a rather simple eigenvalue equation whose roots, under appropriate conditions, give lower bounds which are competitive with the Ritz upper bounds. The input for the theory is the Ritz eigenvalues and their variances; there is no need to compute the full matrix of the squared Hamiltonian. Along the way, we present a Cauchy–Schwartz inequality which underlies many aspects of lower bound theory. We also show that within the matrix Hamiltonian theory used here, the methods of Lehmann and our recent self-consistent lower bound theory (J. Chem. Phys.2020,115, 244110) are identical. Examples include implementation to the hydrogen and helium atoms.


INTRODUCTION
One of the outstanding challenges especially in ab initio quantum chemistry is obtaining lower bounds to atomic and molecular energies, which are as accurate as the upper bounds obtained with the Courant−Fischer theorem from the Ritz variational method. 1 Lower bound methods abound, starting with Temple's seminal expression derived in 1928. 2 Landmarks in the derivation of lower bounds are Weinstein's lower bound of 1934 3,4 and Lehmann's optimization of Temple's lower bound presented in 1949−50. 5,6 Especially Lehmann's expression has turned out to be quite accurate in different settings, however not so for Coulombic systems, as exemplified by computations on the He 7−10 and Li 11 atoms.
In the past few years, we have presented an improvement of Temple's formula for lower bounds of eigenvalues of Hermitian operators Ĥ1 2−15 and like Lehmann's theory it can become as accurate as the Ritz upper bound estimates. 16 However, these results were derived through explicit use of a Lanczos basis set, 17 which depends on a Krylov space, 18 in which one creates basis vectors by repeated application of the Hamiltonian operator. This does not work for unscreened Coulomb potentials 19 because the Coulomb singularity causes the third and higher moments of the Hamiltonian to diverge.
In this paper, we provide an answer to this challenge. Assuming an L dimensional subspace of the Hilbert space (as is standard in the Ritz variational theory), we present a lower bound expression, which depends only on the Ritz eigenvalues and associated variances. In contrast to Lehmann's method, there is no need to compute the matrix representing Ĥ2 in the chosen space. The theory utilizes Lehmann's approach as well as our recent results. To distinguish the present theory from previous ones, we will refer to it as the Pollak−Martinazzo (PM) lower bound theory using the abbreviation PM theory.
In Section 2, we review the known lower bound theories and their simplification when using Lanczos basis sets. Especially for Lehmann's theory, we derive a new and simplified expression of the Lehmann eigenvalues, which replaces the necessity of knowing the full Ĥ2 matrix with the need to know only the variances of the Ritz eigenstates which, in turn, thanks to the Lanczos construct, are available from the Hamiltonian matrix only. We then construct in Section 3 an L + 1-dimensional auxiliary Hamiltonian matrix whose diagonal elements are the L Ritz eigenvalues. The additional dimension is taken to be such that the L + 1 dimensional matrix is parametrically dependent on one of its eigenvalues, which can be set at will, and in practice is set formally to one of the true (unknown) eigenvalues of the Hamiltonian operator. In this construct, the standard deviations associated with the Ritz eigenvalues turn out to be coupling elements, which couple the approximate eigenfunctions associated with the Ritz eigenvalues to the exact eigenfunction associated with the chosen exact eigenvalue. In Sections 3.2 and 3.3, we proceed to show that the resulting eigenvalue expression is consistent with both the Lehmann and our previous selfconsistent lower bound theory when applied to the L + 1dimensional auxiliary problem. The reader interested only in the new lower bound theory and its implementation can skip these at a first reading.
To apply the theory in practice, it is necessary to determine a lower bound to one of the roots of the eigenvalue expression. We show that this is readily obtained by considering how the roots of the equation change with increasing dimensionality. As a result, there is a parallel to the so-called Lehmann pole 20 such that the Lehmann pole, which is employed in the "standard" Lehmann lower bound theory, turns out to be also the pole needed to construct lower bounds through the new lower bound expression. It cannot be over stressed that the "PM" theory derived in this paper depends only on the Ritz eigenvalues and their associated variances. Eigenfunctions associated with the Ritz eigenvalues are only needed for computation of the associated variances but there is no need to compute the full Ĥ2 matrix.
In Section 4, we apply the resulting theory to the He atom using a scaled Schrodinger basis set 25 and to the hydrogen atom using a Gauss Hermite basis set. The resulting lower bounds are superior to estimates based on Temple's expression and the "standard" Lehmann theory, which is not only less accurate but also considerably more expensive since it is based on computation of the L dimensional matrix of Ĥ2. We end with a discussion of the advantages, future challenges, and possible pitfalls in application of the new method to more complex atoms and molecules.

SHORT REVIEW OF PREVIOUS LOWER BOUND
THEORIES 2.1. Framework. We start with a Hamiltonian operator, whose eigenvalues and eigenstates are denoted as H j , 1, 2, ...
The eigenvalues are ordered in ascending order, that is, if j ≤ k, then ε j ≤ ε k . The ground state is given the index 1 (rather than 0) to simplify the notation later on. We also assume the existence of a known orthonormal basis set |Ψ j ⟩ = 1, 2, ... such that the Hamiltonian operator may be represented exactly as To simplify, we will assume that all functions and associated overlaps are real; however, this is not essential; the important property is that the operator under study is Hermitian. As in any practical computation, one never has the full Hamiltonian matrix (except for special cases) but rather its representation in a finite basis set, say the first L states spanning a space . Henceforth, in order to simplify notation, we will not indicate the dimensionality L and denote with the projector onto and with that onto its orthogonal complement.
The Hamiltonian projected onto the finite basis set is As discussed in ref 15, this assumption is somewhat less restrictive than the accepted condition for the validity of the Weinstein lower bound, 4 which is that the Ritz eigenvalue λ j is the closest one to the true eigenvalue ε j , that is, To derive the Temple lower bound, we define with each eigenstate in the projected space a "residual energy" λ̅ j such that With this definition, the residual energy may be expressed in terms of the overlaps and exact eigenvalues as where δ jk is the Kronecker delta. It is a matter of straightforward algebra to show from eq 2.14 that (2. 16) Inserting this identity into eq 2.11 and rearranging, one finds the Temple lower bound expression (2. 17) In the form of eq 2.17, the previous unknown overlap a jj has been replaced by the as yet unknown residual energy λ̅ j . However, we have gained something. Consider the ground-state residual energy. For any k ≥ 2, we have, through the initial ordering of the eigenvalues, the property ε k ≥ ε 2 so that a a (1 ) The Temple lower bound for the ground state now takes the well-known form where ε 2 − , which must be greater than λ 1 , is defined as a lower bound to the first excited-state energy. This lower bound may be obtained through a variety of lower bound methods such as the Weinstein, 3 Bazley, 21, 22 Miller, 23 and Marmorino 24 methods. Introduction of the residual energy made it possible to obtain a practical calculable form of Temple's lower bound formula.
If the Ritz eigenvalues converge to the exact energies when increasing the dimensionality of , so will the Weinstein and Temple lower bounds since the variance vanishes for exact eigenstates. However, the convergence will be much slower than the Ritz convergence due especially to the variances. As noted in ref 8, the integrand of a diagonal matrix element of the Hamiltonian squared is always positive so that all errors in the approximate wavefunction add up. This is not the case for the Hamiltonian itself, where positive and negative errors tend to cancel each other out, leading to more rapid convergence. It is this slow convergence of the Weinstein and Temple lower bounds, as exemplified by precise computations on the He atom, 8 which has hindered their usage.
2.3. Self-Consistent Lower Bound Theory. Instead of using the crude separation of the Hilbert space as in eq 2.10, one may use the projector onto the orthogonal complement to and insert it into the Cauchy−Schwartz inequality (eq 2.9). Upon rearranging, this leads to an improved lower bound inequality This result is superior to Temple's lower bound since the expression in square brackets is always greater than unity. However, the overlaps a jk are unknown. As shown in our previous papers, 14,15 the key to turning this result into a practical one is the use of a Lanczos 24) one finds the needed relation Inserting this into eq 2.20 gives a practical improved lower bound expression 26) One also notes that the same considerations lead to an improved Weinstein lower bound. Assuming as before that a jj ≥ 1/2 but using the projection operator one readily finds that where the second line is valid only when using a Lanczos basis set. The implementation of these results and their improved convergence properties have been discussed in some detail in refs 14 and 15. The central drawback is that the expression is derived by using the Lanczos basis set, which does not exist for Coulombic systems due to the Coulomb singularity. It is this challenge which is addressed in this paper.
2.4. Lehmann Theory. The Temple lower bound as expressed in eq 2.17 is based on a particular choice of a basis function, namely, the eigenfunction of the Ritz eigenvalue. Lehmann noticed that one may choose a better linear combination of states in the space by solving a generalized eigenproblem in this space. The Lehmann equation is, in a form that suits best our purposes, where κ is the Lehmann eigenvalue and |Ω⟩ ∈ is the associated Lehmann eigenfunction. The parameter ρ is known as the Lehmann pole and can be any real number but a Ritz eigenvalue for the above eigenproblem to be well-defined. However, for eq 2.28 to provide lower bounds τ = κ + ρ to the Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article first L* ≤ L lowest eigenvalues (as is customarily needed in quantum chemistry calculations), the sample space must be "good enough" such that λ L* ≤ ε L*+1 holds, and ρ must be limited by the condition λ L* ≤ ρ ≤ ε L*+1 . Only under such circumstances will eq 2.28 deliver L* negative Lehmann eigenvalues and these are lower bounds to the first L* eigenvalues of Ĥ. In practice, then, L* is the highest state for which the inequality λ L* ≤ ε L*+1 holds and ρ is a lower bound to ε L*+1 .
To understand the lower bound property of the Lehmann bounds (the τ's obtained from the Lehmann eigenvalues according to κ + ρ), it is useful to introduce if and only if |y⟩ ≡ |Y⟩ = (Ĥ− ρI)|Ω⟩ where |Ω⟩ satisfies eq 2.28, and in turn where κ −1 is the corresponding quotient (a Ritz eigenvalue of G(ρ)). Then, the Courant−Fischer theorem guarantees that the negative values κ −1 are upper bounds to the exact eigenvalues (ε k − ρ) −1 of G(ρ) from above for ε k lower than ρ; if the negative κ's are sorted in order of decreasing magnitude, |κ L* | ≤ |κ L*−1 | ≤ |κ 1 |, then τ n = ρ + κ n is a lower bound to the (L* − n + 1)th eigenvalue left of ρ, that is the lower bounds are ordered as τ k < ε k . To see the connection with Temple's lower bound, one multiplies eq 2.28 with the bra ⟨Ω| to find that The Ritz variational theorem which underlies the Lehmann construct, as in eq 2.31, shows that the Lehmann eigenfunction is the function that maximizes Temple's lower bound.
To summarize thus far, the Lehmann method builds on the matrices of Ĥ2 and Ĥin the space; diagonalization of eq 2.28 gives the lower bound eigenvalues. The condition that ρ ≤ ε L*+1 implies that one needs knowledge of a non-trivial lower bound to the state ε L*+1 , this could be a Weinstein-or a Bazley-related lower bound.
Interestingly, when using a Lanczos basis, one does not need to know the full Ĥ2 matrix in the projected space but only the variances σ k 2 associated with the respective Ritz eigenvalues. To see this, one multiplies eq 2.28 by the bra ⟨Φ k | to find Rearranging and summing over all k from 1 to L gives an eigenvalue equation, valid for the Lanczos construct and one notes expressly that the variances may be obtained from eq 2.23, that is, all the information is in matrix elements of the Hamiltonian only. The challenge then is to obtain similar results also in the case of Coulombic potentials where the Lanczos construct is not possible.
3. LOWER BOUNDS FOR COULOMBIC SYSTEMS 3.1. Hamiltonian Matrix Construct for Lower Bounds. The "Achilles heel" in the simplifications presented in the previous section is the need to create a Lanczos basis with the full Hamiltonian. Of course, for a finite Hamiltonian matrix representation, any power of the matrix is well defined and does not diverge. As before, in the projected space , we assume that the L-dimensional Hamiltonian matrix is diagonal, with known Ritz eigenvalues and associated variances. At this point, we do not discuss how these variances are computed. We then expand the diagonal Hamiltonian matrix with one additional row and column such that it takes the form i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z where ε is for the time being a parameter. Notice that, for the sake of clarity, we use a simplified notation for ; henceforth, it is understood to have the dimension (L + 1) × (L + 1) where L dim = . In this "auxiliary" Hamiltonian matrix, the standard deviations σ j couple the "Ritz states" to the added new state and, because of the Cauchy interlacing theorem, 26 its eigenvalues x k (k = 0, 1,..., L) are interlaced by the λ k 's, that is, x k−1 ≤ λ k ≤ x k . Among these L + 1 new eigenvalues, one will be the energy ε. To see this, we note that the eigenvalue equation is The expression in the square brackets has to vanish and this implies that and clearly one solution is x = ε. The other L eigenvalues 1 x j , j = 1, ..., L are the solutions of the remaining polynomial equation Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Notice the interesting symmetry: if x k is an eigenvalue of ( ) ε  other than ε, then ε is an eigenvalue of x ( ) k  other than x k . As we shall show below, eq 3.4 lies at the heart of PM theory. Lower estimates on the x k poles will lead under suitable conditions to lower bounds to the eigenvalue under consideration. For example, when interested in the ground-state energy ε 1 , we find that one should expect that the first root x 1 ≥ε 2 so that a lower bound on the first excited-state energy ε 2 will give a lower bound to the ground-state energy, provided, of course, that the Ritz eigenvalue for the ground state is lower than the lower bound for the first excited state. Note the formal similarity between eq 3.4 and the Lanczos-based Lehmann equation for the lower bound derived in eq 2.37. It is in this sense that PM theory generalizes Lehmann theory without the need to use a Lanczos basis set.
The eigenvalue eq 3.4 may also be rewritten as from which one finds that 6) or in other words the eigenvalues of ( ) ε  are increasing functions of the energy parameter ε.
Examination of eq 3.6 might suggest that this monotonicity is only in intervals since one encounters an infinity as either x or ε go through a Ritz eigenvalue. In reality, there is no discontinuity as ε comes close to a Ritz eigenvalue; the same will happen to all the roots x of eq 3.4 except one. The result is that the roots can be arranged to define functions x (k) (ε) that are continuous on the whole real axis except for a single pole singularity at a Ritz eigenvalue λ k and which are monotonically increasing in each connected sub-domain (−∞, λ k ), (λ k , +∞). This property is discussed further in detail in the Appendix where it is shown that the above-mentioned singularity is harmless for the method described below.
Monotonocity has far-reaching implications. Let us set ε = ε 1 and consider the limiting situation where the Ritz eigenvalues λ k approach the exact eigenenergies ε k and have thus vanishingly small variances σ k 2 . We then choose the lowest L eigenvalues to construct the space and the auxiliary matrix. In this limit, the matrix ( ) has L eigenvalues matching the Ritz values and one eigenvalue which diverges to +∞ since the (L + 1) th diagonal entry of the auxiliary matrix causes tr ( ) 1 ε  to diverge. Having chosen one of the eigenvalues of ( ) 1 ε  to be ε 1 , then, necessarily, all other roots of the eigenvalue eq 3.4 are such that x k (ε 1 ) → λ k+1 . Moreover, this occurs from below because of the interlacing theorem (λ k ≤ x k (ε 1 ) ≤ λ k+1 ). In addition, since with increasing accuracy each λ k+1 tends to ε k+1 from above the same holds for the roots x k (ε 1 ), that is, ε k+1 ≤ x k (ε 1 ) in this limit. Now, we can revert the argument. Suppose that the calculation is sufficiently converged such that λ k ≤ ε k+1 and that we know a rough lower bound (yet greater than λ k ) to the (k + 1)th exact energy, call it ε k+1 − . We know that there exists an ε such that ( ) ε  has ε k+1 − among its eigenvalues x k (ε). Indeed, thanks to the symmetry of eq 3.4, such an ε is just the lowest eigenvalue of ( ) andthis is the key pointsince x(ε) is monotonically increasing, ε is guaranteed to be left of (i.e., lower than) ε 1 . In other words, we have managed to convert a lower bound to the (k + 1)th energy into a lower bound to the groundstate energy. This is further shown in Figure 1 for the case k = 1. The question of whether this "map" produces a tight (hence useful) bound will be addressed numerically below, where it will be shown to be indeed the case.
The condition x k (ε 1 ) > ε k+1 deserves some further comments since it is the key to obtaining lower bounds to the ground state. For definiteness, let us focus on the case k = 1 considered in the applications of Section 4. As shown in the Appendix, it is the degree of convergence of the ground-state Ritz eigenvalue that determines the closeness of x 1 to λ 2 irrespective of whether or not λ 2 is close to ε 2 . Hence, if λ 1 is reasonably close to ε 1 , then x 1 should be close enough to λ 2 to guarantee that x 1 > ε 2 . As demonstrated for the computations of lower bounds for the He and H atoms, this property is verifiable by considering the dependence of x 1 on the dimensionality L of the computation. If, and this is the typical case, it is a monotonically decreasing function of the dimensionality, then due to the limit that ultimately λ 2 → ε 2 , the eigenvalue x 1 is guaranteed to lie above ε 2 . The fact that ε 1 is not known is not critical since the property of a monotonically decreasing value of x 1 with increasing dimensionality will hold for a range of ε values close to ε 1 . It is this added property that distinguishes PM theory from Lehmann theory when the Lanczos construction is not exploited or not possible. When the latter is used, the ordering x 1 (ε 1 ) > ε 2 always holds and the two theories become equivalent to each other. In fact, if we had x 1 (ε 1 ) < ε 2 , we could choose a pole ρ in eq 2.37 larger than x 1 and yet below ε 2 (x 1 ≤ ρ ≤ ε 2 ) and bound in this way the ground state from below. However, this is clearly impossible by virtue of monotonicity (eq 3.6) since ε 1 < y 1 (ρ) where y 1 is the inverse function of x 1 .
These same considerations can be generalized to excited states. For the sake of clarity, let us focus on the first excited state ε 2 and set ε = ε 2 in the auxiliary matrix. In the limiting situation considered above, the eigenvalues x k (ε 2 ) approach the corresponding Ritz values from below, but now, due to the interleaving theorem and our ordering of the eigenvalues, x 1 (ε 2 ) ≤ λ 1 and x k (ε 2 ) ≤ λ k+1 for 2 < k ≤ L. The lowest rootx 1 (ε 2 ) is a lower bound to the ground state since, according to the above, we know that ε 2 < x 1 (ε 1 ) and we have to move x 1 leftward to match ε 2 . The remaining eigenvalues, on the other hand, are upper bounds to the states higher than ε 2 for the very same Red dots indicate the positions of the exact energy levels ε k and black vertical bars that of the Ritz eigenvalues λ k . Note that x 1 (ε 2 − ) is a lower bound to ε 1 and x 1 (ε 1 ) is an upper bound to ε 2 . reasons given above. Hence, even a crude lower bound ε k+1 − can be converted into a lower bound to ε 2 − . As is usual when using a finite basis set, the quality of the Ritz eigenvalues deteriorates as one "goes up the eigenvalue ladder". There will be some value L* above which the interleaving property of the Ritz and exact eigenvalues is no longer valid. However, the upper bound quality of the roots of the auxiliary Hamiltonian remains up to L*, that is, x L* ≥ ε L* . In practice, then, one can use the highest index (k + 1 = L*) for which ε L* ≥ λ L*−1 holds and use a lower bound ε L* − (yet such that ε L* − ≥ λ L*−1 ) to obtain lower bounds to all the lower lying states. This is the analogue of the pole in Lehmann theory.
The practical implementation of the PM theory parallels the practical implementation of the Lehmann theory. One starts with a valid lower bound to the Lehmann pole. Then, one computes all roots of eq 3.4 which are below the lower bound to the Lehmann pole and these will be lower bounds to the respective eigenvalues.
Apart from the increased lower bound accuracy obtained through the PM method, we note that the computational expense may be lower than the effort involved in computing the "standard" Lehmann lower bound. For the Lehmann equation, one needs the full Ĥ2 matrix. For the PM method, one only needs the variances associated with the Ritz eigenfunctions. This implies that if |Φ j ⟩ is an eigenfunction, one can compute directly diagonal matrix elements of the sort ⟨Φ j |Ĥ2|Φ j ⟩ and there is no need to first compute the Ĥ2 matrix.
In the next section, we will give a numerical example which shows that the present theory gives improved lower bounds for the excited states as compared to any other existing method.

Lehmann Eigenvalue Equation for
the Hamiltonian of Equation 3.1. The result of the previous subsection has the flavor of the Lehmann lower bound and now we will show that indeed it is identical to it provided that one uses in the Lehmann eq 2.28 the matrix Hamiltonian ( ) ε  (eq 3.1) instead of the full Hamiltonian Ĥ. For this purpose, it is expedient to introduce an auxiliary vector |Ψ ⊥ ⟩ (which is only required to be orthogonal to the |Φ k ⟩ vectors, for k = 1, ..., L) and interpret the matrix ( ) ε  as the representation of an operator K̂(ε) acting in the enlarged, L + 1-dimensional space built with the |Φ k ⟩'s and |Ψ ⊥ ⟩. We then repeat the derivation of the Lehmann lower bound expression starting from the eigenvalue eq 2.28, replacing the full Hamiltonian Ĥwith the Hamiltonian of eq 3.1. By definition and The Lehmann eq 2.28 remains as before. Multiplying it by the bra ⟨Φ k | and rearranging gives We then use the notation Inserting this into left-hand side of eq 3.9 leads to the eigenvalue equation Comparing with eq 3.4 we identify κ + ρ with ε and ρ with x. This Lehmann equation is exact for the Hamiltonian K̂(ε) and is valid also for Coulombic systems. In other words, the PM method is identical to the Lehmann method for systems with Lanczos basis functions as may be inferred by comparing the PM equation and the Lehmann eigenvalue equation for Lanczos systems as in eq 2.37. If however one uses the Lehmann eq 2.28 for systems that cannot use the Lanczos construct such as Coulombic systems, then the PM method gives superior results as shall be exemplified below. However, the PM method is identical to Lehmann's equation provided that one uses it with the matrix Hamiltonian rather than the full Hamiltonian.
This then complements the previous proof of how one obtains lower bounds. Let us suppose that, as before, L* is the highest state for which the interleaving property of the Ritz and exact eigenvalues holds. We then use ε L* − as the Lehmann pole. The lower bound property of the Lehmann equation remains valid so that we know that Finally, before closing this subsection, let us mention two further results that are instrumental to the next one. First, we notice that in the enlarged space, the eigenvector of the auxiliary Hamiltonian K̂(ε) with eigenvalue ε takes the simple form as can be readily verified by inspection. Second, with some straightforward algebra, it is possible to recast the solutions for the eigenvalues other than ε as 3.15) and to write the associated eigenvectors in a form that closely parallels eq 3.14. We will make use of these expressions in the next section where we shed light on the relationship between eq 3.4 and our recent findings. 14−16 3.3. Improved Self-Consistent Lower Bound Theory Using the Hamiltonian of Equation 3.1. We may now also show that our previous improvements of Temple's theory as described in refs 14,16 and eq 2.26 are also identical to eq 3.4. Choosing the parameter ε to be the k-th eigenvalue of the Hamiltonian (Ĥ), the eigenfunction of the Hamiltonian K̂(ε) associated with this eigenvalue is denoted as |ε k ⟩ and the Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article eigenvectors associated with the remaining roots of eq 3.4 as |x j ⟩ (see eqs 3.14 and 3.16, respectively). Similar to the development in Section 2 of this paper, we may rewrite each Ritz eigenvalue as where the residual energy of the Hamiltonian K̂(ε) becomes We then note that for the extended Hilbert space of the Hamiltonian K̂(ε), we have the identity 3.20) so that with our construct  Putting this all into eq 3.21 and rearranging, we get the identity 3.24) and this has the same form as eq 2.26 without invoking a Lanczos basis set. To turn this into a practical expression, it is necessary to estimate the residual energy as defined in eq 3.18, which differs from the residual energy as defined in eq 2.15. For this purpose, one needs lower bounds to the roots x k and this could follow the same procedure as above using the improved Weinstein lower bounds, which now may be also derived by assuming as in Weinstein theory that ⟨ε k |Φ k ⟩ 2 ≥ 1/2 to find that To summarize this section, we have two main results. The first and most important one is the practical one. Given the "Lehmann pole", we obtain lower bounds to eigenvalues without the need to compute the full Ĥ2 matrix; all one needs are diagonal elements of it. Second, we have demonstrated the identity of the Lehmann lower bound expression with the selfconsistent Temple lower bound expression and both are in principle exact. It remains to show that this methodology gives lower bounds for Coulombic systems that are superior to those obtained through the "standard" Lehmann methods. This is demonstrated in the next section.  with r 1 , r 2 the distances of the electrons from the nucleus and the variational parameter α was chosen in all the computations as the value α = 27/16 which minimizes ⟨Ψ 1 |Ĥ|Ψ 1 ⟩. The basis set was constructed using the scaled Schrodinger approach of Nakatsuji. 25 The scaling function (with r 12 the distance between the electrons) was chosen to be g r r r The third function will then be gH gH gH 1 and one continues in this fashion to build up the basis set. Note explicitly that although constructed similarly to the Lanczos algorithm, this procedure does not lead to a tridiagonal representation of the Hamiltonian. In our computation, due to our use of Maple, we were limited to small dimensionality. Even for the seven-dimensional computation, the second lowest eigenvalue is higher than ε 3 so that the lower bound computation is limited to the ground state. The same holds true of course for the Lehmann lower bound where even at L = 7, one finds only one negative eigenvalue which gives the lower bound to the ground-state energy.
4.2. Lower Bounds to the Ground-State Energy of He. The lower-bound property is based on the observation that the solutions of eq 3.4 are monotonically decreasing functions of the dimensionality. This is shown for the specific case of the Helium atom with our chosen basis set in Figure 2 where we plot the dependence of x 1 (ε 1 ) on the dimensionality of the basis set as it changes from 3 to 7. The monotonic decrease is the same when changing the argument in the vicinity of the ground-state energy.
One does not need to know the exact ground-state energy to ascertain that the eigenvalue decreases with increasing dimensionality. As described in Section 3, it is the observation that the eigenvalue x 1 is "trapped" between the Ritz eigenvalue λ 2 and the exact eigenstate ε 2 , which allows us to replace x 1 in eq 3.4 with ε 2 or a lower bound to it such that the resulting lowest root of eq 3.4 will be a lower bound to the ground-state energy.
To test the new theory and compare it with Lehmann, we use the following strategy. For the "standard" Lehmann computation, we choose the best possible Lehmann pole: ε 2 = −2.14597405. This is of course an idealization, typically if one knows one state, one knows the other and certainly if ε 2 is known then so is ε 1 . In a "realistic" scenario, the Lehmann pole for the first excited state would be given by a Weinstein-or Bazley-type lower bound, but for the sake of understanding the new theory without adding in other sources of approximate values, we make this choice. Similarly, the value of ε 2 was used to obtain the Temple lower bound as well as the PM lower bound derived from x 1 (ε 2 ). The resulting lower bounds as well as the Ritz upper bounds are shown in Figure 3. One notes the essential improvement of the PM lower bound (solid blue line), which becomes competitive with the Ritz upper bound when the dimensionality reaches 7. This is further exemplified in Figure 4 where the gap ratio of the lower bound to the upper bound is plotted as a function of dimensionality for the Lehmann (upper blue dashed line) and PM (lower solid red line) lower bounds. At L = 7, the PM gap ratio is 1.03. One notes that in the range 3 ≤ L ≤ 7 the Ritz upper bound and the Temple and Lehmann lower bounds hardly change. The basis set we chose is not optimal in this sense, but critically, as the dimension is increased, the Ritz upper bounds to the excited states improve significantly, as may be seen for the first excited state in Figure 2. It is this improvement in the Ritz eigenvalues and variances of the excited states which leads to the significant improvement of the PM lower bound with dimensionality as shown in Figure 3.

H Atom.
The energy levels of the hydrogen atom are known analytically, yet it serves as a good "playground" for studying lower bounds for this simplest of Coulombic systems. Figure 2. Dependence of the eigenvalue x 1 (ε 1 ) on the dimensionality L of the basis set used for computation of the Helium atom energies (solid red line). Shown also is the second Ritz eigenvalue λ 2 as a function of L (long dashed purple line); it is always larger than x 1 (ε 1 ) and both slowly converge toward the energy ε 2 (shown as the horizontal dashed−dotted blue line) from above. This demonstrates that one may use the energy ε 2 as a lower bound to x 1 (ε 1 ) and thus obtain lower bounds to the ground-state energy. The ground-state wavefunction is 4.13) and the ground-state energy is 1 2 To test the lower bound expressions, we used the normalized antisymmetric harmonic oscillator basis set and H (k) (r) is the kth order Hermite polynomial. This allows us to readily set up the Hamiltonian and Hamiltonian squared matrices and so test the various lower bound theories.
As in the case of the Helium atom, we show in Figure 5 that x 1 (ε 1 ) is a monotonically decreasing function of the dimensionality, ultimately going down to ε 2 . Here, we plot both the difference between the second Ritz eigenvalue and x 1 (ε 1 ) (top, brown long dashed line) and the difference between x 1 (ε 1 ) and the second excited-state energy (−1/8) (bottom red solid line) as functions of the dimensionality. The dotted− dashed blue line is 0, which should be the limit of both lines with increasing dimensionality. For all values, the pole x 1 (ε) is in between the Ritz eigenvalue λ 2 and the exact energy ε 2 .
Then, we compute as a function of L the Temple, Lehmann, and PM lower bounds, in all of them using ε 2 as the Lehmann pole energy. The results are shown in Figure 6 and one notices that PM theory is again superior to the other lower bounds. In Figure 7, we plot the gap ratios for the Lehmann and PM lower bounds; at its worst, the PM gap ratio is ca. 23 and one sees that it improves significantly with the increasing dimensionality of the basis set, reaching 3.05 when L = 70.
With our choice of basis set, for L ≥ 12 the second eigenvalue has the property that ε 2 ≤ λ 2 ≤ ε 3 so that from L = 12, using ε 3 as the "Lehmann pole", one will get lower bounds for the ground and first excited state. This is shown in Figure 8 where the two horizontal lines are the ground and first excited-state energies while the lower dashed line is the PM lower bound to the ground state and the solid red line the PM lower bound to the second state. Comparing with Figure 6, one notes that the ground-state lower bound here is not as good as the one obtained with ε 2 as the Lehmann pole. The reason is quite clear; the Ritz eigenvalue λ 2 converges more rapidly than λ 3 so that the lower bound to x 2 (ε 1 ) given by ε 3 is worse than the lower bound of ε 2 and compared to x 1 (ε 2 ). However, as may be seen from the plot, one is getting a rather "decent" lower bound for the second state −0.1314734 at L = 70 as compared to −0.125. Not as good as the Ritz upper bound (−0.124108), the gap ratio at L = 70 is 36.45, but the PM lower bound for the first excited state is much  between the Ritz upper bound for the second state and the eigenvalue x 1 (ε 1 ) while the lower (red) solid line shows the difference x 1 (ε 1 ) − ε 2 between the eigenvalue x 1 (ε 1 ) and the exact second-state energy ε 2 . Both lines are always positive, demonstrating that the second-state energy is indeed a lower bound to the eigenvalue x 1 (ε 1 ).  However, there remain difficult challenges ahead. Obtaining a Ritz upper bound is easy and straightforward. All that is needed is to construct the Hamiltonian matrix and diagonalize it. Deriving lower bounds is more complex. Even within the present simplified framework, one still needs to verify that the x k 's are monotonically decreasing functions of the dimensionality of the basis set used and it is necessary to compute the standard deviations associated with the Ritz eigenvalues. This implies the need to compute not only the Ritz eigenvalues but also their eigenfunctions. Perhaps, and this will be considered in future computations, it is not necessary to obtain the full orthogonal diagonalization matrix but only those of the first few dozen Ritz states to construct the matrix Hamiltonian and still obtain "good" lower bounds. But it will still be necessary to obtain the associated eigenfunctions of these lowest lying states, increasing the numerical cost of the computations.

DISCUSSION
However, the real challenge is not obtaining the eigenfunctions. A critical element in the theory is obtaining an accurate lower bound to the so-called "Lehmann pole". In the present paper, we took the easy road, by using the known excited-state energies for the helium and hydrogen atoms. As stressed in the paper, the reason we did this was to provide a fair and unbiased comparison of the different lower bound theories. We saw that the present PM theory is superior to any other, yet the quality of the lower bound depends critically on the choice of the Lehmann pole. A "standard" methodology would be to use the Weinstein lower bound. For Hubbard-like Hamiltonians, we have shown 14 that this choice is sufficient for obtaining tight lower bounds. The same is true for atoms as long as one is considering the ground state. However, the rapid reduction of the level spacing between excited states, as is the case in hydrogen, helium, and lithium, 11 presents a serious challenge. Although we have shown how to improve upon the Weinstein lower bound, as for example in eq 2.27, the implementation is based on the assumption that the diagonal overlap matrix element squared is greater than 1/2. The challenge is to know when this assumption holds. In the "standard" Weinstein theory, one must consider the corresponding Ritz eigenvalue and know that it is the closest to the true eigenstate of the Hamiltonian under consideration. This is especially difficult when eigenstates come close together as for excited electronic states of atoms. Here, one must show that the Ritz eigenvalue is closest to the eigenvalue of the matrix Hamiltonian under consideration. These eigenvalues do not necessarily bunch together and this is an advantage; however, one does need to construct an objective criterion which would enable knowing that the overlap condition is valid, and this is not a trivial task.
Another challenge is understanding the choice of the basis set. Even the Ritz upper bound depends on the choice of the basis  lower (orange) line is the PM lower bound for the ground state, and the solid (red) line is the PM lower bound for the second state. The (blue) dashed−dotted line shows the second Ritz eigenvalue, which bounds the second state from above. As noted in the text, the lower bound shown here for the second state is superior to the same obtained from the "standard" Lehmann theory.
set. Different basis sets could give different Ritz upper bounds and PM lower bounds. For example, even in the present application of the theory to He, we made a specific choice of the exponent α in the initial wavefunction. In principle, for a given dimension of the basis set, one could vary α to minimize the Ritz eigenvalue for the ground state. One could also consider maximizing the PM lower bound for the ground state via variation of α. The two variations need not give the same value of the parameter. Different values imply different basis sets. The question of the "best" basis set remains open to both analytical as well as numerical research. We focus in this Appendix on the behavior of the eigenvalues of the auxiliary matrix ( ) ε  when ε approaches a Ritz eigenvalue λ s . We have already shown in the main text that ε is always an eigenvalue of the matrix; here, we show that when ε approaches λ s , the spectrum of the matrix ( ) ε  comprises the full set of Ritz eigenvalues plus a value tending to ±∞ for ε → λ s ∓ , respectively. The finite eigenvalues are smooth functions of ε that can be arranged to define solutions which are smooth and monotonically increasing on the whole real axis except for a pole at a Ritz eigenvalue. We use the pole position for their labeling, that is x (k) (ε), which is such a solution that is singular (only) for ε = λ k .
Let q ε (x) be the characteristic polynomial of the auxiliary matrix is the characteristic polynomial of the L × L upper left block of the matrix and p L (k) (x) = p L (x)/(λ k − x) is a polynomial of degree L − 1. Henceforth, the Ritz eigenvalues are supposed to be non-degenerate; otherwise, additional care is required in the analysis below but the result remains essentially unchanged. Furthermore, as in the main text, we use the notation λ k to denote the Ritz eigenvalues sorted in increasing order, λ 1 < λ 2 <···< λ L .
For x ≈ λ m , upon expanding p L (x) around x = λ m and noticing that the required derivative is −p L (m) (λ m ), one finds /σ s 2 > 0. Since they are always continuous when ε is not a Ritz eigenvalue and (as shown in the main text), they are also monotonic; they are in fact continuous functions that keep increasing when ε sweeps the real axis until they move beyond the highest eigenvalue, λ L . Suppose this happens for ε = λ k−1 , then for ε approaching the next Ritz eigenvalue λ k , the solution tends to +∞ and, for ε > λ k , it increases again starting from −∞. We call this solution x (k) (ε) using the pole position for its labeling. Note that the divergence occurs for one solution at a time since the interlacing property guarantees that for ε not a Ritz eigenvalue, there exists one (and only one) eigenvalue in each of the intervals (taking into account that ε is always an eigenvalue). Stated differently, we have the following theorem.

Theorem
If the Ritz values λ k are non-degenerate; the eigenvalues of the auxiliary matrix ( ) ε  other than ε can be arranged to define functions x (k) (ε) which are smooth and monotonically increasing on the whole real axis except for a pole at ε = λ k . Close to the latter, we have x (k) (ε) = x L (ε) for ε < λ k and x (k) (ε) = x 1 (ε) for ε > λ k , and beyond that point x (k) takes the next higher (lower) root every time ε crosses a Ritz eigenvalue to its right (left).
The functions defined in this way are monotonic in each connected domain, that is, separately for (−∞, λ k ) and (λ k , +∞), but this is harmless for the method described in the main text. To see this, we first notice that, thanks to the symmetry between x and ε mentioned in the main text, the inverse of the function x (k) (ε) is some other function x (l) (ε). The specific value of l can be identified by inspecting the limiting value that x (k) attains for ε → ±∞ since this gives the Ritz eigenvalue that makes the inverse singular (hence its label). We have x (k) (ε) → λ L+1−k ∓ for ε → ±∞, hence l = L + 1 − k. In other words, for any u v , ∈ , we have The situation is illustrated in Figure 9 where we have plotted the solutions x (k) 's for a model three-state problem, along with the eigenvalue ε (red line), as functions of ε. The figure makes clear that discontinuity is a marginal problem and arises only when one insists on labeling the eigenvalues according to their magnitude: each function takes a different x k in each interval (λ m , λ m+1 ), m = 1, 2,..., L − 1. Noteworthily, the symmetry between the x (k) 's and their inverse is apparent from the symmetry of the graph with respect to the main diagonal.
In the main text, in order to find a lower bound to the ground state, we have used k = L, hence the function x (L) and its inverse x (1) , as illustrated schematically in Figure 10. Notice that on the left of λ 1 x (L) (ε) is the lowest lying root (apart from ε) and the same happens for x (1) (ε) when ε ∈ (λ 1 , λ 2 ).    Figure 9 showing in detail the region [λ 1 , λ 2 ] relevant for calculating ground-state lower bounds according to the procedure described in the manuscript using a lower bound ε 2 − to the first excited state. The red stars mark the exact eigenenergies on both the x and the y axes (colors and symbols as in Figure 1 of the main text).