Graph Expansions and Graphical Enumeration Applied to Semiclassical Propagator Expansions

S. A. Fulling

Mathematics Department, Texas A&M University, College Station, Texas, 77843–3368, U.S.A.

Talk presented at the International Conference on Symmetry in Nonlinear Mathematical Physics, Kiev, July, 1995

In recent years T. A. Osborn and his coworkers at the University of Manitoba have extensively developed the well known connected graph expansion and applied it to a wide variety of problems in semiclassical approximation to quantum dynamics [2, 5, 7, 19, 21, 22, 26, 27]. The work I am reporting on attempts to make their work more concrete by applying it to particular quantum systems. The mathematical tools used are some graph theory and some computation, numerical and symbolic.

Today I shall consider only the simplest case, the time-dependent Schrödinger equation with a time-independent scalar potential:

ih̄ \frac{∂ψ}{∂t} = Hψ, H = − \frac{h̄^{2}}{2m}∇^{2} + λV (x).

The propagator K solves the initial-value problem:

ψ(t, x) = ∫_{R^{d}} dy K(t, x, y)ψ(0, y).

When V ≡ 0 we have the free propagator,

K_{0}(t, x, y) = (\frac{m}{2πih̄t})^{d/2} e^{im|x−y|^{2}/2h̄t}.

For any time-independent Hamiltonian, K(t, x, y) is the integral kernel of the operator e^{−itH/h̄}. Therefore, it depends on only two mathematically independent parameters:

\begin{aligned}\frac{tH}{h̄} & = − \frac{h̄t}{2m}∇^{2} + \frac{λt}{h̄} V (x) \\ & ≡ −\frac{1}{2}A∇^{2} + BV. \end{aligned}

This observation is essential in understanding what is meant by “semiclassical expansion”, a term otherwise beset by ambiguity and obscurity.

A digression on Planck’s constant

From a mathematical point of view, h̄ is the constant that appears in the definition of the (inverse) Fourier transform,

ψ(x) = (2πh̄)^{−d/2} ∫_{R^{d}} dp e^{ip·x/h̄} ϕ(p).

• The experimentalist’s h̄ is a constant of nature, equal to 1.054 × 10^{−27} erg-sec.

• The pure theorist is accustomed to choosing units in which h̄ = 1, so that the transform

is

ψ(x) = (2π)^{−d/2} ∫_{R^{d}} dp e^{ip·x} ϕ(p).

• The numerical analyst is more likely to choose h̄ = \frac{1}{2π}, so that the algebra of the discrete Fourier transform comes out simpler;

ψ(x) = ∫_{R^{d}} dp e^{2πip·x} ϕ(p).

• But the asymptotic analyst’s h̄ is none of these; it is an arbitrary small positive parameter: h̄ ↓ 0.

Often it is said that a semiclassical approximation is accurate when h̄ is small. But whatever can that mean, if one can choose units to give h̄ whatever value one likes? Clearly, what must be small is one or more ratios of h̄ to other physical quantities with the same dimensions, and it is those other quantities that determine whether a situation is semiclassical or not. And here is where one needs to make some careful distinctions.

We are studying

\frac{tH}{h̄} = −\frac{1}{2}A∇^{2} + BV.

It involves

an adiabatic parameter, A = \frac{h̄t}{m} , and a coupling constant, B = \frac{λt}{h̄} .

Note that

AB = \frac{λt^{2}}{m} , \frac{A}{B} = \frac{h̄^{2}}{mλ} .

One can investigate K perturbatively, treating any one of these four parameters as small. Physically —

1. An expansion in B is effectively an expansion in λ (traditional perturbation theory in the coupling constant). Each factor B or λ in such a series is attached to a factor V (or one of its derivatives). The magnitude of B answers the question, “How weak is the applied field?”

2. An expansion in A is effectively one in \frac{1}{m} (i.e., m → ∞), which is a kind of semiclassical expansion (the Wigner–Kirkwood series), but is more properly called an adiabatic expansion. Each factor A or \frac{1}{m} in such a series is attached to two differentiation operations upon V. The magnitude of A answers the question, “How slowly varying is the field?” This interpretation is clarified by the variable change x = \sqrt{A} z, which makes

\frac{Ht}{h̄} = −\frac{1}{2}∇_{z}^{2} + BV (\sqrt{A}z).

3. an expansion in AB is effectively one in t, the small-time expansion, in which both powers and derivatives of the potential are treated as small. This expansion (for more complicated operators H) is the subject of intense study in quantum gravity and index theory, but from the point of view of quantum mechanics it seems overly drastic — insufficiently ambitious.

4. A true semiclassical expansion (h̄ → 0) is a series in A/B. Note that in a sense it assumes strong coupling. This is the higher-dimensional, time-dependent version of the famous WKB expansion.

The key distinction between cases 2 and 4 is that in the former the underlying “classical paths” are simply straight lines, whereas in the latter they are the exact classical trajectories of the mechanical system with Hamiltonian H.

Today I concentrate on the adiabatic, or large-mass, expansion. However, I must emphasize that the connected-graph method is applicable to the WKB expansion as well, and also to Hamiltonians involving electromagnetic and gravitational potentials, possibly time-dependent [2, 19, 22, 27]. For technical reasons I assume that V is a C^{∞} function.

The striking fact is that the Wigner–Kirkwood series can be found to all orders in closed form, and its coefficients have a simple combinatorial structure. The derivation [7, 21, 26] starts with a conventional Feynman–Dyson expansion in B; this yields Feynmantype graphs representing convolutions of K_{0} factors and potentials. Then a further expansion in A expands these integrals as more local functionals of the derivatives of the potential to all orders. Under appropriate technical conditions, the resulting series is rigorously asymptotic [33].

The series can be recast into the form K = K_{0}e^{L}, where

L(x, y) ∼ ∑_{j=1}^{∞} ∑_{k=0}^{∞} A^{k}B^{j} ∑_{G_{jk}^{conn}} L_{jk}[G].

Here G_{jk}^{conn} is the set of all connected multigraphs with loops, G, with j vertices and k lines (links plus loops), and L_{jk}[G] is a certain expression involving j factors V (z_{i}) and k operators ∇_{z_{i}} · ∇_{z_{i′}} . A corollary is that

e^{L} ∼ 1 + ∑_{j=1}^{∞} ∑_{k=0}^{∞} A^{k}B^{j} ∑_{G_{jk}^{all}} L_{jk}[G]

— the same formula except that the condition of connectedness is dropped (and a null graph included at the beginning).

Here is an example of a connected multigraph with 3 vertices, 4 links, and 1 loop:

•  ^{◯} • {⌢ \atop ⌣} •

Each vertex represents a V. Each line represents a “contracted” (in the sense of tensor calculus) pair of differentiations acting on those potentials. Let us write out the term denoted by a simpler example graph:

L_{21}[• •] = \frac{1}{2} ∫_{0}^{1} dξ_{1} ∫_{0}^{1} dξ_{2} ξ_{<}(ξ_{>} − 1) ∇V (z_{1}) · ∇V (z_{2}),

where z_{i} = y +ξ_{i}(x−y) (the parametrized line segment from y to x) and ξ_{<} = min(ξ_{1}, ξ_{2}), ξ_{>} = max(ξ_{1}, ξ_{2}).

Connectedness of a graph with few lines implies that it cannot have many vertices: j ≤ k − 1. Therefore, in ∑_{j,k,G_{jk}^{conn}} A^{k}B^{j}L_{jk}[G] only finitely many terms contribute to each order in A! This miracle allows what started as a coupling-constant expansion to be rearranged into a large-mass expansion.

The bad news is that since A = \frac{h̄t}{m} , this expansion is nonuniform in t. No matter how small \frac{h̄}{m} is, the approximation will be accurate only for sufficiently small t. Because \frac{A}{B} = \frac{h̄^{2}}{mλ} is independent of t, an expansion in that parameter would be better. However, it can’t be obtained by rearranging the existing series, since that would require summing all terms with a fixed value of k − j, and there are infinitely many such.

Nevertheless, a WKB (\frac{A}{B} ) series can be constructed, but it requires integration along the exact orbits of the underlying classical-mechanical system, instead of along straight lines. In fact, Osborn and Molzahn have shown [20, 24, and especially 25] that the A series can be recovered from that one by expanding the exact orbits about the straight orbits of the free system; the quantity ξ_{<}(ξ_{>} −1) appearing in the L_{jk} is the Green function solving this classical perturbation problem!

It is not surprising that such an approximation becomes bad at times so large that the true orbits have begun to intersect each other (caustics). Many other researchers have had things to say about this problem much more profound than I have to offer [e.g., 6, 9, 11, 14, 17, 31]. I would, however, like to emphasize a point that evidently is not generally appreciated.

A generic semiclassical approximation is of the form

ψ = A e^{iS/h̄} + O(h̄),

where S is the classical action of a path, and A is the −\frac{1}{2} power of a certain determinant that vanishes when the classical paths intersect. When one is calculating an eigenfunction or a resolvent kernel, the fact that A diverges on the caustic (or “turning point”) is quite properly interpreted as a breakdown of the WKB approximation there. (The true solution must be smooth and finite on the caustic.) The story for time-dependent situations is different, because the Schrödinger equation is not hypoelliptic — its solutions don’t have to be smooth for all time. For the harmonic oscillator, V = \frac{1}{2}ω^{2}x^{2}, the WKB solution for the propagator is known to be exact:

K(t, x, y) = K_{WKB}  (for 0 ≤ ωt ≤ π)

= \sqrt{\frac{mω}{2πih̄ sin (ωt)}} exp \{\frac{imω}{2h̄ sin (ωt)}[(x^{2} + y^{2}) cos (ωt) − 2xy]\}

(known as Mehler’s formula). This means that the infinity that appears at ωt = Nπ is genuine, not an error in the approximation!

This raises the issue of what happens in the propagator for more general potentials. Zelditch [36–37] has proved that if the free Hamiltonian or the harmonic oscillator is perturbed by a smooth function of compact support (the support condition being fundamental!), the singularity set of the exact propagator doesn’t change. The semiclassical interpretation of this theorem seems to be that, although the perturbation may cause paths to cross in the free case, and cause the caustic to spread out to new places in space-time in the oscillator case, nevertheless “most” paths pile up at exactly the same caustic as before. Further progress in rigorously relating the singularities of the propagator to the large-time qualitative behavior of the classical orbits has been made recently by Craig, Kappeler, and Strauss [4].

Davin Potts [30], an undergraduate student working under my direction, programmed the calculation of the first few terms in the Wigner–Kirkwood series and compared them with the exact propagator. Algebraically, these terms agree with a direct expansion of the Mehler formula. Numerically, we found excellent agreement with the exact formula for ωt ≤ 2, but terrible results at ωt = 3 and beyond. This is no surprise, since the harmonic oscillator has a particularly drastic caustic at ωt = π: all the orbits refocus at one point then! (It should be noted that for this system, because the potential is a homogeneous function, the only independent parameter is ωt; there is no issue of choosing the mass, etc., to have physically realistic values.)

We also considered the inverted harmonic oscillator, V = −\frac{1}{2}ω^{2}x^{2}, for which there is no caustic. As expected, the results [30] are noticeably better than for the harmonic oscillator at the same ωt, but far from perfect as soon as ωt becomes comparable to π.

The Wigner–Kirkwood expansion has an obvious calculational advantage over the more difficult WKB approximation, so it is important to determine whether its domain of applicability can be extended. Two gambits are under investigation (both of which have been previously discussed by Makri and Miller [15–16]); (1) turning the power series into a modified Padé approximant [35] that can reproduce the singularity of a true propagator;

(2) iterating the action of an approximate short-time propagator to get a good mediumtime approximation (equivalent to discretization of the Feynman path integral). I hope that these methods can yield reasonable approximations over a large number of periods of oscillation of a system. Each step of iteration requires an integration over configuration space. I believe that this computationally expensive step can have its frequency minimized by applying the method known in numerical analysis as “Richardson extrapolation” [e.g., 1, 12, 32]. (Incidentally, according to Richardson himself [32], part of the early development of that method as applied to differential equations was done here in Kiev [3].) This project has barely begun, so there are no results to report.

I should emphasize that this trouble does not arise in the case of a pure gravitational field, H = −\frac{h̄^{2}}{2m} g^{μν}∇_{μ}∇_{ν} [e.g., 18]. Here there is no distinction between h̄ expansion and \frac{1}{m} expansion. Admittedly some of the computational simplicity of the latter is now lost, since one must integrate along geodesics instead of straight lines [23–24]. But solving for all geodesics on a manifold is much easier than solving for all classical orbits with a nongravitational force included. (Because of the “principle of equivalence”, the manifold of paths has a smaller effective dimension in the purely gravitational case.)

To close, I want to return to the simple Schrödinger equation (in arbitrary dimension, however) and to the title of my talk. Recall that the effect of the potential is expressed by the graph sum

e^{L} ∼ 1 + ∑_{j=1}^{∞} ∑_{k=0}^{∞} A^{k}B^{j} ∑_{G_{jk}^{all}} L_{jk}[G].

For example, in the case j = 2, k = 1, y = x there are three terms:

• •^{◯}  and • •^{◯}  both yield − \frac{i}{24} V ∇^{2}V ,

• • yields − \frac{i}{24} |∇V |^{2}.

This illustrates a shortcoming of the formulation as stated so far: Isomorphic labeled graphs give identical contributions, so to reduce the formula to a sum of linearly independent terms, we must sum over unlabeled graphs instead. (I.e., graphs that differ only in the order of their vertices are equivalent.)

The classification of graphs by isomorphism is a notorious combinatorial problem, believed to be either NP-complete or just slightly less bad [13, 34]. (On the other hand, Goldberg [8] has shown that the enumeration of (simple) graphs can be done by an algorithm with “polynomial delay”.) For the relatively low-degree graphs of practical interest, however, the problem is tractable. The fundamental combinatorial method known as Pólya’s Theorem [10, 28, 29] can be used to classify graphs by vertex spectrum and simultaneously by line spectrum. That is, we can calculate how many multigraphs with j vertices and k lines contain exactly 2 vertices that have exactly 3 lines attached to them and also have exactly 1 pair of vertices with exactly 2 lines between them, etc. It is then a simple matter to find all the graphs in each such category, and thereby compile a complete catalog of all unlabeled multigraphs of size (j, k).

For example, here is a table of all multigraphs with 2 vertices and 3 lines. The columns are labeled by the vertex spectra, and the rows are labeled by the line spectra. Note that one category has two representatives, and several categories have none.

This graph-enumeration work is in progress in collaboration with I. Borosh of Texas A&M University, with some help by a student, A. da Conturbia. We hope to have a

6  51  42  33
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3  _{•}^{◯}_{◯}_{◯} _{•}  _{•} _{•} {⌢ \atop ⌣}
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21  _{•}^{◯}_{◯} _{•}  _{•}^{◯}_{◯} _{•}^{◯}
• •^{◯}{⌢ \atop ⌣}
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

111 _{•}^{◯} _{•}^{◯}

publication by the end of this summer. I believe that with this final step we can claim to have a complete (and elegant) solution of the \frac{1}{m} asymptotics of the Schrödinger equation. In view of recent work of Molzahn and Osborn and coworkers [2, 5, 19, 22, 27], the same methods apply to h̄ asymptotics, to more complicated potentials (e.g., electromagnetic), and to semiclassical solution of the Heisenberg equations of motion.

Work with T. Osborn and F. Molzahn supported by a NATO Grant for Scientific Cooperation. Work with D. Potts supported by a Faculty Research Minigrant from Texas A&M University and by the Undergraduate Fellows Program of TAMU. Computations done on a Sun 3/60 workstation with Mathematica software, provided by a grant from Sun Microsystems extended by Wolfram Research.

I am grateful to J. N. Lyness for pointing me toward reference [12] and to V. V. Kornyak and the librarians of the Institute of Mathematics of the Ukrainian Academy of Sciences for tracking down reference [3]. (Honest scholarship requires me to add that [3] appears to be only marginally relevant to Richardson extrapolation as a computational technique. Bogolubov and Krylov were concerned with deriving an error estimate for finite-difference approximations to solutions of boundary-value problems. In a footnote Krylov thanks Richardson for suggesting the problem while the two were crossing the Atlantic to a conference in Canada. Richardson’s citation in [32] (reported by Joyce in [12]) reads, “The deferred approach to the limit has also been considered by N. Bogolouboff and N. Kryloff in a recent paper [[3]], in Russian.” The document is actually in Ukrainian; one can speculate that Richardson did not know precisely which part of the shipboard conversation had borne fruit in it.)

References

1. F. S. Acton, Numerical Methods that Work, Harper and Row, New York, 1970.

2. A. O. Barvinsky, T. A. Osborn, and Yu. V. Gusev, J. Math. Phys. 36, 30–61 (1995).

3. N. Bogolouboff and N. Kryloff, On the Rayleigh’s Principle in the Theory of the Differential Equations of the Mathematical Physic and upon the Euler’s Method in the Calculus of Variations, Ukr. Akad. Nauk, Trudi Fyz.-Mat. Viddilu, Vol. 3, No. 3, 1926 (in Ukrainian).

4. W. Craig, T. Kappeler, and W. Strauss, Microlocal Dispersive Smoothing for the Schrödinger Equation, preprint, Brown University.

5. R. A. Corns, J. Phys. A 27, 593–607 (1994).

6. J. B. Delos, Adv. Chem. Phys. 65, 161–214 (1986).

7. Y. Fujiwara, T. A. Osborn, and S. F. J. Wilk, Phys. Rev. A 25, 14–34 (1982).

8. L. A. Goldberg, Efficient Algorithms for Listing Combinatorial Structures, Cambridge University Press, Cambridge, 1993.

9. M. C. Gutzwiller, J. Math. Phys. 10, 1004–1020 (1969).

10. F. Harary and E. Palmer, Graphical Enumeration, Academic Press, New York, 1973.

11. E. J. Heller and S. Tomsovic, Phys. Today 46 (7), 38–46 (1993).

12. D. C. Joyce, SIAM Rev. 13, 435–490 (1971).

13. L. Kučera, Combinatorial Algorithms, Hilger, Bristol, 1990.

14. R. G. Littlejohn, J. Math. Phys. 31, 2952–2977 (1990).

15. N. Makri, Computer Phys. Commun. 63, 389–414 (1991).

16. N. Makri and W. H. Miller, J. Chem. Phys. 90, 904–911 (1989).

17. V. P. Maslov and M. V. Fedoriuk, Semi-Classical Approximation in Quantum Mechanics, Reidel, Dordrecht, 1981.

18. S. Minakshisundaram and Å. Pleijel, Canad. J. Math. 1, 242–256 (1949).

19. F. H. Molzahn, J. Phys. A 25, 4913–4940 (1992); erratum ibid. 26, 2275–2276 (1993).

20. F. H. Molzahn and T. A. Osborn, J. Math. Phys. 27, 88–99 (1986).

21. F. H. Molzahn and T. A. Osborn, J. Phys. A 20, 3073–3094 (1987).

22. F. H. Molzahn and T. A. Osborn, Ann. Phys. (N.Y.) 230, 343–394 (1994).

23. F. H. Molzahn, T. A. Osborn, and S. A. Fulling, Ann. Phys. (N.Y.) 204, 64–112 (1990).

24. F. H. Molzahn, T. A. Osborn, and S. A. Fulling, Ann. Phys. (N.Y.) 214, 102–141 (1992).

25. T. A. Osborn and F. H. Molzahn, Phys. Rev. A 34, 1696-1707 (1986).

26. T. A. Osborn and F. H. Molzahn, in Forty More Years of Ramifications: Spectral Asymptotics and Its Applications, ed. by S. A. Fulling and F. J. Narcowich (Discourses in Mathematics and Its Applications, No. 1), Texas A&M University Department of Mathematics, College Station, 1993, pp. 199–236.

27. T. A. Osborn and F. H. Molzahn, Ann. Phys. (N.Y.) 241, 79–127 (1995).

28. K. R. Parthasarathy, Can. J. Math. 20, 40–47 (1968).

29. G. Pólya and R. C. Read, Combinatorial Enumeration of Groups, Graphs, and Chemical Compounds, Springer, New York, 1987.

30. D. M. Potts, Large Mass Approximations in Quantum Physics and A Bridge to Quantum Chemistry, Undergraduate Research Fellow Thesis, Texas A&M University, 1995.

31. P. J. Reynolds, ed., Special Issue on Quantum Chaos, J. Stat. Phys. 68, Nos. 1/2 (1992).

32. L. F. Richardson, Phil. Transac. Roy. Soc. London A 226, 299–349 (1927).

33. A. Saksena, T. A. Osborn, and F. H. Molzahn, J. Math. Phys. 32, 938–955 (1991).

34. S. Skiena, Implementing Discrete Mathematics, Addison–Wesley, Redwood City, 1990.

35. D. F. Styer, Computer Phys. Commun. 61, 374–386 (1990).

36. A. Weinstein and S. Zelditch, Bull. Amer. Math. Soc. 6, 449–452 (1982).

37. S. Zelditch, Commun. Math. Phys. 90, 1–26 (1983).