Acta Polytechnica https://doi.org/10.14311/AP.2022.62.0063 Acta Polytechnica 62(1):63–79, 2022 © 2022 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague ORTHONORMAL POLYNOMIAL PROJECTION QUANTIZATION: AN ALGEBRAIC EIGENENERGY BOUNDING METHOD Carlos R. Handy Texas Southern University, 3100 Cleburne St., Houston, Texas, U.S.A. correspondence: carlos.handy@tsu.edu Abstract. The ability to generate tight eigenenergy bounds for low dimension bosonic or ferminonic, hermitian or non-hermitian, Schrödinger operator problems is an important objective in the computation of quantum systems. Very few methods can simultaneously generate lower and upper bounds. One of these is the Eigenvalue Moment Method (EMM) originally introduced by Handy and Besssis, exploiting the use of semidefinite programming/nonlinear-convex optimization (SDP) techniques as applied to the positivity properties of the multidimensional bosonic ground state for a large class of important physical systems (i.e. those admitting a moments’ representation). A recent breakthrough has been achieved through another, simpler, moment representation based quantization formalism, the Orthonormal Polynomial Projection Quantization Bounding Method (OPPQ-BM). It is purely algebraic and does not require any SDP analysis. We discuss its essential structure in the context of several one dimensional examples. Keywords: Eigenvalue bounding methods, hermitian and non-hermitian linear operators. 1. Introduction The eigenvalue bounding problem for linear ordinary differential equations, or linear partial differential equations (i.e. LODE/LPDE) has been an active area of research for many decades. In the context of quantum physical systems, as represented by the multidimensional Schrödinger equation, − ℏ2 2m ∇2Ψ(−→r ) + V (−→r )Ψ(−→r ) = EΨ(−→r ), (1) the generation of upper bounds to the individual dis- crete state energies is readily obtainable through such well known methods as that of Rayleigh-Ritz (RR) [1]. The challenge has been to find an equally effective lower bound method. A well known lower bound method is that associated with Temple [2]; however, its convergence rate is slow. Nevertheless, it has served as a spring board for other more effective lower bound formulations. Among these is the work by Marmorino et al. [3], and more recently, that of Martinazzo and Pollak [4]. The latter are able to improve upon the convergence rate of Temple’s lower bound formulation. Despite these successes, two important facts remain. The first is that all the above bounding methods require the use of two different bounding formula- tions. One for the upper bounds (RR), another for the lower bounds. That is, they do not define a uni- fied theoretical framework for simultaneously gener- ating lower and upper bounds. Additionally, these methods are based on a Hilbert space representation for quantum systems, dependent on the existence of hermitian hamiltonians. They are of little relevance for bounding the real/complex eigenenergies of non- hermitian systems, particularly those corresponding to PT-symmetry breaking systems [5–8]. In this work we present a novel approach that can generate tight bounds for the discrete states of bosonic or fermionic, low dimension, systems, regard- less if they are hermitian or not. It is referred to as the Orthonormal Polynomial Projection Quantiza- tion Bounding Method (OPPQ-BM), as developed by Handy [9]; and based, in part, on a related method, the OPPQ-Approximation Method by Handy and Vrinceanu [10, 11]. Its general structure is outlined in the following sections, through representative one di- mensional systems, both hermitian and non-hermitian. We outline the full OPPQ-BM theory within the con- text of the one dimensional, double well, sextic anhar- monic oscillator; and then demonstrate the existence of the key structures necessary for its implementation to the PT symmetry breaking problem with potential V (x) = ix3 + iax. Beyond the theoretical interest in bounds, they are also of practical importance for delicate systems where conventional computational methods may yield widely varying results. That is, the availability of tight bounds allows one to discriminate between competing theories. One famous problem of this type corresponds to the Quadratic Zeeman (QZM) effect for superstrong mag- netic fields. This problem was analyzed through many different types of computational methods, resulting in a wide range of values for the most challenging state to compute: the ground state binding energy. This was reviewed by Le Guillou and Zinn-Justin (LG-ZJ) in the context of their order dependent conformal transformation analysis [12]. Using novel, Moment Problem [13] related, computational methods, Handy, Bessis, et al. [14, 15]. were able to confirm the accu- racy of the LG-ZJ analysis, by computing sufficiently 63 https://doi.org/10.14311/AP.2022.62.0063 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Carlos R. Handy Acta Polytechnica tight bounds to the ground state binding energy. Sub- sequent studies by Kravchenko et al. [16], and the more recent work by Schimerczek and Wunner [17], developed a different formulation that yielded vastly improved estimates (not bounds). The work by Handy extended OPPQ-BM to the QZM problem, yielding bounds that significantly improved upon, or were com- petitive with, the estimates by Kravchenko et al. The QZM problem is an example of an important class of problems for which bounding methods are highly relevant. These are classified as singular per- turbation – strongly coupled systems. Such systems involve quantum particles subjected to very strong forces, over relatively short length scales. This is in- herently a multiscale problem, in keeping with the objectives of wavelet analysis [18, 19], etc. Our original immersion into the eigenenergy bound- ing problem was through the study of strong coupling- singular perturbation type systems, such as QZM. The natural framework for regulating these systems is through the use of a non-local, extensive, representa- tion. The power moments provide such a representa- tion. As such, the types of systems studied here are those for which the Schrödinger operator configuration space problem can be transformed into a moments’ equation counterpart for the power moments of the bound state solutions. We refer to these as MER type systems (i.e. those admitting a moment equation representation). This will become clear below. This condition can be relaxed, and essentially imposed on systems that do not admit such MER formulations, making the underlying OPPQ-BM principles applica- ble to many different types of systems. The details of this expanded analysis will be communicated shortly. 1.1. Singular-perturbation strong coupling problems and moment representations It is not widely appreciated that power moments are relevant in understanding the multiscale structure of most systems. Thus, consider the scaling transform of a given wavefunction (i.e. signal, where we have simplified the notation to that of one dimension): SΨ(a, b) ≡ 1 aν ∫ +∞ −∞ dx S(x − b a )Ψ(x), (2) where Lima→0+SΨ(a, b) = Ψ(b), (3) and ν ≡ ∫ dx S(x) ̸= 0. Physicists are biased in favor of attaining an analyt- ical understanding of problems. Therefore the natural question to ask is, what is the analytical dependence in the inverse scale (i.e. 1 a ) for this scaling transform, if Ψ is a bounded, L2 state? It will become clear from Eq. (4), that the power moments of the bound state solution, determine the analytic structure of the inverse scale expansion. These considerations underlie the analysis by Handy [20]. Alternatively, engineers are more oriented towards computational capabilities. If so, then it readily fol- lows that the a → 0 limit can be replaced by the in- tegral ∫ ∞ 0 dα∂α ( αS ( α(x − b) )) , where α ≡ 1 a , which after a convolution substitution gives one the Contin- uous Wavelet Transform [19]. For bound state configurations, if lim|x|→∞ Ψ(|x|)S(|x|eiθ) = 0, exponentially, for arbitrary θ, then the scaling transform becomes ana- lytic in the inverse scale, 1 a . Under this assumption, the scaling transform’s analytic expansion depends on the moments [20] SΨ(a, b) = 1 aν ∞∑ j=0 σj j!aj j∑ p=0 ( j p ) (−b)j−pµ(p), (4) where ∂j xS(0) ≡ σj , and µ(p) ≡ ∫ +∞ −∞ dx xpΨ(x). (5) 2. The moment equation representation The natural extension of the above considerations is to study linear quantum systems whose differential form is transformable into a moment equation. Thus consider the sextic anharmonic oscillator potential, where the physical parameters have been re-scaled: −ϵ2∂2 xΨ(x) + (mx2 + gx6)Ψ(x) = EΨ(x). (6) The nature of physical quantum systems is that the discrete states decay exponentially, and therefore have finite power moments. The unphysical solutions be- come exponentially unbounded in one or both asymp- totic directions, therefore their power moments are infinite. We can multiply the above equation by xp and inte- grate by parts, assuming the underlying wavefunction is that of a discrete state. We then obtain the moment equation representation (MER): gµ(p + 6) = − mµ(p + 2) + Eµ(p) + p(p − 1)ϵ2µ(p − 2), p ≥ 0. (7) This homogeneous MER expression is a finite differ- ence equation of effective order 1 + ms where ms = 5. That is, for any E parameter value, the first six power moments {µ0, µ1, . . . , µ5} (i.e. µ(ℓ) ≡ µℓ) are the initialization moments, or missing moments, and gen- erate all the other power moments through closed form, energy dependent coefficients: µ(p) = ms∑ ℓ=0 ME(p, ℓ) µℓ, p ≥ 0. (8) 64 vol. 62 no. 1/2022 Algebraic Eigenenergy Bounding Method If the coupling strength is large, g >> 1, the nat- ural inclination is to attempt some kind of singular perturbation analysis involving expansions around the kinetic energy term; or, alternatively, a large pertur- bative expansion/resummation analysis. In configurations space, kinetic energy expansions become singular (i.e. expanding in ϵ2) because the order of the differential equation abruptly changes from zero to two. However, the order of the MER relation does not change as its kinetic energy counter- part is set to zero. This is one simple evidence that the MER transformation regulates singular perturba- tion expansions (i.e. kinetic energy expansions). That is, singular perturbation expansions in the moments’ representation are better behaved. 3. Generating eigenenergy bounds within a moments’ representation There are two methods for generating bounds within a moment equation representation (MER). The first method, referred to as the Eigenvalue Moment Method (EMM), was developed by Handy and Bessis [14, 15], and is based on the Moment Problem [13]. Its theoretical-computational structure is based on what is now referred to as semidefnite programming (SDP) [21, 22]. As such, the SDP based formulation of Handy and Bessis is the first use of such meth- ods for quantum operators [22]. Its computational implementation was done through the use of linear programming [23], since SDP algorithms were not known in the 1980s. The second moment representation bounding for- mulation is that presented in this work, OPPQ-BM. Unlike EMM, OPPQ-BM is applicable to the low lying discrete states of any system, hermitian, or non- hermitian, bosonic, or fermionic, provided it admits a moment equation representation (MER). EMM is ap- plicable only for the multidimensional bosonic ground state. For systems admitting both EMM and OPPQ-BM, it is our belief that at its basic level, EMM pro- duces faster converging bounds (as shown in this work through the analysis of the sextic anharmonic oscilla- tor in Eq. (6)); however, if one optimizes the selection of the reference/weight function, then OPPQ-BM can yield significantly faster converging results. The EMM formulation involves sophisticated, non- linear, convex optimization analytical tools. How- ever, OPPQ-BM is purely algebraic (i.e. eigenvalues, eigenvectors, and algebra). Given the power of Math- ematica, with unlimited precision, it can produce spectacular results. It is important to stress that EMM produces (when applicable) tight bounds ab initio. The OPPQ-BM is empirical. If the convergence of a particular parame- ter is numerically observed, then one can confidently generate bounds for the physical energies. We demon- strate this. 4. The eigenvalue moment method The bosonic ground state must be a positive (nonneg- ative) configuration [24], Ψgr(−→r ) ≥ 0. If the corre- sponding Schrödinger equation is transformable into MER form, then one can impose the Moment Problem positivity theorems and constrain the power moments, and in turn the energy and missing moments. This is the EMM-methodology. This was done, several decades ago, by Handy and Bessis (HB) [14, 15]. The Eigenvalue Moment Method (EMM), achieves geometric convergence rates for the bounds to the ground state energy. The only limitation is that it can only be applied to multidimensional bosonic systems, and then only to the ground state. As previously noted, it was used to solve the Quadratic Zeeman (QZM) problem [15]. The bounding of bosonic ex- cited states, through a moment representation, could only be realized more recently through application of OPPQ-BM [9]. One can extend the EMM quantization philosophy to the probability density of one dimensional hermitian Schrödinger operators. This is because the probability density will satisfy a third order, linear, differential equation (LODE). If the LODE representation ad- mits a MER formulation, then one can generate tight bounds to any discrete state. This is because all the discrete states are associated with nonnegative, L2, configurations. If the underlying Sturm-Liouville problem is non- hermitian, then the one dimensional Schrödinger equa- tion can be transformed into a fourth order LODE for the probability density. This approach was used by Handy [8] in precisely computing the a-parameter regimes where the system V (x) = ix3 + iax, violated PT symmetry (i.e. PT symmetry breaking). These re- sults were confirmed through a faster, moment based, estimation procedure [25], that led to a more detailed understanding for the onset of PT symmetry break- ing. These results significantly improved upon the predictions by Delabaere and Trinh [26], based on asymptotic analysis. We revisit this problem in this work. The multidimensional probability density, S(−→r ) ≡ Ψ∗(−→r )Ψ(−→r ) will not generally satisfy a linear partial differential equation; therefore, no MER relation can be generated, and EMM cannot be applied. Never- theless, through OPPQ-BM we can circumvent this difficulty. 65 Carlos R. Handy Acta Polytechnica 5. The orthonormal polynomial projection quantization formalism 5.1. The OPPQ non-orthogonal basis expansion Let us expand the discrete state wavefunction in terms of some appropriate, complete, non-orthonormal basis, {Bn(x) ≡ Pn(x)R(x)}: Ψ(x) = ∞∑ n=0 cnPn(x)R(x). (9) We require a weighted polynomial basis, involving orthonormal polynomials relative to some appropriate positive weight R(x) > 0 (we adopt the one dimen- sional notation for simplicity): ⟨Pm|R|Pn⟩ = δm,n. (10) The weight is real, as are its orthonormal polynomials, involving real polynomial coefficients: Pn(x) = n∑ j=0 Ξ(n) j xj . (11) For non-hermitian systems with complex bound state wavefunctions, it is the projection coefficients, {cn} that become complex. The basis {Bn(x)R(x)} will be complete, but non- orthogonal (i.e. ⟨Bm|Bn⟩ ≠ 0, if m ̸= n). We can rewrite Eq. (10) as m∑ j1=0 n∑ j2=0 Ξ(m) j1 ω(j1 + j2)Ξ(n) j2 = δm,n, (12) where ω(j1 + j2) = ∫ dx xj1+j2R(x), the Hankel mo- ment matrix of the weight, Wi,j ≡ ω(i + j). Knowledge of the Hankel moment matrix allows us to generate the orthonormal polynomials through the Cholesky decomposition method, which involves decomposing the positive Hankel matrix into the form W = CC†. Let êj correspond to a unit vector in the j-th component. We then solve for −→Ξ (j) = (C†)−1êj . (13) This generates the coefficient vector for Pn(x), or −−→ Ξ(n). The projection coefficients are obtainable from the MER relation for the power moments of Ψ(x): µ(p) ≡ ∫ ℜ dx xpΨ(x). (14) Assume that the corresponding MER relation exists µ(p) = ms∑ ℓ=0 ME(p, ℓ)µℓ. (15) It then follows that cn = ⟨Pn|Ψ⟩, = n∑ j=0 Ξ(n) j µ(j) = n∑ j=0 Ξ(n) j ms∑ ℓ=0 ME(j, ℓ)µℓ, (16) cn(−→µ ) = ms∑ ℓ=0 Λ(n) ℓ (E)µℓ, (17) where Λ(n) ℓ (E) = n∑ j=0 Ξ(n) j ME(j, ℓ). (18) The MER relation suggested in Eq. (15) is a ho- mogeneous relation for the power moments. Some normalization condition needs to be imposed: C(−→µ ) = 1. (19) For one dimensional systems, the natural normaliza- tion is the unit (nonlinear) normalization: |−→µ |2 = 1. However, it need not be chosen as such. Alternative choices [9] are linear normalizations such as µ0 = 1 or µ0 + µ1 = 1, etc. 5.2. The OPPQ quantization condition The OPPQ quantization condition requires that the weight be chosen so that the following positive inte- gral is bounded for discrete solutions and infinite for unphysical solutions: I[Ψ, R] = ∫ ℜ dx |Ψ(x)|2 R(x) , (20) = ∞∑ j=0 |cj(E, −→µ )|2, (21) where we have assumed Ψ is a bounded discrete state. More generally, we note that depending on the asymptotic behavior of the physical or unphysical (i.e. unbounded) solutions, and the chosen asymptoic form for the weight, the integral in Eq. (20) will satisfy the quantization condition: I[Ψ, R] = { finite, ⇐⇒ E = Ephys and −→µ = −→µ phys ∞, ⇐⇒ E ̸= Ephys or −→µ ̸= −→µ phys. (22) The OPPQ quantization condition essentially be- comes a shooting method in the E × −→µ , 1 + ms, pa- rameter space (after imposing a normalization). This is the essence of the OPPQ-BM bounding procedure. The focus of the remaining OPPQ formalism is to re- duce this shooting method to a minimization problem in the energy parameter space. 66 vol. 62 no. 1/2022 Algebraic Eigenenergy Bounding Method 5.2.1. Selection of the OPPQ weight The above formalism has two significant advantages. The first is that it tells us that the weight should not be chosen so that it asymptotically vanishes much faster than the asymptotic form of the physical solu- tions. Thus, given a physical, discrete state, Ψ, let C(x) satisfy the asymptotic relation: Lim|x|→∞ Ψ(x) R(x) = C(x). (23) We then want ∫ dx|Ψ(x)C(x)| < ∞. (24) In the works by Handy and Vrinceanu [10, 11], which introduce a particular version of the OPPQ formalism (i.e. to be referred to here to as the OPPQ- Approximation Method (OPPQ-AM)), they argued and demonstrated that the fastest convergence to the discrete states is associated with weights that mimic the asymptotic form of the desired discrete state: Lim|x|→∞C(x) = const. (25) The constant can be finite (i.e. zero), but not infinite in a manner that violates Eq. (24). 5.2.2. Use of the ground state as a weight Another advantage of the above formalism is that we can take R(x) = Ψgr(x) even if we do not know the functional form for the ground state. As long as one does not require that the discrete state wavefunction be reconstructed (i.e. one is only interested in the eigenenergies) then the only infor- mation required for the ground state is that its power moments be known accurately. This then allows us to generate the corresponding orthonormal polynomi- als, allowing for the generation of rapidly converging bounds to the discrete states. We demonstrate this in this work (i.e. Table 8). One can determine the power moments of the ground state wavefunction either through EMM or OPPQ-BM. Note that for bosonic systems, this is an excellent use of EMM, since the ground state is usually the only state that can be determined. Using the ground state as a weight usually yields the fastest convergence. 6. The OPPQ-approximation method (OPPQ-AM) From Eq. (21) and the OPPQ quantization condition in Eq. (22) it follows that the physical energy and missing moment values must satisfy lim→∞cn(Ephys, −→µ phys) = 0. (26) Since cn(E, −→µ ) = −→ Λ (n) (E) · −→µ , from Eq. (17), the 1 + ms linear equations ms∑ ℓ2=0 Λ(N−ℓ1) ℓ2 (E) µℓ2 = 0, (27) 0 ≤ ℓ1 ≤ ms, can be used to approximate the physical energies through the determinantal secular equation Det ( Λ(N−ℓ1) ℓ2 (E) ) = 0. (28) This defines the OPPQ-Approximation Method (OPPQ-AM). As indicated earlier, Handy and Vrinceanu [10, 11] noted significant improvement in the convergence rates when the weight, R, mimics the asymptotic form of the desired physical states. The OPPQ representation is very robust, and no convergence irregularities emerge so long as the weight does not decrease much faster than the asymptotic form of the physical states. Provided this is satisfied, the OPPQ-AM formalism will always converge to the true physical energies; however, the rate of conver- gence depends on the asymptotic properties of the chosen weight. Despite these impressive results, there is no guaran- tee that, for hermitian systems, the energies generated through Eq. (28) will be real. Spurious, small imagi- nary parts may be produced from Eq. (28), that vanish in the N → ∞ limit. By way of contrast, the OPPQ-Bounding Method to be described below, will always generate real energies for hermitian systems. In addition, as its name sug- gests, converging bounds to the discrete state energies can be generated. We make the last observation more explicit. As will be shown below, OPPQ-BM generates, to each or- der, an energy dependent function, LN (E). The local minima, ∂ELN (E) = 0 will approximate the physi- cal energies. These local minima can, essentially, be bounded, through converging lower and upper bounds. We refer to the local minima in LN (E) as the OPPQ- BM estimates, in order to distinguish them from the OPPQ-AM estimates from Eq. (28) and the OPPQ- BM generated bounds. 7. The OPPQ-bounding method We outline the structure of OPPQ-BM for one space dimension problems. The major difference between one dimensional and multidimensional MER type sys- tems is that one dimensional problems have a fixed number of missing moments: ms < ∞. Multidimensional problems involve an infinite hi- erarchy of missing moment subspaces of increasing dimension. That is, ms → ∞. We develop the 1- space dimension OPPQ formalism in a manner that extends to multidimensions. How the normalization prescription is chosen, plays an important role in the formalism [9]. 67 Carlos R. Handy Acta Polytechnica 5 10 15 20 E 35 Log 10 (�n(E)) Figure 1. Log10(λN (E)) for sextic anharmonic oscil- lator system in Eq. (47); N = 20, 40, 100, 120. The OPPQ quantization condition in Eq. (22) is dependent on the positive, (essentially) increasing, sequence, defined by the partial sums: I[Ψ, R] = LimN→∞SN (E, −→µ ), (29) where SN (E, −→µ ) = N∑ j=0 |cj(E, −→µ )|2, (30) = ⟨−→µ |PN (E)|−→µ ⟩, (31) where PN (E) is an energy dependent, positive matrix (if N ≥ ms) of dimension (1 + ms) × (1 + ms): PN (E) = N∑ j=0 (−→ Λ (j) (E) )∗−→ Λ (j) (E), (32) involving the sum over dyad matrix expressions. For non-hermitian systems, the “bra” missing mo- ment vector in Eq. (31) requires the complex conju- gate expression for positive norms on complex vector spaces. It trivially follows, by definition, that 0 < SN (E, −→µ ) < SN+1(E, −→µ ) < . . . < I[E, −→µ ]. (33) The OPPQ quantization condition in Eq. (22) tells us that the physical energy and corresponding missing moments correspond to (Ephys, −→µ phys) points within the E × −→µ parameter space where the functional I[Ψ, R] has a local minimum. Also, for fixed Ephys, the corresponding physical missing moment values are those corresponding to a global minimum in the missing moment space. Therefore, to order N we can focus on the global minimum within the constrained (i.e. normalized) missing moment space: LN (E) ≡ Inf−→µ {SN (E, −→µ )|Cnorm(−→µ ) = 1}, (34) where some convenient normalization has been adopted, C(−→µ ) = 1. 5 10 15 20 E 1 2 3 4 5 Log 10 (�n(E)) Figure 2. Log10(λN (E)) for sextic anharmonic oscil- lator even states in Eq. (65); N = 4, 6, 8, 10, 12; based on the weight RL(ξ) = ξ− 1 2 exp(−ξ). 5 10 15 20 E 2 4 6 8 10 Log 10 (�n(E)) Figure 3. Log10(λN (E)) for sextic anharmonic oscil- lator even states in Eq. (65); N = 4, 6, 8, 10, 12; based on the weight RG(ξ) = ξ− 1 2 exp(−ξ2/2). An important result is that LN (E) < LN+1(E). (35) This trivially follows from Eqs. (33-34). An immediate consequence is that the counterpart to the quantization condition in Eq. (22) now becomes simpler: lim N→∞ LN (E) = { finite ⇐⇒ E = Ephys, ∞ ⇐⇒ E ̸= Ephys. (36) Combining this with Eq. (35) we obtain: 0 < LN (E) 0} < B(U) phys. We then deter- mine the energy interval whose endpoints satisfy λN (E(L) phys;N ) = B(U) phys = λN (E(U) phys;N ). These become the lower and upper bound estimates for that physical energy. The above bounding analysis is initiated in Table 4 for the ground and second excited state of the sextic anharmonic oscillator. Since the use of the RA weight yields very rapid convergence, we see that the coarse upper bounds, B(U) are easily determined. Using these we can generate the bounds for the ground state and first excited state, as given in Table 3. We note that the OPPQ-BM bounds generated in Table 3 (we could have continued tightening the bounds), for the ground state, used 125 power moments giving us bounds at the 18th decimal place. In Eq. (69) we quote the EMM bounds obtained on the basis of 61 power moments. The accuracy is at the 33rd decimal place. However, in Eq. (70), using an optimal (Stieljes representation) weight, the OPPQ-AM converged to 50 decimal places, using only 45 power moments. 71 Carlos R. Handy Acta Polytechnica N E (L) 0;N E (U) 0;N E (L) 2;N E (U) 2;N 25 -0.524852943468474 -0.521778245307431 5.368449680815753 5.384374689469648 50 -0.523268749143883 -0.523268495112775 5.374969291406036 5.374970726288986 75 -0.523268622134517 -0.523268622120587 5.374970008804361 5.374970008875728 100 -0.523268622127553 -0.523268622127552 5.3749700088400432 5.3749700088400468 125 -0.52326862212755223948 -0.52326862212755223934 5.3749700088400449937 5.3749700088400449945 B(U) = −0.876 B(U) = −0.990 Table 3. OPPQ-BM upper and lower bounds for the ground and second excited states, using RA = e−x4/4. N ∂Eλ(E(min) 0;N ) = 0 Log10 ( λ(E(min) 0;N ) ) ∂Eλ(E(min) 2;N ) = 0 Log10 ( λ(E(min) 2;N ) ) 25 -0.523315367444853 -0.877054685910968 5.37640170043752 -0.991603081677466 50 -0.523268622128326 -0.877042389134943 5.37497000884746 -0.991530473233468 75 -0.523268622127552 -0.877042389134815 5.37497000884004 -0.991530473233135 100 -0.523268622127552 -0.877042389134815 5.37497000884004 -0.991530473233135 B(U) = −0.876 B(U) = −0.990 Table 4. OPPQ-BM Analysis for determining coarse upper bounds, B(U), for the ground and second excited states, using RA = e−x4/4 . 9.2. EMM-Ψ2 Another MER representation for the sextic anhar- monic oscillator is possible by working with the prob- ability density, S(x) = Ψ2(x). It is easy to show that the probability density for real potentials satisfies a third order LODE. For the sextic anharmonic os- cillator problem, this 3rd order LODE yields a MER relation of order 3 (i.e. ms = 2). We do not give the details of this analysis, since the following MER representation offers the easiest OPPQ implementa- tion. The value of EMM-Ψ2 is that we can bound the discrete states and use the results to gauge the effectiveness of OPPQ. This is referenced below. 9.3. EMM e− x4 4 Ψ(x) The third MER representation is obtained through the contact transformation, Φ(x) = exp(−x4 4 )Ψ(x). (63) Since this involves a factor identical to the dominant WKB asymptotic form for the physical states, the MER representation for Φ will involve fewer missing moments (i.e. none, after a normalization) than the MER for Ψ, as given in Eq. (48), involving five (5) missing moments (after imposing a normalization). This is desirable since the lower the missing moment order, the faster the convergence of either OPPQ or EMM. We provide the details of both approaches below. 9.4. OPPQ analysis of the (ms = 0) sextic anharmonic double well oscillator The double well anharmonic problem of interest is that for the potential V (x) = x6 − 4x2, where x ∈ ℜ. The physical solutions must die off, asymptotically, according to the dominant WKB expression Ψ(x) ∼ exp(− x4 4 ). If we work with the contact transformation in Eq. (63) Φ(x) = exp(− x4 4 )Ψ(x), we note that the discrete states remain normalizable and exponentially bounded, in the Φ representation. Unphysical Ψ con- figurations (i.e. non-normalizable due to their expo- nentially unbounded form in one or both asymptotic directions) map into non-normalizable Φ configura- tions. The EMM formalism works in either represen- tation, precisely because of this. We note that the power moments for exponentially bounded configura- tions exist; whereas they become infinite (or do not exist) for unphysical configurations. The Φ configurations must satisfy the differential equation Φ′′(x) + 2x3Φ′(x) + (7x2 + E)Φ(x) = 0. (64) The Hamburger moment (i.e. µ(p) ≡ ∫ ℜ dx xpΦ(x)) equation becomes (2p − 1)µ(p + 2) = p(p − 1)µ(p − 2) + Eµ(p). The even parity states will admit a MER for the even order power moments, µ(2ρ) = u(ρ) or (4ρ − 1)u(ρ + 1) = 2ρ(2ρ − 1)u(ρ − 1) + Eu(ρ). (65) The odd parity states µ(2ρ + 1) ≡ ν(ρ) will satisfy the MER : (4ρ + 1)ν(ρ + 1) = 2ρ(2ρ + 1)ν(ρ − 1) + Eν(ρ). (66) We note that the even order Hamburger moments satisfy: µ(2ρ) = u(ρ) = ∫ +∞ 0 dξ ξρ−1/2Φ( √ ξ), where x2 = ξ; however, the odd order Hamburger moments satisfy µ(2ρ + 1) = ν(ρ) = ∫ +∞ 0 dξ ξρΦ( √ ξ). The importance of these relations is that the respective 72 vol. 62 no. 1/2022 Algebraic Eigenenergy Bounding Method N E0 E2 E4 E6 E8 10 -0.5234534028 5.354166238 16.82887249 20 -0.5232677208 5.375031586 16.79298386 31.63171200 49.94322409 30 -0.5232686293 5.374969341 16.79536744 31.74329539 49.45632020 40 -0.5232686220 5.374970019 16.79534604 31.74254727 49.51247497 50 -0.5232686221 5.374970009 16.79534686 31.74254902 49.51148774 60 -0.5232686221 5.374970009 16.79534683 31.74254987 49.51149946 70 -0.5232686221 5.374970009 16.79534683 31.74254984 49.51150043 80 -0.5232686221 5.374970009 16.79534683 31.74254984 49.51150038 90 -0.5232686221 5.374970009 16.79534683 31.74254984 49.51150038 100 -0.5232686221 5.374970009 16.79534683 31.74254984 49.51150038 Table 5. OPPQ-BM Estimates (i.e. ∂EλN (E(min) phys;N ) = 0), Φ representation, (Eq. 65) RL(ξ) = ξ− 1 2 exp(−ξ). moments, in the context of Eqs. (65, 66) allow us to implement a Stieltjes moment analysis through EMM. We will only apply OPPQ on the even Ψ configu- rations, for simplicity. However, since the even order power moments are the moments of Υ(ξ) ≡ Φ( √ ξ)√ ξ , the OPPQ expansion must be relative to this con- figuration. Thus, the relevant OPPQ expansion will be Υ(ξ) = ∞∑ n=0 cnPn(ξ)R(ξ), (67) where Φ’s power moments will satisfy Eq. (65). Since the physical configurations behave, asymptot- ically, as Ψ(x) ∼ exp(− x4 4 ), the transformed expres- sions behave as Φ(ξ) ∼ exp(− ξ2 2 ). Since the transformed system involves a Stieltjes configuration supported on the nonnegative axis, one might take the weight to be the exponential function, R̃L(ξ) = exp(−ξ), with Laguerre polynomials; or the Gaussian R̃G(ξ) = exp(−ξ2/2), restricted to the non- negative real axis. However, for the quantization integral in Eq. (20) to apply, particularly with regards to the generation of bounds, we need to take into account the ξ− 1 2 that is inherent to the transformed, even order, power moments. Based on the previous arguments we will work with RL(ξ) = ξ− 1 2 exp(−ξ) and RG(ξ) = ξ− 1 2 exp(−ξ2/2). The corresponding power moments for the weights, in order to generate the orthonormal polynomials are obtained as follows. The power moments wL(p) ≡ ∫ ∞ 0 dξ ξpRL(ξ) or wL(p) = Γ[p + 1 2 ] satisfy the recursion relation wL(p + 1) = (p + 1/2)wL(p); whereas wG(p) ≡∫ ∞ 0 dξ ξpRG(ξ), or wG[p] = 2 2p−3 4 Γ[ 2p+1 4 ], satisfy the recursion relation wG(p + 2) = (p + 1/2)wG(p). In their original work on the OPPQ-Approximation Method (i.e. Eq. (28)), Handy and Vrinceanu imple- mented OPPQ-AM with both of these types of weights (i.e. using their counterparts along the entire real axis, x ∈ ℜ). They showed that superior (faster converging) results were obtained for weights that emulated the asymptotic form of the physical states. The enhanced efficiency of working with RG(ξ) in- stead of RL(ξ) is evident in Figure 3 compared to Figure 2 (i.e. the functions for the former are converg- ing much faster around the physical energies). The numerical results for the ms = 0 sextic an- harmonic oscillator problem are given in the follow- ing section. We only give numerical results based on using RL(ξ). Instead of working with RG(ξ), we will use the actual ground state wavefunction (i.e. its power moments), Rgr(ξ). The reason is that both have the same asymptotic behavior. That is, Rgr(ξ) = Φgr(ξ) ∼ RG(ξ); therefore, upon solving for the power moments of the ground state wavefunction (including the energy), we can use it to generate its orthonormal polynomials, and quantize the excited states, through OPPQ. To be able to do this requires sufficient accuracy in the determination of the ground state power moments, since these must generate a positive Hankel matrix for the weight (before implementation of a Cholesky analysis). Although this can be done through OPPQ, EMM inherently works with positive matrices and therefore is more efficient for doing this type of analy- sis. 9.5. OPPQ results for the ms = 0, sextic anharmonic oscillator The OPPQ-Φ formulation is a zero missing moment problem (i.e. ms = 0) and converges very fast. In Figure 2 we plot the nested sequence λN (E), for N = 4, 6, 8, 10, 12. This calculation was done based on the RL(ξ) weight. A similar result holds for the RG(ξ) weight as depicted in Figure 3. It clearly reveals the faster convergence afforded by a weight that emulates the asymptotic form of the physical states. In Table 5 we give the OPPQ-BM energy esti- mates (i.e. the local minima, ∂EλN (E(min) phys.;N ) = 0) for the low lying, even parity, discrete state ener- gies, based on Eq. (65), using the inferior weight RL(ξ) = ξ− 1 2 exp(−ξ). Table 5 only cites results (to ten significant figures) for N ≤ 100. Not shown in Table 5 are the OPPQ-BM estimation results for the ground state at N = 120: 73 Carlos R. Handy Acta Polytechnica N ∂ESN (E(min) N ) = 0 Log(SN (E(min) N )) cN (EN ) = 0 E (L) N E (U) N 20 49.94322409213197 7.982552803 49.4426432537 48.893 50.857 30 49.45632020275004 7.997267757 49.4986057623 49.366 49.546 40 49.51247497421839 8.002499541 49.5123070944 49.504 49.520 50 49.51148774312169 8.002657778 49.5114731270 49.510 49.512 60 49.51149945726214 8.002665930 49.5115007312 49.51140 49.51159 70 49.51150043370514 8.002667120 49.5115003992 49.51148 49.51151 80 49.51150037650008 8.002667337 49.511500375879294 49.5114990 49.5115017 90 49.51150037729806 8.002667346 49.511500377394135 49.5115001949 49.5115005597 100 49.51150037738496 8.002667348 49.511500377382972 49.5115003518 49.5115004029 110 49.51150037737999 8.002667348 49.511500377379789 49.5115003736 49.5115003811 120 49.51150037737991 8.002667348 49.511500377379918 49.5115003768 49.5115003780 B = 8.01 Table 6. OPPQ Results for E8: RL = e−ξ√ ξ . Egr = −.52326862212755224. This result lies within the bounds generated through an EMM-Φ analysis (based on the first 29 power moments): −.5232686221275616 < Egr < −.5232686221275495. (68) We can infer that the EMM bounds in Eq. (68) pre- dict the exact answer to 13 decimal places. This result, based on 29 moments, surpasses the OPPQ-BM esti- mation analysis in Table 5 (N = 30) of −0.5232712343, at the fourth-fifth decimal place. Thus EMM-Φ, by this comparison, is approximately three times more efficient. Having said this, we remind the reader that whereas EMM involves sophisticated analysis (i.e. nonlinear convex optimization), OPPQ is purely algebraic, and implementable to arbitrary accuracy through algebraic software, such as Mathematica. Additionally, EMM is, in practice, only applicable to the (multidimensional) bosonic ground state; whereas OPPQ applies to any multidimensional bosonic or fermionic (low dimension) system, including non-hermitian systems. 9.5.1. Generating bounds for E8 We re-examine the results in Table 5 for E8, the slow- est converging energy. As with other methods, such as Rayleigh-Ritz (RR), the manifest convergence of the results (i.e. the N → 120 sequence in Table 5) is no guarantee of the accuracy of the apparent limit in predicting the true energy. Whereas EMM produces converging bounds from first principles; OPPQ defines a procedure by which one can generate bounds pro- vided a certain parameter is empirically determined, specifically the coarse upper bound, B(U). Below we describe the OPPQ-BM bounding procedure in detail, although the same underlies the results in Tables 4 and 3, for the ms = 3 MER formulation in Eq. (48). The only advantage of the current problem is that it is a zero missing moment problem, and therefore easier to implement. In Table 6 we provide the OPPQ-BM eigenenergy estimate for E8 (the second column). The third col- umn contains the increasing, convergent, positive se- quence from which a coarse upper bound, B(U), is to be empirically determined. Using this coarse upper bound, we can generate arbitrarily tight bounds (i.e. the last two columns), as N → ∞. Note that we continued generating these bounds in Table 7. This entire procedure is based on the assumption that the manifest convergence of the third column in Table 6 is correctly bounded, from above, by the empirically determined coarse upper bound, B(U). The fourth column contains the OPPQ-AM estimate. It is important to appreciate that the coarseness of B(U) has nothing to do with the tightness of the bounds, in principle, assuming one can generate high OPPQ expansion orders. Below we compare the accuracy of the OPPQ- BM formalism to order N − 120, both with re- spect to the OPPQ-BM energy estimate (E8 = 49.5115003774, from Table 6) and the OPPQ-BM bounds (49.5115003768 < E8 < 49.5115003780, also from Table 6), as compared to the bounds generated by an EMM analysis. The results are very good and consistent. It is important to note that the OPPQ-BM estimate obtained at lower order, does not have to lie within the OPPQ-BM bounds obtained at higher order. Thus the entry in the second column in Table 6, corre- sponding to OPPQ-BM energy estimate (i.e. the local minima) of E = 49.51148774312169 (i.e. N = 50), lies outside of the bounds generated in the last two columns, for N > 80. However, all OPPQ-BM esti- mates must lie within the bounds calculated at lower order. 9.6. Comparison with EMM bounds Generally, the EMM analysis will be more efficient in generating bounds than OPPQ-BM. That is, fewer moments will be required to generate the same level of tightness of the bounds. However, this depends 74 vol. 62 no. 1/2022 Algebraic Eigenenergy Bounding Method N E (L) 8 E (U) 8 150 49.511500377377302482 49.511500377382545933 160 49.511500377379459378 49.511500377380389040 170 49.511500377379839286 49.511500377380009132 180 49.511500377379908253 49.511500377379940166 190 49.511500377379921131 49.511500377379927287 200 49.511500377379923601 49.511500377379924818 210 49.511500377379924086 49.511500377379924333 220 49.511500377379924184 49.511500377379924235 230 49.511500377379924204 49.511500377379924215 240 49.511500377379924208 49.511500377379924210 250 49.511500377379924209048830 49.511500377379924209558953 260 49.511500377379924209246851 49.511500377379924209360932 270 49.511500377379924209290914 49.511500377379924209316869 280 49.511500377379924209300890 49.511500377379924209306893 290 49.511500377379924209303186 49.511500377379924209304597 300 49.511500377379924209303723 49.511500377379924209304060 B = 8.01 Table 7. Bounds for E8 based on chosen B(U). N E0 E2 E4 E6 E8 1 -0.5232686221275529 5 -0.5232686221276021 5.37847969429 16.78604192264 10 -0.5232686221275522 5.37497046837 16.79542276464 31.74535880954 49.5248113799 20 -0.5232686221275523 5.37497000884 16.79534683270 31.74254983720 49.5115003841 30 -0.5232686221275523 5.37497000884 16.79534683270 31.74254983711 49.5115003772 Table 8. OPPQ-BM Estimates (i.e. ∂ESN (E) = 0) for V (x) = x6 − 4x2, ms = 0, R(x) = Φgr(x) based on the first 61, EMM generated, ground state power moments {ugr(p ≤ 60)}. on the choice of weight, as the following case exem- plifies. Thus, with regards to the ms = 2 missing moment formulation given earlier, the EMM bounds for the ground state −0.523268623844284 < Egr < −0.523268619253327, were based on the first 29 power moments. The same (approximately) level of tight- ness was achieved with OPPQ-BM using more than 50 power moments, based on the results in Table 3. Thus, in this example EMM is vastly superior to OPPQ-BM. However, for the next example the situation signifi- cantly improves for OPPQ-BM over EMM. For the ms = 0 formulation being considered, using the first 62 Stieljes power moments (i.e. {u(p)|p ≤ 61}), EMM achieves the bounds: − .52326862212755223941616949719078449 < EEMM < −.52326862212755223941616949719078395 (69) Using the weight RL the OPPQ-AM estimate achieves this level of accuracy on the basis of approximately 220 power moments; whereas the weight RG gener- ates an OPPQ-AM estimate that surpasses the EMM accuracy based only on the use of 45 moments (OPPQ- AM): EOP P Q−AM = −0.523268622127552239416169497 19078406116564771630604 (70) 9.7. Using the ground state as a weight In Table 8 we implement OPPQ-BM on a representa- tion that uses the generated moments of the unknown ground state R(x) = Ψgr(x). As expected, the con- vergence is very fast in comparison with results based on RL. 9.8. EMM results for the probability density Application of EMM to the probability density (i.e. EMM-Ψ2, an ms = 2 problem) yields the bounds −.52326866 < E0 < −.52326857, 1.0057681 < E1 < 1.0057685, 5.3749699 < E2 < 5.3749701, 10.5725845 < E3 < 10.5725855, 16.795339 < E4 < 16.795351, 23.883886 < E5 < 23.883961, 31.74217 < E6 < 31.74323, 40.301 < E7 < 40.305, 49.506 < E8 < 49.533 all to the same moment expansion or- der of 28. We also note that an EMM-Ψ{ms = 2} bounding formulation on the ground state yields com- parable bounds using the first 28 power moments: −.5232686237 < Egr < −.5232686193 . With regards 75 Carlos R. Handy Acta Polytechnica to E8, extension of EMM-Ψ2{ms = 2} analysis to the first 48 moments (i.e. {u(p)|0 ≤ p ≤ 47}) yields the bounds 49.5115003768 < E8 < 49.5115003798. These serve to further confirm some of the results previously cited. 10. A PT-symmetry breaking, non-hermitian, problem Consider the non-hermitian system −∂2 xΨ(x) + (ix3 + iax)Ψ(x) = EΨ(x), (71) which is known to break PT symmetry for negative values of the ‘a’ parameter (a < −2.611809356, and for a < −5.375879629) from Refs. [8, 25]. The case a = 0 has purely real eigen-energies, as computed by Bender and Boettcher [5], and theoretically confirmed by Dorey et al. [6]. The EMM-Ψ∗Ψ formulation also provided strong numerical evidence for this [27]. States that are PT symmetric satisfy Ψ∗(−x) = Ψ(x), and have real energies, E ∈ ℜ. States that break this symmetry have complex conjugate pairs for the energies (i.e. Ψ∗(−x) ̸= Ψ(x), with complex energies E∗ and E, respectively). We work with the MER-Ψ representation for which the complex power moments on the real axis satisfy µ(p + 3) = − aµ(p + 1) − iEµ(p) − ip(p − 1)µ(p − 2), p ≥ 0. (72) The recursion relation for the energy dependent gen- erator coefficients is ME(p + 3, ℓ) = − aME(p + 1, ℓ) − iEME(p, ℓ) − ip(p − 1)ME(p − 2, ℓ), (73) for p ≥ 0, 0 ≤ ℓ ≤ 2, and ME(ℓ1, ℓ2) = δℓ1,ℓ2 . This is an ms = 2 representation of order 3. Since the asymptotic form of the physical states goes as Ψ(x) ∼ exp(− 2 5 |x| 5 2 ), we can use the Gaussian weight R(x) = exp(−x2/2). The OPPQ-BM formalism ensues as before: Ψ(x) = ∞∑ j=0 cjPj(x)R(x), (74) involving real weights and corresponding orthonormal polynomials. The complex projection coefficients are given by cj(E, µ0, µ1, µ2) = ⟨Pj(x)|Ψ⟩, (75) or cj(E, µ0, µ1, µ2) = ms=2∑ ℓ=0 Λ(j) ℓ (E) µℓ, (76) where the coefficients depend on the orthonormal polynomial coefficients and the ME(p, ℓ)’s: Λ(j) ℓ (E) = j∑ η=0 Ξ(j) η ME(η, ℓ). (77) The OPPQ-BM quantization condition is I[Ψ, R] = ∫ ℜ dx Ψ∗(x)Ψ(x) R(x) = { finite ⇐⇒ E = Ephys and −→µ = −→µ phys, ∞ ⇐⇒ E ̸= Ephys or −→µ ̸= −→µ phys. (78) We adopt the vector notation for the missing moments −→µ ≡ (µ0, µ1, µ2). Substituting the OPPQ representation we obtain: I[Ψ, R] = ∞∑ j=0 c∗ j (E, −→µ )cj(E, −→µ ). (79) Once again, the focus is on : SN (E, −→µ ) = ms∑ ℓ1=0 ms∑ ℓ2=0 µ∗ ℓ1 ( N∑ j=0 Λ(j) ℓ1 (E)∗Λ(j) ℓ2 (E) ) µℓ2 , (80) involving the positive definite matrix ( PN (E) ) ℓ1,ℓ2 = N∑ j=0 Λ(j) ℓ1 (E)∗Λ(j) ℓ2 (E). (81) If we adopt a unit missing moment vector nor- malization, then the focus for OPPQ-BM quantiza- tion is on the behavior of the smallest eigenvalue, λN (E) ≡ Smallest Eigenvalue ( PN (E) ) , viewed as a function of the real and imaginary parts of the energy variable. The focus is on determining the lo- cal minima: ∂Er λN (E(min) N,phys) = ∂Ei λN (E(min) N,phys) = 0. These derivatives can be obtained algebraically. We outline the essential steps for algebraically gen- erating the partial derivatives. First of all, we have λN (E) = ⟨−→µ |PN (E)|−→µ ⟩, involving the normalized missing moment (lowest) eigenvector. This is not an analytic function in E, so we must work with the partial derivatives with respects to Er,i. Thus, we need ∂Er,iλN (E) = ⟨−→µ |∂Er,iPN (E)|−→µ ⟩. We can gen- erate the required matrix expression from Eq. (80) or ∂Er,i PN (E) = ∑N j=0 ( ∂Er,i Λ(j)(E)∗ ) Λ(j)(E) + c.c.. From Eq. (77), so long as the Ξ’s do not depend on the energy parameter, we need ∂Er,i ME(η, ℓ). How- ever, this can be obtained recursively through the expression in Eq. (72). Our immediate interest is on the behavior of λN (E) over the two dimensional complex energy plane, as the parameter ‘a’ is varied, and N = 40. In Figures 5-12 we show the convergence of two PT- breaking (complex conjugate) energies (eventually be- coming real) as we increase the a-parameter through 76 vol. 62 no. 1/2022 Algebraic Eigenenergy Bounding Method Figure 5. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −3.30 Figure 6. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −2.90 Figure 7. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −2.70 Figure 8. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −2.50 Figure 9. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −2.30 Figure 10. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −2.10 77 Carlos R. Handy Acta Polytechnica Figure 11. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = −1.70 Figure 12. Log10(λ40(E)) for PT-breaking regime, −3.30 ≤ a ≤ −1.70 for V (x) = ix3 + iax, a = 0.00 −3.30, −2, 90, −2.70, −2.50, −2.30, −2.10, −1.70 and 0. These parameter values were chosen since it is known that PT-symmetry breaking occurs at ac1 = −2.611809356 (as well as at ac2 = −5.375879629) [8, 25]. Figure 7 shows the behavior of the physical energies (i.e. defined by the local minima), for a = −2.70 ≈ ac1 , the critical value for the onset of symmetry breaking. For this same a = −2.70, Figure 4 shows the behavior of successive λN (E) surfaces for N = 20, 30, 40, 50. Within the two dimensional graphical renderings, we see similar behaviors to that in Figures 1-3. We can im- plement the same OPPQ-BM bounding analysis used previously, for hermitian systems, to bound the real and imaginary parts of the discrete state real/complex energies. The nesting of the respecive λN (Er, Ei) in Figure 4 shows the viability of the previous OPPQ-BM analysis for bounding the real and imaginary parts of the complex-plane energies. This is the focus of a future work. We note that for a = 0, the local minima in Fig- ure 12 correspond, approximately (due to the low order, N = 40) to the PT-symmetric states with ener- gies E0 = 1.15626707198811, E2 = 4.10922875280956, E3 = 7.56227385497590 and E4 = 11.31442182025857 (not shown). These were determined by EMM [8] and OPPQ-AM [10]. 11. Conclusions We have demonstrated the effectiveness of a new eigenenergy bounding procedure implementable for multidimensional discrete states regardless of the her- mitian or non-hermitian character of the associated Schrödinger operator. The discussion centered on several one dimensional systems of this type. The extension to multidimensions have been given else- where [9]. The approach advocated here is purely algebraic. Acknowledgements The author is greatful to the organizers of the 2021 Confer- ence on Analytic and Algebraic Methods In Physics XVIII, particularly Dr. Miloslav Znojil and Dr. Vit Jakubsky, and the support from the Doppler Institute of the Czech Academy of Sciences. Additional remarks from Dr. Daniel Bessis, Dr. John R. Klauder, Dr. Jerzy Cioslowski, and Dr. A. Turbiner are greatly appreciated. Technical support from Dr. Maribel Handy is also greatly appreciated. References [1] W. Ritz. Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen Physik. Journal für die reine und angewandte Mathematik 1909(135):1–61, 1909. https://doi.org/10.1515/crll.1909.135.1. [2] G. Temple. The theory of Rayleigh’s principle as applied to continuous systems. Proceedings of the Royal Society of London Series A, Containing Papers of a Mathematical and Physical Character 119(782):276–293, 1928. http://www.jstor.org/stable/95047. [3] M. G. Marmorino, A. Almayouf, T. Krause, D. Le. Optimization of the temple lower bound. Journal of Mathematical Chemistry 50(4):833–842, 2012. https://doi.org/10.1007/s10910-011-9927-z. [4] R. Martinazzo, E. Pollak. Lower bounds to eigenvalues of the schrödinger equation by solution of a 90-y challenge. Proceedings of the National Academy of Sciences 117(28):16181–16186, 2020. https://doi.org/10.1073/pnas.2007093117. [5] C. M. Bender, S. Boettcher. Real spectra in non-Hermitian Hamiltonians having PT symmetry. Physical Review Letters 80:5243–5246, 1998. https://doi.org/10.1103/PhysRevLett.80.5243. [6] P. Dorey, C. Dunning, R. Tateo. Supersymmetry and the spontaneous breakdown of PT symmetry. Journal of Physics A: Mathematical and General 34(28):L391–L400, 2001. https://doi.org/10.1088/0305-4470/34/28/102. 78 https://doi.org/10.1515/crll.1909.135.1 http://www.jstor.org/stable/95047 https://doi.org/10.1007/s10910-011-9927-z https://doi.org/10.1073/pnas.2007093117 https://doi.org/10.1103/PhysRevLett.80.5243 https://doi.org/10.1088/0305-4470/34/28/102 vol. 62 no. 1/2022 Algebraic Eigenenergy Bounding Method [7] C. M. Bender. Making sense of non-Hermitian Hamiltonians. Reports on Progress in Physics 70(6):947–1018, 2007. https://doi.org/10.1088/0034-4885/70/6/r03. [8] C. R. Handy. Generating converging bounds to the (complex) discrete states of P 2+iX3+iαX Hamiltonian. Journal of Physics A: Mathematical and General 34(24):5065–5081, 2001. https://doi.org/10.1088/0305-4470/34/24/305. [9] C. R. Handy. Exact Christoffel-Darboux expansions: A new, multidimensional, algebraic, eigenenergy bounding method. Physica Scripta 96(7):075201, 2021. https://doi.org/10.1088/1402-4896/abf67e. [10] C. R. Handy, D. Vrinceanu. Orthogonal polynomial projection quantization: a new Hill determinant method. Journal of Physics A: Mathematical and Theoretical 46(13):135202, 2013. https://doi.org/10.1088/1751-8113/46/13/135202. [11] C. R. Handy, D. Vrinceanu. Rapidly converging bound state eigenenergies for the two dimensional quantum dipole. Journal of Physics B: Atomic, Molecular and Optical Physics 46(11):115002, 2013. https://doi.org/10.1088/0953-4075/46/11/115002. [12] J. Le Guillou, J. Zinn-Justin. The hydrogen atom in strong magnetic fields: Summation of the weak field series expansion. Annals of Physics 147(1):57–84, 1983. https://doi.org/10.1016/0003-4916(83)90067-2. [13] J. A. Shohat, J. D. Tamarkin. The Problem of Moments. American Mathematical Society, Providence, RI, 1963. [14] C. R. Handy, D. Bessis. Rapidly convergent lower bounds for the Schrödinger-equation ground-state energy. Physical Review Letters 55:931–934, 1985. https://doi.org/10.1103/PhysRevLett.55.931. [15] C. R. Handy, D. Bessis, G. Sigismondi, T. D. Morley. Rapidly converging bounds for the ground-state energy of hydrogenic atoms in superstrong magnetic fields. Physical Review Letters 60:253–256, 1988. https://doi.org/10.1103/PhysRevLett.60.253. [16] Y. P. Kravchenko, M. A. Liberman, B. Johansson. Exact solution for a hydrogen atom in a magnetic field of arbitrary strength. Physical Review A 54:287–305, 1996. https://doi.org/10.1103/PhysRevA.54.287. [17] C. Schimeczek, G. Wunner. Accurate 2d finite element calculations for hydrogen in magnetic fields of arbitrary strength. Computer Physics Communications 185(2):614–621, 2014. https://doi.org/10.1016/j.cpc.2013.09.023. [18] A. Grossmann, J. Morlet. Decomposition of Hardy functions into square integrable wavelets of constant shape. SIAM Journal on Mathematical Analysis 15(4):723–736, 1984. https://doi.org/10.1137/0515056. [19] C. Handy, R. Murenzi. Moment-wavelet quantization: a first principles analysis of quantum mechanics through continuous wavelet transform theory. Physics Letters A 248(1):7–15, 1998. https://doi.org/10.1016/S0375-9601(98)00645-8. [20] C. R. Handy. Singular-perturbation–strong-coupling field theory and the moments problem. Physical Review D 24:378–383, 1981. https://doi.org/10.1103/PhysRevD.24.378. [21] S. Boyd, L. Vandenberghe. Convex optimization. Cambridge University Press, New York, 2004. [22] J. B. Lasserre. Moments, positive polynomials and their applications. Imperial College Press, London, 2010. [23] V. Chvátal. Linear programming. Freeman, New York, 1983. [24] M. Reed, B. Simon. Methods of modern mathematical physics. Vol. 1-3. Academic Press, New York, 1972. https: //doi.org/10.1016/B978-0-12-585001-8.X5001-6. [25] C. R. Handy, D. Khan, X.-Q. Wang, C. J. Tymczak. Multiscale reference function analysis of the PT symmetry breaking solutions for P 2+iX3+iαX Hamiltonian. Journal of Physics A: Mathematical and General 34(27):5593–5602, 2001. https://doi.org/10.1088/0305-4470/34/27/309. [26] E. Delabaere, D. T. Trinh. Spectral analysis of the complex cubic oscillator. Journal of Physics A: Mathematical and General 33(48):8771–8796, 2000. https://doi.org/10.1088/0305-4470/33/48/314. [27] C. R. Handy. Generating converging eigenenergy bounds for the discrete states of the ix3 non-Hermitian potential. Journal of Physics A: Mathematical and General 34(19):L271–L277, 2001. https://doi.org/10.1088/0305-4470/34/19/102. 79 https://doi.org/10.1088/0034-4885/70/6/r03 https://doi.org/10.1088/0305-4470/34/24/305 https://doi.org/10.1088/1402-4896/abf67e https://doi.org/10.1088/1751-8113/46/13/135202 https://doi.org/10.1088/0953-4075/46/11/115002 https://doi.org/10.1016/0003-4916(83)90067-2 https://doi.org/10.1103/PhysRevLett.55.931 https://doi.org/10.1103/PhysRevLett.60.253 https://doi.org/10.1103/PhysRevA.54.287 https://doi.org/10.1016/j.cpc.2013.09.023 https://doi.org/10.1137/0515056 https://doi.org/10.1016/S0375-9601(98)00645-8 https://doi.org/10.1103/PhysRevD.24.378 https://doi.org/10.1016/B978-0-12-585001-8.X5001-6 https://doi.org/10.1016/B978-0-12-585001-8.X5001-6 https://doi.org/10.1088/0305-4470/34/27/309 https://doi.org/10.1088/0305-4470/33/48/314 https://doi.org/10.1088/0305-4470/34/19/102 Acta Polytechnica 62(1):63–79, 2022 1 Introduction 1.1 Singular-perturbation strong coupling problems and moment representations 2 The moment equation representation 3 Generating eigenenergy bounds within a moments' representation 4 The eigenvalue moment method 5 The orthonormal polynomial projection quantization formalism 5.1 The OPPQ non-orthogonal basis expansion 5.2 The OPPQ quantization condition 5.2.1 Selection of the OPPQ weight 5.2.2 Use of the ground state as a weight 6 The OPPQ-approximation method (OPPQ-AM) 7 The OPPQ-bounding method 8 The OPPQ-bounding method 9 The sextic anharmonic double well potential 9.1 EMM- 9.2 EMM-2 9.3 EMM e-x44(x) 9.4 OPPQ analysis of the (ms = 0) sextic anharmonic double well oscillator 9.5 OPPQ results for the ms=0, sextic anharmonic oscillator 9.5.1 Generating bounds for E8 9.6 Comparison with EMM bounds 9.7 Using the ground state as a weight 9.8 EMM results for the probability density 10 A PT-symmetry breaking, non-hermitian, problem 11 Conclusions Acknowledgements References