Table Of ContentBand Tunneling through Double Barrier in Bilayer Graphene
Hasan A. Alshehaba,b, Hocine Bahloulia,b, Abderrahim El Mouhafidc and Ahmed Jellalb,c∗
aPhysics Department, King Fahd University of Petroleum & Minerals,
Dhahran 31261, Saudi Arabia
bSaudi Center for Theoretical Physics, Dhahran, Saudi Arabia
4
cTheoretical Physics Group, Faculty of Sciences, Choua¨ıb Doukkali University,
1
0
PO Box 20, 24000 El Jadida, Morocco
2
n
a
J
0
2
]
l
l
a
h Abstract
-
s
e By taking into account the full four band energy spectrum, we calculate the transmission prob-
m
ability and conductance of electrons across symmetric and asymmetric double potential barrier
t. with a confined interlayer potential difference in bilayer graphene. For energies less than the inter-
a
m layer coupling γ1, E < γ1, we have one channel for transmission which exhibits resonances, even
for incident particles with energies less than the strength of barriers, E < U , depending on the
- j
d
double barrier geometry. In contrast, for higher energies E >γ , we obtain four possible ways for
n 1
o transmission resulting from the two propagating modes. We compute the associated transmission
c probabilitiesaswellastheircontributiontotheconductance, studytheeffectofthedoublebarrier
[
geometry.
1
v
7
2
4
5
.
1
0
4
PACS numbers: 73.22.Pr; 73.63.-b; 72.80.Vp
1
: Keywords: bilayer graphene, double barrier, transmission, conductance.
v
i
X
r
a
∗[email protected], [email protected]
1 Introduction
In the last few years, graphene [1], a two dimensional one atom thick sheet of carbon, became a hot
research topic in the field of condensed matter physics. Its exceptional, electronic, optical, thermal,
and mechanical properties have potential future applications. For example, its thermal conductivity
is 15 times larger than that of copper and its electron mobility is 20 times larger than that of GaAs.
In addition, it is considered as one of the strongest materials with a Young’s modulus of about 1 TPa,
and some 200 times stronger than structural steel [2]. The most important application of graphene is
to possibly replace silicon in IT-technology. But the biggest obstacle is to create a gap and control
the electron mobility in graphene taking into account the so called Klein tunneling, which makes
the task more complicated [3,4]. However, one can create an energy gap in the spectrum in many
different ways, such as by coupling to substrate or doping with impurities [5,6] or in bilayer graphene
by applying an external electric field [7,8].
Bilayer graphene is two stacked (Bernal stacking [9]) monolayer graphene sheets, each with honey-
comb crystal structure, with four atoms in the unit cell, two in each layer. In the first Brillouin zone,
the tight binding model for bilayer graphene [10] predicted four bands, two conduction bands and
two valance bands, each pair is separated by an interlayer coupling energy of order γ ≈ 0.4 eV [11].
1
At the Dirac points, one valance band and one conduction band touch at zero energy, whereas the
other bands are split away from the zero energy by γ [12]. Further details about band structure and
1
electronic properties of bilayer graphene can be found in the literature [13–21]. Tunneling of quasi-
particles in graphene, which mimics relativistic quantum particles such as Dirac fermions in quantum
electrodynamics (QED), plays a major role in scattering theory. It allows to develop a theoretical
framework, which leads to investigate different physical phenomena that are not present in the non
relativistic regime, such as the Klein-paradox [3,4].
In monolayer graphene, there were many studies on the tunneling of electrons through different
potentialshapes[22–25]. Whilethestudyoftunnelingelectronsinbilayergraphenehasbeenrestricted
to energies less than the interlayer coupling parameter γ so that only one channel dominates trans-
1
mission and the two band model is valid [26–29]. Recently, tunneling of electrons in bilayer graphene
has been studied using the four band model for a wide range of energies, even for energies larger than
γ [30]. New transmission resonances were found that appear as sharp peaks in the conductance,
1
which are absent in the two band approximation.
Motivated by different developments and in particular [30], we investigate the band tunneling
through square double barrier in bilayer graphene. More precisely, the transmission probabilities and
conductanceofelectronswillbestudiedbytackingintoaccountthefullfourbandenergyspectrum. We
analyze two interesting cases by making comparison between the incident energies E and interlayer
coupling parameter γ . Indeed, for E < γ there is only one channel of transmission exhibiting
1 1
resonances,evenforincidentparticleswithE lessthanthestrengthofbarriersU (E < U ),depending
j j
on the double barrier geometry. For E > γ , we end up with two propagating modes resulting from
1
four possible ways of transmission. Subsequently, we use the transfer matrix and density of current
to determine the transmission probabilities and then the corresponding conductance. Based on the
physical parameters of our system, we present different numerical results and make comparison with
significant published works on the subject.
1
The present paper is organized as follows. In section 2, we establish a theoretical framework using
the four band model leading to four coupled differential equations. In section 3, by using the transfer
matrix at boundaries together with the incident, transmitted and reflected currents we end up with
eight transmission and reflection probabilities as well as the corresponding conductance. We deal with
two band tunneling and analyze their features with and without the interlayer potential difference, in
section 4. We do the same job in section 5 but by considering four band energy and underline the
difference with respect to other case. In section 6, we show the numerical results for the conductance
and investigate the contribution of each transmission channel. Finally, we briefly summarize our main
findings in the last section.
2 Theoretical model
Inmonolayergraphene,theunitcellhasinequivalentatoms(usuallycalledAandB).Bilayergraphene
on the other hand is a two stacked monolayer graphene (Bernal stacking) and hence has four atoms
in the unit cell. The relevant Hamiltonian near the K point (the boundary of the Brillouin zone), can
be found using the nearest-neighbor tight binding approximation [17]
V+ v π† −v π† v π
F 4 3
vFπ V+ γ1 −v4π†
H = (1)
−v π γ V− v π†
4 1 F
v π† −v π v π V−
3 4 F
where v = γ03a ≈ 106 m/s is the fermi velocity of electrons in each graphene layer, a = 0.142 nm is
F (cid:126) 2
the distance between adjacent carbon atoms, v = vFγ3,4 represent the coupling between the layers,
3,4 γ0
π = p +ip ,π† = p −ip arethein-planmomentaanditsconjugatewithp = −i(cid:126)∂ . γ ≈ 0.4eV
x y x y x,y x,y 1
istheinterlayercouplingtermandV+,V− arethepotentialsonthefirstandsecondlayer,respectively.
The skew parameters , γ ≈ 0.315 eV and γ ≈ 0.044 eV have negligible effect on the band structure
3 4
at high energy [12,13]. Recently, it was shown that even at low energy these parameters have also
negligible effect on the transmission [30], hence we neglect them in our calculations.
Under the above approximation and for double barrier potential configuration in Figure 1 our
U
4
U
2
1 2 3 4 5
a b c d x
Figure 1: The parameters of a rectangular double barrier structure.
previous Hamiltonian (1) can be written as follows in each potential region where we define regions
as follows: j = 1 for x ≤ a, j = 2 for a < x ≤ b, j = 3 for b < x ≤ c, j = 4 for c < x ≤ d and j = 5 for
2
x > d so that in the j-th region we have
V+ ν π† 0 0
j F
Hj = νF0π Vγj+ Vγ1− ν 0π† (2)
1 j F
0 0 ν π V−
F j
We define the potential on the first and second layer by V± = U ±δ , where U is the barrier strength
j j j j
and δ is the electrostatic potential in the j-th region
j
0, j = 1
U2±δ2, j = 2
V± = 0, j = 3 (3)
j
U4±δ4, j = 4
0, j = 5
(U ,δ ) and (U ,δ ) are the barrier potential and the electrostatic potential in regions 2 and 4, re-
2 2 4 4
spectively.
The eigenstates of (2) are four-components spinors ψj(x,y) = [ψj ,ψj ,ψj ,ψj ]†, here † denotes
A1 B1 A2 B2
the transpose of the row vector. For a double barrier we need to obtain the solution in each regions
as shown in Figure 1. Since we have basically two different sectors with zero (1, 3, 5) and nonzero
potential (2, 4), a general solution can be obtained in the second sector and then set the potential
V± to zero to obtain the solution in the first sector. To simplify the notation, let us introduce the
j
length scale l = (cid:126)vF ≈ 1.76 nm as well as E −→ E and V −→ Vj. Since the momentum along the
γ1 j γ1 j γ1
y-direction is a conserved quantity, i.e [H,p ] = 0, and therefore we can write the spinors as
y
ψj(x,y) = eikyy[φj ,φj ,φj ,φj ]T (4)
A1 B1 A2 B2
As usual, to derive the eigenvalues and the eingespinors we solve H ψ = E ψ . Then, by replacing
j j j j
by (2) and (4) we obtain
U +δ lπ† 0 0 φj φj
j j (cid:126) A1 A1
(cid:126)l0π Uj +1 δj U −1 δ l0π† φφjBj1 eiky = Ej φφjBj1 eiky (5)
j j (cid:126) A2 A2
0 0 lπ U −δ φj φj
(cid:126) j j B2 B2
This gives four coupled differential equations
(cid:20) (cid:21)
d
−il +ky φj = ((cid:15) −δ )φj (6)
dx B1 j j A1
(cid:20) (cid:21)
d
−il −ky φj +φj = ((cid:15) −δ )φj (7)
dx A1 A2 j j B1
(cid:20) (cid:21)
d
−il +ky φj +φj = ((cid:15) +δ )φj (8)
dx B2 B1 j j A2
(cid:20) (cid:21)
d
−il −ky φj = ((cid:15) +δ )φj (9)
dx A2 j j B2
3
where k is the wave vector along the y-direction and we have set (cid:15) = E −U . It is easy to decouple
y j j j
the first equations to obtain
(cid:20) d2 (cid:21)
+(ks)2 φj = 0 (10)
dx2 j B1
where the wave vector along the x-direction is
(cid:34) (cid:35)1/2
ks = −k2+ (cid:15)2j +(δj)2 +s1(cid:113)(cid:15)2(1+4(δ )2)−(δ )2 (11)
j y l2 l2 j j j
where s = ± denotes the propagating modes which will be discussed latter on. Now for each region
one can end up with corresponding wave vector according to Figure 1. Indeed, for regions 1, 3 and 5
we have V± = 0 and then we can obtain
j
(cid:20) (cid:15)2 (cid:15) (cid:21)1/2
ks = −k2+ +s (12)
0 y l2 l2
with (cid:15) = (cid:15) = (cid:15) = (cid:15) , as well as the energy
1 3 5
1 (cid:20) (cid:113) (cid:21)
Es = ± −s+ 1+(2lks)2+(2lk )2 (13)
± 2 0 y
However generally, for any region we can deduce energy from previous analysis as
(cid:34) (cid:114) (cid:35)1/2
(cid:15)s = ±√1 1+2l2(cid:2)(ks)2+k2(cid:3)+2δ2−s 1+4l2(cid:104)(ks)2+k2(cid:105)(cid:16)1+4δ2(cid:17) . (14)
±,j j y j j y j
2
The corresponding energy spectrum of the different regions is shown in Figure 2.
2 3.5 3.5
a b c
3.0 3.0
1
2.5 2.5
2.0 2.0
ΓE1 0(cid:144) ΓE1 1.5(cid:144) ΓE1 1.5(cid:144)
1.0 1.0
(cid:45)1
0.5 0.5
(cid:45)2 0.0 0.0
(cid:45)1.0 (cid:45)0.5 0.0 0.5 1.0 (cid:45)1.0 (cid:45)0.5 0.0 0.5 1.0 (cid:45)1.0 (cid:45)0.5 0.0 0.5 1.0
ky nm(cid:45)1 ky nm(cid:45)1 ky nm(cid:45)1
(cid:72) (cid:76) (cid:72) (cid:76) (cid:72) (cid:76)
Figure 2: Energy as function of the momentum k . (a): V = V = V = 0. (b): δ = 0.2γ ,
y 1 3 5 2 1
U = (0.5,1.2,1.8)γ (red, dashed green, dotted blue). (c): δ = 0.4γ , U = (0.8,1.4,2)γ (red,
2 1 4 1 4 1
dashed green, dotted blue). The dashed horizontal lines in (b) and (c) represent the heights of the
barriers U and U , respectively.
2 4
Associated with each real ks, the wave vector of the propagating wave in the first region, there
0
are two right-going (incident) propagating mode and two left-going (reflected) propagating mode. For
γ > E > 0, k+ is real while k− is imaginary, and therefore the propagation is only possible using k+
1 0 0 0
mode. However when E > γ , both k± are real and then the propagation is possible using two modes
1 0
4
k+ and k−. In Figure 3 we show these different modes and the associated transmission probabilities
0 0
through double barrier structure.
Figure 3: Schematic representation of different modes as well as the corrresponding transmission and
reflection probabilities.
Figure4presentstwodifferentcases. (a): asymmetricdoublebarrierstructureforU < U ,δ = δ
2 4 2 4
and (b): another symmetric one for U = U , δ = δ . It is interesting to note that the Ben results [30]
2 4 2 4
can be recovered from our results by considering the case (b) and requiring b = c in our double barrier
structure. We notice that different channels of transmission and reflection in Figure 3 can be mapped
into all cases in Figure 4 since they are related to the band structure on the both sides of the barriers.
However, the effect of the different structure of the two barriers should appear in the transmission and
reflection probabilities.
a b
V,E V,E
x,ky x,ky
Figure 4: Scheme represents the bandsinside and outside thebarriers for the same interlayer potential
difference. (a): asymmetric for U < U . (b): symmetric for U = U .
2 4 2 4
The solution of (10) can be written as a linear combination of plane waves
φjB1 = a1eikj+x+a2e−ikj+x+a3eikj−x+a4e−ikj−x (15)
where a (m = 1,2,3,4) are coefficients of normalization. Substituting (15) into (6 -9) we obtain the
m
rest of the spinor components:
φjA1 = a1A+−eikj+x−a2A++e−ikj+x+a3A−−eikj−x−a4A−+e−ikj−x (16)
φjA2 = a1ρ+eikj+x+a2ρ+e−jkj+x+a3ρ−eikj−x+a4ρ−e−ikj−x (17)
φjB2 = a1ζ++eikj+x−a2ζ−+e−ikj+x+a3ζ+−eikj−x−a4ζ−−e−ikj−x (18)
5
where we have set
(cid:104) (cid:105)
As± = l(k(cid:15)js±−iδky), ρs = ((cid:15)j −δj)1− l2 ((εkjs−)2δ+)2ky2 , ζ±s = (cid:15)(cid:15)j +−δδjρsAs±. (19)
j j j j j j
Now, we can write the general solution
ψj(x,y) = G M (x)C eikyy (20)
j j j
in terms of the matrices
A+ −A+ A− −A− eikj+x 0 0 0 a
− + − + 1
1 1 1 1 0 e−ikj+x 0 0 a2
Gj = ρ+ ρ+ ρ− ρ− ,Mj = 0 0 eikj−x 0 ,Cj = a3 (21)
ζ++ −ζ−+ ζ+− −ζ−− 0 0 0 e−ikj−x a4
Since we are using the transfer matrix, we are interested in the normalization coefficients, the compo-
nents of C, on the both sides of the double barrier. In other words, we need to specify our spinor in
region 1
φ1A1 = δs,1A+−eik0+x−r+sA++e−ik0+x+δs,−1A−−eik0−x−r−sA−+e−ik0−x (22)
φ1B1 = δs,1eik0+x+r+se−ik0+x+δs,−1eik0−x+r−se−ik0−x (23)
φ1A2 = δs,1ρ+eik1+x+r+sρ+e−ik0+x+δs,−1ρ−eik0−x+r−sρ−e−ik0−x (24)
φ1B2 = δs,1ζ++eik0+x−r+sζ−+e−ik0+x+δs,−1ζ+−eik0−x−r−sζ−−e−ik0−x (25)
as well as region 5
φ5 = ts A+eik0+x+ts A−eik0−x (26)
A1 + − − −
φ5 = ts eik0+x+ts eik0−x (27)
B1 + −
φ5 = ts ρ+eik0+x+ts ρ−eik0−x (28)
A2 + −
φ5 = ts ζ+eik0+x+ts ζ−eik0−x (29)
B2 + + − +
Since the potential is zero in regions 1, 3 and 5, we have the relation
G M (x) = G M (x) = G M (x). (30)
1 1 3 3 5 5
We will see how the above results will be used to determine different physical quantities. Specifically
we focus on the reflection and transmission probabilities as well as related matters.
3 Transmission probabilities and conductance
Implementing the appropriate boundary condition in the context of the transfer matrix approach, one
can obtain the transmission and reflection probabilities. Continuity of the spinors at the boundaries
gives the components of the vector C which are given by
δ ts
s,1 +
Cs = r+s , Cs = 0 (31)
1 δ 5 ts
s,−1 −
rs 0
−
6
where δ is the Kronecker delta symbol. The coefficients in the incident and reflected regions can
s,±1
be linked through the transfer matrix M
Cs = MCs (32)
1 5
which can be obtained explicitly by applying the continuity at the four boundaries of the double
barrier structure (Figure 1). These are given by
G M (a)C = G M (a)C (33)
1 1 1 2 2 2
G M (b)C = G M (b)C (34)
2 2 2 3 3 3
G M (c)C = G M (c)C (35)
3 3 3 4 4 4
G M (d)C = G M (d)C (36)
4 4 4 5 5 5
Now solving the above system of equations and taking into account of the relation (30), one can find
the form of M.
Then we can specify the complex coefficients of the transmission ts and reflection rs using the
± ±
transfer matrix M. Since we need the transmission T and reflection R probabilities and because the
velocity of the waves scattered through the two different modes is not the same, it is convenient to
use the current density J to obtain the transmission and reflection probabilities.
J = ν Ψ†α(cid:126)Ψ (37)
F
to end up with
|J | |J |
tra ref
T = , R = (38)
|J | |J |
inc inc
where α(cid:126) is a 4×4 diagonal matrix, on the diagonal 2 Pauli matrices σ . From (31) and (38), we show
x
that the eight transmission and reflection probabilities are given by [32]
k± k±
Ts = 0 |ts |2, Rs = 0 |rs|2 (39)
± ks ± ± ks ±
0 0
These expressions can be explained as follows. Since we have four band, the electrons can be scattered
between them and then we need to take into account the change in their velocities. With that, we
find four channels in transmission and reflection such that k± is also given by (12). More precisely, at
0
low energies (E < γ ), we have just one mode of propagation k+ leading to one transmission T and
1 0
reflection R channel through the two conduction bands touching at zero energy on the both sides of
the double barrier. Whereas at higher energy (E > γ ), we have two modes of propagation k+ and
1 0
k− leading to four transmission T± and reflection R± channels, through the four conduction bands.
0 ± ±
Since we have found transmission probabilities, let see how these will effect the conductance of our
system. This actually can be obtained through the Landauer-Bu¨ttiker formula [31] by summing on
all channels to end up with
G(E) = G Ly (cid:90) +∞dk (cid:88) Ts(E,k ) (40)
02π y n y
−∞ s,n=±
where L is the width of the sample in the y-direction and G = 4 e2, the factor 4 is due to the valley
y 0 h
and spin degeneracy in graphene.
7
1.0 1.0
0.8 0.8
0.6 0.6
T T
0.4 0.4
0.2 0.2
a b
0.0 0.0
0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8
E Γ E Γ
1 1
1.0 1.0
(cid:144) (cid:144)
0.8 0.8
0.6 0.6
T T
0.4 0.4
0.2 0.2
c d
0.0 0.0
0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0
E Γ E Γ
1 1
Figure 5: Transmission fo(cid:144)r normal incidence with b = b = 10 nm, and(cid:144) ∆ = 0 (blue dashed),
1 2
∆ = 5 nm (green), ∆ = 10 nm (red). (a): U = U = 0.4 γ . (b): U = 0.4 γ and U = 0.6 γ . (c,d):
2 4 1 2 1 4 1
the same parameters as in (a,b), respectively, but with δ = δ = 0.2 γ .
2 4 1
The obtained results will be numerically analyzed to discuss the basic features of our system and
also make link with other published results. Because of the nature of our system, we do our task by
distinguishing two different cases in terms of the band tunneling.
4 Two band tunneling
Barbier[27]investigatedthetransmissionandconductanceforsingleandmultipleelectrostaticbarriers
with and without interlayer potential difference and for E < γ , however the geometry dependance
1
of the transmission was not done. In this section, we briefly investigate the resonances resulting from
the available states in the well between the two barriers and how they influence by the geometry of
the system.
For a normal incidence and for δ = δ = 0 the transmission amplitude is shown in Figure 5a for
2 4
different values of the distance ∆ between the barriers. The dashed blue curve is for a single barrier
with (∆ = 0) and with width (b +b = 20 nm), we note that the transmission is zero and there are no
1 2
resonances in this regime of energy (E < U = U ). Unlike the case of the single barrier, the double
2 4
barrierstructurehasresonancesintheabovementionedrangeofenergy. Thesefulltransmissionpeaks
can be attributed to the bound electron states in the well region between the barriers. In agreement
with [33], the number of these resonances depends on the distance between the barriers. Indeed, for
8
Figure 6: Density plot of transmission probability, for δ = δ = 0, versus (a): E and ky for
2 4
U = U = 0.4 γ , and b = b = ∆ = 10 nm, (b): E and ∆ with ky = 0 and b = b = 10 nm, (c):
2 4 1 1 2 1 2
E and ky with U = 0.4 γ , U = 0.6 γ and b = b = ∆ = 10 nm, (d): E and ∆ with U = 0.4 γ ,
2 1 4 1 1 2 2 1
U = 0.6 γ , ky = 0 and b = b = 10 nm. White and black dashed lines represent the band inside
4 1 1 2
and outside the first barrier, respectively.
∆ = 5 nm we have one peak in the transmission amplitude, increasing the distance allows more bound
states to emerge in the well, and for ∆ = 10 nm there are two peaks (green and red curves in Figure
5, respectively). Figure 5b shows the same results in 5a but with different height of the two barriers
such that U = 0.4 γ and U = 0.6 γ . We see that the asymmetric structure of the double barrier
2 1 4 1
reduces those resonances resulting from the bound electrons in the well between the two barriers. For
δ = δ = 0.2 γ , we show the transmission probability by choosing U = U = 0.4 γ in Figure 5c and
2 4 1 2 4 1
for U = 0.4 γ , U = 0.6 γ in Figure 5d. For single barrier, there are no resonant peaks inside the
2 1 4 1
induced gap which is not the case for the double barrier as clarified in Figure 5c.
Figures 6a,6c present a comparison of the density plot of the transmission probability as a function
of the transverse wave vector k of the incident wave and its energy E between different structure
y
of the double barrier with U = U = 0.4 γ and U = 0.4 γ < U = 0.6 γ , respectively, and for
2 4 1 2 1 4 1
δ = δ = 0 in both. For non-normal incidence in Figure 6a (k (cid:54)= 0) we still have a full transmission,
2 4 y
even for energies less than the height of the barriers, which are symmetric in k . Those resonances are
y
reduced and even disappeared in Figure 6c due to the asymmetric structurer of the double barrier. In
Figures 6b,6d we show the density plot of transmission probability, for normal incidence, as a function
of ∆ and E for the same parameters as in Figure 6a and 6c, respectively. We note that the number
of resonances in Figure 6b, due to the bounded electrons in the well between the barriers, increases
as long as the distance is increasing. They are very sharp for the low energies and become wider at
higher energies. In contrast to Figure 6d and as a result of the asymmetric structure of the double
9
