Schrödinger's Equation
(2-dimensional version)

Terry Robb
<tdr@vaxc.cc.monash.edu.au>

This Notebook solves the two-dimensional time-dependent Schrödinger equation. Because some of the animations take up a large amount of space, I have only included the initial frame in all but the first experiment. Email to me <tdr@vaxc.cc.monash.edu.au> if you want the other animations, or if you want to see other experiments done.

Each of the numerical experiments took about 10 minutes to generate, with most of the time being taken up in rendering the animations.

See the Notebook on the one-dimensional time-dependent Schrödinger equation (also by me) for further comments relating to this application. This Notebook just contains the input and output lines, and no comments.

For yet more information about other symbolic–numeric applications involving the InterCall package and Mathematica just email a single line containing the words find InterCall to mathsource@wri.com.

Setup

<<InterCall.m

InterCall["!schroedinger"];

AddDefault[
   SCHROED2D[$PSI->$PSI,$DT,$T->$T,
             $DX->($B-$A)/($N-1),$A,$B,$N->{$N1,$N2},
             $WRK:>Null,$TMP:>Null,
             $VFUNC->$VFUNC,$LINEARBC->$LINEARBC,
             $N1->ROWS[$PSI],$N2->COLS[$PSI],
             $M->Max[$N1,$N2] ],
   S[Z[$N1,$N2],R,R,R[2],R[2],R[2],I[2],Z[$M],Z[$M],
         RF[R,R],S[Z,Z,Z,Z,Z,Z,R,R[2],I],I,I,I],
   "schroed.o"];

Import[{SCHROED2D}];

AddDefault[CLAMP->Function[{a0,a1,ac,b0,b1,bc,t,x,k},
             a0=1; a1=0; ac=0; b0=1; b1=0; bc=0;],
           S[Z,Z,Z,Z,Z,Z,R,R[2],I] ];
SetDefault[schroed2d,
           $A->{0,0}, $B->{1,1}, $LINEARBC->CLAMP]

dx={0.01,0.01}; dt=dx.dx/4;

plotit[] := Block[{psi=Peek[$PSI],t=Peek[$T]},
   ListDensityPlot[-Abs[psi]//Transpose, Mesh->False,
            PlotRange->{-7,0}, Epilog->drawpot,
            PlotLabel->PaddedForm[t/dt,{3,1}]];1; ];

Experiment 1: Scattering

Initial condition

k = {64.0, 64.0}; a = 0.08^2; x0 = {0.25, 0.25};
psi0=Table[Block[{x={x1,x2}-x0}, Exp[I*k.x-x.x/(2*a)]],
           {x1,0,1,dx[[1]]},{x2,0,1,dx[[2]]}];
psi0=psi0/Sqrt[Plus@@(Plus@@(psi0*Conjugate[psi0]))*dx[[1]]*dx[[2]]];

Potential function

u0=4*k.k; {xp1,xp2}={0.5,0.5}; rpsqr=0.08^2;
potential = Function[{x,t},
   If[(x[[1]]-xp1)^2+(x[[2]]-xp2)^2 > rpsqr, 0, u0]];

drawpot=Circle[{xp1,xp2}/dx,Sqrt[rpsqr]/dx];

Animation

schroed2d[0, $PSI->psi0,$T->0,$VFUNC->potential];

plotit[];

Do[Do[schroed2d[dt],{5}]; plotit[], {15}]

Experiment 2: Double slit

Initial condition

k = {64.0, 64.0}; a = 0.08^2; x0 = {0.25, 0.25};
psi0=Table[Block[{x={x1,x2}-x0}, Exp[I*k.x-x.x/(2*a)]],
           {x1,0,1,dx[[1]]},{x2,0,1,dx[[2]]}];
psi0=psi0/Sqrt[Plus@@(Plus@@(psi0*Conjugate[psi0]))*dx[[1]]*dx[[2]]];

Potential function

u0=4*k.k; dl = 0.025; dw=0.08; dh=0.15;
potential = Function[{x,t},
   If[ x[[1]]+x[[2]]> 1+dl ||
       x[[1]]+x[[2]]< 1-dl ||
       (x[[1]]-x[[2]]< dh && x[[1]]-x[[2]]> dw) ||
       (x[[1]]-x[[2]]<-dw && x[[1]]-x[[2]]>-dh),
      0, u0] ];

px1= 0.00;         py1= 1-dl;
px2=(1-dl+ -dh)/2; py2=(1-dl- -dh)/2;
px3=(1+dl+ -dh)/2; py3=(1+dl- -dh)/2;
px4= 0+dl;         py4= 1.00;
qx1= 1-dl;         qy1= 0.00;
qx2=(1-dl+dh)/2;   qy2=(1-dl-dh)/2;
qx3=(1+dl+dh)/2;   qy3=(1+dl-dh)/2;
qx4= 1.00;         qy4= 0+dl;
rx1=(1-dl+ -dw)/2; ry1=(1-dl- -dw)/2;
rx2=(1+dl+ -dw)/2; ry2=(1+dl- -dw)/2;
rx3=(1+dl+dw)/2;   ry3=(1+dl-dw)/2;
rx4=(1-dl+dw)/2;   ry4=(1-dl-dw)/2;
drawpot={Line[{{px1,py1}/dx,{px2,py2}/dx,
               {px3,py3}/dx,{px4,py4}/dx}],
         Line[{{qx1,qy1}/dx,{qx2,qy2}/dx,
               {qx3,qy3}/dx,{qx4,qy4}/dx}],
         Line[{{rx1,ry1}/dx,{rx2,ry2}/dx,
               {rx3,ry3}/dx,{rx4,ry4}/dx,
               {rx1,ry1}/dx}]};

Animation

schroed2d[0, $PSI->psi0,$T->0,$VFUNC->potential];

plotit[];

Do[Do[schroed2d[dt],{5}]; plotit[], {15}]

Experiment 3: Interaction

Initial condition

ka = { 64.0,64.0}; aa = 0.08^2; xa0 = {0.25, 0.25};
kb = {-64.0,64.0}; ab = 0.08^2; xb0 = {0.75, 0.25};
psi0=Table[Block[{xa={x1,x2}-xa0,xb={x1,x2}-xb0},
                 Exp[I*ka.xa-xa.xa/(2*aa)]+
                 Exp[I*kb.xb-xb.xb/(2*ab)]],
           {x1,0,1,dx[[1]]},{x2,0,1,dx[[2]]}];
psi0=1.4*psi0/Sqrt[Plus@@(Plus@@(psi0*Conjugate[psi0]))*dx[[1]]*dx[[2]]];

Potential function

potential = Function[{x,t}, 0];

drawpot={};

Animation

schroed2d[0, $PSI->psi0,$T->0,$VFUNC->potential];

plotit[];

Do[Do[schroed2d[dt],{5}]; plotit[], {17}]

Experiment 4: Stirring

Initial condition

a = 0.08^2; x0 = {0.5, 0.5};
psi0=Table[Block[{x={x1,x2}-x0}, Exp[-x.x/(2*a)]],
           {x1,0,1,dx[[1]]},{x2,0,1,dx[[2]]}];
psi0=psi0/Sqrt[Plus@@(Plus@@(psi0*Conjugate[psi0]))*dx[[1]]*dx[[2]]];

Potential function

u0=20*k.k; {xp1,xp2}={0.5,0.5}; rpsqr=0.04^2;
omega=N[2Pi/(75*dt)]; amp=0.16;
potential = Function[{x,t},
   If[(x[[1]]-xp1-amp*Cos[t*omega])^2+
      (x[[2]]-xp2)^2 > rpsqr, 0, u0]];

drawpot:={RGBColor[1,0,0],
          Circle[{xp1+amp*Cos[t*omega],xp2}/dx,
                 Sqrt[rpsqr]/dx]};

Animation

schroed2d[0, $PSI->psi0,$T->0,$VFUNC->potential];

plotit[];

Do[Do[schroed2d[dt],{5}]; plotit[], {15}]

Experiment 5: Tunnelling

Initial condition

k = {64.0, 64.0}; a = 0.08^2; x0 = {0.25, 0.25};
psi0=Table[Block[{x={x1,x2}-x0}, Exp[I*k.x-x.x/(2*a)]],
           {x1,0,1,dx[[1]]},{x2,0,1,dx[[2]]}];
psi0=psi0/Sqrt[Plus@@(Plus@@(psi0*Conjugate[psi0]))*dx[[1]]*dx[[2]]];

Potential function

u0=2*k.k; dl = 0.025;
potential = Function[{x,t},
   If[x[[1]]+x[[2]]>1+dl || x[[1]]+x[[2]]<1-dl, 0, u0] ];

drawpot={Line[{{0.,1-dl}/dx,{1-dl,0.}/dx}],
         Line[{{dl,1.00}/dx,{1.00,dl}/dx}]};

Animation

schroed2d[0, $PSI->psi0,$T->0,$VFUNC->potential];

plotit[];

Do[Do[schroed2d[dt],{5}]; plotit[], {15}]

Comments

Numerics

The SCHROED2D routine is based around a 2D Crank-Nicholson solver for parabolic PDE's. Specifically it solves the equation (* see Footnote)

i dpsi/dt = -(d^2psi/dx^2 + d^2psi/dy^2) + V*psi

If you are interested in the numerical method, then the 'Numerical Recipes' book is a reasonable place to get a basic description and some references.

*Footnote:
[Units are chosen so that speed of light c=1, Plank's constant hbar=1, and mass==1/2. An older version of the SCHROED2D routine had mass==1, and hence solved
   i dpsi/dt = -1/2*(d^2psi/dx^2 + d^2psi/dy^2) + V*psi
but the physics is identical if you take account of a change in potential V and time t.]

Timings

The implementation is unconditionally stable (ie DT can be large and the solution won't blow up), and is second order accurate in both space and time. Thus halving DX and DT will lead to a four-fold improvement in accuracy. However the computational time required is proportional to the number of grid points, and since the problem is 2D then halving DX will lead to a four-fold increase in computation, and halving DT will double that again.

Arguments

The arguments to SCHROED2D are as follows:

PSI - the complex wave-function. [complex array]
        On input PSI should contain the wave-function for time T.
        On output PSI will contain the wave-function for time T+DT.
  DT -  the time-step to evolve the wave-function by. [real number]
  T -   the time. [real number]
        On input T should contain the current time.
        On output T will contain the updated time T+DT.
  DX -  the spatial resolution of the region. [real 2-vector].
  A -   the lower left hand corner of the region. [real 2-vector].
  B -   the upper right hand corner of the region. [real 2-vector].
  N -   the number of grid points: Note DX==(B-A)/(N-1). [integer 2-vector].
  WRK - work space, for internal use only. [complex vector]
  TMP - work space, for internal use only. [complex vector]
  VFUNC -    the potential function V(x,t) with x a 2-vector. [real function]
  LINEARBC - the boundary conditions (see below). [subroutine]

Boundary conditions

The boundary condition subroutine LINEARBC required by SCHROED2D needs some further explanation. Boundary conditions that SCHROED2D can handle must be linear, ie a0*psi + a1*psi' = ac  at the lower boundary, and b0*psi + b1*psi' = bc  at the upper boundary.

These are general enough to include simple clamping conditions of psi=0 on the boundaries, which essentially keeps the wave-function in a rectangular box. Also the condition psi'=0 could be used to let the wave-function pass through the boundaries. Anyway one must supply a subroutine that can compute a0, a1, ac, and b0, b1, bc at any given time t. Time-dependent boundary conditions are thus also possible.

The arguments to the boundary condition subroutine must be {a0,a1,ac,b0,b1,bc,t,x,k} with {a0,a1,ac,b0,b1,bc} being output from the subroutine, and {t,x,k} being input to the subroutine.

If k==1, then {a0,a1,ac} should be set to values appropriate for the left edge of the region, and {b0,b1,bc} for the right edge. The second component of x will give the position along the left or right edge that the values are needed.

If k==2, then {a0,a1,ac} should be set to values appropriate for the bottom edge of the region, and {b0,b1,bc} for the top edge. The first component of x will give the position along the bottom or top edge that the values are needed.

examples

For example, the Mathematica function
     Function[{a0,a1,ac,b0,b1,bc,t,x,k}, a0=1; a1=0; ac=0; b0=1; b1=0; bc=0;]
can be used to give simple clamping conditions, ie psi=0, on all four boundaries.

Another example, which lets the wave-function pass through all the boundaries would be:
     Function[{a0,a1,ac,b0,b1,bc,t,x,k}, a0=0; a1=1; ac=0; b0=0; b1=1; bc=0;]
This sets psi'=0 on each boundary. (Here the "psi'" means the derivative of psi with respect to the normal component to the wall.)

A more complicated example which forces psi to be Sin[t]*Sin[Pi*y] on the right wall (x=1) and zero on the other three walls would be:
     Function[{a0,a1,ac,b0,b1,bc,t,x,k}, (* x[[1]] is x and x[[2]] is y *)
        If[k==1, (* x walls *)
              (* left wall:    *) a0=1; a1=0; ac=0; 
              (* right wall:  *) b0=1; b1=0; bc=Sin[t]*Sin[Pi*x[[2]]];],
        If[k==2, (* y walls *)
              (* lower wall: *) a0=1; a1=0; ac=0;
              (* upper wall: *) b0=1; b1=0; bc=0;] ]

Periodic boundary conditions are specified by setting a0=0 and a1=0, but this is not implemented yet.

Experiments

Values for DX and DT were chosen so as to give high resolution and accuracy.

The following point may be obvious but is often forgotten: For a test run just use a small grid, say 20 X 20, as that will be fast. If you are happy with the results then go back and do a high resolution run, say 100 x 100.

Schroedinger's equation is diffusive, so after a while any wave-packet will eventually diffuse away. Thus wave-packets of initial form Exp[I*k*x-x^2/(2a)] (with "a" the initial spread) will naturally spread out over time. After time t the spread should become Sqrt[a^2+(2*t)^2].

In some of the experiments the speed and size of the initial wave-packet has been chosen so that in the time it takes to cross the region of interest it won't diffuse too much. Other cases are worth looking at too, and one can confirm the formula Sqrt[a^2+(2*t)^2].

Settings

About Peek

Peek is a low-level InterCall command that is useful for debugging and some other situations. If an object (such as a large array, or a function etc) is passed as an argument to an external routine, then a special copy of that object is created outside Mathematica. To pull that copy back to Mathematica you can use Peek. The InterCall manual documents three commands called: Peek, Poke and Purge.

Peek and SCHROED2D

Normally one never needs to use Peek, so why is it used in this application?

Well the reason is to do with efficiency. The SCHROED2D routine was written so that it would evolve an array (the wave-function PSI) forward a single time-step DT. A sensible InterCall default setting would thus seem to be
   AddDefault[ SCHROED2D[$PSI->In:>Out, $DT->In, ...
which says that PSI will be supplied (In) by the user, and then its updated value will be returned (Out). This would work fine, but it means that PSI (potentially a very large array) gets exchanged between Mathematica  and schroed2d *every* time-step. And because we need to call schroed2d many times with a small time-step before producing a single frame, the above AddDefault setting will result in a lot of unnecessary communication.

The solution is to instead use the setting,
   AddDefault[ SCHROED2D[$PSI->$PSI, $DT->In, ...
which tell's InterCall that SCHROED2D will already have a copy of the latest value of PSI so it won't need to be copied again.

OK, but how do we get the first copy of PSI (ie. the initial condition) off to schroed2d? And how do we get the latest copy back to Mathematica so we can plot it?

The answer to the second question is to use the Peek command. Peek[$PSI] will pull back the current value of PSI to Mathematica's end. Apart from simple debugging this is probably the only useful application of Peek.

The answer to the first question is to use an option at the end of the very first call to schroed2d. Thus to pass the initial condition psi0 to schroed2d, we could use
   schroed2d[0, $PSI->psi0]
where the option $PSI->psi0 will override the default setting of $PSI->$PSI that has been set with the AddDefault command. The first argument in schroed2d is 0 because we don't want it to evolve PSI just yet, so we set DT to 0 for the first call. Subsequent calls are then, say
   schroed2d[0.0001]
which will evolve PSI by an amount DT=0.0001 of time. Smaller time-steps are better because they lead to greater accuracy in the numerical scheme used; but of course it takes longer to get to a given time T.

InterCall

InterCall was developed by Terry Robb. For more information on InterCall send email to intercall_forward@wri.com. An order form is available via MathSource – email the line send 0202-587-0011 to mathsource@wri.com for details. Information can also be posted or faxed to Analytica below.

Analytica                     Fax: +61 (9) 386 5666
P.O. Box 522
Nedlands 6009
Perth WA
AUSTRALIA.