Execrises (12)

Assigned: April 10
Due: April 17

Setup your Notebook (in-Class)

Copy the Crank-Nicolson formula from the class notebook, 7-QuantumIntro.nb, that advances the wavefunction in unit time steps. 
Also, copy the initial condition function that creates an initial wave-packet of a given size and momentum.
These functions are called cnStep[…] and initial[…][x].
(If you wish, you may also want to copy the function plotψ[…].)

Solution

cnStep[ψ_List, v_List] := Module[{sparseA,r,len},
len = Length[ψ];
sparseA = SparseArray[{{i_,i_} :> 1 + I (1.0 + v[[i]]) 0.5, {i_,j_}/; Abs[i-j] ==  1 -> - I 0.25}, {len, len}, 0.0];
r =  ListConvolve[{I 0.25,1.0 - I 0.5, I 0.25},ψ, 2, 0.0] - I ψ v 0.5;
(* fixed boundaries *)
sol = LinearSolve[sparseA,r];
sol]

initial[k_,x0_,s_][x_] := Exp[I k x] Exp[-((x-x0)/(2 s))^2]/Sqrt[ Sqrt[2 π s^2]];

plotψ[ψ_List,options___]:=Module[{p1,p2,p3,p4},p1=ListPlot[Re[ψ],PlotStyle->RGBColor[0,0,1],Joined->True];p2=ListPlot[Im[ψ],PlotStyle->RGBColor[1,0,0],Joined->True];p3=ListPlot[Abs[ψ],Joined->True];p4=ListPlot[-Abs[ψ],Joined->True];Show[p1,p2,p3,p4,options]]

A stationary start (in-Class)

Examine quantum wave-function dynamics when the initial condition has no initial momentum. This is determined by initializing your solution with k = 0.0, or initial[0.0,x0,s0][x].

Solution

nGrid = 100;ψ0 = (Table[initial[0.0,nGrid/2,8][x], {x, 1, nGrid}] // N);
v0 = Table[0.0, {x, 1, nGrid}];
plotψ[ψ0]

Here, I introduced an auxiliary function, called tenSteps[ψ], that returns a time-advanced wave-function resulting from 10 time steps. This reduces the memory requirements for subsequent graphics.

tenSteps[ψ_List] := Nest[cnStep[#,v0]&, ψ, 10];

(* 100 x 10 = 1000 steps! *)(sol = NestList[tenSteps, ψ0, 100];)//Timing

ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->{0,0.3}]

ListPlot3D[Abs[sol],PlotRange->{0,0.3},Mesh->False]

plotψ[sol[[80]]]

Although the particle starts in the middle of the box and although it is not "moving" in the classical sense, over time, the probablity of finding the object is no longer at the center. In fact, the probability oscillates in time between the center and the sides.

Most likely position (in-class)

Compute the "most likely position" of the object as a function of time for the initial condition above and for the initial condition used in class.
The most likely position is defined as the sum, (∑_(x=1)^nGrid)(SuperStar[ψ]) x ψ, where (SuperStar[ψ])is the complex conjugate of ψ.

avgx[ψ_List] := Plus @@ (Conjugate[ψ] Table[x, {x, 1, nGrid}] ψ) // Chop;

Solution

avgx[sol[[1]]]

avgx[sol[[100]]]

ListPlot[avgx/@sol,PlotRange->{1,100}]

nGrid = 100;ψ0 = (Table[initial[0.6,nGrid/2,8][x], {x, 1, nGrid}] // N);
v0 = Table[0.0, {x, 1, nGrid}];
plotψ[ψ0]

sol = NestList[tenSteps, ψ0, 100];

ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->{0,0.3}]

ListPlot[avgx/@sol,PlotRange->{1,100}]

Most likely momentum (in-class)

For the wave-function computed in class, determine the "most-likely momentum" as a function of time.

The most likely momentum can be calculated in several ways. One easy method is to compute the sum (∑_(x=1)^nGrid)(SuperStar[ψ])(- I (∂_x) ψ). This can be written as

avgk[ψ_List] := - I Chop[ Plus @@ (Conjugate[ψ] ListConvolve[{0.5, 0, -0.5}, ψ, 2, 0.0])];

The command ListConvolve[…] is used to find the central-difference approximation to the first derivative in x. (∂_x) ψ≈ (ψ(x + 1) - ψ(x - 1)) / 2.

Solution

avgk[sol[[1]]]

ListPlot[avgk/@sol]

Momentum Probability

Another way to view the momentum is to transform the solution to a "momentum representation". This is simply the Fourier transform (along the x direction) of the solution, ψ(x, t) -> ψ(k, t).

As described in Mathematica's help manuals, Fourier[…] returns the spectrum in a mixed order. The first nGrid/2 points correspond to k = 0 to -π. The next points correspond to π back down towards 0. For this reason, we need t o re-arrange the spectrum:

rearrangeFourier[p_List] := Reverse[p[[Range[nGrid/2 + 1, nGrid]]] ~ Join~ p[[Range[1,nGrid/2]]]]

kTest=0.6;ListPlot[Abs[Fourier[Table[(E^(I kTest x)),{x,1,nGrid}]]],PlotRange->All]

ListPlot[Abs[rearrangeFourier[Fourier[Table[(E^(I kTest x)),{x,1,nGrid}]]]],PlotRange->All]

solk = rearrangeFourier[Fourier[#]]& /@ sol;

ListDensityPlot[Abs[solk],Mesh->False,ColorFunction->Hue,PlotRange->All]

ListPlot3D[Abs[solk],PlotRange->All,Mesh->False]

Notice that zero momentum is at position nGrid/2. The particle bounces back and forth. After a few bounces, however, it becomes difficult to know whether the particle is moving to the right or to the left!

Heisenberg's Uncertainty

In this problem, you are to examine Heisenberg's Uncertainty Principle. If you specify the location of a particle with high accuracy, you can not know its momentum very well. On the other hand, if you allow uncertainty in a particle's location, then you may specify its momentum more accurately.

Part 1

In order to illustrate this, we'll look at particles initialized with "zero" momentum. We'll first examine the time-dependence of the wave-function for a "narrow" (well-localized) initial condition and next for a broad (not very-well localized) initial condtion.

Let these two initial conditions be called, ψNarrow and ψBroad:

nGrid = 100;
v0 = Table[0.0, {x, 1, nGrid}];ψNarrow = N[Table[initial[0.0,nGrid/2,2][x], {x, 1, nGrid}]];
ψBroad = N[Table[initial[0.0,nGrid/2,16][x], {x, 1, nGrid}]];

Advance the wave-function from both of these initial conditions using cnStep[…]. Describe in words the differences between the two solutions. What do these differences say about the uncertainty relation?

Part 2

Another way to illustrate the uncertainty relation is to look at both the "width" of the particle's probability in both coordinate space, or x-space, and in "momentum space", or k-space. Call the width in coordinate space Δx. Call the width in momentum space Δk. The uncertainty principle states that Δx Δk is a constant.

Below, two functions are defined that compute the width in coordinate space and the width in momentum space:

deltaX[ψ_List] := Module[{avg},
avg = Plus @@ Chop[Conjugate[ψ] Table[x, {x, 1, nGrid}] ψ];
Sqrt[Plus @@ Chop[Conjugate[ψ] Table[(x - avg)^2, {x, 1, nGrid}] ψ]]];
deltaP[ψ_List] := Module[{p,avg},
p = rearrangeFourier[Fourier[ψ]];avg = Plus @@ Chop[Conjugate[p] Table[k, {k, 1, nGrid}] p];
Sqrt[Plus @@ Chop[Conjugate[p] Table[(k - avg)^2, {k, 1, nGrid}] p]]];

These two functions work in the same way, except deltaP[…] uses the Fourier transform of ψ while deltaX[…] uses ψ(x).

{deltaX[ψNarrow], deltaP[ψNarrow], deltaX[ψNarrow] * deltaP[ψNarrow]}

{deltaX[ψBroad], deltaP[ψBroad], deltaX[ψBroad] * deltaP[ψBroad]}

Your assignment: Use the functions defined above and plot Δx versus Δk for various initial conditions. Use the initial[…] function to change the width of ψ keeping the mean momentum at zero. ListPlot your results.