A microscopic model of resonant double-barrier tunneling in a quantum system

Shanna Crankshaw
Department of Physics, University of Florida
Gainesville, Florida 32611

A simple model of resonant electronic transport through a double-barrier

structure was developed. For each region of the system, analytic solutions to
Schrdingers equation were obtained, and these solutions formed the basis of the
Matlab code written to output plots of resonant peaks based on varying barrier
parameters. Resonances in transmission probability and in the electron density
inside the quantum well were then considered for a number of these different
parameters. Also, equilibrium and non-equilibrium cases were modeled for the
same set of parameters; the observed behavior was somewhat dissimilar. As an
application of the model to device physics, a brief discussion of the concepts
governing resonant tunneling diodes is also included.

A familiar fact in the realm of physics is that on very small size scales, some of the

tidiness of the classical, macroscopic world collapses into that of quantum mechanics.

This certainly holds true for sizes on the order of magnitude of fractions of nanometers

the size of a Bohr radius. In this quantum domain, we consider the wave nature of the

electron to describe its interaction with a series of potential barriers. The model

developed is a simple one: a one-dimensional, piecewise-constant system consisting of

two square barriers. Nonetheless, this model does demonstrate a number of interesting

features and may be extended to any length series of quantum wells.

My project was to develop code in Matlab to simulate the behavior of electrons

tunneling through a double-barrier structure. Two features form the basis of the

behavior of the model: transmission probability and well density. From a consideration

1
of the transmission probabilities for a number of different electron energies, we calculate

the current through the two barriers. From a comparable consideration of density as a

function of incident energy, we calculate the total density of electrons inside the quantum

well.

BACKGROUND & THEORY

To introduce the situation, first consider the typical problem of a particle in a one-

dimensional box. The behavior of the particlesay, an electronin this situation follows

from the solution of Schrdingers equation in one dimension, which gives rise to n

discrete energy levels with n2 spacing. As long as this idealistic box has infinite walls, the

electrons wavefunction will be zero outside the box and purely sinusoidal inside. If,

though, this box is modified so as to have walls of finite height, the particle has a chance

to escape from, or tunnel out of, the confining potentials. High potentials relative to the

energy of the particle, and, likewise, wide walls compared to the width of the well,

correspond to less chance of escape, and hence a longer lifetime of the boxs bound

eigenstate.

Now, instead of picturing a particle sitting inside such a potential well, consider a

coherent beam of electrons projected

upon two such barriers of finite

height and width, with the energy of

the electrons lower than that of the

barriers (see Fig. 1). In this one-

FIGURE 1. Diagram of the series of potential barriers
dimensional picture, when the defining the parameters used.

2
incident flux encounters a potential barrier, it has only two options: reflection or

transmission. The probability of transmission through the series of barriers peaks

sharply at particular incident energies. The exact energy of these resonant modes

depends on the height and width parameters1,2 describing the barriers and on the length

of the well between them.

Obtaining the transmission coefficients themselves involves matching the

solutions, for each region of interest, to the one-dimensional, time-independent

Schrdinger equation in position space:

D2 d 2
2m dx 2 + V ( x ) ( x ) = E ( x ) (1)

where m is the mass of the electron, and D = 1055

. 1027 in cgs units [1]. We are assuming

a piecewise constant potential in this model. In those regions where V(x)=0, the equation

simplifies to

D 2 d 2
= E
2m dx 2 (2)
d 2 2mE
2 = = k 2
dx D2

2mE
where k = . We thus have a second-order differential equation to which the two
D2

linearly independent solutions are purely oscillatory, and hence the general solution may

be written as ( x ) = A sin( kx ) + B cos( kx ) , or, more conveniently for our purposes,

( x ) = Aeikx + Be ikx . Similarly, for a nonzero constant potential of magnitude V0,

Schrdingers equation becomes

1 All Matlab code is included in Appendix B

2 See Appendix B, parameters.m

3
 D 2 d 2
= ( E V0 )
2m dx 2 (3)
d 2 2m( E V0 )
2 = = q 2
dx D2

2m(V0 E )
with q = , assuming E < V0 as is the case for the barriers here [1]. For this
D2

case, the two linearly independent solutions are rising and falling exponentials and the

corresponding general solution is ( x ) = Ce qx + De qx . Considering all the regions of the

given model, then, and assuming that Vm < E , we have

Aeikx + Be ikx , x < a

qx qx
Ce + De , a < x < b
(4)
( x ) Feikm x + Ge ikm x , b < x < c ,
px px
He + Ie , c < x < d
Jeikx + Ke ikx , x > d

2m(V1 E ) 2m(V2 E )
here defining q = and p = so as to explicitly allow the left
D2 D2

and right potential barriers to have different magnitudes.

The last idea that we need to develop in order to calculate the transmission

coefficients through our series of barriers is that of probability current. Briefly stated, we

can think of the probability of finding an electron in a particular spatial region changing

due to a probability flow, or current, entering or leaving that region. This current we

define as 3

D *
jx = * (5)
2mi x x

where * ( x) is the complex conjugate of the wavefunction (x) .

3 For a derivation of the probability current from Townsend, see Appendix A

4
 jtrans , x > d
Finally, we define the transmission coefficient as T = [1]. Since in our
j inc , x< a

model the potential before and after the series of barriers is the same--namely, V=0 for

both x < a and x > d --T simplifies to a ratio of the squared moduli of the corresponding

coefficients from the original wavefunctions:

2
H
T= 2
for electrons incident from the left (+k states), and (6)
A

2
B (7)
T= 2
for electrons incident from the right (-k states).
K

These unknown variables become numerically meaningful by satisfying the necessary

boundary conditionsnamely, that the wavefunctions and their first derivatives be

continuous across each interface. These conditions generate two equations per boundary,

or four per barrier. Thus, for our system of two barriers, we must thus satisfy eight

boundary equations. To solve such systems analytically poses a straightforward but

algebraically daunting undertaking, even for a single barrier. This, then, is where Matlab

enters the picture.

The attached Matlab m-file resonant.m4 calculates the probability of

transmission for a single input energy. To demonstrate the n2-dependence of the

resonant energy levels, we plot these transmission coefficients for many different

energies and note the location of the peaks. Generally speaking, if the two barriers are

symmetrichaving equal widths and equal heightsthen the probability of transmission

through both barriers at resonance approaches 100%. The higher the barriers are, the

more like delta functions the transmission peaks become. Another way to think of this

5
feature is in terms of filters: tall barriers are more selective of the energy of electrons that

they will let through than lower barriers. Another general characteristic is that the first

resonant peak occurs at a lower energy for wider wells than for narrow ones (see Fig. 2).

To understand this statement, consider that the uncertainty in position of a particle in a

very narrow well is less than the uncertainty for one in a wide well. To compensate for

this and obey the

constraints of

Heisenbergs notorious

uncertainty principle,

the allowed

momentaand hence

energiesof the more

confined particle are

greater than the

corresponding
FIGURE 2. a) Transmission probability vs. incident energy for
momenta and energies asymmetric barriers; b) for symmetric barriers; c) for tall barriers d) for a
narrow well.
of a particle in a wider

well. As one physics professor once described, A particle does not like to live in a

narrow potential well. [2]

A peak in transmission probability means that the number of electrons tunneling

through the barriers achieves a relative maximum. We define the minimum incident

energy as min and the maximum incident energies (one for the +k and one for the -k

states) L and R (also known as the electrochemical potentials, or quasi-Fermi levels).

4 See Appendix B

6
By integrating all the probabilities of transmission over this energy range of interest

( min < E < L , R ), we obtain the current through the structure:

L ,R
J l ,r = min
T ( E )dE (8)

If the transmission peaks were true delta functions, the current would be a step function,

with one step at every resonance. Electrons that penetrate the potential barriers at these

resonant energies populate the space inside the well; therefore, a maximum in the density

of electrons in the well occurs for every maximum probability in transmission. We may

2
in fact call this density the probability density, and define the quantity dx ( x ) as the

probability of finding the electron between x and dx , where ( x ) describes the

electrons wavefunction inside the well. For a single energy, then, we obtain the well

density5 n( E ) by integrating from x = b to x = c :

c 2
n( E ) = b
dx ( x ) (9)

Repeating this procedure for a number of k values reveals peaks in the same locations as

the transmission resonances, as expected.

The shape of the transmission and density peaks themselves fit Lorentzian

functions, of the form

2
L( E ) = (10)
( )
2
( E Eres ) + 2
2

where , describing the half-width of the peak, is a parameter related to the "escape

frequency" of an electron in the well [3]. This frequency is greatest for electrons closest to

the relevant electrochemical potential L or R , and so the peaks are widest for these

7
higher energies. Thus, as Fig. 1 demonstrates for the transmission peaks, the broadest

peaks correspond to the highest resonant modes.

Any nonzero density of electrons inside the well implies a Coulomb interaction

between the bound particles, generating a repulsive potential that interacts with the

existing well depth. Before we consider the interaction itself, though, we must calculate

the density contribution from every k value incident upon the system. We determine this

total density6 by taking the integral of all the individual n values with respect to energy,

up to the given electrochemical potential:

(11)
N= min
n( E )dE

The well potential Vm rises linearly as a function of this density according to the

expression

Vm = V0 + c N (12)

where c is a constant that may characterize the properties of the particular material or

source of the potential barriers, and V0 is the well height before considering any

Coulomb interactions. Accounting for this dynamic potential, however, complicates the

situation in a nontrivial manner. As the depth of the well decreases, the electronic

wavefunctions change and resonance locations increase accordingly. As the energy of

each resonant level increases, the highest modes will simply pop out the top of the

well, and the wavefunctions will adjust to the newly perceived potentials (see Fig. 3).

This alteration of the density of electrons inside the new well changes the well depth,

which changes the effective density, and so the picture continues. The issue then

5 See Appendix B, resdens.m

6 See Appendix B, nvsvm.m

8
 FIGURE 3. As the height of the potential well increases, the locations
of the resonant levels increase accordingly.

becomes a matter of determining a self-consistent model for this systemthat is, finding

its stable points of operation.

RESULTS

We had expected that as the well potential increases and pushes out resonant modes, the

electron density would simply decrease. The observed case proved somewhat less

obvious than this. If both the + k and the k electrons are filled to the same chemical

potential (from Fig. 1, if L = R ), then there is a true Fermi level E f for the system

namely, E f = L = R . The location of the transmission peaks in this double-barrier

system depends on the parameters of both the barriers; switching the parameters of the

first and second barriers has no effect on the tunneling probability. Thus, if L = R ,

9
exactly the same number of electrons tunnels through the barriers from the left as from

the right, so we obtain no net current through the system.

If we plot the individual +k or -k currents individually, or the total, integrated

density N as a function of Vm , we find that the levels actually increase initially, then

decrease smoothly as expected (see Fig. 4). This feature we explain by recalling two

previously mentioned characteristics: firstly, that the width of the resonant peaks (for T

and for n alike) increases as the

electron's energy increases, such

that the highest modes have the

broadest peaks; secondly, that as

the depth of the potential well

decreases, the locations of

resonant energies increase

accordingly. As we integrate T

and n to obtain J and N,

respectively, we thus obtain

current and density contributions

from broader and broader peaks

as Vm increases. When the

increasing height of the well

finally does push out a resonant FIGURE 4. Plots for a 6.5 well. (a) The total density
decreases smoothly as a function of the well potential. (b)
The current for left-going and right-going electrons is
level completely, we equal, so there is no net current at equilibrium.

10
correspondingly observe a significant drop in the current and total density. This

increase-then-decrease behavior repeats until no resonances remain in the well, at which

time the levels simply decrease.

For the non-equilibrium case, we cannot discuss a true Fermi level of the system,

but we can discuss the local quasi-Fermi levels F + = L and F = R [3]. When these

two are not equal, the behavior of the integrated density and currents proves slightly

more complicated than at equilibrium. Namely, the number of bumps in the curve

increases with well length, a direct result of the fact that there we only have local, quasi-

Fermi levels and not one true one for this non-equilibrium situation. For each bulge in

the equilibrium plots, we find two for non-equilibrium, corresponding to the differing

contributions from the two unequal electrochemical potentials.

The total number of protrusions, for equilibrium and non-equilibrium alike,

depends on the length of the quantum well (see Fig. 5). This result follows from the

previous consideration that longer well lengths correspond to modes of lower energies,

such that the number of resonant peaks in the longer wells is greater than that for the

shorter ones even as Vm increases. The self-consistent operating point of the system

occurs where we can satisfy Eq. (12)graphically, where plot of N vs. Vm and Vm vs. Vm

would intersect.

In all this Matlab modeling, we did not come away wholly unscathed from the

specter of numerical inconsistencies. In particular, a close inspection of Fig. 5 reveals

very small oscillations most pronounced when Vm is close to zero. These pesky

oscillations are a consequence of the sharpness of the first resonant peak and the

finiteness of the energy increment over which we integrate to obtain the current and

11
FIGURE 5. The behavior of the density and current becomes more complex as length of the quantum
well increases. Plots of density and current shown are for well lengths of (a,b) 4, (c,d) 6.5, and (e,f) 9 angstroms.

total density. Matlab works very well with matrices and not as well with continuous

variables, such as we would ideally have the allowed incident energies. We use the

trapezoidal method of integration in the actual Matlab calculations. However, since the

12
first resonant peak is extremely sharp--sharpest for low Vm values--we sample a relatively

small number of points in this first peak. The height of this peak itself thus does not

consistently and accurately represent the true theoretical value, and we see an oscillation

in the height of the peaks with a period exactly equal to the increment of energy used set

in the code parameters. Smaller energy increments, though, correspond to more accurate

sampling of the resonant peaks; therefore, the amplitude and range of the oscillations

becomes markedly smaller for even a change of a factor of five in the energy step size.

On the other hand, the computational speed slows down just as noticeably.

APPLICATIONS

So far, we have discussed a number of features of this double-barrier system

without really considering the actual, physical source of the potentials, or the possible

applications of such a system. We now turn our considerations in this direction and to a

mesoscopic size scale to discuss briefly resonant tunneling diodes, whose basic principles

of operation rely on the type of model we have just described. The operation of most

mesoscopic conductors depends on the formation of a very thin conducting layer of

electrons formed at a pn-junction like AlGaAs/GaAs. In layers of this nature, electrons

flow laterally along the heterojunctions [3]. In resonant tunneling diodes, however,

current flows transversely to the thin layer, in a vertical configuration (see Fig. 6).

13
 FIGURE 6. Resonant tunneling device. A GaAs layer a
few nanometers thick is sandwiched between two AlGaAs
barrier layers of similar thickness. Source: Supriyo Datta,
Electronic Transport in Mesoscopic Systems. (Cambridge
University Press, New York, 1995), p. 247.

We may think of the structure as a filter permitting only electrons of a particular energy

transmit, as mentioned earlier. When we apply a bias voltage, we lower the resonant

energies relative to the incident energies of the electrons from the emitter (here, the top

layer of GaAs). If we continue to increase the applied voltage past some threshold level

VT , the resonant level falls below the range of the emitters incident energiesin other

words, the energy level of the resonance is below that of the conduction band in the

emitter. Since the filter energy is lower than that of the conduction band, a bias

exceeding this threshold voltage effects a sharp drop in the observed current. At this

point, the system exhibits a negative differential resistance [3].

CONCLUSIONS

Although the model we considered in this project comprises a number of

significant physical simplifications, we have thoroughly understood it and have

demonstrated numerically some general, known features of double-barrier tunneling,

such as the behavior of the resonances, as well as some less-anticipated ones like the

14
shape of the N vs. Vm curves. Further efforts to relax some of the assumptions used in the

development of this model certainly suggest a route of future code development.

ACKNOWLEDGEMENTS

I would like to extend many thanks to Dr. Selman Hershfield, my research

mentor, for all of his efforts this summer. Thanks are also in order for Drs. Kevin

Ingersent and Alan Dorsey, the director and associate director of the physics Research

Experience for Undergraduates (REU) program at the University of Florida, and Darlene

Latimer, whose coordination efforts facilitated this program. I also make a grateful

mention of the National Science Foundation for its sponsorship of these summer REU

programs.

15
APPENDIX A: DERIVATION OF THE PROBABILITY CURRENT
Townsend

*
To see how the probability current arises, consider , the time rate of change of the
t

probability density. Using the time-dependent Schrdinger equation, we see that the

time derivative of the wavefunction is given by

( x, t ) 1 D 2 2 ( x, t )
= + V ( x) ( x, t )
t iD 2m x 2

Therefore

* ( x, t ) 1 D 2 2 * ( x, t )
= + V ( x) * ( x, t )
t iD 2 m x 2

and

* *
= + *
t t t
D *2 2
* D 2 2
= + V ( x)
* + + V ( x )
iD 2m x 2 iD 2m x
2

D * * D
2 2 2 2
= +
iD 2m x 2 iD 2m x 2

This equation can be expressed in the form

* j
= x
t x

where

D *
jx = *
2mi x x

This reflects a local conservation lawof probability. For example, if we integrate

between x = a and x = b , we obtain

16
 d b
dt a
dx * = j x (b, t ) + j x (a, t )

If the net probability of finding the particle between a and b increases, it does so because

of a net probability current flowing into the region either at a [positive current flows in

the positive x direction and hence j x (a, t ) > 0 means inward flow at a ] or at b [negative

current means current in the negative x direction and hence j x (b, t ) < 0 means inward

flow at b ]. Thus the probability in a region of space increases or decreases because there

is a net probability flow into or out of that region [1].

17
APPENDIX B: MATLAB M-FILES

parameters.m
resonant.m
resdens.m
resonantcomp.m
lengthen.m
nvsvm.m

18
%parameters.m
%created 29/06/01
%last modified 26/07/01

%Defines barrier parameters for resonant tunneling problem

%cgs units

hbar=(6.6260755/(2*pi))*10^-27; %units ergs

qe=1.60217733e-12; % (qe) ergs = 1eV
Emin=1e-14;
Eincr=1e-14;

a=0e-8; %x placement of first edge of first barrier

w1=.53e-8; %width of first barrier
L=9e-8; %length of well
w2=1.5*.53e-8; %width of second barrier
b=a+w1; %x placement of second edge of first barrier
c=b+L; %x placement of first edge of second barrier
d=c+w2; %x placement of second edge of second barrier

V1=15*qe; %height of first barrier

V2=15*qe; %height of second barrier

m=9.1093897*10^-28; %mass of electron

19
%resonant.m
%created 29/06/01
%last modified 26/07/01
%written by Shanna Crankshaw

%finds the transmission coefficient for one incident energy on a series

of two
% square barriers
%finds the scattering matrix

k=sqrt(2*m*E/hbar^2);
q=sqrt(2*m*(V1-E)/hbar^2);
p=sqrt(2*m*(V2-E)/hbar^2);

m1=[exp(i*k*a) exp(-i*k*a); i*k*exp(i*k*a) -i*k*exp(-i*k*a)];

m2=[exp(q*a) exp(-q*a); q*exp(q*a) -q*exp(-q*a)];
m3=[exp(q*b) exp(-q*b); q*exp(q*b) -q*exp(-q*b)];
m6=[exp(p*c) exp(-p*c); p*exp(p*c) -p*exp(-p*c)];
m7=[exp(p*d) exp(-p*d); p*exp(p*d) -p*exp(-p*d)];
m8=[exp(i*k*d) exp(-i*k*d); i*k*exp(i*k*d) -i*k*exp(-i*k*d)];

%m23=m2*inv(m3);
m23=[cosh(q*w1) (-1/q)*sinh(q*w1); -q*sinh(q*w1) cosh(q*w1)];

if E>Vm
km = sqrt(2*m*(E-Vm)/hbar^2);
%m4=[exp(i*km*b) exp(-i*km*b); i*km*exp(i*km*b) -i*km*exp(-i*km*b)];
%m5=[exp(i*km*c) exp(-i*km*c); i*km*exp(i*km*c) -i*km*exp(-i*km*c)];
%m45=m4*inv(m5);
m45 = [cos(km*L) (-1/km)*sin(km*L); km*sin(km*L) cos(km*L)];

elseif E<Vm
km = sqrt(2*m*(Vm-E)/hbar^2);
%m4=[exp(km*b) exp(-km*b); km*exp(km*b) -km*exp(-km*b)];
%m5=[exp(km*c) exp(-km*c); km*exp(km*c) -km*exp(-km*c)];
%m45=m4*inv(m5);
m45 = [cosh(km*L) (-1/km)*sinh(km*L); -km*sinh(km*L) cosh(km*L)];

else
m45=[1 -L; 0 1];
end

%m67=m6*inv(m7)
m67=[cosh(p*w2) (-1/p)*sinh(p*w2); -p*sinh(p*w2) cosh(p*w2)];

M=inv(m1)*m23*m45*m67*m8; %[A;B]=M*[J;K]

M1=[1 -M(3); 0 -M(4)];

M2=[0 M(1); -1 M(2)];
S=inv(M2)*M1; %Scattering matrix relates incoming to outgoing amplitudes
%[B;J]=S*[A;K]

T=abs(S(2))^2;

20
%resdens.m
%created 24/7/01
%last modified 26/7/01
%written by Shanna Crankshaw

%finds the electron density n(k) inside the well for one value of Vm

%m32=m3*inv(m2);
m32=[cosh(q*w1) (1/q)*sinh(q*w1); q*sinh(q*w1) cosh(q*w1)];

if E>Vm
km = sqrt(2*m*(E-Vm)/hbar^2);
m4=[exp(i*km*b) exp(-i*km*b); i*km*exp(i*km*b) -i*km*exp(-i*km*b)];
mu = inv(m4)*m32*m1;
Mu = mu*[1; S(1)];
F=Mu(1); G=Mu(2);
n=L*(abs(F))^2+L*(abs(G))^2+((F*conj(G))/(2*i*km))*(exp(2*i*km*c)-
exp(2*i*km*b))
-((conj(F)*G)/(2*i*km))*(exp(-2*i*km*c)-exp(-2*i*km*b));

elseif E<Vm
km = sqrt(2*m*(Vm-E)/hbar^2);
m4=[exp(km*b) exp(-km*b); km*exp(km*b) -km*exp(-km*b)];
mu = inv(m4)*m32*m1;
Mu = mu*[1; S(1)];
F=Mu(1); G=Mu(2);
(2*km))*(exp(-2*km*c)-exp(-2*km*b));(2*km))*(exp(2*km*c)-exp(2*km*b))-
(abs(G)^2/

else
m4=[b 1; 1 0];
mu = inv(m4)*m32*m1;
Mu = mu*[1; S(1)];
F=Mu(1); G=Mu(2);
n=(1/3)*(abs(F)^2)*(c^3-b^3)+(1/2)*(F*conj(G)+conj(F)*G)*(c^2-
b^2)+L*abs(G)^2;
end

21
%resonantcomp.m
%written by Shanna Crankshaw
%last modified 26/07/01
%creates density vectors for inside the well
%calculates current

format long;
parameters

energyTl=[]; %initializes energy vectors

energyTr=[];
Tl=[]; %initializes transmission vectors
Tr=[];
nl=[]; %initializes density vectors
nr=[];

for E = Emin:Eincr:muL
resonant;
resdens;
energyTl = [energyTl E];
Tl = [Tl T];
nl = [nl n];
end

V1l=V1; V2l=V2; w1l=w1; w2l=w2;

V2r=V1l; V1r=V2l; w2r=w1l; w1r=w2l;
V1=V1r; V2=V2r; w1=w1r; w2=w2r; %interchanges w1<->w2 and
V1<->V2;

clear V1l V2l w1l w2l V1r V2r w1r w2r

a=0; b=a+w1; c=b+L; d=c+w2; %defines parameters for -k states

for E=Emin:Eincr:muR
resonant;
resdens;
energyTr = [energyTr E];
Tr = [Tr T];
nr = [nr n];
end

energynl=energyTl;
energynr=energyTr;
lengthen; %makes energynl and energynr the same length
n=nl+nr;
parameters %resets parameters for calculation of next +k state

22
%lengthen.m
%last modified 18/07/01
%written by Shanna Crankshaw

%makes vectors energynl and energynr the same length

%pads the shorter density vector with zeros

if length(nl)>length(nr)
nright = 0*nl;
nright(1:length(nr)) = nr;
for count=length(nr)+1:length(nl)
nright(count) = 0;
end
nr = nright;
clear nright
energynr = energynl;
elseif length(nl)<length(nr)
nleft = 0*nr;
nleft(1:length(nl)) = nl;
for count=length(nl)+1:length(nr)
nleft(count) = 0;
end
nl = nleft;
clear nleft
energynl = energynr;
else
nl=nl;
nr=nr;
end

23
%nvsvm.m
%last modified 26/07/01
%written by Shanna Crankshaw

%calculates the total integrated density for many different Vm values

%calculates the total integrated current for many different Vm values

parameters %sets barrier parameters

muL=Emin+500*Eincr %defines electrochemical potentials
muR=Emin+400*Eincr

Vmvalues = 10*qe*(0:.0001:.33); %defines vector of Vm (well hgt) values

dens = zeros(length(Vmvalues,(max(muL,muR)-Emin)/incr+1);

for counter = 1:length(Vmvalues),

Vm=Vmvalues(counter);
counter
resonantcomp;
Nlvalues(counter) = trapz(energynl,nl); %integrates n for +k states
Nrvalues(counter) = trapz(energynr,nr); %integrates n for k states
Jlvalues(counter) = trapz(energyTl,Tl); %calc. current from +k states
Jrvalues(counter) = trapz(energyTr,Tr); %calc. current from +k states
dens(counter,: ) = n; %defines density matrix for single Vm
end

N=Nlvalues+Nrvalues; %total density

J=Jlvalues+Jrvalues; %net current

24
REFERENCES

[1] John S. Townsend, A Modern Approach to Quantum Mechanics. (University Science

Books, Sausalito, CA, 2000), p. 147-188.

[2] Paul Stevenson, Introductory Quantum Mechanics: A Second Course (Lecture notes for
Physics 312). (Rice University, Houston, TX, 2001), p. 1.2.

[3] Supriyo Datta, Electronic Transport in Mesoscopic Systems. (Cambridge University

Press, New York, 1995).

25