I am trying to follow an algorithm that is described in Elementary Quantum Mechanics in 1D.
I want to compute eigen-energies and functions in bound states in the basic case in rectangular potential well shown on page 99 as a simple example. I am doing my code in MATLAB. And I can't get the same solution as shown in page 100. I have tried inspecting the algorithm and it seems ok to me, so I am suspecting whether I'm doing something wrong with the Units or scaling or something like that...but I can't seem to see the problem in there either!
The algorithm is inspected throughout pages 13-23. I have written the algorithm to be computed in general cases, and I have checked it with the Gaussian well too, it produces same failure results...
In summary, what I want to achieve is to get correct results for just one potential well. The algorithm in that case can be summed to this. In the general case in which you have a potential function approximated with many step functions, the transfer matrix would be:
$$T_{0, N+1} = E^{-1}(V_0; a_0)K^{-1}(V_0) M_1 M_2 \cdots M_N K(V_{N+1})E(V_{N+1}; a_{N+1})$$
where the matrices $E^{-1}, K^{-1}, E, K$ and $M$'s are given in tables for cases where the input energy $E$ that is being calculated is lesser or higher than the potential energy of the step in the specific region. (table 2 page 11 and the M matrices are given in table 3 page 16) I would love to use latex to write them here as well, but when I write something other than one-line latex code it just doesn't show (I'd appreciate a tip on that too...?)
The author of the book decomposes the equations make them simpler for computing. Generally, what needs to be done is to calculate the $M$ matrices and then pre-multiply and post-multiply them with the boundary conditions of the potential function. $M$ matrices are dependent on the value of the potential step in the according well and its width.
In the most easiest case of one rectangular potential well the regions are left and right of the well. The potential energy outside the well is 20eV while inside is 0V. The with of the well is 8 Angstroms (or 0.8nm).
I need more reputation to post more than two links, and I think this book I posted is really good for reference, so I can't post what the solution looks like but tell you that it's on page 100.
So, getting back to the problem in hand. I would like to calculate the energy eigenvalues in bound states. The criteria for bound states is that the energies that are being calculated need to be less then the energy of the potential well on it's left and right boundary region.$E<V_L,V_R$ where $V_L=V_R=20eV$ in this case. The MATLAB code makes vector E with values slightly above 0 and slightly less than 20eV.
So, since the problem for rectangular well is greatly simplified it turns out the only thing you need to compute is this (page 99):
$$ t_{11}(E) = \cos(ka) + 0.5 * \left(\frac{\kappa}{k} - \frac{k}{\kappa}\right)\sin(ka) $$
where $k = \sqrt{2mE/\hbar}$, in general case it's $E-V_i$ for it's the case on the region where $E<V_i$ but in this case the value of the well is 0.
And $\kappa = \sqrt{2m(V_{L/R}-E)/\hbar}$
Here's the code:
%main.m
close all
clear all
clc
% parameters for constructing the rectangular potential well
N=101; %number of points
xMin = -0.8e-9; % nm
xMax = -xMin;
x1 = abs(xMin)/2;
U0 = 0; % well depth
%
U = ones(N,1)*20;
x = linspace(xMin,xMax, N);
% creating the potential function U
for cn = 1 : N
if abs(x(cn)) <= x1, U(cn) = U0; end; %if;
end;
% Test U
% plot(x,U,'LineWidth',3);
% stepWidth_and_their_Pots(U, x) function returns the width dx of each steps
% through out the potential function U as well as their corresponding potential value Ustep
[dx, Ustep] = stepWidth_and_their_Pots(U, x);
numE = 50; % scope of input energy levels E to compute
% computes the transfer matrix element t11(E) and plots it
V = transferMatrix(U, dx, Ustep, numE);
function [dx, Ustep] = stepWidth_and_their_Pots(U, x)
flag = [];
j = 1;
for i = 1:length(U)-1
if U(i) ~= U(i+1)
flag(j) = i+1;
if i == length(U)-1
break;
else j = j+1;
end
end
end
j = length(flag);
dx = [];
i = 1;
for k = j:-1:1
if k-1 == 0
break;
else
dx(i) = x(flag(k)) - x(flag(k-1));
i = i+1;
end
end
Ustep = [];
for i = 1:length(flag)-1
Ustep(i) = U(flag(i));
end
end
function V = transferMatrix(U, dx, Ustep, numE)
%constants
hbar = 1.055e-34; % J.s
q = 1.602e-19; % C
m0 = 9.1093897e-31; % kg
m = 0.067*m0;
E = linspace(min(U) + 1e-3, U(1)-1e-3, numE);
V = [];
for j = 1:length(E)
kappa = sqrt(2*m/hbar^2 * q * (U(1) - E(j))); % U(1)=U(N) and both < E
for i = 1:length(dx)
if E(j) > Ustep(i)
k = sqrt(2*m/hbar^2 * q * (E(j) - Ustep(i)));
Ms = [cos(k*dx(i)), -(1/k)*sin(k*dx(i)); k*sin(k*dx(i)), cos(k*dx(i))];
elseif E(j) < Ustep(i)
k = sqrt(2*m/hbar^2 * q * (Ustep(i) - E(j)));
Ms = [cosh(k*dx(i)), -(1/k)*sinh(k*dx(i)); -k*sinh(k*dx(i)), cosh(k*dx(i))];
else
k = 0;
Ms = [1, -dx(i); 0, 1];
end
end
t11 = 0.5 * (Ms(1, 1) + Ms(2, 2)) - 0.5*(kappa*Ms(1,2) + Ms(2,1)/kappa);
%T = [exp(kappa*xMin) 0; 0, exp(-kappa*xMin)] * 0.5 * [1, -1/kappa; 1, 1/kappa] ...
%* Ms * [1, 1; -kappa, kappa] * [exp(-kappa*xMax) 0; 0, exp(kappa*xMax)];
%test = cos(k*dx) + 0.5*(kappa/k - k/kappa)*sin(k*dx);
V(j) = t11;
end
plot(E, V)
end
Can someone help?
UPDATE EDIT: I checked my code and found an error in the mass unit...it was 0.0067 instead of 0.067, which explains one dimension flaw. But I still don't get the correct values. I have checked my units again and again, and I just can't seem to understand where's the error.
so I have: $$\hbar = 1,055 * 10^{-34} Js$$ $$q = 1,602 * 10^{-19} C$$ $$m_0 = 9,1093897 * 10^{-31} kg$$ $$m = 0,067*m_0 = 6,1033 * 10^{-32} kg$$
Length scale is in $nm$ and my energies are in $eV$ Following the formula for $k$ or $\kappa$ I do get that their units are in $m^{-1}$ and that their power is $10^{8}$ which is good since that is multiplied with $nm$ so the values in the matrices $Ms$ should be ok. But I'm still not getting the same results as in the book (page 100)
Here's an update picture and I have updated my code for that stupid mistake in mass $m$.
I'm still not getting the correct results...
$$\begin{align} eqns \end{align}$$
. You don't need to put your information in a table, you can just put the cases as bullet, for example. $\endgroup$ – nicoguaro♦ Aug 20 '15 at 14:58