Solutions (13)

Assigned: April 23
Due: April 30

Setup your Notebook (in-Class)

Copy the Crank-Nicholson formula from the class notebook, 13-QuantumPotentials.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]

cnStep[a_List,b_List,v_List][ψ_List] :=  Module[{sparseA,r,len},(* a faster, pre-defined version of cnStep *)
len = Length[ψ];
sparseA = SparseArray[{Band[{1,2}] -> a, Band[{1,1}] -> b, Band[{2,1}] -> a},len];
r =  ListConvolve[{I 0.25,1.0 - I 0.5, I 0.25},ψ, 2, 0.0] - I ψ v 0.5;
LinearSolve[sparseA,r]]

(* apply cnStep multiple times *)cnStep[n_Integer][ψ_List, v_List] := Nest[cnStep[#,v]&, ψ, n];
cnStep[n_Integer,a_List, b_List, v_List][ψ_List] := Nest[cnStep[a,b,v][#]&, ψ, n];

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->Blue,Joined->True,options];p2=ListPlot[Im[ψ],PlotStyle->Red,Joined->True,options];p3=ListPlot[Abs[ψ],Joined->True,options];p4=ListPlot[-Abs[ψ],Joined->True,options];
Show[p1,p2,p3,p4,options]]

Tunneling

In this problem, you are to compute the reflection and transmission coefficient through a potential barrier. When the barrier than the kinetic energy of the incident particle, the particle "tunnels" through the barrier! Of course, this is impossible for a classical particle. Tunneling is a quantum mechanical process that is important to many phenomena, from practical semiconductor devices to understanding natural radioactivity.

Part 1: Defining the Potential Barrier (In Class)

Define your potential barrier as used in the Section entitled "Potential Barriers (2)". Demonstrate that your notebook properly computes the evolution of the quantum wavefunction as it interacts with the barrier.

Solution

nGrid=200;
v0=+1.0;
v=Table[If[-10≤x≤10,v0,0.0],{x,-(nGrid/2)+1,(nGrid/2)}];
a=Table[-I 0.25,{i,1,nGrid-1}];
b=Table[1.0+I (1.0+v[[i]]) 0.5,{i,1,nGrid}];
ListPlot[v,PlotLabel->"Potential Barrier = "<>ToString[v0],Joined->True]

k0 = 0.6;
Print["Particle Energy = ", k0^2/2];ψ0 = Table[initial[k0,nGrid/4,8][x], {x, 1, nGrid}];plotψ[ψ0, PlotRange-> All]

sol=NestList[cnStep[a,b,v][#1]&,ψ0,200];
ListDensityPlot[Abs[sol]^2,Mesh->False,ColorFunction->Hue,PlotRange->All,PlotLabel->"V0 = "<>ToString[v0]<>"; E0 = "<>ToString[((k0^2)/2)]]

Part 2: Observing Tunnelling

Change the height of the potential barrier, and observe how the changes in the wavefunction solution. Can you observe the conditions when quantum tunneling is significant? Can you see changes in the reflected wavefunction? What happens when the potential barrier height is lower than the incident energy of the particle?

Solution

v0=+0.2;
v=Table[If[-10≤x≤10,v0,0.0],{x,-(nGrid/2)+1,(nGrid/2)}];
a=Table[-I 0.25,{i,1,nGrid-1}];
b=Table[1.0+I (1.0+v[[i]]) 0.50,{i,1,nGrid}];
ListPlot[v,PlotLabel->"Potential Barrier = "<>ToString[v0],Joined->True]

sol=NestList[cnStep[a,b,v][#1]&,ψ0,200];
ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->All,PlotLabel->"V0 = "<>ToString[v0]<>"; E0 = "<>ToString[((k0^2)/2)]]

v0=+0.1;
v=Table[If[-10≤x≤10,v0,0.0],{x,-(nGrid/2)+1,(nGrid/2)}];
a=Table[-I 0.25,{i,1,nGrid-1}];
b=Table[1.0+I (1.0+v[[i]]) 0.5,{i,1,nGrid}];
ListPlot[v,PlotLabel->"Potential Barrier = "<>ToString[v0],Joined->True]

sol=NestList[cnStep[a,b,v][#1]&,ψ0,200];
ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->All,PlotLabel->"V0 = "<>ToString[v0]<>"; E0 = "<>ToString[((k0^2)/2)]]

Part 3: Quantifying Tunnelling

Integrate the probability of finding the particle within each of the three regions as a function of time: (1) to the left of the barrier, (2) within the barrier, and (3) to the right of the barrier. (Of course, the sum of all three must equal to unity!) Examine the variation of these probabilities for several values of the barrier height.

Solution

leftProb[ψ_] := Plus @@ (Abs[ψ[[Range[1,89]]]]^2);
rightProb[ψ_] := Plus @@ (Abs[ψ[[Range[111,200]]]]^2);
barrierProb[ψ_] := Plus @@ (Abs[ψ[[Range[90,110]]]]^2)

Dimensions[sol]

ListPlot[{leftProb/@sol,rightProb/@sol,barrierProb/@sol, (leftProb/@sol) + (rightProb/@sol)+(barrierProb/@sol)},Joined->True,PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]},{RGBColor[0,1,0]}},PlotLabel->"Incident, Barrier, Transmitted Probabilities",AxesLabel->{"Δt",""}]

Notice : the time required for the particle to reflect (classically) from the barrier and return to it's starting point is:

bounceTime = 2 (90-50)/k0 // Ceiling

Other examples:

v0=+0.2;
v=Table[If[-10≤x≤10,v0,0.0],{x,-(nGrid/2)+1,(nGrid/2)}];
a=Table[-I 0.25,{i,1,nGrid-1}];
b=Table[1.0+I (1.0+v[[i]]) 0.5,{i,1,nGrid}];sol=NestList[cnStep[a,b,v][#1]&,ψ0,200];
ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->All,PlotLabel->"V0 = "<>ToString[v0]<>"; E0 = "<>ToString[((k0^2)/2)]]

ListPlot[{leftProb/@sol,rightProb/@sol,barrierProb/@sol},Joined->True,PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]},{RGBColor[0,1,0]}},PlotLabel->"Incident, Barrier, Transmitted Probabilities"]

tunnelProb[v0_,bt_:bounceTime] := Module[{},
v = Table[If[-10 ≤ x ≤ 10,v0,0.0], {x, -nGrid/2 + 1, nGrid/2}];
a = Table[-I 0.25, {i, 1, nGrid-1}];
b = Table[1.0 + I(1.0 + v[[i]]) 0.5, {i, 1, nGrid}];sol = NestList[cnStep[a,b,v][#]&, ψ0, 200];
Part[rightProb /@ sol,bt]
];

Compute a tunneling transmission plot…

ListLogPlot[{Table[{(v0/((k0^2)/2)),tunnelProb[v0]},{v0,0.05,0.30,0.01}],Table[{(v0/((k0^2)/2)),tunnelProb[v0,150]},{v0,0.05,0.30,0.01}]},PlotLabel->"Transmission Probability for E0 = "<>ToString[((k0^2)/2)],AxesLabel->{"v0/E0",""},AxesOrigin->{(0.05/((k0^2)/2)),0.0001},Joined->True, GridLines -> Automatic]

When the barrier is higher than the energy of the particle, then the tunneling probabilty becomes exponentially small as the barrier height get bigger.

1 D Crystal

What about a crystal?

nGrid=300;
v0=+1.0;
vBarrier[i_Integer] := If[(Abs[i- nGrid/2] < nGrid/6) && ((Mod[i,4] == 0) || (Mod[i,4] == 1)),v0,0.0]
v=Table[vBarrier[x],{x,1,nGrid}];
a=Table[-I 0.25,{i,1,nGrid-1}];
b=Table[1.0+I (1.0+v[[i]]) 0.5,{i,1,nGrid}];
ListPlot[v,PlotLabel->"Potential Barrier = "<>ToString[v0],Joined->True]

Energy Levels

mMatrix = Table[Which[j == i, 1+ v[i], j == i+1, -0.5, j == i-1, -0.5, True, 0], {i,1,nGrid}, {j, 1, nGrid}];

{values, vectors} = Eigensystem[Evaluate[mMatrix /. v[i_] ->vBarrier[i]// N]];

ListPlot[values,PlotStyle->{PointSize[0.02]},AxesLabel->{"Mode","Energy"}]

modePlot[i_Integer]:=ListPlot[((Abs[vectors[[nGrid-i+1]]])^2),PlotRange->{{0,nGrid+1},{0,0.1}},Joined->True,Frame->True,Ticks->None,FrameTicks->None,DisplayFunction->Identity, FrameLabel -> "Mode = "<>ToString[nGrid-i+1]<>" Energy = "<>ToString[values[[nGrid-i+1]]]]

Manipulate[modePlot[n], {n,1,100,1}]

Scattering

k0 = Sqrt[2 0.2932];
Print["Particle Energy = ", k0^2/2];ψ0 = Table[initial[k0,nGrid/4,8][x], {x, 1, nGrid}];plotψ[ψ0, PlotRange-> All]

ListPlot[v,PlotLabel->"Potential Barrier = "<>ToString[v0],Joined->True]

sol=NestList[cnStep[a,b,v][#1]&,ψ0,200];
ListDensityPlot[Abs[sol],Mesh->False,ColorFunction->Hue,PlotRange->All,PlotLabel->"V0 = "<>ToString[v0]<>"; E0 = "<>ToString[((k0^2)/2)]]

leftProb[ψ_] := Plus @@ (Abs[ψ[[Range[1,100]]]]^2);
rightProb[ψ_] := Plus @@ (Abs[ψ[[Range[200,300]]]]^2);
barrierProb[ψ_] := Plus @@ (Abs[ψ[[Range[101,199]]]]^2)

Dimensions[sol]

ListPlot[{leftProb/@sol,rightProb/@sol,barrierProb/@sol, (leftProb/@sol) + (rightProb/@sol)+(barrierProb/@sol)},Joined->True,PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]},{RGBColor[0,1,0]}},PlotLabel->"Incident, Barrier, Transmitted Probabilities",AxesLabel->{"Δt",""}]

Notice : the time required for the particle to reflect (classically) from the barrier and return to it's starting point is:

bounceTime = 2 (100-50)/k0 // Ceiling

tunnelProb[e0_] := Module[{ψInitial,bTime},
k0 = Sqrt[2 e0];
(*Print["Particle Energy = ", k0^2/2];*)
bTime = 2 (100-50)/k0 // Ceiling;
ψInitial = Table[initial[k0,nGrid/4,8][x], {x, 1, nGrid}];sol = NestList[cnStep[a,b,v][#]&, ψInitial, 400];
Part[rightProb /@ sol,bTime]
];

tunnelProb[0.29]

Table[{e,tunnelProb[e]}, {e, 0.1, 1.1,0.05}]

ListLogPlot[%,PlotLabel->"Transmission Probability vs. Energy",AxesLabel->{"E0",""},Joined->True, GridLines -> Automatic]

trapProb[e0_] := Module[{ψInitial,bTime},
k0 = Sqrt[2 e0];
(*Print["Particle Energy = ", k0^2/2];*)
bTime = 2 (100-50)/k0 // Ceiling;
ψInitial = Table[initial[k0,nGrid/4,8][x], {x, 1, nGrid}];sol = NestList[cnStep[a,b,v][#]&, ψInitial, 400];
Part[barrierProb /@ sol,bTime]
];

Table[{e,trapProb[e]}, {e, 0.15, 0.60,0.01}];
ListLogPlot[%,PlotLabel->"Trapping Probability vs. Energy",AxesLabel->{"E0",""},Joined->True, GridLines -> Automatic]