1. science
2. mathematics
3. algebra

Northern Arizona University, Vol. I, No. 1, 1-9
Fall 2014
Tunneling in Double Barriers
Congwei Wu12 *
Abstract
In this research project, the main focus was on the rectangular symmetric double barrier system.
There exists an analytic solution for this system hence it is possible to verify numerical results. We
applied multiple methods to compute the numerical values of transmission coefficients for different
initial energies. As a result, a comparison between analytic solution and numerical solutions is given.
Keywords
Double Barriers Tunneling, Monte Carlo, Transfer matrix, WKB approximation
1 Department
of Physics and Astronomy, Northern Arizona University, Flagstaff, AZ
of Mathematics and Statistics, Northern Arizona University, Flagstaff, AZ
*Corresponding author: Congwei Wu, [email protected]
2 Department
Contents
1
Introduction
1
Physics Model
1
1.1 Quantum Tunneling . . . . . . . . . . . . . . . . . . . 1
Challenge to classical physics • Quantum rectangular
potential barrier
1.2 Symmetric Double Barriers . . . . . . . . . . . . . 2
Derivation • Tunneling coefficient • Analytical solution
2
Methods
3
2.1 Method of Solving Non-Linear Equations . . 3
Separation • Metropolis Monte Carlo
2.2 Method of Transfer Matrix . . . . . . . . . . . . . . 4
Background • Transfer matrix • Transmission coefficient •
Probability density function f (x)
2.3 WKB Method with Numerov Scheme . . . . . 6
WKB method • Scheme of Numerov’s method • Algorithm
3
Results and Discussion
7
3.1 Transmission Coefficient . . . . . . . . . . . . . . . 7
devices, spontaneous DNA mutation, and Radioactive
decay. For the application of semiconductor devices,
tunnel diode is one example that uses the quantum tunneling.
For a more complicated structure such as double barriers, it is almost impossible to derive analytically solutions, and many numerical methods are being used. For
instance, the Wentzel-kramers-Brillouin (WKB) method
and Numerov methods are two approximations to solving the Schr¨odinger equation. Meanwhile, new methods such as transfer matrix for solutions of wave functions, and Metropolis Monte Carlo method for estimating
ground-state energy have been developed rapidly. Hence
there is a tendency for computational statistical methods to be involved in the research of quantum physics.
In addition, a circumstance of physicists who become
Quantum Bayesianisms also indicates evidence of combination, although there is still a divergence between
Frequentists and Bayesians in the area of statistics.
Simulation settings • Results and comparison
3.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 8
4
1. Physics Model
Conclusion
9
1.1 Quantum Tunneling
Appendix
9
1.1.1 Challenge to classical physics
References
9
Introduction
Quantum tunneling was developed in the 20th Century and it has widely areas of applications, such as
Scanning tunneling microscope (STM), semiconductor
In classical physics, an object can never be passing
through a physical barrier. It is also hard to imagine that
someone would pass through a wall with some probability. However, in the world of quantum mechanics, there
is some probability for particles to transmit through a
barrier, given that the kinetic energy of those particles
is lower than the required potential energy. So this type
Tunneling in Double Barriers — 2/9
of randomness gives us difficulties not only in physics
interpretation, but also in the methods involving partial
differential equations (PDE).
1.1.2 Quantum rectangular potential barrier
In quantum mechanics, the rectangular potential barrier is a one-dimensional problem which consists of
the tunneling phenomenon. This type of problem is required for solving the one-dimensional time-independent
Schr¨odinger equation. In addition, under a multiple barriers situation, an analytic solution to the Schr¨odinger
equation is not be able to obtained. In this study we are
focusing on the double-barrier problem.
1.2 Symmetric Double Barriers
In this research study, we are focusing on a symmetric structure of double barrier problem, and are proposing
statistical methods to quantum mechanics. As introduced
before, solutions to many problems require physicists to
solve the Sch¨odinger equation in order to obtain a closeform of the wave function. The wave function consists
of the fundamentals of the physical system.
d 2 ψ(x) 2m
+ 2 [E −V (x)]ψ(x) = 0
dx2
h¯
1.2.2 Tunneling coefficient
If E < V0 , then
ψI = A1 eik0 x + A2 e−ik0 x , −∞ < x < 0
ψII = B1 ek1 x + B2 e−k1 x , 0 < x < a
ψIII = C1 eik0 x +C2 e−ik0 x , a < x < b
k1 x
ψIV = D1 e
ψV = F1 e
To avoid tedious calculations, the derivations of transmission coefficient is given as follows. Since the xcoordinate does not have an impact on the derivations,
we constructed the double barriers question as: The bar-
(1)
where m is the mass of the particles, h¯ is the Planck constant, E is the total energy, V (x) is the potential energy,
and ψ(x) is the wave function. Due to the conservation
of energy, the potential V (x) is constant for all regions.
The Schr¨odinger equation can be written as five different
wave functions depending on different regions.
ik0 x
1.2.1 Derivation
−k1 x
+ D2 e
(2)
, b<x<c
, c<x<∞
Also, the wave numbers k are related to the energy via

q
 k = 2m(V − E)/¯h2 , x ∈ (0, a) ∪ (b, c)
1
0
q
k=
 k = 2mE/¯h2 ,
otherwise
0
V0
I
So the time-independent Schr¨odinger equation for
the wave function ψ(x) is
(3)
II
III
IV
a
c
b
Figure 1. General View of Double Barriers
0
V
x
riers (in yellow) are positioned between x = 0 and x = a,
and between x = b and x = c, respectively. The barriers
divide the one dimensional space in five parts. Denote
them as:
Table 1. Table of Regions
Name of Region
Region I
Region II
Region III
Region IV
Region V
Range of x
−∞ < x < 0
0<x<a
a<x<b
b<x<c
c<x<∞
The index on coefficients A, B,C, D and F denotes the
direction of the velocity vector, where 1 represents the
left and 2 represents the right. When E < V0 , k1 becomes imaginary and the wave function is exponentially
decaying within the barrier.
To determine the coefficients A, B,C, D and F, we
applied the boundary conditions of the wave function at
x = 0, x = a, x = b and x = c. Also note the derivative
of the wave function has to be continuous everywhere.
ψI (0) = ψII (0) ⇒ ψI0 (0) = ψII0 (0)
0
ψII (a) = ψIII (a) ⇒ ψII0 (a) = ψIII
(a)
0
0
ψIII (b) = ψIV (b) ⇒ ψIII
(b) = ψIV
(b)
(4)
0
ψIV (c) = ψV (c) ⇒ ψIV
(c) = ψV0 (c)
Substitute all ψ’s by the corresponding wave functions,
the boundary conditions yield the restriction on the coef-
Tunneling in Double Barriers — 3/9
2.1 Method of Solving Non-Linear Equations
ficients
2.1.1 Separation
A1 + A2 = B1 + B2
B1 e
k1 a
−k1 a
+ B2 e
ik0 a
= C1 e
+C2 e
A numerical solution can be given by solving the
system of non-linear equations (5) with Monte Carlo
simulation. The Euler’s formula yields:
−ik0 a
C1 eik0 b +C2 e−ik0 b = D1 ek1 b + D2 e−k1 b
eix = cos(x) + i sin(x)
D1 ek1 c + D2 e−k1 c = Feik0 c
A1 (ik0 ) − A2 (ik0 ) = B1 (k1 ) − B2 (k1 )
B1 k1 ek1 a − B2 k1 e−k1 a = C1 ik0 eik0 a −C2 ik0 e−ik0 a
C1 ik0 eik0 b −C2 ik0 e−ik0 b = D1 k1 ek1 b − D2 k1 e−k1 b
D1 k1 ek1 c − D2 k1 e−k1 c = Fik0 eik0 c .
(5)
1.2.3 Analytical solution
To simplify this type of problem, we treated the
system as consisting of two symmetric double barriers.
That is, widths of the two barriers are the same, namely
a − 0 = a = c − b. By Burstein and Lundqvist [1], the
closed-form of transmission coefficient can be written
as:
2
1
Ttrans = cosh(γl) + iα sinh(γl) · e−i2k0 L
2
−2
1 2
+ β sinh2 (γl)−2 ,
(6)
4
(8)
One important property of the Euler’s formula is that
it separates the imaginary parts from the real parts. So
we use the Euler’s formula rewrite (5) in terms of a
linear combination of a real part and an imaginary part.
That is, express (5) in terms of the form of (Creal ) +
(Cimaginary )i. For each part of this combination, we solve
for the coefficient independently for Creal and Cimaginary .
So (5) can be written as: (i) for all imaginary parts
A2 k0 = 0
(C1 −C2 )k0 cos(k0 a) = 0
(C1 +C2 ) sin(k0 a) = 0
(C1 −C2 )k0 cos(k0 b) = 0
(9)
(C1 +C2 ) sin(k0 b) = 0
F sin(k0 c) = 0
Fk0 cos(k0 c) = 0,
and (ii) for all real parts
where
A2 − (B1 + B2 ) = 0
l = a−0 = a = c−b
(B1 − B2 )k1 = 0
L = b−a
k0
γ
α= −
k0
γ
γ
k0
β= +
k0
γ
r
2m
γ=
(V0 − E).
h¯ 2
k1 (B1 ek1 a − B2 e−k1 a ) + (C1 +C2 )k0 sin(k0 a) = 0
(7)
Therefore, we obtained the closed-form solution of the
transmission coefficient, which can be applied for verifying our numerical results. Note that the analytic solution
is generally different for a different quantum tunneling
problem, but numerical methods are basically the same.
Hence this analytical solution for double barrier problem would help us to determine the best method among
various numerical methods.
2. Methods
B1 ek1 a + B2 e−k1 a − (C1 +C2 ) cos(k0 a) = 0
k1 (D1 ek1 b − D2 e−k1 b ) + (C1 +C2 )k0 sin(k0 a) = 0
D1 ek1 b + D2 e−k1 b − (C1 +C2 ) cos(k0 a) = 0
D1 ek1 c + D2 e−k2 c − F cos(k0 c) = 0
D1 k1 ek1 c − D2 k1 e−k1 c + Fk0 sin(k0 c) = 0.
(10)
So now our question is converted into a problem that
solves two sytsems of non-linear equations, using the
Monte Carlo simulation. In other words, solving a system of non-linear equations as follows

f1 (x1 , x2 , ..., xn ) = 0



 f2 (x1 , x2 , ..., xn ) = 0
(11)
..

.



fn (x1 , x2 , ..., xn ) = 0,
Tunneling in Double Barriers — 4/9
with a algorithm. We assume the true value of xi , i =
1, 2, ..., n is in some interval (ai , bi ). Each interval does
not necessarily be identical but for simplicity they would
be set as the same. By evaluating a large number of
points from this interval, we will eventually obtain a
sequence of points xi∗ , i = 1, 2, ..., n, which would satisfy
our equations within a given error ε. The error ε is
usually very small, and in this case we picked ε = 0.01.
2.1.2 Metropolis Monte Carlo
In statistics, the Metropolis-Hastings algorithm is a
Markov chain Monte Carlo (MCMC) method for obtaining a sequence of random samples from our interested
probability distribution [2], but it is difficult to sample directly. Statistically speaking, Markov chain is a
stochastic process in which every state of the system
only depends upon the previous state. Namely, to obtain
the information of n-th state, all we need is the information of the n − 1-th state of the system. Therefore, the
sequence obtained by MCMC can be used to estimate
some statistics that we are interested in, such as expected
value of some variable. For our purpose, we could also
apply MCMC to sample from the probability distribution (or equivalently square of the wave function) and
compute the ground state energy [3]. Since this study
focused mainly on tunneling, the ground state energy
would not be computed.
For each xi , there is a corresponding interval I =
0
(xi − α, xi0 + α) that covers the true value xi∗ . Then we
let Xi ∼ Unif(xi0 − α, xi0 + α) for the i-th variable. An
algorithm is given below: (for our case, n = 6 since we
set A1 = 1 and F2 = 0)
1. Define a detector function Ξ as
s
n
Ξ=
∑ [ fi (x1 , x2 , ..., xn )]2 .
(12)
4. If Ξnew < Ξold , then this trial is meaningful and
we let Ξold = Ξnew and let xi0 = xij . Because the
trail value xij is closer to the true value compared
with xi0 . Otherwise, go back to step 3 and select
another point xij+1 .
5. Since the trial is good, and we obtain a better Ξ,
move to the next variable xi+1 . That is, repeat step
∗ . If i = n, then set i = 1.
3 but for the point xi+1
Therefore, we applied this algorithm to either equation
(9) or (10), then we obtained a numerical values of the
coefficients A2 , B1 , B2 ,C1 ,C2 , D1 , D2 , and F, given that
E,V0 , a, b, c.
2.2 Method of Transfer Matrix
2.2.1 Background
The method of transfer matrix has been developed
rapidly over the last decades. As Ricco and Azbel mentioned [4], transfer matrix method can be applied to any
types of barrier case and extremely useful when it is hard
to obtain an analytical solution.
2.2.2 Transfer matrix
The Schr¨odinger equation for one dimensional time
independent particle with energy V0 is
d 2 ψ(x) 2m
+ 2 (E −V (x))ψ(x) = 0
dx2
h¯
(13)
where the m represents mass of the particle. We divided
any barrier into N number of intervals, and each interval
the solution of the Schr¨odinger equation can be written
as
ψn (x) = An eikn x + Bn e−ikn x
(14)
i=1
Note for the sequence of solutions {xi∗ }, we should
have Ξ = 0. So a value of Ξ detects whether a
sequence of {xi } is our solution.
2. Pick a starting point xi0 for each i, i = 1, 2, ..., n.
Let α be some relative big number, for example,
10000. Calculate Ξ, denote as Ξold . Define an
interval (xi0 − α, xi0 + α) for which we assume that
the true sequence of solutions {xi∗ } is in this ndimensional space.
3. For the j-th sampling, we randomly select a point
xij from (xi0 − α, xi0 + α), then plug it into (11). So
(11) gives a new value of Ξ. Denote as Ξnew .
where xn−1 < x < xn and kn =
Denote
q
2m(E −V (xn ))/¯h2 .
Knxn = eikn xn
Kn∗xn = e−ikn xn
xn
Kn+1
= eikn+1 xn
∗xn
Kn+1
= e−ikn+1 xn ,
(15)
For every two consecutive intervals, namely (xn , xn+1 )
and (xn+1 , xn+2 ), the left-hand side of wave functions
Tunneling in Double Barriers — 5/9
and their derivatives along with the wave functions is
ψn (xn ) = An Knxn + Bn Kn∗xn
ψn0 (xn ) = An ikn Knxn − Bn ikn Kn∗xn
∗x
x
n+1
n+1
ψn+1 (xn+1 ) = An+1 Kn+1
+ Bn+1 Kn+1
x
∗x
0
n+1
n+1
ψn+1
(xn+1 ) = An+1 ikn+1 Kn+1
− Bn+1 ikn+1 Kn+1
(16)
By the boundary conditions at x = xn+1 , we obtain
ψn (xn+1 )
ψn+1 (xn+1 )
ψn (xn+1 − δ )
=
=
0 (x
ψn0 (xn+1 )
ψn+1
ψn0 (xn+1 − δ )
n+1 )
(17)
where the width of interval is defined as δ = xn+1 − xn .
Denote
ik x
e n+1 n+1
eikn+1 xn+1
Λn (xn+1 ) =
ikn eikn xn+1 −ikn e−ikn xn+1
ik x
e n+1 n+1
eikn+1 xn+1
Λn+1 (xn+1 ) =
ikn+1 eikn+1 xn+1 −ikn+1 e−ikn+1 xn+1
(18)
Then
ψn (xn+1 )
A
= Λn (xn+1 ) n
ψn0 (xn+1 )
Bn
ψn+1 (xn+1 )
An+1
= Λn+1 (xn+1 )
0 (x
ψn+1
Bn+1
n+1 )
So we obtain
ψn (xn+1 )
ψn+1 (xn+1 )
= (Tn+1 ) 0
ψn0 (xn+1 )
ψn+1 (xn+1 )
(19)
So (20) and (22) together yield
ψ0 (x0 )
ψN (xN )
= T1 T2 ...TN 0
ψ00 (x0 )
ψN (xN )
where T j,k , j, k ∈ {1, 2} are elements of the transfer matrix T . So we rewrite (23) and thus
ψ0 (x0 )
ψN (xN )
Ta Tb ψN (xN )
=T
=
ψ00 (x0 )
ψN0 (xN )
Tc Td ψN0 (xN )
(25)
which gives
ψN (xN )
−1 ψ0 (x0 )
=T
ψN0 (xN )
ψ00 (x0 )
(26)
2.2.3 Transmission coefficient
By the definition of transfer matrix T , we can cite
find the transfer coefficient. Note that the wave functions
for region I and region V are defined as:
ψ0 (x0 ) = A1 eik0 x0 + A0 e−ik0 x0
(27)
ψN (xN ) = AN eikN xN
(28)
Then by (27), (28) and (25) we have
ik x
A1 e 0 0 + A0 e−ik0 x0
Ta Tb ψN (xN )
=
ψ00 (x0 )
Tc Td ψN0 (xN )
(29)
After tedious matrices calculation, it yields that
(20)
where Tn+1 is derived by Tn+1 = [Λn (xn+1 )][Λn (xn )−1 ],
namely
1
cos(kn+1 δ )
− kn+1
sin(kn+1 δ )
Tn+1 =
kn+1 sin(kn+1 δ )
cos(kn+1 δ )
(21)
So by (21) we define the matrix Tn as
cos(kn δ ) − k1n sin(kn δ )
Tn =
kn sin(kn δ )
cos(kn δ )
Then define the transfer matrix as
1
N sin(kn+1 δ )
cos(kn+1 δ )
− kn+1
T=∏
cos(kn+1 δ )
n=1 kn+1 sin(kn+1 δ )
(24)
AN =
2ik0 ei(k0 x0 −kN xN )
ik0 (Ta + ikN Tb ) + Tc + ikN Td
2
2ik0 ei(k0 x0 −kN xN )
Ttrans = |AN | = ik0 (Ta + ikN Tb ) + Tc + ikN Td 2
(31)
Rre f = 1 − Ttrans .
(22)
(23)
(30)
(32)
2.2.4 Probability density function f (x)
In order to derive the general probability density function for each n state, then we divide the total length of x
into N of lengths. By (20), let n + 1 = N then we have
ψN−1 (xN−1 )
ψN (xN )
= TN 0
,
(33)
0
ψN−2
(xN−2 )
ψN (xN )
Tunneling in Double Barriers — 6/9
Obviously, (23), (24) and (25) indicate that the matrix sequence of [ψn (xn )ψn0 (xn )]T is a geometric matrix
sequence, which can be written as:
ψn (xn )
AN eikN xN
n ψN (xN )
=
T
=
(34)
ψn0 (xn )
ψN0 (xN )
ikN AN eikN xN
where
1
cos(kn+1 δ )
− kn+1
sin(kn+1 δ )
T =∏
,
cos(kn+1 δ )
n=0 kn+1 sin(kn+1 δ )
n
Ta Tbn
n
.
(35)
T = n
Tc Tdn
N
n
So the wave function for each state of n is
ψn (xn ) = Tan An eikN xN + ikN Tbn AN eikN xN .
(36)
Then we obtain the probability density function f (xn ) for
each n state, for which the wave function is normalized
|ψn (xn )|2
.
f (xn ) = R ∞
2
−∞ |ψn (xn )| dxn
(37)
R
2.3 WKB Method with Numerov Scheme
2.3.1 WKB method
WKB approximation is a half-classical computational
method that can be used to solve a linear partial differential equations with varying coefficients [3]. The idea is
to use some exponential functions to replace the wave
functions in the Schr¨odinger’s equation, then find an
asymptotic series as an approximation of the exponential
function. Finally we take the limit of the asymptotic
series, and obtain an approximation of the wave function. Since we focus on the transmission coefficient in
this study, the WKB method was applied to find the
transmission coefficient instead of the wave functions.
In general, for a differential equation
dn y
dn−1 y
dy
+
a
x
+ ... + an−1 + an y = 0, (38)
1
n
n−1
dx
dx
dx
where ε > 0 is a very small number. We assume that
the solution can be expressed as an asymptotic series,
namely
y(x) ≈ exp
1 ∞ n
δ Sn (x) .
∑
δ n=0
φ10 = φ2
2m
φ20 = 2 (V (x) − E)φ1 .
h¯
(40)
(41)
We write the wave function as
∞
Since the presence of −∞
|ψn (xn )|2 dxn (constant), the
probability density function f (xn ) is normalized, which
has a range from 0 to 1.
ε
n
Then we substitute exp δ1 ∑∞
n=0 δ Sn (x) into the differential equation. Let δ → 0, then starting with n = 0
solve for every single Sn (x). Note most cases the series
is not convergent, so when n is sufficiently large, δ n Sn (x)
would increase. Therefore, once δ n Sn (x) increases, we
stop substituting n, and hence an approximate solution
yields.
For our symmetric double barrier, the differential
equation is the Schr¨odinger equation (1), where our variable of interest is the wave function ψ(x). To apply the
WKB method, first we write the Schr¨odinger equation
(1) as a system of two first-order linear ordinary differential equations with φ1 = ψ and φ2 = dψ/dx. That
is,
(39)
ψ = A(x)eiφ (x) .
(42)
Assume the potential is changing slowly, then d2 φ /dx2 →
0, the solution of the Schr¨odinger equation is
R
1
ψ = A p e±i/¯h |p|dx
|p|
p
k = 2m[E −V0 ].
(43)
(44)
So the asymptotic forms is a power series
∞
A=
∑ h¯ n An (x).
(45)
n=0
2.3.2 Scheme of Numerov’s method
Numerov’s method is a numerical method to solve ordinary differential equations [3], which is especially useful for the one-dimensional time-independent Schr¨odinger
equation. Combined together with the WKB method,
the error of our approximation could be reduced to
very small. There also exists less accurate methods for
solving the Schr¨odinger equation, such as Runge-Kutta
method. However, the basic algorithm for both methods is the same. As long as a minimization scheme of
the Numerov’s method is well designed, then it should
give a relatively close or useful numerical solution to
the double-barrier system. In our design scheme, the
Numerov’s method is mainly used to reduce the error of
the last term of the asymptotic series within the WKB
approximation.
Tunneling in Double Barriers — 7/9
From the Taylor expansion of g(x) about a point x0 :
N∈N
g(x) = g(x0 ) +
∑
k=1
(x − x0 )k g(k) (x0 )
k!
(46)
Then let the distance between x to x0 be h = x − x0 , and
then x = h + x0 . So we rewrite the Taylor expansion of
g(x) as:
So (52) yields the Numerov’s method, and the reminder
O(h6 ) is ignored in calculation. In other words, the error
in each time step is O(h6 ), which decreases with steps,
and the first term of this sequence of error is O(h4 ).
Therefore, we derived the Numerov formula for determining ψn , for n = 2, 3, ..., given ψ0 and ψ1 as two initial
values.
2.3.3 Algorithm
g(x0 + h) = g(x0 ) +
N∈N k (k)
h g (x
∑
0)
(47)
k!
k=1
So (47) means the point x0 of the Taylor expansion is
taking a step forward by an amount of h. That is, the
original expansion point of center is x0 , and we shift
this point to x0 + h. Likewise, we can also derive an
expansion of backward by taking a step of h back:
N∈N
g(x0 − h) = g(x0 ) +
∑
k=1
(−h)k g(k) (x0 )
k!
(48)
For computational purpose, the length h is small, so each
time the expansion point is shifted slightly, forward or
backward. Note
gn+1 = g(xn + h) = g(xn ) +
N∈N k (k)
h g (xn )
∑
k=1
N∈N
gn−1 = g(xn − h) = g(xn ) +
∑
k=1
k!
(49)
(−h)k g(k) (xn )
.
k!
(50)
The sum of (49) and (50) gives
00
gn−1 + gn+1 = 2gn + h2 gn +
h4 4
g + O(h6 ) (51)
12 n
where O(h6 ) denotes terms with order 6 or higher than
order 6. Namely, O(h6 ) denotes the reminder of the
Taylor expansion. Then we obtain a differential equation
as follows:
Applying the Numerov’s method to the double-barrier
problem requires an adjusted scheme, since each boundary condition contains complex numbers. Therefore, a
scheme designed for double-barrier is listed as follows:
1. In order to separate a complex parameter A, we let
2 2
A = Ar + iAi , given that a particle energy E = h¯2mk .
Then we guess two initial values ψ0 and ψ1 .
2. The transmission coefficient can be written as [5]
T3 (C4 −C3 ) C3 +C4 2
Ttrans = +
(54)
4
T3 1 T1
(55)
C3 = ( − )eiT2
T1
4
1
T1
(56)
C4 = ( + )e−iT2
4 T1
Z x4
|k|
T1 = exp(−
dx)
(57)
¯
x3 h
Z x3
|k|
dx
(58)
T2 =
¯
x2 h
Z x2
|k|
T3 = exp(−
dx).
(59)
¯
x1 h
Therefore, we obtained a numerical result of the transmission coefficient.
3. Results and Discussion
3.1 Transmission Coefficient
3.1.1 Simulation settings
Since our study project is focusing on the tunneling in
double barrier, it is necessary to obtain the transmission
h4
coefficient for our model. Note there are three numerical
h2 an gn = 2gn − gn−1 − gn+1 + g4n + O(h6 )
methods, so we set up various numerical systems such
12
(52) that every method can be evaluated under different circumstances. Additionally, the closed-form solution for
The solution of (52) is
the symmetric double barrier system is obtained, so our
numerical simulated environment would be also built on
2
h2
(2 − 5h6 an )gn − (1 + 12
an−1 )gn+1
gn+1 =
+ O(h6 ) this symmetric system. From (7) we knew that l reph2
1 + 12
an+1
resents the width of each barrier, and L represents the
(53) distance between two barriers.
Tunneling in Double Barriers — 8/9
Table 2. Table of Simulation Settings
Circumstances
V0 , in eV l, in m L, in m
I: Basic
1
10−10 10−10
II: Huge Barriers
1
10−9
10−10
III: Long Distance 1
10−10 20−9
IV: High Potential 2
10−10 10−10
reach the same level as in the basic model. For example,
at V0 = 1eV , it demands E = 0.20eV for a Ttrans ≈ 33%,
while at V0 = 2eV , it demands E = 0.40eV for the Ttrans
be able to pass 36%.
Second we applied the each numerical method to the
setting I, the basic model:
Table 4. Table of Ttrans from Different Methods
There are four settings that base on the symmetric
system, which were consisted of four distinct circumstances. The basic circumstance is a possibly realistic
model in real world, and the other three are extremal
situations in different ways. The huge barriers model has
two large barriers with widths of 10−9 meters and the distance between two barriers for the long distance model
is large, also reaching 2 × 10−9 meters. Note that the
total width of the systems for both models are the same,
so it is reasonable that the transmission coefficients for
both models be the same. But due to the existence of the
double barriers, the transmission coefficients might be
somehow different. The last model is called the high potential model, which has a relative high potential energy
for both barriers. Thus if given same kinetic energy, the
transmission coefficient for this model would be likely
lower than the basic model, because it requires higher
energy for particles to be tunneling.
3.1.2 Results and comparison
First we calculated the exact values for all four settings via (7). The results are as follows:
Table 3. Table of Exact Ttrans for All Settings
Results
E
0.10
0.20
0.30
0.40
0.50
I
0.01068
0.33145
0.41192
0.82572
0.73648
Settings
II
III
0.00032 0.01052
0.17954 0.33029
0.25698 0.40974
0.56232 0.82003
0.51041 0.72966
IV
0.00120
0.08061
0.20568
0.36557
0.48706
It is obvious that the long distance model has a higher
transmission coefficient than the huge barriers model. Although the total widths of both models are the same, the
width of barrier has the main effect on the transmission
coefficient. Compared the data between basic model and
long distance model, it is apparent that the Ttrans of long
distance model is very close to the Ttrans of the basic
model. So it indicates that the distance between two barriers does not have a big impact on the Ttrans . Finally, for
the high potential model, E must be higher to be able to
Results
E
0.10
0.20
0.30
0.40
0.50
Methods for Setting I
Monte Carlo Transfer Matrix
0.01564
0.01395
0.33149
0.35142
0.40190
0.42007
0.83555
0.83134
0.74647
0.75164
WKB
0.00002
0.13415
0.36287
0.52334
0.39787
As we expected, the Monte Carlo method is the most
accurate method among all three numerical methods.
For each energy level, the Monte Carlo method gives
almost the same result as the true value. In addition, the
transfer matrix method also yields relative close result,
compared to the true values. Note that as long as the
simulation time sufficiently large, both Monte Carlo and
transfer matrix methods could be very close to the true
values, within error less than 10−10 . In our study, we
constricted the computing time to about 120 minutes,
so both methods did not converge with a very small
error. Through this comparison, it is obvious that the
convergent speed of Monte Carlo is a little faster than
the transfer matrix’s. The reason is probably for the
optimization scheme for the Monte Carlo method, and
hence this method has a higher efficiency.
However, the WKB method seems to be the worst
among all numerical methods. Compared to the true
values, the WKB method is not only underestimated
the Ttrans , but also it has a relative big error, which is
over 40% of the estimate. Even though the Numerov’s
method came with the WKB to eliminate the error, the
results are still not acceptable. The reason for such poor
performance is expected to be the divergent asymptotic
series. As mentioned before, usually the asymptotic
series we chose to approximate is not convergent, and
thus it would give a relative large error because of the
divergence.
3.2 Discussion
Overall, the Monte Carlo method with the scheme
of optimization is the best method among all numerical
methods. A visualized figure is provided as follows:
Tunneling in Double Barriers — 9/9
Figure 2. Figure of data
Note that the estimated line of the Monte Carlo is
almost the same as the true value line, so is the line
of transfer matrix. Therefore, these two methods are
actually both good, even though the Monte Carlo method
is slightly better in estimation. In addition, it is obvious
that the WKB method is off by a large error.
References
[1]
E. Burstein and S. Lundqvist. Tunneling phenomena
in solids. Plenum Press, 1:9–12, 1969.
[2]
R.E. Salvino and F.A. Buot. Self-consistent monte
carlo particle transport with model quantum tunneling dynamics: Application to the intrinsic bistability
of a symmetric double-barrier structure. Journal of
Applied Physics, 72:5975–5981, 1992.
[3]
T. Pang. An introduction to computational physics.
Cambridge University Press, 2:110–115, 2006.
[4]
B. Ricco and M.Ya. Azbel. Physics of resonant unneling. the one-dimensional double-barrier case. Physical Review Letter B, 29:2, 1984.
[5]
A. Dutt and S. Kar. Smooth double barriers in
quantum mechanics. American Journal of Physics,
78:1352–1360, 2010.
4. Conclusion
In this study, for the symmetric rectangular double
barrier system, we derived the analytical solution, and
also introduced three numerical methods to obtain the
transmission coefficient. A simulation study is given
and the Monte Carlo and transfer matrix are the two best
numerical methods.
In the future study, the goal is to develop a better
scheme for the WKB method in order to obtain a less
error estimation. As a classical method, WKB still plays
an important role in the area of quantum mechanics so it
is necessary to have a better performance in estimation.
Appendix