Third-order ladder operators, generalized Okamoto and exceptional orthogonal polynomials
TThird-order ladder operators, generalized Okamotoand exceptional orthogonal polynomials
V. Hussin ∗ , , I. Marquette † , and K. Zelaya ‡ D´epartement de Math´ematiques et de Statistique, Universit´e de Montr´eal, Montr´eal H3C 3J7, QC,Canada School of Mathematics and Physics, The University of Queensland, Brisbane, QLD 4072, Australia
Abstract
We extend and generalize the construction of Sturm-Liouville problems for afamily of Hamiltonians constrained to fulfill a third-order shape-invariance condi-tion and focusing on the “ − x/
3” hierarchy of solutions to the fourth Painlev´etranscendent. Such a construction has been previously addressed in the literaturefor some particular cases but we realize it here in the most general case. The corre-sponding potential in the Hamiltonian operator is a rationally extended oscillatordefined in terms of the conventional Okamoto polynomials, from which we identifythree different zero-modes constructed in terms of the generalized Okamoto poly-nomials. The third-order ladder operators of the system reveal that the completeset of eigenfunctions is decomposed as a union of three disjoint sequences of solu-tions, generated from a set of three-term recurrence relations. We also identify alink between the eigenfunctions of the Hamiltonian operator and a special family ofexceptional Hermite polynomial.
Nonlinear equations have played a fundamental role in understanding the dynamics ofsome physical models, even in cases where the governing physical laws are defined interms of linear equations [1]. Particularly, in quantum mechanics, the nonlinear Riccatiequation [1, 2] has found several interesting applications in the study of exactly-solvablemodels. Such a task is achieved through the factorization method [3–9], a technique in-trinsically related to the Darboux transformation [10, 11]. The latter is also known as ∗ [email protected] † [email protected] ‡ Corresponding author: [email protected] a r X i v : . [ m a t h - ph ] F e b upersymmetric quantum mechanics (SUSY QM) because of the mathematical equiva-lence with the supersymmetric construction of Witten [12–14] in the potential theorymodel. For stationary quantum systems, that is, time-independent Hamiltonians, the cor-responding dynamical law reduces to a Sturm-Liouville eigenvalue problem in the form ofsecond-order differential equation in the spatial coordinates, the solutions of which can beaddressed in terms of either hypergeometric or confluent hypergeometric functions [15,16]for some particular interactions such as the harmonic oscillator, hydrogen atom, interac-tion between diatomic (Morse potential) and polyatomic molecules (Rosen-Morse poten-tial), transparent potential interactions (P¨oschl-teller potential), among others. In thisregard, the Darboux transformation has led to an outstanding progress in the study ofexactly solvable quantum models, where the models previously mentioned can be usedas a departing point to generate new families of potentials with spectrum on demand.Moreover, the latter does not restrict to Hermitian Hamiltonians, and non-Hermitianmodels can be constructed after imposing the non-self-adjoint condition on the result-ing Hamiltonians [17, 18]. The realitiy of the spectrum is preserved [19, 20] in systemswith either broken and unbroken P T -symmetry, extending the conventional systems with
P T -symmetry and real spectrum [21–24].On the other hand, the Painlev´e transcendents form a family of ordinary nonlinearequations that have become a topic under intensive research within both the mathemati-cian and physicists communities. Their solutions cover a wide range of mathematicalmodels with well-defined monodromy and also emerge naturally in several physical prob-lems in the classical and quantum regimes. The Painlev´e transcendents are characterizedby a family of six nonlinear equations P I - P VI defined in terms of complex parameters,whose solutions, in general, cannot be expressed in terms of classical functions [25, 26].Nevertheless, for some specific values of the parameters, a seed function can be usedto generate a complete hierarchy of solutions through the B¨acklund transformation [27],which can be thought as a nonlinear counterpart of the recurrence relations. In par-ticular, the fourth Painlev´e equation can be taken into a Riccati equation with the ap-propriate choice of the parameters [28], where it solves a “simpler” nonlinear equationinstead. Moreover, the fourth Painlev´e equation has also brought new results in thetrend of orthogonal polynomials, where new families were discovered through the hier-archies of rational solutions in terms of the generalized Okamoto, generalized Hermiteand Yablonskii-Vorob’ev polynomials [29]. The Painlev´e transcendents have also foundinteresting applications in the study of physical models in nonlinear optics [30], quantumgravity [31], and SUSYQM [32–36].In particular, the fourth Painlev´e equation arises quite naturally in third-order shape-invariant SUSYQM [37, 38], where the parameters of the Painlev´e transcendent define theeigenvalues of the new Hamiltonians and the respective intertwining operators serve atthe same time as ladder operators. Strinkingly, the so-constructed intertwining operatorsare not in general factorizable in terms of first-order operators. Thus, the results ob-tained in this way generalize those of [39]. It is worth to notice that higher-order ladderoperators have been also studied, in a different way, in the context of supersymmetric(SUSY) partners for the stationary oscillator in both the Hermitian [40, 41] and non-2ermitian regimes [42]. The Hermitian construction is particularly interesting from themathematical point of view, for it has been exploited to construct and study new familiesof exceptional orthogonal polynomials [43, 44], that is, a set of orthogonal polynomialswith some degrees absent from the polynomial sequence. Such exceptional polynomialsfind applications in the construction and study of minimal surfaces [45].In this manuscript, we develop a thorough study of quantum models that satisfy athird-order shape-invariant condition so that the corresponding potential becomes a ra-tional extension of the harmonic oscillator written in terms of the Okamoto polynomials.The latter is achieved through the use of the fourth Painlev´e equation, from which the so-lutions of the “ − x/
3” hierarchy of rational solutions is exploited. Although some workshave been published in the past for a wide range of solutions to the fourth Painkev´eequation, the rational case was studied only for some specific values of the paramaters α and β that define the fourth Painlev´e transcendent [36, 46, 47]. Thus, we provide themost general construction in terms of the “ − x/
3” hierarchy. Our construction leads tothree sequences of solutions, where the three eigenfunctions annihilated by the annihila-tion operator (zero-modes) are written in terms of the generalized Okamoto polynomials.The rest of the eigenfunctions (higher-modes) are computed after identifying a three-termrecurrence relation for each sequence. Such higher-modes are determined iteratively afterfixing the corresponding zero-mode as the initial condition. Interestingly, the latter leadsto a relationship between the higher-modes and the family exceptional Hermite polynomi-als H λ ,n defined by the particular double partition λ ≡ λ k = (1 , , , , · · · , k, k ). Thisprovides us with a set of three different three-term recurrence relations that generatethe set of exceptional Hermite polynomials { H λ k ,n } ∞ n =0 . In previous works [48–50], sev-eral constructions for the exceptional Hermite polynomials have been identified throughhigher-order recurrence relations (of order higher than three), even for the double par-tition λ k . Thus, our approach brings a new simple construction for such exceptionalpolynomials in the form of three-term recurrence relations, a property akin to classicalorthogonal polynomials.The manuscript is structured as follows. In Sec. 2, we briefly discuss the construc-tion of third-order shape invariant Hamiltonians and their relation to the fourth Painlev´eequation, with a particular emphasis to the “ − x/
3” hierarchy of rational solutions. Forsuch a case, the eigenvalues and zero-modes are determined in Sec. 3 in general. In Sec. 4,we introduce an algorithm to construct the higher-modes from the zero-modes through aset of three-term recurrence relations. In this form, the complete spectral information canbe determined. Interestingly, in Sec. 5, we establish a relation between the third-ordershape-invariant potential previously constructed and a higher-order Darboux-Crum trans-formation, which leads to a Wronksian representation of both the potential and eigenfunc-tions. The latter allows defining a three-term recurrence relation for a particular familyof exceptional orthogonal Hermite polynomials. In App. A, we summarize some usefulB¨acklund transformations for the fourth Painlev´e transcendent that lead to the identitiesused throughout the manuscript. In App. B, with aid of the Wronskian representation,we compute the explicit action of the creation operator on the eigenfunctions. Lastly, inApp. C, we summarize some basic properties of the Darboux-Crum (higher-order SUSY)3ransformations.
For completeness, in this section, we briefly discuss the construction of families of sta-tionary Hamiltonians satisfying the third-order shape-invariant condition. Those modelswere initially identified in [51], recently studied in [36, 46, 52] for the stationary case, andin [47] for the generalized time-dependent regime. In such cases, the construction of thecorresponding potentials has been handled with generality. Nevertheless, the completeset of eigenfunctions has not been identified explicitly, and only the zero-modes havebeen presented in some particular cases. Thus, we focus on the general construction ofrationally extended potentials determined in terms of the Okamoto polynomials.Let us consider the eigenvalue equation related to an unknown Hamiltonian H of theform Hφ n = E n φ n , H := − d dx + V ( x ) , (1)with φ n ( x ) the eigenfunctions and E n the respective eigenvalues. The real-valued potential V ( x ) is unknown and determined from a shape-invariant condition HA † = A † ( H + 2 λ ) , HA = A ( H − λ ) , λ > . (2)In the latter, A and A † are known as the intertwining operators . In this particular case,we are interested in the third-order shape-invariant condition [37] so that the intertwiningoperators, in coordinate representation, are represented as A := d dx + A ( x ) d dx + A ( x ) ddx + A ( x ) , (3) A † := d dx + A ( x ) d dx + [2 A (cid:48) ( x ) − A ( x )] ddx + [ A ( x ) − A (cid:48) ( x ) + A (cid:48)(cid:48) ( x )] , (4)with A i ( x ) real-valued functions, for i = 0 , ,
2, to be determined from Eq. (2), along with f (cid:48) ( x ) ≡ df /dx and f (cid:48)(cid:48) ( x ) ≡ d f /dx .Also, it is worth to remark that Eq. (2) also implies that A and A † are annihilationand creation operators, respectively, where the parameter λ dictates the growth rate ofthe energy eigenvalues E n , which increase or decrease by 2 λ units. Such a parametercan also be absorbed after a suitable reparametrization of the spatial coordinate. Thus,without loss of generality, we consider λ = 1.Although the intertwining operators A and A † given in Eqs. (3)-(4) have the mostgeneral form, we consider a convenient factorization so that the Painlev´e transcendentemerges naturally. The latter was first introduced by Cannata et al. [52], and it is usuallyknown as an irreducible factorization of the form A = M † Q , A † = Q † M , (5)4here the sets { Q, Q † } and { M, M † } contain first-order and second-order operators, re-spectively, given by M † := ddx − G ( x ) ddx + B ( x ) , M = ddx + 2 G ( x ) ddx + B ( x ) + 2 G (cid:48) ( x ) , (6) Q † := ddx + W ( x ) , Q = − ddx + W ( x ) , (7)with W ( x ), B ( x ), and G ( x ) real-valued functions to be determined. In analogy to therelations given in Eq. (2), the new operators in Eqs. (6)-(7) define an alternative set ofintertwining relations that are not of the shape-invariant type. To see the latter, let usconsider an auxiliary Hamiltonian of the form H := − d dx + V ( x ) , H ϕ (1) n = E (1) n ϕ (1) n , (8)with V ( x ) the corresponding potential, together with the eigenfunctions ϕ (1) n and eigen-values E (1) n . We thus impose the following intertwining relations: HQ † = Q † ( H + 2) , HQ = Q ( H − , (9) HM † = M † H , HM = M H , (10)such that their combined action bring us back to the original shape-invariant condition (2),with λ = 1.Now, substituting (6)-(7) into (9)-(10), we obtain a set of differential equations involv-ing the potentials V ( x ) and V ( x ) in terms of the unknown functions B , G , and W . Thestraightforward calculations show that V ( x ) = V ( x ) + 2 W (cid:48) + 2 , V ( x ) = − W (cid:48) + W − , (11) V ( x ) = V ( x ) − G (cid:48) , V ( x ) = 2 G + G (cid:48) − B + γ , (12) B = G − G (cid:48) − G (cid:48)(cid:48) G + ( G (cid:48) ) G + d G , (13)where d, γ ∈ R are integration constants. After combining (11)-(12) we obtain the twoadditional relations W = − (2 G + x ) , B = W (cid:48) − W + G (cid:48) + 2 G + γ + 2 , (14)which combined with (13) leads to G (cid:48)(cid:48) = ( G (cid:48) ) G + 6 G + 8 xG + 2[ x − ( γ + 1)] G + d G . (15)In this form, the fourth Painlev´e [27, 29, 53] transcendent w (cid:48)(cid:48) = ( w (cid:48) ) w + 32 w + 4 xw + 2( x − α ) w + βw , (16)5s recovered from (15) through the reparametrization G ( x ) = w ( x )2 , α = γ + 1 , β = 2 d . (17)Although we are only interested in the construction of V ( x ), we can also provide anexplicit expression for V ( x ). From (11) and (14), we get the potentials in terms of thesolutions w ( x ) to the fourth Painlev´e equation (16) as V ( x ) = x − (cid:0) w (cid:48) − xw − w (cid:1) − , (18) V ( x ) = x + ( w (cid:48) + 2 xw + w ) − . (19)The latter reveals that V ( x ) and V ( x ) are regular functions , that is, free of singularities,as long as w ( x ) is a regular function as well. From now on, we focus on the properties of w ( x ) that guarantee the desired regularity.The fourth Painlev´e transcendent (16) has been extensively studied in the litera-ture [25, 27, 29]. Following some previous reports on the construction of third-ordershape invariant quantum models [35, 47], we particularly focus on the rational solutionsof Eq. (16). Such a case has been studied for only some specific values of the parameters α and β , where a connection with a two-step Darboux-Crum transformation was foundin [35]. So far, a generalization for arbitrary values of the parameters is still absent. Thecomplete set of rational solutions to the fourth Painlev´e equation are classified into threehierarchies [27, 29], namely the “ − /x hierarchy”, “ − x hierarchy”, and the “ − x/ generalized Hermitepolynomials [29], which have been already considered in previous works to obtain regularpotentials for some particular rational extensions of the harmonic oscillator [36, 46].In this work, we focus on the “ − x/
3” hierarchy of rational solutions that lead to reg-ular physical models. Rational solutions in such hierarchy are divided into three differenttypes of solutions of the form [26, 27] w [1] m,n ( x ) = − x ddx ln Q m +1 ,n Q m,n , α = 2 m + n , β = − (cid:18) n − (cid:19) , (20) w [2] m,n ( x ) = − x ddx ln Q m,n Q m,n +1 , α = − m − n , β = − (cid:18) m − (cid:19) , (21) w [3] m,n ( x ) = − x ddx ln Q m,n +1 Q m +1 ,n , α = n − m , β = − (cid:18) m + n + 13 (cid:19) , (22)where Q m,n ≡ Q m,n ( x ) stands for the generalized Okamoto polynomials [54], computedfrom the nonlinear recurrence relations Q m +1 ,n Q m − ,n = 92 (cid:2) Q m,n Q (cid:48)(cid:48) m,n − ( Q (cid:48) m,n ) (cid:3) + (cid:2) x + 3(2 m + n − (cid:3) Q m,n , (23) Q m,n +1 Q m,n − = 92 (cid:2) Q m,n Q (cid:48)(cid:48) m,n − ( Q (cid:48) m,n ) (cid:3) + (cid:2) x + 3(1 − m − n ) (cid:3) Q m,n , (24)6 Q k ( x )0 11 12 2 x + 33 8 x + 60 x + 90 x + 1354 64 x + 1344 x + 9360 x + 30240 x + 56700 x + 170100 x + 1275755 1024 x + 46080 x + 817920 x + 7603200 x + 41731200 x + 155675520 x + 493970400 x +1886068800 x + 5304568500 x + 5304568500 x + 3978426375 Table 1: Conventional Okamoto polynomials Q k ( x ) computed from (23) for m = k , n = 0,and the initial conditions Q = Q = 1.for m, n = 0 , , · · · , together with the initial conditions Q , = Q , = Q , = 1 and Q , = √ x . Moreover, with the aid of the B¨acklund transformations, the rationalsolutions (20)-(22) can be cast into the new equivalent forms (see App. A for details) w [1] m,n ( x ) = − √ Q m +1 ,n − Q m,n +1 Q m,n Q m +1 ,n , (25) w [2] m,n ( x ) = − √ Q m +1 ,n Q m − ,n +1 Q m,n +1 Q m,n , (26) w [3] m,n ( x ) = − √ Q m +1 ,n +1 Q m,n Q m +1 ,n Q m,n +1 . (27)Among the generalized Okamoto polynomials, only the Okamoto polynomials Q k ≡ Q k, have zeros distributed outside of the real line. For details, see the zeros distributionin [29] and the analysis provided in Sec. 5.2. Thus, from Eqs. (20)-(22), regular rationalsolutions to the fourth Painlev´e equation ˜ w k ( x ) are obtained only if˜ w k ( x ) ≡ w ( x ) ≡ w [1] k, ( x ) = − x ddx ln Q k +1 Q k , α ≡ ˜ α = 2 k , β ≡ ˜ β = − / , (28)for k = 0 , , · · · . Interestingly, with the aid of the B¨acklund and Schlesinger transforma-tions [28] (see App. A for details), along with the identities (25)-(27), the potential V ( x )defined in (18) takes the form V ( k ) ( x ) ≡ V ( x ) = x − (cid:0) ˜ w (cid:48) k − x ˜ w k − ˜ w k (cid:1) − x − Q k +2 Q k Q k +1 + 4 k + 1 , (29)and the recurrence relation (23) leads to the equivalently expression V ( k ) ( x ) ≡ x − d dx ln Q k +1 + 4 k − . (30)Eqs. (29)-(30) are the most general forms of the rationally extended potentials generatedfrom the “ − x/
3” hierarchy. In particular, using the Okamoto polynomials presentedin Tab. 1, we recover the potential V ( k =1) ( x ) reported in [35], and V ( k =2) ( x ) reportedin [47]. Besides the generality of the potential in (29), we can explicitly determine theeigenfunctions φ n and eigenvalues E n in (1). Details are presented in the following sections.7 Mapping operators: zero-modes and eigenvalues
As discussed in the previous section, the third-order intertwining operators A and A † have been factorized in terms of the second-order and first-order operators introducedin (5) to explicitly show the appearance of the fourth Painlev´e equation and its relationto the potential V ( x ). Nevertheless, to study the spectral information of H , it is moreconvenient and simple if we further decompose the second-order operators as M † = M † M † and M = M M , where the first-order operators M , and M † , are of the form M := − ddx + W ( x ) , M := − ddx + W ( x ) , (31) M † := ddx + W ( x ) , M † := ddx + W ( x ) , (32)where the real-valued functions W , ( x ) are computed from the factorization M † = M † M † after comparing term by term with (6), and using (14). The particular factorizationdiscussed above corresponds to the reducible case discussed in [38,51]. The straightforwardcalculations lead to [46, 47] W ± ( x ) ≡ W ( x ) = − G + G (cid:48) ± √− d G , W ± ( x ) ≡ W ( x ) = − G − G (cid:48) ± √− d G , (33)which, with the aid of (17), take the more convenient form W ± ( x ) = ˜ F − k ± (cid:113) − β w + x , W ± ( x ) = − ˜ F + k ± (cid:113) − β w + x , (34)where ˜ F ± k = ˜ w (cid:48) k ± (2 x ˜ w k + ˜ w k ), with ˜ w k given in Eq. (28). In this form, the B¨acklundtransformations introduced in (A-1), combined with w = ˜ w k , allows us to simplify (34)into W +1 = x ddx ln Q k, Q k +1 , W − = x ddx ln Q k +1 , − Q k +1 , (35) W +2 = x ddx ln Q k Q k, , W − = x ddx ln Q k Q k +1 , − . (36)Notice that W − and W − depend on generalized Okamoto polynomials Q m, − that con-tain a negative index, which are determined iteratively from the recurrence relation (24)with n = 0. For instance, Q , − is computed from Q , = Q , = 1, Q , − from Q , and Q , , Q , − from Q , and Q , , and so on. In Sec. 2, it was shown that the factorization operators M , M † , Q , and Q † define inter-twining relations between H and the auxiliary Hamiltonian H . Similar relations hold for8 H + 2 λH + 2 λ H + 2 λQ † A † = Q † M = Q † M M M = M M M M Figure 1:
Third-order shape-invariant Hamiltonian H . The arrows indicate the intertwining relationamong the different Hamiltonians H , H , and H . For instance, the arrow on top indicates HA † = A † ( H + 2 λ ), as given in (2). The direction of arrows is inverted by using the adjoint relations. the factorized operators introduced in (31)-(32). In analogy to the procedure followed inthe previous section, we introduce another auxiliary Hamiltonian H , together with thecorresponding eigenvalue equation, H := − d dx + V ( x ) , H ϕ (2) n = E (2) n ϕ (2) n , (37)with ϕ (2) n and E (2) n the eigenfunctions and eigenvalues, respectively, so that H intertwineswith H and H through M , M † , M , and M † as H M = M H , H M † = M † H , (38) H M = M H , HM † = M † H . (39)Clearly, the combined action of (38)-(39) brings us back to the initial intertwining relationintroduced in (10). Thus, for clarity, the action of the intertwining relations constructedso far is depicted and summarized in Fig. 1, where it is also evident that A † and A are ladder operators in H . Furthermore, Fig. 1 reveals that the elements of the set { M , M † , M , M † , Q, Q † } define mapping operators among the vector spaces generated bythe eigenfunctions of H and the auxiliary Hamiltonians H and H . We thus get M † : H → H , M : H → H ,M † : H → H , M : H → H ,Q † : H → H , Q : H → H . (40)In the latter, the vector spaces are defined as H = Span { φ n } ∞ n =0 , H = Span { ϕ (1) n } ∞ n =0 , H = Span { ϕ (2) n } ∞ n =0 . (41)Before concluding this section, it is worth noticing that the relations (38)-(39) can befurther exploited by recalling that M , and M † , are first-order operators, and thus canbe exploited to get [47] H = M † M + (cid:15) , H = M M † + (cid:15) ,H = M † M + (cid:15) , H = M M † + (cid:15) , (42)9or some real constants (cid:15) and (cid:15) determined after substituting Eqs. (31)-(32) and Eqs. (34)into (42). After some calculations we obtain (cid:15) = γ − √− d = 2 k − , (cid:15) = γ + √− d = 2 k − . (43) φ ( k )0; j With the operations depicted in Fig. 1, we determine the spectral information of H .Although the eigenvalues of the auxiliary Hamiltonians H and H can be obtained in asimilar way, the latter is not relevant for the present work and thus will be omitted. Letus recall that A = M † M † Q is an annihilation operator, which is also of third-order, andthus three different eigenfunctions are simultaneously annihilated by it, that is, Aφ ( k )0; j = M † Qφ ( k )0; j = M † M † Qφ ( k )0; j = 0 , j = 1 , , . (44)Henceforth, the eigenfunctions φ ( k )0; j will be called zero-modes . Those functions are com-puted from the mappings defined in (40) along with the factorization of the Hamiltonians H , H , and H . From (44) we distinguish three different cases: • φ ( k )0;1 : The zero-mode φ ( k )0;1 ∈ H is annihilated by Q , that is, Qφ = 0. • φ ( k )0;2 : Q does not annihilate the zero-mode, Qφ ( k )0;2 = f (cid:54) = 0, instead we have M † f =0. From (42) it follows that a function annihilated by M † is a zero-mode of H , say f ≡ ϕ (1)0 ∈ H , which can be mapped into the unnormalized zero-mode φ ∈ H through the action Q † ϕ (1)0 . • φ : Neither Q nor M † annihilate the zero-mode, M † Qφ = f (cid:54) = 0, instead wehave M † f = 0. From (42) we conclude that a function annihilated by M † is azero-mode of H , say f ≡ ϕ (2)0 ∈ H . Thus, the respective unnormalized zero-mode φ ∈ H is determined with the mapping Q † M ϕ (2)0 .In Fig. 2, we depict a diagram summarizing the action of the mapping operatorsrequired to compute the zero-modes. The straightforward calculations lead to φ ( k )0;1 ( x ) = N e (cid:82) dxW ( x ) ,φ ( k ) ± ( x ) = N (cid:0) W − W ± (cid:1) e − (cid:82) dxW ± ( x ) ,φ ( k ) ± ( x ) = N (cid:16) ± √− d + ( W − W ± )( W ± + W ± ) (cid:17) e − (cid:82) dxW ± ( x ) , (45)with the respective eigenvalues E = 0 , E ± = γ + 2 ∓ √− d , E ± = γ + 2 ± √− d . (46)10 H H φ ( k )0;3 φ ( k )0;2 φ ( k )0;1 Q ϕ (1)1 ϕ (1)0 M † ϕ (2)0 M † Q † Q † M M Figure 2:
Mappings between the zero-mode eigenfunctions φ ( k )0; j , for j = 1 , ,
3, and the eigenfunctionsof the auxiliary Hamiltonians H and H . The dotted-black lines denote the null-vector. Notably, the zero-mode φ is easily determined in terms of ˜ w k through (14) and W = − ( ˜ w k + x ) = − (cid:18) x ddx ln Q k +1 Q k (cid:19) . (47)We thus get φ ( k )0;1 ( x ) = N (cid:32) e − x Q k +1 (cid:33) Q k , E = 0 , (48)where the Okamoto polynomials Q k are presented in Table. 1 for several values of k . Itis worth recalling that Q k are nodeless functions for x ∈ R . Thus, the zero-mode φ isboth a regular and nodeless function. Thus, φ is the ground eigenfunction of H , with E the lowest eigenvalue among all the sequences. Similarly, we compute the remainingzero-modes from (36), leading to φ ( k )+0;2 ( x ) = N +0;2 (cid:32) e − x Q k +1 (cid:33) Q k +1 , , φ ( k ) − ( x ) = N − (cid:32) e − x Q k +1 (cid:33) Q k +2 , − , (49) φ ( k )+0;3 ( x ) = N +0;3 (cid:32) e − x Q k +1 (cid:33) Q k +2 , − , φ ( k ) − ( x ) = N − (cid:32) e − x Q k +1 (cid:33) Q k +1 , , (50)where N ± and N ± stand for the normalization constants. From (49)-(50), it is clear that φ ( k ) ± = φ ( k ) ∓ and E ± = E ∓ . Henceforth, we can freely fix, without loss of generality, thezero-modes and eigenvalues with the “+” superscript. To illustrate our results, we depictthe behavior of the zero-modes φ ( k )0; j and the corresponding potential V ( k ) ( x ) in Fig. 3 forseveral values of k .So far, we have explicitly determined three regular zero-modes. Moreover, it is imme-diate to realize that each zero-mode converges asymptotically to zero, lim x →±∞ φ j → j = 1 , ,
3. Thus the zero-modes become elements in the space of square-integrable11 a) V ( k ) ( x ) (b) k = 2 (c) k = 4 Figure 3: (a) Potential V ( k ) ( x ) given in (29) for k = 2 (solid-blue), k = 3 (dashed-green),and k = 4 (dotted-red). The zero-modes φ (solid-blue), φ (dashed-green), and φ (dotted-red) are depicted for k = 2 (b) and k = 4 (c). k Q k, ( x )0 11 √ x x + 12 x − √ x (16 x + 192 x + 504 x − x + 7680 x + 80640 x + 362880 x + 453600 x − x − x − x +80372255 √ x (4096 x + 245760 x + 5990400 x + 77414400 x + 569721600 x + 2246952960 x +1600300800 x − x − x − x − x +3258331201125) Table 2: Generalized Okamoto polynomials Q m,n computed from (23)-(24) for m = k , n = 1, and the initial conditions Q , = 1 and Q , = √ x . k Q k, − ( x )0 2 x + 31 √ x x − √ x (4 x − x + 240 x − x − x − x + 60755 √ x (256 x + 6144 x + 34560 x − x − x − x − x + 63149625) Table 3: Generalized Okamoto polynomials Q m,n computed from (23)-(24) for m = k , n = −
1, and the initial conditions Q , − = 2 x + 3 and Q , − = √ x .12unctions, that is, φ j ∈ L ( dx ; R ) such that( φ ( k )0; j , φ ( k )0; j ) := (cid:90) ∞−∞ | φ ( k )0; j ( x ) | < ∞ , j = 1 , , , (51)with the inner product defined for any two eigenfunctions φ and ˜ φ as( ˜ φ, φ ) := (cid:90) ∞−∞ dx [ ˜ φ ( x )] ∗ φ ( x ) . (52)In some cases, the eigenfunctions Ψ ( k )0; j annihilated by the creation operator A † maylead to physical solutions as well. However, in this particular case it is straightforward torealize that any function Ψ j , such that A † Ψ ( k )0; j = 0, diverges asymptotically at x → ±∞ and thus Ψ ( k )0; j (cid:54)∈ L ( dx ; R ), for j = 1 , ,
3. See [47] for additional details. Therefore, nophysical eigenfunctions of A † are extracted from the “ − x/
3” hierarchy, and the spectrumis solely determined from the zero-modes (48)-(50). φ ( k ) n ; j Given that A † is a creation operator in H , the higher-modes φ n ; j are then computedthrough the iterated action on the zero-modes φ j , up to a normalization constant,through φ n ; j ∝ A † n φ j , for j = 1 , ,
3. From (2), it follows that the eigenvalues areincreased by two units after each iteration. In this form, the set of eigenvalues are de-termined straightforwardly from E = 0, E = 2 k + 2 /
3, and E = 2 k + 4 /
3, leadingto E n ;1 = 2 n , E n ;2 = 2 k + 2 n + 2 / , E n ;3 = 2 k + 2 n + 4 / . (53)Notice that none of the eigenvalues overlaps; thus, the respective eigenfunctions shouldbe, in general, all different. Since the creation operator A † is represented by a third-orderdifferential operator, any attempt to determine closed expressions for the higher-modesbecomes a challenging task.In this section, we construct the higher-modes through an alternative mechanism thatrelies on a three-term recurrence relation, similar to the usual construction for classicalorthogonal polynomials [55, 56]. To this end, we notice that, among all the zero-modes φ ,j , there is a common factor e − x /Q k +1 ( x ) that does not depend on the index j . Wethus regard such a factor as the weight function µ k ( x ) in the set of eigenfunctions of H so that φ ( k ) n ; j ( x ) = µ k ( x ) N n,j P ( k ) n ; j ( x ) , µ k ( x ) := e − x Q k +1 ( x ) , j = 1 , , , (54)with N n ; j the normalization constant, and P ( k ) n ; j ( x ) some functions to be determined fromthe following constraints: P ( k )0;1 ( x ) := Q k ( x ) , P ( k )0;2 ( x ) := Q k +1 , ( x ) , P ( k )0;3 ( x ) := Q k +2 , − ( x ) . (55)13n this form, for n = 0, we recover the zero-modes φ j obtained in Sec. 3.2.The weightfunction µ k ( x ) has support over the real line, x ∈ R . It is worth to notice that the spaceof solutions H can be decomposed as the direct sum of the subspaces H = V ⊕ V ⊕ V , V j = Span { φ n ; j } ∞ n =0 , j = 1 , , , (56)where the ladder operators A and A † act as endomorphisms on the elements of the sub-spaces V j , that is, A : V j → V j and A † : V j → V j for all j = 1 , ,
3, see Fig. 4.We can exploit the latter fact to obtain additional properties by working in eachvector subspace. In particular, from the eigenvalue equation (1) we get a linear second-order differential equation for each set of functions { P ( k ) n ; j } ∞ n =0 . It is useful to notice thatthe weight function µ k satisfy µ (cid:48) k µ k = − x − ddx ln Q k +1 , µ (cid:48)(cid:48) k µ k = − − d dx ln Q k +1 + (cid:18) x ddx ln Q k +1 (cid:19) , (57)from which we determine the differential equation for the functions P ( k ) n ; j as d P ( k ) n ; j dx − (cid:18) x Q (cid:48) k +1 Q k +1 (cid:19) dP ( k ) n ; j dx + (cid:18) Q (cid:48)(cid:48) k +1 Q k +1 + 2 x Q (cid:48) k +1 Q k +1 − k E n ; j (cid:19) P ( k ) n ; j = 0 , (58)for j = 1 , ,
3, and n = 0 , , · · · . Now, we present another of the main results of this work, that is, the proper constructionof a linear three-term recurrence relation for the polynomials P ( k ) n ; j , for j = 1 , ,
3, given inEq. (54). Such recurrence relations are fundamental in the theory of classical orthogonalpolynomials [55, 56], a property contained within the
Favard’s theorem [57], and in thiscase are achieved by exploiting the action of the third-order ladder operators A and A † on the elements of V j , for j = 1 , ,
3, instead of the whole space H . Thus, let us considerthe action of the ladder operators A and A † on any arbitrary element of V j . Given thatthe ladder operators are mutually adjoint, we obtain A † φ ( k ) n ; j = C n +1; j φ ( k ) n +1; j , Aφ ( k ) n ; j = C n ; j φ ( k ) n − j , (59)where C n ; j is a proportionality constant computed from( φ ( k ) n ; j , AA † φ ( k ) n ; j ) = | C n +1; j | . (60)To determine the latter, we use the factorization AA † = M † M † QQ † M M along with theintertwining relations (38)-(39) and (42). We thus get AA † = ( H − (cid:15) )( H − (cid:15) )( H + 2),14 V V V ≡ ⊕ ⊕ k k + 22 k + k + k + k + φ ( k )0;1 φ ( k )1;1 kφ ( k ) k ;1 k + 2 φ ( k ) k +1;1 A † A k + φ ( k )0;2 k + φ ( k )1;2 A † A k + φ ( k )0;3 k + φ ( k )1;3 A † A Figure 4: (Color online) Eigenfunctions φ ( k ) n of H arranged according to the increasing eigenvalues E n ,together with the subspace decomposition H ≡ V ⊕ V ⊕ V introduced in (56). The eigenfunctions φ ( k ) n ; j , for j = 1 , , , and n = 0 , , · · · , are given in (54), and computed through the three-term recurrencerelation (70). and, with the aid of (43), we obtain the proportionality constants as C n +1;1 := (cid:115) ( n + 1) (cid:18) n − k + 23 (cid:19) (cid:18) n − k + 13 (cid:19) , (61) C n +1;2 := (cid:115) ( n + 1) (cid:18) n + 23 (cid:19) (cid:18) n + k + 43 (cid:19) , (62) C n +1;3 := (cid:115) ( n + 1) (cid:18) n + 43 (cid:19) (cid:18) n + k + 53 (cid:19) . (63)From the latter, together with Eq. (59), the normalization constants N n ; j are determinedin terms of the normalization constants N j of the corresponding zero-mode as N n +1; j := N j n (cid:89) p =0 C p +1; j , n = 0 , , · · · , j = 1 , , . (64)Moreover, from the explicit action of the ladder operators A and A † , combined with15he differential equations obtained in Eq. (58), we determine the following set of equations: G n ; j ( x ) dP ( k ) n ; j dx + (cid:20) − Q (cid:48) k Q k G n ; j ( x ) + E n ; j w [1] k, (cid:21) P ( k ) n ; j = L n ; j P ( k ) n − j , (65) G n +1; j ( x ) dP ( k ) n ; j dx − (cid:20)(cid:18) x Q (cid:48) k, Q k, (cid:19) G n +1; j ( x ) + (cid:18) − k + E n ; j (cid:19) w [3] k, (cid:21) P ( k ) n ; j = − ˜ L n ; j P ( k ) n +1; j , (66)where G n ; j ( x ) := 29 Q k +2 Q k Q k +1 − E n ; j ≡ Q k +1 , Q k +1 , − Q k +1 + 23 + 2 k − E n ; j , (67) L j := 0 , L n +1; j := 1 , ˜ L n ; j := ( C n +1; j ) , n = 0 , , · · · , j = 1 , , . (68)By subtracting Eq. (65) with Eq. (66), and after performing some calculations, we get˜ L j P ( k )1; j ( x ) = (cid:20) − w [2] k, ( x ) G j ( x ) + E j w [1] k, ( x ) G j ( x ) G j ( x ) + w [3] k, ( x ) (cid:18) − k + E j (cid:19)(cid:21) P ( k )0; j ( x ) , (69)together with the three-term recurrence relation˜ L n +1; j P ( k ) n +2; j ( x ) + G n +2; j ( x ) G n +1; j ( x ) P ( k ) n ; j ( x )= (cid:20) − w [2] k, ( x ) G n +2; j ( x ) + E n +1; j w [1] k, ( x ) G n +2; j ( x ) G n +1; j ( x ) + w [3] k, ( x ) (cid:18) − k + E n +1; j (cid:19)(cid:21) P ( k ) n +1; j ( x ) , (70)for j = 1 , , n = 0 , , · · · , where w [ p ] k, are given in Eqs. (20)-(22) for p = 1 , , P ( k )0; j are given in (55). In terms of finite-difference calculus, wenotice that Eq. (70) is a second-order finite-difference equation that requires two boundaryconditions to determine P ( k ) n ; j uniquely. However, for n = 0, we have L j = 0 and thus P ( k )1; j is determined exclusively from P ( k )0; j , for j = 1 , ,
3. The latter means that we onlyrequire P ( k )0 ,j as the initial condition for each polynomial sequence, the remaining terms arethen computed recursively. To illustrate our results, we present in Tables 4-5 the explicitform of P ( k ) n,j for several values of k and n . Such terms are computed using the generalizedOkamoto polynomials given in Tables 1-3.From the three-term recurrence relation in (70), it seems that P ( k ) n ; j takes the form of arational function. Nevertheless, after computing explicitly some terms, it follows that P ( k ) n,j leads indeed to polynomials. Such a computation is presented in Tables 4-5, for severalvalues of n and k = 1 ,
3. Clearly, the latter is by no means a proof about the polynomialstructure of P ( k ) n ; k , and a more detailed analysis is required. This is a fact discussed in thefollowing sections. 16 P ( k =1) n ;1 ( x ) / D ( k =1) n ;1 P ( k =3) n ;1 / D ( k =3) n ;1 x + 60 x + 90 x + 1351 x (2 x + 9) x (16 x + 480 x + 3528 x + 7560 x + 8505)2 8 x − x − x + 81 64 x + 2496 x + 35280 x + 211680 x +510300 x + 510300 x − x (16 x − x + 1944 x + 4860 x − x (128 x + 4416 x + 58464 x + 390960 x +1370520 x + 2143260 x − x − x − x + 73872 x − x − x + 2755620 x − x + 8448 x − x − x − x − x − x +165337200 x + 310007250 x − n P ( k =1) n ;2 ( x ) / D ( k =1) n ;2 P ( k =3) n ;2 / D ( k =3) n ;2 x + 12 x − x + 7680 x + 80640 x + 362880 x +453600 x − x − x − x + 80372251 x (cid:0) x − x − x + 567 (cid:1) x (512 x +3840 x − x − x − x + 24766560 x + 106142400 x +445798080 x + 208967850 x − x − x + 10080 x − x + 25515 2048 x − x − x +13996800 x +131155200 x − x − x − x − x + 56421319500 x +197474618250 x − x (64 x − x + 121680 x − x +1326780 x + 7960680 x − x (4096 x − x + 4773888 x +85432320 x − x − x + 58205226240 x +430768316160 x − x +1715208112800 x − x − x + 21708102677625)4 256 x − x + 1587456 x − x + 283551840 x − x − x + 4728643920 x − x − x + 102371328 x − x − x +383543424000 x + 1297481898240 x − x − x +351830022223680 x − x +5263973698183200 x + 955156517815500 x − x + 2149102165084875 Table 4: Polynomials P ( k ) n ; j computed through the three-term recurrence relation (70) for j = 1 , n = 0 , , , ,
4. The constants D ( k ) n ;1 and D ( k ) n ;2 are integer scaling factors.17 P ( k =1) n ;3 ( x ) / D ( k =1) n ;3 P ( k =3) n ;3 / D ( k =3) n ;3 x (4 x − x (256 x + 6144 x + 34560 x − x − x − x − x +63149625)1 16 x − x + 360 x + 3240 x − x − x − x − x +18506880 x + 155675520 x + 314344800 x +808315200 x − x − x +17050398752 x (32 x − x + 20592 x − x − x + 280665) x (2048 x − x − x +27530496 x +66624768 x − x − x + 6870679200 x +1717669800 x + 378746190900 x +405799490250 x − x − x + 376992 x − x +15717240 x + 23575860 x − x +22733865 8192 x − x + 19353600 x +102574080 x − x − x + 391975718400 x +884443795200 x − x − x − x +14608781649000 x + 295827828392250 x − x (256 x − x + 2193408 x − x + 661620960 x − x +20096736660 x − x (16384 x − x +149409792 x − x − x + 819467228160 x − x − x +22612029854400 x + 1482481454126400 x − x +12480073465860000 x − x − x + 44226260344641375) Table 5: Polynomials P ( k ) n ; j computed through the three-term recurrence relation (70) for j = 3 and n = 0 , , , ,
4. The constants D ( k ) n ;3 are integer scaling factors. So far, we have identified the potential V ( k ) ( x ) written in terms of the Okamoto polynomi-als and constrained to the third-order shape-invariant condition (2). In contradistinctionto previous works, the spectral problem was determined in general for any value of thearbitrary parameter k = 0 , , · · · . Moreover, a mechanism to generate the higher-modesthrough a three-term recurrence relation, for arbitrary k , was introduced. In this section,for completeness, we discuss the connection from our previous results to the higher-orderSUSY QM construction. It is worth to mention that, for the construction discussed in [35],it was found that the resulting potential was related to a two-step SUSY construction ofthe harmonic oscillator. Such a case corresponds to k = 1 in our construction. We thuscomplement such a discussion, and generalize it for any values of k .Before proceeding, we explore the asymptotic behavior of the rational term Q k +2 Q k Q k +1 inthe potential V ( k ) ( x ) obtained in (29), where Q k ≡ Q k, ,. It is well-known [28, 54] that18 m,n is a polynomial of degree d m,n = m + n + mn − m − n , and thus the rationalpart of the potential V ( k ) is a polynomial of degree 2( k + k + 1) and 2 k ( k + 1) in thenumerator and denominator, respectively. The latter implies that the whole potential V ( k ) ( x ) behaves asymptotically as a polynomial of degree two for x → ±∞ . Thus, fromOblomkov’s theorem [58], it follows that V ( k ) ( x ) can be constructed from the harmonicoscillator through several rational Darboux transformations.The generalized Okamoto polynomials admit a Wronskian representation obtainedfrom the relationship between the Schur polynomials and the τ function defining theToda-type equations (23)-(24). See [59–61] for further details. Such a Wronskian represen-tation allows us identifying the corresponding elements required to perform the Darbouxtransformation to generate V ( k ) ( x ). Thus, from [28] we have Q m,n = C m,n Wr(Ψ , · · · , Ψ m +3 n − , Ψ , · · · , Ψ n − ) , n = 1 , , · · · , m = 0 , , · · · , (71)where the proportionality constant C m,n is fixed so that (23)-(24) are fulfilled, and Wr( · )denotes the usual Wronskian determinantWr( f , f , · · · , f n ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f f · · · f n f (cid:48) f (cid:48) · · · f (cid:48) n ... ... . . . ... f ( n − f ( n − · · · f ( n − n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , f ( n ) ≡ d n fdx n . (72)The functions ψ n ( x ) are determined from the generating function [28] e xξ +3 ξ = ∞ (cid:88) r =0 Ψ r ( x ) ξ r , Ψ r ( x ) := 13 r/ r ! H r (cid:18) x √ (cid:19) , (73)with H r ( z ) = i r H r ( iz ) and H r ( z ) the pseudo-Hermite and Hermite polynomials [53], re-spectively. Notice that C m,n does not modify the form of V ( k ) ( x ) as such a constantvanishes inside the logarithmic derivative. Moreover, C m,n can be absorbed by the zero-modes normalization factor, and thus its explicit form is not required throughout the restof the manuscript.The Wronskian representation for the conventional Okamoto polynomials Q m ≡ Q m, must be handled with caution, as the length of the partition defining the corresponding τ function used in (71) is not the appropriate one, see [59]. In such a case we get Q m = C m Wr(Ψ , · · · , Ψ m − ) , m = 0 , , · · · . (74)Notice that fixing n = 0 in (71) does not lead to the right representation given in (74).On the other hand, the generalized Okamoto polynomials Q m,n possess a symmetry [28]of the form Q m,n ( x ) = exp( − i πd m,n ) Q n,m ( ix ), with d m,n = m + n + mn − m − n thedegree of the Okamoto polynomial. This symmetry allows us getting the alternativerepresentation Q m,n = D m,n Wr( ψ , ψ , · · · ψ n +3 m − , ψ , ψ , · · · Ψ m − ) , (75)19here ψ r ( x ) := ( − i ) r Ψ r ( − ix ) = 13 r/ r ! H r (cid:18) x √ (cid:19) , (76)and D m,n is a proportionality constant. In this form, we have a Wronskian representa-tion in terms of pseudo-Hermite polynomials (71) and one in terms of Hermite polyno-mials (75). Additional representations have been recently reported in terms of pseudo-Wronskians , which are composed of a combination of Hermite and pseudo-Hermite polyno-mials [62]. During the rest of this manuscript, we focus mostly on the representation (75)as it reproduces the right results for n = 0 and n = − V ( k ) ( x )and the higher-order Darboux transformations. To this end, let us introduce the followingshort and convenient notation: W { i , ··· ,i m } = Wr ( ψ i , · · · , ψ i m ) , W { i , ··· ,i m } = Wr (Ψ i , · · · , Ψ i m ) , (77)together with˜ W { i , ··· ,i m } = Wr (cid:16) ˜ ψ i , · · · , ˜ ψ i m (cid:17) , ˜ W { i , ··· ,i m } = Wr (cid:16) ˜Ψ i , · · · , ˜Ψ i m (cid:17) , (78)˜ ψ r := e − x ψ r = e − x r/ r ! H r (cid:18) x √ (cid:19) , ˜Ψ r := e x Ψ r = e x r/ r ! H r (cid:18) x √ (cid:19) . (79)With this new notation, the conventional Okamoto polynomials Q k +1 can be rewritten as Q k +1 = C k +1 W { , ··· , k − , , ··· , k − } = C k +1 e k x ˜ W { , ··· , k − , , ··· , k − } , (80) Q k +1 = D k +1 W { , ··· , k − } = D k +1 e − k x ˜ W { , ··· , k − } , (81)where we have used the Wronskian identity Wr( gf , · · · , gf n ) = g n Wr( f , · · · , f n ) [63].Therefore, the potential V ( k ) ( x ) given in (30) takes the two equivalent forms V ( k ) ( x ) ≡ x − d dx ln ˜ W { , , ··· , k − , , ··· k − } − , (82) V ( k ) ( x ) ≡ x − d dx ln ˜ W { , , ··· k − } + 2 k − . (83)These expressions reveal that the third-order shape invariant construction (2) is relatedto a higher-order SUSY transformation. In particular, (82) shows that H given in (2) isrelated to a state-deleting SUSY transformation of the harmonic oscillator (See App. C fordetails) with the physical energy levels { E ( osc ) n } n ∈S d deleted, where S d = { , , , , · · · k − , k − } . Alternatively, from the Wronskian representation (83), we obtain a state addingconstruction in which the non-physical energy levels { E ( osc ) − n } n ∈S a have been added, with S a = { , , · · · k − } . In this form, we verify the well-known equivalence between thestate-adding and the state-deleting Darboux transformations already noticed by Felder etal. [64]. 20astly, using (71) and the action of the creation operator A † , the higher-modes can bealternatively determined from the Wronskian representation, leading to (see App. B fordetails) φ ( k ) n ;1 ( x ) ≡ µ k ( x ) (cid:113) ˜ N n ;1 W I ( k ) ∪{ n } , φ ( k ) n ;2 ( x ) ≡ µ k ( x ) (cid:113) ˜ N n ;2 W I ( k ) ∪{ n +3 k +1 } ,φ ( k ) n ;3 ( x ) ≡ µ k ( x ) (cid:113) ˜ N n ;3 W I ( k ) ∪{ n +3 k +2 } , (84)with I ( k ) = { , , , , · · · , k − , k − } , µ k ( x ) the weight function given in (54), and˜ N n ; j the corresponding normalization constants, which are different to the ones obtainedin (64). We thus have two different forms to represent the higher-modes. Furthermore,after comparing (54) with (84), it is straightforward to realize that P ( k ) n ;1 ( x ) = (cid:115) ˜ N n ;1 N n ;1 W I ( k ) ∪{ n } ( x ) , P ( k ) n ;2 ( x ) = (cid:115) ˜ N n ;2 N n ;2 W I ( k ) ∪{ n +3 k +1 } ( x ) ,P ( k ) n ;3 ( x ) = (cid:115) ˜ N n ;3 N n ;3 W I ( k ) ∪{ n +3 k +2 } ( x ) , (85)for n, k = 0 , , · · · . From the previous Wronskian representation, it is clear that thefunctions P ( k ) n ; j computed from (70) are indeed polynomials. The Wronskian representation obtained in the previous section provides us with an ad-ditional relation that allows relating the polynomials emerging from the higher-modeswith the family of orthogonal exceptional Hermite polynomials [48]. To this end, let usintroduce the non-decreasing sequence of integers, 0 ≤ λ ≤ · · · ≤ λ (cid:96) , together with thepartition and double partition λ := ( λ , · · · , λ (cid:96) ) , λ := ( λ , λ , · · · , λ (cid:96) , λ (cid:96) ) , (86)respectively. The exceptional Hermite polynomials are defined as [48] H ( λ ) n ( x ) ≡ H λ ,n ( x ) := Wr[ H ν , H ν +1 , · · · H ν (cid:96) , H ν (cid:96) +1 , H n ] , n (cid:54)∈ { ν , ν + 1 , · · · , ν (cid:96) , ν (cid:96) + 1 } , (87)where H m := H m ( x ) and the indexes inside the Wronskian are related to the elements ofthe partition through ν j = λ j + j − , ≤ ν < · · · < ν (cid:96) , j = 1 , · · · (cid:96). (88)From the latter, it is clear that the polynomials P ( k ) n ; j becomes a particular case of theexceptional Hermite polynomials (87). It is worth noticing that the argument of the Her-mite polynomials inside the Wronskian (85) is different to that of the Hermite polynomials21n the Wronskian (87). Nevertheless, a simple reparametrization of the form x → √ x solves the issue, leading to the relation H ( λ k ) σ n ; j ( x ) ≡ D ( k ) n ; j W I ∪{ σ n ; j } ( √ x ) = D ( k ) n ; j (cid:115) N n ; j ˜ N n ; k P ( k ) n ; j ( √ x ) , n = 0 , , · · · , (89)where we have used the index notation σ n ;1 := 3 n , σ n ;2 := 3 n + 3 k + 1 , σ n ;3 = 3 n + 3 k + 2 . (90)with the double partition λ k := (1 , , , , · · · , k, k ), or equivalently the partition λ k :=(1 , , · · · , k ), so that ν = { ν p } kp =1 = { , , , , · · · , k − , k − } , for k = 0 , , · · · . Onthe other hand, the proportionality constant D ( k ) n ; j becomes D ( k ) n ; j := 3 k ( k +1)4 + k + σn ; j (cid:32) k (cid:89) p =1 ν p ! (cid:33) ( σ n ; j !) . (91)Recursion formulas have been previously constructed for the families of exceptionalHermite polynomials [48–50] through the use of higher-order recurrence relations. In par-ticular, for any family of general exceptional Hermite polynomial H ( λ ) n defined in termsof an ( (cid:96) + 1)-order Wronskian (87), there is a (2 (cid:96) + 3)-order recurrence relation, with( (cid:96) + 1) initial conditions, that generates the complete sequence of exceptional polynomi-als. See [48] for more details. From the latter, it follows that the particular sequence λ k defining H ( λ k ) σ n ; j in (89) should be computed from a (2 k + 3)-order recurrence relation.Nevertheless, from the recursion formulas (70) obtained for the polynomials P ( k ) n ; j , we canintroduce a set of new linear and third-order recurrence relations for the exceptional Her-mite polynomials defined by the three disjoint sets S H := { H ( λ k ) σ n ;1 } ∞ n =0 , S H := { H ( λ k ) σ n ;2 } ∞ n =0 ,and S H := { H ( λ k ) σ n ;3 } ∞ n =0 . That is, the rescaled exceptional Hermite polynomials definedin (89) satisfy a three-term recurrence relation akin to that of orthogonal polynomials.Thus, we can reduce the problem of the (2 k + 3)-order recurrence relation to threethird-order recurrence relations, and the number of initial conditions are reduced from( k + 1) to only one for each sequence. As previously mentioned, the eigenvalue problem associated with the Hamiltonian H in (1) is equivalent to a Sturm-Liouville problem. In consequence, the eigenfunctionsshould satisfy the well-known oscillation theorem [2]. In this section, we identify thenumber of zeros of the generalized Okamoto polynomials Q m,n and the polynomials P ( k ) n ; j in each sequence of solutions through the Wronskian representation. Such zeros are dis-tributed over the support of the weight function µ k ( x ), which in this case corresponds tothe real line. Studies in this regard have been addressed by Felder [64] for Wronksians22 m,n n n + n t = n + 2 n + Q m, n n nQ m, n +1 n + m n + m ) Q m +1 , n n nQ m +1 , n +1 n n + m ) + 1Table 6: Number of zeros at x = 0 ( n ), x > n + ), and the total number of zeros in x ∈ R ( n t ) of the generalized Okamoto polynomials Q m,n . n φ ( k ) n ;1 φ ( k ) n ;2 φ ( k ) n ;3 ≤ n ≤ k n n + k + 1 3 n + k + 2 n > k n − k Table 7: Total number of zeros in x ∈ R of the eigenfunctions φ ( k ) n ;1 , φ ( k ) n ;2 , and φ ( k ) n ;3 .defined in terms of sequences of contigous Hermite polynomials, and recently in [65] forgeneral sequences. In the latter, the authors conjectured that a Wronskian of the formWr( H ξ , · · · , H ξ (cid:96) ) characterized by the partition ˜ λ = (˜ λ , · · · , ˜ λ (cid:96) ), with ˜ λ j = ξ j − j + 1, hasthe following number of real zeros [65]:(a) A zero of multiplicity n = d ξ ( d ξ +1)2 at x = 0. With d ξ = p − q , where p and q arethe number of odd and even elements, respectively, in the sequence { ξ , · · · , ξ (cid:96) } .(b) n + simple positive zeros, where n + := 12 (cid:32) (cid:96) (cid:88) j =1 ( − (cid:96) − j ˜ λ j − (cid:12)(cid:12) d ξ + ( (cid:96) − (cid:98) (cid:96) (cid:99) ) (cid:12)(cid:12) (cid:33) (92)(c) The same number of simple negative zeros, n − = n + ,where 0 ≤ ξ < · · · ξ ell and 0 ≤ ˜ λ ≤ · · · ≤ ˜ λ (cid:96) .In particular, considering the generalized Okamoto polynomials representation (71),with m, n = 0 , , · · · and m + n = 2 , , · · · , we identify (cid:96) = n +2 m − n t = n + 2 n + as shown in Tab. 6. The latter reveals that the conventionalOkamoto polynomials Q k ≡ Q k, are nodeless functions, as stated in Sec. 2.On the other hand, we can identify the zeros of the polynomials P ( k ) n ; j , which alsoare the zeros of the eigenfunctions φ ( k ) n ; j . This is achieved with the aid of the Wronskianrepresentation (85) and (92). After some calculations, we obtain the total number ofzeros of the eigenfunctions shown in Tab. 7. From the latter, it is clear that φ (0)0;1 is alwaysa nodeless function for all n and k , and thus it corresponds to the ground state of thesystem. Moreover, the zero-modes φ ( k )0;2 and φ ( k )0;3 correspond to the ( k + 1) and ( k + 2)excited states of the system, respectively. 23 Conclusions
The identification of third-order shape-invariant Hamiltonians associated with rational so-lutions of the fourth Painlev´e transcendent in the “ − x/
3” hierarchy has been developedin the most general case so that the resulting potentials are free of any singularities. Thiscondition leads to a family of potentials written in terms of the conventional Okamotopolynomials and three zero-modes constructed as the product of generalized Okamotopolynomial times a weight function with support on the real line. In this form, the re-sults of this manuscript extend and generalize some previous statements made in earlierworks [35, 47]. Interestingly, the existence of third-order ladder operators leads natu-rally to three sequences of eigenfunctions whose eigenvalues do not overlap. Although theeigenvalues of the shape-invariant Hamiltonian are not equidistant, the three independentsequences define sets of eigenvalues that are indeed equidistant within each sequence. Thelatter is a fundamental property exploited throughout the manuscript to properly identifya second-order differential equation, along with a second-order finite-difference equationfor the eigenfunctions within each sequence. The coefficients of the finite-difference equa-tion are such that only one initial condition per sequence, fixed as the correspondingzero-mode, is required to determine the higher-modes through finite iterations. This is aproperty akin to that of classical orthogonal polynomials [55, 56]. The calculations pre-sented throughout the manuscript were achieved with the aid of the new set of identitiesintroduced in (25)-(27), determined from the appropriate B¨acklund transformations.On the other hand, a direct relation between our results and that of higher-orderDarboux-Crum (higher-order SUSY) transformations [66] of the harmonic oscillator isestablished. In this form, it is found that the potentials constructed from the third-ordershape-invariant condition and the “ − x/
3” hierarchy of rational solutions are equivalentto 2 k state-deleting Darboux transformation. Such a equivalence leads simultaneously toa relation between the higher-modes and the family of exceptional Hermite polynomials,defined by the appropriate partition λ k . Remarkably, the latter provides us with a setof three different three-term recurrence relations for the exceptional Hermite polynomialslabeled by λ k , which, in general, have been identified in previous works [49, 50] throughhigher-order recurrence relations. Therefore, the separation into independent sequences ofeigenfunctions has allowed identifying simple recursion formulas for the exceptional Her-mite polynomials defined by the partition λ k that, to the best of the authors’ knowledge,have been unnoticed.It is worth noting that the existence of third-order ladder operators has been one ofthe key properties exploited throughout the text, and thus a further generalization ofladder operators of a higher order can be considered to address more general systems. Iftrue, the latter could lead to simple recurrence relations for a broader family of excep-tional Hermite polynomials defined by more general partitions. Moreover, the potentialsand solutions related to the third-order shape-invariant Hamiltonian admit general non-rational solutions that have not been exploited in detail. Those are problems that deservespecial attention by themselves, and they will be discussed elsewhere.24 cknowledgments I. Marquette was supported by Australian Research Council Future Fellowship FT180100099.V. Hussin acknowledges the support of research grants from NSERC of Canada. K. Ze-laya acknowledges the support from the Mathematical Physics Laboratory of the Centrede Recherches Mat´ematiques, through a postdoctoral fellowship. K. Zelaya also acknowl-edges the support of Consejo Nacional de Ciencia y Tecnolog´ıa (Mexico), grant numberA1-S-24569.
A B¨acklund transformations
The B¨acklund transformation is one of the most well-known techniques used to generatenew solutions by departing from an already known solution, usually called seed solution .B¨acklund transformations are given in the form of nonlinear recurrence relations. In thisform, we can iterate such a transformation so that we generate different solutions atevery step. Let w ≡ w ( x ; α , β ) be a seed solution to (16) with parameters α = α and β = β , we can then generate eight new solutions through the following B¨acklundtransformations [28] w ± := F − ∓ √− β w , w ± := − F +0 ∓ √− β w , (A-1) w ± := w + 2(1 − α ∓ √− β ) w F +0 ± √− β , w ± := w + 2(1 + α ± √− β ) w F − ∓ √− β , (A-2) F ± := w (cid:48) ± (2 xw + w ) , (A-3)where w ± i ≡ w ± i ( x ; α ± i , β ± i ) are solutions to (16) with α = α ± i and β = β ± i , for i = 1 , , , α and β through α ± = 14 (cid:16) − α ± (cid:112) − β (cid:17) , β ± = − (cid:18) α ± (cid:112) − β (cid:19) , (A-4) α ± = − (cid:16) α ± (cid:112) − β (cid:17) , β ± = − (cid:18) − α ± (cid:112) − β (cid:19) , (A-5) α ± = 32 − α ∓ (cid:112) − β , β ± = − (cid:18) − α ± (cid:112) − β (cid:19) , (A-6) α ± = − − α ± (cid:112) − β , β ± = − (cid:18) − − α ± (cid:112) − β (cid:19) . (A-7)In particular, if we use any of the rational solutions (20)-(22) corresponding to the“ − x/
3” hierarchy we obtain another rational solutions within the same hierarchy. Thisenables us to compute some new identities for the generalized Okamoto polynomials.For instance, using w ≡ w [1] m,n into (A-1)-(A-3), and comparing w +1 with w +4 − w , we25btain the identities − xQ m +1 ,n Q m,n +1 + 3 (cid:0) Q m +1 ,n Q (cid:48) m,n +1 − Q (cid:48) m +1 ,n Q m,n +1 (cid:1) = −√ Q m,n Q m +1 ,n +1 , (A-8) Q m +1 ,n +1 Q (cid:48) m,n − Q (cid:48) m +1 ,n +1 Q m,n = −√ m + 3 n + 1) Q m +1 ,n Q m,n +1 , (A-9) whereas comparing w +2 with w +3 − w leads to − xQ m,n +1 Q m,n + 3 (cid:0) Q m,n +1 Q (cid:48) m,n − Q (cid:48) m,n +1 Q m,n (cid:1) = −√ Q m +1 ,n Q m − ,n +1 , (A-10) Q m +1 ,n Q (cid:48) m − ,n +1 − Q (cid:48) m +1 ,n Q m − ,n +1 = −√ m − Q m,n +1 Q m,n . (A-11) On the other hand, after setting w ≡ w [2] m,n and comparing w +1 with w +4 − w we obtainanother set of identities of the form − xQ m,n Q m +1 ,n + 3 (cid:0) Q m,n Q (cid:48) m +1 ,n − Q (cid:48) m,n Q m +1 ,n (cid:1) = −√ Q m +1 ,n − Q m,n +1 , (A-12) Q m +1 ,n − Q (cid:48) m,n +1 − Q (cid:48) m +1 ,n − Q m,n +1 = √ n − Q m,n Q m +1 ,n . (A-13) These identities allow us to cast the rational solutions of the “ − x/
3” hierarchy into theform presented in (25)-(27).
B Action of the ladder operators on the zero-modesand higher-modes
In this section we show that the action of A † on the zero-modes and higher-modes indeedleads to a ladder operation. To this end, let us consider the zero-mode φ in terms ofthe Wronskian representation (71), together with the set indexes I := { , , · · · , k −
2; 2 , , · · · k − } , I := { , , · · · , k −
5; 2 , , · · · k − } . (B-1)In addition, it is useful to introduce a particular identity, obtained from higher-order SUSYQM [9, 40, 66]. Be H = − d dx + V ( x ) a Hamiltonian operator such that the eigenvalueequation H F ν = E ν F ν holds for ν ∈ R , along with E ν and F ν the eigenvalues andeigenfunctions, respectively. Then, the Wronskian identity [63]Wr (cid:2) Wr( F i , F i , · · · , F i n ) , Wr( F i , F i , · · · , F i n , F i n +1 , F i n +2 ) (cid:3) ∝ Wr( F i , · · · , F i n , F i n +1 ) Wr( F i , · · · , F i n , F i n +2 ) , (B-2)is fulfilled, where the eigenfunctions F i m are not necessarily square-integrable solutions.The Wronskian representation of the Okamoto polynomial Q k can be rewritten byexploiting the properties of the functions ψ n given in (73), such as the identity ddx ψ n = ψ n − . Moreover, we have ψ = 1 and the Wronskian representation of Q k in (71), whichleads to Q k ≡ C k W I = C k W { , ··· , k − , ··· , k − } = C k W { , ··· , k − , , ··· , k − } = C k W { , , ··· , k − , , ··· , k − } = C k W { , ··· , k − , ··· , k − , } = C k W I ∪{ } , (B-3)26here in the latter we have used the notation introduced in (77). Thus, the zero-mode φ ( k )0;1 takes the following form φ ( k )0;1 = N e − x Q k Q k +1 = N e − x W I ∪{ } W I . (B-4)We now proceed with the action of the creation operator A = Q † M M . First, the actionof M is determined through M φ ( k )0;1 = N (cid:20) − x ddx ln W I ∪{ } W I ∪{ k − } (cid:21) e − x W I ∪{ } W I , (B-5)which simplifies using the Wronskian ˜ W { i , ··· ,i n } , introduced in (78)-(79), resulting in M φ ( k )0;1 = N ddx ln ˜ W I ∪{ } ˜ W I ∪{ k − } e − x W I ∪{ } W I . (B-6)Then, using the identity (B-2), and after some straightforward calculations we finally get M φ ( k )0;1 = N e − x W I ∪{ k − , } W I ∪{ k − } . (B-7)Following the same steps as above, the subsequent action of M and Q † is M M φ ( k )0;1 = N e − x W I ∪{ } W I = N e − x W I ∪{ , } W I . (B-8)and φ ( k )1;1 ∝ Q † M M φ ( k )0;1 ∝ N e − x W I ∪{ } W I . (B-9)respectively. The previous procedure can be easily iterated finitely many times as neces-sary, and after n iterations we get φ ( k ) n ;1 ∝ ( A † ) n φ ( k )0;1 ∝ e − x W I ∪{ n } W I . (B-10)We proceed in a similar way to determine the higher-modes related to the other twosequences by starting from the zero-modes φ ( k )0;2 and φ ( k )0;3 . We thus obtain φ ( k ) n ;2 ∝ ( A † ) n φ ( k )0;2 ∝ e − x W I ∪{ n +3 k +1 } W I , φ ( k ) n ;3 ∝ ( A † ) n φ ( k )0;3 ∝ e − x W I ∪{ n +3 k +2 } W I , (B-11)which are the higher-modes, up to a normalization constant.27 Higher-order SUSY transformation
Let us consider the dimensionless rescaled Hamiltonian for the harmonic oscillator H ,together with the respective eigenvalue equation H ≡ − d dx + V ( x ) , V ( x ) = x − , H φ (1) n = E n φ (1) n , (C-1)where the eigenfunctions φ (1) n and eigenvalues E (1) n are respectively given by φ (1) n ( x ) = e − x (cid:112) n n ! √ π H n (cid:18) x √ (cid:19) , E (1) n = 2 n . (C-2)Now, from the iterative action of the intertwining operators B j and B † j , constructed suchthat H j := B † j B j + (cid:15) j , H j +1 := B j B † j + (cid:15) j = B † j +1 B j + (cid:15) j +1 , (C-3)with (cid:15) j a real constant, known as the factorization energy , and the Hamiltonians of theform H j := − d dx + V j ( x ) , (C-4)where the intertwining are operators defined as B j := ddx + β j ( x ) , B † j := − ddx + β j ( x ) , (C-5)and the superpotentials β j are determined from the Riccati equation − [ β j ] (cid:48) + [ β j ] = V j − (cid:15) j , [ β j ] (cid:48) + [ β j ] = V j +1 − (cid:15) j , V = x . (C-6)The new potential partners are V j +1 ( x ) = x ddx j (cid:88) k =1 β k . (C-7)Following the conventional linearization procedure of the Riccati equation, it is straight-forward to obtain the set of linear equations β ( p ) r := − ddx ln u ( p ) r , − d dx u ( p ) r + V p ( x ) u ( p ) r = (cid:15) r u ( p ) r , (C-8)with u ( p ) r the seed functions associated with the factorization energy (cid:15) r of the p -th partnepotential. The superpotentials in (C-5) correspond to the particular case β ( r ) r ≡ β r . Fromthe latter, we can write the j -th partner potential V j +1 ( x ) in term of the eigenfunctionsΦ n ≡ φ (1) n of the initial Hamiltonian h as V j +1 ( x ) = x − d dx ln W (Φ (cid:15) , Φ (cid:15) · · · Φ (cid:15) j ) , Φ (cid:15) j ( x ) ≡ φ (1) (cid:15) j ( x ) . (C-9)28 eferences [1] D. Schuch, Quantum Theory from a Nonlinear Perspective, Riccati Equations inFundamental Physics , Springer, Cham, 2018.[2] E.L. Ince,
Ordinary differential equations , Dover Publications, New York, 1956.[3] E. Schr¨odinger, A Method of Determining Quantum-Mechanical Eigenvalues andEigenfunctions,
Proc. Roy. Irish Acad. Sect. A: Math. Phys. Sci. (1940) 9.[4] E. Schr¨odinger, Further Studies on Solving Eigenvalue Problems by Factorization, Proc. Roy. Irish Acad. Sect. A: Math. Phys. Sci. (1940) 183.[5] E. Schr¨odinger, The Factorization of the Hypergeometric Equation, Proc. Roy. IrishAcad. Sect. A: Math. Phys. Sci. (1941) 53.[6] L. Infeld, and T.E. Hull, The Factorization Method, Rev. Mod. Phys. (1951) 21.[7] B. Mielnik, and O. Rosas-Ortiz, Factorization: Little or great algorithm?, J. Phys.A: Math. Gen. (2004) 10007.[8] B. Mielnik, Factorization method and new potentials with the oscillator spectrum, J. Math. Phys. (1984) 3387.[9] B. Mielnik, L.M. Nieto, and O. Rosas-Ortiz, The finite difference algorithm for higherorder supersymmetry, Phys. Lett. A (2000) 70.[10] V.B. Matveev and M.A. Salle,
Darboux transformation and solitons , Springer-Verlag,Berlin, 1991.[11] M.G. Darboux, Sur une proposition relative aux ´equations lin´eaires,
Comptes RendusAcad. Sci. (1882) 1456.[12] E. Witten, Dynamical breaking of supersymmetry, Nucl. Phys. B (1981) 513.[13] F. Cooper, A. Khare, and U. Sukhatme,
Supersymmetry in Quantum Mechanics ,World Scientific, Singapore, 2001.[14] G. Junker,
Supersymmetric Methods in Quantum and Statistical Physics , Springer,Berlin , 1996.[15] R. Milson, On the Liouville transformation and exactly-solvable Schr¨odinger equa-tions,
Int. J. Theor. Phys. Ann. HenriPoincar´e (2011) 397.[17] O. Rosas-Ortiz, O. Casta˜nos, and D. Schuch, New supersymmetry-generated complexpotentials with real spectra, J. Phys. A: Math. Theor. (2015) 445302.2918] F. Cannata, G. Junker, and J. Trost, Schr¨odinger operators with complex potentialbut real spectrum, Phys. Lett. A (1998) 219.[19] Z. Blanco-Garc´ıa, O. Rosas-Ortiz, and K. Zelaya, Interplay between Riccati, Ermakovand Schr¨odinger equations to produce complex-valued potentials with real energyspectrum,
Math. Meth. Appl. Sci (2019) 4925.[20] K. Zelaya, S. Cruz y Cruz and O. Rosas-Ortiz, On the construction of non-HermitianHamiltonians with all-real spectra through supersymmetric algorithms, In: Geomet-ric Methods in Physics XXXVIII. Trends in Mathematics; P. Kielanowski, A. Odzi-jewicz, E. Previato (Eds.) , pages 283-292, Birkh¨auser, Basel, 2020.[21] A. Contreras and V. Jakubsk´y, Photonic systems with two-dimensional landscapes ofcomplex refractive index via time-dependent supersymmetry,
Phys. Rev. A (2019)053812.[22] M. Znojil, F. Cannata, B. Bagchi, and R. Roychoudhury, Supersymmetry withouthermiticity within PT symmetric quantum mechanics, Phys. Lett. B (2000) 284.[23] B. Bagchi, S. Mallik, and C. Quesne, Generating complex potentials with real eigen-values in supersymmetric quantum mechanics,
Int. J. Mod. Phys. A (2001) 2859.[24] F. Correa, V. Jakubsk`y, and M. S. Plyushchay, PT-symmetric invisible defects andconfluent Darboux-Crum transformations, Phys. Rev. A (2015) 023839.[25] V.I. Gromak, I. Laine, and S. Shimomura, Painlev´e Differential Equations in theComplex Plane , Walter de Gruyter, Berlin, 2002.[26] F. Marcell´an, and W. Van Assche (Eds.),
Orthogonal Polynomials and Special Func-tions: Computation and Applications , Springer-Verlag, Berlin, 2006.[27] A.P. Bassom, P.A. Clarkson, and A.C. Hicks, B¨acklund transformations and thesolution hierarchies for the fourth Painlev´e equation,
Stud. Appl. Math. (1995) 1.[28] P.A. Clarkson, The fourth Painlev´e transcendent, In Differential Algebra and RelatedTopics II , Li Guo and W.Y. Sit (eds.), World Scientific, Singapore, 2008.[29] P.A. Clarkson, The fourth Painlev´e equation and associated special polynomials,
J.Math. Phys. (2003) 5350.[30] M. Florja´nczyk, and L. Gagnon, Exact solutions for a higher-order nonlinearSchr¨odinger equation, Phys. Rev. A (1990) 4478.[31] A.S. Fokas, A.R. Its, and A.V. Kitaev, Discrete Painlev´e equations and their appear-ance in quantum gravity, Commun. Math. Phys. (1991) 313.[32] D. Berm´udez and D. J. Fern´andez C, Supersymmetric quantum mechanics andPainlev´e IV equation,
SIGMA (2011) 025.3033] D. Berm´udez and D.J. Fern´andez C., Solution Hierarchies for the Painlev´e IV Equa-tion, In: Geometric Methods in Physics. Trends in Mathematics , P. Kielanowski, S.Ali, A. Odzijewicz, M. Schlichenmaier, T. Voronov (eds); Birkh¨auser, Basel, 2013.[34] D. Bermudez, D.F. Fern´andez, and J. Negro, Solutions to the Painlev´e V equationthrough supersymmetric quantum mechanics,
J. Phys. A: Math. Theor. (2016)335203.[35] I. Marquette, Superintegrability with third order integrals of motion, cubic algebras,and supersymmetric quantum mechanics. I. Rational function potentials, J. Math.Phys. (2009) 012101.[36] I. Marquette, Superintegrability with third order integrals of motion, cubic algebras,and supersymmetric quantum mechanics. II. Painlev´e transcendent potentials, J.Math. Phys. (2009) 095202.[37] L. ´E. Gendenshteˆın, Derivation of exact spectra of the Schr¨odinger equation by meansof supersymmetry, JETP Lett. (1983) 356.[38] A. Andrianov, F. Cannata, M. Ioffe, and D. Nishnianidze, Systems with higher-ordershape invariance: spectral and algebraic properties, Phys. Lett. A (2000) 341.[39] U.P. Sukhatme, C. Rasinariu, and A. Khare, Cyclic shape invariant potentials,
Phys.Lett. A (1997) 401.[40] D.J. Fern´andez, and V. Hussin, Higher-order SUSY, linearized nonlinear Heisenbergalgebras and coherent states,
J. Phys. A: Math. Gen. (1999) 3603.[41] D.J. Fern´andez, V. Hussin, and O. Rosas-Ortiz, Coherent states for Hamiltoniansgenerated by supersymmetry, J. Phys. A: Math. Gen. (2007) 6491.[42] O. Rosas-Ortiz, and K. Zelaya, Bi-Orthogonal Approach to Non-Hermitian Hamil-tonians with the Oscillator Spectrum: Generalized Coherent States for NonlinearAlgebras, Ann. Phys. (2018) 26.[43] A.B.J. Kuijlaars and R. Milson, Zeros of exceptional Hermite polynomials,
J. Approx.Theor. (2015) 28.[44] D. Gomez-Ullate, Y. Grandati, and R. Milson, Spectral Theory of Exceptional Her-mite Polynomials, arXiv:2012.02354 [math.CA].[45] V. Chalifour and A.M. Grundland, General solution of the exceptional Hermite dif-ferential equation and its minimal surface representation,
Ann. Henri Poincar´e (2020) 3341.[46] I. Marquette, and C. Quesne, Connection between quantum systems involving thefourth Painlev´e transcendent and k-step rational extensions of the harmonic oscillatorrelated to Hermite exceptional orthogonal polynomial, J. Math. Phys. (2016)052101. 3147] K. Zelaya, I. Marquette, and V. Hussin, Fourth Painlev´e and Ermakov equations:quantum invariants and new exactly-solvable time-dependent Hamiltonians, J. Phys.A: Math. Theor. (2021) 015206.[48] D. G´omez-Ullate, Y. Grandati, and R. Milson, Rational extensions of the quantumharmonic oscillator and exceptional Hermite polynomials, J. Phys. A: Math. Theor. (2014) 015203.[49] D. G´omez-Ullate, A. Kasman, A.B.J. Kuijlaars, and R. Milson, Recurrence relationsfor exceptional Hermite polynomials, J. Approx. Theor. (2016) 1.[50] S. Odake, Recurrence relations of the multi-indexed orthogonal polynomials,
J. Math.Phys. (2013) 08350.[51] A.A. Andrianov, M.V. Ioffe, F. Cannata, J.-P. Dedonder, Second order derivativesupersymmetry, q deformations and the scattering problem, Int. J. Mod. Phys. A (1995) 2863.[52] A. Andrianov, F. Cannata, M. Ioffe, D. Nishnianidze, Systems with higher-ordershape invariance: spectral and algebraic properties, Phys. Lett. A (2000) 341.[53] F.W.J. Olver, et al. (eds.),
NIST Handbook of Mathematical Functions , CambridgeUniversity Press, New York, 2010.[54] K. Okamoto, Studies on the Painlev´e equations III. Second and fourth Painlev´eequations, P II and P IV , Math. Ann. (1986) 221.[55] T. S. Chihara,
An Introduction to Orthogonal Polynomials , Dover Publications Inc.,New York, 1978.[56] G. Szeg¨o,
Orthogonal Polynomials , American Mathematical Society, New York, 1959.[57] J. Favard, Sur les polynomes de Tchebicheff,
Comptes Rendus de l’Acad´emie desSciences , (1935) 2052.[58] A.A. Oblomkov, Monodromy-free Schr¨odinger operators with quadratic increasingpotentials, Theor. Math. Phys. (1999) 1574.[59] M. Noumi and Y. Yamada, Symmetries in the fourth Painlev´e equation and Okamotopolynomials,
Nagoya Math. J. (1999) 53.[60] K. Kajiwara and Y. Ohta, Determinant structure of the rational solutions for thePainlev´e IV equation,
J. Phys. A: Math. Gen. (1998) 2431.[61] P.A. Clarkson, Special polynomials associated with rational solutions of the defocus-ing nonlinear Schr¨odinger equation and the fourth Painlev´e equation, Euro. J. App.Math. (2006) 293. 3262] D. G´omez-Ullate, Y. Grandati, and R. Milson, Durfee Rectangles and Pseudo-Wronskian Equivalences for Hermite Polynomials, Stud. App. Math. (2018) 596.[63] S. Odake, Wronskian/Casoratian identities and their application to quantum me-chanical systems,
J. Phys. A: Math. Theor. (2020) 365202.[64] G. Felder, A.D. Hemery, and A.P. Veselov, Zeros of Wronskians of Hermite polyno-mials and Young diagrams, Physica D (2012) 2131.[65] M. ´A. Garc´ıa-Ferrero and D. G´omez-Ullate, Oscillation Theorems for the Wronskianof an Arbitrary Sequence of Eigenfunctions of Schr”odinger’s Equation,
Lett. Math.Phys. (2015) 551.[66] M.M. Crum, Associated Sturm-Liouville systems,
Q. J. Math.6