Complex trajectory method in time-dependent WKB
aa r X i v : . [ qu a n t - ph ] J u l Complex trajectory method in time-dependent WKB
Yair Goldfarb , Jeremy Schiff , , and David J. Tannor Department of Chemical Physics, Department of Mathematics,The Weizmann Institute of Science,Rehovot, 76100 IsraelNovember 1, 2018
Abstract
We present a significant improvement to a time-dependent WKB (TDWKB) formulation developedby Boiron and Lombardi [JCP , 3431 (1998)] in which the TDWKB equations are solved alongclassical trajectories that propagate in the complex plane. Boiron and Lombardi showed thatthe method gives very good agreement with the exact quantum mechanical result as long as thewavefunction does not exhibit interference effects such as oscillations and nodes. In this paper weshow that this limitation can be overcome by superposing the contributions of crossing trajectories.We also demonstrate that the approximation improves when incorporating higher order terms inthe expansion. These improvements could make the TDWKB formulation a competitive alternativeto current time-dependent semiclassical methods. On sabbatical leave from Dept. of Mathematics, Bar-Ilan University, Ramat Gan 52900, Israel.
PACS numbers: . INTRODUCTION The difficulty in performing quantum mechanical calculations of multi-dimensional systemshas stimulated an intensive and ongoing effort in the last three decades to develop numericaltools based on semiclassical mechanics. In this context, we refer to semiclassical mechanicsas the derivation of a quantum mechanical wavefunction or propagator via propagation ofclassical (or classical-like) trajectories. From a physical point of view, semiclassical methodstry to evade the non-locality imbedded in quantum mechanics. Mathematically speaking,semiclassical methods aim at casting the time-dependent Schr¨odinger equation (TDSE),which is a PDE, in terms of ODEs related to classical equations of motion. This trans-formation has significant computational advantages that can ease the inherent difficulty ofmulti-dimensional quantum calculations.The WKB method[1, 2, 3] can be considered as the first of the semiclassical methods. Itsdate of birth almost coincides with the publication of the Schr¨odinger equation in 1926, andvirtually every standard text book in quantum mechanics has a description of the method.The basic idea of the WKB method is to recast the wavefunction as the exponential of afunction and then replace the exponent with a power series in ~ . The WKB method isordinarily applied to the time-independent Schr¨odinger equation and provides for a goodapproximation to the eigenstates as long as one is not too near a classical turning point. Itis only natural that as part of the effort to develop time-dependent semiclassical methods,a time-dependent version of the WKB method would be explored. Surprisingly little workhas been done in this direction[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. A decade ago,Boiron and Lombardi[17] developed a complex trajectory version of time-dependent WKB(TDWKB), which we refer to as CTDWKB. In conventional WKB the leading order termin the phase of the wave function is taken to be O ( ~ − ) and the leading order term in theamplitude is taken to be O ( ~ ). In contrast, the CTDWKB formulation treats the amplitudeand phase on an equal footing. The price to pay for this procedure is that the resultingclassical trajectories propagate in the complex plane. The benefits are that the results aresuperior to standard TDWKB and no singularities are encountered during the integrationof the equation of motion.The CTDWKB equations of motion can be solved analytically and yield the exact wave-function for an initial Gaussian wavepacket in a potential with up to quadratic terms. The2rst-order method was tested numerically by Boiron and Lombardi for scattering of a Gaus-sian wavepacket from a potential barrier. They showed that the method produced verygood results as long as the wavefunction did not exhibit interference effects in the form ofoscillations or nodes[17]. In this paper we present a simple modification to CTDWKB thatprovides an accurate description of oscillations in the wavefunction. We show that complexclassical trajectories, similar to real classical trajectories, can cross in configuration space.By superposing the contributions from two or more crossing trajectories, interference effectsare obtained. We take CTDWKB a step further in another direction by showing that the ap-proximation generally improves when incorporating additional terms in the series expansion.Since the WKB expansion is an asymptotic series, this observation is non-trivial.Two other semiclassical formulations that incorporate complex trajectories should bementioned in relation with CTDWKB. The first is the Generalized Gaussian WavepacketDynamics (GGWPD) developed by Huber, Heller and Littlejohn [18, 19]. One may showthat for an initial Gaussian wavepacket the equations of motion of GGWPD are de factoidentical to the equations of the first-order approximation of CTDWKB. However the GG-WPD has no generalization to arbitrary initial wavefunctions and no systematic way toincrease the accuracy of the approximation. On the other hand, in reference [18], Huber andHeller appreciate the importance of multiple complex trajectories in obtaining interferencephenomena. Here we incorporate the idea of crossing complex trajectories into the moregeneral CTDWKB formulation.The second formulation that is closely related to CTDWKB is Bohmian Mechanics withComplex Action (BOMCA)[20, 21]. CTDWKB and BOMCA begin with the same ansatzof substituting of an exponential function exp[ iS/ ~ ] into the TDSE. Similar to CTDWKB,the BOMCA formulation incorporates equations of motion that propagate along complextrajectories. The first-order equations of motion of BOMCA are identical to the equations offirst-order CTDWKB. The differences between the two formulations are: (1) The equationsof motion in BOMCA are for the coefficients of spatial derivatives of the phase. In CTDWKBthe equations of motion are for the coefficients of an ~ Taylor expansion of the phase and theirspatial derivatives. (2) Incorporating higher order terms of the CTDWKB approximationdoes not effect the the results for lower order terms since each equation of motion dependsonly on lower terms of the expansion. This is not the case with BOMCA where each equationof motion depends on both lower and higher terms resulting in a backward feedback. (3) A3esult of the last difference is that in CTDWKB the equations of motion of the trajectoriesremain classical whereas in BOMCA, the inclusion of higher orders of the approximationaffect the complex trajectories by adding a “quantum force” that yields quantum trajectories.This paper is organized as follows. In section II we formulate TDWKB and the CTD-WKB. Our derivation is more compact than the Boiron-Lombardi derivation and demon-strates how to obtain the equations of motion for higher orders of the expansion in a simplemanner. In Section III we apply the formulation to a Gaussian initial wavepacket propagat-ing in a quartic double-well potential. We demonstrate that superimposing the contributionsof crossing trajectories leads to interference effects and that incorporating higher order termsin the expansion improves the approximation. Section IV is a summary and concluding re-marks. Following Boiron and Lombardi we will refer to the CTDWKB method in the bodyof the paper as the complex trajectory method (CTM) for short.
II. FORMULATIONA. Time-independent vs. Time-dependent WKB
For simplicity we present the one-dimensional version of the CTM derivation. The general-ization to multi-dimensions can be performed in a straightforward manner. The conventionalWKB derivation begins by inserting the ansatz ψ ( x ) = exp (cid:20) i ~ S ( x ) (cid:21) , (2.1)into the time-independent Schr¨odinger equation, where ~ is Planck’s constant divided by2 π . The end result is 12 m (cid:18) dSdx (cid:19) + V ( x ) − i ~ m d Sdx = E, (2.2)where m is the mass of the particle, V ( x ) is the potential energy and E is the eigenvalue. Ifwe assume that S ( x ) can be expanded asymptotically as a polynomial in ~ S ( x ) = S ( x ) + ~ S ( x ) + ~ S ( x ) + ... = ∞ X j =0 ~ j S j ( x ) , (2.3)then, by substituting the last equation into eq.(2.2) and equating powers of ~ , a series ofcoupled ODEs are obtained for the S j ’s. 4he time-dependent WKB begins by inserting the ansatz[4, 5] ψ ( x, t ) = exp (cid:20) i ~ S ( x, t ) (cid:21) , (2.4)into the time-dependent Schr¨odinger equation, i ~ ∂ t ψ = − ~ m ∂ xx ψ + V ( x, t ) ψ, (2.5)The result is the quantum Hamilton-Jacobi equation[4, 5] ∂ t S + 12 m ( ∂ x S ) + V = i ~ m ∂ xx S, (2.6)where the LHS of the equation is in the form of the classical Hamilton-Jacobi equation.Equation (2.6) is formally exact since no approximation has been introduced. In TDWKBformulation we insert into eq.(2.6) a time-dependent version of eq.(2.3) S ( x, t ) = ∞ X j =0 ~ j S j ( x, t ) . (2.7)The result is ∞ X j =0 ~ j ∂ t S j + 12 m ∞ X j ,j =0 ~ j + j ∂ x S j ∂ x S j + V = i m ∞ X j =0 ~ j +1 ∂ xx S j . (2.8)By equating terms having the same powers of ~ we obtain the classical Hamilton-Jacobiequation for S ( x, t ) ∂ t S + 12 m ( ∂ x S ) + V = 0 , (2.9)and equations of motion for S n ( x, t ) , n ≥ ∂ t S n + ∂ x S m ∂ x S n = i m ∂ xx S n − − m n − X j =1 ∂ x S j · ∂ x S n − j . (2.10)Conveniently, each equation depends only on lower order terms. The next step in TDWKBis to convert eqs.(2.9) and (2.10) into a set of ODEs by looking at the evolution of S , S , . . . along classical trajectories, as described in the next section. B. Integrating along classical trajectories
As we mentioned earlier, the first term in the ~ power expansion, S , obeys the Hamilton-Jacobi equation (eq.(2.9)). This equation is an alternative formulation of Newton’s second5aw of motion in terms of an action field. The emergence of classical trajectories in theTDWKB equations provides the incentive to solve eqs.(2.9) and (2.10) by integrating alongsuch trajectories.The link between the Hamilton-Jacobi equation and classical trajectories is demonstratedby defining the velocity field v ( x, t ) ≡ ∂ x S ( x, t ) m (2.11)and considering the trajectories defined by dxdt = v ( x, t ) . (2.12)By taking the spatial partial derivative of eq.(2.9), using the definition of the Lagrangiantime derivative ddt ≡ ∂ t + dxdt ∂ x , and applying eq.(2.11) we obtain the equation of motion forthe velocity along a trajectory as Newton’s second law dvdt = − ∂ x Vm . (2.13)Hence, the trajectories defined are simply classical trajectories.Inserting eq.(2.9) in the Lagrangian time derivative of S yields dS dt = ∂ t S + ∂S m ∂ x S = 12 mv − V, (2.14)where we recognize the equation of motion for the action along a classical trajectory. Notingthat v is a mere dummy variable, we summarize the equations of motion for the zeroth orderterm of TDWKB, S dxdt = ∂ x S m , (2.15) d ( ∂ x S ) dt = − ∂ x V, (2.16) dS dt = 12 m ( ∂ x S ) − V. (2.17)We turn to the higher order terms in the series S n , n ≥
1. Recognizing the LHS of eqs.(2.10)as the Lagrangian time derivative of S n , we can write dS n dt = i m ∂ xx S n − − m n − X j =1 ∂ x S j · ∂ x S n − j . (2.18)These equations do not constitute a closed set of ODEs since they depend on partial deriva-tives such as ∂ xx S n . We close the set of equations by deriving equations of motion for the6artial derivatives on the RHS of eq.(2.18) ( ∂ xx S n − , ∂ x S j and ∂ x S n − j ). We demonstrate theprocess by deriving equations of motion for S and S . Inserting n = 1 in eq.(2.18) yields dS dt = i m ∂ xx S . (2.19)An equation of motion for ∂ xx S is obtained by taking a second spatial partial derivative ofeq.(2.9), ∂ xxt S + 1 m (cid:2) ∂ x S · ∂ xxx S + ( ∂ xx S ) (cid:3) + ∂ xx V = 0 , (2.20)and rewriting it as d ( ∂ xx S ) dt = − m ( ∂ xx S ) − ∂ xx V. (2.21)This equation is derived in reference [17] by a cumbersome finite difference scheme. Itis equivalent to eq.(2.9d) of reference [19] where the equation appears in the context ofGGWPD. Note that an equation of motion for any order of spatial derivatives of S can bederived in a similar fashion by taking consecutive spatial derivatives of eq.(2.20) and thengrouping the Lagrangian time derivative terms. Equations (2.19) and (2.21) provide a closedset of equations of motion for S .Inserting n = 2 into eq.(2.18) yields dS dt = i m ∂ xx S − m ( ∂ x S ) . (2.22)The equations of motion for ∂ x S and ∂ xx S are obtained by first inserting n = 1 in eq.(2.10).We then derive two equations by taking a first and a second spatial partial derivative of theresult. By grouping the Lagrangian time derivatives of ∂ x S and ∂ xx S in each of the twoequations separately we obtain d ( ∂ x S ) dt = i m ∂ xxx S − m ∂ x S · ∂ xx S , (2.23) d ( ∂ xx S ) dt = i m ∂ xxxx S − m ∂ x S · ∂ xxx S − m ∂ xx S · ∂ xx S . The last equations depend in turn on ∂ xxx S and ∂ xxxx S . As mentioned earlier, the equationof motion for these terms can be obtained by additional spatial derivatives of eq.(2.20), aprocess that yields d ( ∂ xxx S ) dt = − m ∂ xx S · ∂ xxx S − ∂ xxx V, (2.24) d ( ∂ xxxx S ) dt = − m (cid:2) ∂ xx S · ∂ xxxx S − ∂ xxx S ) (cid:3) − ∂ xxxx V. S . Thescheme we described for S and S can be extended to any of the higher order terms in theexpansion. Note that incorporating higher order terms S n in the TDWKB approximationdoes not affect the classical trajectories associated with S , defined by eqs.(2.15) and (2.16).We now turn to the source of the distinction between conventional TDWKB and CTM. C. Initial conditions and complex classical trajectories
In conventional TDWKB the initial wavefunction is “divided” between S ( x,
0) and S ( x, ψ ( x,
0) = A ( x ) exp[ iφ ( x )] = exp (cid:20) i ~ S ( x,
0) + S ( x, (cid:21) , (2.25)where A ( x ) and φ ( x ) are the initial amplitude and phase respectively, both taken to be real.The phase is related to the zero-order term S and the amplitude to the first-order correctionterm S according to S ( x,
0) = ~ φ ( x ) , S ( x,
0) = − i ln[ A ( x )] , (2.26)and S n ( x,
0) = 0 for n ≥
2. Note that the initial conditions specified by eqs.(2.26) yieldclassical trajectories that propagate on the real axis since S and its spatial derivatives arereal quantities (see eqs.(2.15) and (2.16)). In contrast, in CTM the amplitude and phaseare treated on an equal footing with far-reaching consequences. The initial wavefunction isspecified by S ( x, S ( x,
0) = − i ~ ln[ ψ ( x, , S n ( x,
0) = 0 , n ≥ . (2.27)Since S is generally complex and since the initial velocity v ( x, ≡ ∂ x S ( x, /m , thetrajectories propagate in the complex plane even if the initial positions are on the real axis( ℑ [ x (0)] = 0). This observation requires us to look at the analytic continuation of thewavefunction in the complex plane and find ways to extract the wavefunction on the realaxis. D. Complex root search and superposition
One of the benefits of conventional TDWKB and CTM compared with BOMCA, is that thetrajectories obey the classical equations of motion and are independent of the order of the8hase expansion we incorporate in the final wavefunction. But the fact still remains thatfor an arbitrary initial position x (0) ∈ C and an arbitrary final propagation time t f the finalposition x ( t f ) is complex and yields an “analytically continued” wavefunction at x ( t f ) ψ [ x ( t f ) , t f ] ≈ exp ( i ~ N X j =0 ~ j S j [ x ( t f ) , t f ] ) , (2.28)where the non-negative integer N is the order of the approximation. References [17, 19,20] include discussions of root search algorithms for the derivation of initial positions thatreach the real axis at a given time. We will not describe all the details here but will juststate the central idea. The complex root search exploits the assumption that the mapping x (0) x ( t f ) is analytic. This property allows for an iterative process that detects the initialpositions that correspond to real final positions. As demonstrated in references [18, 19] andin section III A, for an arbitrary potential and final time, the mapping is only locally analytic.Generally, more than one initial position ends at a final position (whether real or complex).This makes the search for trajectories that end on the real axis more complicated but it hasan important advantage in terms of interference effects.Our main observation is that the contribution of multiple trajectories in CTM can accu-mulate to an interference pattern. For simplicity we make the following assumption. Supposethat L trajectories end at final time t f on real position x ( t f ) Then the final wavefunction isapproximated by a superposition of contributions ψ [ x ( t f ) , t f ] ≈ L X l =1 exp (cid:26) i ~ S l [ x ( t f ) , t f ] (cid:27) , (2.29)where each trajectory (denoted by the index l ) is associated with a phase S l [ x ( t f ) , t f ] S l [ x ( t f ) , t f ] = N X j =0 ~ j S lj [ x ( t f ) , t f ] , (2.30)that is calculated by the CTM equations of motion. In section III we show that this as-sumption is too simplified and does not hold at all times and all positions. For example,for positions associated with a tunneling part of the wavefunction, only one of the multipletrajectories should be taken into account. A partial discussion on the superposition of con-tributions from complex trajectories appears in reference [19] in the GGWPD context. In aforthcoming paper [22] we will explain an alternative derivation of the CTM in which theneed to include multiple trajectories for certain times and positions becomes apparent9 II. NUMERICAL RESULTS
In this section we examine numerically the CTM formulation allowing for the superpositionof complex trajectories. For ready comparison the physical system we choose is identicalto the one studied by Boiron and Lombardi (reference [17] section IVB). The potentialconsidered is a quartic double-well V ( x ) = 1 . × − ( x − x ) . (3.1)The initial wavefunction is a Gaussian wavepacket ψ ( x,
0) = exp (cid:20) − α ( x − x c ) + i ~ p c ( x − x c ) + i ~ γ (cid:21) , (3.2)where α = 1, x c = 0, p c = 5, γ = − i ~ ln( α π ) and we take m = ~ = 1 (all quantities aregiven in atomic units). The initial conditions for the terms in the ~ power-expansion of thephase are S ( x,
0) = iα ~ ( x − x c ) + p c ( x − x c ) + γ = ix + 5 x + γ , (3.3) ∂ x S ( x,
0) = 2 iα ~ ( x − x c ) + p c = 2 ix + 5 , (3.4) ∂ xx S ( x,
0) = 2 iα ~ = 2 i, (3.5) ∂ jx S ( x,
0) = 0 , j ≥ , (3.6) ∂ jx S k ( x,
0) = 0 , j ≥ , k ≥ , (3.7)where ∂ jx S k ≡ ∂ j S k ∂x j .In section III A we analyze the first order approximation of CTM ( N = 1 , S = S + ~ S )and the properties of the trajectories. Section III B is dedicated to the next order of theapproximation ( N = 2 , S = S + ~ S + ~ S ). We omit an analysis of N = 0 since it is wellpresented in reference [17] and only yields poor results. A. First Order approximation, N = 1 The first order approximation of CTM requires the solution of eqs.(2.15), (2.16), (2.17),(2.19) and (2.21). The first two equations define the complex classical trajectories and thenext three equations yield S and S . We start by analyzing the complex classical trajec-tories. As mentioned above, the mapping x (0) x ( t f ) is not one-to-one. For the quartic10otential, we found that three initial positions are mapped to every real final position at t f >
0. For short time scales this observation can be supported analytically. For general po-tentials or for longer time scales than we present here, more than three initial positions mightlead to the same final position[19, 21]. In figures 1(a) and 1(b) we plot complex classicaltrajectories for t f = 3 and t f = 6 respectively. The initial positions of the trajectories canbe divided into three groups referred to as branches [19]. One group of the initial positionsis called the real branch and the other groups are called the secondary branches. The realbranch is characterized by the property that it includes the initial position of a trajectorythat propagates solely on the real axis. We refer to this trajectory as the real trajectory.It can be readily verified that for a Gaussian initial wavefunction there is only a single realtrajectory that initiates at x (0) = x c (see eqs.(2.16), eq.(2.15) and (3.4)). In fig.1(b) we de-pict the real trajectory explicitly. The secondary branches are defined simply as the groupsof initial positions that do not belong to the real branch. Generally, the branches might beinfinitely long curves in the complex plane. We will use the term branches to refer to thelocus of initial positions that leads to final positions where the wavefunction is significantlydifferent from zero. Hence, the branches are curves of finite length in the complex plane,although clearly there is some arbitrariness to their length.In fig.1(a) we see that at short time scales the secondary branches are centered far fromneighborhood of the real axis. We can show analytically that for small times t f the initialpositions that comprise the real branch obey | x (0) | = O ( t f ) whereas the secondary branchesobey | x (0) | = O ( t f ). Note that the linear dependence of the initial momentum on position(eq.(3.4)) allows trajectories with initial positions far from the real axis to reach a real finalposition in a short time. Unlike the secondary branches, the real branch is centered in thevicinity of the real axis at all times. The initial position x (0) = x c is a fixed point of thereal branch and prevents the real branch (recall that this is the locus of initial positions)from “straying” from the neighborhood of the real axis as the final time t f is increased. Atintermediate times (time scales comparable to the time of the collision of the wavefunctionwith the barrier, 4 < ∼ t f < ∼
7) secondary branch (1) reaches the vicinity of the real axis(fig.1(b)) and at longer time scales it continues in the direction of the positive imaginaryaxis. As we demonstrate below, the proximity of secondary branch (1) to the real axis isclosely related to the size of its contribution to the final form of the wavefunction and itsrole in interference effects. Secondary branch (2) does not reach the vicinity of the real axis11or any of the time scales specified below. The contribution of this branch to the absolutevalue of the final wavefunction (eq.(2.29)) is negligible (in the order of 10 − ). Hence, fromhere on we ignore secondary branch (2) and refer to secondary branch (1) as the secondarybranch.As we mentioned in section II D, the existence of more than one branch motivates theattempt to superpose the contributions of the real branch and secondary branch in the finalwavefunction ψ [ x ( t f ) , t f ] = ψ R + ψ S ; ψ R = exp (cid:20) i ~ S Real (cid:21) , ψ S = exp (cid:20) i ~ S Sec (cid:21) , (3.8)where S Real and ψ R are the phase and wavefunction associated with the real branch, and S Sec and ψ S correspond to the secondary branch. In figures 2(a), 2(b) and 2(c) we comparethe exact wavefunction with the numerical results obtained by applying CTM using a two-branch superposition. The figures indicate that when the wavefunction does not exhibitoscillations, the contribution of the real branch is sufficient to obtain a good approximationto the wavefunction. But at intermediate times, when the wavefunction exhibits interferenceeffects, the contribution of both branches must be included. This last last observation appliesin the spatial range up to the classical turning point ( x ≃ t f = 6 of each individual branch and their superposition. Starting from the vicinity of x ≃
22, we observe an exponential increase of ψ S . For x > ∼
23 we have a discontinuity of theapproximation, as we discard the contribution of the secondary branch and include just thereal branch. A description of this divergence appears in reference [19] in the context of theGGWPD formulation.It is interesting to compare the time-dependence of the real and secondary branch con-tributions to the final approximation. A qualitative measure of the contribution of eachbranch is given by the imaginary part of the phase since | ψ R | = (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) i ~ S Real (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) = exp (cid:20) − ℑ ( S Real ) ~ (cid:21) , (3.9)and a similar relation applies for ψ S and ℑ ( S Sec ). In figures 4(a) and 4(b) we plot ℑ ( S Real ) and ℑ ( S Sec ) respectively for a series of final propagation times. We see that the secondary branchhas a significant magnitude only at intermediate times. This observation coincides well with12
60 −40 −20 0 20 40 60−30−20−100102030
Real(x) I m ag ( x ) t f =3 Real BranchSecondary Branch (2) (a)
Secondary Branch (1) −5 0 5 10 15 20 25 30−2−1012
Real(x) I m ag ( x ) t f =6 Secondary Branch (1) Real Branch (b)
Real Trajectory
FIG. 1: Complex classical trajectories with initial positions marked as circles and real final positions(marked as pluses) at (a) t f = 3 and (b) t f = 6. The trajectories arise from an initial Gaussianwavepacket propagating in a quartic double-well potential. The Gaussian is centered at x = 0and has positive initial momentum (the physical parameters are given in the text). In plot (a) wedemonstrate that each final position arises from three initial positions. The initial positions aredivided into a real branch and two secondary branches. The real branch is defined as incorporatinga trajectory that remains on the real axis at all times. The real trajectory is specifically indicatedin plot (b). x | ψ | Exact QM| ψ R | | ψ R + ψ S | (1) (2) (3) (4) (5) (a) x | ψ | (9) (8) (7) (6) (b)
16 18 20 22 2400.050.10.150.20.25 x | ψ | (8)(7) (6) (c) FIG. 2: A comparison between the exact quantum wavefunction and CTM ( N = 1) with a two-branch superposition. The comparison is at a series of final propagation times specified by thenumbers in the parentheses. The plots arise from an initial Gaussian wavepacket centered at x = 0with a positive average momentum, propagating in a quartic double-well potential (the parametersare given in the text). (a) Initially right-propagating wavefunction; (b) the reflected wavefunction;(c) a zoom on a section of (b). For t f = 5 in (a) and t f = 6 in (b) and (c) we plot the results bothfor just the real branch | ψ R | and for the combination of branches | ψ R + ψ S | . The interferencepattern obtained by superposing the contributions is clearly observed. the need to include the contribution of the secondary branch to the final wavefunction onlyat these times. The exponential growth of ψ S that is observed in fig.3 is also apparent infig.4(b), in the negative parts of the graphs for t f = 5 and t f = 6. The divergent magnitudeof ψ S is in contrast to the finite magnitude of ψ R that is observed in fig.4(a). A discontinuity14 x | ψ | t f =6, N=1 Exact QM| ψ R || ψ S || ψ R + ψ S | FIG. 3: CTM approximation at t f = 6 for N = 1. The contributions of each branch to thewavefunction are depicted by plotting | ψ R | , | ψ S | and | ψ R + ψ S | . Note the exponential increase of ψ S begins around x ≃
22. For x ≃
23; we discard the contribution of the secondary branch andinclude just the real branch, leading to a discontinuity of the CTM approximation. in the derivative of ℑ ( S R ) at t f = 5 and t f = 6 is also observed. This discontinuity appearsslightly prior to the points where the contribution of the secondary branch begins to diverge.A close inspection of the complex trajectories at t f = 5 and t f = 6, reveals an interestingproperty of the real trajectory: the real trajectory acts as a boundary between two “regimes”of complex trajectories comprising from the real branch. This can be seen in fig.5, where thetrajectories that initiate from ℑ [ x (0)] > x -axis at values lowerthan the real trajectory while trajectories with ℑ [ x (0)] < x -axis at values higher than the real trajectory. These two regimes correspondto the two legs of the “v”-shaped graph of ℑ ( S Real ) in fig.4(a): the trajectories arising frominitial positions with ℑ [ x (0)] > ℑ [ x (0)] < . Second Order approximation, N = 2 In this section we analyze the effect of incorporating S in the CTM approximation. Inaddition to the five equations that are needed for obtaining the complex trajectories, S and S , we need to solve eqs.(2.22), (2.23) and (2.24). In fig.6(a) we depict the approximatewavefunction for N = 2 at t f = 6. Comparing the N = 2 result with the N = 1 resultplotted in fig.3, we conclude that other than an interval in the neighborhood of x ≃ . N = 2 result (dashed line) lies on top of the exact result (solid line) and is significantlybetter then the N = 1 result. For x > ∼
23, where we incorporate solely the real branchcontribution, the improvement in the approximation is graphically evident from the plots.For x ≤
22 we calculated the relative error between the absolute value of the approximationsand the exact wavefunction using all the data points depicted in figs.3 and 6(a). The resultsare presented in fig.6(b). For N = 1 the mean relative error is 0 .
34% while for N = 2 themean relative error is 0 . ψ R . In the vicinity of x ≃ .
5, the N = 2 results are worse than the N = 1results; moreover, in the N = 2 case ψ S as well as ψ R exhibits a discontinuity. IV. SUMMARY
In this paper we have presented a formulation of complex time-dependent WKB (CTDWKB)that allows the incorporation of interfering contributions to the wavefunction. The centralidea in CTDWKB presented by Boiron and Lombardi[17] is to include both the amplitudeand the phase in the lowest order term of the conventional time-dependent WKB method.The rationale behind this substitution is to treat the phase and the amplitude on equalfootings in the limit ~ →
0. The benefits of the method are twofold. Firstly, CTDWKBexhibits accuracy superior to the conventional TDWKB[17]. Secondly, no singularities ap-pear in the integration of the equations of motion. The method has two main drawbacks .First, the trajectories that emerge obey the classical equations of motion but propagate inthe complex plane (due to complex initial conditions), requiring analytic continuation of thequantum wavefunction. The second drawback is that the reconstruction of the wavefunctionon the real axis requires a root search process. This process can be eased by exploiting theanalytic mapping between initial and final position.16e have incorporated into the CTDWKB method the possibility of contributions frommultiple crossing trajectories. Boiron and Lombardi claim (section V in reference[17]) thatthey use the root search procedure “excluding de facto such double contributions”, althoughthey appreciate the benefit that double contributions have in the GGWPD formulation. Aswe have demonstrated here, considering double contributions allows description of interfer-ence effects that are missing in the Boiron-Lombardi formulation of CTDWKB. Moreover,we have showed how to derive higher orders terms of the approximation in a straightforwardmanner. This process was applied for the derivation of a second order term in the CTDWKBapproximation. The results for N = 2 were better than for N = 1 except for a small intervalin the vicinity of the classical turning point. It was also observed that even though thereare no singularities in the integration of the CTDWKB equations of motion, a singularityappears in the real branch ψ R at intermediate times. For N = 2 an irregularity also appearsin the part of the wavefunction associated with the secondary branch ψ S . We demonstratedthat when a singularity appears in ψ R (at intermediate times), the real trajectory acts asthe boundary between two groups of trajectories associated with the real branch. Each ofthese groups contributes to a different side of the singularity.The CTDWKB formulation has several issues that require more comprehensive study.The most critical issue is to give an analytic explanation of the need to include the con-tributions from multiple classical trajectories (with zero relative phase) and why in somecases these contributions diverge. This will be dealt with in our forthcoming publication[22]. Some insight into the analytic structure of the complex classical trajectories was givenin reference [19] in the context of GGWPD; however, we believe that a more general under-standing of this structure is yet to be developed. This structure presumably is relevant tothe question of when the CTDWKB formulation converges to the exact quantum mechani-cal result. We saw that in most parts of configuration space N = 2 performed better than N = 1, but in other parts of configuration space, where there were singularities, N = 2performed worse. What determines the position and time-dependence of these singularitiesin ψ R at intermediate times? What is the relation between the singularities in CTDWKB vs.conventional time-dependent WKB? Is there any fundamental limitation on the time scalefor which the method is accurate? Since WKB plays such a central role in quantum mechan-ics in general and in semiclassical mechanics in particular, we believe that these questionsare of great general interest. The developments described in this paper together with the17nswers to some of the above questions could make the time-dependent WKB formulationa competitive alternative to current time-dependent semiclassical methods.We wish to acknowledge David Kessler and Uzi Smilansky for useful discussions. Thiswork was supported by the Israel Science Foundation (576 / [1] G. Wentzel, Z. Phys. , 518 (1926).[2] H. A. Kramers, Z. Phys. , 828 (1926).[3] L. Brillouin, CR Acad. Sci, Paris , 24 (1926); L. Brillouin, J. Phys., , 353 (1926)[4] W. Pauli, Die allgemeine Prinzipien der Wellenmechanik , Encyclopedia of Physics, Vol. 5,Springer, Berlin, 1958.[5] K. Gottfried,
Quantum Mechanics, Volume I: Foundations (W. A. Benjamin, New York, 1966)[6] Byung Chan Eu, W. J. Chem. Phys. , 2531 (1972).[7] H. J. Korsch, R. Mohlenkamp, Phys. Lett. A , 110 (1978).[8] S. M. Blinder, Chem. Phys. Lett. , 288 (1987).[9] L. Raifeartaigh, A. Wipf, Found. Phys. Lett. , 307 (1987).[10] M. P. A. Fisher, Phys. Rev. B , 75 (1988).[11] R. Rubin, Junior paper (Princeton, 1990).[12] M. Burdick, H. J. Schmidt, J. Phys. A: Math. Gen. , 579 (1994).[13] C. Sparber, P. A. Markowich, N. J. Mauser, arXiv:math-ph/0109029 (2002).[14] A. S. Sanz, F. Borondo, S. Miret-Art´es, J. Phys.:Condens. Matter , 6109 (2002).[15] Jeong Ryeol Choi, Int. J. Theo. Phys. , 947 (2004).[16] P. Bracken, arXiv:math-ph/0608011v2 (2006).[17] M. Boiron, M. Lombardi, J. Chem. Phys. , 3431 (1998).[18] D. Huber, E. Heller, J. Chem. Phys. , 5302 (1987).[19] D. Huber, E. Heller, R. G. Littlejohn, J. Chem. Phys. , 2003 (1988).[20] Y. Goldfarb, I. Degani, D. J. Tannor, J. Chem. Phys. , 231103 (2006).[21] Y. Goldfarb, D. J. Tannor, submitted.[22] J. Schiff, Y. Goldfarb, D. J. Tannor, in preparation.
10 20 30012345 ℑ ( S R ea l ) x (9) (5)(6) (4)(7)(8) (a) ℑ ( S S e c ) x (9) (8) (6) (5)(4)(7) (b) FIG. 4: (a) ℑ ( S Real ) and (b) ℑ ( S Sec ) are depicted at a series of final propagation times (given inthe parentheses). The results are limited to the spatial interval for which the absolute value of theexact wavefunction is significantly larger than zero. The imaginary part of the phase allows for aqualitative estimate of the contribution of each branch to the probability | ψ R + ψ S | , see eq.(3.9).Figure (b) shows that ℑ ( S Sec ) drops below ∼ Real(x) I m ag ( x ) Real Branch, t f =5 Initial positionsFinal positionsTrajectoriesReal Trajectory
FIG. 5: Complex classical trajectories that correspond to the real branch at t f = 5. The realtrajectory acts as a boundary between two “regimes” of the complex trajectories. Initial positionswith ℑ [ x (0)] < ℑ ( S Real ) at intermediate times (fig.4(a)). x | ψ | t f =6, N=2 Exaxt QM| ψ R || ψ S || ψ R + ψ S | (b)
14 16 18 20 2210 −4 −3 −2 −1 x R e l a t i v e E rr o r ( % ) N=1N=2 (b)
FIG. 6: (a) The second order ( N = 2) CTM approximation is depicted for t f = 6. A discontinuityappears at x ≃ . ψ R and ψ S . (b) The relative error between the absolute value of theexact quantum wavefunction and the CTM approximation for N = 1 and N = 2, based on thedata in fig.3 and fig.6(a). A comparison of the relative errors indicates a clear improvement whentaking an additional order in the CTM approximation.= 2, based on thedata in fig.3 and fig.6(a). A comparison of the relative errors indicates a clear improvement whentaking an additional order in the CTM approximation.