Transfer Matrix Method in a rectangular potential well

Question

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...

Answer 1

I found the mistake in my code...

It's the mass. The author of the book used only the free particle mass while I (taken with some other things I need to calculate) forgot that I was using the mass in GaAs and since I didn't know what to expect for GaAs, I thought I was getting something wrong...

The good news is the algorithm is ok.

so, just use m0 for m in the code and you'll get a nice result, like in the book (page 100)

the code is not fault proof for other potentials though, it still needs testing and tweaking! (!!!)

but it works for a one step square potential at least.

I think I tweaked my code a bit to work on every step potential, if you wish to find bound states:

I'm pinning in an extra potential function I used for testing.

%main.m

close all
clear all
clc
format long e

N=111;        
xMin = -5e-9; %nm  
xMax = -xMin;       
x = linspace(xMin,xMax, N);

%U = square_potential_well(xMin, x, N);

a = 2e-10;
b = 5e-10;
c = 8e-10;
d = 10e-10;
U = random_step_pot(x,a,b,c,d,N);

[dx, Ustep] = step_width_and_potential(U, x);

V = tmm(U, dx, Ustep, 200);

---

function U = square_potential_well(xMin, x, N)

x1 = abs(xMin)/2;         % polovina sirine
U0 = 0; %dubina jame

% inicijalizacija matrice koja predstavlja potencijalnu energiju jame
U = ones(N,1)*20;

% kreiranje pravougaone jame sirine 2*x1
for cn = 1 : N
   if abs(x(cn)) <= x1, U(cn) = U0; end; %if;
end;

% Test U
figure(1)
plot(x,U,'LineWidth',3);

end

---

function U = random_step_pot(x,a,b,c,d, N)

U = ones(N,1);
Umax = 20;
U = U*Umax;


for cn = 1:N
    if x(cn) >= a && x(cn)<=b
      U(cn) = 0;
    elseif x(cn)>b && x(cn)<=c
      U(cn) = -8;
    elseif x(cn)>c && x(cn)<=d
      U(cn) = 16;
    end
end 

figure(1)
plot(x,U, 'LineWidth', 2, 'g')
xlabel('x [nm]');
ylabel('U(x) [eV]'); 

end

---

function [dx, Ustep] = step_width_and_potential(U, x)

    flag = [];
%    flag(1)=1;
    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

    % break nekad ne radi zato dodajem j vrednost flag koju i treba da ima na kraju
    j = length(flag);

    dx = [];
    i = 1; %reinicijalizacija
    for k = 2:j
            dx(i) = x(flag(k)) - x(flag(k-1));
            i = i+1;
    end

    Ustep = [];
    for i = 1:length(flag)-1 
        Ustep(i) = U(flag(i));
    end

end

---

function V = tmm(U, dx, Ustep, countE)

    %konstante
    hbar = 1.054475e-34;     % Js
    q = 1.602e-19;         % C naelektrisanje cestice...
    m0 = 9.1093897e-31; %kg masa slobodnog elektrona
    m = 0.067*m0;    


    %k = [];
    %kappa = [];
    E = linspace(min(U) + 1e-3, U(1)-1e-3, countE); % countE je broj ispitivanih energija
    V = [];
    %Ms = zeros(2,2,length(E));

% za slucaj pravougaone potencijalne jame za svaku energiju ce biti ispunjen samo prvi if uslov i jednom ce se prolaziti kroz for jer postoji samo jedan jama
    for j = 1:length(E)
        kappa = sqrt(2*m/hbar^2 * q * (U(1) - E(j)));
        M = eye(2);
        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 %if

            M = M*Ms;

        end


        t11 = 0.5 * (M(1, 1) + M(2, 2)) - 0.5*(kappa*M(1,2) + M(2,1)/kappa);


        V(j) = t11;       
    end

    figure(2)
    plot(E, V, 'r', 'LineWidth', 2);
    xlabel('E[eV]');
    ylabel('t_{11}(E)');
    %axis([0,20,-2.5,1.5]);
    %print('t11', '-dpng');
end

---

just remember that if you wish to test the results from the book, work with the free particle mass $m_0$ by changing the value for mass $m$ in tmm.m