A Non-relativistic Approach to Relativistic Quantum Mechanics: The Case of the Harmonic Oscillator

A recently proposed approach to relativistic quantum mechanics (Grave de Peralta, Poveda, Poirier in Eur J Phys 42:055404, 2021) is applied to the problem of a particle in a quadratic potential. The methods, both exact and approximate, allow one to obtain eigenstate energy levels and wavefunctions, using conventional numerical eigensolvers applied to Schrödinger-like equations. Results are obtained over a nine-order-of-magnitude variation of system parameters, ranging from the non-relativistic to the ultrarelativistic limits. Various trends are analyzed and discussed—some of which might have been easily predicted, others which may be a bit more surprising.


Introduction
The harmonic oscillator (HO) is a fundamental problem in quantum mechanics, which serves as an elementary description of physical phenomena as varied as the vibrational modes of molecules and solids, or the modes of the electromagnetic field. The most prominent feature of the quantum HO is its equally spaced energy levels, bringing a direct linkage to the concept of quantum field excitation.
Being an exactly solvable model, the quantum HO is discussed in every introductory quantum mechanics textbook [1][2][3][4]. For definiteness, we take the solutions of the non-relativistic quantum HO problem to be the eigenstate solutions of the following ordinary differential equation-i.e., the stationary Schrödinger equation: Here, m is the mass of the particle, k is the curvature of the potential and is the momentum operator, as per usual notation.
From Eq. (1) the energy of the nth level will be with = √ k∕m . Eq. (3) describes a linear relation between the energy and the quantum number n-a unique property of a non-relativistic particle trapped in a quadratic well, hearkening back to the linear relation between energy and action in the classical HO. In principle, n-and therefore E n -extends upwards ad infinitum.
On the other hand, the equipartition theorem ensures that infinite energy implies infinite kinetic energy (since kinetic and potential energy expectation value contributions are equal)-which in turn implies a finite n beyond which relativistic effects cannot be neglected. What happens then? Whatever the answer, this simple argument makes it clear that the large n limit should be regarded as one path towards the relativistic limit.
The relativistic quantum HO problem has been addressed previously [5][6][7] using the time-independent Klein-Gordon equation for a particle in a gauge potential [8,9]: Note that here and throughout this work, we assume units where c = 1 , which allows us to remove the symbol "c" from the formulation. In effect, this leaves m as the relativistic parameter. More specifically, since HO energies increase with decreasing m, the m → 0 limit provides another path to the relativistic limit.
Returning to the Klein-Gordon equation, the important observation here is that Eq. (4) is equivalent to the non-relativistic Schrödinger equation (e.g., Eq. 1) for a particle with effective energy (E 2 ∕2m − m∕2) , moving in an effective potential, Note that in the non-relativistic or large-m limit, m ≈ E , and so Eq. (5) approaches the familiar quadratic form of Eq. (1). (1) More generally, Eq. (5) describes a potential energy that is unbounded from below as x → ±∞ . As a consequence, the positive-energy solutions of the Klein-Gordon oscillator correspond to resonances [5][6][7]. Moreover, even when the anharmonic x 4 term in Eq. (5) is a small perturbation [5], such that tunneling through the barrier can be neglected, the aforementioned isospectral (linear) feature of the HO spectrum is destroyed, leading to a decrease in spacing between levels with increasing energy or n [5]. Of course, the general Klein-Gordon solutions are well known to exhibit some other undesirable features such as negative energy solutions, that will not be dwelt on further here.
Instead, we point out some less well publicized features of the Klein-Gordon HO solutions that will be of interest for the present work. To begin with, the solutions as described above are in close analogy with the behavior of a classical relativistic particle in a quadratic potential. In the classical case, the particle may also be thought to move in an effective potential analogous to Eq. (5)-but with the height of the barrier increasing with the particle energy, such that the motion remains always bounded and periodic [10]. Furthermore, the period of the oscillation increases with energy, due to time dilation along the world line-indicating that the behavior of a particle in a quadratic potential loses its "harmonic" character in the relativistic regime [10].
Finally, we point out another interesting feature, also implied by the above discussion-this being that the effective potential itself depends on the energy value, E. The classical ramifications of this have already been discussed; quantum mechanically, this implies that the non-relativistic-like form of Eqs. (4) and (5) actually describes a nonlinear Schrödinger equation-for which standard self-consistent or mean field solution methods are well established, e.g. in electronic structure [11,12].
In this work, instead of working directly with the Klein-Gordon oscillator equation, we use the "square root" form, as follows [8,9]: Note that the relativistic rest energy has now been subtracted from both sides, in order to result in an energy E that can be directly compared with the usual non-relativistic result of Eq. (3). Equation (6) above (but generalized for arbitrary four-vector potentials) is also often encountered in textbooks, as a covariant relativistic eigenvalue equation. However, it is often dismissed as "impractical", due to the squareroot operator. That said, Eq. (6) is evidently closely related to the so called "spinless Salpeter equation", a relativistic model widely used to describe hadrons as bound states of constituent quarks [13]. In order to avoid the problematic square-root operator, the solutions of this equation have been obtained using non-relativistic-like effective hamiltonians [13,14] or the auxiliary fields formalism [15], with which the present method bears some resemblance.
Recently, we have developed some practical tools for solving Eq. (6) directly-both (numerically) exactly, and approximately [16], for arbitrary scalar potentials. In particular, the previously developed approximate methods rely on a nonlinear Schrödinger or iterative mean-field-type approach, not unlike that described above in the Klein-Gordon context. Such methods require that a conventional Schrödinger equation be solved multiple times. In the present approach, we develop new methods that-at least in the HO context-allow us to the "cut to the chase," by bypassing such multiple rounds of explicit Schrödinger solution.
In the context of exact numerical solutions for the relativistic quantum HO, analytical Hamiltonian matrix elements can be obtained, if one uses a particle-in-a-box basisset representation. Moreover, with this choice, the sign ambiguity of the square root operator completely disappears, as the kinetic energy matrix becomes diagonal and positive-definite. Numerically exact, quadrature-free, positive-energy-and above all, bound-solutions are thus readily obtainable.
In this work, we shall use precisely this approach to compute the energy levels and wavefunctions for the low-lying relativistic quantum HO eigensstates-with high numerical accuracy, in a manner that admits rigorous characterization of basis set truncation error (the only source of error). Note that this approach is rigorously variational. The specific problems considered range from the non-relativistic to the ultrarelativistic, by varying the mass parameter over the very broad interval, 10 −6 ≤ m ≤ 10 3 (in units where ℏ = k = 1).
The exact numerical solutions as obtained above then serve as a benchmark to validate the analytical approximations [16]. As might be expected, it will be shown that exact and approximate relativistic quantum solutions approach Eq. (3), and each other, in the non-relativistic limit where E ≪ m . Conversely, in the ultrarelativistic or m → 0 limit, a radically different, highly nonlinear energy spectrum is observed, for which two simple, but approximate, analytic expressions are derived.
The remainder of this paper is organized as follows. In Sect. 2, the numerical methods used to solve Eq. (6), both approximate and exact, are presented in detail. These are then applied to the HO system, respectively, in Sects. 3 and 5-with the latter section also including comparisons across m and n values. Finally, some concluding remarks are presented in Sect. 6.

The Methods
For the conceptual details of this approach in its general time-dependent formulation, the reader is referred to [16]. Here, we restrict consideration of the derivations to the time-independent or stationary solutions only.
The approach starts by rewriting Eq. (6) in the form which shall be called the stationary "Poirier-Grave de Peralta" (PGP) equation, following Ref. [16]. In Eq. (7) above, ̂ is the dimensionless square-root operator, and V corresponds to the "scalar" time-like component of the full electromagnetic four-potential-taken here to be the quadratic HO potential. With this interpretation, Eq. (7) becomes fully Lorentz covariant [8], although it does involve a square-root operator.
One strategy for avoiding the square root was suggested by Poirier (see the Appendix of Ref. [16]), who noted that the kinetic energy operator of Eq. (7) can be rewritten as a continued fraction, leading to the following appealing equation: Eq. (9) above is referred to as the "Poirier equation." It is exact, and furthermore, avoids the square root sign ambiguity, regardless of the choice of representational basis set (i.e., particle-in-a-box or otherwise). It also suggests a converging sequence of numerical approximations, obtained by truncating the continued fraction at successive orders [16,17]-which, as will be demonstrated in a future paper, can be made to converge extremely rapidly.
Alternatively, the formal similarity between the operator p 2 ∕(1 +̂)m and the nonrelativistic kinetic energy operator p 2 ∕2m , led Grave de Peralta to propose the following equation: Here, the troublesome square-root operator (8) is simply substituted by the corresponding Lorentz parameter: Equation (10) is the so-called "general Grave de Peralta" (gGP) equation. Its similarity with the Schrödinger equation allows standard methods, as used for well-known non-relativistic quantum mechanical problems, to be applied directly in the relativistic domain as well [18][19][20][21][22]. This requires, however, that some specific value be chosen for the parameter . Different values yield different solutions-hence the use of the term "general." Note that regardless of how is chosen, the resultant gGP equation is in general an approximation only; thus, its solution, gGP is not equal to the corresponding exact solution, , of Eq. (7).

Poveda Parametrization and the Grave de Peralta-Poveda-Poirier Equation
Although many different schemes could in principle be adopted for choosing a suitable value for the gGP approximation, one of the most natural must surely be the expectation value of the ̂ operator itself: This is what we call the "Poveda parametrization" [16]. In practice, Eq. (12) cannot be used to determine , since the exact is not known.
To make progress, we replace the exact in Eq. (12) with the specific solution of the gGP equation that corresponds to the resultant value of itself. This choice of gives rise to the "Grave de Peralta-Poveda-Poirier" (GPPP) approximation, GPPP , satisfying the following GPPP equation or self-consistency condition: The GPPP equation is a non-linear Schrödinger equation-expressed by the fact that the Hamiltonian operator itself depends on the solution wavefunction. Such equations may be solved using an iterative self-consistent procedure of the type used routinely, e.g., in Hartree-Fock calculations [11,12]. This method was already recently presented [16] and will not be described further here.
On the other hand, for the present HO application, there is a much faster, more direct way to obtain exact solutions to Eq. (13), which does not involve converging an iterative sequence. This stems from the realization that regardless of the choice of , all gGP solution wavefunctions, gGP , are simply variations of the familiar non-relativistic Hermite-polynomial-based solutions, but with rescaled masses, m → = (1 + ) m∕2 . Thus, where n (x; ) are the usual non-relativistic solutions from Eq. (1), but with m replaced with .
Whereas for gGP, -and therefore -can take on any value, to obtain the exact GPPP value, we need merely solve the following equation: The expectation value in Eq. (15) above is particularly easy to evaluate in the momentum representation-as the matrix representation of ̂ becomes diagonal and analytical, and the corresponding HO ̃n(p; ) forms are of course also well known analytically. Note that the resultant GPPP and GPPP values are state-specific, meaning that they vary with the quantum state index n.
In Sect. 5 we will present exact ground state HO solutions of the GPPP equation (Eq. 13)-as obtained using the approach described above. For the most part however, we shall concentrate more on a new analytical approximation to GPPP, to be described in Sect. 2.2-with results presented in Sect. 3.

Self-consistent (SC)-GPPP Approximation
Consider the following rearrangement of Eq. (6): Equation (16) above certainly holds for the exact relativistic solution, . However, if GPPP is a good approximation, then the relation should also hold approximately for GPPP as well. By replacing with GPPP in Eq. (16), multiplying both sides on the left by ⟨ GPPP | , and finally, making use of Eq. (12), we obtain the approximate relation: Next, we substitute the Eq. (17) formula for GPPP into Eq. (13), and again multiply on the left by ⟨ GPPP | to obtain Now then, a key advantage of Eq. (18) is that it provides relativistic energy levels directly in terms of expectation values-not of the square root operator, as in Eq. (15), but of the much simpler non-relativistic Hamiltonian components, V and p 2 . Of course, for all gGP solutions of the HO, analytic expressions for these expectation values are well known. We are therefore motivated to define a new approximation, wherein GPPP in Eq. (18) is replaced with a different gGP solution (i.e, different value) such that both Eqs. (17) and (18) are now satisfied exactly.
The above defines what we call the "self-consistent" (SC)-GPPP approximation. It is socalled because the SC-GPPP gGP values may be solved for in self-consistent fashion [16]-similar to, but more straightforwardly than, the method outlined at the end of Sect. 2.1. This can be achieved as follows. First, Eq. (18) is replaced with the following exact SC-GPPP expression, for which all expectation values are presumed to be taken with respect to SC−GPPP : This leads at once, via recursive expansion, to a continued fraction form: Of course, Eq. (20) is highly reminiscent of Eq. (9), except that all operator components are replaced with expectation values. Alternatively, we can "close" the continued fraction series, in order to generate an expection-value-based approximation to the original exact square root operator form of Eq. (6): Thus, in the SC-GPPP approach, the relativistic energy levels manifest the correct classical relativistic relation (through expectation values) with the particle's potential energy and linear momentum squared.
In the specific case of the HO system, moreover, these quantities may all be obtained analytically, once the correct SC-GPPP value of has been established. We examine this issue further in Sect. 3.

Analytical SC-GPPP Results for the HO System
From Eqs. (10) and (14), all solutions, n (x; n ) , of the HO gGP Hamiltonian must have energies given by where an n subscript has been added to , to reflect its state-specific nature (Sect. 2.1). Now, from the GPPP formulation, the value of n (and therefore n ) is obtained from the expectation value of ̂ , as in Eq. (15).
Substituting Eq. (15) into (22) then leads to Next, we apply the SC-GPPP condition-i.e., the exact equality version of Eq. (17)-to replace Eq. (23) above with . Now, considering that all n (x; n ) are generic HO solutions, it is clear that the expectation value of the potential energy is one half of the total energy, so that Substituting into Eq. (24) then yields At this point, it is straightforward to convert Eq. (26) into the following cubic equation, with the understanding that the desired solution corresponds to a real positive root: Equation (26) can be easily solved-even analytically! Numerically, one could start with the large-m or non-relativistic limit result, i.e. E n ≈ E n , and solve for E n iteratively. Amusingly, this gives rise to the following continued-fraction-like expansion, known as a "nested radical": The form of Eq. (28) makes it clear that E n < E n always, with the discrepancy increasing with increasing n and or decreasing m. As noted elsewhere [13], the non-relativistic kinetic energy p 2 ∕2m represents an upper bound of the relativistic one √ m 2 + p 2 − m , which can be easily verify by squaring the inequality √ m 2 + p 2 − m ≥ 0 . In this sense a Schrödinger eigenvalue, for a given potential, would be an upper bound of the corresponding eigenvalue of Eq (6), with the same potential. Moreover, the approximated SC-GPPP energy of Eq. (21) tend to overestimate the exact value, considering the fundamental inequality: . In the ultrarelativistic limit, (E n ∕E n ) → 0 , but the precise manner in which this occurs is quite illuminating to consider. From Eq. (27), we see that in the m → 0 limit, Equation (29) above leads to two notable conclusions about the ultrarelativistic limit. First, instead of a linear relationship, the energy levels E n exhibit a sublinear n 2∕3 power law scaling in n. This implies that the energy level spacing decreases with n with an inverse cube root dependence ( n −1∕3 ). Second, whereas the E n scale with mass as m −1∕2 , the ultrarelativistic limit E n become entirely mass-independent. Though the above predictions may perhaps seem surprising, they also emerge from a relativistic phase space analysis of the classical action, as will be shown in Sect. 4, wherein an exact analytical formula is also derived. It is worth noting that the above n 2∕3 behavior of the ultrarelativistic energy of a particle in a quadratic potential, was previously obtained using the previously mentioned effective hamiltonian [13,14] and auxiliary fields methods [15]. Figure 1 shows the three lowest SC-GPPP energy levels for the ℏ = k = 1 HO system, as a function of m, with m ranging from 10 −4 to 10 3 . These curves were computed as solutions to Eq. (27). The corresponding non-relativistic E n curves are also indicated. Note that the latter diverge as m −1∕2 in the m → 0 limit, as discussed. The E n curves, on the other hand, approach finite upper bounds-i.e., mass-independent constants, as predicted by Eq. (29). Moreover, the level spacing decreases with increasing n, also exactly as predicted.

Approximate WKB Results, and Exact Analytical Results for m → 0
In the standard Wentzel-Kramers-Brillouin (WKB) semiclassical theory [23][24][25], approximate energy levels are obtained through the classical actionenergy relationship, S(E), where S is the phase space volume contained within the classical Hamiltonian contour of energy E. With S(E) known, the discrete eigenenergies, E n , are then obtained from the usual half-integer quantized action values, S n = 2 ℏ(n + 1∕2) . This relation is presumed to hold even when the Hamiltonian is relativistic. Though the WKB approximation of the Schrödinger equation is a common topic of many quantum mechanics textbooks, its application in the relativistic domain is lacking. To the author's best knowledge, the first generalization of the WKB method to relativistic Hamiltonian was presented in ref. [26], where the WKB probability distributions for Hamiltonians involving arbitrary kinetic energy are derived. In particular the kinetic term √ m 2 + p 2 and its m → 0 limit were discussed in [26], in the context of the HO.
A bit later in this section, we will derive an analytic WKB energy eigenvalue expression based on the classical version of Eq. (6), for arbitrary m and k values. First, however, we find it instructive to consider the ultrarelativistic m → 0 limit. In this limit, Eq. (6) reduces to Note the absolute value bars that appear around p ! Note also that the parameter m no longer appears anywhere in this expression. These details are very important.
From Eq. (30), the classical Hamiltonian contours for energy E are found to be This is easy to integrate over x, between the two turning points at x tp = ± √ 2E∕k , to obtain the requisite classical action-energy relationship, The semiclassical quantization condition applied to Eq. (32) then yields Note that Eq. (33) above is identical to Eq. (29), except with a slightly different (and smaller) prefactor-which is reduced from 2 to about 1.66608. In any event, the above WKB analysis further confirms the validity of the E n ∝ m 0 n 2∕3 scaling law in the m → 0 limit.
The above n 2∕3 scaling law will look very familiar to those schooled in uniform semiclassical approximation theory [23][24][25]. In this approach, exact analytical solutions for linear-potential systems are used in the vicinity of classical turning points, in order to "regularize" semiclassical solutions between the tunneling and classically allowed regions of space. Of course, those analytical solutions are well-known to be Airy functions. On the "classically allowed" side of the turning point, the Airy functions oscillate, with varying phase that exhibits the 2/3-power scaling law in the semiclassical limit [23][24][25].
But what does this have to do with our current, relativistic application? When viewed from a Fourier-transformed perspective, Eq. (30) clearly becomes a non-relativistic Schrödinger equation with "mass" m mom = (1∕k) and "potential energy" V mom (p) =| p | . This is not quite a linear potential, but almost; in particular, the momentum-space solutions are comprised of piecewise Airy functions on either side of p = 0 [27]. Moreover, this "V-shaped" momentum potential clearly gives rise to discrete bound states, whose boundary conditions at p = 0 imply that this boundary point must either correspond to an Airy function extremum or zero-corresponding, respectively, to the even-or odd-n solutions. Such analysis leads to the following, remarkably simple exact analytic formula in the m → 0 limit, where a n denotes the location of the corresponding Airy extremum or zero, as appropriate. Table 1 presents WKB and exact energy levels for the ultrarelativistic or m → 0 Hamiltonian of Eq. (30), for several low-lying energy levels, as computed using Eqs. (33) and (34), respectively. As is typical for WKB approximations, the worst error is for the ground state energy-although even here, the WKB prediction overestimates by less than 10%. Error magnitudes decrease very quickly thereafter with increasing n-as is also typical for WKB. A bit more interesting is the fact that this decrease is not entirely monotonic, and also changes sign with even/odd nno doubt due to the alternation from Airy function extrema to zeros. All in all, the WKB results are seen to comprise an excellent approximation to the exact results.
We conclude this section with a derivation of the WKB energy levels for the general, or arbitrary-mass, HO system. The semiclassical procedure followed is exactly

Exact Numerical Results, and Comparisons
As described in Sect. 1, we have also performed exact numerical calculations of the exact relativistic quantum HO eigenstate solutions of Eq. (6). This can easily be achieved, despite the square-root operator, if a particle-in-a-box basis set is used-in terms of which the kinetic energy matrix representation becomes diagonal (and positive-definite). Likewise, the matrix representation of the potential energy-though not diagonal-is nevertheless exact, leading to perfectly variational convergence with increasing basis size. In this sense, the present method appears to be more straightforward, at least for the quadratic potential, than the previous methods based on matrix [28,29] or mesh [30,31] representations of the spinless Salpeter equation. The particle-in-a-box basis functions used are as follows: Here, L denotes the width of the box, which is placed symmetrically such that x ranges over −L∕2 ≤ x ≤ L∕2 . Individual basis functions are labeled by the index, N = 1, 2, … . Note that upper case basis labels, N, are used to distinguish from the lower case n labels used for HO eigenstates and energy levels.
In terms of the above basis, the matrix representation of the kinetic energy of Eq (6) is diagonal, as discussed, with individual matrix elements, K NN , given by The corresponding potential energy matrix representation has diagonal elements, and off-diagonal elements, For ℏ = k = 1 , and a set of mass values lying in the interval 10 −6 ≤ m ≤ 10 3 , the eigenstates of the matrix H NN � = K NN � NN � + V NN � were computed using Mathematica [32] and a fortran code of the Jacobi eigenvalue algorithm [33]. For all reported energy levels, the box width L and the basis size N max were increased independently, until numerical differences were smaller than 10 −6 . Convergence criteria and tables of low-lying numerically converged energy levels, as used in the figures that follow, may be found in the Supplemental Material. Note that in the increasingly "anharmonic" relativistic limit, numerical convergence becomes increasingly difficult.
In order to simplify the exposition that follows, the exact numerical values obtained as described above are referred to as the "exact numerical" results, whereas those obtained with the formalism presented in Sect. 3 are called the "self-consistent GPPP" (SC-GPPP) values. However, note that the latter are not the exact GPPP solutions from the GPPP Eq. (13), but rather, those obtained using the approximate self-consistent procedure, computed using Eq. (27). Figure 2 shows the eleven lowest-lying energy levels as computed using the exact numerical method (open shapes)-in units of ℏ , and for four different m values. Also indicated are the corresponding SC-GPPP energies from Eq. (27) (solid shapes). Finally, the thin solid line indicates the non-relativistic or large-m limit.
Several important trends are evident. First, all relativistic energies lie below their non-relativistic counterparts, as they should. Second, these discrepancies increase as either: (a) the mass m decreases; (b) the excitation n increases. With regard to (b), this corresponds to a transition from a linear n to n 2∕3 dependence, vide Eq. (29) and the discussion just below. All of these trends are consistent with expectations.
In comparing the exact numerical and SC-GPPP energies, similar trends also hold, in the sense that the latter systematically overestimate the former-although these two relativistic predictions are generally in quite good agreement with one another, and certainly closer to each other than to the non-relativistic values. Also, the discrepancy trends (a) and (b) above also still hold true.
A similar comparison to that of the other two figures is also provided in Fig. 3. Here, we present absolute energies (left vertical axis) along with energies rescaled by ℏ , for the three lowest-lying energy levels. Thus, the figure allows a detailed comparison for the entire mass interval. First note that in units of ℏ the constancy and even spacing of the energy levels are lost when the oscillator enters the (ultra) relativistic domain. This is a consequence of the increasingly weak dependence of the energy on particle mass for m → 0 . In addition to exact numerical results, the SC-GPPP results from Eq. (27) are also presented in Fig. 3. These are found to be always an overestimate. Not surprisingly, the discrepancy increases with increasing n and decreasing m. What is quite remarkable, however, is what happens in the ultrarelativistic limit-i.e., energy levels represented by all sets of curves become completely independent of particle mass-exactly as predicted in Sect. 3. Also shown in Fig. 3 are the values from the WKB theory, computed using Eq. (36). The agreement between the WKB and exact results is remarkable and improves with increasing n, as expected. This effect can be seen more quantitatively in Table 2, showing a comparison of relativistic HO ground state energies, E 0 , for different mass values, m, as computed using different methods. Note that in addition to the exact numerical method, and the SC-GPPP method of Eq. (27), we also present exact GPPP results in this table, obtained using the new method described at the end of Sect. 2.1. These are the only exact GPPP results presented in this paper. Note that for both the GPPP and SC-GPPP calculations, rescaled mass values, 0 , are also readily available-from Eq. (15).
The table clearly indicates that for a given method and quantity, all values approach constants in the ultrarelativistic limit, m → 0 . The rescaled mass values are particularly illuminating. In every case, 0 > m , as required. However, for the larger two mass values, the relative increase is small, whereas in the ultrarelativistic limit, finite values on the order of 0 ≈ 0.2 are approached, even as m itself becomes arbitrarily small. Another interesting observation here is that the approximate SC-GPPP solutions are actually closer to the exact numerical results than are the exact GPPP results-at least for this ground state comparison. Moreover, the relative difference between the errors of the two approximations is great in the m → 0 limit, but all but vanishes for m → ∞-for reasons that currently elude us. Note that all of the methods described here allow for calculation of eigenstate wavefunctions as well as energy levels. Figure 4 indicates wavefunction probability densities for several different ℏ = k = m = 1 HO excited eigenstates, as computed for the non-relativistic, exact numerical relativistic, and SC-GPPP methods. The exact numerical and SC-GPPP eigenstates appear to be more localized in the potential well, due to the larger rescaled mass of the relativistic particle. As expected, this occurs to an extent that is dramatically increased with increasing n. In addition, keeping with SC-GPPP's "intermediate" role between non-relativistic and exact relativistic solutions (as also exhibited in the energy spectrum), we find that the SC-GPPP compression is less pronounced than that of the exact numerical solutions.
The shape of the density plots is also worth commenting on. Given that SC−GPPP is in fact equivalent to a mass-rescaled non-relativistic HO eigenstate, the form of its density is identical to that of the latter-except compressed. This can be seen clearly in the figure. On the other hand, density plots for the exact numerical solutions are qualitatively different, in that all of the interior peaks are of nearly equal height-as opposed to the familiar Hermite polynomial form, which drops down gently as one approaches the origin. This highly interesting feature of the exact relativistic case can be accounted for using the so-called "bipolar formalism," [34] as follows.
First, it is helpful to think of the actual eigenstate solutions as standing waves, each comprised of two equal and oppositely-moving traveling waves, ± , with = + + − . The traveling wave densities, 2 ± , are equal and smoothly varying, and correspond to the "envelope" of the standing wave. For example, if is  a sinusoidal particle-in-a-box state, then the corresponding ± are the associated right-and left-traveling plane waves, whose density is uniform throughout space. More generally, the envelope density function 2 ± (x) in the classically allowed region of space is nearly inversely proportional to the local classical velocity field, v(x). For the classical non-relativistic HO system, the kinetic energy-and therefore, velocity-is of course greatest at x = 0 , and then tapers off as one approaches the classical turning points. The classical relativistic HO system is similar-except that near x = 0 , the velocity can be much smaller than in the non-relativistic case, and is, of course, limited by a maximum finite value, c. Relativity thus serves to flatten the variation in the velocity field v(x)-and thereby, also, the corresponding envelope density function.

Conclusions
The quantum HO problem is one of a handful of fundamental applications taught in every quantum mechanics course. Yet, a practical relativistic generalization has seemingly eluded us for decades. According to conventional wisdom, one is forced into the dilemma of having to choose either second-order operators in time (e.g. Klein-Gordon) leading to negative energy solutions and far worse paradoxes, or square-root operators whose solution is deemed "impractical" except in the potential-free case.
The present study is part of an ongoing effort to break through this deadlock by dismantling the second tine of the above fork-i.e., by rendering treatment of the square-root operator practical, and even Schrödinger-like. Several exact and approximate methods have been proposed by us, with some of each implemented in the present work. Indeed, this study represents the first time that such exact relativistic methods have been applied to a system with an external potential, to our knowledge. That said, there is clearly a connection with earlier solution methods for the spinless Salpeter equation, widely used to explore model potentials of hadrons as bound state of quarks.
In particular, exact numerical solutions of the covariant stationary state equation of Eq. (6) were obtained for a relativistic particle operating in a quadratic HO potential. By working in a particle-in-a-box basis representation, the square root operator contribution poses no special numerical difficulty, and conventional non-relativistic numerical eigensolver methods may be employed.
The exact numerical relativistic solution eigenenergies are found to be lower than the corresponding non-relativistic energies, with the discrepancy increasing in the relativistic limit of decreasing mass, m, or increasing excitation n-with the latter necessarily implying a nonlinear spectrum. Thus, in units of ℏ the well known constancy and even spacing of the HO energy levels are no longer observed at the (ultra)relativistic domain. These trends can be accounted for using simple intuitive arguments. Other, much more intriguing trends-such as the m 0 n 2∕3 scaling of spectral frequencies that emerges in the m → 0 limit, and the uniformity of the relativistic density function peak heights-require more subtle but still largely straightforward analyses.
We have also employed two natural approximation methods in this work-e.g., GPPP, and SC-GPPP. In both cases, the methods are fully general, but we here exploit some simplifications that arise only in the HO case. Whereas SC-GPPP is actually a further approximation to GPPP, it counterintuitively provides the more accurate results, at least for the relativistic HO ground state energies. In any event, both approximations are quite good, and certainly a lot closer to the exact results throughout the very broad range of mass values considered ( 10 −6 ≤ m < 10 3 ) than are the non-relativistic results. Finally, the WKB semiclassical theory for the HO problem was extended to the (ultra)relativistic domain, and approximate and exact semiclassical expressions for the energy levels of a relativistic HO were derived. The energy levels computed with the relativistic WKB formulas are in very good agreement with the other methods here presented.
We authors feel that the results offered here-presenting a continuous unified description of the HO as it undergoes a nine-order-of-magnitude transition from the non-relativistic to ultrarelativistic regimes-may be of strong interest to a range of scientific researchers, for reasons both practical and fundamental. In addition to this-and of no less import-is the pedagogical value that this approach may also provide. It could be of great benefit to undergraduate and graduate physics students, struggling to understand even the basic consequences of merging Quantum Mechanics with the Special Theory of Relativity-as indeed, so many of their forebears before them have also struggled.