Numerical solutions of the time-dependent Schrodinger equation in two dimensions

5.9 MB·

by  Wytse van Dijk

Numericalsolutionsofthetime-dependentSchro¨dingerequationintwodimensions
Wytse van Dijk,1,2 Trevor Vanderwoerd,1 and Sjirk-Jan Prins1
1DepartmentofPhysics, RedeemerUniversityCollege, Ancaster, Ontario L9K1J4, Canada
2DepartmentofPhysicsandAstronomy,McMasterUniversity,Hamilton,OntarioL8S4M1,Canada∗
(Dated:January30,2017)
ThegeneralizedCrank-Nicolsonmethodisemployedtoobtainnumericalsolutionsofthetwo-dimensional
time-dependentSchro¨dingerequation. Anadaptedalternating-directionimplicitmethodisused,alongwitha
high-orderfinitedifferenceschemeinspace. Extracarehastobetakenfortheneededprecisionofthetime
development. Themethodpermitsasystematicstudyoftheaccuracyandefficiencyintermsofpowersofthe
spatialandtemporalstepsizes.Toillustrateitsutilitythemethodisappliedtoseveraltwo-dimensionalsystems.
7
1
0
I. INTRODUCTION cision rapidly with increasing M. Besides these three ap-
2
proachesthereareotherapproximationsofthetime-evolution

n
operator, e.g., the exponential split-operator method [13] or
The determination of accurate numerical solutions of the
a
the iterative Lanczos reduction [14]. Like time development
J time-dependent Schro¨dinger equation is an ongoing enter-
7 prise. Thequantumwaveequationisfundamentaltotheun- thespatialintegrationcanbeachievedindifferentways,e.g.,
bydifferenttypesoffinitedifferencingorbythepseudospec-
2 derstanding of nonrelativistic atomic and subatomic systems
andphenomena.Consequentlyitoccursinadiversityofphys- tralfastFouriertransformapproach.

] icalsystems. Ideallyanalyticsolutionsareavailable,butmost Sincemany ofthecalculations referredtohave beendone
h
inonespatialdimension,inthispaperweconsiderthegener-
realisticsituationsaretoocomplextoyieldsuchsolutions.
p
alized Crank-Nicolson with two spatial dimensions. A num-
- Inthelastfewyearsanumberofimprovementshavebeen
p ber of articles have appeared recently that describe methods
madetoyieldmoreaccuratesolutionswithgreaterefficiency.
m of solving the two-dimensional time-dependent Schro¨dinger
Thetypeofmethodoftendependsontheproblemathand,i.e.,
o dimensionality, time dependence of the interaction, short- or equation,includingthosewithtime-dependentpotentialsand
c long-time behaviour, etc. The “method of choice” for some nonlinearterms. See,forexample,Refs.[15–20]. Anumber
. oftheseusetheCayleyformforthetimeevolutionoperator.
s years is the Chebyshev polynomial expansion of the time-
c We wish to employ the higher order Pade´ form in order to
evolution operator with (inverse) Fourier transformations to
i enhance the efficiency of the approach. Given the two spa-
s deal with the spatial development as time progresses [1, 2].
y tial dimensions, we pursue an alternating-direction implicit
More recently the Pade´ approximant representation of the
h schemewhichrequiresonlysolvingone-dimensionalimplicit
time-evolution operator is exploited [3–7]. This approach is
p
problemsforeachtimestep. Differentapproacheshavebeen
[ unitary,stable,andallowsforsystematicestimateoferrorsin
terms of powers of the temporal and spatial step sizes. The suggested, such as the use of multigrid partitioning [21], but
1 it is our intention to present one that provides the user with
two approaches have been shown to have comparable effi-
v another efficient alternative. Clearly the method chosen will
cacy [5, 8]. Gusev et al. [9–11] have recently given an im-
7 dependonthecontext.
proved and extended application of the method discussed by
3
InsectionIIweformulatethetimedependenceoftheprob-
1 Puzyninetal.[3]. Theydealwiththemoregeneralproblem
lem. Section IIIis adescriptionof thespatial integration. A
8 of a time-dependent Hamiltonian. Using a truncated Mag-
0 nus expansion with additional transformations, they are able number of applications are discussed in Sec. IV, and Sec. V
1. toobtainstableandefficientsolutionswhichareaccurateup presentsconclusionsandadiscussionofthework.
0 tosixth-orderinthetimestep.
7 Generally the various approaches involve time evolution
II. ACCURATETIME-EVOLUTIONSCHEME
1 andintegrationoverspace. Thusthereareanumberofways
:
v of dealing with the time evolution. Crank-Nicolson approx-
Wesolvethetwo-dimensionaltime-dependentSchro¨dinger
i imates the exponential time-evolution operator by a Cayley
X equation
form which retains unitarity, but is correct only to low order
r in time advance [12]. The Chebyshev polynomial expansion (cid:18) (cid:19)
a ∂
canleadtohighaccuracyevenoversignificanttimeintervals. H(cid:98) −i(cid:126) Ψ(x,y,t)=0, (2.1)
∂t
It is not explicitly unitary. The generalized Crank-Nicolson
approximatestheevolutionoperatorwitha[M/M]Pade´ ap- where
proximant, factorized into M factors of Cayley form. This
(cid:126)2 ∂2 (cid:126)2 ∂2
form is unitary and has a truncation error of O[(∆t)2M+1], H(cid:98) =− − +V(cid:98)(x,y)
where ∆t is the temporal step size. This improves the pre- 2m∂x2 2m∂y2 (2.2)
=K(cid:98)x+K(cid:98)y+V(cid:98)(x,y),
startingwithaninitialwavefunction
∗[email protected] Ψ(x,y,0)=Φ(x,y). (2.3)
2
Thetime-evolutionoperatorofthesystemgivesanexpres- DefiningΨ(n+s/M) ≡ Θ(cid:98)(sM)Ψ(n+(s−1)/M),wecansolvefor
sionforthewavefunctionatatimeintermsofthewavefunc- Ψ(n+1)iterativelystartingwithΨ(n+1/M) =Θ(cid:98)(M)Ψ(n),then
1
tionatanearliertime,i.e., Ψ(n+2/M) =Θ(cid:98)(M)Ψ(n+1/M),andsoon.
2
Ψ(t+∆t)=e−iH(cid:98)∆t/(cid:126)Ψ(t), (2.4) Let us start with the basic substep of the procedure in go-
ing from Ψ(n+(s−1)/M) to Ψ(n+s/M), which we label below
genericallyasΨ0andΨ+,respectively. Wethenwrite
where∆tisthetimeadvance,andwherewehavesuppressed
thespatialcoordinatesxandy inthewavefunction. Wewill
employ the factorized [M/M] Pade´ approximant along with (cid:32) (cid:33)
1+(iH(cid:98)∆t/(cid:126))/z
the alternating-direction implicit method [22]. In keeping Ψ+ = Ψ0, (2.8)
with the expansion of the time-evolution operator discussed 1−(iH(cid:98)∆t/(cid:126))/z¯
inRef.[4],theoperatoriswrittenas
or
M
e−iH(cid:98)∆t/(cid:126) = (cid:89)Θ(cid:98)(M)+O[(∆t)2M+1], (2.5)
s
(cid:16) (cid:17) (cid:16) (cid:17)
s=1 1−(iH(cid:98)∆t/(cid:126))/z¯ Ψ+ = 1+(iH(cid:98)∆t/(cid:126))/z Ψ0, (2.9)
where
Θ(cid:98)(M) ≡ 1+(iH(cid:98)∆t/(cid:126))/zs(M), (2.6) wherez isthegenericzs(M). SinceH(cid:98) = K(cid:98)x +K(cid:98)y +V(cid:98),we
s 1−(iH(cid:98)∆t/(cid:126))/z¯s(M) write
andz(M),s = 1,...,M aretherootsofthenumeratorofthe (cid:20) ∆t(cid:21)
[M/Ms]Pade´approximantofez;thez¯(M)arethecorrespond- 1−i(K(cid:98)x+K(cid:98)y+V(cid:98))(cid:126)z¯ Ψ+
s
(2.10)
ing complex conjugates. Since Ψ(n+1) = e−iH(cid:98)∆t/(cid:126)Ψ(n) (n (cid:20) ∆t(cid:21)
referstothetimetn =n∆t,n=0,1,...),wewrite = 1+i(K(cid:98)x+K(cid:98)y+V(cid:98))(cid:126)z Ψ0
M
(cid:89)
Ψ(n+1) = Θ(cid:98)(M)Ψ(n). (2.7) sothat
s
s=1
(cid:18) ∆t(cid:19)(cid:18) ∆t(cid:19) (∆t)2
1−iK(cid:98)x(cid:126)z¯ 1−iK(cid:98)y(cid:126)z¯ Ψ++K(cid:98)xK(cid:98)y (cid:126)2z¯2 Ψ+
(2.11)
(cid:18) ∆t(cid:19)(cid:18) ∆t(cid:19) (∆t)2 ∆t ∆t
= 1+iK(cid:98)x(cid:126)z¯ 1+iK(cid:98)y(cid:126)z¯ Ψ0+K(cid:98)xK(cid:98)y (cid:126)2z¯2 Ψ0+iV(cid:98) (cid:126)z¯Ψ++iV(cid:98) (cid:126)zΨ0.
InkeepingwithPeacemanandRachford[22],wedefineΨ(cid:101) bytheequation
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
∆t ∆t ∆t ∆t
1−iK(cid:98)x(cid:126)z¯ Ψ(cid:101) = 1+iK(cid:98)y(cid:126)z Ψ0+i 1−iK(cid:98)x(cid:126)z¯ V(cid:98) (cid:126)zΨ0. (2.12)
WeinsertthisexpressionintoEq.(2.11)toobtain
(cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
∆t ∆t ∆t ∆t ∆t ∆t ∆t
1−iK(cid:98)x(cid:126)z¯ 1−iK(cid:98)y(cid:126)z¯ Ψ+ = 1+iK(cid:98)x(cid:126)z 1−iK(cid:98)x(cid:126)z¯ Ψ(cid:101) − 1+iK(cid:98)x(cid:126)z 1−iK(cid:98)x(cid:126)z¯ iV(cid:98) (cid:126)zΨ0
(2.13)
(∆t)2 (∆t)2 ∆t ∆t
−K(cid:98)xK(cid:98)y (cid:126)2z¯2 Ψ++K(cid:98)xK(cid:98)y (cid:126)2z2 Ψ0+iV (cid:126)z¯Ψ++iV(cid:98) (cid:126)zΨ0.
(cid:18) (cid:19)
∆t
OperatingonEq.(2.13)withtheinverseof 1−iK(cid:98)x(cid:126)z¯ ,weget
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
∆t ∆t ∆t ∆t
1−iK(cid:98)y(cid:126)z¯ Ψ+ = 1+iK(cid:98)x(cid:126)z Ψ(cid:101) − 1+iK(cid:98)x(cid:126)z iV(cid:98) (cid:126)zΨ0
(cid:18) ∆t(cid:19)−1 (∆t)2 (cid:18) ∆t(cid:19)−1 (∆t)2
− 1−iK(cid:98)x(cid:126)z¯ K(cid:98)xK(cid:98)y (cid:126)2z¯2 Ψ++ 1−iK(cid:98)x(cid:126)z¯ K(cid:98)xK(cid:98)y (cid:126)2z2 Ψ0 (2.14)
(cid:18) ∆t(cid:19)−1 ∆t (cid:18) ∆t(cid:19)−1 ∆t
+ 1−iK(cid:98)x(cid:126)z¯ iV(cid:98) (cid:126)z¯Ψ++ 1−iK(cid:98)x(cid:126)z¯ iV(cid:98) (cid:126)zΨ0.
3
The inverse operators are expanded, but we must make sure that the expansions are correct to O[(∆t)2M] since the overall
expansion(2.5)isofthatorder.
Simplifyingandkeepingtermsupto(∆t)2M andassumingM >1,weobtain
(cid:18) ∆t ∆t(cid:19) (cid:18) ∆t(cid:19) (∆t)2
1−iKy(cid:126)z¯ −iV (cid:126)z¯ Ψ+ = 1+iKx(cid:126)z Ψ(cid:101) + (cid:126)2z2 KxVΨ0
2(M(cid:88)−1)(cid:18) ∆t(cid:19)l (∆t)2 2(M(cid:88)−1)(cid:18) ∆t(cid:19)l (∆t)2
− iKx(cid:126)z¯ KxKy (cid:126)2z¯2 Ψ++ iKx(cid:126)z¯ KxKy (cid:126)2z2 Ψ0 (2.15)
l=0 l=0
+2M(cid:88)−1(cid:18)iK ∆t(cid:19)liV ∆tΨ++2M(cid:88)−1(cid:18)iK ∆t(cid:19)liV ∆tΨ0+O(cid:2)(∆t)2M+1(cid:3).
x(cid:126)z¯ (cid:126)z¯ x(cid:126)z¯ (cid:126)z
l=1 l=1
TheM =1case,forwhichz =−2,resultsintheequation
(cid:18) (cid:19) (cid:18) (cid:19)
∆t ∆t ∆t
1+iK(cid:98)y 2 +iV(cid:98) 2 Ψ+ = 1−iK(cid:98)x 2 Ψ(cid:101). (2.16)
This equation is a typical implicit equation with the Cayley form. We solve Eq. (2.15) iteratively so that Ψσ+1 → Ψ+ as
σ =0,1,2,... increasesintheequation
(cid:18) ∆t ∆t(cid:19) (cid:18) ∆t(cid:19) (∆t)2
1−iK(cid:98)y(cid:126)z¯ −iV(cid:98) (cid:126)z¯ Ψσ+1 = 1+iK(cid:98)x(cid:126)z Ψ(cid:101) + (cid:126)2z2 K(cid:98)xV(cid:98)Ψ0
2(M(cid:88)−1)(cid:18) ∆t(cid:19)l (∆t)2 2(M(cid:88)−1)(cid:18) ∆t(cid:19)l (∆t)2
− iK(cid:98)x(cid:126)z¯ K(cid:98)xK(cid:98)y (cid:126)2z¯2 Ψσ+ iK(cid:98)x(cid:126)z¯ K(cid:98)xK(cid:98)y (cid:126)2z2 Ψ0 (2.17)
l=0 l=0
2M(cid:88)−1(cid:18) ∆t(cid:19)l ∆t 2M(cid:88)−1(cid:18) ∆t(cid:19)l ∆t
+ iK(cid:98)x(cid:126)z¯ iV(cid:98) (cid:126)z¯Ψσ+ iK(cid:98)x(cid:126)z¯ iV(cid:98) (cid:126)zΨ0.
l=1 l=1
WestarttheiterationwithsettingΨσ=0 = Ψ0. WhenΨσ+1 orderdifferencingmethodisinprinciplefasterthanaspectral
and Ψσ are sufficiently close we stop. Note that we need to methodsinceitscalesasthe’bandedness’timesthesizeofthe
calculateΨ(cid:101) onlyonceforeachsequenceofiterations.Wefind grid,O(bN),ratherthanasO(NlogN)”. Inthecaseoftwo-
thatthisapproachcangiveaccurateresults; typicallyaround dimensional systems using the alternating-direction implicit
sixiterationsarerequiredforpreciseresults. Thisprocesshas approachN isreplacedbyN2,whereasbisunchanged. For
toberepeatedforeachoftheM stepsneededtoachieveafull the purpose of this work we therefore use finite differences.
timestepadvance. One could choose the traditional three-point expression for
There is an alternative approach to solving Eq. (2.15) for the second-order partial derivative. There are however more
Ψ+. ThetermsontherightsideinvolvingΨ+ canbemoved precisemethods. ForinstancetherecentNumerovrecentap-
to the left side and one solves a linear system of equations proach[7]givesmuchhigheraccuracy,asdoesthehigh-order
upon the discretization of the spatial variables. However, as compactfinitedifferenceapproachinRefs.[15,16]. Thetra-
weshowinthenextsection, thekineticenergyoperatorsare ditional approach is O(h2), where h is the spatial step size,
bandeddiagonalmatrices,andthoseoperatorsraisedtosome whereasthehigh-ordercompactmethodisO(h4),andtheNu-
powerwouldresultinmatriceswiththesizeofthebandsin- merovalgorithmisO(h5).Theadvantageoftheseapproaches
creased. Asaresultthegainsinefficiencyofabandedmatrix isthattheyleadtothree-pointformulaswhichmaybeconve-
formulationarelost. nientwhencrossingadiscontinuityofthepotentialorconsid-
eringanadaptivespatialgrid[7].
As in earlier work [4] we consider formulas which allows
onetochooseanarbitraryorderofh.Foraspatialgrid(inone
III. SPATIALINTEGRATION
dimension)withstepsizeh, thesecondderivativeoff(x)is
expandedas
Thenumericalspatialintegrationofthepartialdifferential
equation (2.17) can be done in a number of ways. Two ap- k=r
proaches often considered are the spectral decomposition of f(cid:48)(cid:48)(x)= 1 (cid:88) c(r)f(x+kh)+O(h2r), (3.1)
h2 k
the spatial (kinetic energy) operator or the finite-difference
k=−r
representation of this operator. The relative merits are dis-
cussedbytheauthorsofRef.[23]. Theypointoutthata“low- wherethec(r)arerealconstants,obtainedfrommakingseries
k
4
expansions of the functions f(x±kh). A similar technique A. Errors
isusedbyWangandShaoforthekineticenergyoperatoract-
ingonthewavefunctionofatwo-dimensionalstationarystate The truncation error of the series expansion of the wave
problem[24]. Inanotherarticlethesameauthorssuggestan functionintimeandspacecanbeexpressedas
expansionoftheform[25]
e=e(M)+e(r) =C (∆t)2M+1+C h2r, (4.1)
k=r 1 2
f(cid:48)(cid:48)(x)= (cid:88) a(r)f(cid:48)(cid:48)(x+kh)
k where we are considering the error of the real quantity
k=−r |Ψ(x,y,t)| and C and C are real positive numbers related
(k(cid:54)=0) (3.2) 1 2
tothe(2M +1)thpartialderivativewithrespecttotandthe
k=r
+ 1 (cid:88) b(r)f(x+kh)+O(h4r). (2r)thderivativewithrespecttoxory,respectively. Forsim-
h2 k plicity we assume h = h = h and a xy symmetry of the
x y
k=−r
wavefunction. Ifexactanalyticsolutionsareavailabletheer-
In one dimension the discretized kinetic energy is expressed ror can be calculated by comparison. If that is not the case,
as a banded diagonal matrix with bandwidth of 2r +1, just a good estimate of the error can be made by comparing the
likeinthecaseofEq.(3.1). Thusitseemsthatwithvirtually solutionforparticularM andrtotheoneobtainedwhenone
the same effort the calculation gives much more accurate re- orbothoftheM andrareincreasedbyunity[5].
sults. Acomparisonofthetwoexpansions[5]showsthatfor Tomakecomparisonsofthenumericallyobtainedsolutions
smallervaluesofr thecalculationisindeedmuchmoreeffi- to analytic solutions in cases where the latter are known, we
cient,howeverforlargerrtheaccuracydecreases.Thekinetic definetheerrore2suchthat
energyoperatorresultingfromEq.(3.1)canbemadestrictly
(cid:90) xJ(cid:90) yM
diagonallydominant,whereasthediagonaldominanceofthe e2 = |Ψ(x,y,T)−Ψ (x,y,T)|2 dydx. (4.2)
2 exact
kinetic energy matrix from Eq. (3.2) becomes compromised x0 y0
whenrgoesbeyondten. Inthispaperweuseexpansion(3.1)
Theerrore , whichisaEuclidean/(cid:96) -vectornorm, isamea-
forthekineticenergyoperatorsK(cid:98)xandK(cid:98)y. sureofthea2ccuracyofthewavefunc2tionanditsphase. Alter-
We consider a rectangular domain in space [x0,xJ] × nativelysomeauthorshaveusedthe(cid:96) -vectornorm
[y0,yM] ⊂ R2, which we partition uniformly in each direc- ∞
tion,sothatwithhx =(xJ−x0)/Jandhy =(yM−y0)/M, e =max|Ψ −Ψ (x ,y ,t)| for t=T. (4.3)
xj = x0 + jhx,j = 0,...,J and ym = y0 + mhy,m = ∞ j,m j,m exact j m
0,...,M. The time is also partitioned over the time interval
In the case that no exact solution is available, one can make
from0toT intoN subintervals, sothat∆t = T/N andthe
anestimateoftheerrorbycomparingasolutionobtainedwith
intermediate times are t = n∆t, where n = 0,1,...,N.
n
particular values of M and r to the solution obtained with
The equations we need to solve are typically of the type
M +1andr+1,e.g.,
Eqs. (2.12) and (2.17). If we let Ψ(x ,y ) = Ψ , then
j m j,m
(V(cid:98)Ψ)j,m =Vj,mΨj,mand (cid:90) x(cid:90)J yM
η2 = |Ψ(M,r)(x,y,T)−Ψ(M+1,r+1)(x,y,T)|2 dydx.
2
(K(cid:98)xΨ)j,m =(cid:18)−2(cid:126)m2 (cid:19)h12 k(cid:88)=r c(kr)Ψj+k,m (3.3) Sinceitxt0uryn0soutthate andη areverysimilar,evaluatin(g4.η4)
x k=−r 2 2 2
providesamethodtoestimatetheaccuracyintheabsenceof
for 0≤j+k ≤J.
ananalyticsolution[5,7].
Since in our applications the wave functions are zero near
ThereisasimilarrelationforK(cid:98)yΨexceptthatthesummation theboundaryofthedomain,wecanusethesimplerectangle
is over the second index of Ψ . Thus in Eq. (2.17), for
j,m+k rule for integration. The corrections to higher order polyno-
example,therightsideiscompletelyspecified,buttheΨσ+1
mialapproximationsareallintermsofevaluationsoftheinte-
on the left side needs to be found. This equation is really a
grandneartheendpoints,butsincethewavefunctioniszero
linearsystemofequationswithabandeddiagonalcoefficient
there,onegetsveryaccurateintegralswiththesimplequadra-
matrixovertheindexm. Itcanbesolvedforeachj toobtain
ture[26].
Ψσ+1. It is the strength of the alternating-direction implicit
schemethatcalculationsarereducedtoone-dimensionalones.
B. Example1:solvabletwo-dimensionalpotential
IV. IMPLEMENTATION
Considerthepotential
In this section we consider four examples in which the (cid:126)2
method outlined previously is applied in order to investigate V(x,y)=− (3−2tanh2x−2tanh2y) (4.5)
2m
its accuracy and efficiency. We will also demonstrate the
feasibility of calculating wave functions with more complex with(cid:126)=2m=1,seeFig.1. Thispotentialhasbeenusedby
structure(severalpeaksandvalleys)astheyevolveintime. severalauthorsasonewhichtestsnumericalmethods[15,16,
5
10-2 M = 1
M = 2
10-3 M = 3
M = 4
2 10-4
1
V(x,y) - 10 e2 1100--65
-2
-3 10-7
10
-5 x 0 5 -5 0 y5 10-8
10 -10 10-9
0 5 10 15 20
r
FIG. 1. (Colour online) The potential function of example 1, FIG.2. (Colouronline)Theerrorse ofthenumericalwavefunc-
2
Eq.(4.5). tionsforpotential(4.2).TheparametersareT =1andN =100.
27]. We choose the domain [−20,20]×[−20,20] and solve
In order to determine the efficiency of the calculation we
theproblemfromt=0tot=T =1withthesolutiononthe
obtaintheCPUtimewhentheerroriscloseto,butlessthan,
boundaryequaltozero[27]. Theexactsolutionis
10−6. For a particular M we adjust the spatial step size by
it choosing r until we reach a minimum CPU time. Similarly
ie
Ψ(x,y,t)= , (4.6) foraparticularr, wechoseM toyieldminimumCPUtime.
2coshxcoshy
The results are listed in Tables I and II. In Ref. [4] we give
which is also used to determine the initial wave function. It an estimate for the CPU time as a function of r for the one-
shouldbenotedthatwavefunction(4.6)issquareintegrable dimensionalcalculation.Intwo-dimensionsweexpectthebe-
and is an energy eigenstate with energy −(cid:126)2/(2m). It de- haviour to be similar since errors in x and y integrations are
scribes a bound state at threshold; the energy spectrum at similarandadditive,especiallyinsymmetriccases. Thus,
higher energies is a continuum and the corresponding wave
functionsareunbound. (cid:16) (cid:17)−1/(2r)
In our numerical calculation we allow M = 1,...,6 and CPUtime∝ e(r) r. (4.7)
r = 1,...,20, with∆t = 0.01(or100timesteps)andwith
J = M = 200. The e are plotted as a function of r for
2 ToobtainaformulafortheCPUtimeasafunctionofM when
variousvaluesofM inFig.2.
r andtheerrorareconstant, wetaketherelationoftheerror
We compare this calculation to that of Ref. [27] since the
andM tobe
othertwocalculations[15,16]aredoneoveramuchsmaller
spatialdomainwithmuchsmallertimeintervals.Inourcalcu-
lationwefindgenerallythate (cid:38)e . InRef.[27]thequoted e(M) ∝(∆t/Mν)2M+1, (4.8)
2 ∞
errorsareapproximatelye ≈10−4. Figure2showsthatthe
∞
errorinthecalculationisreducedsignificantlywhenonegoes where ν is a number less than unity. In the one-dimensional
fromM = 1toM = 2and3. WhenM > 3theresultsare case we showed that ν ≈ 1, but Table II shows that a better
identical to those of M = 3. Given that earlier calculations relationship, especiallyforlargerM, hasν < 1. SinceT =
referred to are basically M = 1 calculations with spatial er- N∆tisaconstantwewrite
rorsoftheorderofh3orh5,thismethodresultsinsignificant
improvementinaccuracy.
(cid:16) (cid:17)−1/(2M+1)
The errors of this example for M ≥ 3 saturate at e ≈ CPUtime∝NM2 ∝ e(M) M2+ν. (4.9)
2
1.5×10−9. Inthe“gullies”ofpotential(4.5)themagnitude
ofthewavefunctionislargerthanelsewhere. Inthegulliesat
the boundary of the computational space it is approximately The CPU times as a function of r and M are shown in
10−9. Thenumericalcalculationassumesthatthewavefunc- Fig.3.Theestimatesoftheerrorsareshownassolidlines.For
tion is zero outside the computational domain. The discrep- theCPUtimeasafunctionofM wehaveestimatedν =1/2.
ancybetweenthenumericalwavefunctionandtheexactone Such an estimate seems reasonable in light of the fact that
outside the computational space is the source of the residual theiterativepartoftheprocedureincreasesthetime, andthe
error. numberofiterationsvarywiththevalueofM.
6
TABLEI.SummaryofCPUtimeasrisvariedassumingaconstant
errorforexample1. ThefixedparametersareM = 5, N = 100, 1000
and∆t=0.01. CPU(r)(numerical)
r J h CPU(s) e2(10−7) 800 CPCUP(MU)((rn)(ufmoremriuclaal))
3 290 0.1379 897 8.27
CPU(M)(formula)
4 200 0.1000 396 8.94 e (s) 600
5 165 0.2424 317 9.97 m
6 150 0.2667 302 9.27 PU ti 400
C
7 140 0.2587 306 9.63
8 135 0.2963 317 8.39 200
9 130 0.3077 331 9.96
10 130 0.3077 364 7.59 0
0 5 10 15 20
11 125 0.3200 353 8.27
r, M
12 123 0.3253 392 6.37
13 123 0.3253 424 4.98
14 123 0.3253 457 3.94
FIG.3.(Colouronline)CPUtimeasfunctionofrandM forcalcu-
15 122 0.3279 466 9.75
lationwithpotential(4.5).
16 122 0.3279 535 7.78
17 122 0.3279 577 7.14
18 120 0.3333 581 9.33 where
19 120 0.3333 629 8.62
(cid:18) α2β (cid:19)1/2 e−i(n+1/2)θ
20 120 0.3333 693 7.98 ψ (α,β;x,t)= √
n π2nn! f1/4 (4.12)
−ξ2/2+iT
×H (ξ)e .
TABLEII.SummaryofCPUtimewhenM isvariedforexample1. n
Thefixedparametersarer=6,J =150,andh=0.2667.
ThevariousquantitiesinEq.(4.12)aredefinedasfollows:
M N ∆t CPU(s) e (10−7)
2 (cid:112)
α= mω/(cid:126), f =α4cos2ωt+β4sin2ωt
1 2500 0.0004 398 10.00
2 130 0.0077 86.0 9.36 ξ =β[α2(x−Acosωt)−ksinωt]/f1/2
3 30 0.0333 76.8 9.74 T =α2/(2f)(cid:8)[(β4−α4)x2−k2+β4A2]sinωtcosωt
4 25 0.0400 78.4 9.35 +2[α4kxcosωt+β4A(ksinωt−α2x)sinωt](cid:9)
5 20 0.0500 92.6 9.27
(cid:18)β2sinωt(cid:19) (cid:18)ωt+π(cid:19)
6 18 0.0556 108 9.28
θ =arctan +2πν, ν =int .
7 15 0.0667 128 9.28 α2cosωt 2π
(4.13)
8 14 0.0714 142 9.28
9 12 0.0833 153 9.27
Thewavefunction(4.12)isthepulsatingandoscillatingwave
10 11 0.0909 181 9.27 function of a particle subject to a one-dimensional harmonic
12 10 0.1000 215 9.27 oscillator characterized by ω or α. The initial (t = 0) wave
15 8 0.1250 266 9.27 function is the nth energy state of the particle subject to an
17 7 0.1429 303 9.27 oscillator characterized by β, rather than α, displaced from
20 6 0.1667 353 9.27 theoriginbyamountAandwithamomentum(cid:126)k. Thefunc-
tion H (ξ) is the nth-order Hermite polynomial. This wave
n
function provides a wave packet with more fluctuation than
C. Example2:oscillatingandpulsationharmonicoscillator thetraditionalcoherentwavepacket;forinstance,ithasnodes
wavefunctions whichtravelwiththepacketandwhoseoccurrencespreadand
contractintime.
Consider the potential function for the two-dimensional As an initial study of the accuracy of the method we con-
anisotropicharmonicoscillator, siderthesimplestcaseofanisotropicoscillatorwiththeinitial
state the ground state. The values of the parameters are give
V(x,y)= 1m(ω2x2+ω2y2). (4.10) inTableIII.InFig.4theerrorasafunctionofM, theorder
2 x y ofthediagonalPade´ approximant,isdisplayed. Notethatthe
lowervaluesofM,especiallyforlargerr givenoresultsbe-
Ananalyticsolutionforsuchapotentialis[28]
cause the convergence of the iterative part of the calculation
Ψ(x,y,t)=ψ (α ,β ;x,t)ψ (α ,β ;y,t), (4.11) is not achieved. By decreasing dt convergence can again be
nx x x ny y y
7
TABLEIII.ParametersusedfortheerrorcalculationsofFigs.4and
5.
100 M = 30
(cid:126)=m=1, nx =ny =0 kx =ky =0 M = 15
x0 =y0 =−15, xJ =yM =15 αx =βx =αy =βy =1 10-2 MM == 54
t =2π, N =100, dt=2π/N A =A =2 M = 3
max x y 10-4 M = 2
J=M=100 M = 1
formula
e2 10-6
attained, but in Fig. 4 we keep dt constant . The horizontal 10-8
plateauxarenotcompletedsincethereisnochangeintheer-
rorasM isfurtherincreased. Thecalculationsaredonewith 10-10
double-precisionfloating-pointarithmetic. Weachieveaner-
rorlessthan10−11 forr = 30andM ≥ 6. Theerrorscould 10-12
0 5 10 15 20 25 30
befurtherreducedbyincreasedcomputationalprecision.
r
100 r = 1 FIG.5. (Colouronline)Theerrorsofthenumericalwavefunctions
r = 5 asafunctionofrfordifferentvaluesofM.Theparametersusedare
10-2 r = 10 listedinTableIII.
r = 15
r = 20
10-4 r = 25
r = 30
whereξ issomevalueinthedomainofx. Wehaveassumed
e2 10-6
that x dependence of the wave function of this example is a
10-8 Gaussianandthatthemaximumvalueofitanditseven-order
derivativesoccurwhentheargumentiszero[29,p. 933],i.e.,
10-10
(−1)r(2r)!
10-12 G(x)=e−x2/2, G(2r)(0)= 2r(r!) . (4.15)
0 2 4 6 8 10 12 14
M
ThelastfactorontherightsideofEq.(4.14)isanadjustment
togivereasonableagreementwiththedata. Itamountstoan
effectivespatialstepsizewhichissmallerbyafactorof2.5.
FIG.4. (Colouronline)Theerrorsofthenumericalwavefunctions The shape of the solid curve is very sensitive to the form of
astheyrelatetotheorderofthePade´approximantforvariousorders thisfactor.
ofthespatialexpansion.TheparametersusedarelistedinTableIII.
We plot the progression of the oscillating and pulsating
wavepacketasnumericallydeterminedinFig.6. Theparam-
Forthesamemodelweplottheerrorasafunctionofr in
etersusedforthiscalculationarelistedinTableIV.Theerror
Fig. 5 for a number of values of M. The lower pattern of
e inthecalculationrangesfrom3×10−7afteronetimestep
dotsisonethatisobtainedforeachM valueuptoaparticu- 2
larvalueofr,atwhichtheerrorisconstantasr isincreased
further. These horizontal plateaux only extend to a certain
TABLEIV.Parametersusedforcalculationoftheoscillating,pulsat-
point after which the iterative procedure becomes unstable.
ingwavefunctionshowninFig.6.
TheplateauxareclearlyvisibleforM =4and5,andthebe-
ginningscanbediscernedforM = 1,2,3. Theinstabilityof M =4, r=14, nx =2, ny =1 (cid:126)=m√=1, kx =0, ky =5
thecalculationdoesnotmeanthatwecannotobtainresultsin x0 =y0 =−15, xJ =yM =15 αx = 2, αy =1
those regions. In this calculation dt is the same in all cases; t =2π, N =557, dt=2π/N β =2α , β =2α
max x x y y
whereinstabilitysetsinasmallerdtwillrestorestability. The J=M=200 A =−5, A =0
x y
factthatthecurvessuperimposeontheleftcanbeseenfrom
Fig.4whereforeachM ≥ 6theerrorsconvergeforr suffi-
to5×10−3after557steps.Notethattheoscillatingfrequency
cientlylarge.
isω = 2ω . Thismeansthatwhenthemotionhasexecuted
OnFig.5wehaveplottedasolidlinewhichisanestimate x y
a complete cycle in the x direction it has only gone through
oftheerrorobtainedbyconsideringthetruncationerrorofthe
halfacycleintheydirection. Thusthepacketstartsat(-5,0),
expansioninxory,i.e.,
travelsalongaquarter-ellipticalpathto(5,5), thenbackto(-
e2(r) =C2h2r ≈myax(cid:12)(cid:12)(cid:12)(cid:12)(21r)!∂∂x2r2rΨ(x,y,T)|x=ξ(cid:12)(cid:12)(cid:12)(cid:12)h2r×2.512r, 5fr,e0q),uetonc(i5e,s-5in),eaanchddtoiretchteioinniatiraelfpoouirnttim(-e5s,0th).ecTohrerepspuolsnadtiinngg
(4.14) oscillatingfrequencies.
8
0.4 0.4
|Ψ| 0.2 0 0.1 0.2 |Ψ| 0.2
0 0
10 10
5 5
-10 -5 x 0 5 10-10-5 0 y -10 -5 x 0 5 10-10-5 0 y
0.4 0.4
|Ψ| 0.2 |Ψ| 0.2
0 0
10 10
5 5
-10 -5 x 0 5 10-10-5 0 y -10 -5 x 0 5 10-10-5 0 y
FIG.6.(Colouronline)Oscillatingandpulsatingwavepacketintwodimensionsofexample2.Theplotintheupperleftpanelrepresentsthe
wavefunctionattimet = 0,intheupperrightpanelatt = π/2,inthelowerleftpanelatt = π,andinthelowerrightpanelatt = 3π/2.
Seeanimationhere.
D. Example3:freewavepacket with
ForthefreewavepacketweconsidertheHermite-Gaussian α[(z−A)−kτ]
ξ = √ , θ =arctan(α2τ),
wavefunctionofRef.[28],
1+α4τ2
(4.18)
(cid:18) α (cid:19)1/2
N(α)= √ and τ =(cid:126)t/m.
Ψ(x,y,t)=ψ(α ,n ;x,t)ψ(α ,n ;y,t), (4.16) π2nn!
x x y y
where The travelling wave packet will have nodes whose distribu-
tion,ifthereismorethanonenode,spreadintime.Themodel
issimilartothatofGalbraithetal.[30]inwhosecalculation
ψ(α,n;z,t)=√Nn(α) exp(cid:18)i(z−A)2(cid:19)e−inθ nx =ny =0andαx =αy. Theparametersweusearegiven
1+iα2τ 2τ inTableV.
(cid:18) ξ2 ξ2 (cid:19) Thefreewavepacketattimest=10−3andt=6.3×10−3
×H (ξ)exp − −i
n 2 2α2τ is shown in Fig. 7. The separation of the peaks of the wave
(4.17) functionastimeprogressesisclearlyevident.
9
25 25
20 20
0 5 10 0 5 10
|Ψ| 15 |Ψ| 15
10 10
5 5
0 0
-2 2 -2 2
-1 1 -1 1
x 0 1 -1 0 y x 0 1 -1 0 y
2 -2 2 -2
FIG. 7. (Colour online) The movement and dispersion of a free wave packet in example 3 shown at t = 10−3 in the left panel and at
t=6.3×10−3intherightpanel.Seeanimationhere.
TABLEV.Parametersusedforthecalculationsofthefreepacketin TABLEVI.Parametersusedforthecalculationsofthewavepacket
Fig.7. passingthroughasingleslitasshowninFigs.8and9.
(cid:126)=2m=1, n =2, n =1 k =k =32 (cid:126)=2m=1, n =n =0 k =32,k =0
x y x y x y x y
x0 =y0 =−2.5, xJ =yM =2.5 αx =αy x0 =y0 =−2, xJ =yM =2 αx =αy
N =80, dt=10−4 A =A =0.25 N =79, dt=10−4 A =−0.25,A =0
x y √ x √ y
J=M=200, r=20, M =10 α =12.5 2 J=M=200, r=14, M =6 α =12.5 2
x x
E. Example4:single-slitdiffraction can study different slit configurations and shapes [31] using
thismethod.
Thewavenatureofelectronshasbeenstudiedandobserved
insemiconductornanostructures. Endohetal.[31]havecon-
siderednumericalsimulationsofthepassageofsuchelectrons
through narrow constrictions. Recent experiments observed 3.5×10-10
controlledelectrondiffractionforbothsingle-anddouble-slit 3.0×10-10
configurations[32,33].
The single slit in the barrier is obtained by introducing a 2.5×10-10
potential 2y,t)| 2.0×10-10
V(x,y)=V0f(x)[f(0)−f(y)], (4.19) Ψ(x,0 1.5×10-10
|
wheref(x)isthedifferenceoftwoFermifunctions 1.0×10-10
1 1 5.0×10-11
f(x)= − . (4.20)
−µ(x+x ) −µ(x−x )
1+e w 1+e w
-1.5 -1 -0.5 0 0.5 1 1.5
In the calculation we choose µ = 100, x = 0.05, and
w y
V0 = 1000. Initially the free wave f√unction (4.16) (with
n = n = 0 and α = α = 12.5 2) impinges on the
x y x y
slitanddiffracts. Theparametersofthesingle-slitcalculation
aregiveninTableVI. FIG.8. (Colouronline)Theprobabilitydensityofthewavepacket
Weplottheprobabilitydensity|Ψ(x0,y,t0)|2asafunction diffractedbyasingleslitontheplanewherex=x0 =1.25333ata
timet=0.0057.
of y when x = 1.253333 and t = 0.0057 in Fig. 8. The
0 0
graphhasaremarkablesimilaritytotheFraunhoferdiffraction
intensity.Theslitisnotofuniformwidthandhencewecannot InFig.9weplottwosnapshotsofthewavepacketpassing
compare parameters. The calculation does indicate that one throughtheslit. Sincemostofthepacketisreflected,wemul-
10
5 5
4 4
|Ψ| 3 1 2 3 |Ψ| 3 1 2 3
2 2
1 1
0 0
2 2
1 1
0 y 0 y
-2 -1 x 0 1 2 -2-1 -2 -1 x 0 1 2 -2-1
FIG.9.Thediffractionofawavepacketpassingthroughasingleslit,indicatedbyframedoutlinesonthegraphs.Thetimeintheleftpanelis
t=0.0040andinrightpanelist=0.0079.Seeanimationhere.
tiplytheamplitudeofthediffractedpacketinthefigurebyten O(h2r+∆t2M+1)whereM andrarepositiveintegers.
in order to make the packet’s shape in the region beyond the The examples demonstrate that this method is capable of
slitmorevisible. accuratesolutionsevenwhenthereissignificantfluctuationof
thewavefunction. Themethodsdescribedinthispaperallow
onetoobtainarealistictheoreticalanalysisofthediffraction
V. DISCUSSION experiments that have been done recently. Given the recent
attentiontothetwo-dimensionalnonlinearSchro¨dingerequa-
We have shown that the generalized Crank-Nicolson tionandthePitaevskiiequation,afutureprojectistoexpand
methodcombinedwiththealternating-directionimplicitpro- themethodofthispapertosuchsystems,aswellasthosewith
cedureisapracticalapproachtothedeterminationofnumer- time-dependentinteractionsorsourceterms. Thisineffectis
ical solutions of the two-dimensional Schro¨dinger equation. ageneralizationoftheworkdoneearlieronone-dimensional
The method allows one to study the efficiency and accuracy systems[28].
systematically as a functions of powers of the temporal and Furthermoreitremainstosystematicallyinvestigatetherel-
spatial step sizes. Equation (2.15) is basic to the solution. It ativeefficiencyandaccuracyoftheapproachofthispaperto
canbesolvedindifferentways,butwechoosetouseaniter- othermethodsthathavebeenusedorproposed. Ageneraliza-
ativeapproachwhichmeansonemustfindsolutionsoflinear tiontothreeorhigherspatialdimensionsandtheintroduction
systems of equations whose coefficients form banded diago- of transparent boundary conditions are further natural exten-
nalmatrices. Sincethenumberofiterationsislow,thealter- sionsofthiswork.
nativenoniterativeapproachleadstolesssparsematricesand
acorrespondinglylessefficientprocedure.
Anumberofauthors[15,16,34]haveconsideredalternat- ACKNOWLEDGMENTS
ingdirectionimplicitcompactfinitedifferenceschemeswhich
giveerrorsoforderO(h5+∆t3).Althoughtheyincludenon- We are grateful to the Natural Sciences and Engineering
linearequationsintheiranalysis,theydiscuss,amongothers, Research Council of Canada for the Undergraduate Student
example1ofthispaperasatestcase.Ourschemegiveserrors ResearchAwardsgiventoS.-J.P.(2013)andT.V.(2016).
[1] H.Tal-EzerandR.Kosloff,J.Chem.Phys.81,3967(1984). [6] W. van Dijk and F. M. Toyama, Phys. Rev. E 90, 063309
[2] C.Leforestier,R.H.Bisseling,C.Cerjan,M.D.Feit,R.Fries- (2014).
ner,A.Guldberg,A.Hammerich,G.Jolicard,W.Karrlein,H.- [7] W.vanDijk,Phys.Rev.E93,063307 (2016).
D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, J. Comp. [8] M. Forma´nek, M. Va´nˇa, and K. Houfek, in Numerical Anal-
Phys.94,59(1991). ysisandAppliedMathematics,InternationalConference2010,
[3] I.Puzynin,A.Selin, andS.Vinitsky,Comp.Phys.Comm.123, editedbyT.E.Simos,G.Psihoyios, andC.Tsitouras(Ameri-
1(1999);126,158(2000). canInstituteofPhysics,2010)pp.667–670.
[4] W. van Dijk and F. M. Toyama, Phys. Rev. E 75, 036707 [9] A.A.Gusev,O.Chuluunbaatar,S.I.Vinitsky, andA.G.Abra-
(2007). shevich,Math.Mod.andGeom.2,33(2014).
[5] W. van Dijk, J. Brown, and K. Spyksma, Phys. Rev. E 84, [10] A. Chuluunbaatar, V. L. Derbov, A. Galtbayar, A. A. Gusev,
056703(2011). M.S.Kashiev,S.I.Vinitsky, andT.Zhanlav,J.Phys.A:Math.