{VERSION 5 0 "IBM INTEL NT" "5.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 "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 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 0 0 1 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 }{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 24 "Discrete Wigner Function " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 213 "We u se the discrete Fourier transform to carry out Wigner transforms of wa vepackets that represent more complicated situations than those avail able by means of symbolic integration, and explored in the worksheet \+ " }{TEXT 19 10 "Wigner.mws" }{TEXT -1 235 ". We are interested, in par ticular, in the Wigner distribution for a tunneling wavepacket as calc ulated by numerical techniques. First, we test the calculation on a si mple Gaussian wavepacket for which the analytic answer is available." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "restart; a:=1; p0:='p0';" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 102 "psi := unapply(sqrt(a)*1/Pi^(1/4)*exp(-1/2*(x-I*a^2* p0)^2/(a^2+t*I)-1/2*a^2*p0^2)/(sqrt(-a^2-I*t)),x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "assume(t>0,p,real);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 58 "int(evalc(simplify(abs(psi(x))^2)),x=-infinity ..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "rho:=simpli fy(expand(psi(x+s/2)*conjugate(psi(x-s/2))));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 59 "Our task is now to perform a discrete FT over the vari able " }{TEXT 19 1 "s" }{TEXT -1 27 " to map it into a momentum " } {TEXT 19 1 "p" }{TEXT -1 44 ". This is to be carried out for a number \+ of " }{TEXT 19 1 "x" }{TEXT -1 53 " values, and in this way a distribu tion over momenta " }{TEXT 19 1 "p" }{TEXT -1 15 " and positions " } {TEXT 19 1 "x" }{TEXT -1 101 " will be obtained. Of course, we have to pick numerical values for the average momentum and for time." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "p0:=1;" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 30 "rho2:=simplify(subs(t=2,rho));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 82 "We borrow the code from the quantum mecha ncis example at the end of the worksheet " }{TEXT 19 7 "FFT.mws" } {TEXT -1 61 " to generate the mesh, and to carry out the DFT. We re-la bel " }{TEXT 19 1 "x" }{TEXT -1 4 " to " }{TEXT 19 1 "s" }{TEXT -1 65 " in the code in order to avoid confusion with the current use of " } {TEXT 19 1 "x" }{TEXT -1 139 " as the position in the Wigner distribut ion argument. We are particularly interested in small-size FT explorat ions. The idea is that small " }{TEXT 19 1 "s" }{TEXT -1 44 "-meshes c entered around the actual value of " }{TEXT 19 1 "x" }{TEXT -1 303 " s hould suffice. This becomes an important issue for multi-dimensional c alculations (beyond the scope of Maple, but doable in Fortran). The re ason why we have to be careful is that the DFT assumes periodicity, an d a small mesh can mean that we work with an artificial modification o f the density matrix " }{TEXT 19 3 "rho" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "ds:=0.8; n:=5; N:=2^n;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dp:=evalhf(2*Pi/(N*ds));" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "p_max:=N*dp/2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "sv:=array(1..N): pv:=array(1..N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 153 "for i from 1 to N do:\nif i < N/ 2+1 then\nsv[i]:=evalhf((i-1)*ds): sv[N-i+1]:=evalhf(-i*ds):\npv[i]:=e valhf((i-1)*dp): pv[N-i+1]:=evalhf(-i*dp): end if:\nod:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 15 "Let us try one " }{TEXT 19 1 "x" }{TEXT -1 13 " value first:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "xval :=2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "A small technical wrinkle : for " }{TEXT 19 3 "s=0" }{TEXT -1 29 " (the first value, i.e., for \+ " }{TEXT 19 5 "sv[1]" }{TEXT -1 91 ") the expression for the density m atrix is undefined (0/0); therefore, we don't substitute " }{TEXT 19 7 "s=sv[i]" }{TEXT -1 81 ", but rather take the limit. To figure this \+ out required a little debugging step." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "yR:=array(1..N): yI:=array(1..N):" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 129 "for i from 1 to N do: yR[i]:=evalf(limit(subs(x=xv al,Re(rho2)),s=sv[i])); yI[i]:=evalf(limit(subs(x=xval,Im(rho2)),s=sv[ i])); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "plot([seq([sv[ i],yR[i]],i=1..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "plot([seq([sv[i],yI[i]],i=1..N)],style=point,color=bl ue);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_s:=ds*add(yR[i]^2 +yI[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "evalhf (FFT(n, var(yR), var(yI)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_p:=dp*add(yR[i]^2+yI[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "The imaginary part should vanish; is this true?" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "dp*add(yI[i]^2,i=1..N);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 127 "Well, almost. Perhaps this repres ents a good test of the accuracy of our truncation of the transform to a low number of points." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "yRn:=array(1..N): sq:=evalf(1/sqrt(2*Pi)): for i from 1 to N do: yRn[ i]:=ds*sq*yR[i]: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "N_ p:=dp*add(yRn[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "plot([seq([pv[i],yRn[i]],i=1..N)],style=line);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "Now we set up a loop to step through " }{TEXT 19 1 "x" }{TEXT -1 1 ":" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "Tno rm:=0: for k from 1 to 25 do: xval:=evalf(-4+8/20*(k-0.5));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "yR:=array(1..N): yI:=array(1..N):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 129 "for i from 1 to N do: yR[i]:=evalf (limit(subs(x=xval,Re(rho2)),s=sv[i])); yI[i]:=evalf(limit(subs(x=xval ,Im(rho2)),s=sv[i])); od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_s:=d s*add(yR[i]^2+yI[i]^2,i=1..N);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "e valhf(FFT(n, var(yR), var(yI)));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "yRn:=array(1..N): sq:=evalf(1/sqrt(2*Pi)): for i from 1 to N do: yRn[ i]:=ds*sq*yR[i]: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 87 "N_p:=dp*ad d(yRn[i]^2,i=1..N); print(\"xval,Ns,Np=\",xval,N_s,N_p); Tnorm:=Tnorm+ 8/20*N_p;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 102 "PL[k]:=plot([seq([pv[ i],yRn[i]],i=1..N)],style=line,title=cat(\"x= \",trunc(10*xval),\"/10 \")): od: Tnorm;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "plots[d isplay]([seq(PL[i],i=1..Nk)],insequence=true);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 319 "We observe that the Wigner distribution resembles t he analytical result displayed as a surface plot over x and p. It is p ositive at all values displayed (as it should be for a Gaussian wavepa cket). The total norm preservation is very good. The use of the FFT wi th few points results in a coarse mesh of momentum values." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 257 11 "Exercise 1:" }} {PARA 0 "" 0 "" {TEXT -1 117 "Explore the parameter choices for the FF T that give reasonable results. What is more important: a fine resolut ion in " }{TEXT 19 1 "s" }{TEXT -1 24 ", or a sufficient range?" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 206 "In order to display the results for a wavepacket with tunneling we first have \+ to construct a wavepacket with such properties. One alternative would \+ be to use a FFT-based propagation scheme as introduced in " }{TEXT 19 13 "Ehrenfest.mws" }{TEXT -1 387 ". This integrator works for arbitrar y potentials. For piecewise constant potentials we can use a plane-wav e decomposition of the wavepacket. Let us develop this approach. While it is rather restrictive in applicability, it offers a link between w hat is usually done in quantum mechancis texts, namely the stationary \+ scattering theory approach with the time-dependent wavepacket approach ." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 79 "We a ssume that the potential barrier is different from zero only in the ra nge " }{TEXT 19 9 "0 < x < 1" }{TEXT -1 25 ", and that it has height \+ " }{TEXT 19 2 "V0" }{TEXT -1 72 " there, and is zero everywhere else. \+ For an asymptotic wavenumber value " }{TEXT 19 2 "k0" }{TEXT -1 79 " w e can then calculate the plane wave state. The particle mass is chosen to be " }{TEXT 19 3 "m=1" }{TEXT -1 201 " in our units. We have to co nstruct solutions to the Schroedinger equation from plane-wave segment s, i.e., we calculate energy eigenstates. The scattering is assumed to proceed from left to right. For " }{TEXT 19 3 "x<0" }{TEXT -1 138 " f or a given energy there can be a right-travelling, and a left-travelli ng plane wave (the latter describing the backscattered flux). For " } {TEXT 19 9 "0 < x < 1" }{TEXT -1 34 " we need a superposition, and for " }{TEXT 19 3 "x>1" }{TEXT -1 164 " we need only a right-travelling w ave. The waves have to be matched with continuous derivative at the tw o points where the potential changes in a step-like fashion." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "V0:=2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "q0:=k0->sqrt(k 0^2-2*V0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 119 "PW:=k->piece wise(x<0,exp(I*k*x)+ R*exp(-I*k*x), x>=0 and x<=1,A*exp(I*q0(k)*x) + B *exp(-I*q0(k)*x), x>1, T*exp(I*k*x));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 59 "In order to find an acceptable solution for a given energ y " }{TEXT 19 9 "E = k^2/2" }{TEXT -1 51 " we need to determine the co mplex-valued constants " }{TEXT 19 10 "R, T, A, B" }{TEXT -1 1 "." }} {PARA 0 "" 0 "" {TEXT -1 199 "The squared magnitudes of R and T corres pond to the reflection and transmission probabilities, because the inc ident flux has been normalized to unity (amplitude of 1 for the right- travelling wave at " }{TEXT 19 3 "x<0" }{TEXT -1 2 ")." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "k0:=1.;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 11 "E0:=k0^2/2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 129 "We are testing at an energy which corresponds to one half of the \+ barrier height. Our first attempt to formulate the condition at " } {TEXT 19 3 "x=0" }{TEXT -1 11 " reveals a " }{TEXT 19 12 "bug in maple " }{TEXT -1 22 "8. There should be no " }{TEXT 19 1 "x" }{TEXT -1 80 " left in the expression after taking the limit. We should be able to s ubstitute " }{TEXT 19 3 "x=0" }{TEXT -1 65 ", however, given that it i s included in the piecewise definition." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "C1:=limit(PW(k0),x=0,l eft)=simplify(limit(PW(k0),x=0,right));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "C1:=limit(PW(k0),x=0,left)=simplify(subs(x=0,PW(k0))) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "C2:=limit(PW(k0),x=1,r ight)=simplify(subs(x=1,PW(k0)));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 26 "Our patience with Maple's " }{TEXT 19 9 "piecewise" }{TEXT -1 86 " function has reached its limit! Maple developers have to return to th e drawing board." }}{PARA 0 "" 0 "" {TEXT -1 92 "Let us see whether an other construct available in Maple will help. Let us define explicitly :" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 127 "PW:=k->if x<0 then exp(I*k*x)+ R*exp(-I*k*x) elif x> =0 and x<=1 then A*exp(I*q0(k)*x) + B*exp(-I*q0(k)*x) else T*exp(I*k*x ) fi;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "C1:=limit(PW(k0),x =0,left)=simplify(limit(PW(k0),x=0,right));" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 128 "Ouch! So let us simply define the three adjacent piece s as separate expressions. We label them as parts 1, 2, 3 for the regi ons " }{TEXT 19 3 "x<0" }{TEXT -1 2 ", " }{TEXT 19 5 "01" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 34 "PW1:=k->exp(I*k*x)+ R*exp(-I*k*x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "PW2:=k->A*exp(I*q0(k)*x) + B*exp(-I *q0(k)*x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "PW3:=k->T*exp (I*k*x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 150 "This allows us to av oid the contentious issue in Maple of taking the limit on a piecewise \+ function. For graphing we will want a single function later." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "PW:=k->piecewise(x<0,PW1(k), x>=0 and x<=1,PW2(k), x>1, PW3(k));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "C1:=simplify(subs(x=0,PW1(k0)))=simplify(subs(x=0,PW2 (k0)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "C2:=simplify(sub s(x=1,PW2(k0)))=simplify(subs(x=1,PW3(k0)));" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 23 "Define the derivatives:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "PW1p:=unapply(simplify(diff(PW1(k),x)),k);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "PW2p:=unapply(simplify(diff( PW2(k),x)),k);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "PW3p:=una pply(simplify(diff(PW3(k),x)),k);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "Formulate the continuity of the derivative:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 63 "C1p:=simplify(subs(x=0,PW1p(k0)))=simplify(sub s(x=0,PW2p(k0)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "C2p:=s implify(subs(x=1,PW2p(k0)))=simplify(subs(x=1,PW3p(k0)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "sol:=solve(\{C1,C2,C1p,C2p\},\{A,B, T,R\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "The solution works whe n " }{TEXT 19 2 "k0" }{TEXT -1 83 " is entered as a floating point num ber. For some reason it fails in exact mode for " }{TEXT 19 4 "k0=1" } {TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 47 "To use the coefficients \+ we assign the solution." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 " assign(sol);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "plot([Re(PW (k0)),Im(PW(k0))],x=-10..10,color=[red,blue]);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 7 "PW(k0);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "Now that we know how this works, let us assemble solutions for a rang e of wavenumber values:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 " k0:=0.01: k_max:=2.5: N_k:=80: dk:=evalf((k_max-k0)/N_k);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "for i_k from 1 to N_k do: k_i:=eval f(k0+(i_k-1)*dk);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "A:='A': B:='B' : R:='R': T:='T':" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "C1:=simplify(s ubs(x=0,PW1(k_i)))=simplify(subs(x=0,PW2(k_i)));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "C2:=simplify(subs(x=1,PW2(k_i)))=simplify(subs(x=1,PW 3(k_i)));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "C1p:=simplify(subs(x=0 ,PW1p(k_i)))=simplify(subs(x=0,PW2p(k_i)));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "C2p:=simplify(subs(x=1,PW2p(k_i)))=simplify(subs(x=1, PW3p(k_i)));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "sol:=solve(\{C1,C2, C1p,C2p\},\{A,B,T,R\});" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "assign(s ol): kv[i_k]:=k_i: WF[i_k]:=PW(k_i): od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "We can graph an animation of the basis state solutions:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "#PL:=[seq(plot([Re(WF[i]) ,Im(WF[i])],x=-5..5,color=[red,blue]),i=1..N_k)]:" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 36 "#plots[display](PL,insequence=true);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "Now superimpose those plane-wave s cattering solutions:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 231 "#I P:=(f,g)->evalf[6](Int(simplify(evalc(conjugate(f)*g)),x=-infinity..0, method=_Gquad))+evalf[6](Int(simplify(evalc(conjugate(f)*g)),x=0..1,me thod=_Gquad))+evalf[6](Int(simplify(evalc(conjugate(f)*g)),x=1..infini ty,method=_Gquad));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 83 "The above \+ definition leads to slow calculation, therefore, we use a simpler meth od:" }}{PARA 0 "" 0 "" {TEXT -1 84 "The real axis is discretized and t runcated for the purpose of calculating integrals." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "x_min:=-35: x_max:=35: N_i:=200: dx:=evalf( (x_max-x_min)/N_i);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "IP:= (f,g)->dx*add(evalf(subs(x=x_min+(i-0.5)*dx,conjugate(f)*g)),i=1..N_i) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "a:=3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "phi0:=(2/Pi)^(1/4)/sqrt(a)*exp(-(x+10)^2/ a^2)*exp(1.*I*x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "int(si mplify(evalc(abs(phi0)^2)),x=-infinity..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "IP(phi0,phi0);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 68 "On the other hand, the initial state is completely locali zed in the " }{TEXT 19 3 "x<0" }{TEXT -1 68 " region, so we should be \+ able to get away with one half of the work:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "IP0:=(f,g)->dx*add(evalf(subs(x=x_min+(i-0.5)*dx,c onjugate(f)*g)),i=1..N_i/2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "IP0(phi0,phi0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "c_ n:=[seq(IP0(WF[n],phi0),n=1..N_k)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "plot([seq([k,abs(c_n[k])^2],k=1..N_k)],style=point); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 361 "The above result us smooth. \+ For different initial wavepackets (e.g., narrower ones, which have a b igger momentum spread) the mesh of wavenumber values included may requ ire modification so that the above collection of coefficients indicate s both a reasonable resolution (almost 'smooth' curve), and that the m agnitude of the highest coefficients approaches zero." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 366 "We form a superpositi on quite naively: for a discrete set of scattering states chosen for s ome arbitrary set of wavenumber values we do not have a guarantee of o rthogonality. In the limit of an infinite real axis these states acqui re delta-function normalization. We will have to normalize the wavepac ket represented in terms of the scattering states later manually." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "psi:=add(c_n[i]*(dk)*WF[i]*e xp(-I*0.5*kv[i]^2*t),i=1..N_k):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 114 "The following calculation allows to visualize the time evolution \+ of the wavepacket. It takes some time to compute." }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 93 "PLt:=[seq(plot(evalc(abs(subs(t=1.0*i_t,psi) )^2),x=-15..10,title=cat(\"t= \",i_t)),i_t=0..30)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "plots[display](PLt,insequence=true);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "IP(subs(t=0,psi),subs(t=0,ps i));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "IP(subs(t=10,psi),subs(t=10 ,psi));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "No:=1/sqrt(Re(%));" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 83 "This normalization constant is use d below to have a properly normalized wavepacket." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 117 "Now we implement the Wig ner distribution calculation. The wavepacket is ready for calculation \+ at any reasonable time " }{TEXT 19 1 "t" }{TEXT -1 94 ". First, howeve r, we test the Wigner function calculation by using the Gaussian initi al state." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "#phi:=No*subs( t=10,psi):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "phi:=phi0:" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 231 "Now map the wavefunction to a m esh on which the Fourier transforms can be carried out. First define a Fourier mesh by lifting code from the above section. The next line se rves only as a reminder to explain the relationship between " }{TEXT 19 1 "x" }{TEXT -1 36 " and the Fourier transform variable " }{TEXT 19 1 "s" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 " #rho:=expand(subs(x=x+s/2,phi)*conjugate(subs(x=x-s/2,phi))):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "ds:=0.8; n:=6; N:=2^n;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dp:=evalhf(2*Pi/(N*ds));" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "p_max:=N*dp/2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "sv:=array(1..N): pv:=array(1..N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 153 "for i from 1 to N do:\nif i < N/ 2+1 then\nsv[i]:=evalhf((i-1)*ds): sv[N-i+1]:=evalhf(-i*ds):\npv[i]:=e valhf((i-1)*dp): pv[N-i+1]:=evalhf(-i*dp): end if:\nod:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 112 "In order to be able to carry out Wigner \+ distribution calculations we need the wavefunction on a mesh of spacin g " }{TEXT 19 6 "dx=0.4" }{TEXT -1 97 ", chosen to be commensurate wit h the above Fourier transform mesh. This means that the choice of " } {TEXT 19 1 "x" }{TEXT -1 90 "-values will be restricted. This mesh has to be sufficiently large so that for all points " }{TEXT 19 1 "x" } {TEXT -1 60 " of interest the Wigner transform can be carried out by t he " }{TEXT 19 1 "s" }{TEXT -1 47 "-variable while exploring the neigh bourhood of " }{TEXT 19 1 "x" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 144 "N_x:=512: dx:=0.4: xv:=array(1..N_x): phi_x:=ar ray(1..N_x): for i from 1 to N_x do: xv[i]:=(i-N_x/2)*dx: phi_x[i]:=ev alf(subs(x=xv[i],phi)): od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "PL1:=plot([seq([xv[i],Re(phi_x[i])],i=1..N_x)],color=red):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "PL2:=plot([seq([xv[i],Im(phi_x[i])] ,i=1..N_x)],color=blue):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "plots[display](PL1,PL2,view=[-20..10,-0.6..0.6],style=point);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "dx*add(abs(phi_x[i])^2,i=1.. N_x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "The wavepacket is normal ized." }}{PARA 0 "" 0 "" {TEXT -1 89 "Now the problem of carrying out \+ the Wigner transform becomes one of bookkeeping: we pick " }{TEXT 19 1 "x" }{TEXT -1 40 "-values on the mesh, and read ranges of " }{TEXT 19 1 "s" }{TEXT -1 65 "-values into the appropriate arrays for the FT. Let us pick some " }{TEXT 19 1 "x" }{TEXT -1 14 "-values first." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "Np:=20: i_x_arr:=[seq(210+i* 2,i=1..Np)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "seq(xv[i_x_ arr[i]],i=1..nops(i_x_arr));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "L et's do a single " }{TEXT 19 1 "x" }{TEXT -1 7 "-value:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "xval:=xv[230];" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 105 "A technical point: while we assemble the array fo r the Fourier transform we have to keep in mind how the " }{TEXT 19 1 "s" }{TEXT -1 71 "-array is filled. In the previous application of the FFT algorithm the " }{TEXT 19 1 "s" }{TEXT -1 68 "-values were substi tuted. Here we just want to use look-up from the " }{TEXT 19 5 "phi_x " }{TEXT -1 43 " array. The next line serves as a reminder:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "#sv[i]:=evalhf((i-1)*ds): sv [N-i+1]:=evalhf(-i*ds):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 " yR:=array(1..N): yI:=array(1..N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 249 "i_x:=230; for i from 1 to N/2 do: rho:=evalc(phi_x[i_x+(i-1)]*con jugate(phi_x[i_x-(i-1)])); yR[i]:=evalf(Re(rho)); yI[i]:=evalf(Im(rho) ); rho:=evalc(phi_x[i_x-i]*conjugate(phi_x[i_x+i])); yR[(N-i+1)]:=eval f(Re(rho)); yI[(N-i+1)]:=evalf(Im(rho)); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "plot([seq([sv[i],yR[i]],i=1..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "plot([seq([sv[i],yI[i]],i =1..N)],style=point,color=blue);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_s:=ds*add(yR[i]^2+yI[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "evalhf(FFT(n, var(yR), var(yI)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_p:=dp*add(yR[i]^2+yI[i] ^2,i=1..N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "The imaginary part should vanish; is this true?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "dp*add(yI[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "yRn:=array(1..N): sq:=evalf(1/sqrt(2*Pi)): for i from 1 to N do: yRn[i]:=ds*sq*yR[i]: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "N_p:=dp*add(yRn[i]^2,i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "plot([seq([pv[i],yRn[i]],i=1..N)],style=line);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "Now that we know how it works, we \+ can put it into a loop, and try a time-evolved wavefunction:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "phi:=No*subs(t=10,psi):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 144 "N_x:=512: dx:=0.4: xv:=arra y(1..N_x): phi_x:=array(1..N_x): for i from 1 to N_x do: xv[i]:=(i-N_x /2)*dx: phi_x[i]:=evalf(subs(x=xv[i],phi)): od:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 70 "PlR:=plot([seq([xv[i],Re(phi_x[i])],i=1..N_x)] ,color=red,style=point):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "PlI:=plot([seq([xv[i],Im(phi_x[i])],i=1..N_x)],color=blue,style=point ): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "plots[display](\{PlR ,PlI\},view=[-15..20,-0.5..0.5]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 36 "We adjust the computational mesh of " }{TEXT 19 1 "x" }{TEXT -1 37 "-values at which the Wigner function " }{TEXT 19 6 "f(x,p)" } {TEXT -1 46 " is calculated to span interesting regions in " }{TEXT 19 1 "x" }{TEXT -1 7 "-space." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "Np:=35: i_x_arr:=[seq(220+i*2,i=1..Np)];" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 39 "seq(xv[i_x_arr[i]],i=1..nops(i_x_arr));" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 21 "In the loop over the " }{TEXT 19 1 "x" }{TEXT -1 23 "-mesh we verify in the " }{TEXT 19 5 "Tnorm" } {TEXT -1 86 " variable that the overall norm of the Wigner distributio n approximately equals unity." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "Tnorm:=0: for k from 1 to Np do: i_x:=i_x_arr[k]: xval:=xv[i_x]; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "yR:=array(1..N): yI:=array(1..N ):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 239 "for i from 1 to N/2 do: rho: =evalc(phi_x[i_x+(i-1)]*conjugate(phi_x[i_x-(i-1)])); yR[i]:=evalf(Re( rho)); yI[i]:=evalf(Im(rho)); rho:=evalc(phi_x[i_x-i]*conjugate(phi_x[ i_x+i])); yR[(N-i+1)]:=evalf(Re(rho)); yI[(N-i+1)]:=evalf(Im(rho)); od :" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "N_s:=ds*add(yR[i]^2+yI[i]^2,i= 1..N);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "evalhf(FFT(n, var(yR), va r(yI)));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "yRn:=array(1..N): sq:=e valf(1/sqrt(2*Pi)): for i from 1 to N do: yRn[i]:=ds*sq*yR[i]: od:" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 128 "N_p:=dp*add(yRn[i]^2+ds*sq*yI[i]^ 2,i=1..N); print(\"xval,Ns,Np=\",xval,N_s,N_p); Tnorm:=Tnorm+(xv[i_x_a rr[2]]-xv[i_x_arr[1]])*N_p;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 102 "PL[ k]:=plot([seq([pv[i],yRn[i]],i=1..N)],style=line,title=cat(\"x= \",tru nc(10*xval),\"/10\")): od: Tnorm;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "plots[display]([seq(PL[i],i=1..Np)],insequence=true); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "plots[display](PL[20]); " }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 258 11 "Exercise 2:" }}{PARA 0 "" 0 "" {TEXT -1 328 "Explor e the Wigner distribution using the cuts as shown above for different \+ times, and describe the behaviour. Negative parts in the distribution \+ are apparently strong when the backscattered and forward-moving parts \+ of the wavepacket meet and interfere, but also apparently behind the b arrier, where no interference takes place." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 259 11 "Exercise 3:" }}{PARA 0 "" 0 " " {TEXT -1 256 "Vary the wavepacket parameters (and the potential para meters), and explore the Wigner distribution. For wavepackets moving r apidly over the barrier the regions of negative Wigner functions shoul d be small. What happens for scattering from attractive wells?" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 168 " The Maple program can also be used to read in wavepackets produced wit h Fortran or C programs which employ an FFT-based time integrator (suc h as the one illustrated in " }{TEXT 19 13 "Ehrenfest.mws" }{TEXT -1 98 "). To explore Wigner distributions for such situations can make up an interesting student project." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "102 1 0" 20 }{VIEWOPTS 1 1 0 3 2 1804 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }