Electron transport, interaction and spin in graphene and graphene nanoribbons

Linköping Studies in Science and Technology Dissertation No. 1469

Electron transport, interaction and spin in graphene and graphene nanoribbons

Artsem Shylau

Department of Science and Technology Linköping University, SE-601 74 Norrköping, Sweden

Norrköping 2012

Electron transport, interaction and spin in graphene and graphene nanoribbons

© 2012 Artsem Shylau

Department of Science and Technology Linköping University, SE-601 74 Norrköping, Sweden

ISBN 978-91-7519-816-3 ISSN 0345-7524

Printed by LiU-Tryck, Linköping 2012

To my parents, Valentina and Alexander.

Abstract

Since the isolation of graphene in 2004, this novel material has become the major object of modern condensed matter physics. Despite of enor- mous research activity in this field, there are still a number of fundamental phenomena that remain unexplained and challenge researchers for further investigations. Moreover, due to its unique electronic properties, graphene is considered as a promising candidate for future nanoelectronics. Besides experimental and technological issues, utilizing graphene as a fundamental block of electronic devices requires development of new theoretical meth- ods for going deep into understanding of current propagation in graphene constrictions.

This thesis is devoted to the investigation of the effects of electron- electron interactions, spin and different types of disorder on electronic and transport properties of graphene and graphene nanoribbons.

In paper I we develop an analytical theory for the gate electrostatics of graphene nanoribbons (GNRs). We calculate the classical and quan- tum capacitance of the GNRs and compare the results with the exact self- consistent numerical model which is based on the tight-binding p-orbital Hamiltonian within the Hartree approximation. It is shown that electron- electron interaction leads to significant modification of the band structure and accumulation of charges near the boundaries of the GNRs.

It’s well known that in two-dimensional (2D) bilayer graphene a band gap can be opened by applying a potential difference to its layers. Calcula- tions based on the one-electron model with the Dirac Hamiltonian predict a linear dependence of the energy gap on the potential difference. In paper II we calculate the energy gap in the gated bilayer graphene nanoribbons (bGNRs) taking into account the effect of electron-electron interaction. In contrast to the 2D bilayer systems the energy gap in the bGNRs depends non-linearly on the applied gate voltage. Moreover, at some intermediate gate voltages the energy gap can collapse which is explained by the strong modification of energy spectrum caused by the electron-electron interac- tions.

Paper III reports on conductance quantization in grapehene nanorib-

bons subjected to a perpendicular magnetic field. We adopt the recursive

Green’s function technique to calculate the transmission coefficient which

is then used to compute the conductance according to the Landauer ap- proach. We find that the conductance quantization is suppressed in the magnetic field. This unexpected behavior results from the interaction- induced modification of the band structure which leads to formation of the compressible strips in the middle of GNRs. We show the existence of the counter-propagating states at the same half of the GNRs. The over- lap between these states is significant and can lead to the enhancement of backscattering in realistic (i.e. disordered) GNRs.

Magnetotransport in GNRs in the presence of different types of disorder is studied in paper IV. In the regime of the lowest Landau level there are spin polarized states at the Fermi level which propagate in different directions at the same edge. We show that electron interaction leads to the pinning of the Fermi level to the lowest Landau level and subsequent formation of the compressible strips in the middle of the nanoribbon. The states which populate the compressible strips are not spatially localized in contrast to the edge states. They are manifested through the increase of the conductance in the case of the ideal GNRs. However due to their spatial extension these states are very sensitive to different types of disorder and do not significantly contribute to conductance of realistic samples with disorder. In contrast, the edges states are found to be very robust to the disorder. Our calculations show that the edge states can not be easily suppressed and survive even in the case of strong spin-flip scattering.

In paper V we study the effect of spatially correlated distribution of im- purities on conductivity in 2D graphene sheets. Both short- and long-range impurities are considered. The bulk conductivity is calculated making use of the time-dependent real-space Kubo-Greenwood formalism which allows us to deal with systems consisting of several millions of carbon atoms. Our findings show that correlations in impurities distribution do not signifi- cantly influence the conductivity in contrast to the predictions based on the Boltzman equation within the first Born approximation.

In paper VI we investigate spin-splitting in graphene in the presence of

charged impurities in the substrate and calculate the effective g-factor. We

perform self-consistent Thomas-Fermi calculations where the spin effects

are included within the Hubbard approximation and show that the effective

g-factor in graphene is enhanced in comparison to its one-electron (non-

interacting) value. Our findings are in agreement to the recent experimental

observations.

Popul¨ arvetenskaplig sammanfattning

Anda sedan isoleringen av grafen ˚ ¨ ar 2004, har detta nya material blivit den viktigaste f¨orem˚ alet f¨or den moderna kondenserade materiens fysik.

Trots enorm forskning inom detta omr˚ ade finns det fortfarande ett an- tal grundl¨aggande fenomen som f¨orblir of¨orklarade och utmanar forskare f¨or vidare unders¨okningar. Dessutom, p˚ a grund av dess unika elektron- iska egenskaper, anses grafen vara en lovande kandidat f¨or framtida na- noelektronik. F¨orutom experimentella och teknologiska fr˚ agor, kan grafen anv¨andas som ett grundl¨aggande block av elektroniska komponenter som kr¨aver utveckling av nya teoretiska metoder f¨or att f¨ordjupa f¨orst˚ aelsen av str¨om utbredning i nanostrukturer av grafen.

Denna avhandling till¨agnar ˚ at utredningen av effekterna av elektron-

elektron v¨axelverkan, spin och olika typer av oordning inom elektroniska

och transport egenskaper hos grafen och grafen nanoremsor.

Acknowledgments

Despite of my primary background in applied sciences I always wished to be a theoretical physicist, because I was always impressed by the fact that human mind is able to unravel Nature secrets just with the use of a piece of paper and a pencil (or a computer nowadays). This dream got fulfilled in Sweden where I spent four years as a PhD student at ITN LiU doing my research in theoretical physics. Tack, Sverige!

I would like to thank a lot of people who surrounded me during this period.

First of all, I would like to thank Prof. Igor Zozoulenko for his great supervision, significant contribution to my scientific development and for motivation when it was really needed.

I am thankful to the research administrator Ann-Christin Nor´en and Elisabeth Andersson for their administrative help during my research work.

I had a fruitful collaboration with my colleagues from Germany, Dr.

Hengyi Xu and Prof. Thomas Heinzel, and my polish colleague Dr. Jaros law K los.

Besides the science there were a lot of other things happening in my life.

I met a lot of nice people here, some of them eventually became my friends.

First of all, I am grateful to Olga Bubnova for our inspiring philosophical discussions, pleasant lunch time and just for being a good friend. I am also thankful to Lo¨ıg, Anton, Brice, Sergei, Julia and especially Taras for a good company and funny talks.

All this time I kept in touch with my belarusian friends, Sergei and Alex. We had a good time together whenever I went home, to my lovely Minsk.

Finally, I would like to thank my wife Marina for her love, understanding and patience; my parents and my sister’s family for their love and support which I feel every day no matter where I am.

Artsem Shylau

Norrk¨oping, August 2012

List of publications

Publications included in the thesis

1. A. A. Shylau, J. W. Klos, and I. V. Zozoulenko, Capacitance of graphene nanoribbons, Phys. Rev. B 80, 205402 (2009).

Author’s contribution: Implementation of the self-consistent nu- merical model and all numerical calculations, partial contribution to development of the analytical model. Preparation of all the figures.

Initial draft of the paper.

2. Hengyi Xu, T. Heinzel, A. A. Shylau, I. V. Zozoulenko, Interactions and screening in gated bilayer graphene nanoribbons, Phys. Rev. B 82, 115311 (2010); selected as ”Editor’s Suggestion”.

Author’s contribution: Development of the analytical model. Dis- cussion and analysis of all the obtained results.

3. A. A. Shylau, I. V. Zozoulenko, Hengyi Xu, T. Heinzel, Generic suppression of conductance quantization of interacting electrons in graphene nanoribbons in a perpendicular magnetic field, Phys. Rev.

B 82, 121410(R) (2010).

Author’s contribution: All numerical calculations and preparation of all the figures. Discussion of the results and writing an initial draft.

4. A. A. Shylau and I. V. Zozoulenko, Interacting electrons in graphene nanoribbons in the lowest Landau level, Phys. Rev. B 84, 075407 (2011).

Author’s contribution: All numerical calculations and preparation of all the figures, discussion of the results and writing an initial draft.

5. T. M. Radchenko, A. A. Shylau, and I. V. Zozoulenko, Influence of correlated impurities on conductivity of graphene sheets: Time- dependent real-space Kubo approach, Phys. Rev. B 86, 035418 (2012).

Author’s contribution: Contribution to implementation of the nu-

merical model, performing initial numerical calculations, discussion

of all results.

6. A. V. Volkov, A. A. Shylau, and I. V. Zozoulenko, Interaction- induced enhancement of g-factor in graphene, arXiv:1208.0522v1 [cond- mat.mes-hall], submitted to PRB.

Author’s contribution: Contribution to development of the model, discussion of all results, writing a part of the manuscript, supervision of the numerical calculations.

Relevant publications not included in the thesis

1. J. W. Klos, A. A. Shylau, I. V. Zozoulenko, Hengyi Xu, T. Heinzel, Transition from ballistic to diffusive behavior of graphene ribbons in the presence of warping and charged impurities, Phys. Rev. B 80, 245432 (2009).

2. T. Andrijauskas, A. A. Shylau, I. V. Zozoulenko, Thomas-Fermi

and Poisson modeling of the gate electrostatics in graphene nanorib-

bon, Lith. J. of Phys. 52, 63 (2012).

Contents

1 Introduction 1

2 Quantum transport 3

2.1 Kubo formalism . . . . 4

2.2 Kubo-Greenwood formula . . . . 6

2.3 Landauer approach . . . . 7

2.4 S-matrix technique . . . . 9

3 Electron-electron interactions 11 3.1 The many body problem . . . 11

3.2 Hartree-Fock approximation . . . 12

3.3 Density-functional theory . . . 14

3.4 Kohn-Sham equations . . . 15

3.5 Thomas-Fermi-Dirac approximation . . . 16

3.6 Hubbard model . . . 17

4 Electronic structure and transport in graphene 19 4.1 Basic electronic properties . . . 19

4.2 Graphene nanoribbons . . . 24

4.3 Warping . . . 27

4.4 Bilayer graphene . . . 28

4.5 Dirac fermions in a magnetic field . . . 31

5 Modeling 35 5.1 Tigh-binding Hamiltonian and Green’s function . . . 35

5.2 Recursive Green’s function technique . . . 38

5.2.1 Dyson equation . . . 38

5.2.2 Bloch states . . . 40

5.2.3 Calculation of the Bloch states velocity . . . 42

5.2.4 Surface Green’s function . . . 43

5.2.5 Transmission and reflection . . . 44

5.3 Real-space Kubo method . . . 45

5.3.1 Diffusion coefficient . . . 45

5.3.2 Transport regimes . . . 46

5.3.3 Time evolution . . . 48

5.3.4 Chebyshev method . . . 49

5.3.5 Continued fraction technique . . . 51

5.3.6 Tridiagonalization of the Hamiltonian matrix . . . 53

5.3.7 Local density of states . . . 54

6 Summary of the papers 55 6.1 Paper I . . . 55

6.2 Paper II . . . 56

6.3 Paper III . . . 57

6.4 Paper IV . . . 58

6.5 Paper V . . . 59

6.6 Paper VI . . . 59

Chapter 1 Introduction

Inside every pencil, there is a neutron star waiting to get out.

To release it, just draw a line.

(New Scientist, 2006) In 2004 the researchers from Manchester University, Kostya Novoselov, Andre Geim and collaborators, reported on experimental isolation of graphene [1], a pure 2D crystal consisting of carbon atoms arranged in a honey-comb lattice. For a long time before graphene had been considered only by the- oreticians as a basic block used to build theory for graphite [2] and carbon nanotubes [3]. Its existence was doubt since the theory predicted that perfect 2D crystals are not thermodynamically stable [4]. The discovery of graphene triggered a great scientific interest in this field, as a result graphene is one of the most extensively studied object in a modern con- densed matter physics [5].

Due to its specific lattice structure graphene posses a number of unique electronic properties which make this material interesting for both theoreti- cians, experimentalists and engineers. One of the most important feature of graphene is its linear energy spectrum. This kind of spectrum is known from high-energy physics where it corresponds to massless particles like neutrino. Relativistic-like dispersion relation is responsible for such effects as Klein tunneling - unimpeded penetration of particle through the in- finitely large potential barrier [6]. The experimental discovery of graphene had led to emergence of a new paradigm of ’relativistic’ condensed-matter physics and provided a way to probe quantum electrodynamics phenomena [4]. That is why graphene is sometimes called ”CERN on a desk”.

Being massless fermions electrons in graphene propagate at extremely

high velocities, only 300 times smaller than the velocity of light. This

makes graphene the best known conductor with mobility up to 200,000

cm ² V ⁻¹ s ⁻¹ at room temperatures. Graphene subjected to a perpendicular

magnetic field exhibits the anomalous quantum Hall effect [7]. This is a re-

sult of unusual spectrum quantization and the presence of the 0’th Landau

level, which is equally shared between electrons and holes. Also fractional quantum Hall effect has been experimentally observed in graphene [8].

The range of possible applications of graphene is very broad. With its

high mobility graphene is considered as the main candidate for a future

post-silicon electronics [9]. It is particularly interesting to use graphene in

transistors operating at ultrahigh radio frequencies [10]. Due to its high

optical transmittance (≈ 97.7%) graphene is proposed to be used as a flexi-

ble transparent electrode in touchscreen devices [11]. Also graphene posses

a broad spectral bandwidth and fast responce times, which makes this

material attractive for optoelectronics and, in particular, phototransistors

[12].

Chapter 2

Quantum transport

Electrical transport is a non-equilibrium statistical problem. In principle, one can solve the time-dependent Schr¨odinger equation

i~ ∂|Ψ(~r, t)i

∂t = ˆ H|Ψ(~r, t)i (2.1)

to find a many-body state of the system |Ψ(~r, t)i at any time and then calculate the expectation value of the current operator

I = ˆ Z

S ˆj(~r, t) · d~S. (2.2)

In practice, however, the calculation of the many-body state is an unfea- sible task. Moreover, the wave function representing the state provides a detailed information about the system which is often redundant for deter- mination of transport properties. Hence, one needs to introduce a number of approximations in order to simplify the above problem. Prior doing this it is useful to formulate viewpoints underlying quantum transport theories [13]:

Viewpoint 1: The electrical current is a consequence of an applied electric field: the field is the cause, the current is the response to this field.

Viewpoint 2: The electrical current is determined by the boundary conditions at the surface of the sample. Charge carriers incident to the sam- ple boundaries generate self-consistently an inhomogeneous electric field across the sample. Thus the field is a consequence of the current.

In this chapter two approaches for calculation of transport properties are discussed, namely, the Kubo formalism and the Landauer approach.

The first one belongs to the viewpoint 1, while the latter belongs to the

viewpoint 2.

2.1 Kubo formalism

Kubo formula relates via linear response the conductivity to the equilibrium properties of the system. First we define the expectation value of the current density operator

h~ji = Tr h

~jρ i

, (2.3)

where ρ is a statistical operator

ρ = Z ⁻¹ e ⁻

, Z = Tr h e ⁻

i

. (2.4)

The Hamiltonian, H, describing the system can be split into two parts H = H ₀ + δH,

ρ = ρ ₀ + δρ, (2.5)

where H ₀ corresponds to the system in a global cannonical equilibrium described by the statistical operator ρ 0 . External electric field introduces perturbation in the system which is described by the terms δH and δρ.

Here we focus on the case of uniform electric field [14]. It’s assumed that the field is applied at time t = −∞ and reaches adiabatically its steady value at t = 0

δH = lim

α→0

h e ~ E~re ^{(−iωt+αt)} i

. (2.6)

Substituting ρ from Eq.(2.5) into Eq.(2.3), we have h~ji = Tr h

~j(ρ 0 + δρ) i

= Tr h

~jδρ i

, (2.7)

where we took into account that Tr h

~jρ 0

i = 0, since there is no current in the system before applying the external field. If we substitute ρ given by Eq.(2.5) into the Liouville equation, i~ ˙ρ = [H, ρ], which describes the dynamics of a closed quantum systems, then the change of δρ in time can be found as

i~δ ˙ρ = [H ₀ , δρ] + [δH, ρ ₀ ], (2.8) where we used that i~ ˙ρ ₀ = [H ₀ , ρ ₀ ] and neglected the term [δH, δρ]. Let us now pass to the interaction picture representation of time-dependence of an operator

δρ = e ⁻

∆ρe

. (2.9) Substituting Eq.(2.9) into the left part of Eq.(2.8), we arrive at

i~∆ρ = e

[δH, ρ ₀ ]e ⁻

= lim _α→0 h

e ^{(−iωt+αt)} e

[e~r, ρ 0 ]e ⁻

E ~ i

. (2.10)

Kubo formalism 5

Both δρ and ∆ρ satisfy the following conditions: have the same value at time t = 0, δρ(0) = ∆ρ(0), and equal to zero at t = −∞, δρ(−∞) =

∆ρ(−∞) = 0. Thus, integrating Eq.(2.10), we get δρ(t = 0) = 1

i~ lim

α→0

Z 0

−∞

dte ^{(−iωt+αt)} e

[e~r, ρ ₀ ]e ⁻

E. ~ (2.11) Substituting Eq.(2.11) into the expression (2.7) for the expectation value of the current density, we obtain

h~ji = Tr h

~jδρ i

= 1 i~ lim

α→0

Z 0

−∞

dte ^{(−iωt+αt)} Tr h

~je

[e~r, ρ ₀ ]e ⁻

E ~ i . (2.12) The conductivity tensor is defined as

 j x

j _y



=

 σ xx σ xy

σ _yx σ _yy

  E x

E _y



, (2.13)

which allows to deduce from Eq.(2.12) the components of the tensor σ µν = lim

α→0

Z 0

−∞

dte ^{(−iωt+αt)} K µν , (2.14)

where

K µν = 1 i~ Tr h

j ~ µ e

[er ν , ρ 0 ]e ⁻

i

. (2.15)

Equation (2.15) can be written in another form if we use the following relation [14]:

[r _ν , ρ ₀ ] = ρ ₀ Z 1/k

T

0

dλe ^λH

[H ₀ , r _ν ]e ^−λH

. (2.16) Taking into account that [H ₀ , r _ν ] = −i~ ˙r ν and −e ˙r ν = j _ν , we get

[er ν , ρ 0 ] = i~ρ 0

Z 1/k

T 0

dλe ^λH

j ν e ^−λH

. (2.17) Finally, using Heisenberg representation for time-dependence of the current

~j(t) = e

~je ⁻

, (2.18) we obtain the Kubo formula for conductivity which is formulated in terms of the current-current response function

K _µν = Z 1/k

T

0 dλhj µ (0)j _ν (t − i~λ)i. (2.19)

2.2 Kubo-Greenwood formula

Equations (2.14) and (2.19), which constitute the Kubo formula for con- ductivity are very important, since they reflect the underlying physics of the linear responce theory. However, these equations are not suitable for a particle calculations and one needs to work out another form of the con- ductivity formula. As in the previous section we start with splitting the Hamiltonian, H, into two parts H = H ₀ + δH corresponding to the sys- tem in the equilibrium and the perturbation respectively. Hence, if the full Hamiltonian is defined as

H = 1

2m (~p + e ~ A) ² + eφ (2.20) and the perturbation is incorporated in the vector potential, ~ A = ~ A ₀ + ~ A _ext , we get for δH:

δH = e ~ A _ext · ~v. (2.21)

If we assume the time dependence of the external electric field to be ~ E(t) = Ee ~ ^−iωt , then using the relation ~ E = − ^{d ~} _dt ^A , we arrive at

δH = e

iω E · ~v. ~ (2.22)

The expectation value of the current operator, hji = Tr h ρˆj i

, can be written as

hji = X

kk

hk|δρ|k ^′ ihk ^′ |ˆj|ki = V ² Z Z

dEdE ^′ g(E)g(E ^′ )hE|δρ|E ^′ ihE ^′ |ˆj|Ei, (2.23) where g(E) is a density of states and we used Tr[ρ ₀ ˆj] = 0. If we assume δρ(t) = δρ · e ^iωt and use Eq.(2.8), then we get

hE ^′ |δρ|Ei = f _{F D} (E ^′ ) − f F D (E)

E ^′ − E − ~ω − i~α hE ^′ |δH|Ei, (2.24) where α is a small constant and we used that ρ 0 |Ei = f F D (E)|Ei [14]. Sub- stituting Eq.(2.24) into Eq.(2.23) and recalling that ~j = − _V ^e ~v, we obtain the expectation value of the current operator

hji = V ² Z Z

dEdE ^′ g(E)g(E ^′ )  e iω

  − e V



× f _{F D} (E ^′ ) − f F D (E)

E ^′ − E − ~ω − i~α hE ^′ | ~ E · ~v|EihE|~v|E ^′ i. (2.25) This equation allows us to derive the components of the conductivity tensor

σ _ij (ω) = Re



− e ² iω V

Z Z

dEdE ^′ g(E)g(E ^′ )

× f F D (E ^′ ) − f F D (E)

E ^′ − E − ~ω − i~α hE ^′ |~v i |EihE|~v j |E ^′ i



, (2.26)

Landauer approach 7

which can be further simplified using the relation Re



α→0 lim 1 i

1

E ^′ − E − ~ω − i~α



= πδ(E ^′ − E − ~ω) (2.27) and performing integration over E ^′

σ _ij (ω) = − e ² πV ω

Z

g(E)g(E + ~ω)hE + ~ω|~v i |EihE|~v j |E + ~ωi

× [f F D (E + ~ω) − f F D (E)]dE. (2.28) If one is interested in DC conductivity, Eq.(2.28) should be considered in a limit ω → 0, lim ω→0 f

(E+~ω)−f

(E)

~ ω = ^∂f

_∂E ^(E) , σ _ij (ω) = e ² πV ~

Z

|g(E)| ² hE|~v i |EihE|~v j |Ei



− ∂f _{F D} (E)

∂E



dE. (2.29) Finally, if we consider the case of a very low temperature T → 0, the deriva- tive of the Fermi-Dirac function can be substituted by the delta function

∂f

(E)

∂E = δ(E − E F ), which simplifies the integration over E. Thus, re- calling that g(E) = Tr[E − H], we obtain the relation for the conductivity tensor

σ ij = e ² π~

V Tr [v i δ(E F − H)v j δ(E F − H)] . (2.30) This equation is a starting point of the numerical real-space time-dependent Kubo method which will be described in details in the Chapter 5.

2.3 Landauer approach

Landauer approach belongs to the viewpoint 2, i.e. a constant current is forced to flow through a scattering system and the asking question is what the resulting potential distribution will be due to the spatially inho- mogeneous distribution of scatters [15]. Calculation of the current in the Landauer approach requires to divide formally the system into three parts, namely, the perfect leads and the scattering region, as depicted on Fig.(2.1).

The leads, in turn, are connected to the infinite reservoirs which represent the infinity and contain many electrons in a local equilibrium character- ized by the Fermi-Dirac distribution function. The basic idea behind this approach is that the electron has a certain probability to transmit through the scattering region [16]. Hence, the current carrying by an electron in a state with a wave-vector k is

J _k = ev(k)T (k) (2.31)

Sample (scattering region) Left

lead

Right

Reservoir lead Reservoir

a)

b)

0 L

Figure 2.1: a) Schematic illustration of the system which consists of the scattering region connected to perfect leads. The leads itself are coupled to microscopic reservoirs. b) Potential profile. The left and right leads have a constant potential µ L and µ R , respectively, equal to the chemical potential in the reservoirs.

with T (k) being a transmission probability. Full current supplied by the left lead is a sum over all states

I L = 2 X

k

J k f F D (E(k), µ L ), (2.32) where the factor 2 is due to spin-degeneracy, µ _L is a chemical potential in the left lead and

f _{F D} (E(k), µ) = 1 1 + exp 

E(k)−µ k

T

 (2.33)

is the Fermi-Dirac distribution function. Hence, we have I L = 2e X

k

v(k)T (k)f F D (E(k), µ L ) =

"

X

k

→ 1 2π

Z dk

# (2.34)

= 2e 2π

Z _∞

0 v(k)T (k)f F D (E(k), µ L )dk. (2.35) Changing integration variables by dk = _dE ^dk dE and using the expression for the group velocity v = ¹ _~ ^dE _dk , we get

I L = 2e h

Z _∞

U

T (E)f F D (E, µ L )dE. (2.36) Similarly, the current supplied by the right lead is

I _R = − 2e h

Z ∞ U

T (E)f F D (E, µ _R )dE. (2.37)

S-matrix technique 9

Sum of both contributions gives the net current I = I _L + I _R = 2e

h Z ∞

U

T (E) [f F D (E, µ _L ) − f F D (E, µ _R )] dE. (2.38) The potential drop between the reservoirs is

eV = µ _L − µ R . (2.39)

In the case of very low bias, the Fermi-Dirac functions can be expanded in the Taylor series

f _{F D} (E, µ _L ) − f F D (E, µ _R ) ≈ −eV ∂f _{F D} (E, µ)

∂E , (2.40)

which results in I = 2e

h Z ∞

U

T (E)



−eV ∂f F D (E, µ)

∂E



dE. (2.41)

This equation allows to calculate conductance of the system G = I

V = 2e ² h

Z ∞ U

T (E)



− ∂f F D (E, µ)

∂E



dE. (2.42)

At a very low temperature the derivative of the Fermi-Dirac distribution function can be replaced by the Dirac delta function δ(E −µ) which reduces Eq.(2.42) to

G = 2e ²

h T (µ). (2.43)

This equation shows that the conductance of the perfect conductor (i.e.

T = 1) is finite and thus the resistance (G ⁻¹ ) is non-zero. The following explanation can be used [16]: in the contacts (reservoirs) the current is carried by infinitely many transverse modes, however inside the conductor only few modes supply the current. It leads to redistribution of the current among current-carrying modes which results in the interface resistance.

2.4 S-matrix technique

As it was shown in the previous section, the current (or conductance) can be formulated in terms of the transmission function T . The powerful method to calculate T is the scattering matrix technique. Scattering matrix (or S- matrix) relates the outgoing amplitudes b = (b ₁ , b ₂ , ..., b _n ) to the incident amplitudes a = (a 1 , a 2 , ..., a n ), see Fig.(2.1),

b = S(E)a, S(E) =

 r(E) t

(E) t(E) r

(E)



. (2.44)

S matrix has a 2N · 2N dimension, where N is a number of transmission channels. The transmission probability now equals to

t _m←n (E) = |S mn | ² (2.45)

Prior to calculation of the S-matrix, one can determine its general prop- erties. S-matrix must be unitary, that is a consequence of current conser- vation: the incoming electron flux P

n |a| ² must be equal to the outgoing flux P

n |b| ²

b ⁺ b = a ⁺ a,

a ⁺ (1 − S ⁺ S)a = 0, (2.46)

S ⁺ S = I.

Moreover S-matrix is also a symmetric matrix, S = S ^T . This fact reflects the time-reversal symmetry of the Schr¨odinger equation, H = H ^∗ . A non- zero magnetic field breaks the time-reversal symmetry. In this case, we have S _{B ~} = S ^T

− ~ B .

S-matrix and Green’s function S-matrix can be expressed in terms of Green’s function. Outside the scattering region solution of the Schr¨odinger equation has the form of plane waves ψ _n (~r) ≈ e ^ik

^z , where we assume that the system with a cross-section area A is uniform in x and y directions, ~r = (~ρ, z). Each plane wave corresponds to a scattering channel n characterized by a transverse momenta ~q _n and longitudinal momenta k _n with the energy E = (1/2m)(k _n ² + q ² _n ). If we define the Green’s function G(E) = (E + iη − H) ⁻¹ with matrix elements between scattering channels m and n as

G _mn (z, z ^′ , E) = A ⁻¹ Z

d~ρ Z

d~ ρ ^′ exp (−i~q m ~ρ) exp (−i~q n ρ ~ ^′ )h~r|G(E)|~r ^′ i.

(2.47) Then, the transmission coefficient, following Fisher and Lee [17], can be calculated as

t _mn = −i~ √

v _m v _n G _mn (z, z ^′ , E) exp [−i(k m z − k n z ^′ )], (2.48)

where z and z ^′ are taken outside the scattering region, i.e. z > L and

z ^′ < 0, see Fig.(2.1), and v _n = k _n /m is the velocity in channel n. This

relation is very important, since it shows the connection between different

transport formalisms.

Chapter 3

Electron-electron interactions

With advent of quantum mechanics, the physical laws which govern par- ticles motion and interactions between particles became known. However the exact analytical solution is possible only for a system consisting of two particles. A typical piece of solid consists of approximately 10 ²³ particles.

Even if it would be possible to write down all the differential equations required to describe this system, the solution of these equation is an un- feasible task in principle. The problem of finding the solution arises from electron-electron interaction which makes the motion of particles correlated and couples corresponding differential equations. Therefore it is of great importance to develop approximated methods which provide a simplified form of the electron-electron interaction and reduces the number of equa- tions needed to be solved.

3.1 The many body problem

The Hamiltonian of a many-body system of interacting particle is written as

H = ˆ T ˆ _e + ˆ T _n + ˆ V _e−e + ˆ V _e−n + ˆ V _n−n

= − X

i

~ ²

2m _e ∇ ² i − X

I

~ ²

2M _I ∇ ² I + 1 4πεε ₀

X

i>j

e ²

|~r i − ~r j | (3.1)

+ 1

4πεε ₀ X

I>J

Z _I Z _J e ²

| ~ R _I − ~ R _J | − 1 4πεε ₀

X

i,I

Z _I e ²

|~r i − ~ R _I | ,

where the first two terms describe the kinetic energy of electrons and nuclei.

The last three terms result from the Coulomb interaction between electrons, electron-nuclei and nuclei-nuclei respectively. The Hamiltonian acts on the wave-function Ψ({~r i }, { ~ R I }) which depends on the position of all electrons and nuclei in the system

HΨ({~r ˆ i }, { ~ R _I }) = EΨ({~r i }, { ~ R _I }). (3.2)

An essential simplification of Eq.(3.2) can be done with the use of the Born- Oppenheimer approximation which neglects the coupling between the nu- clei and electronic motion. In thermodynamic equilibrium electrons move much faster than nuclei, since M ≫ m e , which allows to treat nuclei as stationary particles and neglect the kinetic term ˆ T _n in the Hamiltonian.

Hence one can deal only with the electronic part, ˆ H e , of the full Hamil- tonian which corresponds to the system of interacting electrons moving in the effective potential produced by nuclei

H ˆ _e = ˆ T _e + ˆ V _e−e + ˆ V _e−n + ˆ V _n−n , (3.3) H ˆ e Φ({~r i }) = EΦ({~r i }). (3.4) Even though electronic wave-function Φ({~r i }) still depends on the positions of nuclei, ~ R _I are just parameters of Eq.(3.3) and the number of differential equations needed to be solved is greatly reduced.

3.2 Hartree-Fock approximation

The Born-Oppenheimer approximation significantly simplifies the problem of interacting particles by eliminating the coupling between electrons and nuclear motion. However determination of the exact solution of the many- particle electronic wave-function is still not feasible. The basic idea of the Hartree-Fock approximation is to substitute the system of interacting electrons by the motion of single electrons in the average self-consistent field generated by all the other electrons in the system.

The Hamiltonian of the many-particle system is given by H = ˆ X

k

p ² 2m _e + X

k

V (~r _k )+ e ² 8πεε ₀

X

kk

1

|~r k − ~r k

| = X

k

H ˆ _k + X

kk

H ˆ _kk

, (3.5)

where the term V (~r _k ) = P

I V (~r _k − ~ R _I ) describes interaction between k- electron with all nuclei in the system located at { ~ R _I }. The operator ˆ H _k is a one-particle operator, while ˆ H _kk

depends on the position of two particles.

The simplest way to construct the many-particle wave-function is to write down it in the form of a product of single-particle wave-functions, φ _k (~r), which have to be determined,

Φ({~r k }) = Y

k

φ(~r _k ), hφ i (~r)|φ j (~r)i = δ ij . (3.6)

Let us calculate an expectation value of the energy [14]

E = hΦ| ˆ H|Φi = X

k

hφ k | ˆ H k |φ k i + e ² 8πεε ₀

X

kk

 φ k φ k

1

~r _k − ~r k

    φ k φ k

.

(3.7)

Hartree-Fock approximation 13

According to the variational principle the closer values of φ _k to the exact solution the smaller value of the energy, thus

δ E − X

k

E _k (hφ k |φ k i − 1)

!

= 0, (3.8)

where E _k are Lagrange parameters. Changing φ _i → φ i + δφ _i and keeping terms linear in respect to δφ i , we arrive at the Hartree equation

"

− ~ ²

2m _e ∆ + V (~r) + e ² 4πεε ₀

X

k6=i

Z |φ k (~r k )| ²

|~r k − ~r i | d~r _k

#

φ _i (~r) = E _i φ _i (~r). (3.9) The third term in the brackets has a simple interpretation. If we define the charge density as n(~r) = e P

i |φ i (~r)| ² , then the term U H (~r) = e

4πεε 0

Z n(~ r ^′ )

|~r ^′ − ~r| d~ r ^′ , (3.10) called the Hartree term, describes the Coulomb interaction between the i-th electron located at ~r with all the other electrons in the system.

Since electrons are fermions, the many-particle wave-function must change the sign under the interchange of the coordinates of any two particles. The wave-function given by Eq.(3.6) does not satisfy this condition. In order to construct an antisymmetric wave-function one can use Slater determinant

Φ({~q k }) = 1

√ N!

φ ₁ (~q ₁ ) . . . φ _N (~q ₁ ) .. . . .. .. . φ ₁ (~q _N ) . . . φ _N (~q _N )

, (3.11)

where ~ q _i = {~r i , σ _i } denotes both position and spin of the electron and the factor √ ¹

N ! is used for normalization. Following the same way as before, i.e.

applying the variational principle to the expectation value of the energy, the new form of the wave-function results in the equation [14]



− ~ ²

2m _e ∆ + V (~r)



φ i (~r) + e ² 4πεε ₀

X

k6=i

Z |φ k (~ r ^′ )| ²

|~r ^′ − ~r| d~ r ^′ φ i (~r)

− e ² 4πεε 0

X

k6=i

Z φ ^∗ _k (~ r ^′ )φ _i (~ r ^′ )

|~r ^′ − ~r| d~ r ^′ · φ k (~r) = E _i φ _i (~r), (3.12)

called the Hartree-Fock equation. The additional term in Eq.(3.12) is

known as the exchange interaction. It does not have classical analog and

results from the Pauli exclusion principle. The exchange interaction term

which arises in the Hartree-Fock approximation has a non-local form in

contrast to the Coulomb interaction. This makes calculations more com-

plicated. In the density-functional theory (discussed in Sec.(3.3)) a number

of approximations are used to deduce a local form of the exchange interac-

tion.

3.3 Density-functional theory

Density-functional theory (DFT) is one of the most widely used model- ing method applied in physics and chemistry for calculation of electronic properties of complex systems. The basic idea behind DFT is to describe the system in terms of the electronic density instead of operating with a many-body wave function [18].

Hohenberg-Kohn theorems

In 1964 Hohenberg and Kohn proved two theorems which made the DFT possible [19]. They state that a knowledge of the ground-state density can, in principle, determine all the ground-state properties of a many-body system [18].

Theorem 1: An external potential V _ext (~r) uniquely determines the elec- tronic density for any system of interacting particles.

Proof: Assume that the same electron density n(~r) results from two po- tentials V _ext ¹ (~r) and V _ext ² (~r) differing by more than constant. Obviously, V _ext ¹ (~r) and V _ext ² (~r) belong to distinct Hamiltonians ˆ H ¹ (~r) and ˆ H ² (~r) which produce different wave-functions Ψ ¹ (~r) and Ψ ² (~r). The ground-state state energy associated with the Hamiltonian ˆ H ¹ (~r) is

E ¹ = hΨ ¹ | ˆ H ¹ |Ψ ¹ i. (3.13) According to the variational principle no other wave-function can give lower energy, i.e.

E ¹ = hΨ ¹ | ˆ H ¹ |Ψ ¹ i < hΨ ² | ˆ H ¹ |Ψ ² i (3.14) Since the Hamiltonians differs by the external potentials only, we can write H ˆ ¹ = ˆ H ² + V _ext ¹ − V ext ² , which gives us for the expectation value

hΨ ² | ˆ H ¹ |Ψ ² i = hΨ ² | ˆ H ² |Ψ ² i + Z 

V _ext ¹ − V ext ²

 n(~r)d~r. (3.15) Substituting it into Eq.(3.14) and recalling that E ² = hΨ ² | ˆ H ² |Ψ ² i, we obtain

E ¹ < E ² + Z 

V _ext ¹ − V ext ²

 n(~r)d~r. (3.16) Interchanging labels (1) and (2), we find in the same way that

E ² < E ¹ + Z 

V _ext ² − V ext ¹

 n(~r)d~r. (3.17) Addition of Eq.(3.16) and Eq.(3.17) leads to contradiction

E ¹ + E ² < E ¹ + E ² . (3.18)

Hence the theorem is proved by reductio ad absurdum.

Kohn-Sham equations 15

Theorem 2: The exact ground-state density n(~r) is the global minimum of the universal functional F [n].

Proof: Since the electron density n(~r) uniquely determines wave-function Ψ, the universal functional F [n] can be defined as

F [n] = hΨ| ˆ T + ˆ U _e−e |Ψi. (3.19) For a given external potential V _ext (~r), the energy functional is written as

E[n] = F [n] + Z

V _ext (~r)n(~r)d~r. (3.20) According to the variational principle, it has a minimum only for the ground-state wave-function Ψ. For any other wave-function Ψ ^′ which pro- duces density n ^′ (~r), we get

E[Ψ] = F [n]+

Z

V _ext (~r)n(~r)d~r < F [n ^′ ]+

Z

V _ext (~r)n ^′ (~r)d~r = E[Ψ ^′ ]. (3.21)

3.4 Kohn-Sham equations

In the Kohn-Sham method [20] one considers a system of non-interacting electrons moving in some effective potential v _{ef f} (~r) (which will be defined later)



− ~ ²

2m ∇ ² + v _{ef f} (~r)



φ _i (~r) = ǫ _i φ _i (~r). (3.22) An obtained set of the Kohn-Sham orbitals φ _i determines the electron density

n(~r) = X

i

|φ i (~r)| ² . (3.23)

The energy functional equals to E[n] = T _S [n] + V _{ef f} [n] = T _S [n] +

Z

n(~r)v _ext (~r)d~r + V _H [n] + E _xc [n], (3.24) where the first term,

T s [n] = X N

i=1

Z φ ^∗ _i (~r)



− ~ ² 2m ∇ ²



φ i (~r)d~r, (3.25) is a single-electron kinetic energy functional. The second term describes the potential energy acquired by the charged particles in an external electric field. The Hartree term is given by

V H = e ² 8πεε ₀

Z Z n(~r)n(~r ^′ )

|~r − ~r ^′ | d~rd~ r ^′ . (3.26)

The last term, E _xc [n], arises from the exchange-correlation interaction. It can be shown [21] that the functional given by Eq.(3.24) corresponds to the the effective potential in the form

v ef f (~r) = v ext (~r) + e 4πεε 0

Z n(~ r ^′ )

|~r − ~r ^′ | d~ r ^′ + v xc (~r). (3.27) Equations (3.22), (3.23) and (3.27) constitute the basis of the Kohn-Sham method, and are solved self-consistently for the ground state density and the effective potential.

The explicit form of the exchange-correlation potential v _xc can be de- termined using other approximations. The most widely used are the local spin density approximation [21] and the generalized gradient approximation [22].

3.5 Thomas-Fermi-Dirac approximation

According to the Hohenberg-Kohn theorems discussed in Sec.(3.3), the total energy of the system may be written as

E[n(~r)] = Z

T [n(~r)]d~r + Z

V (~r)n(~r)d~r, (3.28) where T [n(~r)] is a kinetic energy functional of the electron density n(~r) and V (~r) is an external potential. The Thomas-Fermi approximation assumes that the kinetic-energy functional is a local function of the density. This assumption allows to rewrite Eq.(3.28) in the form

µ = T [n(~r)] + V (~r). (3.29) Let us now derive the relation between the potential and charge density in graphene. Taking into account dispersion relation for graphene, E =

±~v F |~k|, and using n = g v g _s R _d~k

(2π)

(v _F = 10 ⁶ m/s and g _v = g _s = 2, see Chapter 4 for details), we get [28]

sgn[n(~r)]~v _F p

πn(~r) + V (~r) = µ. (3.30) An alternative way is to rewrite Eq.(3.28) in the form which directly relates electron density to the external potential [23, 24]

n(~r) = Z

dEρ(E − V (~r))f F D (E, µ), (3.31)

where ρ(E) is a single-electron density of states, calculated in the presence

of homogeneous potential. Figure (3.1) illustrates application of Eq.(3.31).

Hubbard model 17

Figure 3.1: Schematic illustration of Thomas-Fermi model.

Locally the dispersion relation corresponding to the ideal system (with homogeneous external potential) is preserved. Filling up the states lying between the charge neutrality point and the Fermi energy level one obtains local electron density.

Even though the Thomas-Fermi model misses quantum mechanical ef- fects (e.g. quantization), it produces quantitatively similar results in com- parison to more rigorous models and widely used in graphene physics [25, 26, 27, 28].

3.6 Hubbard model

The Hubbard model, originally proposed by John Hubbard in 1963 [29], is the simplest model of interacting particles in a lattice. The interaction is assumed to take place only between particles located at the same site (or atom). Despite of its simplicity rigorous analytical solution is found only for a one-dimensional problem [30].

The Hubbard Hamiltonian in the tight-binding approximation consists of two terms, namely the kinetic energy term and the on-site potential

H = −t ˆ X

hi,ji,σ

(a ⁺ _i,σ a j,σ + h.c.) + U X

i

n _i↑ n _i↓ , (3.32)

where n _i,σ = a ⁺ _i,σ a _i,σ is the occupation number operator. The parameter U = const describes the strength of on-site Coloumb interaction and can be determined using ab-initio calculations or extracted from experimental data.

Equation (3.32) can be rewritten within the mean-field approach using substitutions

n _i↑ = hn i↑ i + (n i↑ − hn i↑ i), (3.33)

n _i↓ = hn i↓ i + (n i↓ − hn i↓ i), (3.34)

where hn iσ i denotes average occupation of spin σ at site i. Hence, we have n _i↑ n _i↓ = n _i↑ hn i↓ i + n i↓ hn i↑ i − hn i↓ ihn i↑ i + (n i↑ − hn i↑ i)(n i↓ − hn i↓ i)

| {z }

≈0

, (3.35)

where the product of two deviations from the average values is assumed to be small. Substituting the result of Eq.(3.35) into Eq.(3.32) we derive Hubbard Hamiltonian in the mean-field approximation

H ˆ ^{M F} = −t X

hi,ji,σ

(a ⁺ _i,σ a _j,σ + h.c.) + U X

i

(n _i↑ hn i↓ i + n i↓ hn i↑ i − hn i↓ ihn i↑ i) .

(3.36)

Despite of its simplicity the Hubbard model was applied to investigate the

properties of different materials. It reproduces a variety of phenomena

observed in solid state physics, such as ferromagnetism, metal-insulator

transition and superconductivity.

Chapter 4

Electronic structure and transport in graphene

4.1 Basic electronic properties y

x A

B a) b)

Figure 4.1: a) Graphene lattice consisting of two interpenetrating triangu- lar sublattices A (red circles) and B (blue circles) with unit vectors ~a 1 , ~a 2

and nearest-neighbours vectors ~δ ₁ , ~δ ₂ , ~δ ₃ . The yellow parallelogram marks unit cell containing two atoms. b) The structure of reciprocal lattice de- fined by unit vectors ~b 1 and ~b 2 . The grey hexagon is the first Brillouin zone.

Real space and reciprocal lattices Carbon atoms in graphene are arranged in a honeycomb lattice shown on Fig.(4.1). This structure can be described [31, 32] as a triangular lattice with unit vectors

~a ₁ = a cc

2 (3, √

3),~a ₂ = a cc

2 (3, − √

3), (4.1)

where a _cc = 0.142 nm is a carbon-carbon distance. Unit cell contains two atoms and has an area S _cell = ^{a √} ₄ ³ , where a = a _cc √

3 = 0.246 nm is a lattice constant. Each point of the sublattice A is connected to its nearest- neighbors by the vectors

δ ~ 1 = a cc

1 2 ,

√ 3 2

!

, ~ δ 2 = a cc

1 2 , −

√ 3 2

!

, ~ δ 3 = a cc (−1, 0) . (4.2)

For a given lattice one can easily build a reciprocal lattice with a unit vectors ~b i defined by the relation ~b i ~a j = 2πδ ij , i.e.

~b ₁ = 2π 3a _cc (1, √

3),~b ₂ = 2π 3a _cc (1, − √

3). (4.3)

The first Brillouin zone, which is the Wigner-Seitz primitive cell of the reciprocal lattice [33], is shown on Fig.(4.1)(b). There are two inequivalent points which are of special interest in graphene physics

K = ~ 2π 3a _cc

 1, 1

√ 3



, ~ K ^′ = 2π 3a _cc

 1, − 1

√ 3



. (4.4)

Dispersion relation Graphene lattice can be considered as two inter- penetrating triangular sublattices A and B which are defined by vectors

R ~ _p,q ^A = p~a 1 + q~a 2 , ~ R ^B _p,q = ~δ 1 + p~a 1 + q~a 2 , (4.5) where p, q are integer numbers.

In a single-electron approximation the tight-binding Hamiltonian for electrons in graphene is given by

H ˆ _tb = −t X

p,q

a ⁺ _p,q b _p,q + a ⁺ _p,q b _p−1,q + a ⁺ _p,q b _p−1,q+1

+ h.c., (4.6)

where a ⁺ _p,q (a _p,q ) and b ⁺ _p,q (b _p,q ) create (annihilate) an electron on sublattices A and B at site ~ R ^A _p,q and ~ R ^B _p,q respectively and t = 2.77 eV is a nearest- neighbor hopping integral.

The wave-function for the lattice can be written in the form

|Ψi = X

p,q

ζ _p,q ^A a ⁺ _p,q + ζ _p,q ^B b ⁺ _p,q

|0i, (4.7)

where ζ p,q ^A(B) is a probability amplitude to find the electron at site ~ R ^A(B) p,q

Substituting Eqs.(4.6),(4.7) into the Schr¨odinger equation, ˆ H _tb |Ψi = E|Ψi,

and calculating the matrix elements h0|a p,q H ˆ _tb a ⁺ _p,q |0i and h0|b p,q H ˆ _tb b ⁺ _p,q |0i

Basic electronic properties 21

with use of the commutation relation, one arrives to the system of difference equations

−t ζ p,q ^B + ζ _p−1,q ^B + ζ _p−1,q+1 ^B

= Eζ _p,q ^A , (4.8)

−t ζ p,q ^A + ζ _p+1,q ^A + ζ _p+1,q+1 ^A

= Eζ _p,q ^B . The states ζ _p,q ^A and ζ _p,q ^B can be written in the Bloch form

ζ _p,q ^A = ψ _p,q ^A e ^{i~k ~ R}

, ψ _p,q ^A = ψ ^A _p+1,q = ψ _p+1,q−1 ^A , (4.9) ζ _p,q ^B = ψ ^B _p,q e ^{i~k ~ R}

, ψ ^B _p,q = ψ ^B _p−1,q = ψ _p−1,q+1 ^B .

Substituting Eq.(4.9) and Eq.(4.2) into Eq.(4.8) and omitting indexes (p, q), one gets

−tφ(~k)ψ ^B = Eψ ^A , (4.10)

−tφ ^∗ (~k)ψ ^A = Eψ ^B , where

φ(~k) ≡ e ^{i~k ~ δ}

+ e ^{i~k ~ δ}

+ e ^{i~k ~ δ}

, (4.11) or in a matrix form

H ˆ

 ψ ^A ψ ^B



= E

 ψ ^A ψ ^B



, ˆ H ≡ 0 −tφ(~k)

−tφ ^∗ (~k) 0

!

. (4.12)

In order to obtain dispersion relation, one needs to determine the eigenval- ues of the matrix ˆ H, which are calculated using the relation det | ˆ H − ˆ IE| = 0,

E(~k) ² = t ² |φ(~k)| ² , E(~k) = ±t

v u

u t 1 + 4 cos ²

√ 3 2 a _cc k _y

! + 4 cos

√ 3 2 a _cc k _y

! cos  3

2 a _cc k _x

 . (4.13)

Dirac equation The spectrum of graphene, given by Eq.(4.13), is sym-

metric in respect to energy E = 0. If the Fermi energy coincides with

this point (E _F = 0), i.e the states are occupied only up to zero energy,

it corresponds to the case of electrically neutral graphene. One is usually

interested in electronic properties close to a charge-neutrality point. There

are six points in the k-space where the energy equals to zero. These points

are at the corners of the first Brillouin zone, see Fig.(4.1). Only two of

them, ~ K and ~ K ^′ , are inequivalent.

0

a) 3 b)

Figure 4.2: a) Dispersion relation calculated using Eq.(4.13). The energy is given in units of the hopping integral. b) Close to the charge neutrality point dispersion relation is linear and has a form of a cone determined by the Dirac equation (see Eq.(4.23) below).

Let’s expand the function φ(~k) in the Hamiltonian of Eq.(4.12) near the K-point ~

~k = ~ K + ~q, (4.14)

where ~q is some small (|~q| < | ~ K|) vector having origin at ~ K. Substituting Eq.(4.14) in φ(~k), one gets

φ( ~ K + ~q) = X 3

i=1

e ^{i ~ K~ δ}

e ^{i~ q~ δ}

= e ⁱ

e ^{i~ q~ δ}

+ e ⁱ⁰ e ^{i~ q~ δ}

+ e ⁻ⁱ

e ^{i~ q~ δ}

. (4.15)

Considering the continuum (low energy) limit (a _cc → 0) [32], one can ex- pand exponents in Taylor series

a lim

→0 e ^{i~ q~ δ}

≃ 1 + i~q~δ i . (4.16) After performing some straightforward algebra, we arrive at

φ(~k) = 3a _cc 2

"

−

√ 3 2 − i 1

2

!

q x − 1 2 + i

√ 3 2

! q y

#

. (4.17)

If we rotate the system of coordinates on angle θ = ^π ₆ by operator

R =

 cos θ − sin θ sin θ cos θ



=

√ 3 2 − ¹ ₂

1 2

√ 3 2

!

, (4.18)

Basic electronic properties 23

Eq.(4.15) is finally reduced to

φ(~q) = − 3a _cc

2 (q _x + iq _y ). (4.19)

Substituting Eq.(4.19) into Eq.(4.12), one gets 3

2 a _cc t

 0 q x + iq y

q _x − iq y 0

  ψ ^A ψ ^B



= E

 ψ ^A ψ ^B



. (4.20)

With the use of the Pauli matrices, ~σ = (σ _x , σ _y ), where σ x =

 0 1 1 0

 , σ y =

 0 −i i 0



, (4.21)

the Hamiltonian can be written in a vector form

H ˆ _K = ~v _F ~σ~q, (4.22)

where v _F = ³ ₂ ^a

_~ ^t ≃ 10 ⁶ m/s is a Fermi velocity. Equation (4.22) is al- gebraically identical to a two-dimensional relativistic Dirac equation with vanishing rest mass known as Weyl’s equation for a neutrino, where the two-component wave function (or spinor) represents pseudo-spin which re- sults from the presence of two sublattices [34, 35]. The eigenenergies of ˆ H _K are

E = ±~v F q. (4.23)

In order to find eigenfunctions Eq.(4.20) can be rewritten in the following way

~ v _F q

 0 e ^{iθ(~ q)} e ^−iθ(~q) 0

  ψ ^A ψ ^B



= E

 ψ ^A ψ ^B



, (4.24)

where θ(~q) = arctan(q y /q x ). Taking into account Eq.(4.23) for eigenener- gies and using normalization condition |ψ ^A | ² + |ψ ^B | ² = 1, one gets

|Ψ K,s ~ (~q)i = 1

√ 2

 e ^{iθ(~ q)/2} se ^−iθ(~q)/2



, (4.25)

where the sign s = ± corresponds to the eigenenergies ±~v F q.

Besides spin-degeneracy each level is double-degenerated due to valley.

One has to operate by a full wave function which includes the contribution from both valleys. The same procedure can be repeated to obtain the effective Hamiltonian and wave function near K ^′ -point

H ˆ _K

= ~v _F ~σ ^∗ ~q, |Ψ K ~

,s (~q)i = 1

√ 2

 e ^−iθ(~q)/2 se ^{iθ(~ q)/2}



. (4.26)

If the wave-vector ~q rotates once around the Dirac point, i.e. θ → θ + 2π, the wave-function acquires an additional phase equals to π, hence Ψ K,s ~ (θ ± 2π) = −Ψ K,s ~ (θ) which is a characteristics of fermions.

Electron’s wave function in graphene has a chiral nature. Helicity can be interpreted as a projection of pseudospin vector on direction of motion and defined by the operator ˆh = ~σ~q/|~q|. It can be easily shown using Eq.(4.22), that eigenvalues of the helicity operator equal to h = ±1. Since ˆh commutes with the Hamiltonian, helicity is a conserved quantity and responsible for such effects as the Klein tunneling [36].

4.2 Graphene nanoribbons

Electronic properties of graphene nanoribbons (GNR) depend on the type of a edges. One can distinguish two types, namely, zig-zag and armchair GNR’s.

-3 -2 -1 0 1 2 3

E/t 0

-π ka x π

a) b)

L

Figure 4.3: a) Structure of an armchair graphene nanoribbon lattice. Each edge of the ribbon is terminated by both A and B atoms. b) Dispersion relation of the ribbon with 10 atoms in transverse direction calculated using the tight-binding Hamiltonian.

Armchair graphene nanoribbons Dispersion relation of armchair GNR can be derived solving the Schr¨odinger equation in the following form [37]

~ v _F



 



0 k x + ik y 0 0

k _x − ik y 0 0 0

0 0 0 −k x + ik _y

0 0 −k x − ik y 0



 





 

 ψ A

ψ _B ψ _A ^′ ψ _B ^′



 

 = E



 

 ψ A

ψ _B ψ ^′ _A ψ ^′ _B



 

 ,

(4.27)

Graphene nanoribbons 25

where ψ _A(B) and ψ _A(B) ^′ are the probabilty amplitudes on the sublattice A(B) for the state near the ~ K and ~ K ^′ points, respectively. The total wave function has the form

Ψ = e ^{i ~ K~ r} Ψ K ~ + e ^{i ~ K}

^{~ r} Ψ ^′ _{K ~} . (4.28) Let us first find a solution near the K point. Substituting in Eq.(4.22) wave vector ~k = 

−i _∂x ^∂ , −i _∂y ^∂ 

, one gets

 0 −i _∂x ^∂ + _∂y ^∂

−i _∂x ^∂ − _∂y ^∂ 0

  ψ ^A ψ ^B



= ǫ

 ψ ^A ψ ^B



, (4.29) where ǫ = _~ ^E _v

F

. Due to translational invariance in ~x-direction the wave function can be written in the form

Ψ _{K ~} (x, y) = e ^ik

^x

 φ ^A (y) φ ^B (y)



, (4.30)

which allows us to reduce the problem to a system of two differential equa- tions

( k x φ B (y) + ^∂φ _∂y

^(y) = ǫφ A (y) k _x φ _A (y) − ^∂φ _∂y

^(y) = ǫφ _B (y) ⇒

( φ _B (y) = ^k _ǫ

φ _A (y) − ¹ _ǫ ^∂φ _∂y

^(y)

∂

φ

(y)

∂y

+ z ² φ _A (y) = 0

(4.31) where z ² = ǫ ² − k ² x . The general solution of the system of equations is a sum of plane waves

 φ _A (y) = Ae ^izy + Be ^−izy

φ _B (y) = ^k

^−iz _ǫ Ae ^izy + ^k

^+iz _ǫ Be ^−izy (4.32) Similar derivation can be done for the wave functions describing the states near K ^′ -point, which gives

 φ ^′ _A (y) = Ce ^izy + De ^−izy

φ ^′ _B (y) = ^−k

_ǫ ^−iz Ce ^izy + ^−k

_ǫ ^+iz De ^−izy (4.33) In order to find the unknown coefficients A and B, one can utilize the boundary conditions. The armchair nanoribbon edge consist of atoms be- longing to both sublattices, see Fig.(4.3)(a), therefore one can expect that both Ψ _A and Ψ _B should vanish at the edges

Ψ _A (0) = Ψ _B (0) = Ψ _A (L) = Ψ _B (L) = 0, (4.34) where Ψ _A(B) is an A(B) component of the total wave function (4.28).

Hence, we have

φ _A (0) + φ ^′ _A (0) = 0 φ _B (0) + φ ^′ _B (0) = 0 e ^iK

^L φ _A (L) + e ^−iK

^L φ ^′ _A (L) = 0 e ^iK

^L φ _B (L) + e ^−iK

^L φ ^′ _B (L) = 0

(4.35)