The SUSY partners of the QES sextic potential revisited

In this paper, the SUSY partner Hamiltonians of the quasi-exactly solvable (QES) sextic potential Vqes(x)=νx6+2νμx4+μ2−(4N+3)νx2 , N∈Z+ , are revisited from a Lie algebraic perspective. It is demonstrated that, in the variable z = x 2, the underlying sl2(R) hidden algebra of V qes(x) is inherited by its SUSY partner potential V 1(x) only for N = 0. At fixed N > 0, the algebraic polynomial operator h(x, ∂ x ; N) that governs the N exact eigenpolynomial solutions of V 1 is derived explicitly. These odd-parity solutions appear in the form of zero modes. The potential V 1 can be represented as the sum of a polynomial and rational parts. In particular, it is shown that the polynomial component is given by V qes with a different non-integer (cohomology) parameter N1=N−32 . A confluent second-order SUSY transformation is also implemented for a modified QES sextic potential possessing the energy reflection symmetry. By taking N as a continuous real constant and using the Lagrange-mesh method, highly accurate values (∼20 s. d.) of the energy E n = E n (N) in the interval N ∈ [ − 1, 3] are calculated for the three lowest states n = 0, 1, 2 of the system. The critical value N c above which tunneling effects (instanton-like terms) can occur is obtained as well. At N = 0, the non-algebraic sector of the spectrum of V qes is described by means of compact physically relevant trial functions. These solutions allow us to determine the effects in accuracy when the first-order SUSY approach is applied on the level of approximate eigenfunctions.


Introduction
In quantum mechanics, quasi-exactly solvable (QES) systems are spectral problems, H ψ = E ψ, for which it is possible to obtain in closed (analytical) form only a finite number of exact solutions while the remaining ones remain unknown.A systematic approach to quasi-exact solvability is purely algebraic [25,24,27].In this case, the Hamiltonian can be rewritten as a constant coefficient quadratic combination in the generators of a Lie algebra.The underlying hidden algebraic structure then leads to non-trivial dynamical features of these systems.For instance, one can mention the appearance of the so-called energy-reflection (ER) symmetry [22] as well as the existence of a generating function for a set of orthogonal polynomials P n (E) in the energy variable E [5].
Interestingly, starting from a given QES model, the generation of new classes of QES problems was pointed out by means of supersymmetric quantum mechanics (SUSY) techniques [21,20,19,14].A key element is the factorization method as described first in Dirac's book [11] and further developed in Ref. [15], whereas complete reviews on supersymmetric quantum mechanics can be found in Ref. [9,4,2,12,16].In particular, in Ref. [14], the SUSY framework is one of the three methods employed to derive new QES potentials.Even more, in the last two decades, much attention has been paid to quasi-solvable models from the point of view of N -fold supersymmetry [1,3] (see also [23] and references therein).
In most cases, the key ideas about QES models are illustrated taking as a prototype the sextic potential V qes (x) = ν 2 x 6 +2 ν µ x 4 +(µ 2 − (4N +3)ν) x 2 .This one-dimensional problem possesses a hidden sl 2 Lie algebra.If the parameter N is a positive integer number, it is possible to find only (N + 1) exact analytical solutions.The exact ground state ψ qes 0 (x; N ) can always be obtained explicitly.Therefore, the first-order SUSY partner Hamiltonian H 1 with potential V 1 (x) can also be constructed for any positive integer value of N .A natural question arises, namely, Does the SUSY partner potential of V qes possess a hidden Lie algebra?.To the best of the authors' knowledge, such a relevant question has not been previously studied in the literature.In this respect, the situation is very different from previous works where the emphasis relies on finding methods for generating QES potentials.Here, we focus on the search for a hidden algebraic structure and special properties of the first-order partner Hamiltonian H 1 and its solutions.
Along these lines, we also explore the sextic QES system combining the SUSY framework and the variational method to analyze the connection, in the non-algebraic sector of the spectrum, between the bosonic and fermionic QES Hamiltonians.Since highly accurate approximate solutions are constructed, the description of the first-order SUSY mechanism at the level of approximate solutions is investigated in detail.
The goal of the present study is threefold.Firstly, for the lowest values of the parameter N = 0, 1, 2, we will derive an algebraic polynomial operator h(x, ∂ x ; N ).As a distinguished feature, the zero modes of h(x, ∂ x ; N ) give the N exact polynomial eigenfunctions of the SUSY partner potential V 1 (x) which, in general, is not a polynomial but a rational function in the variable x 2 .Towards a Lie-algebraic characterization, the analytical properties of V 1 are clearly indicated.For N > 2, the generic expression of h(x, ∂ x ; N ) is presented as well.To overcome the 1-SUSY requirement of the explicit knowledge of the exact (analytical) ground state function ψ qes 0 (x; N ) a confluent secondorder SUSY transformation is applied on a QES sextic potential to generate isospectral QES models.
Secondly, for the three lowest states n = 0, 1, 2 of the QES potential V qes (x), highly accurate values of the energy E n = E n (N ) as a function of N are displayed within the interval N ∈ [−1, 3].Thus, we investigate the system with N being considered as a continuous real parameter.The numerical results are presented with ∼ 20 significant digits; the corresponding calculations are performed using the Lagrange-mesh method.In particular, the critical value N c above which tunneling effects can occur is computed.
Finally, for the case N = 0 where solely the exact ground state of V qes (x) is known, we generate approximate solutions for V 1 (x).They are constructed by means of SUSY techniques allied with the variational method.As a first step, compact variational trial functions ψ trial for the excited states of V qes are designed.They encode relevant physical properties of the system.Afterward, a SUSY transformation acting on these ψ trial will produce approximate eigenfunctions for V 1 .This simple idea allows us to determine the effects in accuracy when SUSY is implemented on the level of approximate eigenfunctions, another interesting practical aspect absent in the literature.

Generalities
We consider the following one-dimensional spectral problem in non-relativistic quantum mechanics: defined on the real line, x ∈ (−∞, ∞).The corresponding Hamiltonian operator for our specific problem is of the form: where m denotes the mass of the particle, and is the quasi-exactly solvable sextic potential, here ν, µ, N are real parameters.At ν = 0, the potential (3) reduces to the exactly-solvable simple harmonic oscillator V es (x) = µ 2 x 2 .Hereafter, we will adopt atomic units ℏ = 1, m = 1.Some remarks are in order: • For ν ̸ = 0 and arbitrary N and µ, N not necessarily an integer number, the potential V qes (3) admits infinitely many bound states, see Fig. 1.No scattering states occur.• At fixed µ and ν > 0, there exists a special value N = µ 2 −3 ν 4 ν above which the potential (3) develops two symmetric degenerate minima.Thus, tunneling effects (instanton-like terms) can take place.For large N → ∞, these two minima are located at , respectively, and • The Hamiltonian ( 2) is invariant under the parity symmetry x → −x, which gives rise to alternating symmetric and anti-symmetric bound states.
• At µ = 0, N and ν ̸ = 0 arbitrary, the energy reflection (ER) symmetry in (2) emerges [22].For instance, the wave functions of the ER symmetric levels are connected to each other by analytic continuation x → i x, ψ E → ψ −E .
• Formally, the most general eigenfunction of (2) reads here c 1 , c 2 are constants of integration, and w(z) = w(q, α, γ, δ, ϵ, z) satisfies the bi-confluent Heun differential equation z w ′′ + w ′ (z 2 ϵ + δ z + γ) + w (α z − q) = 0. We will focus on the case where ψ(x) is normalizable, and one of its factors is a polynomial function in the x−variable.Other solutions are also interesting; see for example [18] where a sextic potential with a centrifugal barrier is considered.
• The Hamiltonian (2) possesses a hidden sl 2 (R) Lie algebra [25].In the particular case when N is a positive integer number, N ∈ Z + , this algebra admits a finitedimensional irreducible representation and, thus, one can find algebraically only (N +1) exact even-parity eigenfunctions and their corresponding energies explicitly.
• As stated in [17], for a positive integer number N , the system (3) is intimately related to the exact cancellation of real and complex nonperturbative saddles to all orders in the semiclassical expansion.
3. The hidden sl 2 Lie algebra of V qes (x) We begin with the well-established result that the Hamiltonian H in (2) can be transformed into a suitable sl 2 Lie algebraic operator [25], [26].To this end, the gauge factor: is introduced.In the Z 2 -invariant variable z = x 2 , one can construct the gauge-rotated Hamiltonian which can be rewritten as a constant coefficient quadratic combination in terms of the sl 2 generators [25] Explicitly, Moreover, if the parameter N takes positive integer values, then the spectral problem possesses (N + 1)-exact polynomial eigenfunctions P j with j = 0, 1, 2, . . ., N .In this case, the (N + 1)-exact solutions of the original Hamiltonian H in (1) take the form All of them correspond to even-parity states, including the ground state, possessing an even number of nodes.It is worth mentioning that odd-parity exact solutions of H (1) appear for positive half-integer values of N = 1 2 , 3 2 . ... They correspond to excited states and can be associated with the existence of the same hidden sl 2 algebra.However, the basic object we will use later in first-order supersymmetric quantum mechanics is precisely the ground state solution.Therefore, the present study does not consider such odd-parity exact states of H (1) except in the part concerning the confluent second-order SUSY transformation.

Exact even-parity polynomial solutions: cases
This section reviews explicit results for the lowest positive integer values of N .They are well-known in the literature and will be used later in the present consideration.For convenience, we take ν = 1 and µ = 1.The QES sextic potential reads see Fig. 1.The corresponding exact solutions of the Schrödinger equation are the following: • At N = 0, one exact solution occurs only.It is given by the nodeless function: with energy It corresponds to the exact ground state of the system.Accordingly, the polynomial factor P (N =0) 0 ≡ 1 is an eigenfunction of the sl 2 Lie-algebraic operator h (8).
• At N = 1, solely two exact solutions appear: with energy and here They correspond to the ground and second excited state, respectively.The polynomial factors P = 2 x 2 − √ 3 + 1 are exact eigenfunctions of the Lie-algebraic operator h (8).
• At N = 2, following the general theory, three exact analytical solutions occur only: with energy Next, where and with eigenvalue They correspond to the ground, second, and fourth excited states.Thus, the corresponding factors P 1 are also eigenpolynomials of the Lie-algebraic operator h (8).

QES sextic potential with energy reflection symmetry: examples with N = 0, 1
As previously mentioned in the introduction, at µ = 0, the Hamiltonian (2) with potential V qes (x) as in (3) possesses the energy reflection (ER) symmetry E ↔ −E [22].For the potentials with ER symmetry, we study a more general QES potential, the two-parametric sextic potential V qes ER (x): where κ = 0, 1 indicates the parity of the algebraic solutions.The case κ = 0 corresponds with the potential V qes (x) in (3) at µ = 0. Let us revisit in detail concrete solutions using the most straightforward examples, namely N = 0 and N = 1.They will be exploited later in the confluent SUSY algorithm, a degenerate case of the second-order supersymmetric quantum mechanics: • At N = 0, the potential, the single eigenfunction, and the corresponding eigenvalue are where ν > 0 and is the energy of the ground (κ = 0) and first-excited (κ = 1) state ψ (N =0) ER , respectively.
• If N = 1, there are two exact analytical solutions for each value of κ, the potential, the analytic eigenfunctions, and its energies are: When κ = 0, the previous expressions correspond to the ground (−) and second excited (+) states whereas at κ = 1 they describe the first (−) and third excited (+) state, respectively.

4.
Exact-and WKB-numerical solutions, case N = 0 For any value of N , by means of the user-friendly LagrangeMesh Mathematica Package [10] (LMMP), one can easily determine the exact-numerical eigenfunctions and eigenvalues.It is worth mentioning that the underlying Lagrange Mesh Method is an approximate variational method simplified by a Gauss quadrature associated with a certain mesh; for further details, see Ref. [10].The high efficiency and direct control of the involved accuracy in arithmetic manipulations (and final results) allow us to obtain highly accurate energies using a personal laptop.For instance, an excessive value with 30 correct significant digits can be achieved in short CPU times.A similar situation occurs for the corresponding wave functions.Therefore, in practice, we will refer to these numerical results as exact solutions, dropping the term numerical.
Just for N = 0, we present in Table 1 the energies E n and the expectation value ⟨ x 2 ⟩ for the lowest ten bound states n = 0, 1, 2, . . ., 9, explicitly.The calculations used a Hermite mesh with 800 mesh points for the energy.A comparison with WKB results is indicated as well.
Table 1.Case N = 0: exact numerical solutions E exact for the sextic potential 1 2 V qes (x, N = 0) defined by (11).Both the energy and the corresponding expectation value ⟨ x 2 ⟩ were obtained using the LagrangeMesh Mathematica Package [10] In Table 2, the energy E = E(N ) as a function of the cohomology parameter N is depicted for the three lowest states n = 0, 1, 2. If N > − 1 2 , the potential V qes (x, N ) (11) has a maximum located at x = 0, and it corresponds to a zero energy value.It implies that for negative energies E < 0, the dynamics takes place in a classically forbidden region.At fixed n, the energy E n (N ) is a smooth decreasing function of the parameter N .In particular, there exists a critical value N = N c such that the ground state energy E 0 (N c ) of the system vanishes.Accordingly, for N > N c , the appearance of tunneling effects (exponentially small instanton-like terms) are expected.Interestingly, if N is a positive (half)integer number, the corresponding instanton-like terms are absent in the exact analytical solutions!This puzzle was resolved in [17] by including complex saddles.It is an interesting open question whether this is valid for non-positive (half)integer values of N .
Table 2. Sextic potential 1  2 V qes (x, N ) defined by (11): the energy E n = E n (N ) in atomic units, as a function of N , for the three lowest states n = 0, 1, 2. As for the critical value (see text) we found N c ≈ 0.73295312615213043.The results were obtained using the LagrangeMesh Mathematica Package [10].

First-order supersymmetry
We start from a certain Schrodinger operator , where some of its eigenfunctions ψ n (x) and eigenenergies E n are known explicitly.Then we suppose the existence of a first-order differential operator where In the last expression, u(x) is called seed solution, and it solves the equation H 0 u = ϵ u, where ϵ is a constant known as factorization energy.We use the notation u ′ (x) = du(x)/dx.The seed solution u is not necessarily a physical solution ψ of H 0 .In the literature α(x) ≡ u ′ (x)/u(x) = (ln u(x)) ′ is called the superpotential.From the intertwining relation (29) it follows the expression It is said that V 1 is the SUSY partner of V 0 .Moreover, from the equation above, the importance of avoiding zeroes of u is clear, so the SUSY partner potential V 1 is regular in the same domain of V 0 .This imposes the condition ϵ ≤ E 0 .
If we introduce the adjoint operator of A + 1 , we can show by direct substitution that the operators A 1 and A + 1 factorize the Hamiltonians as Also, from the intertwining relation and the factorization of H 0 , it can be shown that if ψ k is an eigenfunction of H 0 with eigenvalue E k then is an eigenfunction of H 1 with the same eigenvalue E k .The factor 1/ √ E k − ϵ comes from the normalization condition.
Remark The subindex k in ψ k indicates the energy level, i.e., ψ k is the k-th excited state of H 0 .In contrast, in the transformed wavefunction ϕ k of H 1 , it is not necessarily the k-th excited state, as we will show below.From the factorization of H 1 , we can see that the function annihilated by A 1 is also a solution of the eigenvalue equation of H 1 with eigenvalue ϵ.This function is known as the missing state and its expression is If ϕ ϵ fulfills the boundary conditions, ϵ belongs to the spectrum of H 1 , otherwise it does not.In general, if ϵ 0 = E 0 and u(x) = ψ 0 then Sp(H 1 ) differs from Sp(H 0 ) only on the ground state energy.Below, we present the supersymmetric partner potential V 1 (x, N ) corresponding to V 0 = 1 2 V qes (x, N ) for the lowest cases N = 0, 1, 2, 3.

1-SUSY partner potential
We take V 0 = 1 2 V qes (x, N = 0).In this case, we solely know the exact ground state solution ψ with energy E 0 = 1/2.No exact analytical solutions for the excited states of H 0 are known.Then, we can only choose u = ψ (N =0) 0 . Using (31), we immediately obtain the SUSY partner potential which is in complete agreement with the general Eq.(31) in [14].
Up to an additive constant, V 1 (35) coincides with the sextic potential 1 2 V qes (x, N = −3/2), see (11).Thus, just like V 0 , it is an analytic function with no poles on the real (and complex) domain.Remarkably, the SUSY partner Hamiltonian H 1 with potential (35) still possesses a hidden sl 2 Lie algebraic structure but with a negative value of parameter N .The trivial SUSY solution of For fixed µ = ν = 1 and N = 0, the confining potential V 0 = V qes (x)/2 defined by (3) (blue curve) and its SUSY partner V 1 (x) as in (35).

1-SUSY partner potential
Here V 0 = 1 2 V qes (x, N = 1), for which we know two exact analytical solutions, namely ψ where P (N =1) 0 = 2 x 2 + √ 3 + 1 vanishes only in the complex domain.The SUSY transformation retrieves a rational extension of V qes .For the above potential V 1 , we know only a single exact solution (the first excited state) with energy E It is very common that the SUSY transformation changes the parameters of a potential and adds a finite term.In this case, such a finite part is a rational function in the variable x 2 , with complex poles given by the complex zeros of P (N =1) 0 .Now, let us introduce the gauge factor thus, ϕ 2 = x Γ 1 and the polynomial factor of ϕ 2 is simply x.Clearly, in the x-variable the gauge rotated Hamiltonian Γ −1 1 H 1 Γ 1 will contain non-polynomial coefficients.To obtain an algebraic differential operator with polynomial coefficients we construct the gauge-rotated Hamiltonian = − (2x where λ is a real constant.That way, the spectral problem H 1 ϕ = E ϕ is equivalent to the zero-mode equation h (susy) p(x) = 0 being p = p(x) a function to be determined.In this case, ϕ 2 (x) = p(x) Γ 1 is the solution of the Hamiltonian H 1 .

1-SUSY partner potential
Here V 0 = 1 2 V qes (x, N = 2), for which we know three exact analytical solutions, ψ (N =2) 0 = 2x 4 + 6x 2 + 3.For the above potential V 1 we know two exact solutions: a first excited state with energy E 2, and a third excited state with energy E (N =2) 4 = 9 2 + 2 √ 2, respectively.Now, let us introduce the gauge factor Again, in the x-variable the gauge rotated Hamiltonian Γ −1 2 H 1 Γ 2 will contain non-polynomial coefficients.In order to obtain an algebraic differential operator with polynomial coefficients we build the gauge-rotated Hamiltonian where λ is a real constant.That way, the spectral problem H 1 ϕ = E ϕ is equivalent to the zero-mode equation h (susy) p = 0.By construction, at λ = 9 2 ± 2 √ 2, the kernel of h (susy) | N =2 admits two fifth order polynomials solutions p 5,± , i.e., they obey h (susy) | N =2 p 5,± = 0. Nevertheless, in the x-variable, this operator can not be rewritten as a constant coefficient quadratic combination in the sl 2 generators (7).

1-SUSY partner potentials V 1 with arbitrary integer N > 0
Here V 0 = 1 2 V qes (x, N ), for which we know (N + 1) exact solutions, ψ 0 , ψ 2 . . ., ψ N .Substituting u = ψ 0 in (31) we arrive to the following expression here and [N/2], respectively ‡.For the Hamiltonian H 1 with the above potential V 1 (45), the N known exact solutions are given by (33).In (45), the part containing no complex poles is given by the sextic potential 1 2 V qes (x, N ) with the replacement N → N − 3 2 .Remark: At large distances |x| → ∞, the last two rational terms in (45) vanish.The dominant term behaves as 1 x 2 .Thus, in this limit the asymptotic behaviour x 2 ] of the solutions for V 1 is the same as that occurring for V qes .Using the gauge factor we construct the gauge-rotated Hamiltonian where λ is a real constant and Y N +1 , Z N +1 are polynomial functions in x 2 of the same order (N + 1).Hence, the spectral problem By construction, the kernel of h (susy) | N admits N polynomial solutions p 4N −3 of order (4N − 3), i.e., they obey h (susy) | N p 4N −3 = 0. Nevertheless, in the x-variable, this operator can not be rewritten as a constant coefficient quadratic combination in the sl 2 generators (7).
Moreover, we can recover an intertwining relation for the Hamiltonian h given in (6).From the intertwining relation (29) we can define the SUSY partner of h as h 1 = Γ −1 N hΓ N and the intertwining operator as Also, by construction, if P (z) solves hP (z) = EP (z) then P(z) = A † 1 P (z) solves h 1 P(z) = E P(z).

Second-order supersymmetry: confluent case
Let us present this second-order SUSY transformation as an iteration of two first-order transformations; see more details in [13,8,7].Based on the results of the previous subsection 5.1, from H 1 we will construct a Hamiltonian H 2 = − 1 2 d 2 dx 2 + V 2 (x) using a second intertwining operator A + 2 , namely where the seed function v(x) solves the spectral problem H 1 v = ϵ v.Note that we have used the same factorization energy ϵ of the previous SUSY step; this is the characteristic feature of the confluent SUSY transformation.The SUSY partner potential becomes A simple choice of v comes from (33), v = 1/u; however, this selection leads to V 2 = V 0 .
From the reduction-of-order formula, the general choice is given by where ω 0 is a constant.To simplify notation, we can define Accordingly, the second-order confluent SUSY partner potential of V 0 is Note that there is an intertwining relation between H 0 and H 2 : Moreover, the operator B + and its adjoint B = (B + ) + factorize as follows: The missing state of H 2 becomes and the eigenfunctions of H 0 are mapped to those of H 2 as Below, we display the confluent SUSY partners of the potentials with energy reflection symmetry presented in subsection 3.2.To simplify notation, we set ν = 1.

Even parity with N = 0
To obtain the confluent SUSY partner of where Γ(z) is the Gamma function, and E n (z) is the exponential integral function ∞ 1 e −zt /t n dt.We fixed the lower limit of the integral x 0 = 0 to use the parity properties of u 2 .It is necessary to avoid zeros of ω(x), for this reason V with ER symmetry (blue curve) and its confluent SUSY partner (orange curve).Right: Ground state of this sextic potential (blue curve) and the corresponding to the confluent SUSY partner (orange curve).The parameters are κ = 0, µ = 0, ν = 1, ω 0 = 2 1 4 Γ (5/4) + 0.01.
Remark: For the confluent SUSY partner potential V 2 (x; N = 0) the original ER symmetry remains.However, the wave functions of the ER symmetric levels are not connected by the original analytical continuation x → i x.

An example with N = 1
Since, at N = 1, there are two known exact eigenfunctions of H 0 , there are two different choices for the seed functions u(x).For even parity κ = 0, they are the ground and the V ) with ER symmetry (blue curve) and its confluent SUSY partner (orange curve).Right: First excited state of V 0 (blue curve) and the corresponding of the confluent SUSY partner (orange curve).The parameters are κ = 1, µ = 0, ν = 1, ω 0 = 0.01 + 2 −5/4 Γ(3/4).second excited states.For odd parity κ = 1, they correspond to the first and the third excited states.The confluent algorithm can be applied as in the previous cases.Firstly, we pick the seed function (u = ψ (N =1) ER,± , see ( 28)); so we can construct the function ω(x) defined as in (52).Secondly, it is important to find the domain of ω 0 such that ω(x) is nodeless.Thirdly, it is straightforward to obtain the confluent SUSY partner V 2 (see (53)).Eventually, there will be two eigenfunctions of H 2 associated with the same energies of the initial system H 0 , one of them ϕ ϵ is constructed as in (55) and the second one ϕ n is calculated using (56).As a result, we can build four different Hamiltonians because, for each parity, we have two options for the seed function.
7. SUSY in the non-algebraic sector of V (qes) : case N = 0 Let us consider again the case N = 0.For the QES sextic potential only the ground state function ψ 0 (12) is know analytically.The corresponding 1-SUSY partner potential V 1 (x) can be found immediately.Even though no more exact analytical solutions exist other than ψ 0 , it is important to emphasize that both potentials, V 0 and V 1 , as well as their corresponding (unknown) eigenfunctions and eigenvalues, are still connected by SUSY means.The set of excited states of V 0 can be calculated using numerical and approximate methods.In particular, for one-dimensional systems, a SUSY scheme based on a hierarchy of Hamiltonians has been introduced in previous works to compute the excited states (see [9] and references therein).However, the accuracy of the so obtained solutions is rather limited.
Therefore, it would be worth analyzing the effect of a SUSY transformation acting on the approximate solutions of V 0 .This would generate approximate solutions for V 1 .
The accuracy of approximate solutions can be easily estimated using a direct numerical method.Hence, the concrete question we aim to answer can be formulated as follows: How does the accuracy of approximate solutions change when a SUSY transformation is applied?.
We adopt the variational approach to compute the first excited states of V 0 (x) as a first step.The corresponding trial functions ψ trial (x) are designed on physical grounds and a criterion of simplicity.The accuracy of the obtained solutions, energies, and wave functions is estimated.Afterward, using the operator A + 1 , we calculate approximate solutions ϕ = A + 1 ψ trial for V 1 and determine how the accuracy is modified.
7.1.V 0 = 1 2 V qes (x, N = 0): approximate variational solutions 7.1.1.First excited state For the first excited state of V 0 , we employ the trial function where the c i , (i = 0, 1, 2, . . ., k), are (k + 1)-variational parameters to be determined by the minimization procedure of the energy functional.As a result of calculations, we can always put c 0 = 1.The function (61) possesses the following properties: • the orthogonality condition with the exact ground state function is satisfied identically, ⟨ ψ • it possesses a definite odd parity ψ trial (−x) = − ψ trial (x).
• the node is correctly located at x = 0.

The corresponding energy functional
can be evaluated analytically in terms of Bessel functions.However, we do not present the corresponding (lengthy) expression explicitly.
For the lowest values of k = 1, 2, . . .6, the minimization of (62) gives the results of the energy E (1st) var displayed in Table 3.Using the highly accurate LagrangeMesh Mathematica Package [10] we obtain the corresponding exact numerical result E   7.1.2.Second excited state Similarly, for the second excited state of V 0 = 1 2 V qes (x, N = 0) we use the following trial function where the b i , (i = 0, 1, 2, . . ., k), are (k + 1)-variational parameters to be determined by the minimization procedure of the energy functional and the orthogonality condition.
For the values of k = 4, 5, 6, 7, the energy E    trial (64) is of order 10 −4 or less.In particular, the two nodes of (64) are located at x ≈ ±0.5017 in agreement with the exact numerical result.where the values of parameters c i correspond to the optimal results obtained for ψ (1st) trial previously.The energy obtained using ϕ 1,trial (x) as a variational function with no-free parameters is presented in Table 5.
Next, we compute the approximate first excited state of V 1 acting the operator A + 1  (30) onto the second excited state ψ (2nd) trial defined in (64).Explicitly, where the optimal values of parameters b i are taken from (65).In this case the relative error e r ≡ ϵvar−ϵexact ϵexact is of order ≈ 10 −7 with respect to the exact value ϵ 2 = E (2nd) exact .By construction, the orthogonality condition ⟨ ϕ 2,trial |ϕ 0,trial ⟩ = 0 is fulfilled exactly.

Conclusions
In summary, for the QES sextic potential with integer N > 0 an algebraic polynomial operator h(x, ∂ x ; N ) that governs the N exact polynomial solutions of its 1-SUSY partner V 1 (x) is constructed.These odd-parity solutions are polynomials in the variable x of order (4N − 3), occurring in the form of zero modes.Nor in the x-variable neither in x 2 , the operator h(x, ∂ x ; N ) is sl 2 Lie-algebraic.In the case N = 0, the potential V 1 possesses a sℓ 2 hidden Lie algebra, but no exact solutions occur.
At fixed N > 0, the potential V 1 (x) splits into two additive parts.The first one is polynomial, and it is given by V qes (x) with a different quantized parameter N → N − 3 2 ,

Figure 1 .
Figure 1.At fixed µ = ν = 1, we display the confining potential V qes (x) in (3) for different values of the parameter N .For N > − 1 2 it develops two symmetric degenerate minima.
49) ‡ [a] denotes the integer part of a

7. 2 .
Ground state and first excited state of V 1 (x): approximate SUSY solutions Now, for the 1-SUSY partner potential V 1 (35) we calculate the approximate ground state solution ϕ 1,trial by simply acting the operator A + 1 (30) on the first excited state ψ (1st) trial defined in (61).No further variational minimization is involved.Explicitly, ϕ 0,trial = A + 1 ψ

Figure 11 .
Figure 11.Ground state function ϕ 0 of V 1 (x, N = 0).The exact numerical wavefunction (left), and the local relative error δϕ ≡ ϕ trial −ϕexact ϕexact . Results are displayed in atomic units.For comparison, the WKB calculations E W KB are presented as well.Here ∆E ≡ E W KB −Eexact

Table 5 .
Ground state energy of V 1 (x; N = 0): the energy ϵ 1,var is obtained calculating the expectation value of H 1 on the function ϕ 1,trial = A + 1 ψ