Reconciling Semiclassical and Bohmian Mechanics: II. Scattering states for discontinuous potentials
aa r X i v : . [ qu a n t - ph ] F e b Reconciling Semiclassical and Bohmian Mechanics:II. Scattering states for discontinuous potentials
Corey Trahan and Bill Poirier
Department of Chemistry and Biochemistry, and Department of Physics,Texas Tech University, Box 41061, Lubbock, Texas 79409-1061 ∗ In a previous paper [J. Chem. Phys. +Ψ was presented for stationary bound states Ψ of the one-dimensional Schr¨odinger equation, suchthat the components Ψ and Ψ approach their semiclassical WKB analogs in the large action limit.Moreover, by applying the Madelung-Bohm ansatz to the components rather than to Ψ itself, theresultant bipolar Bohmian mechanical formulation satisfies the correspondence principle. As a result,the bipolar quantum trajectories are classical-like and well-behaved, even when Ψ has many nodes, oris wildly oscillatory. In this paper, the previous decomposition scheme is modified in order to achievethe same desirable properties for stationary scattering states. Discontinuous potential systems areconsidered (hard wall, step, square barrier/well), for which the bipolar quantum potential is found tobe zero everywhere, except at the discontinuities. This approach leads to an exact numerical methodfor computing stationary scattering states of any desired boundary conditions, and reflection andtransmission probabilities. The continuous potential case will be considered in a future publication. I. INTRODUCTION
Much attention has been directed by theoreti-cal/computational chemists towards developing reliableand accurate means for solving dynamical quantummechanics problems—i.e., for obtaining solutions tothe time-dependent Schr¨odinger equation—for molecu-lar systems. Insofar as “exact” quantum methods areconcerned, two traditional approaches have been used:(1) representation of the system Hamiltonian in a fi-nite, direct-product basis set; (2) discretization of thewavefunction onto a rectilinear grid of lattice points overthe relevant region of configuration space. Both ap-proaches, however, suffer from the drawback that thecomputational effort scales exponentially with systemdimensionality.
Recently, a number of promising newmethods have emerged with the potential to alleviatethe exponential scaling problem once and for all. Theseinclude various basis set optimization methods, and build-and-prune methods, such as those based onwavelet techniques. On the other hand, a completely different approachto the exponential scaling problem is to use basis sets orgrid points, that themselves evolve over time. The idea isthat at any given point in time, one need sample a muchsmaller Hilbert subspace, or configuration space region,than would be required at all times—thus substantiallyreducing the size of the calculation. For basis set calcula-tions, much progress along these lines has been achievedby the multi-configurational time-dependent Hartree(MCTDH) method, developed by Meyer, Manthe and co-workers.
More recently, time-evolving grid, or “quan-tum trajectory” methods (QTMs) havealso been developed, and for certain types of systems,successfully applied at quite high dimensionalities.
QTMs are based on the hydrodynamical picture ofquantum mechanics, developed over half a century agoby Bohm and Takabayasi, who built on the earlier work of Madelung and van Vleck. QTMs are inher-ently appealing for a number of reasons. First, they offeran intuitive, classical-like understanding of the underly-ing dynamics, which is difficult-to-impossible to extractfrom more traditional fixed grid/basis methods. In ef-fect, quantum trajectories are like ordinary classical tra-jectories, except that they evolve under a modified po-tential V + Q , where Q is the wavefunction-dependent“quantum potential” correction. Second, QTMs holdthe promise of delivering exact quantum mechanical re-sults without exponential scaling in computational effort.Third, they provide a pedagogical understanding of en-tirely quantum mechanical effects such as tunneling and interference. They have already been used tosolve a variety of different types of problems, includingbarrier transmission, non-adiabatic dynamics, andmode relaxation. Several intriguing phase space gener-alizations have also emerged, of particular rele-vance for dissipative systems.
Despite this success, QTMs suffer from a significantnumerical drawback, which, to date, precludes a com-pletely robust application of these methods. Namely:QTMs are numerically unstable in the vicinity of am-plitude nodes. This “node problem” manifests in severaldifferent ways: (1) infinite forces, giving rise to kinky,erratic trajectories; (2) compression/inflation of trajecto-ries near wavefunction local extrema/nodes, leading to;(3) insufficient sampling for accurate derivative evalua-tions. Nodes are usefully divided into two categories, depending on whether Q is formally well-behaved (“typeone” nodes) or singular (“type two” nodes). For sta-tionary state solutions to the Schr¨odinger equation, forinstance, all nodes are type one nodes. In principle, typeone nodes are “gentler” than type two nodes; however,from a numerical standpoint, even type one nodes willgive rise to the problems listed above, because the slight-est numerical error in the evaluation of Q is sufficient tocause instability.In the best case, the node problem simply results insubstantially more trajectories and time steps than thecorresponding classical calculation; in the worst case, theQTM calculation may fail altogether, beyond a certainpoint in time. Several numerical methods, both “exact”and approximate, are currently being developed to dealwith this important problem. The latter category in-cludes the artificial viscosity and linearized quantumforce methods, both of which have proven to be verystable. While such approximate methods may not cap-ture the hydrodynamic fields with complete accuracy innodal regions, they do allow for continued evolution andlong-time solutions, often unattainable via use of a tradi-tional QTM. The “exact” methods include the adaptivehybrid methods, and the complex amplitude method. In the adaptive hybrid methods, for which hydrodynamictrajectories are evolved everywhere except for in nodalregions, where the time-dependent Schrodinger equationis solved instead to avoid node problems. Although theyhave been applied successfully for some problems, thesemethods are difficult to implement numerically, since notonly must the hydrodynamic fields be somehow moni-tored for forming singularities, but there must also bean accurate means for interfacing and coupling the twocompletely different equations of motion. The complexamplitude method is cleaner to implement, but is onlyexact for linear and quadratic Hamiltonians.In a recent paper, hereinafter referred to as “pa-per I,” one of the authors (Poirier) introduced a newstrategy for dealing with the node problem, based on abipolar decomposition of the wavefunction. The idea isto partition the wavefunction into two (or in principle,more) component functions, i.e. Ψ = Ψ + Ψ . Onethen applies QTM propagation separately to Ψ and Ψ ,which can be linearly superposed to generate Ψ itself atany desired later time. In essence, this works because theSchr¨odinger equation itself is linear, but the equivalentBohmian mechanical, or quantum Hamilton’s equationsof motion (QHEM) are not. In principle, therefore, onemay improve the numerical performance of QTM cal-culations simply by judiciously dividing up the initialwavepacket into pieces.Although bipolar decompositions have been aroundfor quite some time, their use as a tool for cir-cumventing the node problem for QTM calculationsis quite recent. Two promising new exact methodsthat seek to accomplish this are the so-called “counter-propagating wave” method (CPWM), and the “cov-ering function” method (CFM). In the CPWM, thebipolar decomposition is chosen to correspond to thesemiclassical WKB approximation, for which all of thehydrodynamic field functions are smooth and classical-like, and the component wavefunctions are node-free.Interference is achieved naturally, via the superpositionof left- and right-traveling (i.e. positive- and negative-momentum) waves. For one-dimensional (1D) stationarybound states, it can be shown that the resultant bipo-lar quantum potential q ( x ) becomes arbitrarily small inthe large action limit, even though the number of nodes becomes arbitrarily large. (Note: in accord with the con-vention established in Ref. 24, upper/lower case will beused to denote the unipolar/bipolar field quantities). Inthe CFM, the idea is to superpose some well-behavedlarge-amplitude wave, with the actual ill-behaved (nodalor wildly oscillatory) wave, so as to “dilute” the undesir-able numerical ramifications of the latter.This paper is the second in a series designed to ex-plore the CPWM approach, introduced in paper I. Asdiscussed there in greater detail, there are many motiva-tions for this approach, but the primary one is to recon-cile the semiclassical and Bohmian theories, in a mannerthat preserves the best features of both, and also sat-isfies the correspondence principle. For our purposes,this means that the Lagrangian manifolds (LMs) for thetwo theories should become identical in the large actionlimit (Sec. II B). As described above, a key benefit ofthe CPWM decomposition is an elegant treatment of in-terference, the chief source of nodes and “quasi-nodes” (i.e. rapid oscillations) in quantum mechanical systems.An interesting perspective on the role of interference insemiclassical and Bohmian contexts is to be found in arecent article by Zhao and Makri. Whereas paper I focused on stationary bound states for1D systems, the present paper (paper II) and the nextin the series (paper III) concern themselves with sta-tionary scattering states. The CPWM decomposition ofpaper I is uniquely specified for any arbitrary 1D state—bound or scattering—and in the bound case, always sat-isfies the correspondence principle. However, the non- L nature of the scattering states is such that the paper Idecomposition generally does not satisfy the correspon-dence principle in this case. Simply put, the quantumtrajectories and LMs exhibit oscillatory behavior in atleast one asymptotic region (thereby manifesting reflec-tion), whereas the semiclassical LMs do not. This is nota limitation of the CPWM, but is rather due to the fun-damental failure of the basic WKB approximation to pre-dict any reflection whatsoever for above-barrier energies,as has been previously well established. In semi-classical theory, a modification must therefore be made tothe basic WKB approximation, in order to obtain mean-ingful scattering quantities. As discussed in Sec. II B andin paper III, our approach will be to apply a similar mod-ification to the exact quantum decomposition (actually, a reverse modification) such that the correspondence prin-ciple remains satisfied, and the two theories thus recon-ciled, even for scattering systems.It will be shown the modified CPWM gives rise tobipolar Bohmian LMs that are identical to the semi-classical LMs, regardless of whether or not the actionis large. Put another way, this means that the bipolarquantum potentials q effectively vanish , so that the resul-tant quantum trajectory evolution is completely classical .Moreover, the resultant component wavefunctions, Ψ ( x )and Ψ ( x ), correspond asymptotically to the familiar “in-cident,” “transmitted,” and “reflected” waves of tradi-tional scattering theory. Thus, the modified CPWM im-plementation of the bipolar Bohmian approach providesa natural generalization of these conceptually fundamen-tal entities throughout all of configuration space , not justin the asymptotic regions, as is the case in conventionalquantum scattering theory.The above conclusions will be demonstrated for bothdiscontinuous and continuous potential systems, in pa-pers II and III, respectively. Discontinuous potentials—e.g. the hard wall, the step potential, and the squarebarrier/well—serve as a useful benchmark for the modi-fied CPWM approach, because the scattering componentwaves (e.g. “incident wave,” etc.) in this case are well-defined throughout all of configuration space, accordingto a conventional scattering treatment. Although this isno longer true for continuous potentials, the foundationlaid here in paper II can be extended to the continu-ous (and also time-dependent) case as well, as describedin paper III. Additional motivation for the developmentof a scattering version of the CPWM, vis-a-vis the rel-evance for chemical physics applications, is provided inpaper III. Additional motivation for the consideration ofdiscontinuous potentials is provided in Sec. II B of thepresent paper. II. THEORYA. Background
1. Bohmian mechanics
According to the Bohmian formulation, theQHEM are derived via substitution of the 1D (unipolar)wavefunction ansatz,Ψ( x, t ) = R ( x, t ) e iS ( x,t ) / ¯ h (1)into the time-dependent Schr¨odinger equation. For the1D Hamiltonian, ˆ H = − ¯ h m ∂ ∂x + V ( x ) , (2)this results in the coupled pair of nonlinear partial dif-ferential equations, ∂S ( x, t ) ∂t = − S ′ m − V ( x ) + ¯ h m R ′′ R ,∂R ( x, t ) ∂t = − m R ′ S ′ − m R S ′′ , (3)where m is the mass, V ( x ) is the system potential, andprimes denote spatial partial differentiation.The first of the two equations above is the quantumHamilton-Jacobi equation (QHJE), whose last term isequal to − Q ( x, t ), i.e. comprises the quantum potentialcorrection. The second equation is a continuity equation.When combined with the quantum trajectory evolution equations, i.e. P = m dxdt = S ′ ,dPdt = − ( V ′ + Q ′ ) , (4)the continuity equation ensures that the probability [i.e.density, R ( x, t ) , times volume element] carried by indi-vidual quantum trajectories is conserved over the courseof their time evolution.
2. CPWM decomposition for stationary states
In paper I, we derived a unique bipolar decomposition,Ψ( x ) = Ψ + ( x ) + Ψ − ( x ) , (5)for stationary eigenstates Ψ( x ) of 1D Hamiltonians of theEq. (2) form, such that:1. Ψ ± ( x ) are themselves (non- L ) solutions to theSchr¨odinger equation, with the same eigenvalue, E ,as Ψ( x ) itself.2. The invariant flux values, ± F , of the two solutions,Ψ ± ( x ), equal those of the two semiclassical (WKB)solutions.3. The median of the enclosed action, x , equals thatof the semiclassical solutions.There are other important properties of the Ψ ± ( x ), asdiscussed in Sec. I, and in Ref. 24. Nevertheless, theabove three conditions are sufficient to uniquely specifythe decomposition. In the special case of bound (i.e. L )stationary states, the real-valuedness of Ψ( x ) implies thatthe Ψ ± ( x ) are complex conjugates of each other. B. Scattering systems
It is natural to ask to what extent the above analysismay be generalized for scattering potentials. Certainly,Ψ( x ) itself is no longer L , nor even real-valued, andthere are generally two linearly independent solutionsof interest for each E , instead of just one. Condition(1) above poses no difficulty for Ψ ± ( x ), as these com-ponent functions are non- L and complex-valued, evenin the bound eigenstate case. In principle, condition (2)is not difficult either; although the flux value dependson the normalization of Ψ itself, which is not L , cer-tain well-established normalization conventions for scat-tering states exist, that can be applied equally well tosemiclassical and exact quantum solutions. There is noaction median per se for scattering states, as the actionenclosed within the Ψ + ( x ) and Ψ − ( x ) phase space La-grangian manifolds (LMs) is infinite; however,the scattering analog of condition (3) is related to theasymptotic boundary conditions, and it is here that oneencounters difficulty. Moreover, an additional concern israised by the doubly-degenerate nature of the continuumeigenstates, namely: should each scattering Ψ( x ) haveits own Ψ ± ( x ) decomposition, or should there be a sin-gle Ψ ± ( x ) pair, from which all degenerate Ψ( x )’s may beconstructed via arbitrary linear superposition?To resolve these issues, we will adopt the same gen-eral strategy used in paper I, i.e. we will resort to semi-classical theory as our guide, wherever possible. Wewill also exploit certain special features of the scatter-ing problem not found in generic bound state systems,such as the asymptotic potential condition V ′ ( x ) → x → ±∞ (where primes denote spatial differentia-tion), and its usual implications for scattering theory andapplications. The basic WKB solutions are given byΨ sc ± ( x ) = r sc ( x ) e ± is sc ( x ) / ¯ h , (6)where r sc ( x ) = s mFs ′ sc ( x ) and s ′ sc ( x ) = p m [ E − V ( x )](7)The corresponding positive and negative momentumfunctions, specifying the semiclassical LMs, are given by p sc ± ( x ) = ± s ′ sc ( x ). Equations (6) and (7) apply to bothbound and scattering cases; note that for both, Ψ sc ± ( x ) arecomplex conjugates of each other. The asymptotic poten-tial condition ensures that these approach exact quantumplane waves asymptotically, with the usual scattering in-terpretations, i.e. Ψ + ( x ) in the x → −∞ asymptoticregion is the incoming wave from the left (usually takento be the incident wave), Ψ + ( x ) as x → ∞ is the out-going wave from the left (the usual transmitted wave),etc.Insofar as determining the corresponding exact quan-tum solutions Ψ ± ( x ), the procedure described in paperI is still appropriate for bound and semi-bound (i.e. onone side only) states, in that the results satisfy the corre-spondence principle globally, as desired (for semi-boundexamples, consult the Appendix). For true scatteringstates, however, this procedure fails, in the sense that ifΨ + ( x ) is chosen to match the normalization and flux ofΨ sc+ ( x ) in the x → ∞ asymptote, then it will necessar-ily approach a nontrivial linear superposition of Ψ sc+ ( x )and Ψ sc − ( x ) in the x → −∞ asymptote, and vice-versa.There is therefore an ambiguity as to how the corre-sponding quantum Ψ ± ( x )’s should be defined, i.e. whichasymptotic region should be used to effect the correspon-dence. More significantly though, either choice will resultin component functions Ψ ± ( x ) with substantial interfer-ence in one of the two asymptotic regions. This is due topartial reflection of the exact quantum scattering states,which is not predicted by the basic WKB approxima-tion. Thus, in the large action limit, the exact quan-tum solutions manifest large-magnitude quantum poten-tials, q ± ( x ), and rapidly oscillating field functions q ± ( x ), r ± ( x ), and p ± ( x )—exactly the undesirable behavior thatthe CPWM was introduced to avoid—whereas the corre-sponding basic WKB functions are smooth, and asymp-totically uniform.The lack of any partial reflection is a well-understoodshortcoming of the WKB approximation —i.e.,the basic Ψ sc ± ( x ) components, though elegantly con-structed from smooth classical functions r sc ( x ) and s sc ( x ), do not in and of themselves correspond to anyactual quantum scattering solutions Ψ( x ). In light of thebipolar decomposition ideas introduced in paper I, how-ever, our perspective is the reverse one: for any actual quantum Ψ( x ), can one determine an Eq. (5) decompo-sition such that the resultant Ψ ± ( x ) resemble their well-behaved semiclassical counterparts, and is such a decom-position unique? Among other properties, the Ψ ± ( x )LM’s should become identical to the semiclassical LM’s inthe large action limit, so as to satisfy the correspondenceprinciple. Based on the considerations of the previousparagraph it is clear that the paper I decomposition doesnot achieve this goal, when applied to stationary scatter-ing states.We defer a full accounting of these issues—in thecontext of completely arbitrary continuous potentials V ( x )—to paper III, wherein it will be demonstrated howto compute exact quantum reflection and transmissionprobabilities (and stationary scattering states) using onlyclassical trajectories, and without the need for explicitnumerical differentiation of the wavefunction. In thepresent paper, we lay the foundation for paper III, byfocusing attention onto two key aspects whose develop-ment comprises an essential prerequisite.First, as the paper III approach treats V ( x ) as a se-quence of steps, the present paper II will focus exclu-sively on the step potential and related discontinuouspotential systems, for which V ( x ) = const in betweensuccessive steps. Discontinuous potentials are importantfor chemical physics, because they model steep repulsivewells, and are used in statistical theories of liquids. More-over, they hold a special significance for QTM methods,for which they serve as a “worst-case scenario” bench-mark. Indeed, conventional QTM techniques always failwhen applied to discontinuous potentials. To date, Theonly such calculations that have been performed havecomputed the quantum potential from a completely sep-arate time-dependent fixed-grid calculation (the “analyt-ical approach”) rather than directly from the quan-tum trajectories themselves. Even if one could propagatetrajectories for discontinous systems using a traditionalQTM, the trajectories that would be generated would bevery kinky and erratic, and a great many time trajec-tories and time steps would thus be required.Second, since the new Ψ ± ( x ) do not satisfy condition(1), unlike the paper I CPWM decomposition, the timeevolution of these two component functions is clearly notthat of the time-dependent Schr¨odinger equation. More-over, since the | Ψ ± | are constant over time [because | Ψ | itself is stationary, and Eq. (5) is presumed unique], the two Ψ ± ( x, t ) time evolutions must be coupled together .It is essential that the nature of this coupling be com-pletely understood, in order that the present approachmay be generalized to non-stationary state situations—e.g. to wavepacket scattering, as will be discussed infuture publications. The ramifications for QTMs areequally important. Accordingly, the present paper fo-cuses on the QTM propagation of the wavefunction andits bipolar components—with a keen eye towards gener-ality and physical interpretation—even though the statesinvolved are stationary. This approach leads to a peda-gogically useful reinterpretation of “incident,” “transmit-ted,” and “reflected” waves—very reminiscent of ray op-tics in electromagnetic theory—which is applicable muchmore generally than traditional usage might suggest. C. Basic applications
The necessary theory will be developed over the courseof a consideration of various model application systemsof increasing complexity.
1. free particle system
Let us first consider the simplest case imaginable, thefree particle system, V ( x ) = 0. In this case, the exact so-lutions Ψ ± ( x ) = Ψ sc ± ( x ) clearly satisfy the conditions ofSec. II A 2, and the bipolar quantum potentials q ± ( x ) arezero everywhere. Thus, the bipolar decomposition devel-oped for bound states in paper I can be used directly withthis continuum system, requiring only the slight modifi-cation that arbitrary linear combinations of Ψ sc+ ( x ) andΨ sc − ( x ) are to be allowed, in order to construct arbitraryscattering solutions Ψ( x ). For convenience, the linearcombination coefficients will from here on out be directlyincorporated into the amplitude functions, r ± ( x ), andphase functions, s ± ( x ), so that Eq. (5) is still correct.If from all solutions Ψ( x ) one considers only that whichsatisfies the usual scattering boundary conditions (i.e. in-cident wave incoming from the left) then the negative mo-mentum wave Ψ − vanishes, and Ψ( x ) = Ψ + ( x ). Thereis zero reflection, and 100% transmission. Put anotherway, the incident flux, lim x →−∞ j + ( x ), is equal to thetransmitted flux, lim x → + ∞ j + ( x ), where j ± ( x ) = ¯ h im (cid:20) Ψ ∗± ( x ) d Ψ ± ( x ) dx − d Ψ ∗± ( x ) dx Ψ ± ( x ) (cid:21) = (cid:20) p ± ( x ) m (cid:21) r ± ( x ) , (8)[both flux values are equal to F , as in Eq. (7)].In the quantum trajectory description, flux mani-fests as probability-transporting trajectories, which movealong the LMs. For the boundary conditions describedabove, there are only positive momentum trajectories,moving uniformly from left to right with momentum p + ( x ) = √ mE . If a Ψ − ( x ) contribution were present,its trajectories would move uniformly in the opposite di-rection [ p − ( x ) = −√ mE .] Since the two componentsΨ ± ( x ) are in this case uncoupled, the positive and neg-ative momentum trajectories would have no interactionwith each other.
2. hard wall system
We next consider the hard wall system: V ( x ) = n x ≤ ∞ for x >
0. (9)In the x ≤ ± ( x ) components are ex-actly the same as in the free particle case, except thatthe Ψ(0) = 0 boundary condition imposes the additionalconstraints, s − (0) = s + (0) + π mod(2 π ) ; r − (0) = r + (0) . (10)This also results in only one linearly independent solutioninstead of two, i.e. Ψ( x ) ∝ sin( kx ), with k = √ mE/ ¯ h .Regarding the LMs and trajectories, in the x < + ( x )LM trajectories move uniformly to the right, towards thehard wall at x = 0.It is natural to ask what happens when the Ψ + ( x ) LMtrajectories actually reach x = 0. There are two rea-sonable interpretations. The first is that the trajectorieskeep moving uniformly into the x > x > This underscores the fact thatunlike Ψ( x ) itself, the individual Ψ ± ( x ) components perse are unconstrained at the origin—though the Eq. (10)constraint implies a unique correspondence between thetwo. This interpretation also makes it clear that for thehard wall system, the paper I decomposition is essentiallyidentical to the present decomposition, as is worked outin detail in the Appendix.In the second interpretation, the effect of the hardwall at x = 0 is to cause instantaneous elastic reflec-tion of a Ψ + ( x ) LM trajectory momentum, from p = p + = + √ mE to p = p − = −√ mE . Afterwards,the reflected trajectory propagates uniformly backward,along the Ψ − ( x ) LM. In this interpretation, the trajec-tories never leave the allowed configuration space, x ≤ hopping from one LM to the other—not un-like that previously considered, e.g., in the context ofnon-adiabatic transitions. The trajectory hopping in-terpretation is adopted in the present paper, and in pa-per III, but the first interpretation will also be recon-sidered in later publications. Note that for discontin-uous potentials—and indeed more generally —one canregard trajectory hopping as the source of Ψ ± ( x ) inter-action coupling.For the hard wall case, trajectory hopping only mani-fests at x = 0, the sink of all Ψ + ( x ) LM trajectories, andthe source of all Ψ − ( x ) LM trajectories. If these trajec-tories are to be regarded as one and the same via hop-ping, then a unique field transformation for r , s , and allspatial derivatives, must be specified. Fortunately, theunique correspondence between Ψ + ( x ) and Ψ − ( x ) de-scribed above, enables one to do just that. In particular,Eq. (10) specifies the correct transformations for r and s ,as transported by the quantum trajectories. All spatialderivatives of arbitrary orders can then be obtained viaspatial differentiation of Eq. (6)—although in the hardwall case, only the s ′ condition, p − (0) = − p + (0) is rele-vant, because all higher order derivatives are identicallyzero.Since the magnitudes of the p and r fields associatedwith a given quantum trajectory are unchanged as a re-sult of the trajectory hop, Eq. (8) implies that the inci-dent and reflected flux values are the same (apart fromsign), and so the scattering system exhibits 100% reflec-tion and zero transmission (along each LM, the flux isinvariant ). These basic facts of the hard wall system areof course well understood. The point, though, is that wehave now obtained the information in a time- dependent quantum trajectory manner, rather than through theusual route of applying boundary conditions to time- independent piecewise component functions. In otherwords, Eq. (10) now refers to individual quantum tra-jectories , rather than to wavefunctions.This shift of emphasis is very important, and leadsto quite a number of conceptual and computational ad-vantages. For instance, the standard description of thehard wall stationary states would decompose these intoplane wave components interpreted as “incident” and“reflected” waves. This language suggests a process, orchange over time—i.e. a state that is initially incident,at some later time is somehow transformed into a re-flected state. Nothing in the standard description, how-ever, would seem to render transparent the usage of suchterminology, i.e. Ψ( x ) is stationary, and so the reflectedand transmitted components are in fact both present forall times. Of course, a localized superposition of station-ary states, i.e. a wavepacket, may well exhibit such anexplicit transformation over the course of the time evo-lution, as such a state is decidedly non-stationary. In-deed, wavepackets are relied upon by the more rigorousformulations of scattering theory, in order to justify theuse of terms such as “reflected wave,” even in a station-ary context. Such formulations, though certainly legit-imate, seem always to require a clever use of limits, thesubtle distinction between unitary and isometric trans-formations, and other esoteric mathematical tricks.On the other hand, the time-dependent bipolar quan-tum trajectory hopping picture presented above providesa physicality to such language that is immediately ap-parent. Over the course of the time evolution, althoughthe wavefunction as a whole is stationary, each individ-ual trajectory is first incident from the left, then collides with the hard wall, and is subsequently reflected backtowards the left (i.e. towards x → −∞ ). The bipo-lar quantum trajectories are all classical, as the bipolarquantum potentials, q ± ( x ), are zero everywhere exceptat the wall itself. Interference arises naturally from thesuperposition of the two LMs—i.e., from the trajectoriesthat have already progressed to the point of reflecting,vs. those that have not reflected yet. In contrast, sinceΨ( x ) itself exhibits very substantial interference, and aninfinite number of nodes, the traditional unipolar QTMtreatment would be very ill-behaved, i.e. R ( x ) wouldoscillate wildly in the large k limit, and Q ( x ) would benumerically unstable near the nodes. Apart from theseimportant pragmatic drawbacks, the incident/reflectedinterpretation of the quantum trajectories would also belost.The bipolar quantum trajectory description of the hardwall system is very reminiscent of ray optics, as used todescribe the reflection of electromagnetic waves off of aperfectly reflecting surface. Indeed, much can be gainedfrom applying a ray optics analogy to quantum scatteringapplications, especially where discontinuous potentialsare concerned. One can construct a simple gedanken-experiment as follows. Let x L < t = 0, all trajectories on thepositive LM lying to the right of x L are ignored , as is thenegative LM altogether. One then evolves the retainedtrajectories over time, and monitors the contribution thatjust these trajectories make to the total wavefunction. Insome respects, it is as if the point x L were serving as theinitial wavefront for some incoming wave, that at t = 0had not yet reached the hard wall/reflecting surface. Ofcourse, if the actual wave were in fact truncated in thisfashion, then the discontinuity in the field functions atthe wavefront would result in a very non-trivial propaga-tion over time, owing to the high-frequency componentsimplicitly present. For actual waves, the precise natureof the wavefront is known to have a tremendous impacton the resultant dynamics. We avoid such complicat-ing details by always interpreting the “actual wave” tobe the full stationary wave itself, i.e. the truncation isconceptual only.In the ray optics analogy, the above situation is like asource of light located at x L , which is suddenly “turnedon” at t = 0. It takes time for the wavefront to propagateto the reflecting surface, and additional time for the re-flected wavefront to make its way back to x = x L . Priorto the latter point in time, the evolution of the trun-cated electromagnetic wave is decidedly non -stationary;afterwards however, a stationary wave is achieved, atleast within the region of interest, x L ≤ x ≤
0, as thewavefront has by this stage propagated beyond this re-gion. The same qualitative comments apply to the bipo-lar quantum case, although of course the evolution equa-tions are different.A similar prescription may be used to achieve rudi-mentary “wavepacket dynamics,” even in the context ofpurely stationary states. Instead of retaining all initialtrajectories that lie to the left of x L , one retains onlythose that lie within some finite interval. The result-ing time evolution is analogous to a light source that isturned on at t = 0, and then turned off at some later time(prior to when the wavefront arrives at the reflecting sur-face). The initial “wavepacket” has uniform density, andmoves with uniform speed towards the hard wall. Inter-ference fringes then form after the foremost trajectorieshave been reflected onto the negative LM. Eventually,all trajectories within the interval are reflected, at whichpoint interference ceases (the nodes are “healed” ), uni-form density is restored, and the reflected wave travelswith uniform speed in the reverse direction, back towardsthe starting point x L . Qualitatively, this behavior isclearly similar to that undergone by actual wavepacketsreflecting off of barrier potentials. D. More complicated applications
The ideas described above can be easily extended tomore complicated discontinuous potential systems, suchas up- and down-step potentials, and any combination ofmultiple steps, e.g. square barriers and square wells. Inpaper III, they will even be extended to arbitrary contin-uous potentials. In every case, the ray optics analogyfrom electromagnetic theory may also be extended ac-cordingly. This approach provides a useful perspective onglobal reflection and transmission in scattering systems,and in particular, demonstrates how such quantities maybe obtained from a single, universal expression for localreflection and transmission.
1. step potential system—above barrier energies
We next consider the step potential system: V ( x ) = (cid:26) x ≤ V for x >
0, (11)Classically, this system exhibits 100% transmission if thetrajectory energy is above the barrier (i.e.
E > V ), and100% reflection if the trajectory energy is below the bar-rier ( E < V ). Quantum mechanically, all above barriertrajectories are found to exhibit partial reflection andpartial transmission, although there is a general increasein transmission probability with increasing energy. Thebelow barrier quantum trajectories exhibit 100% reflec-tion, as in the classical case; however, they also mani-fest tunneling into the classically forbidden x > x → ±∞ . Incoming trajectories can therefore originate from either asymptote, thus giving rise to two linearlyindependent solutions, Ψ( x ). This is in stark contrastto the hard wall system, for which incoming trajectoriescould only originate from x → −∞ , thus resulting in onlyone linearly independent solution for Ψ( x ).In the standard time-independent picture, one startswith the four piecewise solutions,Ψ A ± ( x ) = e ± ip A x/ ¯ h and Ψ B ± ( x ) = e ± ip B x/ ¯ h , (12)where region A corresponds to x ≤ B to x ≥ p A = √ mE and p B = p m ( E − V ) . (13)Matching Ψ( x ) and Ψ ′ ( x ) boundary conditions at x = 0,and specifying asymptotic boundary conditions for Ψ( x ),then enables a unique determination of the four complexcoefficients A ± and B ± inΨ( x ) = (cid:26) A + Ψ A + ( x ) + A − Ψ A − ( x ) for x ≤ B + Ψ B + ( x ) + B − Ψ B − ( x ) for x ≥ . (14)In general, the solution coefficients depend on the par-ticular stationary solution of interest. For the usual scat-tering convention of an incident wave incoming from theleft (Fig. 1) the solutions are A + = 1 ; A − = R = (cid:18) p A − p B p A + p B (cid:19) B + = T = (cid:18) p A p A + p B (cid:19) ; B − = 0 , (15)where R and T are respectively, reflection and transmis-sion amplitudes. When flux is properly accounted for,the resultant reflection and transmission probababilities(which add up to unity) are given by P refl = | R | ; P trans = (cid:18) p B p A (cid:19) | T | . (16)Note that Eqs. (15) and (16) above are correct for bothan “up-step” and a “down-step”—i.e. for V positive ornegative. We can also apply these equations to the “op-posite” boundary conditions, i.e. to an incident waveincoming from the right, by simply transposing A and B subscripts, and + and − subscripts ( p A and p B are stillpositive). This is important, because any stationary so-lution Ψ( x ) can be obtained as some linear superpositionof left-incident and right-incident solutions.Regarding the time-dependent interpretation, it is ev-ident that upon reaching the step discontinuity, left-incident trajectories must be partially reflected and par-tially transmitted. The trajectory is suddenly split intotwo, one that continues to propagate along the positiveLM for the transmitted B region (i.e. the B + LM) andthe other being instantaneously reflected down to the A − LM. Moreover, since probability carried by individ-ual quantum trajectories is conserved, this splitting
FIG. 1: Component waves for a left-incident stationary eigen-state of the up-step barrier problem with
E > V , as describedin Sec. IV A. Solid and dashed lines represent real and imagi-nary contributions, respectively. The dot-dashed line denotesthe location of the step. must be done in a manner that preserves both probabilityand flux. In other words, the local splitting of the trajec-tory at x = 0 must correspond to Eq. (15), which is nowregarded as a local condition, giving rise to local reflectionand transmission amplitudes, R and T . For the presentstep potential case, these local quantities are directly re-lated to the global P refl and P trans values via Eq. (16).For multiple step potentials (Sec. II D 3), the global ex-pressions above [Eq. (16)] no longer apply; however, alocal, time-dependent trajectory version of Eq. (15) does turn out to be correct.Such an expression, immediately applicable to all sin-gle and multiple step systems, can be written as follows: r refl = (cid:18) p i/r − p trans p i/r + p trans (cid:19) r inc ; s refl = s inc (17) r trans = (cid:18) p i/r p i/r + p trans (cid:19) r inc ; s trans = s inc . (18)In the above equations, “inc” refers to any trajectory,locally incident on some particular step from some par-ticular direction, which spawns both a locally reflectedtrajectory, “refl,” and a locally transmitted trajectory,“trans”. The quantity p i/r is the (positive) momentumassociated with the locally incident/reflected trajectory;similarly, p trans (also positive) is associated with the lo-cally transmitted trajectory. For above-barrier incidenttrajectories, note that the local reflection and transmis-sion amplitudes are both real, thus ensuring the realityof r and s for the spawned trajectories.Returning to the step potential system, the ray opticspicture can once again shed some interesting light. Theoptical analog of the step is an interface between two media with different indices of refraction. Light incidenton such an interface will partially reflect back towardsthe original source, and partially refract forwards intothe new medium. The refraction is completely analogousto the discontinuous change in momentum, ( p B − p A ),that suddenly occurs as one crosses the step (Fig. 3). Inany event, the ray optics gedankenexperiment describedin Sec. II C 2 can also be applied to the step potentialsystem, in order to obtain a particular stationary solutionΨ( x ) with any desired boundary conditions.For instance, suppose one is interesting in constructingthe left-incident wave solution, i.e. that of Eq. (15). At t = 0, only the Ψ A + wave is considered, and only thosetrajectories for which x ≤ x L , as before. As the incidenttrajectories reach the step, two new waves are dynami-cally created from the spawned trajectories: a transmit-ted wave traveling to the right, and a reflected wave trav-eling to the left. A plot of the overall density | Ψ( x, t ) | soobtained will change over time, as the transmitted and re-flected wavefronts propagate into their respective regions(Fig. 6). Eventually, however, these wavefronts will prop-agate beyond the region of interest, i.e. x L ≤ x ≤ x R ,where x R > x ) obtainedwithin the region of interest will be exactly equal to thestationary solution with the desired boundary condition.As in the hard wall case, one can also perform steppotential “wavepacket dynamics” by restricting consid-eration to just those initial Ψ A + trajectories lying withinsome coordinate interval. The wavepacket will propagatetowards the step with uniform density and speed. As thefirst few trajectories hit the step, a uniform transmittedwave will be formed in the B region. In the A region,the sudden appearance of a Ψ A − wave will introduce in-terference wiggles in the overall density plot (althoughno nodes per se , owing to partial reflection only). Even-tually, after all trajectories have progressed beyond thestep, well-separated transmitted and reflected wavepack-ets emerge, propagating in their respective spaces and di-rections. There is no longer any interference in the A re-gion, as the incident wave is now gone, having been com-pletely divided into the two final contributions. P refl and P trans values may be determined via monitors placed at x L and x R , either by integrating probability over time asthe respective wavepackets travel through, or by record-ing the (constant) amplitude values R and T , and apply-ing Eq. (16).
2. step potential system—below barrier energies
The case for which the incident trajectory energies arebelow V requires special discussion. In this case, theclassical LMs and trajectories are confined to the A re-gion only (i.e. to x ≤ B region is classi-cally forbidden. In the language of Sec. II B, these belowbarrier states are therefore semi-bound, implying thatthere is only one linearly independent stationary solu- FIG. 2: Wavefunction plot for a left-incident stationary eigen-state of the up-step barrier problem with
E < V , as describedin Sec. IV A. Solid and dashed lines represent real and imag-inary contributions, respectively, for the analytical solution.Squares and circles denote corresponding numerical results.The shaded box represents the tunneling region. tion, Ψ( x ) which without loss of generality, must be real-valued. This in turn implies that the Ψ ± ( x ) are complexconjugates of each other, as in the bound state case dis-cussed in paper I. Indeed, one option is to simply applythe paper I decomposition to such problems. This ap-proach is discussed in detail in the Appendix, wherein itis shown to provide a natural extension of classical tra-jectories into the tunneling region.On the other hand, the trajectory hopping-based de-composition scheme offers a different, but also very nat-ural means to accomplish the same task—which has theadded advantage that all bipolar quantum potentials van-ish , except at x = 0. The idea is simply to treat all ex-pressions in Sec. II D 1 as being literally correct for thebelow barrier case as well, with the understanding thatthe requisite quantities need no longer be real. In partic-ular, p trans = i ¯ hκ [Eq. (22)] becomes pure positive imag-inary, implying that the transmitted trajectories “turn acorner” in the complex plane, and start heading off in thepositive imaginary direction, with speed ¯ hκ/m (Fig. 5).Along this path, the transmitted wave is an ordinaryplane wave; however, when analytically continued to thereal axis in the x > ◦ clockwise rotationin the complex plane), the familiar exponentially dampedform results (Fig. 2).For the reflected wave, Eq. (17) states that the re-flected “phase” remains unchanged. However, r refl isnow complex , leading to an effective phase shift of 2 δ ,where δ is defined in the Appendix [Eq. (21)]. For lo-calized wavepackets, a physical significance can be at-tributed to this phase shift, in both quantum mechanicsand electromagnetic theory; it is the source of the Goos-H¨anchen effect, a time delay observed in conjunction with total internal reflection. Consequently, in the time-dependent wavepacket context, it may be more appro-priate to associate the phase shift with time , rather thanwith s or r —specifically, with the delay time needed toaccrue sufficient action so as to compensate for the shift.For stationary states, however, such a time delay wouldbe inconsequential, because all trajectories are identicalapart from overall phase. Consequently, we do not con-sider such time delays explicitly in this paper, though wewill return to this issue in future publications.
3. multiple step systems
The most interesting case is that for which there aremultiple discontinuities, occurring at arbitrary locations x k (with k = 1 , , . . . , l ), and dividing up configurationspace into l + 1 regions, labeled A , B , C , etc. In each re-gion, the potential energy has a different constant value,i.e. V ( x ) = V A in region A , etc. From an optics pointof view, this system is analogous to a stack of differentmaterials, each with its own thickness, and index of re-fraction. Our primary focus in this paper will be squarebarrier/well systems for which l = 2, and V A = V C . How-ever, all of the present analysis extends to the more gen-eral case described above.In the standard time-independent picture, the solu-tion is obtained via a straightforward generalization ofEqs. (12), (13), and (14). However, even when compara-ble left-incident boundary conditions are specified as inSec. II D 1—i.e. A + = 1, and (for l = 2) C − = 0—theremaining coefficient values are fundamentally differentfrom those of the single-step case. To begin with, onlythe l ’th step exhibits the characteristics of a (locally) left-incident single-step solution; all other steps involve fournon-zero coefficients, corresponding locally to some su-perposition of left- and right-incident waves. Even moreimportantly, however, the expressions for the coefficientvalues as a function of system parameters in no way re-sembles Eq. (15); in particular, these now depend explic-itly on the x k values, as well as on V A , V B , etc. The sameis also true for the global P refl and P trans expressions, ascompared with Eq. (16).It is this dependence on the other steps that gives riseto the global nature of the time-independent solutions;i.e. the coefficient values at one step depend in principleon the properties of all of the other steps, no matter howfar away these might be located. Consequently, a reflec-tion probability as obtained from the A − value associatedwith the first, k = 1 step, cannot be determined withoutextending the analysis out to the final step at x = x l ,in the standard time-independent picture. On the otherhand, a primary goal of the time- dependent approach isto construct a completely local theory, for which local re-flection and transmission amplitudes associated with anygiven trajectory, as it encounters a given step k , depend only on the properties of the k ’th step (i.e. on x k , and onthe p or V values to the immediate left and right of x k ).0 FIG. 3: Bipolar trajectory plot for the square barrier prob-lem with
E > V , as described in Sec. IV B. One trajectoryin five is indicated in the figure. The black/gray solid linesindicate positive/negative LM trajectories, respectively. Theopen circles represent recombination points. The dashed linesdenote the two barrier edges, x = 0 and x = 1. In fact, from the point of view of the given trajectory, itmust be immaterial whether the potential contains othersteps or not—implying that the correct local relationsfor the spawned trajectories, if they exist at all, must beexactly those already specified in Eqs. (17) and (18).How is it possible that for stationary wavefunctions,whose time evolution is presumably trivial, an inherentlyglobal problem can be converted to a local one, simply byswitching from a time-independent to a time-dependentperspective? This is because of the bipolar decompo-sition, which provides each step with not one, but twosets of incident trajectories, one from the left, and onefrom the right. When there are multiple steps, not onlydoes this result in a non-trivial superposition for the re-sultant locally reflected and transmitted waves, but thetrajectories themselves are subject to multiple spawnings,which effectively enable them to traverse back and forthover the same regions of configuration space an arbitrarynumber of times (Fig. 3). This crucial feature ultimatelygives rise to the rich global scattering behavior observedeven in two-step systems. However, it is wholly missedby any time-independent treatment, even a bipolar one,which can only summarize the net superposition of allleft-traveling and right-traveling waves.We now discuss how the local time-dependent theorydescribed above gives rise to the correct stationary so-lutions, which is readily understood by invoking the rayoptics description introduced earlier. For simplicity anddefiniteness, we consider only the square potential case,which is optically analogous to say, a single pane of glasssurrounded by vacuum. If a single step gives rise to a sin-gle reflection, then two steps, like a pair of mirrors, resultsin an infinite number of reflections. The same is true of apane of glass, within which a single beam of light will bereflected back and forth at the edges an arbitrary number
FIG. 4: Seven snapshots of the superposition wavefunction,Ψ( x, t ), for the
E > V square barrier problem, as computedusing the numerical algorithm of Sec. III. The shaded box rep-resents the barrier region. All units are atomic. The evolvingdiscontinuities found at intermediate times in these curves de-note wavefronts for the newly created reflected/transmittedcomponents. Over time, the magnitudes of these disconti-nuities (i.e. corrections) become arbitrarily small, signifyingthat numerical convergence to the left-incident stationary so-lution has been achieved. of times. Of course, these reflections are not perfect; aportion of the incident flux always escapes as transmis-sion into the surrounding vacuum. Consequently, eachsuccessive internal reflection is exponentially damped, inaccord with Eq. (17).If the globally incident wave is incoming from the left,then at x , there are two contributions to Ψ B + . One con-tribution is the portion of the left-incident Ψ A + wave thatis locally transmitted through the first step. Apart froma phase factor, the resultant B + value would be givenby Eq. (15) if this were the only contribution. However,there is also a contribution that arises from the locallyreflected part of the right-incident wave, Ψ B − . This con-tribution is zero for a single step system, but of coursenon-zero in the multiple step case. Although the sec-ond contributing wave is right-incident, we can still useEqs. (17) and (18) to compute the contribution to Ψ B + ,as discussed in Sec. II D 1. For the second, k = l = 2 stepat x , there are only left-incident waves; consequently,Ψ C + and Ψ B − are obtained from a single source each,i.e. Ψ B + (Fig. 3).The above description refers to the stationary state re-sult, obtained by our gedankenexperiment in the largetime limit only. In practice this result would be achievedin stages. As in the previous examples, we imagine thatat time t = 0, one retains only those trajectories forwhich x ≤ x L < x . This one-sided trajectory restric-tion is somewhat analogous to continuous wave cavityring-down spectroscopy. When the wavefront first hitsthe first interface at x = x , there is partial reflectionand transmission, exactly identical to what would hap-1pen for a single step system. The reflected wavefrontpropagates beyond the left edge of the region of interestat x = x L , and for some time, the reflected amplitudepassing through this left edge is constant. The initiallytransmitted wavefront eventually reaches the second stepat x = x (i.e. the far side of the pane of glass), lead-ing to a second transmission into the C region, and asecond reflection back through the B region. Eventu-ally, the second transmitted wavefront reaches the rightedge of interest at x = x R , after which the transmittedamplitude remains constant for some time.Neither the globally transmitted nor reflected ampli-tudes for the times indicated above, as determined viamonitors at x = x L and x = x R respectively, are cor-rect. However, we have not yet described the steadystate solution. To do so requires an accounting of thesecond reflected wavefront, which eventually reaches the x = x step again, this time incident from the right.The resultant locally transmitted wave becomes an in-stantaneous second contribution to Ψ A − , and the locallyreflected wave plays the same role for Ψ B + . These newcontributions give rise to discontinuities in these waves,that subsequently propagate to the left and right, respec-tively (Fig.4). The new Ψ A − wave discontinuity eventu-ally reaches x = x L , where it is recorded by the monitor,giving rise to a sudden change in the reflection probabil-ity value.The Ψ B + discontinuity propagates to the second step,where it spawns new discontinuities in Ψ C + and Ψ B − .The former constitutes the border between first- andsecond-order transmitted waves, registered at sufficientlylater time by the monitor at x = x R . The latter, second-order Ψ B − wave heads back towards the first step, to giverise to third-order waves, with commensurate disconti-nuities, etc. In principle, this process continues indefi-nitely, resulting over time in global transmitted and re-flected waves of arbitrarily high order. However, Eq. (17)and the relation P refl + P trans = 1 ensure that the resultconverges to a stationary solution exponentially quickly.Moreover, since C − is necessarily zero throughout thisprocess, it is clear that the stationary state that is con-verged to is indeed the one corresponding to the desiredboundary condition of a globally incident wave that isincoming from the left.Note that in an actual optical system as describedabove, the spatial dimensionality is three rather thanone, and the incident wave would usually be taken atsome angle to the normal. If in addition, the beam hasa finite width, then one would observe separate reflectedbeams for each order, of exponentially decreasing bright-ness. The one-dimensional quantum case, however, isanalogous to a normal incident beam, for which all or-ders of reflection are superposed. In addition to provid-ing a pedagogical understanding of the dynamics that isvery much analogous to the optical example provided, thepicture above also suggests a practical numerical methodthat may be used to obtain stationary scattering statesof any desired boundary condition (via superposition of globally left- and right-incident wave solutions, obtainedindependently).Note that the “wavepacket dynamics” version of theray optics analogy may also be applied. In this case,the resultant initial square wavepacket is somewhat rem-iniscent of pulsed wave cavity ring-down spectroscopy. Once the wavepacket has penetrated the middle region B (i.e. the pane of glass), it reflects back and forth betweenthe two edges, with each reflection giving rise to a left- orright-propagating outgoing square wavepacket in region A or C , and a temporary interference pattern in region B . The amplitude of the central wavepacket dissipatesexponentially in time. All of this complicated behavioris indeed qualitatively observed in actual wavepacket dy-namics for such systems, but in the present context, isreconstructed entirely from a single stationary state. III. NUMERICAL DETAILS
In this section, we discuss several remaining issuespertaining to the numerical methods used to generateand propagate the various bipolar component waves, forthe examples discussed in Secs. II D and IV. For allof these examples, the numerical algorithm used corre-sponds to the gedankenexperiment with one-sided trun-cation, i.e. to continuous wave cavity ring-down. Inessence, this consists of just two basic operations: (1)each piece-wise bipolar component of the wavefunction[i.e. Ψ A ± ( x ), Ψ B ± ( x ), etc.] is independently propagatedin time over its appropriate region of space, using thestandard QHEMs and QTMs; (2) whenever a trajectoryreaches a turning point, it is immediately deleted, and re-placed with two new trajectories, spawned in the appro-priate locally transmitted and reflected component LMs.The first operation above, i.e. QTM propagation ofthe wavefunction components, is very straightforward.Note that for simplicity, we have throughout this paperused time-independent expressions for Ψ( x ) and its com-ponents, but in reality these evolve over time—even forstationary states, via ˙ s = ∂s ( x, t ) /∂t = − E . We havetherefore been rather lax in distinguishing Hamilton’sprinciple function from Hamilton’s characteristic func-tion, although from a trajectory standpoint, it is alwaysthe former that is implied. Since each component is sta-tionary in its own right, the time evolution of the hy-drodynamic fields is governed by the quantum station-ary Hamilton-Jacobi equation (QSHJE), rather than theQHJE. Moreover, the fact that the piecewise r is con-stant implies that the component quantum potentials are zero , resulting in classical HJE’s and trajectories. Theseconclusions are trivially correct for the present paper,for which all components are plane waves; however, thearguments also extend to arbitrary continuous potentialsystems. From a numerical perspective, the use of classicaltrajectories offers many advantages over a conventionalQTM propagation. To begin with, the trajectories them-2selves are always smooth if V ( x ) is smooth, resulting infar fewer trajectories and larger time steps than wouldotherwise be the case. Even more importantly, however,since the quantum potential is not required, there is noneed to compute on-the-fly numerical spatial derivativesof the local hydrodynamic fields. Consequently, for agiven component, the trajectories are completely inde-pendent and need not communicate—again resulting infewer of them. Indeed, it is possible to perform an essen-tially exact computation using only a single trajectoryper wavefunction component . This feature is particularlyimportant for the very frontmost trajectory of the initialensemble, which for brief periods at later times, will (viaspawning) come to be the only trajectory to occupy agiven component LM. The subsequent evolution of theselone trajectories does not require the presence of nearbytrajectories.The spawning of new trajectories, i.e. operation (2)above, also bears further discussion. In principle, thisis always achieved via application of Eqs. (17) and (18).For the above barrier case, R and T are both real, ensur-ing the reality of r and s for the spawned trajectories—although in the case of a down step, R <
0, resulting ina negative r refl . This is in accord with the conventionsdiscussed in paper I. However, in this paper, we find itnumerically convenient to adopt the more usual r > r refl , it isreplaced with − r refl , and π is added to s refl . A similar,but more complicated modification is also applied to thebelow-barrier trajectories, for which Eqs. (17) and (18)yield complex amplitudes. In this case, the phase shift is2 δ , as discussed in Sec. II D 2 and the Appendix.For a single step system, the algorithm is now essen-tially complete. At the initial time t = 0, a variablenumber of particles (or synonymously, grid points) aredistributed uniformly along the Ψ A + manifold, to theleft of x = x L . The extent of these points must be largeenough that at the end of the propagation, there are stillΨ A + grid points that have not yet reached x L . The gridspacing is mostly arbitrary, but must be small enoughthat at sufficiently later times, there is always at leastone trajectory per component LM. The propagation isconsidered complete when the reflected and transmittedwavefronts travel beyond x L and x R , respectively.For multiple step systems, the situation is similar, butsomewhat more complex. The primary new feature isthe recombination of wavefunction components arisingfrom two sources, i.e. from two locally incident wavescoming from opposite directions. The present algorithmwould seem to yield two subcomponent wavefunctionsfor every component, each with its own set of trajecto-ries. If left “unchecked,” this would lead to undesirablefurther multifurcations for higher orders/later times. Asimple solution would be to propagate each subcompo-nent long enough that there is at least one trajectoryfor each, then extrapolate the corresponding subcom-ponent wavefunctions to a common position, where anew trajectory is constructed for the superposed com- ponent wavefunction, which is then propagated in lieuof the subcomponent trajectories. This requires dynam-ical fitting (see below), or at the very least, extrapo-lation. Although these numerical operations would bevery stable in the present context, to rule these out al-together as sources of error in Sec. IV, we have adopteda much simpler approach—i.e. the grid spacing is cho-sen such that trajectories from the two component wavesincident on a given step always arrive at the same time(Fig. 3). The corresponding subcomponent wavefunc-tion values are then simply added together when form-ing the spawned trajectory. Adopting once again the r > ± , the corresponding field values are then obtained via r = p Ψ ∗± Ψ ± and s = arctan [Im(Ψ ± ) / Re(Ψ ± )].As discussed in Sec. II D 3, multiple step systems allowfor infinite reflections that perpetually modify | Ψ( x, t ) | ,in principle for all time. In practice however, there isexponential convergence within the region of interest, x L ≤ x ≤ x R , so that one would not run the calcula-tion indefinitely, but only until the desired accuracy isreached. Accurate “error bars” on the computed global P refl and P trans values are conveniently provided by themagnitudes of the most recent discontinuous jumps asrecorded by the monitors at x L and x R . Note that thenumber of digits of accuracy scales only linearly withpropagation time. However, the rate of convergence de-pends on the energy value. Near the barrier height, inparticular, convergence may take quite a long time, asthe exponent is close to zero. For all other energies, onlya few “cycles” should be required, depending on the levelof accuracy desired.If in addition to reflection and transmission probabili-ties, the actual stationary solution over the region of in-terest is also desired, then it is necessary to reconstructΨ( x ). This is obtained from the final grid, after the prop-agation is finished, using a multiple step generalization ofEq. (14). The first step is to reconstruct the componentwavefunctions Ψ A ± ( x ), Ψ B ± ( x ), etc., via interpolationor fitting of the hydrodynamic field values from the cor-responding dynamical grid points onto a much finer com-mon grid (used e.g. for plotting purposes). The secondstep is to linearly superpose the ± components onto theplotting grid, and to assemble the pieces together overthe coordinate range of interest. For the discontinuoussystems considered here, the number of dynamical gridpoints per component can be as small as one—i.e. muchsmaller, even, than the number of wavelengths! To ourknowledge, such performance has never been achievedpreviously by a QTM; however, it does require that theplotting grid be much finer than the dynamical grid, e.g.at least several points per wavelength, in order to ade-quately represent the interference fringes of the the su-perposed solution, Ψ( x ).3 IV. RESULTS
In this section, we apply the numerical algorithm pre-viously described to three different applications: the up-step potential, the square barrier, and the square well.
A. Up-step Potential
The first system considered is the up-step potential,i.e. Eq. (11) with V >
0. Since there are no multiplereflections, this is in principle a trivial application forthe current algorithm; it therefore serves as a useful nu-merical test. Both above barrier (Sec. II D 1) and belowbarrier (Sec. II D 2) energies are considered. We choosemolecular-like values for the constants, i.e. V = 0 . m = 2000 a.u. The left and right edges ofthe region of interest are taken to be x L = − . x R = 1 . t = 0, 51trajectory grid points are distributed uniformly over theinterval − ≤ x ≤ − .
06 a.u.). Thisnumber is far greater than what would be needed for dy-namical purposes, but is chosen so as to avoid construc-tion of a separate plotting grid (Sec. III). The hydrody-namic field functions for the initial Ψ A + ( x ) wavepacketover the above interval are taken to be r ( x ) = 1 a.u. − and s ( x ) = √ mEx .For the above barrier calculation, the energy E =2 V = 0 .
018 hartree was used. The trajectory propa-gation and termination were performed exactly as de-scribed in Sec. III. The real and imaginary parts of allthree resultant wavefunction components [i.e. Ψ A ± ( x )and Ψ B + ( x )] at the final time, t = 550 a.u. are pre-sented in Fig. 1. All three components exhibit the de-sired plane wave behavior, e.g. no interference is evidentwithin a given component. The resultant Ψ( x ) does ex-hibit interference in the A region, however, arising fromthe superposition of Ψ A + and Ψ A − .For the below barrier calculation, the system was givenan energy equal to one half of the barrier height, i.e. E = V / . B is achieved,not through a quantum potential, but via analytic con-tinuation. At sufficiently large time ( t = 1100 a.u.) thefinal wavefunction is reconstructed from the components,i.e. Ψ A ( x ) = Ψ A + ( x ) + Ψ A − ( x ), and Ψ B ( x ) = Ψ B + ( x ).The real and imaginary parts of the reconstructed wave-function are presented in Fig. 2. In the figure, squaresand circles denote the numerical results obtained via thepresent algorithm, whereas the solid and dashed lines rep-resent the well-known analytic solutions. The agreementis essentially exact. Note that the real and imaginaryparts are in phase throughout the coordinate range—i.e., S ( x ) is a constant, so apart from a phase factor, Ψ( x ) isreal. Note also that the tunneling region exhibits thedesired exponential decay. B. Square Barrier
The second system considered is the square barrier.This is a two-step potential ( l = 2), with V A = V C =0, and V B = V >
0. The two steps comprise the leftand right edges of the barrier, at x = 0 and x = w ,respectively. The constants are chosen as follows: V =0 .
018 hartree; m = 2000 a.u.; w = 1; x L = − x R =2 a.u. Initially, 75 trajectory grid points are distributeduniformly over the interval − ≤ x ≤ − .
05 a.u.), which again, is far more than are dynamicallyrequired. The same initial hydrodynamic field functionsare used as in Sec. IV A.Both above barrier (
E > V ) and below barrier ( E 036 hartree. Once again, the trajectorypropagation and termination were performed exactly asdescribed in Sec. III. In order to converge P trans to 10 − ,a propagation time of 3000 a.u. was required. This corre-sponds to 3 complete cycles, i.e. a 3rd-order calculation.Figure 3 is a plot of the quantum trajectories for thiscalculation, in which every fifth trajectory for each ofthe five component wavefunctions is indicated. Trajec-tory spawning at the two steps is very clearly evident,as is recombination of pairs of incident waves (indicatedby circles). On the whole, this figure demonstrates allof the anticipated analogues with ray optics, i.e. paralleltrajectories, reflection and refraction.In Fig. 4, the time evolution of the superposition stateΨ( x, t ) is represented, via snapshots of the real and imag-inary parts at seven different times. At t = 0 a.u., theincident wavefront is located at x = − t = 250a.u., the wavefront has spawned Ψ A − ( x ) and Ψ B + ( x )trajectories; the former gives rise to the kink (really adiscontinuity) somewhat to the left of the first step. By t = 550 a.u. and t = 650 a.u., the Ψ A − ( x ) wavefront hasmoved outside the region of interest, though the Ψ B + ( x )wavefront has not quite reached the second step. Afterit does so, two new wavefronts are propagated along theΨ C + ( x ) and Ψ B − ( x ) LMs (e.g. t = 800 a.u.), the for-mer of which propagates beyond the region of interest by t = 900 a.u. Subsequent discontinuity magnitudes be-come exponentially smaller, so that at sufficiently largetime (i.e. t = 2000 a.u.), the resultant Ψ( x ) has con-verged to the correct stationary solution.For the below barrier case, E = V / . 009 hartree, x L = − . w = 0 . − level of convergence ( t = 1400 a.u., or2 cycles). Fig. 5 indicates how the tunneling dynamicsis achieved. After the wavefront hits the first step, thetransmitted Ψ B + wave is propagated along the imaginaryaxis ( iy ), until the point y = w is reached. When thisoccurs, it is necessary to analytically continue Ψ B + down4 FIG. 5: Schematic of the algorithm used to propagate tra-jectories into the classically forbidden region of the E < V square barrier problem, as discussed in Secs. III and IV B. to the real axis, in order to compute amplitudes for thenew Ψ C + ( x ) and Ψ B − ( x ) trajectories that are spawnedat the second step. The latter component propagates inthe (negative) imaginary direction, along w − iy ′ , until y ′ = w , at which point analytic continuation is once moreapplied (this time for the first step), and the pattern re-peated.Seven snapshots of the superposition density, | Ψ( x, t ) | , are displayed in Fig. 6. The initial density—equal to just the | Ψ A + ( x ) | density—is uniform. Afterthe wavefront encounters the first step, a reflectedΨ A − ( x ) emerges, giving rise to clearly evident inter-ference in the region A . The Ψ B ( x ) = Ψ B + ( x ) waveis at this stage perfectly exponentially damped. Uponencountering the second step, a second contribution,Ψ B − ( x ) emerges; however, this first-order correction isalready extremely small, owing to the large amount oftunneling that has occurred. The global transmittedwave, Ψ C + ( x ), though small, is clearly seen to haveuniform density.In addition to the two detailed trajectory calculationsdescribed above, we computed P refl and P trans for a largerange of w and E values, so as to fully explore (withoutloss of generality) the entire range of the square barrierproblem. The numerical results are presented, and com-pared with known analytical values, in Fig. 7. Twoaspects of this study bear comment. First, for all w and E values considered, the computed P refl and P trans valuesagree with the exact values to within an error compara-ble to that predicted by the level of numerical conver-gence. In particular, the oscillatory energy dependenceis perfectly reproduced. Second, the closer the barrierpeak is approached from either above or below in energy,the longer the time required to achieve a given level ofconvergence, as predicted in Sec. III. In particular, forthe calculations closest to the barrier peak, 5-7 cycleswere required in order to approximately maintain a 10 − convergence of the transmission probability. Although agreater number of particles are required in this case, thisposes no great limitation in practice, since one would pre-sumably never require a calculation precisely at the peak FIG. 6: Seven snapshots of the superposition density, | Ψ( x, t ) | , for the E < V square barrier problem, as com-puted using the numerical algorithm of Sec. III. The shadedbox represents the barrier region. All units are atomic. Notethat interference manifests in the incident (left) region only af-ter some incident trajectories have struck the left barrier edge,causing reflected trajectories to be created (i.e. just prior to t = 195, initially). The most advanced reflected trajectory de-fines the reflected wavefront, manifesting as the left-movingdiscontinuity, e.g. at t = 195 and t = 287.FIG. 7: Transmission and reflection probabilities as a func-tion of energy, for square barrier potentials of three differentwidths, w , as discussed in Sec. IV B. Solid lines denote an-alytical results; open/closed circles denote numerical results,as obtained via algorithm of Sec. III. The vertical dot-dashedlines represent the barrier height, i.e. E = V . energy.5 FIG. 8: Transmission and reflection probabilities as a func-tion of energy, for square well potentials of three differentwidths, w , as discussed in Sec. IV C. Solid lines denote an-alytical results; open/closed circles denote numerical results,as obtained via algorithm of Sec. III. C. Square Well As the final system, we consider the square-well po-tential, i.e. the square barrier but with V < 0. In thescattering state context, there is no tunneling for thissystem, but in other respects it resembles the squarebarrier. From an optics point of view, the square wellcorresponds to a central medium with larger index of re-fraction than its surroundings, whereas the square barriercorresponds to a smaller index of refraction, giving riseto the possibility of total internal reflection (i.e. tunnel-ing). The parameters are as in Sec. IV B, except that V = 0 . 009 hartree, and three different w values are con-sidered: w = 2 a.u., w = 4 a.u., and w = 16 a.u.As the time evolution and trajectory pictures are sim-ilar to those of the previous sections, we focus only onthe P refl and P trans calculations, for which once again, alarge range of energies was considered (0 . < E < . − . The energy-resolvedreflection and transmission probabilities are presented inFig. 8.As in the square barrier case, excellent agreement isachieved with the exact analytical results, i.e. on the or-der of the level of convergence. This is true despite thefact that the square well energy curves are decidedly moreoscillatory than the square barrier curves—particularlyfor wide barriers, for which the E dependence is verysensitive indeed. One particularly important feature ex-hibited by the exact curves is the so-called Ramsauer- Townsend effect, i.e. the phenomenon of 100% trans-mission and zero reflection, even at very low energies.This occurs when sin(2 k B w ) = 0, and may be regardedas a purely quantum mechanical resonance phenomenon.Yet it is reproduced perfectly here in the bipolar decom-position, using classical trajectories . Indeed, the bipolarpicture provides an interesting physical explanation, i.e.the right-incident and left-incident waves of the first stepgive rise to spawned contributions to Ψ A − ( x ) that ex-actly cancel each other out via destructive interference. V. SUMMARY AND CONCLUSIONS As described in paper I, the Schr¨odinger equation islinear, yet the equivalent QHEM—obtained via substitu-tion of the Madelung-Bohm ansatz into the Schr¨odingerequation—are not. This aspect of Bohmian mechanicssuggests that it can be beneficial, both from a pedagog-ical and a computational perspective, to apply a suit-able bifurcation (or “multifurcation”) to the wavefunc-tion prior to applying the QHEM. Indeed, following thepaper I CPWM bipolar decomposition for 1D stationarybound states, the quantum trajectories become morewell-behaved and classical-like, in precisely the limit inwhich there are more nodes, and the usual unipolar calcu-lation breaks down. Moreover, the resultant componentLMs admit a natural physical interpretation in terms ofthe corresponding semiclassical LMs.In the generalization to the stationary scattering statesconsidered here, a somewhat different bipolar decompo-sition is found to be required. The new decomposition isstill unique, at least for discontinuous potentials. How-ever, the resultant components Ψ ± ( x ) are no longer so-lutions to the Schr¨odinger equation in their own right,as a result of which their time evolution is coupled. Thenew scheme—though fundamentally different from theold one—nevertheless bears a correspondence to a mod-ified version of semiclassical theory appropriate for scat-tering systems. Curiously, this semiclassical modificationis not simply a higher order treatment in ¯ h ; if it were,the corresponding exact quantum modification consid-ered here would not exist.In any event, the new decomposition also gives riseto its own physical interpretation, specific to the scat-tering context. In particular, for right-incident bound-ary conditions, the left and right asymptotes of Ψ + ( x )respectively represent incident and transmitted waves,whereas the left asymptote of Ψ − ( x ) represents the re-flected wave [the right asymptote of Ψ − ( x ) approacheszero]. That these conceptually useful asymptotic bipo-lar assignments—found even in the most elementarytreatments of scattering—may be extended throughoutconfiguration space [even for continuous potentials (pa-per III)] represents an important leap forward, especiallyfor QTMs.Another pedagogically and numerically useful devel-opment from this approach is the inherently time-6dependent ray optics interpretation that naturally arises,particularly in the discontinuous potential context. Theray optics approach is anticipated to be a relevant guid-ing force in subsequent generalizations of the presentmethodology, i.e. to continuous, multidimensional po-tential systems, and—as per the discussion in Secs. II C 2and II D 1—for non-stationary wavepacket dynamics.This approach also provides a much simpler trajectory-based explanation of scattering terminology as appliedin a stationary context—e.g. why the “reflected wave” isso-called, despite being present from the earliest times—than those traditionally used. The discontinuous potential applications consideredin this paper—hard wall (Sec. II C 2), step potential(Secs. II D 1, II D 2, and IV A), and square barrier/well(Secs. II D 3, IV B, and IV C)—are significant for severalreasons. To begin with, these are the first discontinuousapplications of a genuine QTM calculation that have everbeen performed, to the authors’ knowledge. The singularderivatives associated with discontinuous potential func-tions would wreak havoc with standard numerical differ-entiation routines. Second, the use of a time- dependent method for stationary, or time- independent , applications,is also significant. Ordinarily, the time dependence ofstationary states is regarded as trivial. In the presentcontext, this is true in a sense for the hard wall and steppotential systems, because the correct answer is “builtin” the method itself. For multiple step systems, how-ever, the dynamical truncated wave approach, i.e. thegedankenexperiment introduced in Sec. II C 2 and furtherdeveloped in later sections, yields decidedly nontrivial re-sults.In particular, the algorithm uses only single step scat-tering coefficients to obtain global scattering quantitiesfor multiple step systems. In effect, the time depen-dent nature of this approach allows computation of global properties using a completely local method. Not onlywere exact quantum results obtained for a full range ofsystem parameters, but the numerical resources neces-sary to achieve this—i.e. the number of trajectories andtime steps—were decidedly minimal. Indeed, the algo-rithm lives up to the promise made in paper I, of per-forming an accurate quantum calculation with fewer tra-jectories than nodes —a prospect virtually unheard of ina unipolar context.In future publications, we will naturally attempt togeneralize the methodology described here and in pa-per III, for the type of multidimensional time-dependentwavepacket dynamics relevant to chemical physics appli-cations. In this context, the scattering version of theCPWM decomposition developed here is an absolutelyessential first step, as reactive scattering is the underpin-ning of all chemical reactions. Additional discussion andmotivation will be provided in paper III. ACKNOWLEDGEMENTS This work was supported by awards from The WelchFoundation (D-1523) and Research Corporation. The au-thors would like to acknowledge Robert E. Wyatt andEric R. Bittner for many stimulating discussions. DavidJ. Tannor and John C. Tully are also acknowledged. Appendix: Bipolar decomposition of semi-boundstates As discussed in Secs. II B, II C 2, and II D 2, semi-bound stationary states in 1D are bounded on one sideonly, as a result of which they are real-valued and singly-degenerate, like bound states. Consequently, they areamenable to the CPWM bipolar decomposition schemeintroduced in paper I. In this appendix, we apply this de-composition to two semi-bound systems: the hard wallsystem, and the below-barrier up-step system. A. Hard wall system From Ref. 24, the most general bipolar decomposi-tion of a hard wall stationary state—corresponding toSec. II A 2 condition (1) only—is found to satisfy − cot [ s ( x ) / ¯ h ] = (cid:18) mF ¯ h (cid:19) (cid:20) − cot( kx ) k + B (cid:21) , (19)where s ( x ) = s + ( x ) = − s − ( x ), and r + ( x ) = r − ( x ) isobtained from s ( x ) via Eq. (7) (without “sc” subscripts).The arbitrary parameters F and B are the invariantflux and median action parameters associated with con-ditions (2) and (3), respectively, although the definitionof F has been changed slightly to account for the scat-tering normalization convention, r + ( x ) = 1. Note that only the semiclassical values for these parameters yieldsa solution that satisfies the correspondence principle inthe large action (i.e. k ) limit. In particular, the choice B = 0 and F = ¯ hk/m yields the desired semiclassicalresult, s ( x ) = ¯ hkx ; all other choices exhibit undesirableoscillatory behavior in r ± ( x ), s ± ( x ), and q ± ( x ). B. Up-step system For the hard wall system considered above—which isjust the special case of the up-step potential in the limit V → ∞ —exact agreement is achieved between semiclas-sical and quantum LM’s in the x < V values—i.e. for general below-barrier up-stepstationary states—there is of course also tunneling intothe forbidden region, which must be accounted for. Thepaper I bipolar decomposition therefore results in LM’sthat span the entire coordinate range −∞ < x < ∞ .7These LMs are given by the following analytical expres-sion: p ( x ) = (cid:26) ¯ hk for x ≤ hκe κx sin 2 δ e κx − e κx cos 2 δ for x > , (20)where δ is given bytan δ = κk = s(cid:18) V E (cid:19) − , (21)and κ = p m ( V − E ) . (22)The p ( x ) LM function of Eq. (20) is continuous every-where, including at the potential discontinuity at x = 0.In the A region, it agrees exactly with the semiclassicalsolution; in the B region, it decays exponentially to zero.The paper I approach thus yields a very natural way toextend trajectories into the tunneling region. Note that the quantum potential in this region is not zero; indeed,it exhibits a discontinuity at x = 0 that exactly balancesthat of V ( x ) itself, so that the bipolar modified potentialis continuous across the step. Unlike the above-barriercase, the paper I solution does not manifest oscillatorybehavior in the large action limit, and so this approachwould at first glance appear to be ideal. There are tworeasons, however, why it is not pursued here. The firstreason is that r ( x ) diverges asymptotically as x → ∞ ,which according to preliminary numerical investigations,appears to lead to numerical instabilities for completelyQTM-based propagation schemes. Second, if the barrierwere to fall off again at larger x values, so that the tun-neling region were finite, then the asymptotic behaviorwould be once again undesirably oscillatory. This wouldbe the case, for example, for the below-barrier energiesof the square barrier system of Sec. IV B. ∗ Electronic address: [email protected] J. M. Bowman, J. S. Bittman, and L. B. Harding, J. Chem.Phys. , 911 (1986). Z. Baˇci´c and J. C. Light, Annu. Rev. Phys. Chem. , 469(1989). B. Poirier and J. C. Light, J. Chem. Phys. , 4869(1999). B. Poirier and J. C. Light, J. Chem. Phys. , 211 (2000). H.-G. Yu, J. Chem. Phys. , 8190 (2002). X.-G. Wang and T. Carrington, Jr., J. Chem. Phys ,101 (2003). R. Dawes and T. Carrington, Jr., J. Chem. Phys. , 726(2004). B. Poirier, J. Theo. Comput. Chem. , 65 (2003). B. Poirier and A. Salam, J. Chem. Phys. , 1690 (2004). B. Poirier and A. Salam, J. Chem. Phys. , 1704 (2004). H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem.Phys. Lett. , 73 (1990). U. Manthe, H.-D. Meyer, and L. S. Cederbaum, J. Chem.Phys. , 3199 (1992). C. L. Lopreore and R. E. Wyatt, Phys. Rev. Lett. , 5190(1999). F. S. Mayor, A. Askar, and H. A. Rabitz, J. Chem. Phys. , 2423 (1999). R. E. Wyatt, Chem. Phys. Lett. , 189 (1999). R. E. Wyatt and E. R. Bittner, J. Chem. Phys. , 8898(2001). R. E. Wyatt and K. Na, Phys. Rev. E , 016702 (2001). R. E. Wyatt, Quantum Dynamics with Trajectories: Intro-duction to Quantum Hydrodynamics (Springer, New York,2005). D. Bohm, Phys. Rev. , 166 (1952). D. Bohm, Phys. Rev. , 180 (1952). T. Takabayasi, Prog. Theor. Phys. , 341 (1954). E. Madelung, Z. Phys. , 322 (1926). J. H. van Vleck, Proc. Natl. Acad. Sci. U.S.A. , 178(1928). B. Poirier, J. Chem. Phys. , 4501 (2004). Y. Zhao and N. Makri, J. Chem. Phys. , 60 (2003). R. E. Wyatt, C. L. Lopreore, and G. Parlant, J. Chem.Phys. , 5113 (2001). E. R. Bittner, J. B. Maddox, and I. Burghardt, Int. J.Quantum Chem. , 313 (2002). D. V. Shalashilin and M. S. Child, J. Chem. Phys. ,10028 (2000). I. Burghardt and L. S. Cederbaum, J. Chem. Phys. ,10303 (2001). I. Burghardt and L. S. Cederbaum, J. Chem. Phys. ,10312 (2001). C. J. Trahan and R. E. Wyatt, J. Chem. Phys. , 7017(2003). A. Donoso and C. C. Martens, J. Chem. Phys. , 6309(2002). E. R. Bittner, J. Chem. Phys. , 6309 (2002). K. H. Hughes and R. E. Wyatt, J. Chem. Phys. , 4089(2004). B. K. Kendrick, J. Chem. Phys. , 5805 (2003). D. K. Pauler and B. K. Kendrick, J. Chem. Phys. , 603(2004). S. Garashchuk and V. A. Rassolov, J. Chem. Phys. ,1181 (2004). K. H. Hughes and R. E. Wyatt, Phys. Chem. Chem. Phys. , 3905 (2003). S. Garashchuk and V. A. Rassolov, J. Chem. Phys. ,8711 (2004). E. R. Floyd, Physics Essays , 135 (1994). M. R. Brown, arXiv:quant-ph/0102102 (2002). D. Babyuk and R. E. Wyatt, J. Chem. Phys. , 9230(2004). B. Poirier, J. Chem. Phys. , (submitted). M. V. Berry and K. V. Mount, Rep. Prog. Phys. , 315(1972). N. Fr¨oman and P. O. Fr¨oman, JWKB Approximation (North-Holland, Amsterdam, 1965). J. Heading, An Introduction to Phase-integral Methods (Methuen, London, 1962). P. R. Holland, The Quantum Theory of Motion (Cam-bridge University Press, Cambridge, 1993). J. B. Keller and S. I. Rubinow, Ann. Phys. , 24 (1960). V. P. Maslov, Th´eorie des Perturbations et M´ethodesAsymptotiques (Dunod, Paris, 1972). R. G. Littlejohn, J. Stat. Phys. , 7 (1992). J. R. Taylor, Scattering Theory (John Wiley & Sons, Inc.,New York, NY, 1972). B. Poirier and T. Carrington, Jr., J. Chem. Phys. , 17(2003). B. Poirier, Found. Phys. , 1191 (2000). J. C. Tully, J. Chem. Phys. , 562 (1971). J. D. Jackson, Classical Electrodynamics , 2nd ed. (JohnWiley & Sons, New York, NY, 1975). L. Brillouin, Ann. Phys. , 177 (1914). J. O. Hirschfelder, A. C. Christoph, and W. E. Palke, J.Chem. Phys. , 5435 (1974). M. D. Wheeler, S. M. Newman, A. J. Orr-Ewing, andM. N. R. Ashfold, JJ. Chem. Soc. Faraday Trans. , 337(1998). S. Gasiorowitz,