Quantum Anharmonic Oscillators: A Truncated Matrix Approach

This study aims at implementing a truncated matrix approach based on harmonic oscillator eigenfunctions to calculate energy eigenvalues of anharmonic oscillators containing quadratic, quartic, sextic, octic, and decic anharmonicities. The accuracy of the matrix method is also tested. Using this method, the wave functions of the anharmonic oscillators were written as a linear combination of some finite number of harmonic oscillator basis states. Results showed that calculation with 100 basis states generated accurate energies of oscillators with relatively small coupling constants, with computation time less than 1 minute. Including more basis states could result in more correct digits. For instance, using 300 harmonic oscillator basis states in a simple Mathematica code in about 8 minutes, highly accurate energies of the oscillators were obtained for relatively small coupling constants, with up to 15 correct digits. Reasonable accuracy was also found for much larger coupling constants with at least three correct digits for some low lying energies of the oscillators reported in this study. Some of our results contained more correct digits than other results reported in the literature.


Introduction
Quantum anharmonic oscillators have long been used to test the power and shortcomings of new approximation techniques proposed to solve the Schrödinger equation of quantum systems. They are also often used in testing computational approaches originally designed for systems with many fermions [1]. Moreover, anharmonic oscillators can be used to represent various challenging potentials, such as the double-well potential, which is very often used in theoretical physics studies [2]. More importantly, they were found to be useful in modeling many phenomena in nuclear physics, solid-state physics, atomic and molecular physics, and laser theory [3][4]. Their importance is mainly due to the anharmonic nature of vibrations of many quantum systems [5], ranging from that of diatomic molecules [6][7] to extended solids [8].
Various analytical and computational approaches have been developed to calculate the energies of the anharmonic oscillators. Some of the approaches include an algebraic method based on the ladder operator [9], analytic quasilinearization method [10], Lie algebra [11,12], the Poincare-Linstedt method [13], multiple-scale perturbation theory [14], Wick's normal ordering technique [15], examination of polynomial solution [16], quantum Monte Carlo method [17], and pertur-bation theory [18]. Many other approaches have also been developed to calculate the energies of the systems.
Some approaches, including the semi-classical and standard perturbation theory, have limited success in calculating the energies of the anharmonic oscillators due to the divergence issues even for small coupling constants [10,[19][20][21][22]. It is also well known that for large coupling, the perturbation expansion becomes worse and results in very poor accuracy [10]. This means that utilizing an effective and simple method with a high level of accuracy is still challenging [2].
A simple truncated matrix approach presented in [23][24][25][26] was successfully applied to obtain accurate ground state and some excited state energies of the helium atom and its isoelectronic series. Recently, energies of linear oscillator and quartic oscillators with potential ( ) = 2 + 4 were determined using this method with a matrix of size 25 for fixed values of = 0.05 and different values of ranging from 0 to 1.9 [27]. This method was also implemented in a simple Matlab code to obtain some energies of a pure quartic anharmonic oscillator with fixed and matrix of size ranging from 7 to 40 [28]. However, in their study, only some energies of the quartic oscillators with one coupling constant were reported. Hence, the accuracy of the method for other anharmonic oscillators with various coupling constants is still unknown. Therefore a comprehensive study utilizing this method to calculate energies of various anharmonic oscillators with a much larger range of coupling constant values is of great importance in testing the accuracy of the method in such systems.
To the best of our knowledge, the matrix method [23][24][25][26] has not been applied to calculate energy eigenvalues of other quantum anharmonic oscilla-tors containing sextic, octic, and decic anharmoni-city. In addition, despite being reported for quartic oscillators in [27,28] and linear oscillators in [27], their results were not enough to test the accuracy of the method for the oscillators since limited values of coupling constants were used.
Hence, it is the main purpose of this research to implement the matrix method based on harmonic oscillator basis states in a simple Mathematica code to obtain accurate ground-state energies and some excited state energies of quadratic, quartic, sextic, octic, and decic anharmonic oscillators with a large range of coupling constants. In addition, the accuracy of the method is tested for the systems.

Computational Details
By treating anharmonicity terms as a perturbation, the Hamiltonian of the anharmonic oscillators can be written Where Ĥ 0 is the Hamiltonian of the harmonic oscillator, which in the study of anharmonic oscillators is written in its dimensionless form as: Using this expression, Ĥ 0 is in the units of ℏω/2, p is in the units of √ℏmω and x is in the units of √ℏ/mω. As a result, energies are in the units of ℏω/2. These units are used throughout this article. It is important to note that some papers defined equation (2) in a slightly different way. Therefore, the conversion of units is vital to compare the energies of the oscillators. in equation (1) is the perturbation parameter which will be set equal to 1 at the end and Ĥ′ is the perturbation term which in this research is of the form: Where = 2, 4, 6, 8, 10 correspond to quadratic, quartic, sextic, octic, and decic anharmonicities, respectively and are the coupling constants defined to be 2 = , 4 = , 6 = , 8 = and 10 = . More explicit forms of equation (3) were shown in equations (3a), (3b), (3c), (3d), and (3e) for the respective oscillators. Figure 1 illustrates the oscillators.
Substituting equation (5) into (4) and multiplying from left by 〈 |, one obtains the reduced form of the Schrödinger equation in N dimension: The elements of the Hamiltonian matrix Hji can be obtained by the standard method, i.e.
H ji = 〈j|Ĥ 0 |i〉 + ∑ c m 10 m=2,even 〈j|x m |i〉 (8) The first term in equation (8) is the well-known matrix elements of the unperturbed harmonic oscillators which have the following form: Furthermore, the second term can be easily calculated using the fact that matrix elements of x are well known from quantum mechanics classes, which using the units here are of the form: The 〈j|x m |i〉 term in equation (8) can be easily and effectively evaluated utilizing a very convenient matrix multiplication offered in Mathematica.
Once the Hamiltonian matrix elements of the anharmonic oscillators shown in equation (8) are determined for a particular matrix size N with particular values of coupling constants, the matrix is then diagonalized to obtain N energy eigenvalues of the oscillators, which was performed using a Mathematica code in this research.

A. Hamiltonian matrix representation
Matrix representation of the Hamiltonian of the anharmonic oscillators was performed by solving equation (8) using equations (9) and (10). This could be easily done "by hand" using very simple algebraic manipulations involving ladder operators as discussed in any undergraduate quantum mechanics class, which can also be easily found in most lecture notes and problem sets of undergraduate quantum mechanics. On the one hand, this "by hand" calculation is important to help undergraduate students understand simple algebraic manipulations involving Dirac notation and creation and annihilation operators. On the other hand, for larger values of in equation (8), this "by hand" type calculation becomes less effective since the same procedure is repeated with no new physical significance observed from such repetition. Additionally, such "by hand" calculation can only be used to diagonalize the Hamiltonian of relatively small size, which (of course) leads to less accurate energy eigenvalues of the systems. It is for these reasons that we used Mathematica to do the job. Only some elements of the Hamiltonian matrix of the form of equation (8) are shown in equation (11) due to the limited space. The Mathematica codes, which can be downloaded from the supplementary material, can be easily adjusted to deal with any order of anharmonicities and matrix of any size. The codes not only determine the Hamiltonian matrix but also calculate the energies of the oscillators.

B. Convergence of Energies as a Function of Basis States
Diagonalizing equation (11) resulted in N energy eigenvalues of the systems. Table 1 shows the first five energies of the anharmonic oscillators using a different number of basis states for pure quadratic ( = 1.0), pure quartic ( = 0.2), pure sextic ( = 0.2), pure octic ( = 0.2), and pure decic anharmonicity ( = 0.2 ). Despite the fact that N energy eigenvalues were obtained, only the first five energies were reported due to the limited space provided. Exact values are taken from [1] for quartic, sextic and octic oscillators and from [15] for the quadratic oscillators.
It is clear from Table 1 that increasing the number of basis states does have a major effect on the accuracy of the energies obtained, as shown by the increase in the number of decimal places that are in agreement with the exact values (the underlined digits). This effect became more significant as one goes to higher order anharmonicities. For instance, when using 5 basis states, the ground state energies of quartic oscillators agree with the exact values to 3 digits while only 2 correct digits are found for the sextic oscillator and 1 correct digit for the octic oscillator.
The exact values for the decic oscillator were not presented in [1] and hence were not compared to our results. It is also obvious that when using 100 and 200 basis states, all the results agree with the exact values for all anharmonicities up to 8 and even to 9 decimal digits, indicating that our results are accurate. It might be true here that adding a basis from 100 to 200 only has a minor effect on the accuracy of our calculation. However, it was found that for much larger coupling constants, this has a more significant effect. Therefore, in the calculation presented in the next section, we used 300 basis states.
It is also obvious from Table 1 that using a relatively small number of basis states resulted in less accurate energy eigenvalues for the oscillators. For example, using N=5, which might be considered as the largest number of basis states that one could use to diagonalize the Hamiltonian matrix using 'by hand' calculation explained in the previous section, errors obtained for the first four energy levels of the quadratic oscillators were about 0.06%, 0.40%, 2.88%, and 8.49% respectively. This increasing trend in the errors as one goes to higher energies strongly suggested that 'by hand' type calculation using Hamiltonian of relatively small size can not be used if one wants to obtain highly accurate energies, especially for higher state energies of the anharmonic oscillators. Therefore, the use of modern computational techniques such as Mathematica is highly recommended to obtain much more accurate energies.  Figure 2 shows the convergence of the energies as a function basis states for the sextic oscillators. As can be seen, from N=20 onwards, the energies converged more rapidly, in agreement with the results of [28] for the quartic oscillators. However, as shown in Table 1, to obtain much more accurate energies, one needs 200 basis states for relatively small coupling constants. It was also found that using 300 basis states for much larger coupling constants significantly improved the accuracy of the calculation.
In terms of the computation time, a larger number of basis states resulted in longer computation time, as expected. For N=5, it took about 0.12 seconds to obtain the energies of the oscillators, but a larger deviation from the exact values was obvious from Table 1 and Figure 2, especially for larger quantum numbers n. Accurate results were obtained with only 20 basis states, which were obtained in 1.62 seconds. Accuracy was significantly improved when using 100 basis states with a computation time of approximately 48.34 seconds. Meanwhile, calculations involving 200 and 300 basis states added more correct digits to the results and took about 4.35 minutes and 8.54 minutes, respectively. These results showed that the truncated matrix method implemented in this research was not only effective and simple, but also accurate. The computation time presented here can be faster or lower depending on the computer used. In this study, AMD E1-1200 APU with Radeon (tm) HD Graphics 1.40 GHz with 2 GB RAM was used in performing the calculation.
By considering the accuracy and numerical efforts needed, it is highly recommended to use 100 basis states for quick calculations, especially for oscillators with relatively small coupling constants. However, if one needs to generate more correct digits, especially for those with relatively large coupling constants, calculation using 300 basis states can be used.

C. Energies of the Quadratic and Quartic Oscillators
The first few energies of the quadratic and quartic oscillators for various coupling constants using N=300 and 30 digit precision calculation are shown in Table 2, compared with corresponding results using different methods in the literature. It was important to note that some coupling constants values taken from literature and presented here may look different from those in the original papers but they are actually equivalent after conversion of units. This is due to the fact that while some papers used the same units as units in this article, some others used slightly different units. As a consequence, some energies are multiplied by 2 to have the same units as ours. Table 2 clearly shows that for relatively small coupling constants, energies from our calculation are in excellent agreement with the most accurate energies in the literature. Moreover, we reported more correct digits in our results compared to some other results in the literature. For much larger coupling constants, our approach results in less accurate energies but is still comparable to other approaches. For instance, using β=2000 for the ground state energy, we have 6 correct digits, which are the same as the results of [30].   [28] used only 40 basis states compared to 300 in this article. Therefore, it is clear that we reported more correct digits, and hence underlining correct digits is irrelevant for this particular case. Table 2 clearly shows that for relatively small coupling constants, energies from our calculation are in excellent agreement with the most accurate energies in the literature. Moreover, we reported more correct digits in our results compared to some other results in the literature. For much larger coupling constants, our approach results in less accurate energies but is still comparable to other approaches. For instance, using β=2000 for the ground state energy, we have 6 correct digits, which are the same as the results of [30].

D. Energies of Sextic, Octic, and Decic Oscillators
Energies of sextic, octic, and decic oscillators for various coupling constants using N=300 are shown in Table 3, Table 4, and Table 5, respectively. The energies are also compared with correspondding results using different methods in the literature. The comparisons show that as in the case of quadratic and quartic oscillators, the energies are very accurate for small coupling constants, and the accuracy decreases as coupling constants increase    [32] but is still comparable to the accuracy of some other results in the literature.

E. Accuracy of The Truncated Matrix Approach
From Table 3, Table 4, and Table 5, it can be seen that, in general, there is a declining trend in the accuracy of the matrix approach as the coupling constants increase. For instance, for the ground state energy of the quartic oscillator, when the coupling constant is 0.2 there were 15 correct digits in our results, which were reduced to 10 and 6 correct digits for coupling constants of 2 and 2000, respectively. This is also true for the case of quadratic, sextic, octic, and decic oscillators. This is expected since the Hamiltonian used in this article was written in terms of an unperturbed Hamiltonian and a perturbation which was assumed to be much smaller than the unperturbed Hamiltonian. Despite this declining trend in accuracy, we still reported reasonable accurate energies for oscillators with large coupling constants.
In addition, coupling constants ranged from -11 to 2000 have been tested in our calculation, and the results showed that at least 3 correct digits were found in our energies compared to the exact values in the literature. This strongly indicates that despite using a simple calculation, our results can produce highly accurate energies for small coupling constants and reasonable accuracy for large ones. This also indicates that the divergence issues experienced by the ordinary perturbation expansion method have not been found in this study.
Finally, some of our results contained more correct digits than some other methods in the literature. For example, using [31] for the exact ground state energy of the quartic oscillator, there are 15 correct digits obtained from our calculation compared to only 2 correct digits reported in [9], 4 correct digits presented in [11], and 6 correct digits in [30].

Conclusion
Accurate energies of quadratic, quartic, sextic, octic, and decic oscillators with various coupling constants were obtained using a simple matrix approach using harmonic oscillator eigenfunctions. Highly accurate energies of all anharmonic oscillators were obtained for relatively small coupling constants. It was also found that the accuracy decreased as the coupling constants increased but a reasonable degree of accuracy for much larger coupling constants was still obtained for relatively low energy levels of the oscillators.