Reconciling Semiclassical and Bohmian Mechanics: III. Scattering states for continuous potentials
aa r X i v : . [ qu a n t - ph ] F e b Reconciling Semiclassical and Bohmian Mechanics:III. Scattering states for continuous 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,such that the components Ψ and Ψ approach their semiclassical WKB analogs in the large actionlimit. The corresponding bipolar quantum trajectories, as defined in the usual Bohmian mechanicalformulation, are classical-like and well-behaved, even when Ψ has many nodes, or is wildly oscillatory.A modification for discontinuous potential stationary stattering states was presented in a secondpaper [J. Chem. Phys. !!! !!!! (2005)], whose generalization for continuous potentials is given here.The result is an exact quantum scattering methodology using classical trajectories. For additionalconvenience in handling the tunneling case, a constant velocity trajectory version is also developed. I. INTRODUCTION
This paper is the third in a series investigating the useof “counter-propagating wave methods” (CPWMs) for solving the Schr¨odinger equation exactly. CPWMsare a particular variant of the more general multipolardecomposition methods, wherein the wavefunction Ψ isdecomposed into two or more components, Ψ j . Thus, fora two-term, or “bipolar” decomposition, Ψ = Ψ + Ψ .This rather trivial-seeming procedure can be quite advan-tageous in the context of quantum trajectory methods(QTMs), i.e. trajectory-based numer-ical techniques for performing exact quantum dynamicscalculations, based on Bohmian mechanics. Conventional QTMs use a single-term or “unipolar” rep-resentation of the wavefunction from which all otherquantities, such as the quantum trajectories themselves,are uniquely determined. Multipolar decomposition, onthe other hand, can lead to radically different QTM be-havior for the individual Ψ j components, owing to thefact that the Bohmian equations of motion are non-linear. As applied to wavepacket dynamics for reactive scat-tering systems, QTMs suffer from a significant and well-known numerical drawback, which to date, precludes acompletely robust application of these methods. Namely,QTMs are numerically unstable in the vicinity of ampli-tude nodes and “quasi-nodes” (i.e. rapid oscillations), owing to singularities in the “quantum potential,” Q ,which together with the classical potential, V , determinesthe quantum trajectories. In the reactive scattering con-text, such behavior is always observed, due to interfer-ence between the incident and reflected waves. On theother hand, if the latter two contributions to the total Ψwere somehow separated, and associated with two differ-ent interference-free Ψ j components, the node problemmight well be circumvented. If, in addition, the Ψ j com-ponent field functions were smooth and slowly-varying,far fewer QTM trajectories and time-steps might be re-quired than for Ψ itself—although the latter could bereconstructed at any desired time simply via linear super- position of the components. These numerical advantagesthus provide significant motivation for consideration ofthe bipolar approach. A much more detailed discussionmay be found in the first two articles of this series, pa-per I (Ref. 1) and paper II (Ref. 2).The most obvious aspect of any bipolar decomposi-tion, including those restricted along the lines of the pre-ceding paragraph, is that it is not unique. Thecovering function method, for instance, treats Ψ asthe difference between two very large-amplitude compo-nents, thus “diluting” the effects of interference. Theone-dimensional (1D) CPWM approach, on the otherhand, regards the bipolar decompositionΨ = Ψ + + Ψ − (1)as a superposition of right- and left-traveling counter-propagating waves, Ψ ± . For stationary states at least,the Eq. (1) decomposition is defined such that the Ψ ± components correspond to semiclassical WKB approxi-mations, Ψ sc ± , in the large-action limit. In addition toproviding pedagogical value (semiclassical and Bohmianmechanics can not be so reconciled in a unipolar context),the semiclassical field functions are typically smooth andslowly-varying, i.e. the semiclassical-like CPWM com-ponents Ψ ± provide the desirable numerical advantagesdescribed in the preceding paragraph.In paper I (Ref. 1), a unique CPWM bipolar de-composition was determined for bound stationary eigen-states of arbitrary 1D Hamiltonians. The resultant Ψ ± field functions are smooth and interference-free, and ap-proach the WKB approximations in the large-action limit(within the classically allowed region of space) as desired.Moreover, the quantum potentials q ± become vanishinglysmall in this limit, so that the bipolar quantum trajecto-ries approach classical trajectories, and only a small num-ber are required for a numerical propagation, regardlessof excitation energy, E . In contrast, the unipolar Ψ ex-hibits motionless trajectories, and an arbitrarily increas-ing number of nodes in the large-action/energy limit. Re-sults were presented for both the harmonic and Morseoscillator potentials.In paper II (Ref. 2), the 1D CPWM ideas were mod-ified somewhat for stationary scattering states of dis-continuous potentials. Although the CPWM decomposi-tion of paper I is uniquely specified for any arbitrary 1Deigenstate—bound or scattering—and in the bound case,always satisfies the correspondence principle, the non- L nature of the scattering states is such that the paper Idecomposition generally does not satisfy correspondenceglobally. As discussed in Sec. II A, this requires thatglobal modifications must be made in order to enable acorrespondence between Ψ sc ± and Ψ ± . These are suchas to lead to substantial differences between the density functions | Ψ sc ± ( x ) | and | Ψ ± ( x ) | , although the trajecto-ries are identical. Moreover, the resultant Ψ ± are foundto correspond to the familiar “incident,” “transmitted,”and “reflected” waves of traditional scattering theory in the appropriate asymptotic limits. From the time-dependent standpoint, reflection was found to be due totrajectory hopping from one CPWM component to theother, as can be naturally understood using a simple rayoptics analogy. The method was applied to several ele-mentary discontinuous potential systems, including thesquare barrier/well.In the present paper (paper III), the Ref. 2 formu-lations are generalized to incorporate both continuous and discontinuous potential systems. Once again, ananalogy is drawn from semiclassical mechanics, albeit a“sophisticated” version less frequently considered(Sec. II A). As in paper II, a time-dependent methodbased on ray optics is developed for computing station-ary states of any desired energy and boundary condi-tions. For the most part, the discussion of the precedingparagraph still applies, but some additional key pointsshould be emphasized. First, the bipolar decompositionnow provides a sensible definition of “incident,” “trans-mitted,” and “reflected” waves throughout all space , notjust asymptotically. Second, the explicit hopping of tra-jectories from one CPWM component to the other is re-placed with a coupling term in the time-evolution equa-tions. Third, the trajectories become completely classi-cal . Finally, an alternative methodology is also devel-oped, based on the use of constant velocity trajectories ,which can be readily applied to barrier tunneling situ-ations. The new methods are found to be remarkablyefficient, accurate, and robust across a diverse range of1D test potentials and system energies (Sec. IV).The paper is organized as follows. A derivation anddiscussion of the time-evolution equations for the CPWMbipolar components Ψ ± , both for classical and constantvelocity trajectories, are presented in Sec. II and the Ap-pendices. Sec. III provides numerical details of the vari-ous bipolar algorithms used to compute stationary states.Results are presented in Sec. IV for four benchmark ap-plications: Eckart barrier; square barrier; uphill ramp;double-Gaussian barrier. For the first three, these arecompared with known analytic solutions. Concluding re-marks, including prospects for future development, maybe found in Sec. V. II. THEORYA. Semiclassical Approximations
Let V ( x ) be the potential energy for a 1D scatter-ing system, and E the energy of some stationary stateΨ( x ) such that E > V ( x ) for all x . In reality, Ψ( x ) al-ways manifests some reflection, even though the energy isabove the potential barrier. However, basic WKB theorypredicts zero reflection in this case, as the classical trajec-tories do not turn around. This is also evident from theform of the two basic WKB counter-propagating wavesolutions, Ψ sc ± ( x ) = r sc ( x ) e ± is sc ( x ) / ¯ h , (2)where r sc ( x ) = s mFs ′ sc ( x ) and s ′ sc ( x ) = p m [ E − V ( x )] , (3) m is the mass, F is the invariant flux (paper I),and primes denote spatial differentiation. Note thatsince scattering solutions are non-square-integrable, thechoice of F is arbitrary; however, throughout this pa-per we shall adopt the usual left-incident wave normal-ization convention, lim x →−∞ r sc ( x ) = 1, so that F =lim x →−∞ s ′ sc ( x ) /m . The positive and negative momen-tum functions, p sc ± ( x ) = ± s ′ sc ( x ), specify the classical tra-jectories and semiclassical Lagrangian manifolds (LMs) associated with the Ψ sc ± ( x ) solutions. Clearly, theclassical trajectories do not change direction, as there areno turning points along the real x axis. Therefore, for anarbitrary linear combination solution,Ψ = α + Ψ sc+ + α − Ψ sc − , (4)the entire reflected wave must be due to the Ψ sc − ( x ) con-tribution. The boundary conditions presume that theleft-traveling wave contribution vanish in the x → ∞ limit, implying that α − = 0 and Ψ( x ) = Ψ + ( x )—i.e.,there is no reflected wave.Various strategies have been developed to deal withthe above difficulty of the basic WKB method. The most common involve analytic continuation into thecomplex plane. Specifically, if no turning point is lo-cated along the real axis, then one is found elsewherein the complex plane. The path of the incident waveis deformed so as to give rise to a reflected wave uponencountering the complex turning point. An analysis ofthe Stokes and anti-Stokes lines that emanatefrom the complex turning point and cross the real x axis,enables one to effectively recast the above procedure interms of purely real-valued x . The net effect (as inter-preted in this paper) is the introduction of coupling fromΨ + = α + Ψ sc+ to Ψ − = α − Ψ sc − , resulting in an α − valuethat changes over x , thus yielding meaningful semiclassi-cal partial reflection probabilities.Alternatively, there is an approach due toBremmer that from the start is formulatedentirely on the real x axis. This approach is preferredfor the present purpose, as it provides a common bipolarfoundation not only for basic WKB and “sophisticated”(i.e. capable of predicting partial reflection) semiclassicalmethods, but also for exact quantum scattering applica-tions, as will be shown in Sec. II B. Bremmer’s idea is tomodel the continuous potential V ( x ) using a collectionof discontinuous steps. Solutions are determined for anarbitrary step size, which is then made infinitesimallysmall.Locally, along any given step, the two “exact”Schr¨odinger solutions (for the model step potential) areplane waves, exp( ± ip k x ), where p k = p m ( E − V k ),and V k is the (constant) potential value along the k ’thstep (Appendix A). A viable global solution, Ψ( x ), mustmatch boundary conditions appropriately at each of thestep edges. In reality, a pure positive momentum localplane wave in one step would be joined to some linearcombination of positive and negative momentum planewaves in the adjacent step, corresponding to partial re-flection and transmission off of the local step. If the re-flection contribution is ignored , e.g. so that only positive momentum local plane waves are involved, then the re-sultant approximate global solution can be easily shownto be just Ψ sc+ ( x ) in the limit of infinitesimal step size.Similarly, Ψ sc − ( x ) is obtained from the negative momen-tum local plane waves. This approximation thus leads tothe uncoupled basic WKB solutions of Eqs. (2) and (3).The above approach clearly demonstrates how “con-tinuous reflection” —i.e. that arising from continu-ous potentials—may be interpreted in more familiar dis-continuous terms. It also suggests that reflection re-quires coupling between positive and negative momentumwavefunction components . This approach has been usedto develop a sophisticated WKB approximation, simi-lar to the Stokes/anti-Stokes approach described above,wherein one-way reflection from Ψ + to Ψ − is retained,but back reflection is ignored. On the other hand, if all reflection is retained, one can derive a coupled pairof first order differential equations for Ψ ± that exactly describes the quantum scattering solution Ψ of Eq. (1).Although the Ψ ± components are not themselves sta-tionary solutions, they do correspond to the familiar in-cident/transmitted/reflected interpretations as discussedin paper II. B. Exact Quantum Dynamics Using ClassicalTrajectories
1. Time-evolution equations
In paper II, an exact quantum, CPWM bipolar de-composition scheme for stationary scattering states waspresented such that a suitable Ψ ± could be constructedfor a solution Ψ with any desired boundary conditions, FIG. 1: Time-dependent analysis of the Ψ ± coupling, ob-tained via trajectory hopping due to local reflection off ofinfinitesimal potential steps. Continuous potential reflectionis thus described in more pedagogical discontinuous terms. using a simple time-dependent approach. The classicaltrajectories and LMs were found to be identical to thoseof the basic WKB approximations, Ψ sc ± ( x )—thus satis-fying the correspondence principle, and giving rise tosmooth, well-behaved, interference-free CPWM compo-nent functions, Ψ ± . Although the method formally ap-plies only to discontinuous potentials, Sec. II A suggeststhat it ought to be extendible to continuous potentialsas well, simply by modeling the latter using infinitesi-mally small steps. This idea is indeed quite straight-forward to implement, as discussed in Appendix A andFig. 1. As in paper II, the result is a trajectory-basedscattering methodology that is both local and exact inits treatment of reflection and transmission, due to thetime-dependent nature of the approach. In contrast, theglobal character of reflection and transmission in con-ventional time-independent exact quantum methods iswell-recognized. Of the time-independent semiclassicalapproximations, the most sophisticated depend on therelative placement of local potential features such as dis-continuities and turning points, although the moreapproximate methods treat such features independently.From Appendix A, the hydrodynamic (Lagrangian)time-evolution equations for the CPWM bipolar compo-nents are found to be d Ψ ± dt = (cid:20) ∓ p ′ m + i ¯ h ( E − V ) (cid:21) Ψ ± ± p ′ m Ψ ∓ , (5)where ± p = ± p m ( E − V ) are the momenta, which de-fine the Ψ ± LM and trajectories. Note that as expected,these are identical to those of Ψ sc ± ( x ), i.e. the trajectoriesare completely classical even though the solutions are ex-actly quantum mechanical (Appendix B). All of the otherdesirable properties, as described in the preceding para-graph and in paper II, are also found to be true.
2. Comparison with unipolar Bohmian mechanics
It is worthwhile to compare Eq. (5) with theusual unipolar evolution equations of Bohmianmechanics.
Both are exact QTMs, but quantumeffects are incorporated in very different ways. Inunipolar Bohmian mechanics, the trajectory evolutionis determined by the modified potential V + Q , wherethe quantum potential correction Q is responsible forall quantum effects. Evaluation of Q requires explicitdouble spatial differentiation of the wavefunction Ψ,which in turn requires specialized numerical techniques(Sec. III B). Also, Q diverges at nodes, and is otherwisenumerically unstable in the presence of interference. Incontrast, the present bipolar QTM scheme avoids inter-ference difficulties as desired, owing to the semiclassicalcorrespondence. Quantum effects do not manifest inthe trajectories themselves (as these evolve completelyclassically), and there is no quantum potential in Eq. (5).So where do quantum effects come from in Eq. (5)?As anticipated in Sec. II A, these must be due to inter-component coupling, i.e. the last term in the equa-tion. Note that the amount of coupling is proportionalto p ′ ∝ V ′ . Thus, the coupling vanishes in the limit that V ′ → no wavefunction spatialderivatives are required in Eq. (5)—a decided advantageover the unipolar approach. In principle, every trajec-tory “splits” into a forward and backward moving tra-jectory pair at every point in space and time (paper II).The forward-moving trajectory remains on the same LMas the source, whereas the backward trajectory “hops”onto the other LM. However, unlike the situation in pa-per II, this splitting, hopping, and subsequent recombin-ing need not be considered explicitly, as it is all implicitlydealt with at the differential equation level. Nevertheless,conceptually, one may regard trajectory hopping as the source of coupling and quantum effects.Note that for scattering systems, V ′ → x . This implies that there isno asymptotic coupling (or trajectory hopping) be-tween the two Ψ ± components—an essential featurefrom the perspective of computing scattering quanti-ties such as global reflection and transmission proba-bilities (which according to normalization and bound-ary conventions discussed in Sec. II A, are obtained re-spectively via P refl = lim x →−∞ | Ψ − ( x ) | and P trans =lim x →∞ [ p ( x ) /p ( − x )] | Ψ + ( x ) | ). Equation (5) also en-sures that the asymptotic Ψ ± ( x, t ) solutions are the de-sired plane waves , and not some arbitrary linear super-position such as a sine wave. Note that in this regard,the hydrodynamic time-derivative aspect of Eq. (5) is es-sential , e.g. the ordinary Schr¨odinger equation would notpreclude asymptotic sine wave solutions.Although Ψ ± and Ψ sc ± ( x ) are identical in the asymp-totic limits (apart from the constant scaling factor α ± ),they differ in the interaction region [i.e. α ± ( x ) de-pends on x in this region], even though the trajecto-ries and LMs are the same throughout. This raisessome interesting questions vis-a-vis the interpretation ofstandard Bohmian quantities in a coupled bipolar con- text, which will be explored more fully in Sec. II B 3.Here, we consider the bipolar quantum potential, whichin the conventional sense would be obtained via q ± = − (¯ h / m )( r ′′± /r ± ) with Ψ ± = r ± exp( is ± / ¯ h ) and r ± and s ± real. But, it is not clear that such a definition shouldapply in the case where the Ψ ± are coupled and do not in-dividually evolve according to the Schr¨odinger equation.Indeed, such a definition would be inconsistent with thetime evolution of s ′± , which in any event is itself incon-sistent with the classical trajectory evolution ( s ′± = ± p ).In this context, it is perhaps more natural to define q ± such that ( V + q ± ) determines the trajectory evolution.According to this definition, q ± = 0 for the present bipo-lar CPWM formulation, even throughout the interactionregion.
3. Additional properties
As in paper II, the desired stationary state solution, asobtained from the Eq. (5) evolution equations, is not ob-served at all times, but only asymptotically in the large t limit. The same ray optics and continuous wave cavityring-down interpretations that apply in paper II also ap-ply here. Thus, at t = 0, one starts with a left-incidentasymptotic plane wave truncated outside the interactionregion.It can be shown (paper II, Appendix B) that as t → ∞ ,within any finite x interval that includes the interactionpotential, the resultant Ψ( x, t ) converges exponentiallyquickly to the correct time-dependent stationary state so-lution with appropriate x boundary conditions. A proofis provided in Appendix B, which, for comparison withthe time-dependent Schr¨odinger equation, relies on theEulerian (partial time derivative) version of Eq. (5), i.e. ∂ Ψ ± ∂t = ∓ pm Ψ ′± + (cid:20) ∓ p ′ m + i ¯ h ( E − V ) (cid:21) Ψ ± ± p ′ m Ψ ∓ . (6)Note that Eq. (6) involves a single spatial derivative ofthe wavefunction, unlike Eq. (5).Equation (6) above can be employed to derive a veryuseful flux relationship, ∂ρ ± ∂t = − j ′± ± p ′ m Re [Ψ + ∗ Ψ − ] , (7)where ρ ± = | Ψ ± | is the density, and j ± = ± ( p/m ) ρ ± isthe flux, defined in terms of the actual classical trajectoryvelocities. Taken individually, the ∂ρ ± /∂t do not obeycontinuity (except when p ′ = 0), implying that the to-tal probability for each Ψ ± component is not conserved.Together, however, they do satisfy a kind of continuityrelation, in that ∂ ( ρ + + ρ − ) /∂t = − ( j + + j − ) ′ , implyingthat the total probability for both Ψ + and Ψ − is con-served . This is an important, nontrivial result—quitedistinct from the usual probability conservation of Ψ it-self. As described in Fig. 2, Eq. (7) in effect states that ∂ρ ± /∂t is equal to the usual negative flux divergence plus FIG. 2: Schematic indicating the flux, probability, and con-tinuity relationships that exist between Ψ ± components forthe CPWM bipolar decompositions considered in this paper.All probability that leaves the Ψ − component must flow intothe Ψ + component. Consequently, R [ ρ + ( x ) + ρ − ( x )] dx isconserved for all time. a coupling term, ± ∂ρ cpl /∂t = ± ( p ′ /m )Re [Ψ + ∗ Ψ − ], rep-resenting the rate at which probability flows from oneCPWM component to the other.It should be stated that the classical trajectory CPWMbipolar decomposition scheme, in the form describedabove, is not viable for computing stationary eigenstatesbelow the potential barrier maximum. Tunneling per se is not the issue, since one can in principle apply Eq. (5)in the tunneling regime using analytic continuation (in amanner similar to that applied in paper II) as has beenconfirmed in numerical tests. However, a difficulty ariseswhen | p | → | p ′ | → ∞ , i.e. in the vicinity of the real x axis turning points. For discontinuous potentials, thisonly occurs for energies near a piece-wise barrier energy V k , manifesting as substantially increased propagationtimes (paper II). For continuous potentials, all energiesbelow the barrier height exhibit turning points along thereal x axis, and are therefore problematic. Specifically,the p ′ terms in Eqs. (5) and (6) lead to numerical insta-bilities near the turning points. C. Exact Quantum Dynamics Using ConstantVelocity Trajectories
1. Motivation
It should be noted that turning point/caustic issuessimilar to those described in Sec. II B 3 are also facedby semiclassical methods. Thus, similar remedies maypresumably be applied here (e.g. complex plane path de-formation), although these are not considered further inthis paper. Instead, we are guided by one of our origi-nal motivations, to develop methods that avoid suchsemiclassical difficulties altogether. To this end, we turnonce again to the ray optics interpretation of the bipolarapproach, as discussed in paper II.Note that in classical optics, the “ray” interpretation ofa given time-evolving electromagnetic field is not neces-sarily unique; one has a certain freedom to define rays asconvenient for a given application. Consider the caseof total internal reflection, for instance, at the bound- ary between two materials with different indices of re-fraction. Quantum mechanically, this corresponds to asingle-step discontinuous potential at an energy belowthe barrier height (paper II). According to the usual rayinterpretation, the incident rays refract to the extent thatthey become parallel to the interface, and therefore donot penetrate at all into the second medium. This pic-ture is physically somewhat incorrect however, in thatit does not capture the exponentially-damped evanes-cent wave.
Quantum mechanically, this correspondsto tunneling into the discontinuous step, which is sim-ply not described using the bipolar classical trajectorymethod as presented in Sec. II B.On the other hand, a simple, alternative ray interpre-tation can be applied to total internal reflection that does predict the evanescent wave correctly. Accordingto this interpretation, the incident ray velocities remain constant across the interface . These rays are reflectedwithin the second medium, giving rise to the evanescentwave, and also providing an explanation for the Goos-H¨anchen phenomenon.
The quantum-mechanicalanalog would be constant velocity bipolar trajectories,for which the incident wave asymptotic velocities remainconstant throughout the interaction region. Presumably,the time evolution equations derived from such a choice oftrajectories would not depend very sensitively on the en-ergy E , even in the vicinity of the barrier height, so thattunneling, and classical turning points/caustics wouldpose no special difficulties.
2. Fr¨oman and Fr¨oman methodology
The semiclassical-like CPWM discussed in Sec. II Band the Appendices does not generalize in any straight-forward manner for trajectories other than classical.However, there is an alternative time-independent for-mulation, conceptually similar to the Bremmer (B) ap-proach but differing in the details, which does allow forsuch a generalization. This approach is due to Fr¨omanand Fr¨oman (F), who use it to define generalized semi-classical approximations, although it can also be usedto derive corresponding exact quantum bipolar decom-positions. The guiding principle of the F approach (asinterpreted here) is that semiclassical solutions shouldsatisfy the invariant flux property —i.e., the left sideof Eq. (3) with s ′ sc ( x ) essentially arbitrary. If s ′ sc = p ischosen classically, then the usual basic WKB solutionsresult for Ψ sc ± ( x ). Curiously, however, the correspondingexact quantum Ψ ± ( x ) are not the same as in Sec. II B.In the F case, the decomposition of Eqs. (1) and (4) isuniquely defined via the time-independent Schr¨odingerequation and the relationΨ ′ = α + (cid:0) Ψ sc+ (cid:1) ′ + α − (cid:0) Ψ sc − (cid:1) ′ (8)= − p ′ p Ψ + ip ¯ h (Ψ + − Ψ − ) . (9)Thus, in comparing Eq. (8) to Eq. (4), the α ± coefficientsbehave in the derivative as if they were x -independentconstants, though, in fact, they are not (except whenthere is no coupling).Although not identical, the F and B bipolar decom-positions are somewhat similar in the classical trajec-tory case [compare Eq. (9) to Eq. (B2)], and both areplagued by similar numerical instabilities near turningpoints, as has been confirmed in numerical testing. How-ever, the main advantage of the F approach is that itcan be applied to arbitrary trajectories, p . In partic-ular, we choose constant velocity trajectories obtainedfrom the left-incident plane waves, i.e. p = √ mE [with lim x →−∞ V ( x ) = 0]. This gives rise to exact quan-tum CPWM components Ψ ± that are smooth and well-behaved everywhere, even near turning points and bar-rier maxima. Note that for constant velocity trajectories,Eqs. (1) and (9) imply that Ψ ± are linear combinationsof Ψ and Ψ ′ .
3. Time-evolution equations
In analogy with Sec. II B, the goal is to develop a time-dependent method to compute the F constant velocityCPWM bipolar decomposition in the large t limit. How-ever, since the F construction of the Ψ ± is radically dif-ferent from that of B, a substantially different approachthan that of Appendix A must be used. We have de-veloped several different derivations, all of which yieldthe same final results. The simplest strategy is to “workbackwards” through Appendix B, but using Eq. (9) forconstant velocity trajectories instead of Eq. (B2). How-ever, a more pedagogical approach may also be employed,as presented below.As in Eq. (8), the basic idea is to presume that thecoefficients α ± act as constants, but with respect to time derivatives rather than spatial derivatives. In particu-lar, it seems more natural here to refer to total ratherthan partial time derivatives, given the trajectory-basednature of the methodology. Less obvious at this stageis the fact that the linear combination (Ψ + − Ψ − ) mustbe used rather than Ψ = Ψ + + Ψ − itself, in order to beconsistent with the time-independent F results. We thusobtain the condition d Ψ + dt − d Ψ − dt = α + d Ψ sc+ dt − α + d Ψ sc+ dt = − iE ¯ h (Ψ + − Ψ − ) + pm Ψ ′ (10)Together with the time-dependent Schr¨odinger equation(converted to total derivative form), as applied to theEq. (1) linear combination, Eq. (10) above gives rise to aunique set of time-evolution equations for d Ψ ± /dt . In theconstant velocity trajectory case, for which Ψ sc ± ( x, t ) ∝ exp h i ( ±√ mEx − Et ) / ¯ h i , these are found to be d Ψ ± dt = i ¯ h ( E − V ) Ψ ± − i ¯ h V Ψ ∓ . (11)
4. Additional properties
Using a procedure analogous to that described at thestart of Appendix B, one can show that for station-ary states, Eq. (11) is in fact consistent with the time-independent F constant velocity CPWM bipolar decom-position. This requires the Eulerian version of Eq. (11),i.e. ∂ Ψ ± ∂t = ∓ pm Ψ ′± + i ¯ h ( E − V ) Ψ ± − i ¯ h V Ψ ∓ . (12)What is not so clear, however, is whether starting witha truncated asymptotic plane wave, one necessarily ap-proaches the stationary solution in the large t limit(Sec. II B 3, paper II). Although this has not yet beenproven, it is a reasonable assumption, given both thecounter-propagating trajectory nature of the method andflux properties similar to the B classical trajectory case(discussed below). Moreover, for all of the numerical ap-plications considered in Sec. IV, exponentially fast con-vergence is in fact observed.The F constant velocity time-evolution equations[Eq. (11)] offer some decided advantages over the classicaltrajectory approach. Since only energy quantities appearon the right hand side, there is no need to resort to an-alytic continuation in order to handle tunneling for thebelow-barrier energies. The equations may therefore beapplied with equal ease throughout the energy spectrum,and in fact, the resultant Ψ ± are qualitatively similarabove, below, and just at the potential barrier maximum(Sec. IV A). Another advantage is that the numericalpropagation does not require spatial differentiation of anykind—not even of the potential, V, to determine forcesdriving the trajectories. Furthermore, Eq. (11) may benumerically implemented as is for discontinuous poten-tials, just as easily as for continuous potentials, withoutthe need for explicit splitting and recombining of trajec-tories as in paper II (Sec. IV B).On the other hand, the F constant velocity approachintroduces a drawback that the classical trajectory meth-ods do not have to contend with when the potential is“asympotically asymmetric”—by which we mean simplythat lim x →−∞ V ( x ) = 0 = lim x →∞ V ( x ) = V , corre-sponding e.g. to an exoergic or endoergic chemical re-action. Whereas in the x → −∞ limit, the evolutionequations are uncoupled, this is not true in the x → ∞ limit if V = 0, in which case the time-dependent Ψ ± are not expected to converge. Many techniques could beused to remedy this situation, e.g. trajectories describedvia smoothly-varying sigmoid (tanh-like) functions ratherthan uniform or classical trajectories. For purposes ofthis paper, we adopt a simpler solution, wherein two F constant velocity CPWM bipolar decompositions areused, one each for the reactant and product regions.Standard boundary condition matching techniques arethen applied to join these together (Sec. III C).Regarding flux properties, Eq. (12) can be used to ob-tain ∂ρ ± ∂t = − j ′± ± V ¯ h Im [Ψ + ∗ Ψ − ] , (13)which should be compared with Eq. (7) for the B classicaltrajectory case. Note that the total combined probabil-ity for both Ψ + and Ψ − is conserved here as well, i.e. ∂ ( ρ + + ρ − ) /∂t = − ( j + + j − ) ′ , and Fig. 2 still applies.The coupling term ∂ρ cpl /∂t = (2 V / ¯ h )Im [Ψ + ∗ Ψ − ] is ofcourse different from that of Sec. II B 3, although in bothcases, the wavefunction inner product cross terms areinvolved. Several additional relations unique to the con-stant velocity case may also be easily derived from theabove equations, e.g. dρ ± /dt = ± ∂ρ cpl /∂t, (14)which states that the probability lost by a positive mo-mentum trajectory is gained by the negative momentumtrajectory directly “beneath” it (same x ), and ρ ′ + = ρ ′− for stationary solutions , (15)which states that the positive and negative componentdensity functions are identical apart from a constant.Note that for asymptotically symmetric potentials, andthe normalization and boundary conditions described inSecs. II A and II B 2, this constant must be the transmis-sion probability itself, P trans . III. NUMERICAL DETAILS
In this section, we discuss the numerical details associ-ated with implementing the bipolar time-evolution equa-tions of Sec. II as a practical and robust algorithm forcomputing stationary scattering states. Although anyboundary conditions may be considered, our emphasisshall be on left-incident solutions, for which Ψ − → x → ∞ . The region of interest is defined via x L ≤ x ≤ x R , a region presumed to include the entire potential in-teraction to the desired level of accuracy. According tothe time-dependent ray optics picture developed in pa-per II, at the initial time, Ψ = Ψ + is a left-incident planewave truncated on the right at x = x L . This initial wavepropagates into the interaction region, wherein it is cou-pled to Ψ − and eventually reaches a stationary state.Outside of the interaction region, the Ψ ± coupling maybe regarded as effectively zero. Thus, one may computereflection and transmission probabilities, P refl and P trans ,simply by evaluating Ψ + ( x R ) and Ψ − ( x L ) at sufficientlylarge times.From a purely numerical point of view, theabove scheme offers many advantages—e.g., the abil-ity to compute specific scattering states without theneed for complex scaling or complex absorbingpotentials that would increase the coordi-nate range unnecessarily. Moreover, the density of gridpoints may be substantially reduced, owing to the fact that the interference-free component functions Ψ ± aregenerally smooth and slowly varying as compared to Ψitself. For the same reason, a larger time-step is also an-ticipated. Most importantly however, since q ± = 0, thereis no need to compute on-the-fly numerical spatial deriva-tives. With regard to parallel computation, therefore,there is no need for trajectories to communicate within aLM, although the coupling requires position-dependentcommunication between positive and negative LM tra-jectories. A. Algorithmic Issues
In order to solve the time-evolution equations numer-ically, two trajectory grids are created, one for eachbipolar component of the total wavefunction. Hereafter,“upper” will be used to describe attributes of the Ψ + component, and “lower” for the Ψ − component (e.g.“lower/upper grid,” etc.). On each grid, the correspond-ing wavefunction component is spatially discretized atthe grid point locations and evolved in time.For the applications considered here, we have foundit convenient to modify the ray optics picture somewhatfrom the form described above and in paper II. First, wedefine an initial negative LM grid, even though Ψ − itselfis zero initially. The idea here is that unlike paper II,we wish to avoid explicit trajectory hopping. Conse-quently, both sets of grid points exist for all time andevolve independently of each other (though the compo-nent wavefunctions do interact). Second, to avoid un-necessary computational effort, no propagation is doneoutside of the interaction region of interest. Instead, atuniform time intervals, new Ψ + trajectories are fed in at x = x L , with an initial Ψ + value consistent with the pos-itive momentum asymptotic plane wave solution. Later,as these upper grid trajectories reach x = x R , they arediscarded. The Ψ − situation is similar, except that theinitial Ψ − ( x R ) value is set to zero, and the trajectoriesare discarded as they reach x = x L .The most substantial difference from the ray optics pic-ture, however, is found in our implementation of the con-stant velocity method, for which the grids are distributeduniformly throughout the interaction region at all times .At the initial time, the upper and lower grids coincide.The initial Ψ + is taken to be the asymptotic plane wavesolution throughout the interaction region, and the initialΨ − is set to zero. The above modification, which we callthe “non-wavefront” approach, is certainly not necessary,and is introduced simply for convenience. Calculationsperformed both ways reveal that both converge exponen-tially to the correct stationary solution. However, theinitial convergence of the non-wavefront calculations isfaster, owing to the fact that the coupling takes effect im-mediately. The numerical algorithm is also easier to im-plement. Consequently, all constant velocity results pre-sented in Sec. IV were obtained using the non-wavefrontcode.As per Sec. II, the upper and lower grids move clas-sically or with constant velocity, as appropriate, andthe Ψ ± contributions evolve in accord with Eq. (5) orEq. (11). Since the grids move in opposite directions,they do not remain commensurate over time. Numericalinterpolation (Sec. III B) is therefore required to computethe coupling contribution from the component wavefunc-tion of one grid to the other. As an alternative to theabove trajectory-based (Lagrangian) approach, we havealso implemented a fixed-grid (Eulerian) bipolar propa-gation algorithm for the constant velocity case. The pri-mary advantage of fixed-grid propagation is that the twogrids remain commensurate for all time, thus avoidingthe need for interpolation. However, this is achieved atthe expense of switching to the Eulerian evolution equa-tions [Eq. (12)], which require explicit numerical spatialdifferentation of the bipolar wavefunction components. B. Low Level Issues
For the numerical propagation of the time-evolutionequations, two fourth-order explicit integrators wereconsidered—the multi-step Adam’s/Bashforth, andRunge-Kutta methods. Although, both integrators arefourth-order accurate, and require approximately thesame CPU time, multi-step integrators can not be triv-ially implemented for the asymptotically asymmetric po-tential case (Sec. III C). Consequently, the basic fourth-order Runge-Kutta method was used for all results pre-sented herein. In the future, the time integration canbe made more accurate by using a time-adaptive Runge-Kutta approach such as the Cash-Carp method. Formost of the calculations performed here, a time-step of∆ = 10 a.u. was found to be sufficient to achieve a com-puted accuracy of 10 − or better (Sec. IV). Generallyspeaking however, the time-step issue is quite problem-dependent, as it can depend on both grid point velocitiesand spacing.As discussed in Sec. III A, the trajectory-based algo-rithms require interpolation of both bipolar componentsin order to compute coupling contributions. For all exam-ples considered here, these interpolates were calculatedvia a polynomial moving least squares (MLS) method. In MLS methods, local low-order polynomial fits are cal-culated about each grid point via a small stencil of sur-rounding neighbor points. To ensure that the fit willexactly represent the function values at the grid pointlocations, the stencil size and polynomial order must beset equal, thus effectively transforming the MLS methodinto a moving-interpolation method. For all calculationsperformed here, five stencil points were used, correspond-ing to a symmetric stencil (for the interior grid points)and quartic polynomial interpolation.For the fixed-grid algorithms, explicit numerical spatialdifferentiation of both CPWM bipolar components mustbe performed. For the present study, these were calcu-lated using centered, fourth-order finite differences for the interior grid points, and one-sided fourth-order finitedifference for the boundary grid points at x L and x R . C. Asymptotically Asymmetric Potentials
As discussed in Sec. II C 4, the constant velocity time-evolution equations will exhibit coupling in the x → ∞ asymptote if V ( x ) is asymptotically asymmetric, i.e. iflim x →∞ V ( x ) = V = 0. To remedy this, we con-struct two sets of solutions, one for the “reactant” region, x ≤ x , and another for the “product” region, x ≥ x ,with x the dividing point. For the former solutions (de-noted via the “L” subscript), Eq. (11) may be used di-rectly, with p L = √ mE . For the product solutions,however (denoted “R”), the evolution equations must bemodified slightly to accommodate the generalized asymp-totic potential condition, d Ψ R ± dt = i ¯ h ( E − V − V ) Ψ R ± − i ¯ h ( V − V ) Ψ R ∓ , (16)since the positive trajectory momentum is now p R = p m ( E − V ). It is clear from Eq. (16) that the cou-pling vanishes as x → ∞ .For all propagation times, the wavefunction and itsfirst derivative must be continuous across x . Remark-ably, we can treat this boundary condition like an elemen-tary plane wave scattering off a step potential (paper II).This is achieved via Eq. (9), in terms of which the twoexact conditions becomeΨ L + + Ψ L − = Ψ R + + Ψ R − p L (cid:2) Ψ L + − Ψ L − (cid:3) = p R (cid:2) Ψ R + − Ψ R − (cid:3) , (17)where Ψ L + = Ψ L + ( x ), etc. In the trajectory-based al-gorithm, the left and right incident wavefunction values,Ψ L + and Ψ R − , are known at any given time, whereas thelocally transmitted values Ψ R + and Ψ L − are unknown.Equation (17) enables one to solve for the two unknowns,and thus to continue propagating trajectories through thedividing point, x .The solutions areΨ R + = (cid:18) p L p L + p R (cid:19) Ψ L + + (cid:18) p R − p L p L + p R (cid:19) Ψ R − Ψ L − = (cid:18) p R p L + p R (cid:19) Ψ R − − (cid:18) p R − p L p L + p R (cid:19) Ψ L + , (18)which do indeed correspond to transmission and reflec-tion coefficients for stationary states of the discontinu-ous step potential (Appendix A and paper II). Numeri-cally, the propagation is implemented as follows. Equa-tions (11) and (16) are used until a trajectory reaches x ,at which point it is replaced with a new trajectory on theother side of x via Eq. (18). This requires a “known”value from the opposite LM. If the grids are not com-mensurate (i.e. trajectories from opposite LMs do notarrive at x at the same time), it is necessary to applyextrapolation to the opposite manifold to determine the“known” value. FIG. 3: Converged positive component density ρ + ( x ) (solidline) and negative component density ρ − ( x ) (dashed line) forthe E = 450 cm − stationary state of the V = 400 cm − Eckart barrier system, as obtained using the classical trajec-tory CPWM bipolar decomposition. Shaded area indicatesthe Eckart potential V ( x ). IV. RESULTS
In this section, four different stationary scattering ap-plications are considered, all with comparable charac-teristic parameters—e.g., the same proton-like mass of m = 2000 a.u., barrier height V = 400 cm − ≈ . × smaller than in paper II), and interac-tion region, [ x L = − , x R = 3 a.u.]. All four systemswere solved using the non-wavefront, F, constant velocitytrajectory-based method, hereinafter referred to as the“constant velocity trajectory” scheme. For the Eckartsystem (Sec. IV A), additional calculations were also per-formed using the constant velocity fixed-grid scheme, andthe wavefront, B, classical trajectory method, now re-ferred to as the “classical trajectory” scheme. A. The Eckart Barrier
The first application considered is the symmetricEckart problem, defined via V ( x ) = V sech( αx ) , (19)with parameter values V = 400 cm − , and α = 3 . P refl and P trans values were computed for various energies E > V , to a convergence error of10 − , and in every case found to match the exact ana-lytical results to within the same error. A density plotof the CPWM bipolar components for a typical solution( E = 450 cm − ≈ .
002 a.u.) is presented in Fig. 3. Formost energies E , these calculations are roughly as effi-cient as the constant velocity calculations described be-low. However, the required propagation time does indeedincrease rapidly as E → V , as expected. In this limit,the classical trajectory bipolar decomposition also be-comes unstable; the small central peaks evident in Fig. 3become increasingly pronounced, eventually developinginto singularities.The second study is a detailed convergence and ef-ficiency comparison between the trajectory and fixed-grid implementations of the constant velocity method(Sec. III A). This study also serves as a benchmark forestablishing the numerical parameter values needed toachieve a 10 − accuracy level. Both schemes were ap-plied to a calculation of the E = V stationary state, i.e.the classical “worst-case scenario,” for a varying numberof grid points, N . A time-step of ∆ = 1 a.u. was used forthe trajectory calculation, and ∆ = 0 . t max = 10 ,
000 a.u. These parameters were sufficientlyconverged as to ensure that their contribution to the to-tal numerical error is insignificant (i.e. 10 − or less). InFig. 4, the resultant CPU times on a 2.60 GHz Pentiumcomputer are presented. Note that despite the ten-foldincrease in the number of time-steps for the fixed-gridcase, the total CPU time is still markedly faster for allgrid sizes considered. This is due both to the searchoperation (to find the nearest opposite LM trajectories)and the interpolation matrix inversions required at eachtime-step by the trajectory scheme.In Fig. 5, the log of the computed P refl and P trans er-rors (relative to known analytical values ) are plottedversus N . Across N , the trajectory-based results areseen to be much more accurate than the fixed-grid re-sults, although the latter are substantially improved byincreasing the density of grid points. For example, inorder to achieve 10 − accuracy for both P refl and P trans ,only N = 25 trajectories were needed, whereas N = 91fixed-grid points were required. This is consistent withthe fixed-grid method requiring explicit spatial differen-tiation, introducing an additional source of numerical er-ror. Note that for this comparison, the CPU times arecomparable—4.3 seconds and 2.7 seconds, respectively.The fixed-grid calculation is slightly faster, although thissituation would be reversed for more accurate calcula-tions (e.g. 10 − ) and/or higher dimensionalities. Theperformance of both methods is in any event remark-able, with the large- N trajectory case representative ofthe most accurate trajectory-based quantum scatteringcalculations performed to date. Further improvementsare also possible, however, as discussed in Sec. III B.The final study was a constant velocity trajectorycalculation of P refl and P trans over a wide range of E FIG. 4: CPU time vs. no. of grid points, N , required to prop-agate trajectory (empty circles) and fixed-grid (filled circles)implementations of the constant velocity CPWM bipolar de-composition to a final system time of t max = 10 ,
000 a.u. Thetrajectory and fixed grid time-steps used were ∆ = 1 a.u. and∆ = 0 . P refl (solidlines) and transmission probabilities P trans (dashed lines) vs.no. of grid points N for the E = 400 cm − stationary state ofthe V = 400 cm − Eckart barrier system, as obtained usingtrajectory and fixed-grid implementations of the constant ve-locity CPWM bipolar decomposition. Errors are relative toknown analytical results (Ref. [47]). values—including those above, below, and at the bar-rier maximum, E = V . The parameters ∆ = 10 a.u., t max = 10 ,
000 a.u., and N = 31 were chosen to cor-respond to a target accuracy of 10 − . The resultingstationary state solution for the E = V case is displayedin Figs. 6 and 7. Note the contrast between the smooth,slowly-varying CPWM bipolar densities ρ ± in Fig. 6 vs.the oscillatory total density ρ of Fig. 7. Nevertheless,as indicated in the latter figure, the numerically recon-structed total density agrees with the analytical result FIG. 6: Converged positive component density ρ + ( x ) (solidline) and negative component density ρ − ( x ) (dashed line) forthe E = 400 cm − stationary state of the V = 400 cm − Eckart barrier system, as obtained using the constant veloc-ity trajectory CPWM bipolar decomposition. Shaded areaindicates the Eckart potential V ( x ).FIG. 7: Total wavefunction density ρ ( x ) for the E = 400 cm − stationary state of the V = 400 cm − Eckart barrier system,as obtained using the constant velocity trajectory CPWMbipolar decomposition (filled circles), and compared with ex-act analytical results (solid line). Shaded area indicates theEckart potential V ( x ). to within the target accuracy of 10 − at all grid points.Note also that ρ + and ρ − are identical apart from a con-stant, as predicted in Sec. II C 4.For all 26 energy values considered, the resulting con-stant velocity bipolar decompositions are qualitativelysimilar to those presented here for E = V . Indeed, inno respect do these calculations seem to depend sensi-tively on the value of E , unlike the classical trajectorycase. Consequently, the same numerical parameters asdescribed above were used for all energies. The resulting1 FIG. 8: Computed reflection probabilities P refl (filled circles)and transmission probabilities P trans (empty circles) vs energy E for the V = 400 cm − Eckart barrier system, as obtainedusing the constant velocity trajectory CPWM bipolar decom-position, and compared with exact analytical results (solidline). computed P refl and P trans values are presented in Fig. 8and compared with analytical values. Once again, allerrors are found to be smaller than 10 − . B. The Square Barrier
Although the primary focus of this paper is continu-ous potentials, it is worth emphasizing that the methodspresented here may be applied to discontinuous poten-tials as well. As a case in point, we reexamine the squarebarrier system introduced in paper II. Note that the clas-sical trajectory scheme for
E > V would yield resultsidentical to those of paper II, as is easily verified fromEq. (5). Instead, we focus on the constant velocity tra-jectory scheme, which can be applied directly withoutany algorithmic modification.For the following square barrier study, we used abarrier height of V = 400 cm − , and barrier edges of x = − x = 1, respectively. As in Sec. IV A, N = 31 initial grid points were distributed uniformlythroughout the interaction region, and the time-step wasdefined as ∆ = 10 a.u. The energy was taken to be E = 450 cm − ≈ .
002 a.u., i.e. slightly above-barrier.The propagation was continued until the computed P refl and P trans values were both converged to less than 10 − ,which was found to require approximately 1000 time-steps.Figure 9 is a density plot of the converged CPWMbipolar component solutions. As in the Eckart case, theresults are well-behaved and interference-free through-out, although there is a discontinuity in the first spatialderivative at the step edges (but not for the total Ψ it- FIG. 9: Converged positive component density ρ + ( x ) (solidline) and negative component density ρ − ( x ) (dashed line) forthe E = 450 cm − stationary state of the V = 400 cm − square barrier system, as obtained using the constant veloc-ity trajectory CPWM bipolar decomposition. Shaded areaindicates the square barrier potential V ( x )FIG. 10: Computed reflection probabilities P refl (filled circles)and transmission probabilities P trans (empty circles) vs energy E for the V = 400 cm − square barrier system, as obtainedusing the constant velocity trajectory CPWM bipolar decom-position, and compared with exact analytical results (solidline). self). For the total density ρ , the computed values andthe analytical results are once again in agreement to 10 − or better, at all grid points.The calculation described above was repeated for a to-tal of 45 different energy values. The resultant converged P refl and P trans values are presented in Fig. 10, as are theexact analytical results. Relative to the latter, all com-puted probability errors are found to be less than 10 − .2 FIG. 11: Converged positive component density ρ + ( x ) (solidline) and negative component density ρ − ( x ) (dashed line) forthe E = 500 cm − stationary state of the V = 400 cm − up-hill ramp system, as obtained using the constant velocity tra-jectory CPWM bipolar decomposition, modified for asymp-totically asymmetric potentials. The dividing point is x = 0.Shaded area indicates the uphill ramp potential V ( x ). C. The Uphill Ramp
The next application considered is an asymptoticallyasymmetric system—the continuous “uphill ramp” de-fined via V ( x ) = V h (cid:16) x α (cid:17)i , (20)with the parameters V = 400 cm − and α = 0 . Note that lim x →−∞ V ( x ) = 0 andlim x →∞ V ( x ) = V , as presumed in Sec. III C.The propagation was performed using the constant ve-locity trajectory scheme described in Sec. III C, with thedividing point chosen to be x = 0. All other computa-tional parameters were identical to those of the previousexamples, except for the energy, which was chosen to be E = 500 cm − ≈ . P refl and P trans were converged toless than 10 − . Figure 11 is a density plot of the re-sultant CPWM bipolar component solutions. There is adiscontinuity at x = 0, which is to be expected giventhe “step-like” nature of the join (the same behavior isobserved in paper II). Importantly, however, this discon-tinuity does not give rise to any numerical problems, sincespatial derivatives are not required in the trajectory inte-gration (though one must be a bit careful with the inter-polation). Away from x , both density plots are smoothand well-behaved.Despite the discontinuities in ρ ± , the computed ρ andits first derivative are continuous, owing to the Eq. (18)relations. This is evident in Fig. 12, a density plot for FIG. 12: Total wavefunction density ρ ( x ) for the E =500 cm − stationary state of the V = 400 cm − uphill rampsystem, as obtained using the constant velocity trajectoryCPWM bipolar decomposition (filled circles), and comparedwith exact analytical results (solid line, Ref. [48]). Shadedarea indicates the uphill ramp potential V ( x ). the converged total Ψ. As is clear from the plot, no visi-ble discrepancies may be observed between the computedand analytically exact solutions. D. The Double-Gaussian Barrier
The previous three example systems serve as usefulbenchmarks for testing and evaluating the new bipolarmethodologies under a wide range of conditions. In par-ticular, all three have known analytical solutions, againstwhich the computed results may be compared. However,we feel it is also worthwhile to consider at least one sys-tem which has not previously been solved. One suchexample, which also presents a qualitatively new varietyof problem, is the symmetric double-Gaussian barrier po-tential, V ( x ) = V n exp h − β ( x − x ) i + exp h − β ( x + x ) io . (21)For the results presented here, the following parameterswere used: V = 400 cm − ; β = 9 a.u.; x = 0 .
75 a.u.Once again, the constant velocity trajectory methodwas employed, with the same computational parametersas described previously. The energy was chosen to bejust at the barrier height, E = V . Figure 13 displaysthe resulting converged CPWM bipolar densities. As inthe previous examples, these are identical apart from aconstant, and are otherwise smooth and slowly-varying.Note that the double-barrier nature of the potential givesrise to interesting features in the bipolar densities notpreviously observed. In particular, despite the fact that V ( x ) is changing in the central well region between the3 FIG. 13: Converged positive component density ρ + ( x ) (solidline) and negative component density ρ − ( x ) (dashed line)for the E = 400 cm − stationary state of the V =400 cm − double-Gaussian barrier system, as obtained usingthe constant velocity trajectory CPWM bipolar decomposi-tion. Shaded area indicates the double-Gaussian potential, V ( x ). two barriers, the bipolar densities are nearly flat acrossthis region, resulting in a well-defined “reaction interme-diate” state between reactants and products. Moreover,the density plots enable one to assign quantitative prob-ability values for the intermediate state.In contrast, the above interpretation and quantitativeassignments would be very difficult, if not impossible,to glean directly from Ψ( x ) itself. This is evident fromFig. 14, a density plot of the total Ψ for the above double-Gaussian barrier calculation, obtained via interpolationand superposition of the two converged bipolar compo-nent solutions, Ψ ± . Note the interference present bothin the reactant region (due to reflected trajectories) and the intermediate region. Thus, not only the intermedi-ate probabilities, but also the reflection probability, aredifficult to read directly from such a plot. V. SUMMARY AND CONCLUSIONS
Scattering applications are of paramount importancefor chemical reactions, because all reactions may be re-garded as scattering events. From a theoretical exactquantum perspective, therefore, multichannel scatteringtheory, both time-dependent and time-independent,will always play an essential role. At the same timehowever, trajectory-based methods also bring much tobear on dynamics, providing great insight into reactiveprocesses, vis-a-vis the determination of which trajecto-ries make it past the barrier to products vs. those thatdo not. Quantum trajectory methods (QTMs) thereforeexhibit great potential promise as a chemical dynamicstool, combining a trajectory-based description with ex- FIG. 14: Total wavefunction density ρ ( x ) (solid line) for the E = 400 cm − stationary state of the V = 400 cm − double-Gaussian barrier system, as obtained using the constant ve-locity trajectory CPWM bipolar decomposition. Shaded areaindicates the double-Gaussian potential, V ( x ). act quantum dynamics. However, all reactive systemsexhibit interference between the incident and reflected(non-reactive) waves, thus causing numerical instabilityproblems for conventional unipolar QTMs. CPWM bipo-lar decompositions offer a natural means of alleviatingthis interference difficulty, by splitting the reactant re-gion Ψ into incident Ψ + and reflected Ψ − components,neither of which exhibits interference on its own. More-over, this splitting can be extended through the inter-action region over to the product region, by which pointΨ + has transformed smoothly into the transmitted wave,and Ψ − has damped to zero.As discussed in Sec. II and the Appendices, the CPWMapproach borrows conceptually from semiclassical scat-tering methods. Indeed, for the first bipolar decompo-sition considered (B, Sec. II B) the bipolar trajectoriesare simply equal to the classical trajectories themselves,and the correspondence principle is satisfied in preciselythe usual WKB limit, V ′ →
0. Three features, however,contribute to render the present approach fundamentallydifferent from basic WKB theory: (1) time-dependentformulation; (2) coupling between Ψ + and Ψ − ; (3) uni-versal, local reflection and transmission formulae (see alsopaper II). Point (3) is what determines point (2), i.e. ifonly local transmission were considered without reflec-tion, then the coupling would vanish and the basic WKBsolutions would result. The combination of (1) and (3)provides a local physical understanding related to the rayoptics picture in electromagnetic theory, and also givesrise to useful flux relations and numerical algorithms. Re-garding the second bipolar decomposition considered, i.e.the F, constant velocity scheme, this was motivated bypractical concerns, but also by an alternative ray opticsdescription (Sec. II C). This can be related to the semi-4classical approach of Fr¨oman and Fr¨oman, and in thatcontext, also satisfies a generalized kind of correspon-dence principle.Note that neither set of evolution equations [Eqs. (5)and (11)] involves a quantum potential; instead, all quan-tum effects manifest through Ψ ± coupling. In both cases,the trajectories themselves are not “context-sensitive,” in that they may be computed independently of the Ψ ± evolution. Moreover, no spatial differentiation of thewavefunction is required, although there may be situa-tions where explicit calculation of one spatial derivativeis numerically advantageous (Sec. IV A). Several othernumerical modifications have also been introduced forconvenience (Sec. III A), resulting in a shift from an ex-act time-dependent Schr¨odinger interpretation of the dy-namics [wherein the exact stationary state is “revealed”over time (paper II)] to what may be regarded as moreof a relaxation approach. Be that as it may, the result-ing algorithms offer a remarkably simple, efficient, andaccurate means of performing reactive scattering calcu-lations of all kinds in 1D (Sec. IV). Indeed, the time-stepsrequired for the benchmark molecule-like systems consid-ered here are orders of magnitude larger than for typicalfixed-grid calculations performed at a comparable level ofaccuracy. Moreover, the converged bipolar solution den-sity plots render the determination of global reflectionand transmission probabilities, as well as probabilitiesfor reaction intermediate states, quite straightforward.In future publications we will continue to generalizethe methodology described here and in paper II, for thetype of multidimensional time-dependent wavepacket dy-namics relevant to real chemical physics applications. Asadditional motivation for the present work, we now sketchhow this might be achieved. First, it is necessary to gen-eralize the stationary state results of this paper and pa-per II for arbitrary time-evolving wavepackets. This isdone initially for the discontinuous step potential, andthen generalized for arbitrary continuous potentials ina manner similar to Appendix A. In fact, much of thegroundwork is already laid, in that the time-dependentframework has already been introduced.The generalization to multidimensional systems is lessstraightforward but can certainly be achieved (such cal-culations have already been performed, as will be re-ported in a future publication). Conceptually at least,many direct chemical reactions can be described usinga single scattering reaction coordinate, plus additional“bound” coordinates. It is natural to consider applyingthe current CPWM bifurcation to the former, and thepaper I bifurcation to each of the latter. However, thetotal number of wavefunction components would then be2 D where D is the number of degrees of freedom. On theother hand, for most time-dependent wavepacket calcu-lations, node formation in Ψ is associated primarily withthe reaction coordinate itself, due to wavepacket reflec-tion off of the reaction profile barrier. Thus, a naturalapproach would be to bifurcate only along the reactioncoordinate. Only two component wavefunctions result, regardless of D . It remains to be seen whether such aprocedure will render QTM calculations possible for ac-tual molecular systems. Nevertheless, it seems very likelythat some such bipolar or multipolar approach will go along way towards ameliorating the infamous node prob-lem, which has thus far severely limited the effectivenessof QTMs in the molecular arena. Acknowledgments
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. Ja-son McAfee is also acknowledged for his aid in convertingthis manuscript to an electronic format suitable for thearXiv preprint server.
APPENDIX A: DERIVATION OF BIPOLARTIME-EVOLUTION EQUATIONS
Following the notation of paper II, Sec. II D 3, wepresume a discontinuous potential, for which x k ( k =1 , , . . . , l ) denote the locations of the l discontinuities.The x k divide configuration space into l + 1 regions, la-beled 0 ≤ k ≤ l . In each region k , the potential has adifferent constant value, V k . The discontinuous systemdescribed above may be used to model any continuous potential system, V ( x ), by defining V k = V ( x mk ) [where x mk = ( x k + x k +1 ) / x k +1 − x k ) → k .As the derivation is a time-dependent one, it is conve-nient to introduce a small (ultimately differential) timeincrement ∆, which is then used to determine the x k as follows. Associated with each region k is the posi-tive classical momentum value, p k = p m [ E − V k ] [i.e. p sc+ ( x mk )]. The x k are chosen such that a particle movingwith momentum p k would traverse the region k in time∆, i.e. ( x k +1 − x k ) = ∆( p k /m ). Consider a trajectorywhich at t = 0, is located at the k ’th region midpoint, x mk ,heading to the right with momentum p k . This is clearlya positive LM trajectory, carrying a contribution of thecomponent wavefunction Ψ + k = Ψ + ( x = x mk , t = 0)(Fig. 1). It is presumed that there are also negative LMtrajectories along which Ψ − is propagated, but for nowthe emphasis is on Ψ + .From t = 0 to t = ∆ /
2, the positive LM trajectorytravels from x = x mk to the discontinuity at x = x k +1 .As per paper II, the propagation is that of a plane wave,i.e. Ψ + k → Ψ + ( x k +1 , ∆ / + k exp (cid:18) i ∆2¯ h (cid:20) p k m − V k (cid:21)(cid:19)
5= Ψ + k exp (cid:18) i ∆2¯ h [ E − V k ] (cid:19) . (A1)At this point, the trajectory splits into two: one thatcontinues in the forward direction along the positive LM,transmitting into region k + 1 with momentum p k +1 ; theother reflected backwards along the negative LM withmomentum − p k (Fig. 1). According to paper II Eqs. (17)and (18), the trajectory bifurcation introduces a factor of2 p k / ( p k + p k +1 ) into the the transmitted Ψ + wave, and( p k − p k +1 ) / ( p k + p k +1 ) into the reflected Ψ − wave.During the remaining time evolution from t = ∆ / t = ∆, the Ψ + trajectory moves to x mk +1 [the midpoint ofthe adjacent ( k + 1)’th region], resulting in an additionalphase factor analogous to that of Eq. (A1). The finalresult isΨ + k → Ψ + ( x mk +1 , ∆) (A2)= Ψ + k exp (cid:18) i ∆¯ h [ E − V k − V k +1 ] (cid:19) (cid:18) p k p k + p k +1 (cid:19) . The reflected trajectory, meanwhile, moves back to theoriginal location at x mk , resulting in the following for Ψ − :Ψ + k → Ψ − ( x mk , ∆) (A3)= Ψ + k exp (cid:18) i ∆¯ h [ E − V k ] (cid:19) (cid:18) p k − p k +1 p k + p k +1 (cid:19) . We thus find that Ψ + k at time t = 0 contributesto both Ψ +( k +1) and Ψ − k at time t = ∆. However,these must be combined with similar contributions fromthe initial Ψ − ( k +1) in order to determine the total finalΨ +( k +1) and Ψ − k (Fig. 1). Following an analysis similarto the above, the negative LM contributions are easilyshown to be the following:Ψ − ( k +1) → Ψ + ( x mk +1 , ∆) (A4)= − Ψ − ( k +1) exp (cid:18) i ∆¯ h [ E − V k +1 ] (cid:19) (cid:18) p k − p k +1 p k + p k +1 (cid:19) . Ψ − ( k +1) → Ψ − ( x mk , ∆) (A5)= Ψ − ( k +1) exp (cid:18) i ∆¯ h [ E − V k − V k +1 ] (cid:19) (cid:18) p k +1 p k + p k +1 (cid:19) . By adding Eqs. (A3) and (A4), we obtain the final ex-pression for Ψ + ( x mk +1 , ∆). Subtracting Ψ + k , dividing by∆, and taking the limit ∆ → d Ψ + /dt . We thus obtainlim ∆ → (cid:20) Ψ + ( x mk +1 , ∆) − Ψ + k ∆ (cid:21) =lim ∆ → ((cid:20) (cid:18) p k − p k +1 p k + p k +1 (cid:19) + i ¯ h ( E − V k − V k +1 ) × (cid:18) p k p k + p k +1 (cid:19) Ψ + k − (cid:20)
1∆ + i ¯ h ( E − V k +1 ) (cid:21) × (cid:18) p k − p k +1 p k + p k +1 (cid:19) Ψ − ( k +1) ) . (A6)Recall that ∆ = ( x k +1 − x k ) m/p k . In the small ∆ limit, V k +1 → V k , p k +1 → p k , and − ( p k − p k +1 ) / ( x k +1 − x k ) → p ′ ( x ) = ∂p ( x ) /∂x . Substituting into Eq. (A6), and re-placing k subscripts with functions of x = x k yields: d Ψ + dt = (cid:20) − p ′ m + i ¯ h ( E − V ) (cid:21) Ψ + + p ′ m Ψ − (A7)A similar analysis applied to Eqs. (A4) and (A5) yields: d Ψ − dt = (cid:20) p ′ m + i ¯ h ( E − V ) (cid:21) Ψ − − p ′ m Ψ + (A8) APPENDIX B: PROOF THAT BIPOLARSTATIONARY SOLUTIONS SATISFYTIME-INDEPENDENT SCHR ¨ODINGEREQUATION
Consider the stationary solutions of Eq. (6), i.e. ∂ Ψ ± /∂t = − ( i/ ¯ h ) E Ψ ± . Substituting in these expres-sions for the time derivatives, and rewriting to obtainexpressions for the spatial derivatives, yields:Ψ ′± = (cid:18) − p ′ p ± ip ¯ h (cid:19) Ψ ± + p ′ p Ψ ∓ (B1)These equations are identical to those of the bipolar time-independent stationary state decomposition described inRef. 25. Adding the Ψ + ′ and Ψ −′ equations togetherresults in Ψ ′ = ip ¯ h (Ψ + − Ψ − ) . (B2)Applying spatial differentiation and substitutingEq. (B1) into the resulting right hand side yieldsΨ ′′ = − ( p / ¯ h )Ψ. Substitution into the Schr¨odingerequation then results in − ¯ h m Ψ ′′ + V Ψ = p m Ψ + V Ψ = E Ψ . (B3)The stationary solutions of Eq. (6) are therefore consis-tent with the time-independent Schr¨odinger equation. ∗ Electronic address: [email protected] B. Poirier, J. Chem. Phys. , 4501 (2004). B. Poirier, J. Chem. Phys. , (in press). D. Babyuk and R. E. Wyatt, J. Chem. Phys. , 9230(2004). R. E. Wyatt,
Quantum Dynamics with Trajectories: Intro-duction to Quantum Hydrodynamics (Springer, New York,2005). 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). D. V. Shalashilin and M. S. Child, J. Chem. Phys. ,10028 (2000). R. E. Wyatt and E. R. Bittner, J. Chem. Phys. , 8898(2001). R. E. Wyatt and K. Na, Phys. Rev. E , 016702 (2001). I. Burghardt and L. S. Cederbaum, J. Chem. Phys. ,10312 (2001). E. R. Bittner, J. B. Maddox, and I. Burghardt, Int. J.Quantum Chem. , 313 (2002). K. H. Hughes and R. E. Wyatt, Phys. Chem. Chem. Phys. , 3905 (2003). E. Madelung, Z. Phys. , 322 (1926). J. H. van Vleck, Proc. Natl. Acad. Sci. U.S.A. , 178(1928). D. Bohm, Phys. Rev. , 166 (1952). D. Bohm, Phys. Rev. , 180 (1952). T. Takabayasi, Prog. Theor. Phys. , 341 (1954). P. R. Holland,
The Quantum Theory of Motion (Cam-bridge University Press, Cambridge, 1993). E. R. Floyd, Physics Essays , 135 (1994). M. R. Brown, arXiv:quant-ph/0102102 (2002). J. R. Taylor,
Scattering Theory (John Wiley & Sons, Inc.,New York, NY, 1972). J. Heading,
An Introduction to Phase-integral Methods (Methuen, London, 1962). N. Fr¨oman and P. O. Fr¨oman,
JWKB Approximation (North-Holland, Amsterdam, 1965). M. V. Berry and K. V. Mount, Rep. Prog. Phys. , 315(1972). 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). M. S. Child,
Molecular Collision Theory (Dover, New York,1996). B. Poirier and T. Carrington, Jr., J. Chem. Phys. , 77(2003). H. Bremmer, Commun. Pure Appl. Math. , 105 (1951). B. Poirier and T. Carrington, Jr., J. Chem. Phys. , 17(2003). J. D. Jackson,
Classical Electrodynamics , 2nd ed. (JohnWiley & Sons, New York, 1975). L. Brillouin, Ann. Phys. , 177 (1914). J. O. Hirschfelder, A. C. Christoph, and W. E. Palke, J.Chem. Phys. , 5435 (1974). see for example, Int. J. Quant. Chem. , (1978), specialissue. W. P. Reinhardt, Ann. Rev. Phys. Chem. , 223 (1982). V. Ryaboy, N. Moiseyev, V. A. Mandelshtam, and H. S.Taylor, J. Chem. Phys. , 5677 (1994). G. Jolicard and E. J. Austin, Chem. Phys. Lett. , 106(1985). T. Seideman and W. H. Miller, J. Chem. Phys. , 4412(1992). U. V. Riss and H.-D. Meyer, J. Phys. B: At. Mol. Phys. , 4503 (1993). J. G. Muga, J. P. Palao, B. Navarro, and I. L. Egusquiza,Phys. Rep. , 357 (2004). G. M. Phillips and P. J. Taylor,
Theory and Applicationsof Numerical Analysis (Academic Press, New York, 1996). W. H. Press et al , in
Numerical Recipes , 1st ed. (CambridgeUniversity Press, Cambridge, England, 1989). B. Fornberg, Math. Comp. , 699 (1988). C. Eckart, Phys. Rev. , 1303 (1930). Z. Ahmed, Phys. Rev. A , 4761 (1993). S. Flugge,