Table Of Content

Accurate numerical solutions of the time-dependent Schrödinger equation W. van Dijk∗ Physics Department, Redeemer University College, Ancaster, Ontario L9K 1J4, Canada and Department of Physics and Astronomy, McMaster University, Hamilton, Ontario L8S 4M1, Canada F. M. Toyama† 7 Department of Information and Communication Sciences, 0 Kyoto Sangyo University, Kyoto 603-8055, Japan 0 (Dated: February 2, 2008) 2 Wepresentageneralizationoftheoften-usedCrank-Nicolson(CN)methodofobtainingnumerical n solutionsofthetime-dependentSchrödingerequation. Thegeneralizationyieldsnumericalsolutions a accuratetoorder(∆x)2r−1inspaceand(∆t)2M intimeforanypositiveintegersrandM,whileCN J employ r=M =1. Wenote dramatic improvement in theattainable precision (circa 10 or greater 2 ordersofmagnitude) along with severalorders ofmagnitude reductionof computational time. The 1 improved method is shown to lead to feasible studies of coherent-state oscillations with additional short-range interactions, wavepacket scattering, and long-time studiesof decaying systems. ] h PACSnumbers: 02.60.-x,02,70.-c,03.67.Lx,03.65.-w p - p m I. INTRODUCTION cesses [9]. However, the methods used in the past, and still employedcurrently,are limited inthat the solutions o c Whereas there are a number of examples of exact an- often degrade after a certain time interval, so that they . reduce to noise. Furthermore for processes in which the s alytic solutions of time-independent problems in quan- c wave function spreads or travels away from the source tum mechanics, such solutions of time-dependent prob- i one often requires such a large number of space steps s lems are few. The analytic solutions of both types that y that the computation becomes prohibitive. do exist tend to provide approximate models to actual h physical systems. No doubt such solutions are instruc- The goal of this study is to improve the existing stan- p [ tive for gaining insight into the behavior of the physi- dard approach by allowing for relatively large step sizes calsystemsthatthey describe. Neverthelessbecause the bothintimeandinspaceandthustoreducethenumber 1 modelsthemselvesareoftenapproximationsandbecause of basic arithmetical calculations while obtaining more v one wishes to describe real systems as precisely as possi- accurate solutions. We have been able to make signifi- 0 ble, one also relies on accurate numerical methods. cantimprovementtooneconventionalapproach,viz.,the 5 Crank-Nicolson(CN) implicit integrationscheme for the 1 The time-dependence of nonrelativistic quantum sys- 1 tems, which is the focus of this paper, has become time-dependent Schrödinger equation. Many years ago 0 theCNapproachwasshowntobesuccessfulinthestudy important in diverse areas of atomic and subatomic 7 of wavepacket scattering in one dimension by Goldberg physics. Examples of these include the study of nu- 0 et al. [8]. In recent years the CN method continues to clear processes such as the decay of unstable nuclei and / s associated phenomena like atomic ionization [1, 2] and be employed for its space- and/or time-development al- c gorithmto study various time-dependent problems. See, bremsstrahlung [3, 4, 5], the study of fundamental pro- i s cessesnecessaryforquantumcomputing[6],thestudyof for example, Refs. [10, 11, 12, 13]. The attractive as- y pect of this method is that the solutionis constrainedto mesoscopic physics or nanophysics devices [7], and the h be unitary at every time step. It is this constraint that p motion of atoms in a trap. A reliable and accurate nu- makesthesolutionstableregardlessoftime-orspace-step : merical determination of the time-dependent wave func- v size. Although the evolution of the solution is unitary, tion such as we discuss in this paper will no doubt be i the wave function is not correct if the step sizes are too X necessary and/or helpful in making advances in the un- large. Theerrorisof ((∆x)2,(∆t)3),where∆xand∆t r derstanding of a variety of quantum processes. O a are the spatial and temporal step size, respectively. In this paper we consider the numerical solution of The method was successfully generalized to two di- thetime-dependentnonrelativisticSchrödingerequation. mensional scattering by Galbraith et al. [14] and, more Much has been learned about basic scattering processes recently,tomultichannelscattering[15]. Furthermoreal- from the numerically generated solutions of traveling ternative methods which are fast computationally were wavepackets as they pass through a potential region [8], introducedbyKosloffandKosloff[16]. Theseinvolvethe as well as the time-evolution of unstable quantum pro- fast Fourier transform of the kinetic energy operator of the Schrödinger equation. Variants of this method were discussed in Ref. [17]. Although this approach is fast ∗Electronicaddress: [email protected] and is able to handle large time intervals in one pass, it †Electronicaddress: [email protected] is not unitary and does require a large number of space 2 intervals. where ImprovedCNalgorithmshavebeendiscussedbyafew authors. Mi¸sucu et al. [5] introduced a seven-point for- ¯h2 ∂2 H = +V(x), (2.2) mula for the second-order spatial derivative with error −2m∂x2 of ((∆x)6) and an improved time-integration scheme witOh an error of ((∆t)5). They claim to obtain two and φ(x) is a given wave function at initial time t . In O 0 orders of magnitude improvement in the time advance. this sectionwe use the standardtime-advance procedure Moyer [18] uses a Numerov scheme for the spatial inte- oftheCNmethod,butgeneralizethespatialintegration. gration method with error of ((∆x)6) but has a one- In Sec. III we generalize the time-evolution procedure. O stage time evolution giving the CN precision in time, The time evolution of the system can be expressed in i.e., with an error of ((∆t)3). Moyer also introduces O termsofanoperatoractingonthe wavefunctionattime transparentboundaryconditionsforunconfinedsystems. t which gives the wave function at a later time t+∆t We havefoundthese veryusefulfor making long-timeor according to the equation large-spaceproblemstractable[7],butwewillnotdiscuss suchboundaryconditions furtherinthis paper. Puzynin iH∆t/¯h et al. [19, 20] indicate how to generalize the time devel- ψ(x,t+∆t)=e− ψ(x,t). (2.3) opment to higher order, but do not discuss spatial inte- iH∆t/¯h gration. Thetime-evolutionoperatore− canbeexpanded In this paper we present a generalizationof the often- togiveaunitaryapproximationoftheoperatorbysetting used CN method of obtaining numerical solutions of the time-dependent Schrödinger equation. The gener- 1 1iH∆t/¯h alization yields numerical solutions accurate to order e−iH∆t/¯h = 1+− 12iH∆t/¯h +O((∆t)3). (2.4) (∆x)2r−1 in space and (∆t)2M in time for any positive 2 integers r and M, while CN employ r = M = 1. By Inserting the approximate form of the operator into appropriate choice of r and M the improvement can be Eq. (2.3), we obtain the equation ofsuchanaturethathithertocomputationallyunfeasible problemsbecomedoable,andsolutionswithlowtomod- est precision can now be obtained extremely accurately. 1 1 1+ iH∆t/¯h ψ(x,t+∆t)= 1 iH∆t/¯h ψ(x,t), In the following we consider the generalization of spa- 2 − 2 (cid:18) (cid:19) (cid:18) (cid:19) tial integration in Sec. II and the generalization of the (2.5) time integration in Sec. III. In Sec. IV we discuss errors with an errorof ((∆t)3). Here we focus on the second- O and a way of dealing with a particular type of bound- order spatial derivative in H of Eq. (2.2) and leave im- ary condition. We study specific examples to illustrate provementswithrespecttothetimederivativetoSec.III. theimprovementofthegeneralizationsoverthestandard We generalize the usual three-point formula and the CN procedure in Sec. V. Some general observations and seven-point formula of Mi¸sicu et al. [5], for the second- conclusions are made in Sec. VI. order derivative to a (2r+1)-point formula. Such a formula has the form II. SPATIAL INTEGRATION k=r 1 y′′(x) y(2) = c(r)y(x+kh)+ (h2r), (2.6) ≡ h2 k O We describe a general procedure for solving the one- kX=−r dimensional time-dependent Schrödinger equation (r) where c are real constants. To obtain the coefficients ∂ k i¯h H ψ(x,t)=0, ψ(x,t )=φ(x), (2.1) c(r) we make expansions ∂t − 0 k (cid:18) (cid:19) 1 1 y(x+kh) = y(x)+(kh)y(1)(x)+ (kh)2y(2)(x)+ + (kh)2r+1y(2r+1)(x)+ (h2r+2) 2! ··· (2r+1)! O 1 ( 1)2r+1 y(x kh) = y(x) (kh)y(1)(x)+ (kh)2y(2)(x) + − (kh)2r+1y(2r+1)(x)+ (h2r+2), − − 2! −··· (2r+1)! O for k =1,2,...,r; y(i) denotes the ith derivative with respect to x. When we add the two equations, the terms with odd-order derivatives cancel, resulting in the equation (kh)2 (kh)4 (kh)2r 2 y(2)(x)+2 y(4)(x)+ +2 y(2r)(x)=y(x+kh)+y(x kh) 2y(x)+ (h2r+2). (2.7) 2! 4! ··· (2r)! − − O 3 Thus we obtain the system of r equations in r unknowns, i.e., y(2k)(x) for k =1,...,r, (h)2 (h)4 (h)2r 2 y(2)(x)+2 y(4)(x)+ +2 y(2r)(x) = y(x+h)+y(x h) 2y(x) 2! 4! ··· (2r)! − − (2h)2 (2h)4 (2h)2r 2 y(2)(x)+2 y(4)(x)+ +2 y(2r)(x) = y(x+2h)+y(x 2h) 2y(x) 2! 4! ··· (2r)! − − . . . (rh)2 (rh)4 (rh)2r 2 y(2)(x)+2 y(4)(x)+ +2 y(2r)(x) = y(x+rh)+y(x rh) 2y(x). (2.8) 2! 4! ··· (2r)! − − Wesolvetheseequationstoobtainy(2)(x). Itisevident The notation includes z(1) which is consistent with that 1 fromthetermsontherightsideofEqs.(2.8)thaty(2)(x) used in the generalization of the time dependence of the has the form of Eq. (2.6) and the coefficients c(r) can wave function discussed in the next section. k be identified. Because the equations (2.8) are invariant The solution ψj,n+1 is obtained by solving the system under the change of h to h, the coefficients satisfy the of linear equations − relation c(r) = c(r) for k = 1,2,...,r. For example, −k k AΨ =A∗Ψ , (2.12) the first sevensets of coefficients (up to the fifteen-point n+1 n formula) are given in Table I. where the matrix A is the (2r+1)-diagonal matrix r k=0 1 2 3 4 5 6 7 d a a a 0 0 1 2 r 1 −2 1 ··· a d a a a  1 1 1 r−1 r  2 −5 4 − 1 ··· 2 3 12 a2 a1 d2 ar−2 ar−1 3 −49 3 − 3 1  . . . ··· . .  18 2 20 90  . . . . .  4 −205 8 −1 8 − 1  . . . . .  72 5 5 315 560   657 −−−2611556883260066800991 175372 −−−25157165 1112857069 −−−1101175028 31191275205 −−1661732 1 A= a0r aar−r1 aarr−−21 ······ adr1 dar+11 ... , 88200 4 24 108 528 3300 30888 84084      dJ−1 a1    TABLEI: The coefficients c(r) up to r=7.  a1 dJ  k  (2.13) where the superscript (r) of the a is assumed. The ma- Let us partition the range of x and t values so that k trix A∗ is the complex conjugate of matrix A. The wave x = x + j∆x, j = 0,1,...,J and t = t + n∆t, j 0 n 0 functionatt ,i.e.,Ψ ,isacolumnvectorconsisting n = 0,1,...,N. The numerical approximation of the n+1 n+1 of the ψ as components, and can be determined if wave function at a mesh point in space and time is de- j,n+1 Ψ is known. The matrix equation (2.12) can be solved noted as ψ ψ(x ,t ) and we set V = V(x ). Using n j,n ≈ j n j j using standard techniques. expression (2.6) in Eq. (2.5), we obtain k=r i¯h∆t i∆t ψ c(r)ψ + V ψ III. TIME ADVANCE j,n+1− 4m(∆x)2 " k j+k,n+1# 2h¯ j j,n+1 k=−r X k=r In this section we extend the work of Puzynin et i¯h∆t i∆t =ψ + c(r)ψ V ψ , (2.9) al. [19, 20]. The basic idea is to replace the exponen- j,n 4m(∆x)2 "k=−r k j+k,n#− 2h¯ j j,n tialoperatorexp( iH∆t)bythediagonalPadéapproxi- X − mant. The [M/M] Padéapproximantof the exponential for j = 0 to J. The indices in the sums may go out of function may be written as range, so we set ψ =0 when j <0 and j >J. Define j,n i¯h∆t b M b≡ 2m(∆x)2, z1(1) ≡−2 and a(kr) ≡ z1(1)ck(r), (2.10) f(z)=ez = a0+a1z+···+aMzM = mX=0amzm , and subsequently b +b z+ +b zM M i∆t/¯h 0 1 ··· M bm′zm′ dj ≡1+a(0r)− z(1) Vj, j =0,1,...,J. (2.11) mX′=0 (3.1) 1 4 where the am and the bm′ are complex constants. It is these equations contain no am and hence can be solved evident that when z =0, a0/b0 =1, which makes one of for the bm′, which in turn can be inserted in the first the coefficients arbitrary. By convention we take b = 1 M equations to obtain the a . The numerator and the 0 m whichimmediatelyfixesa =1. Thereare2M constants denominator of the diagonal Padé approximant of the 0 remaining, which can be found from the known coeffi- exponential function have been studied extensively [21]. cientsoftheseriesexpansionoftheexponentialfunction, When each is factored it is found that the roots of the giving an error term in Eq. (3.1) (z2M+1). The prop- denominator are the negative complex conjugates of the O ertyofPadéapproximantsthatcanbeusedtoadvantage roots of the numerator. Thus the [M/M] Padé approxi- is that, if f(z) is unitary, so is its diagonalPadéapprox- mant of the exponential function leads to imant [21]. muInltigpelynienrgalEwqe. (s3o.l1v)ebfoyrtthheedceoneoffimciineanttosrasmo tahnadtbm′ by ez = M 1−z/zs(M) + (z2M+1), (3.3) s=1 1+z/z¯s(M)! O M ∞ M Y m′=0bm′zm′! i=0cizi!= m=0amzm!, (3.2) where zs(M), s = 1,...,M, are the roots of the numera- X X X tor, and z¯(M) is the complex conjugate of z(M). These s s where the c are known since ez = ∞ zi/i!. Multi- rootscanbefoundtoadesiredprecisionforvirtuallyany i i=0 plying out the sums on the left side of Eq. (3.2), and value of M. We have found them to 17 digit precision equating the coefficients ofz throughPz2M onboth sides, for M up to 20, a sample of which for M =1 to 5, each weobtain2M equationsin2M unknowns. ThelastM of rounded to five decimal places, is given in Table II. M s=1 2 3 4 5 1 −2.00000+i0.00000 2 −3.00000+i1.73205 −3.00000−i1.73205 3 −4.64437+i0.00000 −3.67781−i3.50876 −3.67781+i3.50876 4 −4.20758+i5.31484 −5.79242+i1.73447 −5.79242−i1.73446 −4.20758−i5.31483 5 −4.64935+i7.14205 −6.70391+i3.48532 −7.29348+i0.00000 −6.70391−i3.48532 −4.64935−i7.14205 TABLE II: The roots zs(M) of the numerator of the Padé approximant of theexponential function for M from 1 to 5. WeusethePadéapproximanttoexpressthetime evo- AssumingthatΨ isknown,wedetermineΨ from n n+1/M lution operator. Define the operator Eq. (3.7) which has a form similar to that of Eq. (2.5). We use therefore the same method of Sec. II to ob- iH∆t/¯h 1 tain Ψn+1/M. This is repeated to obtain in succession K(M) − zs(M) , (3.4) Ψn+2/M,Ψn+3/M,...,Ψn+(M−1)/M,Ψn+1. Since the op- s ≡ 1+ iH∆t/¯h eratorsKs(M) commute,theycanbeappliedinanyorder. z¯(M) s so that M IV. DISCUSSION OF ERRORS AND e−iH∆t/¯h = Ks(M)+O((∆t)2M+1). (3.5) BOUNDARY CONDITIONS s=1 Y A. Errors iH∆t/¯h Since Ψ =e Ψ , we write the relation n+1 − n Inthissectionwediscusstheerrorsasafunctionofthe M Ψ = K(M)Ψ . (3.6) orders of the method, i.e., r and M. Let us separate the n+1 s n truncation errors due to the integration over space and s=1 Y those due to integrationovertime. At a giventime t the Defining Ψ K(M)Ψ , we can solve for spatial integration with the rth-order expansion yields a n+s/M s n+(s−1)/M ≡ truncation error Ψ recursively, starting with n+1 Ψ =K(M)Ψ . (3.7) e(r) =C(r)(∆x)2r, (4.1) n+1/M 1 n 5 where C(r) is assumed to be slowly varying with r. Ac- tuallyC(r) = ψ(2r)(x∗,t)/(2r!)forsomex∗ intherange 30 z(M) | | avge of spatial integration, and thus is model dependent. If z(M) we specify an acceptable error, the step size ∆x can be 25 min z(M) adjusted to obtain that error. Since ∆x = (x0 xJ)/J, max − an adjustment of ∆x is equivalent to a change in J. Re- 20 calling that x x is fixed, we obtain 0− J |z(M)| 15 x x e(r) 1/2r 0 J ∆x= − = , (4.2) 10 J C(r) (cid:18) (cid:19) 5 and constant 0 e(r) . (4.3) 0 5 10 15 20 ≈ J2r M We have assumed that C(r) is approximately constant. FIG.1: Theaverage,theminimum andthemaximumvalues The CPU time for the calculation is proportional to the number of basic computer operations in solving the ma- of {|zs(M)|,s=1...M}, as a function of M. trix equation (2.12). This involves elementary row oper- ationsonr 1rowsinJ 1columnstobringthematrix − − to upper triangular form, plus J back substitutions to obtain the solution. Hence r,M increases from 1 through 5 or larger. For increas- r ing M the CPU time continues to decline although the CPU time # operations Jr . (4.4) ∝ ∝ ∝ (e(r))1/2r decrements become smaller at larger M. For increasing r there is a minimum depending on the specified error This form gives a minimum (optimum) CPU time which and beyond the minimum the curve shows a slow in- occurs when crease with increasing r. Superimposed on the curves are the CPU times (as dots) of a model calculation (see lne(r) r . (4.5) Sec.VA),inwhichthenumericalandexactsolutionscan ≈− 2 be compared. Clearly the theoretical trends, including the minimum as a function of r, occur in the computed For the time integration we assume a truncation error example. It should be noted that it “pays” to increase independentofr. Foragivenr the errorduetofinite∆t has a first term in the expansion e(M) =C(M)(∆t)2M+1, (4.6) 50 whereagainC(M) isassumedtobeaslowlyvaryingfunc- e(r)(theory) tion of M. We note that the factor 1 in the numerator e(M)(theory) 2 40 and denominator of Eq. (2.4) is replaced by 1/z(M) in e(r)(calculated) s each of the M factors (3.4) of Eq. (3.5). As M increases e(M)(calculated) the averageover different values of s of zs(M) , which we me 30 denote as z(M), also increases. In fact|z(M)|is a linear U ti avge avge P function of M as is seen in Fig. 1. The effective expan- C 20 (M) sionparametercanbeapproximatedby2∆t/z rather avge than∆tandhenceisproportionalto∆t/M. Thuswecan 10 replace the relation of Eq. (4.6) by e(M) C(M)(∆t/M)2M+1, (4.7) 0 ≈ 2 4 6 8 10 12 14 16 18 20 wheretheconstantC(M) isappropriatelyadjusted. Ifwe r, M take the total time t =N∆t to be fixed, then max FIG.2: ThenormalizedtheoreticalvariationoftheCPUtime foragivenerrorof1.0×10−8. Thecalculated CPUtimesfor 1 theexample of Sec. VA are shown as dots. − CPU time e(M) 2M +1 . (4.8) ∝ (cid:16) (cid:17) In Fig. 2 the curves of the (scaled) CPU times are plot- M indefinitely, whereas there is an optimum value of r ted. Both curves clearly show the sharp decline when which depends on the magnitude of e(r). 6 B. Boundary conditions V. EXAMPLES We consider three systems to which this numerical Below Eq.(2.9) we indicate thatwe setψ =0 when j,n method may be applied. The first two allow us to make j < 0 or j > J, or when j goes out of range. These are acomparisonwiththeexactsolutionandtotestthepre- appropriateboundaryconditionswhenthewavefunction cisionofthenumericalprocedure. Thethirdinvolvesthe and its first r+1 derivatives are zero at the boundaries, time evolution of a quasi-stable quantum process. since it was assumed in the derivation of the method that all these derivatives exist. If however that is not the case, for instance, at the boundary of an (in)finite A. Oscillation of a coherent wavepacket square well or barrier where the second-order derivative doesnotexist,onemustdevisewaysofincorporatingthe The oscillationof a coherent state in the harmonic os- proper boundary conditions. cillator well is described in Ref. [23]. The time evolution One case of importance, which we discuss in the third of such states has been discussed recently in connection example (see Sec. VC) of the paper, is the case of radial with the quantum abacus. (See Ref. [6].) In that case behavior of a partial wave when angular momentum de- there is a point interactionat the center of the oscillator compositionhasbeen done. IntheS-wavecasethe wave well. It is of interest to consider narrow but finite-range function,definedonlyfornonnegativevaluesoftheradial interactions to simulate more realistic physical systems. coordinate, is zero at the origin but the first derivative In orderto test the robustness of the quantum gates one is finite. Muller [22] discusses a related, but not identi- needs a very stable numerical procedure. We test the cal, situation. He considers three-point formulas for the precisionofthenumericalprocedurebyinvestigatingthe Coulomb potential which lead to a radial wave function case without the central interaction, so that the numeri- whichiszerowhentheradialvariableρ=0,buthasfirst cal results can be compared with the exact ones. and second derivatives which are nonzero at ρ = 0. His The time-dependent Schrödinger equation is approach can be adapted to the (2r+1)-point formula of this work. i¯h∂ ψ(x,t)= ¯h2 ∂2 + 1Kx2 ψ(x,t). (5.1) ∂t −2m∂x2 2 We treatthiscasebymakingthe ansatzthatthe wave (cid:18) (cid:19) function behaves like an odd function about the origin We consider the time evolution of the initially displaced and continues in the unphysical region of the negative ground-state wave function radial variable. With this assumption we do not affect thebehaviorofthesystematpositivevaluesoftheradial ψ(x,0)= α1/2e−12α2(x−a)2, (5.2) variable, but all the required derivatives exist. Further- π1/4 more,thewavefunctionatnegativej valuescanbecom- binedwiththeoneswithcorrespondingpositivej values, where α4 = mK/¯h2, ω = K/m, and a is the initial so that the space need not be enlarged but instead the displacement. Theclosedexpressionforthetimeevolved p firstfewmatrixelementsofthematrixAcanbechanged wave function is to account for the boundary condition. This is achieved by replacing A by A′ =A B in Eq. (2.13) where ψ = α1/2 exp 1(ξ ξ cosωt)2 − exact π1/4 −2 − 0 (cid:20) 1 1 i( +ξξ sinωt ξ2sin2ωt) , (5.3) 0 a1 a2 a3 ar−2 ar−1 ar 0 − 2 0 − 4 0 ··· ··· (cid:21) 0 a a a a a 0 0  2 3 4 ··· r−1 r ··· where ξ =αx and ξ =αa. We set ¯h=m=1, ω =0.2, 0 a a a a 0 0 0 0  ... ...3 ...4 ...5 ··· ...r ... ... ... ··· a[n4d0a,4=0]1.0T.hWeepcehriooodseoofuorscsipllaacteiosnucihs tthhaentxT∈=[x100,πx.JW]=e B =0 ar−1 ar 0 0 0 0 0 . a−llow the coherent state to oscillate for eleven periods  ··· ··· beforecomparingthenumericalsolutiontotheexactone. 0 a 0 0 0 0 0 0   r ··· ··· The error is calculated as e using the formula [19] 0 0 0 0 0 0 0 0  2  ··· ···  ... ... ... ... ... ... ... ...  (e2)2 = xJ dx ψ(x,t1) ψexact(x,t1)2, (5.4)  (4.9) Zx0 | − | A hard-core type potential could be dealt with in the where t = 11T for our example. The results including 1 same way. Different forms of boundary conditions are the relative CPU time [34] are displayed in Table III. more complicated to implement, but Ref. [22] suggests In the above tests we have tried to obtain a precision an approachto including such boundary conditions. For better than 10−8. While varying the number of steps the purpose of the radial wave function of a nonsingular for the spatial integration, we kept the number of time potential, the above approachis sufficient. factors per time step constant at 20. Given that the 7 M r ∆t ∆x J e CPU time functionisusuallynotknownandanestimatesuchaswe 2 20 20 π 0.44444 180 6.717×10−9 18.38 have given here would be all that is available. The main 20 15 π 0.38095 210 7.044×10−9 14.23 pointisthatdramaticimprovementsresultboththeoret- ically and computationally when larger values of r and 20 10 π 0.27586 290 7.506×10−9 10.84 M are employed. 20 7 π 0.18182 440 9.353×10−9 10.12 20 5 π 0.09877 810 9.330×10−9 12.83 20 4 π 0.05755 1390 9.871×10−9 18.61 B. Propagation of a wavepacket 20 3 π 0.03810 2100 2.102×10−7 23.67 20 2 π 0.03810 2100 1.624×10−4 18.42 For this example we return to the workof Ref. [8] and 20 1 π 0.03810 2100 1.749×10−1 13.75 consider the main features of that analysis with a view 20 10 π 0.26667 300 5.106×10−9 10.77 of determining the improvement brought about by the 15 10 π/1.5 0.26667 300 5.153×10−9 12.13 generalizations of this paper. This problem was revis- 10 10 π/3 0.26667 300 4.995×10−9 16.16 ited by Moyer [18] to illustrate the efficacy of the Nu- 5 10 π/15 0.26667 300 8.787×10−9 40.42 merovmethodandthe useoftransparentboundarycon- 3 10 π/150 0.26667 300 1.840×10−9 242.9 ditions for the propagation of free-particle wavepackets. 1 10 π/3000 0.26667 300 5.046×10−4 1627 The authors of Ref. [8] consider wavepackets impinging onasquarebarrierandstudytheirbehaviorintime. We TABLE III: Summary of computational time and errors in- consider first free wavepacket propagation (without po- curredbyusingthenumericalintegrationprocedurewhenthe initial wave function is the displaced ground state. The last tential), and second the reflection and transmission of a column indicatesarelativeCPU runtime. Theupperhalfof wavepacketby a smooth potential. the table gives the effects of changing the number of spatial ThuswefirstassumeV(x)=0andtakeasinitialwave steps; the lower half the effects of changing the number of time steps. function ψ(x,0)=(2πσ02)−1/4eik0(x−x0)e−(x−x0)2/(2σ0)2. total space is fixed and spans 80 units, we adjusted the (5.5) number of spatial steps J to give the required precision. (Note thatourσ is thatofRef.[8]dividedby√2.) The 0 We limited (arbitrarily)the maximumnumber of spatial wave function at later time is given by stepsto2100. WithM =20the15-pointformula(r =7) ismostefficient. Whenr<4(lessthan9-pointformula), ψ(x,t)=(2πσ2)−1/4[1+i¯ht/(2mσ2)]−1/2 0 0 × we were unable to reach the precision criterion because of the imposed limit on J. It is clear from the trend however that the efficiency is significantly less for the (x x )2/(2σ )2+ik (x x ) i¯hk2t/(2m) exp − − 0 0 0 − 0 − 0 . lower r values. The 9-point formula is roughly half as 1+i¯ht/(2mσ2) (cid:26) 0 (cid:27) efficient as the 15-point formula. (5.6) The effect of different order time formulas as seen in We use parameters comparable to those of Ref. [8]. thelowerpartofTableIIIisevenmoredramatic. Forthe We set ¯h=1 and m= 1. The coordinate range we take 2 spatialintegrationweusedthe21-pointformula(r =10), is from 0.5 to 1.5 rather than from 0 to 1 since over and varied the time-order formula, i.e., M, from20 to 1. the smal−ler space the normalization of the packet is not We see at least two orders of magnitude improvementin as precise as we require because the tails of the Gaus- computational speed as M is increased over this range. sian are nonzero outside the (0,1) interval. We choose A comparison with the standard CN approach (r = σ = 1/20, k = 50π, ∆t = 2(∆x)2, and allow as much 0 0 M = 1) is instructive. We considered the same sys- time as it takes the packet to travel from x = 0.25 to 0 tem with x0 = 25 and xJ = 25, ∆x = 0.005 and around 0.75. For the final position the numerically cal- − ∆t = 0.5(∆x)2. The standard CN method yielded an culated wave function is compared to the analytic one eexrrpoornoefntei2al=ly7t.o1e×21=0−25.7whe1n0−t3=atT/t4=w1h0iTch. iWnchreearseeads aconmdpeu2teisddreetseurlmtsi.neWd.e IonbsTeravbelethIVatwtheelitsrtasdoitmioenaolf CthNe × the CPU time in Table III is given in seconds, the CPU method (M =1 and r =1) has a low precision and that time required to complete this last calculation exceeded using greater J (smaller ∆x) results in modest gain in 24 hours. precision. By using higher-order time formula one can ThecomputedCPUtimesshowninFig.2exceededthe make significant gain in precision (seven orders of mag- “theoretic” values by increasing amounts as r increased nitude)withnoincreaseincomputationaltimecompared beyond 10. This can be attributed to the approximate with that of Ref. [8]. (Compare the first row to the last nature of the error analysis in which the model depen- two rows of Table IV.) Rows 5 through 9 of Table IV dence of C(r) (and C(M)) was neglected. In this exam- illustrate the transition from less precise to more precise ple a more elaborate analysis could be done since the solutions. It is consistentwith the finding of the authors wave function is known analytically. In practical situa- of Ref. [5] who use an M =2, r = 3 method and obtain tionswhereanumericalmethodisusedtheanalyticwave two orders of magnitude improvement of the results of 8 Refs.[1, 3]. The resultsaresensitive to surprisinglyhigh 11 orders of ∆x and ∆t. Another test using wavepacket scattering to show the 00..88 efficacy of the higher-order approach is scattering from a potential. Rather than using the square barrier of Ref. [8], we consider the repulsive Po¨schl-Teller type po- k)k) 00..66 TT((kk)) R(R( RR((kk)) tential [24, 25] of the form k), k), TTwwpp((kk)) T(T( 00..44 RR ((kk)) wwpp ¯h2 β2λ(λ 1) V(x)= − . (5.7) 00..22 2m cosh2βx 00 Since this potentialdoes nothavediscontinuitiesthe im- 11 11..55 22 22..55 33 provedCN method workswell with it. The transmission kk and reflection coefficients are known analytically. We can also compute them by considering the wavepacket FIG.3: Thetransmissionandreflectioncoefficientsasafunc- tionofk whenβ =1,λ=2.5, m=1,andσ =10. Thesub- Eq. (5.5) incident on the potential. Over a sufficiently 0 script “wp” indicates that the coefficients are obtained from long time the wavepacket will have interacted with the theemergingwavepackets. Thequantitieswithoutsubscripts potential and transmitted and reflected packets emerge are calculated using thetime-independentmethod. andtravelawayfromthe potentialregion. At thatpoint we can calculate the probabilities of the particle repre- C. Long-time behavior of decay of quasi-stable sented by the packet on the left and on the right of the system potential;theseprobabilitiescorrespondtothetransmis- sionandreflectioncoefficientsprovidedthepacketissuf- There are few analytically solvable models of the time ficiently narrow in momentum space. This means that evolution of unstable quantum systems [26, 27]. Real- one needs an initial packet which is wide in coordinate istic systems need to be solved numerically. Long-time space. In our calculation we choose β = 1, λ = 2.5, calculations are required for systems which require both m = 1, and σ = 10. This gives a spread in the inci- 0 nuclear and atomic time-scales such as ionization and dent momentum-space wave packet of σ = 0.05. The k bremsstrahlung due to radioactive decay of the nucleus width in momentum space ofthe reflected andtransmit- of an atom [2, 3, 5]. To study the short-time anoma- ted wavepackets also has this value. The domain of the lous power law behavior of the decay and the long-time x coordinates is from 300 to +300 and the initial posi- − inverse-power law behavior the method of this paper is tion of the packet is at x = 150 to ensure that there 0 − appropriate. This is especially relevant because of the is no overlap of the initial packet and the potential. We recently observed violation of the exponential-decay law findgoodagreementbetweenthetransmissionandreflec- at long times [28]. tions probabilities determined by plane-wave scattering To illustrate the numerical method discussed in this approach and the time-dependent calculation as shown paper as applied to decaying systems let us consider a in Fig. 3. variantofthemodelwithaδ-shellpotential[9],butwith the δ function replaced by a gaussian. Thus in the S partial wave of a spherical system the potential is M r J e CPU time 2 λ 1 1 2000 9.418×10−2 2.20 V(ρ)= exp[ (ρ a)2/w2], (5.8) w√π − − 4000 2.189×10−2 18.04 8000 5.368×10−3 151.92 where ρ is the radial coordinate. This potential reduces 16000 1.336×10−3 1287.8 to the δ-shell potential, V (ρ) = λδ(ρ a) when w δ 2 2 2000 3.018×10−4 5.99 0. For small but finite values of w thi−s potential lea→ds 3 3 2000 1.321×10−6 12.57 to scattering results which are good approximations of 4 4 2000 6.577×10−9 22.64 thoseoftheδ-shellinteraction[35]. Initiallythequantum 5 5 2000 3.648×10−11 37.13 system is in the state 6 6 2000 8.437×10−13 56.05 10 10 440 3.606×10−9 2.12 ψ(ρ,0)= 2/asin(πρ/a). (5.9) 20 20 260 4.542×10−9 2.51 In our example we take h¯p=1, λ=3, m= 1, a=1 and 2 w = 0.10. Using the numerical method of this paper, TABLE IV: Summary of computational parameters used to calculate the propagating free packet and compare it to the including the modification of matrix A as described in analytic wavepacket. Sec. IVB to take care ofthe boundary conditions at ρ= 9 0, we determine the wave function at later times, i.e., themainpacketwhichcorrespondstotheexponentialde- ψ(ρ,t). From that we obtain the nonescape probability, cay at the resonance energy, and the follower,which is a asafunctionoft,P(t)= a ψ(ρ,t)2 dρ,whichisshown smallblipthatstretchesintime(travelsmoreslowly)and 0 | | in Fig. 4. It clearly shows the exponential decay-region is due to energy components in the initial wave function R whichhavelowerenergythattheresonanceenergy. Ifthe maximum spatial coordinate,which is set at 800 for this 0 calculation, were set at a smaller value, say 400, then -1 one observes a fuzziness in the precursor of the right- -2 most wave function. This can be attributed to the fi- 1 nite space in which the wavepacket travels so that the -3 0.8 P(t) -4 P(t) 00..46 fbaosutnpdraercyurasnodrhinatserbfeeernespwairtthiatllhyerweflaevcetefrdonfrtoomf tthheemrigahint 0 g1 -5 0.2 wavepacket. The numerical parameters for this calcula- o l 0 tion are ∆ρ=0.1, ∆t=0.02, r =20 and M =20. -6 0 0.2 0.4 0.6 0.8 1 -7 t -8 VI. REMARKS -9 0 2 4 6 8 10 12 14 ThegeneralizedCNmethodthatwehavepresentedin t this paper gives many orders of magnitude improvement in the precision of the results and several orders of mag- FIG. 4: The nonescape probability as a function of time for nitude in the computational time required to obtain the the interaction with λ = 3, a = 1, and w = 0.10. We also results. Clearly,sincethismethodenhancestheefficiency take¯h=1 and m= 1. 2 of the numerical calculations, it can be a significant tool for studying time-dependent processes. It goes beyond in time, the inverse power-law behavior for long times, the improvement of Ref. [5] in a systematic way. It also thedeviationfromexponentialdecayatshorttimes,and is anadvanceoverthe method ofRef. [18], since the Nu- the transition regions [9]. merov method has an error ((∆x)6), and it is difficult The quadratic short-time behavior is seen in Fig. 4 O to see how it can be generalizedsystematically to higher as is the inverse power law behavior at long times [29]. order spatial errors. The generalized time evolution al- Remarkably the decaying system can be studied in this gorithm can be applied to Moyer’s [18] method. manner for a time exceeding thirty half-lives. In Fig. 5 WehaveappliedthemethodforrandM uptoandin- we plot the square of the absolute value of the wave cluding20forboth. Havingachievedsignificantimprove- function at times t = 5,10,15. Notice that the wave ments with these values of r and M we did not consider larger values although there does not seem to be a practical reason that this cannot be done. Although the ap- 0.08 proach seems to saturate at r around 10 (see Table III), 22 VVVV ((((rrrr)))) thereisnovisiblesaturationinthetime-evolutionpartof 1111 11..55 the problem. One would expect higher orders of spatial 0.06 ||ψψ((rr,,00))||22 11 errors to be significant when the wave function and/or 00..55 the potential has large spatial fluctuations. Even so, at |ψ(r,t)|2 00 r = 7 we are using a 15-point formula for the second- 0.04 00 00..44 00..88 11..22 orderspatialderivative,anditissurprisingthatasmaller rr number of points in the formula is not sufficient for op- timal results in this case. It should be noted that the 0.02 higher order method as discussed in this paper are suit- able only for well-behaved, sufficiently differentiable so- 0 lutions; these occur when the potential function is well 0 20 40 60 80 100 120 140 behaved. As the authors of Ref. [10] point out for sin- r gular functions higher-order methods do not necessarily lead to greater accuracy. FIG.5: Thesquareoftheabsolutevalueofthewavefunction It should be noted that in this paper we consider pri- as a function of ρ at times t = 5,10,15 for the same parameters apply as in Fig. 4. The insert gives the t= 0 graph as marily one-dimensionalsystems, but the method applies well as thescaled potential, V1(ρ)=V(ρ)/8. equally well to partial-wave equations of two- or three- dimensional systems. The study of the decaying quasi- function (packet) has three distinct regions: a precursor stable state is an example of the latter. due to energy components of the initial wave function An interesting avenue to investigate further is the larger than that associated with the exponential decay, impact that this approach may have on two- or 10 three-dimensional systems, where the number of vari- (short-range) potential region and little variation in the ables involved is equal to the dimension. The asymptotic regions. A great savings in computational Peaceman-Rachford-type approach [30], also known as time can be achieved by using different space-step sizes the alternating-direction implicit method, of factoring inthedifferentregions. Oneneedstoinvestigatewhether the approximation of the time evolution operator may such variable step size can be incorporatedin the gener- applyasitdidinRefs.[14,31]ormorerecentlyin,forex- alized spatial integration scheme of this paper. The ap- ample, Refs. [32, 33]. Using the Crank-Nicolson method proachthatdealtwiththediscontinuousfirst-orsecond- the authors of Refs. [14, 31] show that in two dimension order derivative in this paper and Ref. [22] is worth ex- the kinetic energyparts ofthe evolutionoperatorfactor- ploring. We intend to study this in the future. izes. Whether such factorization can be generalized in the spirit of the method of this paper is under investiga- tion. Acknowledgments Calculationsonone-dimensionalmultichannelsystems indicate that this approach also leads to substantially greater efficiencies. We are grateful to Professor Y. Nogami for carefully Preliminarystudies with the Numerov spatialintegra- reading the manuscript and the constructive comments, tion scheme [18] and the generalized time evolution as as well as useful and helpful discussions. One of the au- described in this work indicate that significant improve- thors(WvD)expressesgratitudeforthehospitalityofthe ments occur if one incorporates appropriate changes in DepartmentofInformationandCommunicationSciences the spatial step size for different regions of space. This ofKyotoSangyoUniversitywhere mostofthis workwas is important in the case of discontinuous potentials and completed. He also acknowledges the financial support potentials that have great variation in some region and for this research from the Natural Sciences and Engi- little or no variation in other regions. Furthermore it is neeringCouncilofCanadaandthe JapanSocietyforthe wellknownthatthewavepackethaslargefluctuationina Promotion of Science. [1] W. van Dijk, F. Kataoka, and Y. Nogami, J. Phys. A: [20] I. Puzynin, A. Selin, and S. Vinitsky, Comp. Phys. Math. Gen. 32, 6347 (1999). Comm. 126, 158 (2000). [2] F. Kataoka, Y. Nogami, and W. van Dijk, J. Phys. A: [21] G. A. Baker, Jr. and P. Graves-Morris, Padé Approxi- Math. Gen. 33, 5547 (2000). mants: Part I: Basic Theory (Addison-Wesley Publish- [3] C. A. Bertulani, D. T. de Paula, and V. G. Zelevinsky, ing Company, Reading, Massachusetts, 1981). Phys.Rev.C 60, 031602 (1999). [22] H. G. Muller, Laser Physics 9, 138 (1999). [4] W.vanDijkandY.Nogami, Few-BodySystemsSupple- [23] L. I. Schiff, Quantum mechanics, International series in ment 14, 229 (2003). pure and applied physics (McGraw-Hill Inc., New York, [5] S¸. Mi¸sicu, M. Rizea, and W. Greiner, J. Phys. G: Nucl. 1968), 3rd ed. Part. Phys. 27, 993 (2001). [24] S. Flu¨gge, Practical Quantum Mechanics (Springer- [6] T. Cheon, I. Tsutsui, and T. Fu¨löp, Physics Letters A Verlag, New York,1974). 330, 338 (2004). [25] L. D. Landau and E. M. Lifshitz, Quantum mechanics [7] C.N.Veenstra,W.vanDijk,D.Sprung,andJ.Martorell (Nonrelativistic theory), vol. 3 of Course of theoretical (2006), e-print cond-mat/0411118. physics (Pergamon Press, Oxford, 1977), 3rd ed., trans- [8] A.Goldberg,H.M.Schey,andJ.L.Swartz,Am.J.Phys. lated from Russian byJ. B. Sykesand J. S. Bell. 35, 177 (1967). [26] W. van Dijk and Y. Nogami, Phys. Rev. C 65, 024608 [9] R.G. Winter, Phys. Rev. 123, 1503 (1961). (2002); Phys.Rev.C 70, 039901(E) (2004). [10] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and [27] W.vanDijkandY.Nogami,Phys.Rev.Lett.90,028901 B.P.Flannery,Numerical recipes inC:Theartofscien- (2003). tificcomputing(CambridgeUniversityPress,Cambridge, [28] C. Rothe, S. I. Hintschich, and A. P. Monkman, Phys. 1992), 2nd ed. Rev. Lett.96, 163601 (2006). [11] M. Patriarca, Phys.Rev.E 50, 1616 (1994). [29] L. A. Khaflin,Soviet Physics JETP 6(33), 1053 (1958). [12] X. Qian, J. Li, X. Lin, and S. Yip, Phys. Rev. B 73, [30] D. W. Peaceman and J. H. H. Rachford, J. Soc. Indust. 035408 (2006). Appl. Math. 3, 28 (1955). [13] V. M. Vyas, T. S. Raju, C. N. Kumar, and P. K. Pani- [31] K.C.Kulander,K.S.Devi,andS.E.Koonin,Phys.Rev. grahi, J. Phys.A: Math. Gen. 39, 9151 (2006). A 25, 2968 (1982). [14] I.Galbraith,Y.S.Ching,andE.Abraham,Am.J.Phys. [32] N.H.Shon,A.Suda,andK.Midorikawa,RIKENReview 52, 60 (1984). 29, 66 (2000). [15] W. van Dijk, K. Kiers, Y. Nogami, A. Platt, and [33] K. L. Ishikawa, Phys. Rev.A 70, 013412 (2004). K.Spyksma,J. Phys.A: Math. Gen. 36, 5625 (2003). [34] The CPU timeisarelative measure. Thecalculations of [16] D.KosloffandR.Kosloff,J.Comp.Phys.52,35(1983). TableIIIweredoneonacomputerwithanIntelPentium [17] C. Leforestier, et al., J. Comp. Phys. 94, 59 (1991). 4 Mobile CPU 1.90 Ghz processor. Those of Table IV were done with a AMD Opteron processor 250 running [18] C. A. Moyer, Am.J. Phys. 72, 351 (2004). at 2.4 GHz. The CPU times within a table give a rough [19] I. Puzynin, A. Selin, and S. Vinitsky, Comp. Phys. comparison of the efficiency of the method. CPU times Comm. 123, 1 (1999). in different tables should not becompared.