Quantum Harmonic Oscillator

Quantum harmonic oscillator (QHO) involves square law potential (x 2 ) in the Schrodinger equation and is a fundamental problem in quantum mechanics. It can be solved by various conventional methods such as (i) analytical methods where Hermite polynomials are involved, (ii) algebraic methods where ladder operators are involved, and (iii) approximation methods where perturbation, variational, semiclassical, etc. techniques are involved. Here we present the general outcomes of the two conventional semiclassical approximation methods: the JWKB method (named after Jeffreys, Wentzel, Kramers, and Brillouin) and the MAF method (abbreviated for “ modified Airy functions ” ) to solve the QHO in a very good precision. Although JWKB is an approximation method, it interestingly gives the exact solution for the QHO except for the classical turning points (CTPs) where it diverges as typical to the JWKB. As the MAF method, it enables very approximate wave functions to be written in terms of Airy functions without any discontinuity in the entire domain, though, it needs careful treatment since Airy functions exhibit too much oscillatory behavior. Here, we make use of the parity conditions of the QHO to find the exact JWKB and approximate MAF solutions of the QHO within the capability of these methods.

JWKB method is known to give exact eigenenergies for the QHO, but eigenfunctions fail at and around the classical turning points (CTPs) where f ¼ 0 (or, equivalently, E n ¼ U r ð Þ) in (3) as typical to the JWKB method [1][2][3][4][5][6][7]14]. These discontinuities prevent us from using continuity at the boundaries by equating the JWKB solutions of two neighboring regions directly at the CTPs to find the eigenenergy-dependent coefficients in the general JWKB eigenfunctions (wave functions). It also prohibits the use of normalizability of the eigenfunctions between À∞ and ∞. To surmount the problem, parity conditions of the problem regarding the symmetry of the QHO in the dimensionless form are used, and advanced computational software such as Mathematica can be used to achieve these calculations [3,4,14,24]. Moreover, asymptotic matching is required in the JWKB solutions to maintain the normalizability except for the CTPs as discussed above [4,7,14]. As to the MAF method, it does not exhibit discontinuities at the CTPs, though highly oscillating behavior of the Airy functions requires careful handling in finding their zeros and the parity treatment used in the JWKB solution seems straightforward to be also applicable to the MAF solution of the QHO [3,[19][20][21][22]. Although it was originally suggested in 1931 by Langer in [25], finding zeros of highly oscillatory Airy functions became practical as the advances in computational software and the MAF method became widespread by the 1990s [3,[18][19][20][21][22][23]. In this work, we present the general outcomes of the conventional JWKB and MAF methods as two semiclassical conventional methods and solve the QHO by using the parity condition of the problem in the dimensionless form pedagogically. We will discuss the treatment of parity matching and asymptotic matching in solving the QHO by these semiclassical methods.

Exact solution of the QHO in 1D by the analytic method
The QHO in (3) is a bound-state problem which can be written in 1D for the potential function in 1D (U or, simply, whose solution by various conventional approaches (such as analytical, algebraic, approximation, etc.) is given in any fundamental textbooks, that is, [1][2][3]5] and whose results can be summarized as follows [14]: i. Change of variable in (4) and (5): ii. TISE for the QHO in 1D in dimensionless form: Note that here f ≕ k 2 is a function of λ n E n ð Þ&y and λ n ≥ 0 since E n ≥ 0. Moreover, f λ n ; y ð Þis an even function as shown in Figure 1 (λ is chosen as continuous including the discrete energy values assuming that the eigenenergies have not been found yet).
iii. Exact eigenenergies: iv. Exact eigenfunctions (wave functions) in y: v. By using (6), we have the wave functions in x: We used two different symbols (φ and ψ) to label the functions in two different independent variables (in x and in y, respectively). Exact wave functions in (9) (via exact eigen energies in (8)) are given for even and odd n values in Figures 3 and 4 along with the JWKB solutions for comparison. H n in (9) and (10) is Hermite polynomials with indice n (named after the French mathematician Charles Figure 1.
Rodriguez formula: H n x ð Þ ¼ À1 ð Þ n e x 2 ð Þ d n dx n e Àx 2 ð Þ Generating function: Hn x ð Þt n n! Some of the Hermite polynomials: Recurrence relations: Hermite). Some of the properties of Hermite polynomials are tabulated in Table 1, and calculation of conversion factor β which exhibits a quantization, namely, is given along with the related Mathematica codes in [14].

A review of the JWKB solution of the QHO
2D plot of Figure 1 is schematically given in Figure 2 for the QHO under study (in the dimensionless form) from which we have the following outcomes [14]:

JWKB solution of eigenfunctions (wave functions) of the QHO
Conventional first-order JWKB solution of the QHO given in the normal form in (4) or (7) is as follows: where y t is either of the classical turning points (CTPs: either "y 1 , on the left" or "y 2 , on the right" depending on the region under question) and c J1 &c J2 are arbitrary JWKB constants. Once solution in any region is found (say, φ JII ), solution in the adjacent region (say, φ JIII ) can directly be found by the conventional JWKB connection formulas given in [3,4,14] without calculating it via (15). The integrals here are the definite integrals whose upper and lower values should be chosen as the related turning point (either y 1 or y 2 ) and the variable y should be in the correct ascending integration order. Normally, constant coefficients in the general solutions are determined from normalization by applying the boundary conditions of the eigenvalue problem. However, it is useless since the boundary conditions correspond to the CTPs at which (and also in a narrow region around) the conventional first-order JWKB solutions typically diverge [3,4,7,14]. This might be thought as a violation of continuity requirement of the acceptable wave function properties concerning continuity, but higher order JWKB approximation can fix it. Now, due to the discontinuities at the boundaries between the adjacent regions (such as between I and II or between II and III), the unidirectional JWKB connection formulas (surely, for the first-order JWKB) given in the literature [3,4,7,14] cannot be used to find the constant coefficients in the general solution in (18) and (19). Note that these connection formulas can be used to determine the structure of the JWKB solutions in all regions (I, II, and III), but they cannot be used to find the constant coefficients (which will be a function of eigenenergy) as explained. However, we are fortunately not helpless: (7) is an even function (see Figure 2), we should have even and odd-parity solutions. If we start by considering the exact solutions in (9) and (10) and considering them to be approximate to the JWKB solution (shown with tilde and subscript J), we have the following outcomes [14]: where p&q are positive real constants regarding the even-parity (EP) and oddparity (OP) initial values of the physical system. Remember that we use φ for the xsystem and ψ for the y-system as shown in (16). In finding the constant coefficients, we can take AEp ¼ AEq ¼ 1, and alternating sign can be modified as a parity matching as follows [14]:  (19) where the superscripts (par.m.) and (asy.m.) represent parity matched and asymptotically matched JWKB solutions, respectively. Eqs. (18) and (19) tells that we will take AEp ¼ AEq ! 1 to find the solutions in 0 ≤ y ≤ y firstly by using (21) for the asymptotic matching and then extending it to the second quadrant according to the parity under question. Note that asymptotically matched general JWKB ð Þ 1 solution can be obtained as follows (see [3,4,7,14] for details): φ asy:m: ð Þ J λ n ; y ð Þ ¼φ asy:m: ð Þ JI λ n ; y ð Þ ¼ eitherk 1φ1 λ n ; y , y t1 À Á ork 2φ2 λ n ; y , y t1 À Á φ JII λ n ; y ð Þ¼φ J λ n ; y t1 , y , y t2 À Á φ asy:m: ð Þ JIII λ n ; y ð Þ ¼ eitherk 1φJ1 λ n ; y t2 , y À Á ork 2φJ2 λ n ; y t2 , y À Á 8 > > < > > : (20) so that they exhibit the following asymptotic behaviors:

Even-parity (EP) wave functions
When initial values at x ¼ y ¼ 0 for the EP case in (17), namely (by using (16)), is applied to the JWKB solution in (15), we find the following: where the second complementary solution (in the sine form) has been canceled and calculation of the integral in the cosine term can be calculated by the similar change of variable as in (13) whose result will give η y; 0 ð Þ(see Eq. (18) below and apply η y; 0 ð Þ ¼ η λ n ! y; y ! 0 ð Þ ).
iiÞψ J, II λ n ; y and by using (16) Now, by applying the JWKB connection formula with a small phase term α, we get.
, for λ n , y , ∞, and the asymptotically matched (modified) wave function in region III via (20) and (21) of [3,4,7,14] gives: ψ asy:m: ð Þ IIIλ n ¼ λ n ; y Abbreviations we use for the EP JWKB solutions here (and also for the OP solutions in the next subsection) are as follows [14]: Since we have already calculatedψ J, II β; λ n ; x ð Þandψ asy:m: ð Þ IIIλn ¼ λ n ; yÞ À in the first quadrant (0 ≤ y ≤ λ n ), JWKB solutions in the other regions can be easily written as in (18). JWKB wave functions regarding the EP case are given in Figure 3 along with the exact solutions for comparison.
Again, since we have already obtainedψ J, II λ n ; y ð Þandψ m: ð Þ J, III λ n ; y ð Þin 0 ≤ y ≤ λ n , JWKB solutions in the second quadrant can be written in terms of them as shown in (19). JWKB wave functions regarding the OP case are given in Figure 3 along with the exact solutions for comparison.

General structure of the MAF approximation to the bound-state wave functions
Formal MAF method suggests a solution to the TISE in (7) in terms of Airy functions as follows: where Ai and Bi represent the Airy functions (namely, Ai x ð Þ and Bi x ð Þ are the linearly independent solutions of the Airy differential equation y ″ x ð Þ À xy x ð Þ ¼ 0 in x), a 1 &a 2 are the arbitrary constants which will be found from boundary values, and F&G are the functions to be determined. Note that the first variable λ n is the eigenenergies (constant values quantized by index n) which will also be determined soon. So, for now, we can consider all these functions as one dimensional in only y for simplification (say,ψ J λ n ; y ð Þ≕ψ J y ð Þ, ξ λ n ; y ð Þ≕ ξ y ð Þ, F λ n ; y ð Þ≕ F y ð Þ, etc.). If we choose one of the linearly independent solutions, say F y ð Þ:Ai ξ y ð Þ ½ , to substitute in the TISE in (4), then it gives: Now, with the choice of the last term in (35) as zero, we find the following: Here, the property of the Airy functions, Ai ″ ξ ð Þ ¼ ξAi ξ ð Þ, was used [3,18]. The integral interval in (36) is also chosen tactically in a fashion that it invokes a relationship with the turning point y t (representing the correct order y t1 or y t2 to give a non-imaginary result), and it can be written in a more explicit and conventional form (by also using in our two-variable form here) as follows: where f λ n ; y ≥ 0 ð Þ¼k 2 λ n ; y ð Þ ¼ Àκ 2 λ n ; y ð Þand y t1 &y t2 are the CTPs at the interface of the regions I À II&II À III, respectively. The remaining terms in (44) and (45) are also made zero as follows: Starting from the second term, we have.
where b 1 is some constant, and finally, making the first term in (35) zero (which is the only assumption in the MAF method), we have the following: Or more correctly in two-variable form in our eigenvalue system.
can be thought as a measure of the accuracy of the MAF solution, namely, P λ n ; y ð Þ !0, as MAF solution gets more accurate [18]. The same results would also be obtained if we had chosen the other linearly independent solution, G y ð Þ:Bi y ð Þ, in (34). Consequently, using the results found in (37) and (38), the general solution suggested in (34) can be written explicitly in the standard form of the MAF formula as follows: or more correctly in two variables here in our study.
where c 1 ¼ a 1 :b 1 and c 2 ¼ a 2 :b 2 are the arbitrary constants to be determined from the boundary values as mentioned and ∂ y ξ λ; y ð Þ represents the first derivative of ξ with respect to y. Using the result in (38), the approximation term P λ n ; y ð Þin (40) can be rewritten explicitly as follows:

MAF solution of eigenenergies
For a symmetrical f as in Figure 2, we have even-parity (EP) and odd-parity (OP) MAF wave functions just as in JWKB method, but now it leads to two different MAF quantization formulas with two different MAF universal constants regarding EP and OP solutions as given in [3] and as we study in this section. We again use the symbolism in (9) (φ⇔ψ) and start with the first quadrant, by applying that lim y!∞ ¼ 0 requires c 2 ¼ 0 in (42), namely, ψ MAF, n λ n ; y ð Þ≕ψ Mn λ n ; y ð Þ ¼c 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ∂ y ξ λ n ; y ð Þ p Ai ξ λ n ; y ð Þ ½ , where the denominator can be written in the following form [3]: i. Even-parity (EP) eigenenergies: if we apply the EP formulas of the exact solution in (17), by using (16) where we used ξ 0 , ξ 0 0 , and ξ ″ 0 for simplification, namely, ξ 0 ¼ ξ λ n ; 0 ð Þ, ξ 0 0 ¼ ∂ y λ n ; y ð Þ y¼0 , and ξ ″ 0 ¼ ∂ yy λ n ; y ð Þ y¼0 , respectively. We assumed here φ λ n ; 0 ð Þ¼1, and we will use then parity correction for φ λ n ; 0 ð Þ¼À1 case just as in the JWKB calculations. If we take the derivative of (36), we get.
where the last term vanishes as y ! 0 to give.
Now, by the substitution of ξ 0 ! ÀZ sn , we have.
where the subscripts sn stand for s, symmetrical solution (EP), and n, quantization order (nth quantization), and Z sn is the nth solution of the differential equation in (52) regarding the symmetrical solution. Now, by using the results in (13), we find the MAF quantization formula regarding the symmetrical solution: where ζ sn is the universal MAF constants regarding the symmetrical solution whose values are given in Table 2 along with the JWKB solutions (which are already exact) for some n values in comparison. Note that we used n s to represent the symmetrical (EP) MAF indices in Table 2.
ii. Odd-parity (OP) eigenenergies: similarly, if we apply the OP formulas of the exact solution in (17), by using (16), to the MAF wave functions, we have the following:ψ ii ð Þ ∂ yψ M λ n a ; y (56) where, similarly, the subscripts an stand for a, antisymmetrical (OP), and n, quantization order (nth quantization), and Z an is the nth solution of the equation in (55) regarding the asymmetrical solutions. Similarly, by using the results in (13), we find the MAF quantization formula regarding the antisymmetrical solution: where ζ an is the universal MAF constants regarding the antisymmetrical solution whose values are given in Table 2 along with the JWKB solutions (which are already exact) for some n values in comparison. Note that we used n a to represent the antisymmetrical (OP) MAF indices in Table 2. Similarly, calculation of the constant coefficient in (60) for the antisymmetric boundary values given in (54) with q = 1 gives. O:P: :c 1a ¼Ẽ 7=6 a n 3π ð Þ 1=6 ffiffiffiffiffiffiffi 2β 3 pẼ 4=3 a n 3π ð Þ 2=3 Ai 0 À 1 4Ẽ 4=3 a n 3π ð Þ 2=3 h i À Ai À 1 4Ẽ 4=3 a n 3π ð Þ 2=3 h i n o (66) Since the MAF solutions of both EP and OP solutions are very close to the exact solutions given in Figures 3 and 4, their absolute and relative error graphs with respect to the exact solution are given in Figures 5 and 6. We can also see that there are no discontinuities at the CTPs in the MAF solutions when compared with the JWKB solutions given in Figures 3 and 4.

Conclusion
Here we studied the fundamental outcomes of the two conventional semiclassical approximation methods, namely, JWKB and MAF methods pedagogically, and obtained the solutions of the QHO by these semiclassical methods by using the parity conditions of the expected solutions by using the dimensionless form of the QHO system. We applied the asymptotic matching and parity matching procedure to obtain the correct form of semiclassical solutions. As expected, JWKB solutions diverge at and around the CTPs, whereas MAF solutions do not. As also expected (since being typical), JWKB eigenenergies are exact, whereas MAF eigenenergies are unfortunately not but very accurate as expected from an approximation method. In the MAF method, function p in (40) or in (43) is assumed zero. Indeed, it is very close to zero to give approximate results, and function P in (40) or in (43) can be used as an approximation criterion for the MAF method [3,18]. However, improved MAF methods (IMAF) or perturbation corrections concerning the nonzero P function seem straightforward to improve the accuracy of the MAF solutions as in [3,20,22]. Normally, for an even potential function in the TISE, EP and OP initial values are as given in (17), but due to the conversion factor β in (11) or (16), for the QHO in the dimensionless form (in ψ), we have (22) and (29). In our notation, we have used the notation, φ⇔ψ, where real physical system is in φ and the dimensionless form is in ψ. Since the standard formulation is given according to the real physical systems, JWKB and MAF formulas in the literature such as in [1][2][3][4][5][6][7][19][20][21][22][23] surely correspond to the initial values β ! 1 in our dimensionless form formulation in ψ. Consequently, we hereby present a full JWKB and MAF solution concerning the quantized conversion factor β in (11).