On the Method of Stationary Phase in Calculating the Propagator of the Coulomb–Kepler Problem

A theorem is proven on the method of stationary phase to be applied in the calculation of the propagator from the Green’s function of the Coulomb - Kepler problem, by means of Fourier transformation. The Green’s function is available in position representation in an exact compact form as ingeniously derived by Hostler (1963). The Fourier transform, at first, exists in the distributional sense. In order to avoid series representation in Schwartz space, the integrals are defined away from the real axis of the integration variable and by taking principal values. The theorem deals with integrands with a highly oscillating exponential term. It is shown that the method of stationary phase renders exact results in cases relevant to the problem, whereby the integrand can contain an analytic amplitude factor in addition to the exponential term.


Introduction
The propagator K considered is defined by the unitary operator Θ(t) exp[−iH t∕ℏ] where the Heaviside function Θ(t) = 1 if t > 0 and zero if t < 0 . In position representation an initial wave function 0 (r � ) is propagated to (r, t) by K as In terms of a complete set of eigenfunctions with H = E , one writes By differentiation with respect to t, one immediately finds that K(r, r � , t) obeys the following partial differential equation (PDE) By time Fourier transformation, one defines the energy amplitude K E and subsequently makes use of (2) as follows: For the non-relativistic Hamiltonian H of the Coulomb-Kepler problem, in a milestone achievement, Hostler [1,2] succeeded to exactly determine the Green function G E in a compact form. The Hamiltonian reads In (5), the coupling constant is defined for the electric and gravitational case, respectively; it is positive, except for repellent charges; m denotes the reduced mass of the twobody problem. The Hamiltonian H encompasses both the hydrogen atom and the two-body Kepler problem of celestial mechanics which differ in scale up to 80 decimals. Hostler's solution will be noted further below. We try to derive the (1) (r, t) = ∫ dr � K(r, r � , t) 0 (r � ), t > 0.
(2) K(r, r � , t) = Θ(t) ∑ exp[−iE t∕ℏ] (r) * (r � ).  propagator K from Hostler's Green's function by means of the inverse Fourier transformation as follows: In Sect. 3, the potential-free case as noted in [3] will be recovered. By taking the limiting procedure above, we avoid using series representations in Schwartz space [4], and simultaneously define how to approach or sidestep the essential singularity at E = ∞ . We remark on different notations: In [5], the propagator is termed kernel, also the notation resolvent is sometimes used instead of Green's function in operator form.
By its explicit time dependence, the propagator delivers more immediate information than the Green's function. By knowing K(r, � , t) , one can follow, for instance, the scattering process of a wave packet through the whole time interval, finite or infinite, including the behavior near the scattering center, at variance with the asymptotic scattering description.
Stimulated by Richard Feynman, there have been attempts to calculate the propagator by the method of path integrals; see [6,7] and for a critical study [8]. Also encouraged by Feynman, Hostler's Green's function was used to calculate the propagator in [9]. Recently, a numerical investigation was published [10], where the propagator is represented in a basis of square-integrable Sturmian functions. All these attempts that can be said were only partially successful and so far did not lead to important applications.
The present study is part of a further attempt to calculate the propagator from Hostler's Green' s function. It was motivated by our previous study of the quantum mechanics of rectilinear orbits [11,12]. There, the Kustaanheimo-Stiefel transformation [13] was applied which leads to a model Hamiltonian of four decoupled harmonic oscillators. One could then apply coherent states of the harmonic [14] and of the reversed harmonic oscillator [15,16]. As a shortcoming of the method, the evolution of the initial state is controlled by the eccentric anomaly and not by the Schrödinger time t. The necessary condition that eccentric anomaly and time are in one-one correspondence failed to be true before the simulated rectilinear orbit entered the physically interesting region.
In the next Section, we deform the integration path of the E integral in (6) to a quarter arc circle close to the singularity at E = ∞ . One observes a highly oscillating exponential phase in the integrand. In Sect. 3, as a kind of litmus test, the propagator of the potential-free case is determined from the Green's function: (i) by exact integration of the modified Fourier transformation in terms of the complex error function and its asymptotic form, (ii) by the method of stationary phase. In Sect. 4, the method of stationary phase is generalized to include integrands which will be relevant to the problem of non-zero potential, work which is still in progress. Appendices 1 and 2 contain proofs on the equivalence of the method of stationary phase and the exact calculation of certain improper integrals.

Transformation of the Fourier Integral in Complex Energy Plane
Hostler's Green's function G Ho differs from G E , defined in (4), by a constant factor as To write down the exact solution of the Fourier transformed PDE (3), given by Eq. (8) of [2], we state the symbols there used where x, y are the Lambert coordinates, k is a complex wave number, and a complex number which is proportional to the coupling constant and, via k, depends on the energy E (by the notation of [2], = Ze 2 ∕(4 )) . The functions M and W denote the Whittaker functions as defined in [17], Γ is the Gamma function. The solution, which factorizes in the Lambert coordinates, reads In the potential-free case, the coupling constant = 0 and thus = 0 . The Whittaker functions are then elementary functions as which give rise to the well known free-particle Green's function; see e.g. [18,19], The Green function depends on k = k(E) . We set In (12), the signs of the square roots are chosen such that in the case E > 0 one has outgoing radial waves, whereas, x = r + r � + r 12 , y = r + r − r 12 , if E < 0 , the free-particle Green's function (11) describes evanescent radial waves. In (6), we apply the transformation E → = i∕k(E) , which is well defined, since differs from zero for | E | ≤ R , and > 0.
We assume that the integrand in (6) is analytic within the closed contour of Fig. 1. Then, by Cauchy's theorem, the original integral is equivalent to the integral along the quarter arc of the circle: We have to comment on the analyticity assumption. The closed path lies in the first quadrant with Re( ) ≡ R ≥ 0 and Im( ) ≡ I ≥ 0 , whereby I can get zero only close to the origin where also R gets small. The Γ factor which appears in the integrand (6) through G Ho , depends on through Γ(1 − i ) = Γ(1 − 2 ) ; since I > 0 , poles cannot appear; the poles of Γ(z) lie on the negative real axis at integer values z = −n , n = 0, 1, ... , and are connected with the discrete part of the energy spectrum. As to the Whittaker functions which constitute G Ho , we use their connections with the confluent hypergeometric functions Φ (Kummer function) and U (Tricomi function) as defined in [20]: where z = −ik x or z = −ik y with the Lambert coordinates x and y larger zero. The function Φ(a, b, z) is analytic in all three complex entries, whereas U(a, b, z) is multivalued in z with the branch cut defined along the negative real axis of the complex variable z. But on the contour, Im(k) ≡ k I > 0 which implies that Real(z) > 0 ; k I > 0 follows from the relations k I = I ∕| | 2 > 0 . By (13), dE∕d is analytic within the contour. In concluding, the integrand in (6) is analytic within the closed path and, thus, justifies the application of Cauchy's integral theorem.
In view of (12), we have the following relations on the circle:

Fig. 1 The integration path (E) is exemplified by a parametric plot in the interval
The magnitudes are nondimensionlized. The singularity at (E = ±∞) = 0 is circumvented by the quarter arc of a circle of radius which eventually shrinks to zero. For demonstration, the imaginary part > 0 of the energy E has some finite value The endpoints E 1,2 = ∓R are real numbers. For small , the corresponding angles 1,2 are close to 0 and ∕2 , respectively. The conditions Im[E( 1,2 )] = 0 imply and We remark that 2 is a dimensionless number with both > 0 and > 0 arbitrary small. Now, we blow up the circle arc by the transformation → A where which implies The end points E 1,2 on the -circle correspond to the values In conclusion, the Fourier integral (6) is transformed as follows; with E = −A 2 , k = i A and t > 0 . It is observed that, through the factor exp(iA 2 t∕ℏ) , there is an unbounded variation of the integrand in the limit → 0.

The Free-Particle Propagator
From Eq. (3.113) of [3], we quote the propagator as As compared to the method of Fourier transformation applied below, the above result is quite easily derived from the spectral representation (2).

Applying the Method of Fourier Transformation
In the following, we reproduce the standard result (23) by the method of Fourier transformation (FT), using the free-particle Green's function G (0) E given in (11). The FT method cannot be avoided in the present context of the given Hostler's Green's function. The application of (22) (17) to the potential-free case constitutes a kind of litmus test. In view of (22) and (11), we have to consider with c = t∕ℏ and s = r 12 . The integral F(c, s) can be expressed in terms of the complex error function erf. Subsequently, in the limit → 0 , asymptotic forms of erf are applied with due attention paid to the occurrence of a sign function which is decisive and depends on the phases of the complex boundary values A 1 and A 2 , respectively. Actually, the sign function is mentioned in [21], but we did not find it in standard handbooks of mathematical physics like [20,22,23] nor in the monograph [24]. Furthermore, for the existence of the integral F in the limit → 0 , we introduce an infinitesimal dissipation as We give details in Appendix 1 including a numerical confirmation of the sign function in the asymptotic forms of the error function. As is shown in Appendix 1, the standard result (23) for the free-particle propagator is recovered by the Fourier integral (24).

Applying the Approximation of Stationary Phase
In the integrand of F in (24), the exponent Φ = iA 2 c − sA is highly oscillating near the integration boundaries A 1,2 when comes close to zero. The phase Φ is stationary when d Φ∕dA=0, which is the case at A = A s with We assume, at first, that both c and s are larger than zero, so A s lies on the negative imaginary axis. The integration path is deformed such that it passes through A s along the imaginary axes. With A = A s − i , we have the identity It is elementary to determine the following stationary phase (sp) integral F sp 0 for c > 0 where the last expression coincides with F 0 obtained in (39) by means of the complex error function. The procedure above was inspired by Sec.3 of the monograph [24].
The method of stationary phase is exact in the given case. Apparently, the method picks out here the relevant interval of the integrand around the stationary point whereby the contribution of the remaining interval oscillates to zero.

Generalization of the Method of Stationary Phase
A theorem will be proved for the following integrand F 1 which is generalized by an analytic amplitude function f(A) as

Theorem
With the definitions (30) and (31), the application of the method of stationary phase is exact as follows: The proof is given in Appendix 2.

Summary
The route was outlined to determine the propagator in position representation by means of energy Fourier transformation (FT) from the exact Green's function of the Coulomb-Kepler problem as derived by Hostler [1,2]. The FT was defined as a principal value and in addition deformed in the complex energy plane in order to deal with the singularity at the energy E = ∞ . The method was tested for the case of the potential-free case whereby the standard textbook result was recovered (i) by means of the analytic performance of the FT, and (ii) by means of the method of stationary phase. In the latter case, a theorem was proved which gives exact results for a certain class of integrands which are relevant to the problem with nonzero potential. The latter problem is still in progress.

Apprendix 1: Propagator of Free-Particle by Fourier Transformation
By (24), we have to determine the following integral: In terms of the error function erf, one obtains For small , i.e., for large | A | , the absolute value | Z | is of order 1∕ → ∞.
We invoke the asymptotic approximation of the error function paying attention to the offset (Z) = ±1 which depends on the phase of the complex number Z; see [21], with where, by definition, − < arg(Z) < . With the aid of (32) to (35), one gets where by c → c − i 3∕2 we introduced the infinitesimal dissipation according to (25). One has to keep in mind that both c ≡ t∕� > 0 and s > 0 . In the limit → +0 , the modulus of the first term within the curly bracket of T 1 vanishes proportional to exp(−s∕ ) , and the second term of T 1 proportional to exp(−1∕ √ ) . As a consequence, We are left with the sign function 1,2 = (Z 1,2 ) , which depends on the phases 1,2 of Z 1,2 = Abs[Z 1,2 ] exp(i 1,2 ) . One easily verifies that in the limit → 0 , the phases 1,2 do not depend on the parameters c > 0 and s > 0 . One finds (32) which implies that 1 = −1 and 2 = 1 with the result which leads to the following textbook result [3], already stated in (23), We conclude this appendix by determining the sign functions (Z 1,2 ) numerically. To this end, we calculate erf(Z 1,2 ) for decreasing → n = exp(−n) , n = 1, 2, 3, 4 . Results, by means of the built-in Mathematica routine for the error function [25], are given in Table 1.
From the numbers in Table 1, it is inferred that 1 = −1 and 2 = +1 which is consistent with the asymptotic formula for the error function in [21].

Appendix 2: Proof of Theorem
We define the integrals J n as and take advantage of the related recurrence relation

Lemma:
If one defines then (38) lim where Later on, the range of the variable s can be extended by analytical continuation.

Proof of Lemma
The proof results from induction. By means of exact integration, one finds for first low-order cases that We assume that (45) is true up to some integer N > 1 and show that To this end, we replace in the expression (43) of J sp N the parameter A s by −is∕(2c) and differentiate with respect to s to obtain where the latter identity holds by induction assumption. The square bracket can be related to the Hermite polynomials H n (x) by taking into account the following generating function: We substitute x for s as follows: to obtain by the induction assumption the representation (51) below which implies the following relation:  With the aid of (51) and the recurrence relation H N+1 = 2xH N (x) − 2NH N−1 (x) , Eq. (47) can be written as follows: We make use of the definition of in (49), of 2 = 1 , of 4c 2 = −i , and of (51) to arrive at the desired result which proves the claim (46) of the Lemma. As a consequence, we proved the following:

Theorem
For analytic amplitudes f(A), which can be represented by a power series in A, and for real parameters c, the integral from A 1 to A 2 can be exactly calculated by the method os stationary phase as follows: