In this work, we present a rigorous mathematical scheme for the derivation of the *s*th-order perturbative corrections to the solution to a one-dimensional harmonic oscillator perturbed by the potential *V*_{per}(*x*) = *λx*^{α}, where *α* is a positive integer, using the non-degenerate time-independent perturbation theory. To do so, we derive a generalized formula for the integral $I=\u222b\u2212\u221e+\u221ex\alpha \u2061exp(\u2212x2)Hn(x)Hm(x)dx$, where *H*_{n}(*x*) denotes the Hermite polynomial of degree *n*, using the generating function of orthogonal polynomials. Finally, the analytical results with *α* = 3 and *α* = 4 are discussed in detail and compared with the numerical calculations obtained by the Lagrange-mesh method.

## I. INTRODUCTION

Approximation methods play a crucial role in quantum mechanics since the number of problems that are exactly solvable is small in comparison to those that must be solved approximately. To our knowledge, the hydrogen atom, harmonic oscillators, and quantum particles in some specific potential wells have exact solutions,^{1–4} and two cold atoms interacting through a point-like force in a three-dimensional harmonic oscillator potential^{5} can also be solved analytically. Recently, Jafarov *et al.* reported an exact solution to the position-dependent effective mass harmonic oscillator model.^{6} Because of this, approximation methods have been developed early since the dawn of quantum mechanics. One of the essential approximation methods is the perturbation theory (PT) established by Schrödinger in 1926.^{7} Later, it was immediately used to interpret the LoSurdo–Stark effect of the hydrogen atom by Epstein.^{8} Although the PT is not well convergent at higher-order corrections,^{1–3,9,10} many more efficient approximation methods have been developed to treat quantum-mechanically complex problems,^{11} it is still a paramount and elementary approximation method in quantum physics. For instance, the PT contributes significantly to quantum optics and quantum field theory, as discussed in Ref. 12. For this reason, the PT is usually presented clearly and discussed in detail in quantum mechanics textbooks for undergraduate students.^{1–3,9}

The one-dimensional harmonic oscillator is not only a rich pedagogical example for approximation theories in quantum mechanics^{13–24} but also an excellent candidate for numerical^{25,26} and analytical or algebraic methods,^{6,27–30} owing to its simple calculation and exact solution. Moreover, the one-dimensional harmonic oscillator potential has also played a key role in studies of ultra-cold atomic quantum gases for the last two decades.^{5,31–51} Moshinsky and Smirnov^{52} provided a deep review of the role of quantum harmonic oscillators in modern physics. In standard quantum mechanics textbooks,^{1–3} a one-dimensional anharmonic oscillator is usually presented as an application of PT because the solutions can be constructed based on the exact solution to a one-dimensional harmonic oscillator. The two perturbative potentials that are usually considered are *V*_{per}(*x*) = *λx* and *V*_{per}(*x*) = *λx*^{4}. It is common to present the calculation of the first- and second-order corrections of the energy; however, the corrections to the wave function are usually not provided in detail.^{1–3} Moreover, physicists have also studied the generalized case *V*_{per}(*x*) = *λx*^{2β}, where *β* is a positive integer. The case *β* = 2 has been studied using the WKB method,^{53} the intermediate Hamiltonian,^{54} the Padé approximation,^{55–57} the Heisenberg matrix mechanics,^{58} and the variational perturbation.^{59} In addition, the cases *β* = 3 and *β* = 4 were studied in Refs. 60 and 61, and the considered potential of the harmonic oscillator is *V*(*x*) = *x*^{2} instead of *V*(*x*) = 0.5*x*^{2}. The authors intended to establish efficient methods to solve the problem mathematically, regardless of their physical meaning. Recently, the problem has been extended to sextic (*x*^{6}) and decatic (*x*^{10}) potentials using polynomial solutions^{62,63} and a polynomial perturbative potential.^{64} Interestingly, we realized that in the above-mentioned works, the authors only considered even values of the power of *x*, while the odd cases have not been studied. It is also interesting to note that the wave function was not considered in the above-mentioned articles. Essentially, the applications of the one-dimensional anisotropic oscillator can be found in chemistry, in which the perturbative potential is used to study the vibration in molecules.^{65–69} In addition, the perturbative potential *λx*^{4} has recently been used to model the Brownian motion of particles in optical tweezers.^{70} Consequently, it is necessary to compute the approximated wave function and the energy of a one-dimensional harmonic oscillator perturbed by the potential *V*_{per}(*x*) = *λx*^{α} for arbitrary eigenstates with arbitrary values of *α*.

The goal of this study is to present a systematic and complete treatment of the *s*th-order perturbative corrections to the solution to a one-dimensional harmonic oscillator perturbed by the potential *V*_{per}(*x*) = *λx*^{α}. To achieve this goal, we derived a formula for $I=\u222b\u2212\u221e+\u221ex\alpha \u2061exp(\u2212x2)Hn(x)Hm(x)dx$, with *H*_{n}(*x*) being the Hermite polynomial of degree *n*. Our scheme is based on the so-called generating functions of orthogonal polynomials.^{71} Because the potential depends solely on the spatial coordinate and the states are non-degenerate, the non-degenerate time-independent PT is used to derive the corrections to the wave function and energy. Note that our results can be used for arbitrary eigenstates of a one-dimensional anharmonic oscillator with an arbitrary power coefficient *α*. This is significantly different from previous works, as discussed above.

The remainder of this paper is organized as follows: Sec. II briefly outlines the time-independent PT for non-degenerate states. Section III presents the main results and discussion. Finally, conclusions are presented in Sec. IV. For simplicity, we use atomic units in which *ℏ* = *m* = *ω* = 1 throughout this study. In addition, the notation $Xn,s\alpha $ denotes the *s*th-order perturbation correction to the physical quantity *X* in the state with the quantum number *n* and the power coefficient *α*.

## II. NON-DEGENERATE TIME-INDEPENDENT PERTURBATION THEORY AND THE 2*s* + 1 RULE

This section presents the time-independent PT for non-degenerate states and a recurrence relation to obtain higher-order corrections to the wave function and energy. We followed the procedure of Fernandez.^{9} The Schrödinger equation describing a one-dimensional quantum system is as follows:

where $H\u0302$ is the Hamiltonian operator and *E*_{n} is the eigenvalue corresponding to the eigenfunction *ψ*_{n}. The Hamiltonian can be split into two parts,

with $H\u03020$ being the Hamiltonian operator, whose eigenvalues and eigenfunctions are analytically solvable and satisfy the equation

and $H\u0302\u2032$ being sufficiently small and considered as a small perturbation with parameter *λ*. The Taylor formula is used to expand the energy and the eigenfunction of the Hamiltonian $H\u0302$ as a function of perturbation parameter *λ*,

where *s* is the order of the perturbative correction. Substituting Eq. (4) into Eq. (1), we obtain a recurrence equation expressing the relation between the corrections to the eigenfunction *ψ*_{n,s} and the energy *E*_{n,s} as follows:

The perturbation correction to the wave function, *ψ*_{n,s}, is then expanded as a linear combination of the eigenfunctions of the non-perturbative Hamiltonian, $H\u03020$, as follows:

where *c*_{mn,s} = ⟨*ψ*_{m,0}|*ψ*_{n,s}⟩ is the expanding coefficient. Substituting Eq. (6) into Eq. (5) and then integrating over the whole space after multiplying both sides by $\psi m,0*$ yield

where $H\u0302mk\u2032=\u27e8\psi m,0|H\u0302\u2032|\psi k,0\u27e9$. For non-degenerate states, *E*_{m,0} ≠ *E*_{n,0} if and only if *m* ≠ *n*, we can derive the general expression for the correction to the energy by letting *m* = *n* in Eq. (7); hence, we obtain

The normalization condition ⟨*ψ*_{n}|*ψ*_{n}⟩ = 1 results in a constraint on the correction to the wave function,

Finally, in the case of *m* ≠ *n*, we can derive the expanding coefficients for the correction to the wave function as follows:

and

for the case of *m* = *n*. Substituting *s* = 1 into Eq. (8) derives the formula for the first-order correction to the energy, which is the average of the perturbation potential with respect to the eigenfunction *ψ*_{n,0},

and the expanding coefficient for the first-order correction to the wave function is then derived as follows:

Obtaining the second-order correction to the energy is also straightforward. It is obtained as follows:

These results can be found in standard quantum mechanics textbooks.^{1–3} To derive higher-order corrections to the energy, it is obvious that one can use the recurrence given by Eq. (8). However, there is another way to quickly compute the corrections to the energy. It is called the 2*s* + 1 rule, in which, once we know the *s* order of the correction to the wave function, we are allowed to compute the corrections to the energy up to the 2*s* + 1 order. For the detailed derivation of the rule, one should refer to the textbook.^{9} Below, we list the formula for the third-, fourth-, and fifth-order corrections to the energy used for calculations in Sec. III,

## III. RESULTS AND DISCUSSION

### A. Derivation of the *s*th-order perturbative corrections to the solution to a one-dimensional anharmonic oscillator

The Schrödinger equation describing a one-dimensional harmonic oscillator induced by a perturbation potential *λx*^{α} is

where *λ* is the strength of the external field, which gives rise to the perturbation, and *α* is a positive integer. In the absence of the perturbation, Eq. (18) is the well-known equation describing a one-dimensional harmonic oscillator with a wave function

where $An=12nn!\pi $ is the normalization constant, *n* is the quantum number, and *H*_{n}(*x*) is the Hermite polynomial of degree *n*, and the energy is given by

Since Eq. (18) cannot be analytically solvable, the PT is then chosen to approximate the solutions. As discussed above, the first-order correction of the wave function is given by the following equation:

where

is the expanding coefficient of the first-order correction to the wave function. Making use of Eq. (A12) (see the Appendix), we obtain the following equation:

where *δ*_{m,n} denotes the Kronecker delta, satisfying

By substituting Eq. (23) into Eq. (21), the first-order correction to the wave function for arbitrary states is obtained by

The first-order correction to the energy is then computed by the following equation:

Combining the known wave function of a one-dimensional harmonic oscillator and Eq. (A12), the first-order correction to the energy is obtained by

It is interesting to note that *E*_{n,1} is non-zero only if

owing to the mathematical property of the Kronecker delta. Because *k* and *ℓ* are integers, the above equation has the solutions for even *α* solely. This means that in the case of odd *α*, the first-order correction to energy always equals to zero. Therefore, the anharmonic oscillator does not feel the presence of the external field in this case if only the first-order approximation is considered. Therefore, it is necessary to compute higher-order corrections to the energy. Regarding the second-order correction of energy, it can be computed as follows:

Once again, the integral can be treated by making use of Eq. (A12),

Substituting back into Eq. (30), the general second-order correction to the energy is obtained as follows:

Using Eq. (15), the third-order correction to the energy for arbitrary states could be derived as follows:

The wave function of a one-dimensional anharmonic oscillator with first-order correction is given by the following equation:

and the corresponding energy with third-order correction is as follows:

Obviously, the above results can be used to approximate the wave function and the energy for arbitrary states and arbitrary power *α*. Since the calculation for the higher-order corrections is more tedious for the general case *α*, we constrained our calculation for the general power *α* at this point. In the following, we discuss two particular circumstances in which *α* = 4 and *α* = 3 in detail and then extend the calculation to the second-order correction to the wave function and the fourth- and fifth-order corrections to the energy for these cases.

### B. Typical cases: *α* = 4 and *α* = 3

In this section, we discuss two particular circumstances where *α* = 4 and *α* = 3 as two typical examples of PT. In addition, to validate the significance of the higher-order corrections to the energy obtained in the textbook,^{9} we calculated the second-order correction to the wave function for these two cases and then computed the energy up to the fifth-order approximation. The numerical results obtained by the Lagrange-mesh method^{25,26} were then used as the benchmark to validate the applicable range of the analytically approximated results.

First, let us discuss the case where *α* = 4 explicitly. According to Eq. (25), the first-order correction to the wave function is given by the following equation:

where (*a*)_{n} = *a*(*a* + 1)⋯(*a* + *n* − 1) is the Pochhammer symbol. The first-, second-, and third-order corrections to the energy are obtained by using Eqs. (27), (32), and (33), respectively, with *α* = 4,

The second-order correction can be derived by the 2*s* + 1 rule, which was presented in Sec. II. The calculation shows that

Consequently, the fourth- and fifth-order corrections to the energy are obtained, respectively, as follows:

Next, we compare the analytically approximated formulas with numerical calculations to validate the exactness and determine the applicable range of the approximated formulas. For this purpose, the relative deviation is given by the following equation:

where *E*_{num} and *E*_{ana} are the numerical and analytical results, respectively. As shown in Fig. 1, in the small *λ* regime (0 ≤ *λ* ≤ 0.05), the approximated formulas match well to the numerical calculation. In this regime, the tenth-order corrected formula has the lowest relative deviation, approximately zero. However, the higher-order formulas diverge rapidly as the perturbative parameter increases. Nevertheless, in the large *λ* regime, the first-order corrected formula has, in general, a smaller relative deviation compared to others.

Finally, we explicitly present the corrections to the energy and wave function in the case of *α* = 3. The first- and second-order corrections to the wave function of this case are given, respectively, as follows:

and

Straightforwardly, we obtain

Similarly, the analytically approximated formulas were compared to the numerical calculation. Because the odd-order corrections to the energy are zero, the second- and fourth-order corrections are considered. The results are depicted in Fig. 2. The comparison shows that the second- and fourth-order corrected energies are relatively identical to the numerical results. The deviations of the ground and first-excited state are less than 1%, while those of higher-excited states are larger, ∼5%, which is still acceptable.

## IV. CONCLUSION

In summary, we have provided a mathematical procedure to derive the wave function and energy for any arbitrary states of a one-dimensional anharmonic oscillator with the general perturbation potential *V*_{per}(*x*) = *λx*^{α} using the time-independent non-degenerate PT. Subsequently, the explicit results of two particular cases in which *α* = 3 and *α* = 4 are discussed in detail. In addition, the second-order correction to the wave function and up to fifth-order correction to the energy for these cases were computed. The analytical results were then compared with the numerical solutions obtained using the Lagrange-mesh method. The comparison indicates that in the regime in which the perturbation parameter is small, the analytical results agree well with those obtained numerically and then dramatically diverge as the perturbation parameter increases. The higher-order corrections to the energy diverge rapidly as the perturbative parameter increases. In addition, the fifth-order corrected energy in the case of *α* = 4 is the same as that printed in Ref. 9; however, the procedure is more comfortable to approach. The results in the present article are anticipated to be a useful reference.

## ACKNOWLEDGMENTS

This work was supported by the Ministry of Education and Training of Vietnam (Grant No. B2021-SPS-02-VL). Mr. Tran Duong Anh-Tai acknowledges support provided by the Okinawa Institute of Science and Technology Graduate University (OIST). The authors are thankful to Mr. Mathias Mikkelsen for introducing the Lagrange-mesh method to them.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: THE INTEGRAL $I=\u222b\u2212\u221e+\u221ex\alpha \u2061exp(\u2212x2)Hn(x)Hm(x)dx$

In this work, the Hermite polynomials of degrees *m* and *n* are expanded by generating function,^{71} respectively,

Then, multiplying both sides of the above two equations by *x*^{α} exp(−*x*^{2}) and taking the integral from −∞ to +∞, we obtain

To simplify the calculation, let us introduce a new variable *y* = *x* − (*t* + *z*), and then the right-hand side is rewritten as

The binomial in (A4) is expanded by the binomial formula

and hence,

The integral in Eq. (A6) can be straightforwardly deduced as

therefore solely the integrals with even coefficients *i* = 2*k* are considered

It can be seen in (A3) that the value of the integral is the coefficient of *t*^{n}*z*^{m}; thus, we need to find that expanding coefficient. To do so, we expand

by the Taylor formula and

by the binomial formula. Substituting into (A4), we obtain

## REFERENCES

^{M}perturbation

^{2m}anharmonic oscillators