Scattering from Multiple Potential Barriers

Student Projects in AP1603 investigating quantum mechanics.

All projects must have…
Name:____________
Email:_____@columbia.edu

Introduction

The second general theme of our Mathematica investigations is oscillations. We have seen that simple models of classical oscilations, such as those found on strings, leads naturally to partial differential equations. These equations require careful solution on a computer to avoid numerical instability or inaccuracy. In this project, you are to go "atomic" and use Mathematica to continue your investigations of the numerical solutions to the time-dependent Schrödinger's equation.

I suggest you perform numerical experiments similar to those I presented in our classroom notebooks. Like last quarter, your goal is to produce a complete Mathematica notebook containing a summary of your investigations. I have prepared a general outline as a suggested guide to a notebook on quantum oscillations. If you wish, you can investigate other related topics (but I recommend that you speak with me before you begin programming.) Your notebook should contain Mathematica expressions and graphics which illustrate your solution. Try not to include a large number of repeated expressions. Instead, generate a table or graphic of your results. In all cases, format your notebook and include texual comments and descriptions. Your notebook need not be long. It should be interesting and arrive at a definite conclusion.

Scattering from Multiple Potential Barriers

How does a quantum ball reflect from multiple potential barriers?

For example, consider the following two potential barriers, one with multiple potential hills:

Off[General::"spell1"];
nGrid=200;
xB1=Floor[0.6 nGrid];xB2=Floor[0.7 nGrid];dB=Floor[0.02 nGrid];
v1=Table[Which[xB1≤x≤xB1+dB,0.1,True,0.],{x,1,nGrid}];
v2=Table[Which[xB1≤x≤xB1+dB,0.1,xB2≤x≤xB2+dB,0.2,True,0.],{x,1,nGrid}];
ListPlot[v1,PlotLabel->"Single Potential Barrier",Joined->True,PlotRange->{0,0.2}]
ListPlot[v2,PlotLabel->"Multiple Potential Barriers",Joined->True,PlotRange->{0,0.2}]

For particles moving to the right with initial momentum less than about 0.45, the first barrier height is larger than the kinetic energy of the particle. The particle can tunnel through the barrier, but classically, it is reflected. 
When the second  barrier is added, normal modes can exist. The presense of these modes influence significantly the coupling of quantum particles through the first barrier. 
Question:  As a function of initial momentum (or energy), what is the reflection coeffiecient from a single barrier and a multiple barrier? Can you suggestion a reason for the difference?

To proceed, I suggest the following three steps.

Step 1: Preparing the Notebook

First, you must prepare your notebook so that you are able to compute solutions to Schrodinger's equations on a grid, including a variable potential, V(x). 
Of course, we must use the Crank-Nicolson method to solve the finite-difference approximation to the continuous partial differential equation with good accuracy. We set the grid-spacing to unity in both the time and spatial directions.

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[ψ];
r =  ListConvolve[{I 0.25,1.0 - I 0.5, I 0.25},ψ, 2, 0.0] - I ψ v 0.5;
sparseA = SparseArray[{Band[{1,2}] -> a, Band[{1,1}] -> b, Band[{2,1}] -> a},len];
LinearSolve[sparseA,r]]

cnStep[n_Integer][ψ_List, v_List] := Nest[cnStep[#,v]&, ψ, n]; (* apply cnStep multiple times *)
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];p2=ListPlot[Im[ψ],PlotStyle->Red,Joined->True];p3=ListPlot[Abs[ψ],Joined->True];p4=ListPlot[-Abs[ψ],Joined->True];Show[p1,p2,p3,p4,options]]

Step 2: Example Solutions

Use cnStep[…] to simulate quantum particles bouncing from the barriers as a function of the initial momentum. Remember, you are trying to determine what is the probability that an object will reflect from the barrier.

A v1 and v2 example

a = Table[-I 0.25, {i, 1, nGrid-1}];
b1 = Table[1.0 + I(1.0 + v1[[i]]) 0.5, {i, 1, nGrid}];
b2 = Table[1.0 + I(1.0 + v2[[i]]) 0.5, {i, 1, nGrid}];

Below, we set the initial momentum to be p0 = 0.325, corresponding to an initial energy of 0.325^2/2 = 0.053.

ψ0 = Table[initial[0.325,nGrid/4,14][x], {x, 1, nGrid}];
plotψ[ψ0]

ψ1 = cnStep[a,b1,v1][ψ0];
plotψ[ψ1]

sol1=NestList[cnStep[5,a,b1,v1][#1]&,ψ0,100];
ListDensityPlot[Abs[sol1],Mesh->False,ColorFunction->Hue,PlotRange->All]

sol2=NestList[cnStep[5,a,b2,v2][#1]&,ψ0,100];
ListDensityPlot[Abs[sol2],Mesh->False,ColorFunction->Hue,PlotRange->All]

The probability of finding the particle to the left of the barrier, x ≤ xB1, is given by the following.

probLeft[ψ_List] := Plus @@ Part[Chop[Conjugate[ψ]  ψ],Range[xB1]]

p1=ListPlot[1-probLeft/@sol1,PlotRange->{0,0.4},PlotStyle->{Red}];
p2=ListPlot[1-probLeft/@sol2,PlotRange->{0,0.4},PlotStyle->{Blue}];
Show[p1,p2,AxesLabel->{"time","Transmission Probability"}]

Step 3: Eigensystem Analysis (Optional)

Use Eigensystem to investigate how the energy levels (frequencies) and spatial structures of the normal modes of the box change for the two potential wells. What are the normal modes of v1 and v2? Does the presence of a normal mode influence the reflection coefficient when the particle's kinetic energy equals the normal modes of the multiple well?

Hint: We know that  the normal modes of a square potential well (with steep barriers) are sine and cosine functions. For the multiple well case, v2, the width of the well is xB2 - xB1 - dB = 0.08 nGrid. The lowest energy mode occurs with a momentum, k = π/(0.08 nGrid) = 39.27/nGrid. For nGrid = 200, this gives k = 0.19 and an energy of 0.019. Higher energy levels increase as squares of the integers, 4, 9, 16, …, or 0.076, 0.171, 0.304, …. Clearly, only the lowest two modes are true modes of the multiple well. (Higher energies can not be trapped by a barrier with a height of only 0.1!) When the incident wave function is near a resonance, there will be greater tunneling. The eigenvalue analysis allows us to determine the precise energy corresponding to the modes of the multiple barriers.

Defining the Eigenvalue Matricies

From the optional class notebook, 8-QuantumModes.nb, we define the matrix mMatrix expressing Schrödinger's equation in matrix form.

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}];

Example Eigenvalue Analysis

Once mMatrix is defined, we can find the "normal" modes using the built-in command Eigensystem[…].

vb[i_Integer] := v2[[i]];{values, vectors} = Eigensystem[Evaluate[(mMatrix /. v[i_] ->vb[i])// N]];

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

The "trapped" modes of the multiple well can be seen by plotting the eigenvectors:

modes=Table[ListPlot[((Abs[vectors[[i]]])^2),PlotRange->{{0,nGrid+1},{0,0.10}},Joined->True,PlotStyle->{Hue[((30-i)/5)],Thickness[0.02]},Frame->True,Ticks->None,FrameTicks->None, PlotLabel -> "Mode = "<>ToString[nGrid-i+1]],{i,nGrid,nGrid-19,-1}];

Manipulate[modes[[n]], {n, 1, 19,1}]

{values[[nGrid - 9]], values[[nGrid - 19]]} (* the energies of the first two modes of the multiple well *)

The actual energy values are close to (but less than) expectations.

Step 4: Plot your results

Finally, summarize your results by plotting the transmission (or reflection) coefficient as a function of p0. Remember to look only after the first reflection. If you can, compare the two potential wells v1 and v2.