{VERSION 4 0 "IBM INTEL NT" "4.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 
1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 20 "Ehrenfest's theorems" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 369 "We first
 develop a procedure which propagates a one-dimensional wavepacket in \+
space and time. Then we demonstrate the validity of the theorems, whic
h determine the time evolution of position and momentum expectation va
lues. This allows one to investigate the semi-classical limit in which
 the position and momentum average values satisfy Hamilton's equations
 of motion." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 83 "The first step involves setting up a coordinate- and mo
mentum-space discretization." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "dx:=0.3; n
:=10; N:=2^n;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dp:=evalhf
(2*Pi/(N*dx));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "readlib(F
FT):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "with(LinearAlgebra)
: Digits:=15:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "xv:=Vector
(N): yR:=Vector(N): yI:=Vector(N): pv:=Vector(N):" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 129 "for i from 1 to N do:\nxv[i]:=evalhf((i-N/2
-1/2)*dx); if i < N/2+1 then\npv[i]:=evalhf((i-1)*dp); pv[N-i+1]:=eval
hf(-i*dp); fi; od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "p_max
:=(N/2-1)*dp;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 172 "The maximum ava
ilable momentum value on the mesh puts a limit on the momentum distrib
ution that we can choose for the intial wavepacket! It can be increase
d by reducing dx." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 233 "Now we set \+
up the initial wavepacket: we can choose the initial average position,
 average momentum, and momentum spread. We should also have a choice f
or the shape of the wavepacket. To make matters simple we choose a Gau
ssian shape:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "x0:=-50; p0
:=1; sigma:=1/4;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "fx:=exp
(-sigma^2*(x-x0)^2/2)*exp(I*p0*x)/(Pi/(sigma)^2)^(1/4);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "int(evalc(conjugate(fx)*fx),x=-infi
nity..infinity);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 128 "The initial \+
wavepacket is normalized properly! The width parameter sigma describes
 the momentum spread around the average value." }}{PARA 0 "" 0 "" 
{TEXT -1 38 "We can calculate the average momentum:" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 64 "int(evalc(conjugate(fx)*(-I)*diff(fx,x)),
x=-infinity..infinity);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 21 "The av
erage position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "int(eval
c(conjugate(fx)*x*fx),x=-infinity..infinity);" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 62 "We can also calculate the deviation from the average p
osition:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "sqrt(int(evalc(
conjugate(fx)*(x-x0)^2*fx),x=-infinity..infinity));" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 49 "plot([Re(fx),Im(fx)],x=-75..15,color=[red
,blue]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 257 11 "Exercise 1:" }}
{PARA 0 "" 0 "" {TEXT -1 184 "Vary the parameters that determine the G
aussian wavepacket (particularly the momentum spread), and observe the
 change in the packets appearance (in particular, its spread in positi
on)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "for i from 1 to N d
o: yR[i]:=evalf(subs(x=xv[i],Re(fx))); yI[i]:=evalf(subs(x=xv[i],Im(fx
))); od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 197 "Now we pick a potent
ial function. In principle, the potential should not vary as a functio
n of x in the region where the wavepacket is defined initially, i.e., \+
that region should be potential-free." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 3 "x0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "Vp:=
x->2*(exp((x-100)/50)-exp((x0-100)/50));" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 65 "Vi:=Vector(N):\nfor i from 1 to N do: Vi[i]:=evalf(Vp
(xv[i])); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots
):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 502 "We graph the real and imag
inary parts of the wavepacket's discretized wavefunction displaced by \+
the amount of average initial kinetic energy along with the potential.
 The baseline for the wavefunction plot indicates the classical turnin
g point (where it intersects the potential energy graph). Note that th
e variation in the potential energy is gentle at first (approximately \+
constant force over the extent of the wavepacket), while the variation
 becomes more steep as one approaches the turning point." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "P1:=plot([seq([xv[i],Vi[i]],i=1..N)
],color=green):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "P2:=plot
([seq([xv[i],p0^2/2+yR[i]],i=1..N)],color=blue):" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 55 "P3:=plot([seq([xv[i],p0^2/2+yI[i]],i=1..N)],c
olor=red):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "display([P1,P
2,P3],view=[-100..100,-0.5..1.5]);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 27 "dt:=0.1: dth:=0.5*dt:\nt:=0;" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 29 "T2:=Vector(N): V1:=Vector(N):" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 57 "for i from 1 to N do: T2[i]:=exp(-I*0.5*p
v[i]^2*dth); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "for i f
rom 1 to N do: V1[i]:=exp(-I*Vi[i]*dt); od:" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 460 "Tstep:=proc() local i,s; global y,xv,pv,Vi,Eni,n,
N,yR,yI,T2,V1,t,dt,dx;\nevalhf(FFT(n,var(yR),var(yI)));\nfor i from 1 \+
to N do: y:=T2[i]*(yR[i]+I*yI[i]); yR[i]:=Re(y); yI[i]:=Im(y); od:\nev
alhf(iFFT(n,var(yR),var(yI)));\nfor i from 1 to N do: y:=V1[i]*(yR[i]+
I*yI[i]); yR[i]:=Re(y); yI[i]:=Im(y); od:\nevalhf(FFT(n,var(yR),var(yI
)));\nfor i from 1 to N do: y:=T2[i]*(yR[i]+I*yI[i]); yR[i]:=Re(y); yI
[i]:=Im(y); od:\nevalhf(iFFT(n,var(yR),var(yI))); t:=(t+dt);\nt; end:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Fx:=unapply(-diff(Vp(x)
,x),x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 8 "At time " }{TEXT 258 2 
"t " }{TEXT -1 90 "= 0 we compare the force evaluated at the average p
osition to the force expectation value:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 18 "F0:=evalf(Fx(x0));" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 49 "Fexp:=dx*add((yR[i]^2+yI[i]^2)*Fx(xv[i]),i=1..N);" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 450 "We set up a time loop with an ex
terior loop in which we store graphs of the time evolution of the wave
function and the momentum-space probability at a large time interval. \+
An inner loop steps over dt: in this loop we collect information about
 the position and force operator expectation values. Our aim is to ana
lyze these expectation values to verify Ehrenfest's theorem, as well a
s to compare with classical time evolution for the average position." 
}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 101 "The ti
me loop requires a large amount of CPU time. Be prepared for hours to \+
complete the calculation!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
25 "Nt:=55: EKL:=[]: EPL:=[]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 31 "xL:=[[t,x0]]: FxL:=[[t,Fexp]]: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 23 "for jt from 1 to Nt do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "fo
r it from 1 to 50 do: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "t_i:=Tste
p(); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 146 "x0t:=dx*add((yR[i]^2+yI[i
]^2)*xv[i],i=1..N): xL:=[op(xL),[t_i,x0t]]: Fexp:=dx*add((yR[i]^2+yI[i
]^2)*Fx(xv[i]),i=1..N): FxL:=[op(FxL),[t_i,Fexp]]: " }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "Pre[jt]:=
plot([seq([xv[i],p0^2/2+yR[i]],i=1..N)],color=blue):" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 373 "yRt:=Vector(N): yIt:=Vector(N): for i from 1 to
 N do: yRt[i]:=yR[i]: yIt[i]:=yI[i]: od: evalhf(FFT(n,var(yRt),var(yIt
)));\nEP:=dx*add(Vi[i]*(yR[i]^2+yI[i]^2),i=1..N); EPL:=[op(EPL),[t_i,E
P]]: EK:=add(pv[i]^2/2*(yRt[i]^2+yIt[i]^2),i=1..N)/add(yRt[i]^2+yIt[i]
^2,i=1..N); EKL:=[op(EKL),[t_i,EK]]: Pmom[jt]:=plot([seq([pv[i],yRt[i]
^2+yIt[i]^2],i=1..N)],color=black,style=point):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "B:=a
nimate(Vp(x),x=-100..120,d=-1..1,frames=Nt,color=green):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "A:=display([seq(Pre[j],j=1..Nt)],in
sequence=true,view=[-100..120,-0.1..0.9]):" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 25 "display(A,B,thickness=3);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 352 "The wavepackets momentum distribution evolves as a fun
ction of time: the potential slows down the wavepacket, forces a turn-
around, and then accelerates it to the left. The momentum distribution
 changes as a result of the interaction: it becomes quite narrow as th
e wavepacket approaches the turning point, and then broadens again aft
er the turnaround." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "displ
ay([seq(Pmom[j],j=1..Nt)],insequence=true,view=[-2..2,0..500],symbol=d
iamond);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 117 "Let us check the tim
e evolution of the energy expectation values: kinetic, potential, and \+
total energy are displayed:" }}{PARA 0 "" 0 "" {TEXT -1 310 "According
 to the equations of motion for observables in QM the expectation valu
e of the total energy has to remain constant, as [H,H]=0 and H does no
t depend on timein our problem. The constancy of the energy expectatio
n value serves to demonstrate that the numerical propagation of the wa
vepacket is accurate." }}{PARA 0 "" 0 "" {TEXT -1 141 "We can also see
 how kinetic energy on average is converted into potential energy and \+
back as the wavepacket scatters from the potential wall." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "ETL:=[seq([EKL[i][1],EKL[i][2]+EPL[
i][2]],i=1..Nt)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "plot([
EKL,EPL,ETL],color=[red,blue,green],thickness=3);" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 53 "How small is the kinetic energy at the turning poi
nt?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "EKL[27],EKL[28],EKL[
29];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 167 "We did not resolve the t
ime axis well enough to answer the question whether the average kineti
c energy actually becomes zero! A detailed graph of the neighborhood o
f  " }{TEXT 260 2 "t " }{TEXT -1 167 "= 135-140 suggests that it does \+
not; this makes sense as the wavepacket cannot be a function with zero
 curvature. Nevertheless, it does come quite close to this limit." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 531 "
Now we compare the evolution of the position expectation value with pu
rely classical evolution. In classical evolution the force is evaluate
d only at the position x0(t) to control the further time evolution. Th
e quantum calculation of the average position via Ehrenfest's theorem \+
would rely on a sampling of the force over the extent of the wavepacke
t. At early times we expect the two methods to agree, at later times w
hen the wavepacket spreads and hits a region where the potential varie
s more rapidly we can expect differences." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 25 "X1:=plot([xL],color=red):" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 28 "NE:=diff(x(s),s$2)=Fx(x(s));" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 23 "IC:=x(0)=x0,D(x)(0)=p0;" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 55 "sol:=dsolve(\{NE,IC\},x(s),numeric,output=l
istprocedure);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "x_c:=subs
(sol,x(s));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "x_c(10);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "X2:=plot('x_c(s)',s=0..Nt*dt
*50,color=blue):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "display
(X1,X2,thickness=3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 176 "Note how
 the classical calculation for the position expectation value predicts
 the turnaround at a slightly later time, and a larger value of x0 at \+
which the turnaround occurs." }}{PARA 0 "" 0 "" {TEXT -1 170 "To verif
y that Ehrenfest's theorem holds for the numerical calculation we can \+
check whether the second derivative of x0(t) agrees with the expectati
on value of the force." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "
x2pL:=[]: for i from 2 to nops(xL)-1 do: x2pL:=[op(x2pL),[xL[i][1],(xL
[i+1][2]+xL[i-1][2]-2*xL[i][2])/dt^2]]; od:" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 183 "A:=plot([seq(FxL[20*i-10],i=1..nops(FxL)/20)],col
or=red,style=point,symbol=diamond): B:=plot([seq(x2pL[20*i-1],i=1..nop
s(x2pL)/20)],color=blue,style=point,symbol=cross): display(A,B);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 355 "Quite clearly the numerical solut
ion of the TDSE satisfies Ehrenfest's theorem: otherwise the evaluatio
n of the second derivative of the average position with respect to tim
e would not coincide with the expectation value of the force. We can u
se this as a confirmation of the validity of the numerically obtained \+
result for the position expectation value." }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 282 "On the other hand, we find tha
t the position expectation value x0(t) follows the classical evolution
 only to the point where the wavepacket experiences a substantially di
fferent amount of acceleration in the front versus back parts. This oc
curs shortly before the turnaround point. " }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 295 "An interesting question arises
: suppose we have a classical statistical ensemble. This ensemble cont
ains faster and slower particles. What is the evolution of the particl
e density? The faster particles will have more energy, i.e., a turning
 point further to the right, than the slower particles." }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 259 41 "Classical s
tatistical ensemble simulation" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 145 "We dial up an initial distribution that \+
corresponds to a wavepacket: intially there is no correlation between \+
the particle position and momentum." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 9 "NEns:=50:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
12 "x0,p0,sigma;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 163 "We need a Ga
ussian distribution in position (width 1/sigma) around x0, and a distr
ibution in momentum around p0. This can be achieved with the normal di
stribution." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "Pval:=Vector
(NEns): Xval:=Vector(NEns):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
12 "with(stats):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "xrand:=
stats[random, normald](NEns):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 56 "for i from 1 to NEns do: Pval[i]:=p0+sigma*xrand[i]: od:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "xrand:=stats[random, normald
](NEns):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "for i from 1 to
 NEns do: Xval[i]:=x0+(1/sigma)*xrand[i]: od:" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 56 "The expectation value follows from the ensemble averag
e:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "xC0:=add(Xval[i],i=1.
.NEns)/NEns;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "pC0:=add(Pv
al[i],i=1..NEns)/NEns;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "E
PC0:=evalf(add(Vp(Xval[i]),i=1..NEns)/NEns);" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 40 "EKC0:=add(0.5*Pval[i]^2,i=1..NEns)/NEns;" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 202 "For each testparticle (ensemble m
ember) we would like to propagate the Newton equation. We can calculat
e expectation values as in quantum mechanics, except that they are cal
culated as ensemble averages." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 58 "We have a time loop as before, and a part
icle loop inside." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "t_i:=0:
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "EKCL:=[[t_i,EKC0]]: EPCL:=[[t_i
,EPC0]]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "xCL:=[[t_i,xC0]]: for j
t from 1 to Nt do: t_f:=evalf(t_i+50*dt): for j from 1 to NEns do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "IC:=x(t_i)=Xval[j],D(x)(t_i)=Pval[j
];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "sol:=dsolve(\{NE,IC\},x(s),nu
meric,output=listprocedure);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "x_c
:=subs(sol,x(s)); p_c:=subs(sol,diff(x(s),s));" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 51 "Xval[j]:=x_c(t_f): Pval[j]:=p_c(t_f): od: t_i:=t_f:" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "xC:=add(Xval[i],i=1..NEns)/NEns;
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "EPC:=evalf(add(Vp(Xval[i]),i=1.
.NEns)/NEns); EKC:=add(0.5*Pval[i]^2,i=1..NEns)/NEns;" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 80 "EKCL:=[op(EKCL),[t_i,EKC]]: EPCL:=[op(EPCL),[t_i
,EPC]]: xCL:=[op(xCL),[t_i,xC]]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "
od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "ETCL:=[seq([EKCL[i][
1],EKCL[i][2]+EPCL[i][2]],i=1..Nt)]:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 58 "plot([EKCL,EPCL,ETCL],color=[red,blue,green],thicknes
s=3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 186 "The fact that the kinet
ic energy does not vanish at the turning point for the average positio
n means that the particles do not all reach their respective turning p
oints at the same time." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "
X3:=plot(xCL,color=green): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
30 "display(X1,X2,X3,thickness=3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 905 "The classical statistical mechanics simulation has a substanti
ally lower value for the location of the turning point: we can conclud
e that the wavepacket's tunneling behaviour allows it to penetrate the
 potential wall (which was very evident from watching the real part of
 the wavefunction). It is interesting that the classical calculation f
or the expectation value, which initially coincides with the QM calcul
ation for the wavepacket's true expectation value (justified by Ehrenf
est's theorem) leads to a slightly larger excursion into the classical
ly forbidden regime, and to a delay in the return of the wavepacket. T
his suggests that the tunneling process (as calculated by the numerica
l solution of the TDSE) may involve shorter time scales than classical
 propagation. The question of tunneling time, and its experimental ver
ification in particular, represents an active area of research current
ly." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 
0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }