Next Article in Journal
Leggett–Garg-like Inequalities from a Correlation Matrix Construction
Previous Article in Journal
Defending Many Worlds via Case Discrimination: An Attempt to Showcase the Conceptual Incoherence of Anti-Realist Interpretations and Relational Quantum Mechanics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Laplace Method for Energy Eigenvalue Problems in Quantum Mechanics

1
Department of Physics, Georgetown University, 37th and O Sts. NW, Washington, DC 20057, USA
2
Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany
*
Author to whom correspondence should be addressed.
Quantum Rep. 2023, 5(2), 370-397; https://doi.org/10.3390/quantum5020024
Submission received: 5 February 2023 / Revised: 25 February 2023 / Accepted: 31 March 2023 / Published: 20 April 2023

Abstract

:
Quantum mechanics has about a dozen exactly solvable potentials. Normally, the time-independent Schrödinger equation for them is solved by using a generalized series solution for the bound states (using the Fröbenius method) and then an analytic continuation for the continuum states (if present). In this work, we present an alternative way to solve these problems, based on the Laplace method. This technique uses a similar procedure for the bound states and for the continuum states. It was originally used by Schrödinger when he solved the wave functions of hydrogen. Dirac advocated using this method too. We discuss why it is a powerful approach to solve all problems whose wave functions are represented in terms of confluent hypergeometric functions, especially for the continuum solutions, which can be determined by an easy-to-program contour integral.

1. Introduction

Quantum mechanics problems are typically solved via the differential equation form of the Schrödinger equation in position space and using the Fröbenius method for the generalized series solution. Unfortunately, this method of solution for bound-state problems cannot be directly used to solve continuum problems. In this work, we show how one can employ the Laplace method to solve these problems with contour integrals in the complex plane. This approach naturally allows one to solve both bound-state problems and continuum problems with similar efforts. It may also provide useful ways to experiment with and visualize the bound-state and continuum wave functions using the same numerical methods for both types of solutions. This method does have restrictions of its own, though. It can only solve problems whose wave functions can be represented in terms of confluent hypergeometric functions.
This approach is as old as quantum mechanics itself. It was introduced by Schrödinger when he determined the spectrum of hydrogen in his first wave mechanics paper in 1926 [1]; the method was based on Schlesinger’s famous differential equations textbook [2]. In 1937, Dirac advocated for using this same technique [3]. Oddly, he did not include it in his quantum mechanics textbook [4] (which was revised to its third edition ten years later), so it never was broadly adopted by the physics community. Modern texts that use it include both Landau and Lifshitz [5] and Messiah [6], where it is employed in the appendix to describe the properties of confluent hypergeometric functions and Konishi and Paffuti [7], where they employ it to solve the linear-potential problem, but not other problems solved by confluent hypergeometric functions. Capri [8] also uses it, but not as a general method. He suggests using an integral form (that corresponds to the Laplace method solution) to determine an integral representation for the bound-state wave functions of the Coulomb problem, and then evaluates them by residues.
The Laplace method for solving differential equations is summarized in our earlier paper on Schrödinger’s first solution for hydrogen [9]. We provide a briefer summary here of this method. The basic idea is to find an integrating factor for the differential equation and then solve it using a contour integral constructed from this integrating factor.
The Laplace method works for arbitrary order linear differential equations for Φ ( ξ ) that have constant coefficients a m or linear coefficients b m (in the dependent variable ξ ). These differential equations take the form given by
m ( a m + b m ξ ) Φ ( m ) ( ξ ) = 0 ,
where the ( m ) superscript denotes the mth derivative of the function Φ ( ξ ) with respect to ξ . For quantum mechanics solutions, we concentrate on the m = 2 case because the time-independent Schrödinger equation is a second-order differential equation. The solution to Equation (1) is constructed by introducing integrating factors and is represented in the form
Φ ( ξ ) = γ e ξ z R ( z ) d z ,
with the integral being over a properly chosen contour γ in the complex z plane. The function R ( z ) is determined from the integrating factor mentioned above.
The form of the ansatz for Φ ( ξ ) allows us to compute derivatives by differentiating under the integral sign,
Φ m ( ξ ) = γ e ξ z z m R ( z ) d z
for properly chosen contours, where this procedure is well defined. Plugging this representation into the differential equation yields
γ e ξ z m ( a m + b m ξ ) z m R ( z ) d z = 0 ,
which motivates the definition of two polynomials P ( z ) and Q ( z ) via
P ( z ) = m a m z m and Q ( z ) = m b m z m
that convert the differential equation into
γ e ξ z P ( z ) + Q ( z ) ξ R ( z ) d z = 0 .
If the integrand is the derivative of a complex-valued function that has the same value at the endpoints of the contour γ , then we can immediately obtain the solution to the differential equation. We require that R ( z ) satisfies the following differential equation
P ( z ) R ( z ) = d d z Q ( z ) R ( z ) ,
so that the integrand is a pure differential and that the function V ( z ) is defined by
V ( z ) = Q ( z ) R ( z ) e ξ z
has equal values at the endpoints of the contour (or is single-valued for a closed contour), so that evaluating the integral of the perfect differential yields 0 as required to solve the differential equation. Note that more than one contour can satisfy these conditions; indeed, for an order-m equation, we know that exactly m linearly independent solutions are represented by m inequivalent contours.
The function R ( z ) is then found by integrating Equation (7) after dividing it by Q ( z ) R ( z ) and recognizing the resulting logarithmic derivative. Solving, yields R ( z ) as
R ( z ) = 1 Q ( z ) exp z P ( z ) Q ( z ) d z .
Armed with R ( z ) , we immediately find the solution of the original differential equation to be
Φ ( ξ ) = γ e ξ z 1 Q ( z ) exp z P ( z ) Q ( z ) d z d z ,
where the contour γ needs to be chosen so that the vanishing endpoint condition is fulfilled. In quantum mechanics, we require the solution to satisfy additional properties, such as being always finite, or being square-integrable.
Next, we apply these methods to solve for the bound states of the simple harmonic oscillator in one-, two-, and three-dimensions. We do the same for the Coulomb problem in two- and three dimensions, as well as solving the Morse potential. We do not consider the spherical harmonics problem, or the Pöschl-Teller and Hulthén potentials, because all of these cases require full hypergeometric functions for the wave function. These can be solved by the closely related Laplace transform method, as discussed in the work of Tsaur and Wang [10]. Others have used the Laplace transform method too [11,12,13].

2. Bound States with the Laplace Method

We start with the application of the Laplace method to the bound-state problem. As mentioned above, we examine the simple harmonic oscillator in one-, two- and three-dimensions; the one-dimensional case is treated twice because there are two different forms for the wave function ansatz that are commonly used. We also cover the inverse r potential (Coulomb problem) in two- and three-dimensions. Finally, we discuss the Morse potential.
The potentials for these different problems are summarized in Table 1, where we use μ to denote the mass of the particle (sometimes this is the effective mass of a two-body problem), ω is the oscillator frequency, e is the magnitude of the charge of an electron and of a proton, V 0 > 0 is an energy scale for the Morse potential, and a > 0 is an inverse length for the Morse potential. Since all of these potentials do not have a Schrödinger equation that is in the Laplace form, we must do two additional things to arrive at the Laplace form: first, we construct an ansatz for the wave function, and compute the differential equation for the unknown function in the wave function ansatz; and second, we use a dimensionless independent variable (which sometimes is related to the original variable by a change in functional form). These choices are also summarized in Table 1. The wave function ansatz arises from a number of different strategies. For the one-dimensional simple harmonic oscillator, we have a different form for the even and the odd solutions, in higher dimensions, we separate out the angular and radial degrees of freedom—in two dimensions we use ρ and ϕ for the polar coordinates, while in three-dimensions we use r for the radial coordinate and θ (polar angle to the z-axis) and ϕ (azimuthal angle in the x-y-plane) for the angular coordinates. The radial functions also have a power-law behavior as they approach the origin included in the ansatz. The second method for the simple harmonic oscillator in one dimension removes a Gaussian factor from the wave function. The independent variable (always dimensionless and denoted by ξ ) is proportional to the square of the (radial) coordinate for the simple harmonic oscillators, is linear for the Coulomb problems, and is exponential for the Morse potential. For the second way, we treat the simple harmonic oscillator in one dimension as linear.
Note that the ansatz we use is not the standard one used in textbooks, where one incorporates the asymptotic behavior for large arguments of the wave function as well. Instead, we are following the methodology of Schrödinger in his original paper, where only the behavior at the origin is incorporated into the ansatz.
Plugging the wave function ansatz and the functional form of the independent variable into the Schrödinger equation, yields a differential equation in the Laplace form for all of the examples we are working on in this paper. The results of this exercise are summarized in Table 2. There, you can see that all but the last row (which will be treated separately) have Schrödinger equations that have been transformed into the form:
ξ Φ ( ξ ) + β Φ ( ξ ) + δ λ 2 ξ Φ ( ξ ) = 0 ,
where β , δ , λ R with λ = 1 or 1 2 > 0 in the bound-state cases we consider here (when we later consider continuum solutions, the transformed Schrödinger equation for Φ has the same form, but λ becomes an imaginary number, as we discuss below). Notably, the form of Equation (11) is the same form used to treat the three-dimensional hydrogen atom, as discussed in Refs. [1,9]. Nevertheless, the differences in the parameters in the Laplace form of the final differential equation, along with the different wave function ansatzes, mean that the analysis can be slightly different for some of these cases. We summarize the procedure in as general terms as possible, and show where the cases vary, as needed. Note further that this is not the standard Kummer equation, which arises when one includes the asymptotic behavior of the wave function for a large argument as well. The form we use is in many respects easier to derive and simpler to work with. Of course, we will see the final answers are the same as one obtains with the standard methodology, as must be so.
As an illustration of this procedure, consider the even solutions of the one-dimensional simple-harmonic oscillator, as given in the first row. We use ξ = μ ω x 2 / , so the kinetic energy operator becomes
2 2 μ d 2 d x 2 2 ω ξ d 2 d ξ 2 + 1 2 d d ξ
and the time-independent Schrödinger equation transforms into the equation in the top row of Table 2 after we divide both sides by 2 ω and we put all terms on the left-hand side of the equation. The other rows are derived similarly.
To begin the Laplace method, we construct the P and Q polynomials as described in the introduction. They become
P ( z ) = β z + δ ,
Q ( z ) = z 2 λ 2 = ( z λ ) ( z + λ ) .
The ratio of P / Q is then
P ( z ) Q ( z ) = α + z λ + α z + λ ,
where
α + = β λ + δ 2 λ and α = β λ δ 2 λ .
This can be immediately integrated and exponentiated to form the integrand for the contour-integral form of the solution, which, up to a constant prefactor, is given by
Φ ( ξ ) = γ d z e ξ z ( z λ ) α + 1 ( z + λ ) α 1 ,
where we still need to choose the contour γ . We must choose the contour γ so that the value of the integrand of Equation (17) multiplied by Q ( z ) = ( z λ ) ( z + λ ) , has equal values at the endpoints. Closed contours will always satisfy this. If the contour is not closed, then the function that must be equal at the endpoints is
V ( z ) = e ξ z Q ( z ) R ( z ) = e ξ z ( z λ ) α + ( z + λ ) α
and this must hold for all ξ . Note that we always have α + > 0 and α + + α = β > 0 , but α can take on positive or negative values. In general, it is difficult to find contours where this function will be equal at the endpoints (for arbitrary values of V) if the path is not a closed path. However, special values, such as 0 or , are easier to enforce, because they automatically hold for all ξ ; for V = one must use a careful limiting procedure to ensure the difference of V at the two endpoints actually vanishes.
The selection of these contours, and examining the integral over inequivalent contours is at the heart of the Laplace method. It requires a discussion of some technical details from complex analysis. Specifically, the integrand is generically multivalued when α ± are non-integers. Consider a complex-valued function of the form z α . We can write this as
z α = e α log ( z ) .
The complex logarithm is multivalued in the complex plane; consider
log ( z ) = log | z | e i θ = log ( | z | ) + i θ ,
with z = | z | e i θ written in polar form. Clearly, if we fix | z | and move through θ from 0 2 π , we will return to the same point on the complex plane, but the value of the logarithm will have picked up an additive shift by 2 π i and thereby will not be a continuous (or even well-defined) function. So, we instead write the logarithm function as log ( z ) = Log ( z ) + 2 π i m , m Z , where Log indicates the principal value of the logarithm function, which we take to be the complex logarithm function for 0 < θ < 2 π in this example; note that this choice is a domain that has a cut in the complex plane running from 0 to infinity along the positive real axis. Log ( z ) is then single-valued, but the domain is smaller than the full complex plane due to the cut. This part of the complex plane given by 0 < θ < 2 π is a branch of the logarithm, defined to be a domain in the complex plane on which a function is single-valued. A branch is chosen to be maximal, in the sense that no point can be added to it and still maintain single-valuedness.
The boundary of a branch is called the branch cut. From our definition of the complex logarithm, it is clear that the branch cut will be the ray where θ = 0 (or more completely, where θ = ( 2 n π for n Z ) originating at the origin; that is, the positive real axis plus the origin. In this case, the origin is a branch point of the logarithm function since it is a point that is common to any branch cut one can draw for the logarithm. The branch cut can be any curve that extends from the origin to infinity to have a single-valued logarithm function. Note that the complex plane, C , does not include infinity, and as such there is no “point at infinity;" instead approaching infinity, means that the complex number continues to increase without bound ( | z | ).
Critical to our work is a procedure that determines the phase of a number z in a domain on the complex plane that has a branch cut. In the preceding discussion, the phase θ on the complex plane was defined in the standard convention: the phase is calculated with respect to the origin. That is, to determine the phase of z, we set the origin as the anchor of a vector of length | z | starting from the horizontal (and oriented along the positive real axis), and rotate it counterclockwise until we reach the point z. Then, the phase of z is the angle swept out by this line. When one needs the phase of a point z relative to a reference point ( z r ) , in a domain with a branch cut, one generalizes the process. We draw a path from the reference point to z that does not cross any branch cut. Such a path will always exist for the cases we consider here. Then, we track the phase, relative to the reference point, as we move along the path, until we reach z. The phase is tracked by observing how the angle between the reference point and a line whose endpoint traverses along the path varies from the reference point to z; explicit examples of how this works will be given below. Note that while the path (and hence the tip of the line) never crosses a branch cut, the arrow drawn from the reference point to the tip often does cross a branch cut. Hence, the final phase that is used, may look like we rotated an arrow from the reference point to the final point z across a branch cut. But, as we see above, because we follow a path that does not cross the branch cut, this is fine.
We now return to the complex power function, which becomes
z α = e α Log ( z ) e 2 π i m α .
Clearly, if α Z , then e α Log ( z ) = z α , which is single-valued. Otherwise, it is multivalued. In general, α ± in Equation (17) are not integers. This means that the integrand is not single-valued in general. We must draw branch cuts to restrict the complex plane to a domain where the integrand is single-valued. The contour is then required not to cross any branch cut (but it can have endpoints at branch cuts or branch points).
For a given multivalued function, there can be many inequivalent ways to draw branch cuts (the only requirement is that the integrand is single-valued within the domain that has the branch cuts drawn). From the form of the integrand in Equation (17), we know that the branch points will be determined by the exponents α ± ; that is, ± λ is a branch point if and only if α ± Z , respectively. By construction, α + + α = β , so looking at Table 2, we immediately see two cases: those where β Z always (rows 3, 5, and 6) and those where β Z (all other cases); note the Morse potential can be an integer in special cases. These two cases have different allowed branch cuts. In addition to determining the branch cuts and the allowed contours (which must have V ( z ) be equal at the endpoints for all ξ values), we also require the function Φ ( ξ ) to be nonsingular for all ξ . This condition is different from the conventional requirement of square-integrability of the wave function, but it is the condition needed for the Laplace method. We will find all bound-state wave functions that satisfy this condition are also square-integrable, but we do not prove this.
First, we can eliminate any closed contour that does not cross a branch cut or contain any branch points. This is because such a closed contour always gives a zero integral due to Cauchy’s theorem and the fact that the function is analytic inside the domain defined by the branch. An open contour that goes from a finite point z to another finite point z inside the domain (that does not cross a branch cut) is also ruled out because it will not be possible to satisfy V ( z ) = V ( z ) for all 0 ξ < .
These restrictions then imply that contours must go to infinity, or end at a branch point (technically they could also end anywhere on a branch cut, but similar to the open contour case above, one won’t be able to satisfy the condition V ( z ) = V ( z ) for all allowed ξ in that case). Rather than determine all possible contours next, we now consider the condition that the solution must be finite as ξ . Since there is a factor of e ξ z in the integrand, no contour can go to infinity with Re ( z ) > 0 without the integral diverging as ξ + . Moreover, the contour cannot have an endpoint at z = λ because using the stationary phase method to asymptotically determine the limit of the integral for large ξ , again has an exponential term that will diverge (for details see [9]).
For all cases, we can have a branch cut running from to λ along the negative real axis and from λ to along the positive real axis. When β is an integer, we can instead have the branch cut run along the real axis between λ and λ . These two possibilities are depicted in Figure 1.
We describe next the three possible contours we can have for all 0 < ξ (we will determine the behavior for ξ = 0 last). Two contours can be drawn when we have the branch cuts extending from λ to along the negative real axis and from λ to along the positive real axis. The first one, γ 1 , is the so-called Hankel contour, which starts at just below the real axis, loops around the branch point, and goes back to above the real axis (see Figure 1a). The second possibility, γ 2 , starts at the branch point at λ and goes to (either above or below the real axis, both end up giving the same asymptotic behavior as ξ 0 ; we draw it above the real axis here). The final contour can be drawn when the branch cut extends from λ to λ (which is only possible when α + + α = β is an integer). This contour, γ 3 , encircles the branch cut and the branch points; it is sometimes called the dog-bone contour. Note that a similar γ 2 exists when the branch cut runs from λ to λ , but it behaves the same as the γ 2 that we analyze. A stationary phase analysis of γ 3 shows that it diverges as ξ , so it is ruled out [9].
Now, we need to determine if any of the remaining contours yield a finite result as ξ 0 . First, in the case of the Hankel contour, the phase of the integrand will be different below and above the real axis, since we drew a branch cut along the negative real axis. Hence, Φ ( ξ ) for ξ = 0 for the Hankel contour, because the integrand remains finite but the contour extends an infinite length, and the difference in the phases means those infinite contributions do not cancel. Thus, we can rule out the γ 1 contour. For the contour γ 2 from λ (in the left half plane), if α < 0 , the integral will diverge as z λ . If α 0 , we find that the contour does satisfy the condition that V ( z ) has equal values (in fact V ( z ) = 0 ) at the endpoints of the contour γ 2 for ξ > 0 . But, when ξ = 0 , V ( λ ) = 0 because there is a factor of z + λ raised to a positive power (it is a finite number when α = 0 ), while for z , | V | goes like | z | β . So, our condition that the endpoints have equal values of V for all values of ξ does not hold for ξ = 0 , eliminating γ 2 . This eliminates all possible contours when α ± are not integers.
We still have to consider the case when β Z for cases where this occurs. This does not change any of the prior analyses; it simply restricts us to exclude the cases where the branch cut ran from λ to λ —hence, we have no dog-bone contour to consider when β is not an integer.
Since there are no valid solutions for the general case, we must re-evaluate our prior assumption, that in general α ± Z , to see if we can find another contour in that case. It is clear that we will require at least one of them to be an integer, in general. When β Z , there is no ambiguity, as in this case, if one is an integer, they both must be. When β Z , the analysis is more complex. Note that, as shown in Table 2, these cases are the one- and three-dimensional oscillators and the Morse potential, since a Z , in general. Recalling that in all cases, δ , λ > 0 , we know that α + > α . Moreover, the only way to obtain a new contour is to have one that surrounds a pole, that is a point z at which the integrand is proportional to ( z z ) N for N Z + . Since α + + α = β > 0 and α + > α , we can only have α be a non-positive integer (we can have α = 0 , because in Equation (17), the exponent of ( z + λ ) is ( α 1 ) , so α = 0 will result in a first-order simple pole). That is,
| α | = N , N = 0 , 1 , 2 . . .
We will find that α = 0 corresponds to the ground state, and α < 0 Z will yield higher order poles, corresponding to excited states.
Note that the constraint that α is a nonpositive integer yields the quantization condition for each bound state (see Table 3). Key to our analysis is that integer exponents lead to single-valued powers, so z = λ is no longer a branch point, and we can enclose it with a closed contour (Figure 2). Since this point is a pole of order | α | + 1 , this contour is non-trivial and distinct from the previously analyzed ones. Notably, the cases when β Z are only slightly different; z = λ is still a branch point, and we choose the branch cut to run from λ along the positive real axis, as shown in Figure 2b. Thus, all cases of the form in Equation (17) lead us to the same new contour, a closed contour about a pole of order N + 1 at z = λ . This yields the solution in the form of a closed contour integral (with the contour encircling the point z = λ once in the counter-clockwise direction, as shown in Figure 2)
Φ ( ξ ) = d z e ξ z ( z λ ) α + 1 ( z + λ ) α 1 ,
which we evaluate by residues:
Φ ( ξ ) = 2 π i | α | ! d | α | d z | α | e ξ z ( z λ ) α + 1 z = λ .
From the Rodrigues formula for the associated (alternatively, “generalized”) Laguerre polynomials L N ( b ) ( ξ ) [14],
L N ( b ) ( ξ ) = 1 N ! e ξ ξ b d N d x N e x x N + b x = ξ ,
we can see that the form of Equation (24) is similar. To establish the exact correspondence, we define u = ( z λ ) ξ (the point at which we evaluate the derivative, z = λ , corresponds to u = 2 λ ξ ) and the inverse is z = ( λ ξ u ) / ξ . Re-expressing the right-hand side as a function of u yields
Φ ( ξ ) = 2 π i | α | ! 1 | α | ξ | α | d | α | d u | α | e u e λ ξ 1 α + 1 u ξ α + 1 u = 2 λ ξ = 2 π i ( 1 ) β 1 | α | ! ξ β + 1 e λ ξ d | α | d u | α | e u u α + 1 u = 2 λ ξ = 2 π i ( 1 ) β 1 2 λ β + 1 | α | ! e λ ξ L | α | ( β 1 ) ( 2 λ ξ ) .
Hence, except for a constant prefactor, which will be set via normalization, we find that
Φ ( ξ ) e λ ξ L | α | ( β 1 ) 2 λ ξ .
Because the Laguerre polynomial of order N goes like ξ N as ξ ± , the asymptotic behavior of Φ ( ξ ) is dominated by the exponentially decaying term. Thus, this solution will be finite for all 0 ξ and is the wave function we sought to find.
The last detail we need to work out is how the integer from the quantization condition N = | α | , relates to the conventional principal quantum number n. In all harmonic oscillator cases, we have n = N , while the two-dimensional Coulomb problem has n | m | 1 = N and the three-dimensional Coulomb problem has n l 1 = N . These results are summarized, in terms of the principal quantum number, in Table 3.

3. Examples: Bound States for the Simple Harmonic Oscillator

Having now constructed the general methodology to solve bound-state problems using the Laplace method, we show the concrete details of how the method is used for the two-dimensional simple harmonic oscillator, which we treat in polar coordinates. The Schrödinger equation is solved first by separating variables, that is by letting ψ ( ρ , ϕ ) = R ( ρ ) Φ ¯ ( ϕ ) . The solution for the angular function is
Φ ¯ ( ϕ ) = e i m ϕ , m = ± 1 , ± 2 , . . .
as summarized in Table 1. The angular function Φ ¯ should not be confused with the general function Φ used in our discussion of the Laplace method. From this, we obtain a radial equation:
R ( ρ ) + 1 ρ R ( ρ ) + 2 μ E 2 μ 2 ω 2 2 ρ 2 m 2 ρ 2 .
Now, we define the new independent variable ξ = μ ω ρ 2 , and make the ansatz that R ( ξ ) = ξ | m | 2 Φ ( ξ ) (also summarized in Table 1). Making these substitutions we find the Laplace form of the radial equation as
ξ Φ ( ξ ) + ( | m | + 1 ) Φ ( ξ ) + E 2 ω 1 4 ξ Φ ( ξ ) = 0 ,
as summarized in Table 2. Hence,
β = | m | + 1 , δ = E 2 ω , and λ = 1 2 .
Next, we construct
α + = 1 2 | m | + 1 + E ω and α = 1 2 | m | + 1 E ω .
The quantization condition becomes
N = | α | = n E n = ω 2 n + | m | + 1 .
From Equations (27) and (32), we find the desired solution to the differential equation is
Φ ( ξ ) = e ξ 2 L n | m | ξ
as summarized in Table 3. This then yields the following for the full wave function:
ψ n , m ( ρ , ϕ ) = ρ | m | e μ ω 2 ρ 2 L n | m | μ ω ρ 2 e i m ϕ ,
up to a normalization constant, which has not yet been determined.
The solutions we derived for the one-dimensional simple harmonic oscillator may not look familiar to many. But they actually are the standard form, expressed in terms of a Gaussian multiplied by a Hermite polynomial, once we realize that there is an identity relating associated Laguerre polynomials to Hermite polynomials, given by
H 2 n ( x ) = 4 n n ! L n 1 2 x 2 H 2 n + 1 ( x ) = 2 4 n n ! x L n 1 2 x 2 ,
where H n is the physicist’s form for the Hermite polynomial defined by Rodrigues’ formula
H n ( x ) = 1 n e x 2 d n d x n e x 2 .
The remaining cases follow in a similar fashion and can be constructed by following the summarizing formulas in Table 1, Table 2 and Table 3.
What remains is for us to discuss the second way to solve the one-dimensional harmonic oscillator, which obtains the Hermite polynomials directly from Rodrigues’ formula after determining the residue at the pole. This method is summarized in the last row of Table 1, Table 2 and Table 3 and it uses a different ansatz for the wave function, which also leads to a differential equation in the Laplace form, but a different one from all of the other cases. It is given by
Φ ξ 2 ξ Φ ξ 2 α Φ ξ = 0 ,
where α = 1 2 E ω ; note we must have E 0 (since the Hamiltonian is a positive semidefinite operator in this case), so α 1 2 . In addition, we need to find a finite solution for all < ξ < . We now go through the steps of the Laplace method. First, compute the polynomials
P ( z ) = z 2 2 α
Q ( z ) = 2 z ,
and then exponentiate the antiderivative of the ratio P / Q , which now includes only one power of z raised to a potentially noninteger exponent. This allows us to write the solution to Equation (38) as the contour integral
Φ ( ξ ) = C d z e ξ z z 2 4 z α 1 .
Next, we restrict the exponent of z such that α Z . The branch points of the integrand in Equation (41) are at z = 0 and at infinity. To construct a single branch of the integrand, we choose to draw the branch cut from z = 0 z = along the negative real axis.
We next consider the possible contours as shown in Figure 3. By the same arguments used before, we can eliminate all closed contours that don’t contain a branch point or end at a branch cut and all open contours between any two finite points in the complex plane. Thus, the only candidate contours are ones with at least one endpoint at z = 0 or that have endpoints that go to infinity. As the contour goes to infinity, it must remain inside a cone with an angle of ± π / 4 of the real axis, to be bounded. There is a cone around the positive real axis and also one around the negative real axis.
Thus, the contours to analyze are a Hankel-like contour around the branch cut, which we call γ 4 , a contour from negative infinity to positive infinity (which does not cross the negative real axis), called γ 5 , and a contour from the origin going to infinity with Re ( z ) < 0 , called γ 6 , or from the origin to infinity with Re ( z ) > 0 , called γ 7 ; as long as we remain inside the cones about the real axis as we go to infinity. It is fairly easy to see any other contour with endpoints at 0 and can be deformed into one of these contours or can be mapped to them by taking z z . For example, the integral over γ 7 is converted into a contour from the branch point at 0 that runs to infinity inside the white cone along the positive real axis by transforming z z , which is equivalent (up to a complex phase) to the integral over γ 7 but with ξ ξ .
The first condition we have is that the function V ( z ) has identical values at the endpoints for all ξ . In this case, we have
V ( z ) = e ξ z z 2 4 z α .
It is easy to see that lim z 0 V ( z ) = if α < 0 , while lim z 0 V ( z ) = 0 if α > 0 . The asymptotic behavior of V is dominated by the e z 2 term, when | z | , so we have lim z V ( z ) = 0 when we lie inside the white cone about the real axes. This implies, that we can only have z = 0 as an endpoint, when α > 0 . But, in that case, any integral that has 0 as an endpoint diverges, due to the power-law behavior of the integrand near z = 0 having too negative of an exponent. So no contour can have z = 0 as an endpoint, eliminating γ 6 and γ 7 .
The next condition is to check the integral as ξ . Since the contour lies inside the cone, it seems like it will always be bounded, and should never diverge. But, in the region near z 2 ξ , the integrand actually behaves like e ξ 2 and can give large contributions. The way to evaluate such situations is to perform an asymptotic analysis on the integrals to determine their value for large | ξ | . The standard way to do this is called the steepest-descents approach, which notes that the contributions to the exponential are largest near the maximum of the exponent, which is then described by a simple quadratic near the extremum, yielding a Gaussian integral that can be evaluated exactly. The analysis is straightforward and we describe it carefully, starting with the Hankel contour γ 4 .
The saddle-point approximation allows us to approximate integrals of the form
Φ ( ξ ) = γ d z e h ( z ; ξ ) g ( z ) ,
which corresponds to the integral we need to evaluate for the Laplace method solution. Note that we have two options for how to proceed here. We can pick g ( z ) = z α 1 , or we can pick g ( z ) = 1 via writing z α 1 = e ( α 1 ) ln z and absorbing this term into the exponent h ( z ; ξ ) . Both approaches yield the same final results, but the former is simpler than the latter, because there is only one saddle point in this case, making the analysis easier. This means we have h ( z ; ξ ) = ξ z 1 4 z 2 . Taking the derivative, to find the extrema, we have h ( z ; ξ ) = ξ 1 2 z . Setting this equal to zero, tells us the extremum occurs at z = z 0 = 2 ξ , which is called the saddle point; note that we have h ( z 0 ; ξ ) = ξ 2 . In complex analysis, one direction through the saddle point is a minimum, while the other is a maximum, yielding a saddle-point shape for the exponent near the saddle point z 0 . One must choose the direction for the contour through the saddle point to traverse the maximum, not the minimum. In this case, the maximum direction is along the real axis, which is simple to see, because we have a quadratic for the exponent with a negative real coefficient.
The asymptotic analysis next deforms the contour to go through the saddle point along the maximum direction. When ξ < 0 , this saddle point lies on the branch cut, so we will deform the Hankel contour to pass infinitesimally below it and parallel to the real axis—once below the negative real axis and once above. This yields two contributions for the steepest-descents integral. For the contribution from below the real axis, we parameterize the contour as given by γ z 0 + t near the saddle point, so that
h ( z 0 + t ; ξ ) h ( z 0 ; ξ ) 1 4 t 2 .
Since the integrand decays quickly away from the saddle point, we extend the limits on t to run from to and we approximate g ( z ) g ( z 0 ) . This then gives us the contribution from the saddle point below the real axis to be
d t e h ( z 0 ; ξ ) 1 4 t 2 g ( z 0 ) 2 π ( 2 | ξ | ) α 1 e i π ( α 1 ) e ξ 2 .
The contribution from the saddle point above the negative real axis is similar—it has an overall negative sign, because the contour runs from right to left instead of left to right and the sign of the phase in the exponent is positive because we are above the branch point. The total asymptotic estimate for the integral is then
γ 4 d z e ξ z 1 4 z 2 z α 1 2 α + 1 i π | ξ | α 1 sin π ( α 1 ) e ξ 2 .
Since we assume α is not an integer, the coefficient is nonzero and this gives a leading contribution that goes like e ξ 2 . Looking at Table 1, we see that the full wave function is proportional to e 1 2 ξ 2 Φ , so this solution will go as e 1 2 ξ 2 as ξ , which diverges. So, γ 4 does not yield a finite wave function. Moreover, a similar analysis yields the same asymptotic behavior for γ 5 . Note that this analysis is similar to the Fröbenius analysis when the series does not truncate, and we reject the solution due to the wave function growing as we go to infinity.
So we do not obtain a finite solution from any of the possible contours, and we again must change our assumption that α Z to allow for a new contour. When α Z , the integrand is single-valued, so there no longer is a branch point or branch cut. Consequentially, we can now enclose the origin with a new closed contour. Our only choice is then that α 0 , yielding the quantization condition given by
n = α = E ω 1 2 E n = ω n + 1 2 , n = 0 , 1 , 2 . . .
This determines the energy levels of the one-dimensional oscillator. Now, we write the (unnormalized) solution to the differential equation for this closed contour as
Φ ( ξ ) = d z e ξ z z 2 4 z n 1 , n = 0 , 1 , 2 ,
with a closed contour that encircles the origin. By completing the square in the exponential term in Equation (48), we can re-write this integral as
Φ ( ξ ) = e ξ 2 d z e ξ z 2 2 z n + 1 .
Now, we let u = ξ z 2 , and by making this substitution, we obtain, up to a constant prefactor,
Φ ( ξ ) = e ξ 2 d u 1 n e u 2 ( u ξ ) n + 1 .
This integral can be evaluated by residues about the pole of order n + 1 at u = ξ :
Φ ( ξ ) = 1 n e ξ 2 lim u ξ d n d u n e u 2 = 1 n e ξ 2 d n d ξ n e ξ 2 .
This is precisely the Rodrigues formula for the n th degree Hermite polynomial H n ( ξ ) . That is, up to a constant prefactor,
Φ ( ξ ) = H n ( ξ ) = H n μ ω x ,
which allows us to write the (unnormalized) wave function for the 1D harmonic oscillator as:
ψ ( x ) e μ ω 2 x 2 H n μ ω x .

4. Continuum Solutions with the Laplace Method

There are several quantum systems whose energy eigenstates have energy eigenvalues that lie in the continuum and that we can also treat with the Laplace method; this includes the free particle in two, and three dimensions, the continuum solutions of the Coulomb problem in two- and three-dimensions, and the continuum solutions of the one-dimensional Morse potential. The steps for obtaining a differential equation in the Laplace form are similar to what we already showed above, and in Table 4 we summarize the results for these different models. Note that the substitutions required for the free-particle problems look like the final wave function will diverge at the origin, but we require the function Φ ( ξ ) to have a high-enough order zero at the origin that the final wave function remains finite everywhere. In Table 5, we show the final differential equations obtained by this procedure, which is similar to the bound-state form in Equation (11), but with the sign of the λ ¯ 2 term changed. We have written that term as ( i λ ¯ ) 2 instead of + λ ¯ 2 to simplify the notation that we need in solving the problem. The Morse potential, on the other hand, keeps the form of Equation (11), but some parameters now become complex.
We begin with the free particle and Coulomb problems, each treated in both two and three dimensions, because they all are treated similarly. The differential equation in the Laplace form takes the form
ξ Φ ( ξ ) + β Φ ( ξ ) + δ ( i λ ¯ ) 2 ξ Φ ( ξ ) = 0 ,
for all of these cases. This means that we can write the solution as
Φ ( ξ ) = γ d z e ξ z z i λ ¯ α + 1 z + i λ ¯ α 1 ,
by following the Laplace method. Here, we have
α ± = i β λ ¯ ± δ 2 i λ ¯ .
In each of these four cases, λ ¯ = 1 . For the free-particle problems, we have δ = 0 , so there is only one α .
As with bound states, we restrict the contours over which we integrate by requiring that the wave function be finite everywhere. As discussed in [9], this constraint yields a rotated dog-bone shaped contour which encloses the branch points ± i , as seen in Figure 4. The reason why is that any contour that runs out to infinity will yield an infinite result for Φ ( 0 ) . We do, however, mention that the case of a 3d free particle actually does not need the branch cut, because both α ± are positive integers, and the integrand is not multivalued. In this case, the integral must be taken along the imaginary axis from i to i, connecting the two points where V = 0 for all ξ . Note that this choice could also have been used for the Coulomb potential cases (because V = 0 at the branch points for all ξ as well), but there is value to using the dog-bone contour because it allows us some additional options for evaluating Φ ( ξ ) numerically, as we discuss below; since the solutions corresponding to both choices of the contour are proportional to each other (shown below), we can freely choose either one. The ambiguity in the prefactor is always removed when the final wave function is normalized (but we will not discuss normalization in this work).
It is critical to evaluate the phases properly when determining the integrand. With the branch cut structure we use, the function is single-valued everywhere in the complex plane, except along the branch cut itself, but the function has different values on both sides of the branch cut. We start by picking a reference point, which will be the origin, just to the right of the branch cut, which we call z = 0 + , as shown in Figure 5a. The multivalued function in the integrand is f ( z ) = ( z i ) α + 1 ( z + i ) α 1 and we focus on how to determine the phase consistently for this function. We find that f ( 0 + ) = ( i ) α + 1 ( i ) α 1 , which we evaluate with the standard phases i = e i π 2 and ( i ) = e i π 2 . Then we have f ( 0 + ) = exp i π 2 ( α α + ) . Noting the form of α ± , we find f ( 0 + ) = exp ( π 2 δ ) .
Now, to calculate f ( z ) anywhere in the complex plane, we first draw a path from the reference point 0 + to z that does not cross the branch cut (see Figure 5b,c). Then, we examine how arrows drawn from the upper branch point i to a point on the path rotate as we move along the path from 0 + to z. This determines the change in the phase for the factor ( z i ) . We repeat with an arrow drawn from i to a point along the path, and follow it from 0 + to z. The rotation of the arrow here, again determines the change in the phase for the factor ( z + i ) . Each of those factors will have the change in phase multiplied by the corresponding exponent, and that will determine the phase of f ( z ) . We will show below that using this procedure produces a single-valued function over the entire complex plane. But first, we will use this procedure to convert the integral form of our solution to a single integral that runs along the real axis.
We deform the contour in Figure 4 to be infinitesimally close to the imaginary axis along the vertical portions, and wrapping infinitesimally close to the branch points in the circular portions, we obtain a result expressed as the sum of the contributions from two vertical lines (downward from i i just to the left of the axis and upward from i i to the right of the axis) and from the infinitesimal circular arcs enclosing the two branch points. It is easy to see that the contributions from the circular arcs around the branch points will be zero. Since Re ( α ± ) > 0 , the integral around each arc will go to zero, so the branch points do not contribute to the integral. The integral over the full contour then reduces to the sum of the contributions from the two vertical lines. Note that since they lie on either side of the branch cut, there is a phase difference between the two integrals, and the contributions will not cancel.
To determine this phase difference, we follow the procedure described above (see Figure 5 for a graphical representation). We start with the piece of the contour that runs from i to i along the right-hand side of the imaginary axis. We can draw a path from the reference point at 0 + to any point along this piece of the deformed contour as a straight vertical line. We can immediately see that the arrow from i to z on the contour does not change direction after moving along the path. Neither does the arrow drawn from i to z. This means the phase for f ( z ) is the same as the phase for f ( 0 + ) along the entire piece of the contour. We write z = i y , with y real, and this part of the contour integral becomes the following:
I 1 = 1 1 i d y e i y ξ | y 1 | α + 1 | y + 1 | α 1 e π 2 δ .
The second piece of the contour runs along the left-hand side of the imaginary axis from i to i . We draw a path from the reference point 0 + to z along the contour, by an arc that goes around the upper branch point at i. One can immediately see that an arrow drawn from i to the reference point 0 + , winds by 2 π in the counterclockwise direction as it goes around the path to z. So the change in the phase for the factor ( z i ) is 2 π . The arrow drawn from i to the reference point, will rotate first to the right, and then to the left, but it ultimately accumulates no net phase, so its change in phase is 0. The fact that the final point ended on the left-hand side of the branch point does not determine the change in the phase along the path, it is the net motion of the arrow as we follow along the path that does this. This point can often be misunderstood. It is important to note that our rule for determining the phase of the function does not input by hand a change of phase when crossing the branch cut. Instead, we follow the described procedure to determine the change in the phases. Again, we let z = i y and we note that the integral runs down the imaginary axis, so we find the contribution from this piece of the contour is the following:
I 2 = 1 1 i d y e i ξ y | y 1 | α + 1 | y + 1 | α 1 e 2 i π ( α + 1 ) e π 2 δ .
The factor e 2 i π = 1 can be ignored. The real part of α + is either an integer or a half-odd integer. In the former case, the added factor is 1 and can be ignored, in the latter case it is 1 . The imaginary part of α + adds in a factor of e π δ . Hence, when we combine the two integrals together (and recall that we have to switch the order of the limits in the second integral) we find the total contour integral becomes
Φ ( ξ ) = i e π 2 δ e π 2 δ 1 1 d y e i ξ y | 1 y | α + 1 | 1 + y | α 1 ,
where the minus sign is for when the real part of α + is an integer and the plus sign is for when it is a half-odd integer.
To convert this into a form that is easily expressed in terms of confluent hypergeometric functions, we let y = 1 + 2 x and substitute into the integral to find that
Φ ( ξ ) = i e π 2 δ e π 2 δ 2 α + + α 1 e i ξ 0 1 d x e 2 i ξ x | 1 x | α + 1 x α 1 .
Comparing to the standard integral form of the Kummer function (for Re ( b ) > Re ( a ) > 0 ) (as given by Equation (13).4.1 of Ref. [15])
M ( a , b , z ) = Γ ( b ) Γ ( a ) Γ ( b a ) 0 1 d x e z x x a 1 ( 1 x ) b a 1 ,
with Γ ( x ) the Gamma function, we find that
Φ ( ξ ) = i e π 2 δ e π 2 δ 2 α + + α 1 Γ ( α + ) Γ ( α ) Γ ( α + + α ) e i ξ M ( α , α + + α , 2 i ξ ) .
Note that for E > 0 , the numerical prefactor is never zero, so we can always remove it from further discussion in the summary of the wave functions. It will enter, and is important, when we evaluate the results numerically below. One does need to complete the normalization step for the final wave functions (which is usually done with delta-function normalization), but we will not discuss that further here and instead will only summarize unnormalized wave functions, with the prefactor removed. Note that the result here corrects a sign error in the final result for the continuum wave function in Ref. [9] arising from an inconsistent definition of the phase of the multivalued function.
We comment briefly here on the 3d free particle. This case results in just the I 1 term, which has the same final form as we have for the cases with the “dog-bone” contour (just with a different prefactor). So its result falls into the same category as the other three cases. It is just that we do not need to worry about any phase issues in working with the integrand in this case. Because there is no branch point, we cannot describe that solution with a closed contour, because such an integral always vanishes. Instead, we simply integrate from one “zero” point to the other.
One other point to note is the Kummer relation (Equation 13.2.39 of Ref. [15])
M ( a , b , z ) = e z M ( b a , b , z ) ,
we can show that the unnormalized Φ ( ξ ) , given by e i ξ M ( α , α + + α , 2 i ξ ) , is real for real ξ . In particular, we have
e i ξ M ( α , α + + α , 2 i ξ ) = e i ξ M ( α , α + + α , 2 i ξ ) ,
where we used the facts that α + + α and ξ are both real. Then, if you note that α = α + = ( α + + α ) α , and use the Kummer relation, the right hand side of the equation becomes e i ξ e 2 i ξ M ( α , α + + α , 2 i ξ ) , which is equal to the original function and hence shows that this combination is real. Thus, if we drop the constant prefactor, the unnormalized Φ ( ξ ) can always be chosen to be real-valued. Note as well that at ξ = 0 , we have Φ ( ξ ) = 1 , because M ( a , b , 0 ) = 1 for all cases where the Kummer function is well defined from its power-series, which is the situation we have here. So this choice of contour leads to the correct continuum wave function, as we claimed earlier. A summary of the unnormalized wave functions for all potentials (which have real radial functions) is given in Table 6. These results follow by simply plugging in the explicit values of α ± and noting that there are identities between confluent hypergeometric functions and other functions, such as Bessel functions and spherical Bessel functions (as summarized in Section 13.6 of the NIST Digital Library of Mathematical Functions [15]). We do not show the details for how to carry out that algebra here.
One of the benefits of using contour-integral representations for the continuum wave functions is that it allows us to explore different ways to determine the wave functions numerically. For example, in this case we have three equivalent numerical representations. The first involves a real integral and is given in Equation (62). The second involves our original integral representation in Equation (54), where, for concreteness, we use a circular contour of radius R centered at the origin for evaluating the wave function. The third is to develop a power series representation and then to numerically evaluate the series. This is done by deforming the contour until it has a very large radius, and then extracting the residue at infinity.
We describe how to determine this power series next. The key to completing the calculation is the determination of the Laurent series near the point at infinity. This is most easily determined by approaching infinity along the positive imaginary axis; we add a constant shift by i first, because we know the final result has a factor of e i ξ . So, we let z = i + i / y for y a positive real number near zero. The phase for the factor z i is i π , because the arrow from the reference point wraps by π as we move up the imaginary axis, while the phase from z + i is zero. This gives us an overall phase factor of e i π α + . The integrand (including the change of variables factor i / y 2 ) is then given by
i e i π α + π 2 δ e i ξ e i ξ y 1 y α + + α | 1 2 y | α + 1
(the term e π 2 δ comes from the phase of the function at our reference point). We now want to expand this in a Laurent series for small y. Using the generalized binomial theorem for complex powers, yields
i e i π α + π 2 δ e i ξ m = 0 ( i ξ ) m m ! j = 0 ( α + 1 ) j j ! ( 2 ) j y j m α + α ,
where ( α 1 ) j is the Pochammer symbol for the falling factorial, given by ( α 1 ) ( α 2 ) ( α j ) . If we now perform a contour integral around the point at infinity, the result will be given by the residue, which is determined by the coefficient of the expansion in Equation (66) of the term 1 / y . Since α + + α is an integer, we can immediately determine the residue. It is given by
Residue = i e i π α + π 2 δ e i ξ j = α + + α 1 ( α + 1 ) j j ! ( 2 ) j ( i ξ ) j + 1 α + α ( j + 1 α + α ) ! .
Shifting the summation to start from zero, then gives
Residue = i e i π α + π 2 δ e i ξ ( 2 ) α + + α 1 j = 0 ( α + 1 ) j + α + + α 1 ( j + α + + α 1 ) ! ( 2 i ξ ) j j ! .
The confluent hypergeometric functions are typically expressed in terms of the rising factorials, instead of the falling factorials. Converting between them gives
( α + 1 ) j + α + + α 1 = ( 1 ) j + α + + α 1 ( 1 α + ) ( j + α + + α 1 ) .
This can then be expressed in terms of Gamma functions as
( α + 1 ) j + α + + α 1 = ( 1 ) j + α + + α 1 Γ ( j + α ) Γ ( 1 α + ) = ( 1 ) j + α + + α 1 ( α ) ( j ) Γ ( α ) Γ ( 1 α + ) .
The term in the denominator can be written as
( j + α + + α 1 ) ! = ( α + + α ) ( j ) Γ ( α + + α ) .
This means we have established that
( α + 1 ) j + α + + α 1 ( j + α + + α 1 ) ! = ( 1 ) j + α + + α 1 ( α ) ( j ) ( α + + α ) ( j ) Γ ( α ) Γ ( 1 α + ) Γ ( α + + α ) .
Using Euler’s reflection formula Γ ( z ) Γ ( 1 z ) = π / sin ( π z ) gives our final result:
Residue = e π 2 δ e 2 i π α + 1 4 π 2 α + + α Γ ( α + ) Γ ( α ) Γ ( α + + α ) e i ξ j = 0 ( α ) ( j ) ( α + + α ) ( j ) ( 2 i ξ ) j j ! .
Multiplying by 2 π i to determine the integral via the calculus of residues then produces a result equal to that in Equation (62), provided the Kummer function satisfies
M ( α , α + + α , 2 i ξ ) = j = 0 ( α ) ( j ) ( α + + α ) ( j ) ( 2 i ξ ) j j ! ,
which is the standard definition of this confluent hypergeometric function, as long as α + + α is not a nonpositive integer. This provides a numerical way to evaluate the wave function via a power series expansion.
We can also evaluate the wave function using a simple numerical integration over the two different integral formulas for Φ . The real-valued integral in Equation (60) is straightforward to evaluate, with the caveat that one evaluates the complex exponentials carefully, noting that x a + i b = x a e i b ln x , for example. Care must also be taken when the real part of the exponents is less than one, because those functions are not so easily integrated using traditional integration rules, without changing variables to remove their non-polynomic behavior first.

5. Contour Integral around Circular Path for the Coulomb Problem

For the contour integral around the circular path (see Equation (54)), we need to carefully determine the phase of the multivalued function in the integrand again. We use our standard approach, relating everything to our reference point. With the angle θ measured relative to the positive imaginary axis, the parametrization of the radius R integral is z = R e i ( θ + π 2 ) = R ( sin θ + i cos θ ) ; for concreteness, we will pick R = 2 in the figures. Note that the conventional angle for describing z in the complex plane is measured from the real axis, hence, the polar angle for z is θ + π 2 . The integral then becomes
Φ ( ξ ) = i 0 2 π d θ R e i ( θ + π 2 ) e R ξ ( sin θ + i cos θ ) R e i ( θ + π 2 ) i α + 1 R e i ( θ + π 2 ) + i α 1 × e i ϕ 2 ( α + 1 ) e i ϕ 1 ( α 1 ) e π 2 δ .
Now we need to determine the phases ϕ 1 and ϕ 2 as a function of θ . The graphics in Figure 6 are helpful for this task, and represents the situation when 0 θ < π .
The phase ϕ 1 corresponds to the angle, from which an arrow is drawn from i to 0 + winds when it moves along the indicated path from the reference point 0 + to a point on the circle (Figure 6a). From the angles and sides of the shaded triangle, we can extract the relation between ϕ 1 and θ . We employ the law of cosines
l 2 = R 2 + 1 2 R cos ( π θ ) = R 2 + 1 + 2 R cos θ ,
and the law of sines
sin ϕ 1 R = sin ( π θ ) l = sin θ l ,
and thus find
sin ϕ 1 = R sin θ R 2 + 1 + 2 R cos θ .
For extracting ϕ 1 , we need to take the arcsin of the right-hand side, and be careful about which quadrant to choose ϕ 1 (will be discussed in detail later). Note that the final result looks like we just drew the angle ignoring the branch cut, but we did carefully follow the procedure of traversing a path that does not cross the branch cut, as required for the determination of the angle.
The phase ϕ 2 we determine in a very similar way. This phase ϕ 2 is the winding angle of an arrow drawn from + i to 0 + , which winds along the indicated path from the reference point 0 + to a point on the circle (see Figure 6b). Note that for θ = 0 , the arrow has already flipped, i.e., rotated by π . The corresponding angle in the shaded triangle is thus 2 π ϕ 2 . Again, we employ the law of cosines
l 2 = R 2 + 1 2 R cos θ = R 2 + 1 2 R cos θ ,
and the law of sines
sin θ l = sin ( 2 π ϕ 2 ) R = sin ϕ 2 R
Thus, we find
sin ϕ 2 = R sin θ R 2 + 1 2 R cos θ .
Now we repeat the same procedure for π θ < 2 π . The corresponding plots for determining ϕ 1 and ϕ 2 are shown in Figure 7. In Figure 7a, the angles in the shaded triangle are θ π and 2 π ϕ 1 . The law of cosines
l 2 = R 2 + 1 2 R cos ( θ π ) = R 2 + 1 + 2 R cos θ ,
and the law of sines
sin ( 2 π ϕ 1 ) R = sin ϕ 1 R = sin ( θ π ) l = sin θ l ,
yield
sin ϕ 1 = R sin θ R 2 + 1 + 2 R cos θ ,
which is the same as Equation (78). For ϕ 2 , the geometric relation are shown in Figure 7b. Again, we use the law of cosines
l 2 = R 2 + 1 2 R cos ( 2 π θ ) = R 2 + 1 2 R cos θ ,
and the law of sines
sin ( 2 π θ ) l = sin θ l = sin ϕ 2 R
and obtain
sin ϕ 2 = R sin θ R 2 + 1 2 R cos θ ,
which is the same as Equation (81). Thus, the same formulae for sin ϕ 1 and sin ϕ 2 are valid for all 0 θ < 2 π . It remains to determine correctly ϕ 1 and ϕ 2 , i.e., choosing them from the right quadrants depending on the integration angle θ . Figure 6 and Figure 7 are helpful in this respect. Let us start with ϕ 1 . The phase ϕ 1 (arrow from i to z) rotates from 0 to π / 2 , when θ rotates from 0 to cos 1 1 R . Then, when θ rotates from cos 1 1 R to π , ϕ 1 rotates from π / 2 to π , etc. Overall, this yields the following relations
0 θ < cos 1 1 R : 0 ϕ 1 < π 2 cos 1 1 R θ < π : π 2 ϕ 1 < π π θ < π + cos 1 1 R : π ϕ 1 < 3 π 2 π + cos 1 1 R θ < 2 π : 3 π 2 ϕ 1 < 2 π .
We can repeat the procedure for ϕ 2 , which regards the arrow drawn from + i to z. We get for ϕ 2
0 θ < cos 1 1 R : π ϕ 2 < 3 π 2 cos 1 1 R θ < π : 3 π 2 ϕ 2 < 2 π π θ < π + cos 1 1 R : 2 π ϕ 2 < 5 π 2 π + cos 1 1 R θ < 2 π : 5 π 2 ϕ 2 < 3 π .
In general, when evaluating the integral, we can simply use a trapezoidal rule, dividing the θ interval evenly. While the result is independent of the radius R, the appearance of z ξ in the exponent of the exponential function produces accuracy issues for large R and ξ —this means, for accurate numerical work, we should use as small an R as possible (we found R = 1.1 to be good with 100 000 steps). It also means at some point, the direct numerical integration will fail when ξ is large enough, due to precision issues similar to using the power series to compute e x for large x.
As an example, we plot three results for the radial Coulomb wave function for small ξ and three different energies in the continuum in Figure 8. These results are normalized according to the result from the one-dimensional integral for comparison. The results agree with the exact results, as expressed in terms of the confluent hypergeometric function, but requires no knowledge of that function. It instead requires just a moderate computing exercise.
All the other continuum solutions except for the Morse potential, proceed in a similar fashion and can be evaluated with the circular contour integral. The Morse potential is different for two reasons: (i) the contour is not given as a rotated dog-bone contour about the imaginary axis and (ii) the function must vanish as x . We treat its continuum solutions next.

6. Continuum Solutions of the Morse Potential

The Morse potential differs from the other continuum cases. We see in Table 5 that the Laplace form of the Schrödinger equation is once again of the form of Equation (11), so ± λ R once again, but now β C . This differs from both the previous problems of the form in Equation (11), where β was real, as well as the previous continuum problems, where β was an integer. As a result, the solution for Φ can still be written in the contour integral form given in Equation (17), but the sum of the exponents α + + α 2 = β 2 is not an integer. This means the integrand is not single-valued as | z | , so we must draw the branch cut as two pieces, each running from infinity to one of the branch points ± λ . Hence, we draw the branch cut along the real axis in two pieces: one piece travels from λ to Re ( z ) and the other piece travels from Re ( z ) = to λ . This branch-cut structure does not allow us to enclose the branch points at ± λ with a contour (as we did for previous continuum cases), nor enclose just one of the two branch points (as in the bound-state problems). However, we can now try another possible contour that connects the two points ± λ . For concreteness, we choose it to run through the origin, and we already see a parallel to the other continuum cases, where we took the limiting behavior of a contour running just next to the imaginary axis. For the Morse potential, we have 0 ξ < . Recall that ξ = 2 2 μ V 0 a e a x , so as x we have ξ and as x we have that ξ 0 . Since the Morse potential becomes large and positive for x , we must have ψ ( x ) 0 . Hence, ϕ ( ξ ) must go to zero as ξ . This is a more stringent condition than just having the wave function be bounded as we used previously. This condition eliminates the new contour from λ + λ , because an analysis similar to the previous cases shows that this solution diverges as ξ .
Now, let us instead consider the contour as shown in Figure 1b: the contour running from λ . The integral solution over that contour is
Φ ( ξ ) = λ d z e ξ z z λ α + 1 z + λ α 1 ,
where λ = 1 2 . Now, we make the substitution z = t λ to obtain
Φ ( ξ ) = e ξ 2 ( 1 ) β 1 0 d t e ξ t t + 1 α + 1 t α 1 .
Immediately, we see that this integral will go to zero as ξ because the exponential term in the integrand guarantees a convergent integral that will vanish in the limit. Moreover, this integral closely resembles the integral representation of the Tricomi confluent hypergeometric function U ( a , b , z ) , given by (Equation 13.4.4 in Ref. [15])
U ( a , b , z ) = 1 Γ ( a ) 0 d t e z t t + 1 b a 1 t a 1 .
This means we can write the integral solution to the Morse differential equation as
Φ ( ξ ) = ( 1 ) β 1 Γ ( α ) e ξ / 2 U ( α , β , ξ ) .
One might be concerned that this function is not bounded for ξ 0 . Indeed, a simple power-counting argument shows that the magnitude of the integrand behaves like 1 z when ξ = 0 . But, because the exponent is complex, it will produce oscillations, which can allow the integral to converge and be bounded. To settle this question, we look at the well-known asymptotics of the Tricomi function U. In the limit where the argument ξ goes to zero, the behavior of U is governed by the value of the real part of the second parameter, in our case β . For the Morse potential, Re ( β ) = 1 , so the asymptotic behavior of U as ξ 0 is (Equation 13.2.18 of Ref. [15])
U ( α , β , ξ ) = Γ ( β 1 ) Γ ( α ) ξ 1 β + Γ ( 1 β ) Γ ( 1 α + ) + O ( ξ 2 Re ( β ) ) .
Once again, Re ( β ) = 1 , so the ξ 1 β term will be ξ i Im ( β ) , which has a modulus of 1 ( ξ i Im ( β ) ξ i Im ( β ) = 1 ) and a phase that varies rapidly as ξ 0 . Thus, the Tricomi function does not diverge, but oscillates, as ξ 0 . We will see below that the behavior generally looks like a cosine of x for x .
This contour yields a solution that will be finite everywhere, as well as going to zero as ξ . Thus, it satisfies all of our requirements, and we can write the solution for the unnormalized Morse potential wave function in one dimension ψ ( ξ ) as
ψ ( ξ ) = ξ β 1 2 e ξ 2 U ( α , β , ξ ) ,
Recall that ξ = 2 2 μ V 0 a e a x in this solution. Finally, we note that we can use the so-called Kummer relation (see Equation (13).2.42 of Ref. [15])
U ( a , b , z ) = Γ ( 1 b ) Γ ( a b + 1 ) M ( a , b , z ) + Γ ( b 1 ) Γ ( a ) z 1 b M ( a b + 1 , 2 b , z ) ,
to relate the above solution in terms of the Tricomi function U to the sum of two complex conjugate Kummer functions M, which is the form of the Morse continuum wave function that appears in the literature [16,17]. The divergence in each M as x is canceled exactly by their sum leading to a finite result. We plot the continuum wave functions for some typical values in Figure 9. You can see the behavior is as anticipated. We have a rapid decay for x < 0 , there is a transition region near x = 1 , and then the form is of a constant amplitude sinusoidal oscillation as x increases in the positive direction.
We summarize all of our continuum solution results in two tables. Table 6 shows the results expressed as functions of ξ , while Table 7 shows the results in terms of the original variables. The normalization of these wave functions is subtle and depends on the scheme that will be used, so we do not discuss the issue of normalization here.

7. Conclusions

We have shown how the almost forgotten Laplace method can be employed as a powerful tool to determine the bound and continuum states of all problems that are solved with confluent hypergeometric equation wave functions. The bound-state solutions provide a way to relate the orthogonal polynomials that arise to their Rodrigues formulas, which naturally emerge via the poles in the contour integrals needed to determine the wave functions. Quantization of the energies in the bound states also occurs naturally. For the continuum solutions, we find a simple contour integral that gives us the wave functions, but it requires a precise determination of the phases of the terms that are raised to complex-valued powers—hence it requires a proper understanding of branch cuts and how to evaluate the polar radius and phase of complex numbers on a cut plane. The main result we developed here is a contour-integral representation of the continuum wave functions, which can be easily programmed, and may be a more useful way to represent these functions than the more common approach via confluent hypergeometric functions.
This method also allows for a solution to the linear potential problem (which we did not discuss here). That problem is treated with methods similar to what we developed here in textbooks. In particular, the textbook by Konishi and Paffuti [7] has a nice treatment of this. Curiously, Landau and Liftshitz [5] and Messiah [6] both use the Laplace method to determine the properties of the confluent hypergeometric functions in their appendix, but they do not use the Laplace method for finding solutions to the Schrödinger equation in the main part of the text.

Author Contributions

Conceptualization, J.K.F.; Methodology, J.K.F., Investigation, J.C., A.G. and J.K.F., Formal analysis, J.C., A.G. and J.K.F., Writing—original draft, J.C., A.G. and J.K.F., Writing—Review and editing, J.C., A.G. and J.K.F., Funding acquisition, J.K.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Science Foundation under Grant No. PHY-1915130 and the McDevitt bequest at Georgetown University.

Data Availability Statement

There is no data used in this study.

Acknowledgments

We acknowledge useful discussions with Wesley Mathews Jr.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Schrödinger, E. Quantisierung als Eigenwertproblem (Erste Mitteilung). Ann. Phys. 1926, 384, 361–376. [Google Scholar] [CrossRef]
  2. Schlesinger, L. Einführung in die Theorie der Differentialgleichungen mit einer unabhängigen Variablen; G. J. Göschensche Verlagshandlung: Leipzig, Germany, 1900. [Google Scholar]
  3. Dirac, P.A.M. Complex Variables in Quantum Mechanics. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 1937, 160, 48–59. [Google Scholar]
  4. Dirac, P.A.M. Principles of Quantum Mechanics, 3rd ed.; Clarendon Press: Oxford, UK, 1947. [Google Scholar]
  5. Landau, L.D.; Lifshitz, E.M. Quantum Mechanics: Non-Relativistic Theory, 3rd ed.; Pergamon Press: Oxford, UK, 1977. [Google Scholar]
  6. Messiah, A. Mécanique Quantique; Dunod: Paris, 1959; translated into English as Quantum Mechanics; John Wiley and Sons: North Holland, The Netherlands, 1966. [Google Scholar]
  7. Konishi, K.; Paffuti, G. Quantum Mechanics: A New Introduction; Oxford University Press: Oxford, UK, 2009. [Google Scholar]
  8. Capri, A.Z. Nonrelativistic Quantum Mechanics; Benjamin Cummings: Menlo Park, CA, USA, 1985. [Google Scholar]
  9. Galler, A.; Canfield, J.; Freericks, J.K. Schrödinger’s original quantum-mechanical solution for hydrogen. Eur. J. Phys. 2021, 42, 035405. [Google Scholar] [CrossRef]
  10. Tsaur, G.-Y.; Wang, J. A universal Laplace-transform approach to solving Schrödinger equations for all known solvable models. Eur. J. Phys. 2014, 35, 015006. [Google Scholar] [CrossRef]
  11. Arda, A.; Sever, R. Exact solutions of the Schrödinger equation via Laplace transform approach: Pseudoharmonic potential and Mie-type potentials. J. Math. Chem. 2012, 50, 971–980. [Google Scholar] [CrossRef]
  12. Das, T. A Laplace transform approach to find the exact solution of the N-dimensional Schrödinger equation with Mie-type potentials and construction of Ladder operators. J. Math. Chem. 2015, 53, 618–628. [Google Scholar] [CrossRef]
  13. Chung, W.; Kim, Y.; Kwon, J. Laplace transform method in one dimensional quantum mechanics on the semi infinite axis. J. Math. Chem. 2022, 60, 1080–1088. [Google Scholar] [CrossRef]
  14. Szegö, G. Orthogonal Polynomials; American Mathematical Society: Providence, RI, USA, 1939. [Google Scholar]
  15. Olver, F.W.; Olde, J.; Daalhuis, A.B.; Lozier, D.W.; Schneider, B.I.; Boisvert, R.F.; Clark, C.W.; Miller, B.R.; Saunders, B.V.; Cohl, H.S.; et al. (Eds.) NIST Digital Library of Mathematical Functions; National Institute of Standards and Technology: Gaithersburg, MD, USA, 2022. Available online: http://dlmf.nist.gov/ (accessed on 1 June 2022).
  16. Nicholls, R.W. Continuum Wavefunctions for Morse Molecules. Chem. Phys. Lett. 1981, 79, 317–320. [Google Scholar] [CrossRef]
  17. Matsumoto, A. Generalized matrix elements in discrete and continuum states for the Morse potential. J. Phys. B At. Mol. Opt. Phys. 1988, 21, 2863–2870. [Google Scholar] [CrossRef]
Figure 1. The contours we consider when α ± are not integers. The branch cuts in each case are shown by the dashed red lines. In panel (a), we show a Hankel-like contour, in panel (b) a contour that runs to infinity, and in panel (c) we show the dog-bone contour that surrounds the two branch points when the branch cut is drawn between them.
Figure 1. The contours we consider when α ± are not integers. The branch cuts in each case are shown by the dashed red lines. In panel (a), we show a Hankel-like contour, in panel (b) a contour that runs to infinity, and in panel (c) we show the dog-bone contour that surrounds the two branch points when the branch cut is drawn between them.
Quantumrep 05 00024 g001
Figure 2. Contours that lead to the correct solution for the wave function and provide the quantization condition for the energy. The contour is always a closed contour in the counter-clockwise direction encircling the point z = λ , as shown in panel (a). In some cases, shown in panel (b), a branch point remains at z = λ and the branch cut runs from there to z = along the real axis.
Figure 2. Contours that lead to the correct solution for the wave function and provide the quantization condition for the energy. The contour is always a closed contour in the counter-clockwise direction encircling the point z = λ , as shown in panel (a). In some cases, shown in panel (b), a branch point remains at z = λ and the branch cut runs from there to z = along the real axis.
Quantumrep 05 00024 g002
Figure 3. Four possible contours for solving the one-dimensional simple harmonic oscillator using the second ansatz. The contours must lie within the white “cone” regions as they go to infinity (there is no other restriction on the contours except not to cross the branch cut for finite values of z). In panel (a), we show the Hankel contour γ 4 , which goes around the branch cut, and γ 5 , a contour that runs parallel to the real axis. In panel (b), we show two contours starting from the branch point and running to infinity either below ( γ 6 ) or above ( γ 7 ) the branch cut.
Figure 3. Four possible contours for solving the one-dimensional simple harmonic oscillator using the second ansatz. The contours must lie within the white “cone” regions as they go to infinity (there is no other restriction on the contours except not to cross the branch cut for finite values of z). In panel (a), we show the Hankel contour γ 4 , which goes around the branch cut, and γ 5 , a contour that runs parallel to the real axis. In panel (b), we show two contours starting from the branch point and running to infinity either below ( γ 6 ) or above ( γ 7 ) the branch cut.
Quantumrep 05 00024 g003
Figure 4. Rotated ’dog-bone’ shaped contour for evaluating the contour integral for the continuum wave functions of the Coulomb problem in two and three dimensions.
Figure 4. Rotated ’dog-bone’ shaped contour for evaluating the contour integral for the continuum wave functions of the Coulomb problem in two and three dimensions.
Quantumrep 05 00024 g004
Figure 5. Procedure for determining the phases of the integrand in Equation (55) along the vertical pieces of the contour (we deformed the contour slightly for clarity in the image). (a) We pick a reference point z = 0 + and show arrows from i and i to the reference point. (b) We draw paths from the reference point 0 + to each point z of this piece of the deformed contour γ c (here running upward vertically to the right of the branch cut). The arrows drawn from the upper and lower branch points to z do not change their net direction as they move along the path from 0 + to z, thus the phase for f ( z ) is the same as the phase for f ( 0 + ) on the right side of the branch cut. (c) For reaching the points of the contour on the left side of the branch cut, the arrow drawn from i to z needs to rotate by 2 π , while the net change in direction of the arrow drawn from the lower branch point is zero; note that the arrow along the path from the reference point to z is allowed to cross the branch cut, even though the path never crosses the branch cut.
Figure 5. Procedure for determining the phases of the integrand in Equation (55) along the vertical pieces of the contour (we deformed the contour slightly for clarity in the image). (a) We pick a reference point z = 0 + and show arrows from i and i to the reference point. (b) We draw paths from the reference point 0 + to each point z of this piece of the deformed contour γ c (here running upward vertically to the right of the branch cut). The arrows drawn from the upper and lower branch points to z do not change their net direction as they move along the path from 0 + to z, thus the phase for f ( z ) is the same as the phase for f ( 0 + ) on the right side of the branch cut. (c) For reaching the points of the contour on the left side of the branch cut, the arrow drawn from i to z needs to rotate by 2 π , while the net change in direction of the arrow drawn from the lower branch point is zero; note that the arrow along the path from the reference point to z is allowed to cross the branch cut, even though the path never crosses the branch cut.
Quantumrep 05 00024 g005
Figure 6. Geometry for determining the relationship between the phases ϕ 1 and ϕ 2 and θ . In the figures, the symbol l is used for the unknown length on each triangle, and R is chosen to equal 2, for concreteness; do not conflate l with the quantum number for total angular momentum. This case corresponds to 0 θ π . Panel (a) shows the geometry for the angles from the lower branch point at i , while panel (b) shows the results for the geometry from the upper branch point i.
Figure 6. Geometry for determining the relationship between the phases ϕ 1 and ϕ 2 and θ . In the figures, the symbol l is used for the unknown length on each triangle, and R is chosen to equal 2, for concreteness; do not conflate l with the quantum number for total angular momentum. This case corresponds to 0 θ π . Panel (a) shows the geometry for the angles from the lower branch point at i , while panel (b) shows the results for the geometry from the upper branch point i.
Quantumrep 05 00024 g006
Figure 7. Similar figure as in Figure 6, except here we have π θ 2 π . Again, panel (a) shows the geometry from the lower branch point, while panel (b) shows the geometry from the upper branch point.
Figure 7. Similar figure as in Figure 6, except here we have π θ 2 π . Again, panel (a) shows the geometry from the lower branch point, while panel (b) shows the geometry from the upper branch point.
Quantumrep 05 00024 g007
Figure 8. Plot of the continuum Coulomb wave functions for l = 0 and three different values of E: (a) E = 0.1 ; (b) E = 1 , and (c) E = 10 (all energies are in Hartrees). The power series is shown in red, the contour integral with R = 1.1 in green, and the one-dimensional integral in black. The three approximations lie on top of each other until they start to fail—the power series fails around ξ 20 , while the contour integral fails around ξ 30 . The errors typically occur due to the loss of digits of precision in the expressions being evaluated.
Figure 8. Plot of the continuum Coulomb wave functions for l = 0 and three different values of E: (a) E = 0.1 ; (b) E = 1 , and (c) E = 10 (all energies are in Hartrees). The power series is shown in red, the contour integral with R = 1.1 in green, and the one-dimensional integral in black. The three approximations lie on top of each other until they start to fail—the power series fails around ξ 20 , while the contour integral fails around ξ 30 . The errors typically occur due to the loss of digits of precision in the expressions being evaluated.
Quantumrep 05 00024 g008
Figure 9. Plot of the continuum Morse wave function for three different values of E: (a) E = 0.1 ; (b) E = 1 , and (c) E = 10 . We use 2 a 2 2 μ as the energy unit and 1 a as the length unit. The case we consider is for V 0 = 2 a 2 2 μ . Note how the continuum wave rapidly decays for x < 0 , where the Morse potential becomes large and positive. Because the Morse potential decays to zero exponentially fast, the continuum solution rapidly looks like a simple cosine wave for large positive x with a fixed amplitude. In the region around x = 0 we see a transition between the two behaviors.
Figure 9. Plot of the continuum Morse wave function for three different values of E: (a) E = 0.1 ; (b) E = 1 , and (c) E = 10 . We use 2 a 2 2 μ as the energy unit and 1 a as the length unit. The case we consider is for V 0 = 2 a 2 2 μ . Note how the continuum wave rapidly decays for x < 0 , where the Morse potential becomes large and positive. Because the Morse potential decays to zero exponentially fast, the continuum solution rapidly looks like a simple cosine wave for large positive x with a fixed amplitude. In the region around x = 0 we see a transition between the two behaviors.
Quantumrep 05 00024 g009
Table 1. Quantum-mechanical potentials with bound states that are analyzed in this work. For each, we give the form of the potential, the general form of the wave function where Φ is the unknown part determined by using the Laplace method, and the form for the independent variable ξ used for each problem; in one-dimension, the independent variable is x, in two dimensions it is ρ = x 2 + y 2 , and in three dimensions it is r = x 2 + y 2 + z 2 . Notably, 0 ξ < in all cases but the last (where ξ ). In the table, m is the z-component angular momentum quantum number with eigenvalues m and l is the total angular momentum quantum number with eigenvalue 2 l ( l + 1 ) . Moreover, μ is the mass of the (effective) particle, ω is the angular frequency of the oscillator, e is the magnitude of the charge of an electron, is Planck’s constant, E is the energy of the corresponding energy eigenstate, a is a real constant with units of inverse length, and V 0 has units of energy.
Table 1. Quantum-mechanical potentials with bound states that are analyzed in this work. For each, we give the form of the potential, the general form of the wave function where Φ is the unknown part determined by using the Laplace method, and the form for the independent variable ξ used for each problem; in one-dimension, the independent variable is x, in two dimensions it is ρ = x 2 + y 2 , and in three dimensions it is r = x 2 + y 2 + z 2 . Notably, 0 ξ < in all cases but the last (where ξ ). In the table, m is the z-component angular momentum quantum number with eigenvalues m and l is the total angular momentum quantum number with eigenvalue 2 l ( l + 1 ) . Moreover, μ is the mass of the (effective) particle, ω is the angular frequency of the oscillator, e is the magnitude of the charge of an electron, is Planck’s constant, E is the energy of the corresponding energy eigenstate, a is a real constant with units of inverse length, and V 0 has units of energy.
ProblemPotentialIndependent
Variable
Wavefunc.
Form
1D SHO,
Even
V = 1 2 μ ω 2 x 2 ξ = μ ω x 2 Φ ( ξ )
1D SHO,
Odd
V = 1 2 μ ω 2 x 2 ξ = μ ω x 2 x Φ ( ξ )
2D
SHO
V = 1 2 μ ω 2 ρ 2 ξ = μ ω ρ 2 ρ | m | Φ ( ξ ) e i m ϕ
3D
SHO
V = 1 2 μ ω 2 r 2 ξ = μ ω r 2 r l Φ ( ξ ) Y l m ( θ , ϕ )
2D
Coulomb
V = e 2 ρ ξ = 2 μ E 2 ρ ρ | m | Φ ( ξ ) e i m ϕ
3D
Coulomb
V = e 2 r ξ = 2 μ E 2 r r l Φ ( ξ ) Y l m ( θ , ϕ )
Morse
Potential
V = V 0 e 2 a x 2 e a x ξ = 2 2 μ V 0 a e a x ξ 2 μ E a Φ ( ξ )
1D SHO
Method 2
V = 1 2 μ ω 2 x 2 ξ = μ ω x e μ ω 2 x 2 Φ ( ξ )
Table 2. More details on the solutions for the wave functions using the Laplace method. The 1D SHO can be solved in two ways with the Laplace method, and we treat them both. In the second column, we provide the form of the Schrödinger equation obtained by the substitutions detailed in Table 1. Notably, for all but the last case, the coefficient of the first derivative term is either an integer or a half-odd integer, except in the Morse potential. It is always larger than 1 for all of these cases except the first row. Moreover, a 0 is the reduced Bohr radius ( a 0 = 2 / μ e 2 ).
Table 2. More details on the solutions for the wave functions using the Laplace method. The 1D SHO can be solved in two ways with the Laplace method, and we treat them both. In the second column, we provide the form of the Schrödinger equation obtained by the substitutions detailed in Table 1. Notably, for all but the last case, the coefficient of the first derivative term is either an integer or a half-odd integer, except in the Morse potential. It is always larger than 1 for all of these cases except the first row. Moreover, a 0 is the reduced Bohr radius ( a 0 = 2 / μ e 2 ).
ProblemLaplace Form of the Schrödinger Equation α ± = β λ ± δ 2 λ
1D SHO,
Even
ξ Φ + 1 2 Φ + E 2 ω 1 4 ξ Φ = 0 1 4 ± E 2 ω
1D SHO,
Odd
ξ Φ + 3 2 Φ + E 2 ω 1 4 ξ Φ = 0 3 4 ± E 2 ω
2D
SHO
ξ Φ + ( | m | + 1 ) Φ + E 2 ω 1 4 ξ Φ = 0 2 | m | + 1 2 ± E 2 ω
3D
SHO
ξ Φ + ( l + 3 2 ) Φ + E 2 ω 1 4 ξ Φ = 0 2 l + 3 4 ± E 2 ω
2D
Coulomb
ξ Φ + ( 2 | m | + 1 ) Φ + 2 a 0 2 μ E ξ Φ = 0 2 | m | + 1 2 ± a 0 2 μ E
3D
Coulomb
ξ Φ + 2 ( l + 1 ) Φ + 2 a 0 2 μ E ξ Φ = 0 l + 1 ± a 0 2 μ E
Morse
Potential
ξ Φ + 2 2 μ E a + 1 Φ + 2 μ V 0 a 1 4 ξ Φ = 0 2 μ E ± 2 μ V 0 a + 1 2
1D SHO
Method 2
Φ 2 ξ Φ + 2 E ω 1 Φ = 0 N/A
Table 3. Quantization condition for the Laplace method (by tradition, the principal quantum number n starts from 0 for all cases, except the Coulomb cases, where it starts from | m | + 1 in two dimensions and from l + 1 in three dimensions); N is required to be a nonnegative integer from the quantization condition arising from Laplace’s method. By tracing back through the definitions of Φ and ξ in each case, one obtains the standard wave functions for each problem, up to a normalization constant, that still needs to be determined. The δ in the index of the associated Laguerre polynomial in the last column of the Morse potential satisfies δ = 2 μ V 0 a . All models, except for the Morse potential, have an infinite number of bound states. The Morse potential has a finite number, where we are required to have n < 2 μ V 0 / a . H n ( ξ ) is the nth Hermite polynomial.
Table 3. Quantization condition for the Laplace method (by tradition, the principal quantum number n starts from 0 for all cases, except the Coulomb cases, where it starts from | m | + 1 in two dimensions and from l + 1 in three dimensions); N is required to be a nonnegative integer from the quantization condition arising from Laplace’s method. By tracing back through the definitions of Φ and ξ in each case, one obtains the standard wave functions for each problem, up to a normalization constant, that still needs to be determined. The δ in the index of the associated Laguerre polynomial in the last column of the Morse potential satisfies δ = 2 μ V 0 a . All models, except for the Morse potential, have an infinite number of bound states. The Morse potential has a finite number, where we are required to have n < 2 μ V 0 / a . H n ( ξ ) is the nth Hermite polynomial.
ProblemQuantization
Condition
Energy
Quantization, E n
Form of Φ ( ξ )
1D SHO,
Even
N = n = E 2 ω 1 4 ω 2 n + 1 2 e ξ / 2 L n 1 2 ( ξ )
1D SHO,
Odd
N = n = E 2 ω 3 4 ω 2 n + 1 + 1 2 ξ 1 / 2 e ξ / 2 L n 1 2 ( ξ )
2D
SHO
N = n = E 2 ω 2 | m | + 1 2 ω 2 n + | m | + 1 e ξ / 2 L n ( | m | ) ( ξ )
3D
SHO
N = n = E 2 ω 1 2 2 l + 3 2 ω 2 n + l + 3 2 e ξ / 2 L n l + 1 2 ( ξ )
2D
Coulomb
N = n | m | 1 = a 0 2 μ E | m | + 1 2 2 2 μ a 0 2 n 1 2 2 e ξ L n | m | 1 ( 2 | m | ) ( ξ )
3D
Coulomb
N = n l 1 = a 0 2 μ E l + 1 2 2 μ a 0 2 n 2 e ξ L n l 1 ( 2 l + 1 ) ( ξ )
Morse
Potential
N = n = 2 μ V 0 2 μ E a a 2 2 2 μ n 2 μ V 0 a 2 e ξ / 2 L n ( 2 n 2 δ 1 ) ( ξ )
1D SHO
Method 2
N = n = E ω 1 2 ω n + 1 2 H n ξ
Table 4. Summary for how to convert the Schrödinger equation into the Laplace equation for the five problems with continuum solutions. For each, we give the form of the potential, the general form of the wave function where Φ is the part of the solution that is found by using the Laplace method, and the independent variable ξ used for each problem. Note that E > 0 for these continuum problems. We have m denoting the quantum number for the z-component of angular momentum and l denoting the quantum number for the total angular momentum and Y l m ( θ , ϕ ) the spherical harmonic.
Table 4. Summary for how to convert the Schrödinger equation into the Laplace equation for the five problems with continuum solutions. For each, we give the form of the potential, the general form of the wave function where Φ is the part of the solution that is found by using the Laplace method, and the independent variable ξ used for each problem. Note that E > 0 for these continuum problems. We have m denoting the quantum number for the z-component of angular momentum and l denoting the quantum number for the total angular momentum and Y l m ( θ , ϕ ) the spherical harmonic.
ProblemPotentialIndependent
Variable
Wavefunc.
Form
2D Free
Particle
V = 0 ξ = 2 μ E 2 ρ ρ | m | Φ ( ξ ) e i m ϕ
3D Free
Particle
V = 0 ξ = 2 μ E 2 r r l Φ ( ξ ) Y l m ( θ , ϕ )
2D
Coulomb
V = e 2 ρ ξ = 2 μ E 2 ρ ρ | m | Φ ( ξ ) e i m ϕ
3D
Coulomb
V = e 2 r ξ = 2 μ E 2 r r l Φ ( ξ ) Y l m ( θ , ϕ )
Morse
Potential
V = V 0 e 2 a x 2 e a x ξ = 2 2 μ V 0 a e a x ξ i 2 μ E a Φ ( ξ )
Table 5. Final differential equation and exponents α ± for continuum cases to be solved by the Laplace method. The second column summarizes the final form of the Schrödinger equation obtained by the substitutions detailed in Table 4. Note that because δ = 0 for the free-particle cases, there is only one exponent for those problems. In all problems except for the Morse potential, we have λ ¯ = 1 .
Table 5. Final differential equation and exponents α ± for continuum cases to be solved by the Laplace method. The second column summarizes the final form of the Schrödinger equation obtained by the substitutions detailed in Table 4. Note that because δ = 0 for the free-particle cases, there is only one exponent for those problems. In all problems except for the Morse potential, we have λ ¯ = 1 .
ProblemLaplace Form of the
Schrödinger Equation
α ±
2D Free
Particle
ξ Φ + ( 2 | m | + 1 ) Φ + ξ Φ = 0 | m | + 1 2
3D Free
Particle
ξ Φ + 2 ( l + 1 ) Φ + ξ Φ = 0 l + 1
2D
Coulomb
ξ Φ + ( 2 | m | + 1 ) Φ + 2 a 0 2 μ E + ξ Φ = 0 | m | + 1 2 i a 0 2 μ E
3D
Coulomb
ξ Φ + 2 ( l + 1 ) Φ + 2 a 0 2 μ E + ξ Φ = 0 l + 1 i a 0 2 μ E
Morse
Potential
ξ Φ + 2 i 2 μ E a + 1 Φ + 2 μ V 0 a 1 4 ξ Φ = 0 i 2 μ E ± 2 μ V 0 a + 1 2
Table 6. Summary of the results of the Laplace method for continuum cases in terms of the variable ξ . The solution for Φ ( ξ ) , as defined in Table 5, is expressed in terms of the confluent hypergeometric functions M ( a , b , z ) and U ( a , b , z ) .
Table 6. Summary of the results of the Laplace method for continuum cases in terms of the variable ξ . The solution for Φ ( ξ ) , as defined in Table 5, is expressed in terms of the confluent hypergeometric functions M ( a , b , z ) and U ( a , b , z ) .
ProblemConfluent Hypergeometric Form of Φ ( ξ )
2D Free Particle e i ξ M | m | + 1 2 , 2 | m | + 1 , 2 i ξ
3D Free Particle e i ξ M l + 1 , 2 l + 2 , 2 i ξ
2D Coulomb e i ξ M | m | + 1 2 + i a 0 2 μ E , 2 | m | + 1 , 2 i ξ
3D Coulomb e i ξ M l + 1 + i a 0 2 μ E , 2 l + 2 , 2 i ξ
Morse Potential ( 1 ) β 1 Γ ( α ) e ξ / 2 U ( α , β , ξ )
Table 7. Summary of the results of the continuum cases we solved with the Laplace method in terms of the original independent variable. For the free particle cases, we express the more common form of the confluent hypergeometric function. Here, J n ( x ) is the Bessel function of the first kind, and j n ( x ) is the spherical Bessel function of the first kind.
Table 7. Summary of the results of the continuum cases we solved with the Laplace method in terms of the original independent variable. For the free particle cases, we express the more common form of the confluent hypergeometric function. Here, J n ( x ) is the Bessel function of the first kind, and j n ( x ) is the spherical Bessel function of the first kind.
ProblemUnnormalized Wave function
2D Free Particle J | m | 2 μ E 2 ρ e i | m | ϕ
3D Free Particle j l 2 μ E 2 r Y l m θ , ϕ
2D
Coulomb
ρ | m | e i 2 μ E 2 ρ e i | m | ϕ × M | m | + 1 2 + i a 0 2 μ E , 2 | m | + 1 , 2 i 2 μ E 2 ρ
3D
Coulomb
r l e i 2 μ E 2 r Y l m ( θ , ϕ ) × M l + 1 + i a 0 2 μ E , 2 l + 2 , 2 i 2 μ E 2 r
Morse
Potential
2 2 μ V 0 a e a x i 2 μ E a exp 2 μ V 0 a e a x U i 2 μ E 2 μ V 0 a + 1 2 , 2 i 2 μ E a + 1 , 2 2 μ V 0 a e a x
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Canfield, J.; Galler, A.; Freericks, J.K. The Laplace Method for Energy Eigenvalue Problems in Quantum Mechanics. Quantum Rep. 2023, 5, 370-397. https://doi.org/10.3390/quantum5020024

AMA Style

Canfield J, Galler A, Freericks JK. The Laplace Method for Energy Eigenvalue Problems in Quantum Mechanics. Quantum Reports. 2023; 5(2):370-397. https://doi.org/10.3390/quantum5020024

Chicago/Turabian Style

Canfield, Jeremy, Anna Galler, and James K. Freericks. 2023. "The Laplace Method for Energy Eigenvalue Problems in Quantum Mechanics" Quantum Reports 5, no. 2: 370-397. https://doi.org/10.3390/quantum5020024

Article Metrics

Back to TopTop