Semiclassical propagation of Wigner functions

Semiclassical propagation of Wigner functions

T. Dittrich, E. A. Gómez, L. A. Pachón Departamento de Física, Universidad Nacional de Colombia, Bogotá D.C., Colombia.
CeiBA – Complejidad, Bogotá D.C., Colombia.
July 16, 2019

We present a comprehensive study of semiclassical phase-space propagation in the Wigner representation, emphasizing numerical applications, in particular as an initial-value representation. Two semiclassical approximation schemes are discussed: The propagator of the Wigner function based on van Vleck’s approximation replaces the Liouville propagator by a quantum spot with an oscillatory pattern reflecting the interference between pairs of classical trajectories. Employing phase-space path integration instead, caustics in the quantum spot are resolved in terms of Airy functions. We apply both to two benchmark models of nonlinear molecular potentials, the Morse oscillator and the quartic double well, to test them in standard tasks such as computing autocorrelation functions and propagating coherent states. The performance of semiclassical Wigner propagation is very good even in the presence of marked quantum effects, e.g., in coherent tunneling and in propagating Schrödinger cat states, and of classical chaos in four-dimensional phase space. We suggest options for an effective numerical implementation of our method and for integrating it in Monte-Carlo–Metropolis algorithms suitable for high-dimensional systems.

03.65.Sq, 31.15.Gy, 31.15.Kb

I Introduction

Molecular dynamics, by the spatial and temporal scales it involves, straddles the borderline between quantum and classical behavior. Quantum effects—above all the very concept of chemical reaction—are obviously crucial. Notwithstanding, a good deal of a molecule’s external and even internal motion can be understood on purely classical grounds: precisely the situation for which semiclassical methods have been conceived and are optimally suited. Moreover, recent developments in chemical physics, experimental as well as theoretical, are advancing vigorously in the direction of complex time evolution and high excitations. This faces adiabatic techniques like even time-dependent density-functional theory (TDDFT) RG84 () and time-dependent Hartree GBR82 (); MMC90 () with formidable challenges while it favors semiclassical approaches which do not suffer from limitations in this respect.

Therefore, the modest renaissance semiclassics in the time domain (for a comprehensive review see SG96 ()) is presently enjoying comes as no surprise. It includes approaches in configuration space as well as less established phase-space methods. The well-known shortcomings plagueging semiclassics in coordinate representation (based on WKB approximations BM72 () and the van-Vleck–Gutzwiller propagator Vle28 (); Gut67 (); Lit92 ()), such as divergences at caustics and the root-search problem, have largely been overcome by sophisticated refinements of the original method Mil74 (). Remarkably, most of them already switch internally to mixed representations (combining position with momentum) or full phase space in order to smooth out divergences and to reduce the set of contributing classical solutions to the neighbourhood of a single trajectory, permitting the construction of initial-value representations (IVRs): uniform approximations LS77 (); LM&78 (); Kla86 (), semiclassical IVRs Mil74 (); Mil98 (); TWM01 (); LM06 () (for reviews see Mil01 (); TW04 (); Kay05 ()), Gaussian-wavepacket propagation Hel75 () with its numerous ramifications, notably the Herman-Kluk propagator HK84 (); Gro99 (); BA&01 () and Heller’s cellular dynamics Hel91b (). A similar but independent approach, the dephasing representation Van04 (); Van06 (); LMV09 () employs the shadowing lemma for chaotic dynamics to represent the classical flow through a Planck cell by a single trajectory in a semiclassical approximation for correlation functions.

These methods have gained widespread acceptance in practical applications but tend to be technically cumbersome. Their implementation in ab initio molecular-dynamics simulation with “on the fly” update of the electronic evolution MTM99 (); DM02a () has been proposed CA&09 (). It has to cope, though, with the high computational costs of calculating Hessians for the potential-energy surface of the electron sector—so much so that recourse is sought even at a complete suppression of determinantal prefactors TP09 ().

Semiclassical approximations working directly in phase space promise a fresh, structurally more transparent approach. They are inherently free of the above problems and thus offer a number of tempting advantages: (i) phase-space propagation provides a natural exact IVR by construction, without Gaussian smoothing, (ii) it can be formulated exclusively in terms of canonically invariant quantities (no projection required) that (iii) allow for simple geometrical interpretations, (iv) as far as determinantal prefactors arise, they are also canonically invariant, and (v) being based on the density operator, not on wavefunctions, the extension to non-unitary time evolution is immediate, opening access to decoherence and dissipation ORB09 ().

Yet there is a price to be paid. The two representations mainly considered as candidates for semiclassical phase-space propagation both come with their specific virtues and vices: The Husimi function Hus40 () is based on coherent states KS85 (), hence closely related to semiclassical IVRs and Gaussian-wavepacket propagation. Coherent states possess a natural interpretation in terms of classical probabilities but are overcomplete, contain arbitrary parameters, and, in order to work optimally, require to complexify phase space Kla79 (); BA&01 (). The Wigner function Wig32 (); HO&84 (); Sch01 (), by contrast, provides a parameter-free one-to-one representation of the full density operator but encodes information on quantum coherences as small-scale (“sub-Planckian”) oscillatory fringes Ber77 (); Zur01 (). They prevent a probability interpretation and pose serious problems for their propagation, pointed out more than 30 years ago by Heller Hel76 ().

Owing to these seemingly prohibitive difficulties, the semiclassical propagation of Wigner functions still awaits being explored in its full potentiality. For a long time, it has been almost synonymous to a mere classical propagation of phase-space distributions, discarding quantum effects in the time evolution altogether Mcl83 (); GH95 (). As a first attempt to improve on this so-called classical Wigner model, it suggests itself to include higher-order terms in in the quantum analogue of the Poisson bracket, the Moyal bracket Ber77 (), upon integrating the evolution equation for the Wigner function. Such additive quantum corrections give rise to accordingly modified “quantum trajectories” LS82 (); Lee92 (); Raz96 (); DM01 (); THW03 (). However, they tend to become unstable even for short propagation time MP09 () and suffer from other practical and fundamental problems MS96 (): As we shall point out below, the very concept of propagating along single deterministic trajectories is inadequate even in the semiclassical regime.

A more convincing alternative has emerged in the context of quantum chaos where semiclassical Wigner functions Ber77 (); Hel77 (); Ozo98 () are appreciated as valuable tools, e.g., in the study of scars in wavefunctions and of spectral features related to them Ber89 (); AF98 (); TAO01 (). By applying semiclassical approximations directly to the finite-time propagator, expanding the phase instead of the underlying evolution equation, an appropriate phase-space propagation method for semiclassical Wigner functions has been achieved RO02 (). It features, in a geometrically appealing manner, pairs of (real) classical trajectories and the symplectic area enclosed between them as their classical underpinning—not surprisingly in view of the well-known pairing of paths and diagrams, one forward, one backward in time, as it occurs in the propagation of quantities bilinear in the state vector, particularly the density operator.

While this technique is tailor-made to time-evolve the specific class of semiclassical Wigner functions Ber77 (); Hel77 (); Ozo98 (); TAO01 (), semiclassical approximations to the propagator proper, independent of the nature of initial and final states, have been constructed DVS06 () on basis of phase-space path integrals Mar91 () and of the Weyl-transformed van Vleck propagator Ber89 (). This approach opens the door towards a far broader range of applications. It has already proven fruitful in an analysis of the spectral form factor of classically chaotic systems in terms of phase-space manifolds DP09 ().

Here, by difference, we pretend to demonstrate its viability as a practical instrument for the propagation of molecular systems, including in particular numerical simulations. We therefore focus on its performance in benchmark models established in chemical physics like the Morse oscillator and the quartic double well and study standard tasks such as propagating Gaussians and computing autocorrelation functions. We are confident that the remarkable accuracy and robustness of the scheme and its unexpected numerical stability and efficiency suggest it as a new competitive option in semiclassical propagation.

In Sec. II, we introduce the basic building blocks of our method, specify our two main approaches, via the van-Vleck propagator and via phase-space path integration, and discuss general properties of the semiclassical Wigner propagator, such as its geometrical structure and its asymptotics towards the classical limit and towards weak anharmonicity. This section partially corresponds to a more detailed account of material already published in Ref. DVS06 (). Readers mainly interested in numerical results might proceed directly to Sec. III, dedicated to a broad survey of the performance of semiclassical Wigner propagation for prototypical models of nonlinear molecular potentials, in the computation of dynamical functions like autocorrelations as well as of detailed structures in phase space. We also test our method in particularly challenging situations, involving strong quantum effects, such as the reproduction of coherent tunneling and the propagation of Schrödinger-cat states, or a classically chaotic dynamics. Strategies to optimize the implementation of this method in numerical, including Monte-Carlo-like, algorithms are suggested in Sec. IV. We conclude in Sec. V with an outlook to immediate extensions of our work.

Ii Constructing the semiclassical Wigner propagator

ii.1 Definitions and basic relations

Alternative to its definition as a Weyl-ordered transform of the density operator HO&84 (), the Wigner function can be introduced as the expectation value of phase-space reflection operators Ozo98 () or as expansion coefficient function in a basis of phase-space displacement operators GT88 (), besides other options. For our purposes, the standard definition suffices,


where denotes the density operator and is a vector in -dimensional phase space (we adopt this ordering throughout the paper but assign to the horizontal and to the vertical axis in phase-space plots).

We shall restrict ourselves to unitary time evolution within the dynamical group generated by a time-independent Hamiltonian . The extension to time-dependent Hamiltonians is immediate but will be suppressed to avoid technicalities. The von-Neumann equation for the density operator, , translates into an equation of motion for the Wigner function, , involving the Weyl symbol of the Hamiltonian . For , the Moyal bracket HO&84 (); Sch01 () converges to the Poisson bracket.

The time evolution of the Wigner function over a finite time can be expressed as an integral kernel, the Wigner propagator ,


It forms a one-parameter group, implying in particular the initial condition and the composition law


An important quantity that must not be confused with the Wigner propagator is the Weyl transform of the evolution operator , referred to as the Weyl propagator for short Ozo98 (),


It cannot be used to propagate Wigner functions as it stands, but is related to the Wigner propagator by a self-convolution Mar91 (); RO02 (); DVS06 (),


with . It serves as a starting point for the construction of semiclassical approximations, based on the van-Vleck propagator (Sec. II.2) as well as on phase-space path integration (Sec. II.3).

ii.2 van Vleck approximation for the Wigner propagator

A straightforward route towards a semiclassical Wigner propagator is replacing the Weyl propagator in Eq. (5) by the Weyl transform of the van Vleck propagator Ber89 (); Mar91 (); Ozo98 (). Transformed from the energy to the time domain, it reads Mar91 (); Ozo98 (),


The sum runs over all classical trajectories connecting phase-space points to in time such that (midpoint rule). and are its stability matrix and Maslov index, respectively. The action , with , the Weyl Hamiltonian evaluated on the trajectory (to be distinguished from ) and , the symplectic area enclosed between the trajectory and the straight line (chord) connecting to Ber89 () (chord rule, vertically hashed areas in Fig. 1).

Figure 1: Symplectic areas entering the van-Vleck-based semiclassical Wigner propagator, Eqs. (10, II.2). The vertically hashed areas correspond to the phases of the Weyl propagators (6) according to the chord rule. The symplectic area (slanted hatching) enclosed between the two classical trajectories and the two transverse vectors and determines the phase (II.2) of the propagator. The classical trajectory (dashed) is to be distinguished from the propagation path (bold) connecting the initial argument of the propagator to the final one, .

Substituting Eq. (6) in (5), one arrives at


where indices refer to classical trajectories contributing to the Weyl propagators in Eq. (5). The principal challenge is now evaluating the -integration. As it stands, Eq. (II.2) couples the two classical trajectories , , to one another only indirectly via , the separation of their respective midpoints. This changes as soon as an integration by stationary-phase approximation is attempted, consistent with the use of the van Vleck propagator. Stationary points are specified implicitly by the condition . Combined with the midpoint rule , this implies


Equation (8) constitutes a simple geometrical rule for semiclassical Wigner propagation Ozo98 (): It is based on pairs of classical trajectories , , that need not coincide with one another nor with the trajectories passing through and but must have midway between their respective initial points and likewise for . This condition is trivially fulfilled for identical pairs , as a type of diagonal approximation. By including also non-diagonal terms, i.e., non-identical trajectory pairs, we add classical information contained in third-order terms in the action and partially compensate for the general limitations of the van Vleck approach on the level of wavepackets. This becomes particularly evident in the case of pure initial states : Applying the van Vleck propagator in position representation separately to ket and bra, allowing for off-diagonal terms in the ensuing double sum over classical orbits, and then transforming back to the final Wigner function leads to the same result as the present derivation, albeit through a more tedious and less transparent route.

To complete the Fresnel integral over , we note that


where denotes the symplectic unit matrix Ber89 (). Combined with the determinantal prefactors inherited from the van Vleck propagator, this produces


our main result for the semiclassical Wigner propagator in van Vleck approximation. The phase is determined by


with and . Besides the two Hamiltonian terms it includes the symplectic area enclosed between the two trajectory sections and the vectors and (Fig. 1).

Figure 2: Classical skeleton of the semiclassical Wigner propagator according to Eqs. (10,II.2). A set of initial points (panel a, left target pattern) surrounding the initial point of the propagation path is propagated along classical trajectories (dotted lines), giving rise to a corresponding set of final points (right deformed target pattern). The manifold formed by their midpoints (cone-like structure, panel b) constitutes the “illuminated area”, the support of the semiclassical propagator. Inside this area, two trajectory pairs contribute to each point, one of them elliptic (“upper” shell of the cone, blue), the other hyperbolic (“lower” shell, red). Its boundary gives rise to caustics in phase space where the two trajectory pairs collapse into one. Panel c depicts the generic structure of the underlying phase as a function of the integration variable in Eq. (II.2), an odd cubic polynomial of the phase-space coordinates with two extrema, a minimum and a maximum (blue symbols) and two saddles (red symbols), corresponding to the elliptic and the hyperbolic trajectory pairs, respectively.

In the following we list a number of general features of Eqs. (10,II.2):

  1. Equation (10) replaces the Liouville propagator,


    localized on the classical trajectory initiated in , by a “quantum spot”, a smooth distribution peaked at the support of the classical propagator but spreading into the adjacent phase space and structured by an oscillatory pattern that results from the interference of the trajectories involved.

  2. The propagator (10) does involve determinantal prefactors. However, they do not result from any projection onto a subspace like or and are invariant Ber89 () under linear canonical (affine) transformations Bal80 ().

  3. It deviates from the Liouville propagator only if the potential is anharmonic. For a purely harmonic potential, the two operations, propagation in time and forming midpoints between trajectories, commute, so that all midpoint paths coincide with each other and with the classical trajectory . This singularity restores the classical delta function on , see Sec. II.4.3.

  4. By contrast to position-space semiclassics Lit92 (), the number of trajectory pairs contributing to the summation in Eq. (10) ranges between 0, 1, and 2 following a universal pattern (Fig. 2): Within an “illuminated area”, a sector with its tip on the classical trajectory, the sum contains two trajectory pairs (four stationary points, Fig. 2c). In the “shadow region” outside this sector, stationarity cannot be fulfilled by real trajectories. Along the border there is exactly one solution (two stationary points). Unexpected in phase space, this caustic arises as the projection onto phase space of the manifold of midpoints (Fig. 2b).

  5. The set of trajectory pairs to enter the calculation of the propagator is only restricted by the midpoint rule (8). However, this does not constitute a double-sided boundary condition since every pair fulfilling Eq. (8) for the initial points contributes a valid data point to the propagator: no root-search. The freedom in the choice of trajectory pairs can be exploited to optimize algorithms, see Sec. IV.1.3.

  6. Of the two trajectory pairs contributing to the illuminated area, exactly one is elliptic, associated to a pair of opposite extrema of the action, one is hyperbolic and associated to a pair of saddles (Fig. 2c). They can be distinguished by monitoring , positive for the former and negative for the latter. The two sets form separate sheets of the action (II.2). Within each of them, amplitude and phase of the propagator are slowly varying functions of , facilitating their numerical treatment, see Sec. IV.

  7. The propagator’s oscillatory pattern encodes information on quantum coherences and allows us to propagate the “sub-Planckian oscillations” that characterize the Wigner function, solving the problem of the “dangerous cross terms” pointed out by Heller Hel76 (). See Sec. III.5.

  8. The principle of propagation by trajectory pairs is consistent with the properties of a dynamical group. It translates the concatenation of propagators into a continuation of trajectories, if the convolution (3) is evaluated by stationary phase as well.

  9. Equations (10, II.2) fail if the stationary points approach each other too closely. This is the case near the central peak of the propagator on the classical trajectory and along the caustics mentioned above. The same problem arises in the limits , , and of weak anharmonicity. It can be overcome by means of a uniform approximation to the -integration in Eq. (II.2), see Sec. II.4.2.

  10. As a consequence of the underlying van Vleck approximation, the propagator (10) is not properly normalized. Specifically, the slow decay of its oscillatory tail renders it nonintegrable with respect to its initial or final arguments, see Sec. III.3.1 for quantitative details. The problem is solved by the same uniform approximation as announced above.

ii.3 Path-integral approach

The quality of the van Vleck approximation (10) to the Wigner propagator is limited by its principal ingredient, the van Vleck propagator itself. Just as in configuration-space propagation, a comprehensive solution is provided on the basis of path integrals Fey48 (); Sch81 (). In phase space an analogous theoretical framework is available since the seminal work of Marinov Mar91 (), largely following Feynman: Time is discretized into equidistant sections , where and is the time span to be propagated over. At each , two short-time propagators, composed of the corresponding Weyl propagators according to Eq. (5), are concatenated by means of Eq. (3). In the continuous limit , , the full propagator is obtained as a double path integral comprising two distinct path variables and , inherited, respectively, from Eqs. (3) and (5),


with an action integral


Only is subject to boundary conditions and while is free and may be associated to quantum fluctuations.

For this subsection, we restrict ourselves to a single degree of freedom, , and to standard Hamiltonians , , but admit an explicit time dependence and an arbitrary degree of nonlinearity of the potential. In this case, the Weyl Hamiltonian, , is obtained by replacing the operators , in with corresponding real variables, .

Returning to discrete time, we are faced with evaluating the action


denoting and . To obtain a semiclassical approximation, we expand the action (15) in around a classical trajectory defined implicitly in discrete time by , , and around , . By expanding around a single trajectory, we are sampling less classical information than in the van Vleck approximation based on trajectory pairs. This renders our path-integral approach inferior in some respects detailed below. Note that the action (15) is not only an odd polynomial in , it is even lacking a linear term, hence , since according to Hamilton’s equations of motion, the first term cancels exactly the leading order of a Taylor expansion of with respect to .

Expanding to third order in all variables (there is no -term in anyway), this leaves us with


The integrations over the can now be performed, much like the -integrations in the derivation of the Feynman path integral Fey48 (); Sch81 (), resulting in one-step propagators


where we use shorthands and and define . The factor is analogous to the determinantal prefactors occurring in position-space path integration in that it involves higher derivatives of the potential. Similarly as in the van-Vleck-based Eq. (10), however, it constitutes a mathematically more benign canonical invariant.

Taking into account that the one-step propagators (17) are already written in a frame moving with the classical trajectory, we observe three independent actions of the propagator on the evolving quantum spot: (i), a rigid shift along the classical trajectory, (ii), a rotation (shear) with the linearized phase-space flow around the classical trajectory, if it is elliptic (hyperbolic), and (iii), quantum Airy-function spreading in the negative -direction if is positive and vice versa. Their concatenation is most conveniently performed in Fourier space, where, under certain conditions fulfilled here fn1 (), the convolution (3) reduces to a mere multiplication. With rescaled coordinates of common dimension ,


the transformation of the Wigner propagator to Fourier phase space is given as


Going to the continuum limit, we thus obtain


Here, is the stability matrix of , the subscript selects the second component of , and


Inverting the Fourier transform (19) we recover the Wigner propagator. In the following subsection, we work out this result in detail for the specific case .

ii.4 Weak anharmonicity, short-time and classical limits

ii.4.1 Path-integral approach for weak cubic nonlinearity

Equation (20) permits an analytical evaluation in the case of constant second- and third-order derivatives of the potential, which we will discuss in this subsection. It can arise in different scenarios: To begin with, consider a static binding potential with cubic nonlinearity, e.g., , and a trajectory close to the quadratic minimum at , so that the linearized flow is elliptic. Then and we can also assume . In accordance with Eq. (21), denote . In this case, the stability matrix corresponds to a rotation by ,


and the phase of the Fourier transformed propagator (20) reduces to elementary integrals,

In a reference frame rotating by with respect to , , , the propagator after inverting the Fourier transform takes a particularly transparent form


where we abbreviate , , , and . The corresponding result for the hyperbolic case , e.g., close to a quadratic maximum of a cubic potential, is obtained from Eq. (23) by replacing and .

Equation (23) provides us with a precise geometric outline that complements the features of the Wigner propagator pointed out in the sequel of Eq. (II.2):

  1. The quantum spot formed by the propagator fills a sector in phase space of opening angle (Fig. 3c), rotating by if is the angle coordinate along the classical trajectory. Its nodes form straight lines running parallel to its sidelines, with separations given by Airy-function zeros (Fig. 3).

  2. For short scaled time , the oscillatory tail is quasi-one-dimensional, , pointing in the negative -direction if and vice versa, and scales with time as . It is periodic in with period and invariant under time reversal, , , with respect to .

  3. As expected for a uniform approximation, Eq. (23) resolves the sharp caustics along the outlines of the quantum spot present in the van Vleck approximation (10) (Fig. 3b) into a smooth penumbra (Fig. 3c,d). This restores the correct normalization of the propagator. On the other hand, it restricts the oscillatory pattern to straight nodelines while with Eq. (10), nonlinear deformations owing to higher nonlinearities in the potential can well be reproduced, see Sec. III.2.

  4. The finite weight the propagator (23) achieves in the penumbra zone outside the illuminated region of Eq. (10) can be interpreted as the contribution of complex trajectories. However, we shall not pursue this issue here.

ii.4.2 Perturbation theory from van Vleck approach

A uniform approximation equivalent to Eq. (23) can be obtained through an alternative route which is particularly helpful in order to analyze the asymptotics at weak anharmonicity and at short time and the classical limit for the Wigner propagator: We treat the weak cubic nonlinearity explicitly as a perturbation of an otherwise quadratic Hamiltonian


to obtain expressions for the various ingredients of the van Vleck propagator (10) through a perturbation expansion in .

From the perturbed orbits , we obtain pairs of trajectories with initial points displaced by from , the corresponding centers , and the deviations of these final midpoints from the endpoint of the perturbed trajectory starting in ,


with time-dependent coefficient matrices


Arising from a perturbation expansion of the equations of motion, Eqs. (25, 26) are reliable only up to leading order in , hence their validity is also restricted in time to . They constitute the classical skeleton of the quantum spot and provide the input for its determinantal prefactor and action.

It is, however, more convenient to express these quantities in terms of , the final argument of the propagator relative to the classical trajectory starting in , instead of the initial displacement . By means of a few elementary phase-space transformations that diagonalize simultaneously the two coefficient matrices (26), the pair of quadratic equations (25) can be resolved for as a function of . We thus obtain the propagator in the form


with and . We identify the cosine term in Eq. (II.4.2) as the contribution of the pair of hyperbolic trajectories (saddles of the action, Fig. 3e) and the sine term as the elliptic contribution (extrema, Fig. 3f), cf. remark (vi) after Eq. (II.2) and Fig. 2c.

Unfortunately, Eq. (II.4.2) becomes spurious precisely in the limit for which it has been devised, since in this limit, the stationary points of the action coalesce and the stationary-phase approximation underlying Eq. (II.4.2) breaks down. The problem can be fixed, however, in a straightforward manner: Since the action in Eq. (II.2) is odd in the integration variable and the two known trajectory pairs entering Eq. (II.4.2) mark its four stationary points and (Fig. 2c), corresponding resp. to the elliptic and hyperbolic pairs, we can uniquely reconstruct an effective action that contains only linear and cubic terms in (Fig. 2c) and thus permits an exact evaluation of the -integration. The result (cf. Fig. 3c), replacing Eq. (II.4.2), reads


with and , as above. It follows from Eq. (23) as a short-time approximation, to leading order in . Equation (II.4.2), in turn, is recovered if we replace the Airy functions by their asymptotics for large negative argument AS84 (), .

Figure 3: Different versions of the Wigner propagator for a harmonic potential with weak cubic anharmonicity as a function of the final phase-space coordinates at a scaled time : Wigner propagator in the semiclassical approximation (23) based on the exact quantum result (37) (panel a), the van Vleck approximation (II.4.2) (b), phase-space path integration (c), and short-time version (28) of the same approach (e). Panels d,f are contributions to Eq. (II.4.2) corresponding to hyperbolic (cosine term, d) and elliptic (sine term, f) trajectory pairs, respectively, cf. Fig. 2b,c. Parameter values are , , , . The underlying classical trajectory (red line in all panels) has been launched at . Color code ranges from red (negative) through white (zero) through blue (positive).

Equation (28) has been obtained combining perturbation theory with semiclassics on the level of the van Vleck propagator augmented by a uniform approximation—a route that might appear less systematic and controlled than the access through phase-space path integration that led to Eq. (23). Its virtue, though, lies in the fact that it is based on trajectory pairs, hence readily extends to a higher number of degrees of freedom, while this is far from obvious in the case of path integration. This feature together with its simplicity and the fact that is is adapted to short propagation time suggests Eq. (28) as a suitable choice for on-the-fly molecular-dynamics applications MTM99 (); DM02a ().

ii.4.3 Asymptotics , ,

As pointed out the discussion of Eq. (II.2), the semiclassical Wigner propagator replaces the Liouville propagator by a smooth quantum spot only if the potential is not purely harmonic. It should converge to the Liouville propagator for vanishing anharmonicity. Indeed, invoking the asymptotic expression VS04 () for the Airy functions in Eq. (23), we find for , . Separating and in the arguments of the delta functions, this proves to be equivalent to , the classical Liouville propagator. Considering the dependence of on , , and , this includes the classical limit (for subtleties of this limit for Wigner functions and their propagation, see CF&08 ()) and the limit . Moreover, the powers , , and , respectively, with which , , and appear in , indicate a hierarchy of the corresponding limits: At finite , the smooth quantum spot already collapses to a delta function in the limit of weak anharmonicity, i.e., for a quantum harmonic oscillator. And even at finite , taking into account the different scaling with time of and as they enter Eq. (II.4.2), the lateral () extension of the quantum spot disappears first with while its longitudinal () dimension (in -direction) scales with as it does with .

Iii Numerical results

iii.1 Models

iii.1.1 Morse oscillator

We choose the Morse oscillator for a detailed study of semiclassical Wigner propagation for two reasons: Above all, being prototypical for strongly anharmonic molecular potentials and correspondingly complex dynamics KK04 (); BB06 (); BL07 (); RG09 (), it is a widely used benchmark for numerical methods in this realm DS88 (); SG96 (); Gro99 (); HRG04 (); KRO04 (). Secondly, the Morse oscillator has the remarkable advantage, shared only with very few other potentials, that closed analytical expressions are available not only for its energy eigenstates Mor29 () but even for their Wigner representations DS88 (), a feature that greatly facilitates the comparison with exact quantum-mechanical results LS82 (); DS88 (); HM&98a (); HM&98bc (); FR&00 (); SK06 () and hence provides a solid basis for an objective test of the performance of semiclassical approximations to the propagator.

The one-dimensional Morse potential Mor29 ()


is determined by the depth and inverse width of the potential well. Quantum-mechanically, its spectrum is discrete for and continuous for . The number of bound states is given (up to ) by the parameter , to be interpreted as an inverse Planck’s constant in natural units. Analytical expressions in terms of Laguerre polynomials are available for the eigenfunctions , see Refs. Mor29 (); DS88 (). Applying the Weyl transform (1) to then leads to the corresponding Wigner eigenfunctions DS88 (). Exact solutions for the continuum states, on the other hand, are not known. Numerical evidence indicates, however, that the error caused by their omission in the basis set underlying quantum calculations of the time evolution is acceptable as long as the energy does not come too close to the threshold .

iii.1.2 Quartic double well

The quartic double-well potential is in many respects complementary to the Morse oscillator. From a phenomenological point of view, it is the standard model for the study of coherent tunneling, and therefore constitutes a particularly hard problem for semiclassical propagation. Technically, it is characterized by the absence of a continuum, which facilitates quantum calculations. At the same time, no analytical expressions for its eigenfunctions are known.

To begin with, define the quartic double-well potential as , where is the oscillation frequency near the minima at and their depth. In natural units and , the Schrödinger equation reads


Its only parameter, the dimensionless barrier height , measures approximately the number of tunneling doublets, half the number of eigenstates below the barrier top.

Eigenvalues and eigenfunctions underlying the exact quantum calculations shown in the sequel have been obtained in a basis of harmonic-oscillator eigenstates, with unit frequency and centered at .

iii.2 Propagating delta functions: the propagator as a stand-alone quantity

To our best knowledge, no detailed account of the propagator of the Wigner function as a quantity of its own right has been published to date, except for the recent Ref. DVS06 (). Indeed, it is inaccessible to semiclassical approximation schemes based on Gaussian smoothing. We therefore find it appropriate to illustrate its basic properties with data obtained for the two models featured in this section. Moreover, the interference pattern characterizing the Wigner propagator is a sensitive diagnostic for the performance of the different semiclassical approximations we are proposing. In addition, an analysis of the propagator allows to assess by mere inspection the validity of propagation schemes based on non-classical but deterministic trajectories LS82 (); Lee92 (); Raz96 (); DM01 (); THW03 (), see Sec. III.4.

The naked propagator is equivalent to the time evolution of a delta function as initial Wigner function, which is not an admissible Hilbert-space element. Notwithstanding, since Wigner functions are measurable, e.g., through quantum-state tomography SB&93 (), so is the propagator, under weak conditions fn1 (), by unfolding the final Wigner function from the initial one, which adds legitimacy to considering the propagator as a stand-alone quantity.

Figure 4: Van-Vleck-based semiclassical Wigner propagator (10, II.2) (b) and corresponding exact quantum result (37) (a) as a function of the final phase-space coordinates for the Morse oscillator (29) at intermediate energy at , approximately the period of the underlying classical orbit (red line). Parameters are , , , , the initial point is . Color code as in Fig. 3.

The linear checkerboard structure of the quantum spot appearing in Fig. 3 is owed to the cubic anharmonicity of the potential to which the figure refers (its poor resolution in the quantum result, panel a, is a numerical artefact). In Fig. 4, we contrast it with data for the Morse oscillator, at an intermediate energy where the nonlinearity is already considerable, and after a propagation time of the order of the period of the underlying classical orbit. The pattern of Fig. 3 can still be discerned but is severely distorted. The nodelines are now strongly curved, and the pattern even folds over, giving rise to a phase-space swallow-tail catastrophe. The van-Vleck-based semiclassical propagator (10, II.2) (panel b) captures this complicated structure surprisingly well. The path-integral-based approximation, by contrast, does not reproduce these distortions (not shown).

Figure 5: As Fig. 4 but for the quartic double well (30) with at energy above the barrier at , approximately the period of the underlying classical orbit (red line). The initial point is . Color code as in Fig. 3.

Figure 5 shows a similar comparison for the quartic double well, at an energy above the barrier top. Here, the interference pattern is even more complex than in Fig. 4. While the van-Vleck-based propagator does not agree with the quantum calculation in the intricate fine details of the oscillatory fringes, it correctly reproduces much of the more global features of the distribution. For an account of its performance at energies below the barrier top, see Sec. III.6.

iii.3 Propagating Gaussians

Gaussian initial states deserve their ubiquity not only for their unique physical properties, exemplified by coherent states. In the context of the Wigner representation, Gaussians gain special relevance as they constitute the only admissible Wigner functions that are positive definite and therefore can be interpreted in probabilistic terms Hud74 (). They have achieved a fundamental rôle for semiclassical propagation as they provide a natural smoothing which allows to reduce the time evolution of an entire phase-space region to the propagation along a single classical trajectory.

By difference to Gaussian-wavepacket propagation and semiclassical IVRs, however, the semiclassical Wigner propagator does not involve any smoothing by construction and thus allows to time-evolve arbitrary initial states, Gaussian or not. A more specific question is then how semiclassical approximations to it operate on the space of Gaussian phase-space distributions. We find that they generally transform them into a broader class of functions which, in view of the above theorem, can no longer remain positive definite.

Define Gaussians in phase space Hel91a () by


The covariance matrix controls size, shape, and orientation of the Gaussian centered in . The more specific class of minimum-uncertainty Gaussians, equivalent to Wigner representations of coherent states Sch01 (), is characterized by . In what follows, in two dimensions , we choose , so that


equivalent to a position-space wave function Hel91a () .

iii.3.1 Autocorrelation functions

A standard interface between dynamical and spectral data, autocorrelation functions have a wide range of applications in atomic and molecular physics and provide a robust and easily verifiable assay of the accuracy and efficiency of propagation methods. The overlap for an initial state upon squaring readily translates into an expression in terms of the Wigner representation of that state, , or, involving the propagator explicitly,