3.1  Tight-Binding Method

During the past decade, graphene (Castro Neto et al. 2009; Geim and Novoselov 2007; Novoselov et al. 2004, 2005), a single layer of graphite with a planar honeycomb structure as illustrated in Figure 3.1, has been extensively studied because of its astonishing properties that are mostly attributed to quantum phenomena from 2D confinement effects (Novoselov et al. 2004, 2005; Zhou et al. 2006). In a perfect graphene, each carbon atom forms σ-bonds with its three nearest neighbors (sp² hybridization). The electronic states near the Fermi energy are dominated by the π and π* bands, which are derived from the weakly interacting p_z orbitals. The most important characteristic of the electronic structure of graphene is the degenerate states π and π* at the K point of graphene (hexagonal) Brillouin zone (BZ) (inset of Figure 3.2), making graphene a zero-bandgap semi-metal (Reich et al. 2002; Wallace 1947).

The tight-binding method is an approximation assuming that the wave functions tightly bound to the atoms, so-called “tight-binding,” such that the atomic wave functions can be used as a basis for expanding the crystal wave functions. The method is helpful to investigate how chemical bonding, atomic arrangement, and structural symmetry reflect to the characteristic of electronic structure. This method is the simplest method for computationally calculating band structures. Since each carbon atom forms bonds with its three nearest neighbors, it is assumed that the wave functions of s, p_x, and p_y orbitals are tightly bound to the atoms. The method has been widely used and proven to be efficient for studying this class of materials as summarized in Table 3.1. Generally, suppose that we start with an atomic wave function, for example, p_z orbital, on a sublattice A (Figure 3.1), which is centered at coordinate R_A,

ϕ (r − R_A) is an atomic wave function associated with this atom. It is assumed that this state interacts with an atomic wave function on a sublattice B (Figure 3.1), which is centered at site R_B,

These two orbitals dominantly attribute to π and π* bands near the Fermi level. As long as we consider about the electronic structures near the Fermi level, we can form electronic states that can be used as the basis functions for the crystal wave functions assuming that there are no other orbitals, for example, s, p_x, and p_y orbitals can mix into the states near the Fermi level since the bands that correspond to the dispersion of bonding and antibonding molecular orbital are π and π* bands. Given that there are two atoms in the unit cell at sublattices A and B, we can construct two Bloch states as

where i = A, B, with the summation running over all the N unit cells in the crystal (the vectors R). It can be verified that these states obey Bloch's theorem. Let R′ be another lattice vector,

Then, this state at position r + R is given by

where R′ − R is another lattice vector. This leads to

that satisfies Bloch's theorem. The crystal eigenstates can be expanded in these two basis functions, |k_A〉 and |k_B〉. We can construct the approximate crystal eigenstates as

where the expansion coefficients c_A and c_B are to be determined. The eigenstates |k〉 are normalized by |c_A|² + |c_B|² = 1. We want to find the eigenvalues ∊ (k) that are exactly the number of bands that we can expect at each k-point of a Hamiltonian operator $\widehat{H}$

such that

Inserting the unit operator $\widehat{I}={\sum }_j|k_j⟩⟨k_j|=|k_A⟩⟨k_A|+$

$|k_B⟩⟨k_B|$ in front of |k〉 on both sides of the equation and multiplying by 〈 k_i|, where i = A, B, we can write this equation in the form

In order to solve this system of linear equations, we need to evaluate 〈 k_i|k_j〉

It is obvious that R″ − R′ is another lattice vector. We can define R″ − R′ = R given that

The integral 〈R_i + R′|R_j + R + R′〉 = 〈R_i|R_j + R〉 is independent to the lattice vectors R′ because the relative position between R_i and R_j + R is unchanged under the same translation vector R′. The summation over the lattice vectors R′ will cancel with the factor 1/N because there is no explicit dependence on R′ so that

In the framework of the tight-binding approximation, the overlap integral 〈R_i|R_j + R〉 is nonzero only for the same orbital on the same atom, that is, only for i = j and R = 0, therefore,

The two basis functions,|k_A〉 and |k_B〉, are orthogonal. With this approximation, Equation 3.9 can be written as

with

This system of equations can have nonzero solutions only if its determinant vanishes,

This is known as the secular or characteristic equation. The Hamiltonian matrix $H_{ij}=⟨k_i\left|\widehat{H}\right|k_j⟩$

is N × N matrix, where N is the number of basis functions used for expanding the crystal eigenstates. Alternatively, the eigenvalues can be obtained by directly diagonalizing the Hamiltonian matrix with all the known values of elements H_ij. The eigenvalues of a symmetric matrix are real. The solutions of this equation yield N eigenvalues ∊₁(k), ∊₂(k), …, ∊_N(k), which are exactly the number of bands that we can expect for each value of k.

The Hamiltonian matrix elements H_ij between atomic states can be obtained by

$\begin{array}{c}⟨k_i\left|\widehat{H}\right|k_j⟩=\frac{1}{N}{\displaystyle \sum R^{\prime },R^{″}{e}^{ik\cdot (R^{″}-R^{\prime })}⟨R_i+R^{\prime }\left|\widehat{H}\right|R_j+R^{″}⟩}\\ =\frac{1}{N}{\displaystyle \sum R^{\prime },R{e}^{ik\cdot R}⟨R_i+R^{\prime }\left|\widehat{H}\right|R_j+R+R^{\prime }⟩}\end{array}$

At this point, an important approximation in the framework of the tight-binding method is introduced by taking the Hamiltonian matrix elements to be nonzero only if: (i) the orbital is on the same atom, that is, for i = j and R = 0,

where the diagonal elements ∊_i are referred to as the onsite energies or the energy offset for an atom at site i, or (ii) the orbital is on atoms at nearest-neighbor sites located at c_NN,

that is, R_i = R_j + R + c_NN. The summation running over all R will be nonzero when the nearest neighbors are in the same unit cell (R = 0) or when they are across unit cells (R is a combination of the primitive lattice vectors). The t_ij is referred to as the hopping integral between atoms at sites i and j. A range of nearest-neighbor sites can be included into hopping elements (Equation 3.19) for better description of Hamiltonian matrix but more hopping integrals t_ij have to be included to evaluate the matrix elements. In practice, we can consider ∊_i and t_ij as parameters or we can even consider all the matrix elements as parameters as well. In both cases, these parameters are fitted to reproduce certain properties from experimental data or first-principles calculations. Then these parameters can be used to calculate other properties. Alternatively, these parameters can be used as variables to investigate a change to the characteristic of electronic structure under influence of that parameter, that is, ∊_i, t_ij.

The Hamiltonian matrix elements H_ij in Equation 3.15 can be obtained by another form of Hamiltonian that is defined as

where ∊_i is the site energy for an atom at site i, t_ij is the hopping integral between atoms at sites i and j, and $a_i^{†}$

and a_i are the creation and annihilation operators, respectively, of π-electron at site i. By this definition, $a_i^{†}$ is equivalent to an atomic state|R_i〉, and a_i is equivalent to an atomic state 〈R_i|. j represents the index of summation that is taken over only the nearest neighbors of interest, typically first neighbors, with truncation of farther neighbors. The interactions with farther neighbors, higher-order correction, are ignored for simplicity in most qualitative interpretations. The abbreviation “h.c.” stands for Hermitian conjugate of the term $a_i^{†}{a}_j$ . Tight-binding method provides an insightful of how the interactions between the constituents influence on the electronic structures.

3.2  Electronic Structures of Pristine Graphene and Honeycomb Lattice

The lattice structure of graphene contains two sublattices per primitive unit cell as shown in Figure 3.1. The primitive unit cell is defined by the primitive lattice vectors,

where a is a lattice constant of graphene. With this choice of lattice vectors, the reciprocal primitive lattice vectors are given by

The lattice vectors R are formed by all possible combinations of lattice vectors

where n₁ and n₂ are integers. The lattice vectors connect all equivalent points in space. The two sublattices are usually referred to as sublattice “A” and “B” which are indicated by different colors in Figure 3.1. Given that there are two atoms that are located at the two sublattices in the primitive unit cell. Suppose that we choose the position of sublattice A to be at lattice point R_A, then sublattice B will be located at R_B = R_A + c with the basis vector

The electronic structures of pristine graphene can be qualitatively described by tight-binding model. The Hamiltonian is defined as

For simplicity, let us start by considering only the hopping integral or the interactions from first nearest neighbors. The summation is taken over only the first nearest neighbors. The interactions with farther neighbors are ignored for a moment. In the unit cell as shown in Figure 3.1, sublattice A is surrounded by three neighbors which are located at R_B ≡ R_A + c, R_B − a₁, and R_B − a₁ − a₂. Likewise, sublattice B is surrounded by three neighbors, which are located at R_A, R_A + a₁, and R_A + a₁ + a₂. Therefore, the Hamiltonian in Equation 3.25 can be explicitly written as

where|R_A〉 and|R_B〉 are the π atomic states associated with the atoms at sublattices A and B, respectively. ∊_A and ∊_B are the site energies for an atom at sites A and B, respectively. t_AB is the hopping integral between atoms at sites A and B. To evaluate the Hamiltonian matrix elements in Equation 3.15, we need to use two bases for expansion of the crystal wave functions which obey Bloch's theorem. These bases are given by Equation 3.3,

with the summation running over all the N unit cells in the crystal (the vectors R), since the Hamiltonian matrix as a function of k is defined by

or

The Hamiltonian matrix elements can be evaluated by using,

The first diagonal matrix element is given by

The overlap integral 〈R_i|R_j〉 = δ_ij is nonzero only for the same orbital on the same atom, that is, only for i = j; hence

The element $H_{AA}=⟨k_A\left|\widehat{H}\right|k_A⟩$

is

By using Equations 3.30 and 3.32, the off-diagonal element $H_{BA}=⟨k_B\left|\widehat{H}\right|k_A⟩$

is given by

$\begin{array}{l}⟨k_B\left|\widehat{H}\right|k_A⟩=-t_{BA}{\displaystyle \sum R(\frac{1}{N}⟨R_B|R_B⟩+\frac{1}{N}{e}^{ik\cdot a_1}⟨R_B|R_B⟩}\\ +\frac{1}{N}{e}^{ik\cdot (a_1+a_2)}⟨R_B|R_B⟩),\end{array}$

Similarly, the second diagonal matrix element $H_{BB}=⟨k_B\left|\widehat{H}\right|k_B⟩$

can be evaluated, by first considering

$\begin{array}{c}\widehat{H}|k_B⟩={\displaystyle \sum R\left({\epsilon }_B|R_B⟩\frac{1}{\sqrt{N}}⟨R_B|R_B⟩\right)}\\ -t_{AB}{\displaystyle \sum R\left(|R_A⟩\frac{1}{\sqrt{N}}⟨R_B|R_B⟩\right)}\\ -t_{AB}{\displaystyle \sum R\left(|R_A⟩\frac{1}{\sqrt{N}}{e}^{-ik\cdot a_1}⟨R_B-a_1|R_B-a_1⟩\right)}\\ -t_{AB}{\displaystyle \sum R(|R_A⟩\frac{1}{\sqrt{N}}{e}^{-ik\cdot (a_1+a_2)}⟨R_B-a_1-a_2|R_B}\\ -a_1-a_2⟩),\end{array}$

The element $H_{BB}=⟨k_B\left|\widehat{H}\right|k_B⟩$

is

By using Equations 3.30 and 3.35, the off-diagonal element $H_{AB}=⟨k_A\left|\widehat{H}\right|k_B⟩$

is given by

$\begin{array}{l}⟨k_A\left|\widehat{H}\right|k_B⟩=-t_{AB}{\displaystyle \sum R(\frac{1}{N}⟨R_A|R_A⟩+\frac{1}{N}{e}^{-ik\cdot a_1}⟨R_A|R_A⟩}\\ +\frac{1}{N}{e}^{-ik\cdot (a_1+a_2)}⟨R_A|R_A⟩),\end{array}$

Therefore, the Hamiltonian matrix is

In this matrix, one can readily understand that the diagonal elements H_AA and H_BB are the energies at the atomic site called “onsite energy” or the energy offset for an atom at sites A and B. The H_AB is the interaction from atom at B sublattices to the atom at sublattice A, or vice versa for H_BA. Notice that, the Hamiltonian matrix is a Hermitian matrix that has two real eigenvalues for each value of k. The solutions to this Hamiltonian can be obtained by diagonalizing the Hamiltonian or from the secular equation

which is

having solutions

or

where the multiplication

With $k=k_x\widehat{x}+k_y\widehat{y}$

and the lattice vectors defined by Equation 3.21, we can evaluate product of

For simplicity, we can define

For a simple nearest-neighbor tight-binding model, one obtains the dispersion relation (Reich et al. 2002; Wallace 1947):

If the atom at sublattice A is identical to the atom at sublattice B (∊_A = ∊_B ≡ ∊ and t_AB ≡ t), then the above equation will get reduced to

Here, the atoms at sites A and B are carbon atoms; we can choose the energy offset ∊ = 0. The values of ∊ and t are found by fitting experimental data or first-principles calculations. In practice, ∊ and t are adjusted to reproduce a good description of π bands at the K point. Based on the first-principles calculations, a typical value of ∊ = 0 eV and t = 2.7 eV are used (Reich et al. 2002). With this choice of reciprocal lattice vectors in Equation 3.22, we can plot the electronic structures of graphene, Equation 3.47, as illustrated in Figure 3.2 within the first BZ as shown in the inset. The band structure of graphene is plotted along Γ − K − M − Γ directions. The lines connecting the high symmetry points, Γ − K − M − Γ, are defined by

where 0 ≤ k ≤ 1 and the unit along these paths is in the unit of 2π/a. In the case that the atom at sublattice A is equivalent to the atom at sublattice B, there will be doubly degenerate energy levels at the K point or at the six corners of the hexagonal BZ, so-called Dirac point (Zhou et al. 2006). Among the six Dirac points, only two are independent denoted by K and K′ in the inset of Figure 3.2, where the others are related to K and K′ by a combination of the reciprocal lattice vectors in Equation 3.22. This characteristic is strongly related to the inversion symmetry of crystal structure, where the inversion center is located between the two sublattices, such that the doubly degenerate levels at the K point will always be presented in other atomic species as well, for example, silicene (Şahin et al. 2009). If the atom at sublattice A is different from the atom at sublattice B, the degeneracy is lifted as in the case of hexagonal form of boron nitride (h-BN) (Pease 1952). The band structures of h-BN can be calculated from Equation 3.46 with the onsite energies of boron ∊_B = 3.6 eV and nitrogen ∊_N = −1.0 eV, and the hopping integral t_BN = 2.5 eV, where the values of all parameters are found by fitting first-principles calculations (Jungthawan et al. 2011). The difference between the two sublattices creates a direct bandgap instead of doubly degenerate levels at K point in h-BN. The bandgap of graphene can be created by several approaches (Peng and Ahuja 2008), one of which is to break the equivalence of A and B sublattices. This simple model clearly shows that the equivalence of A and B sublattices is crucial to electronic structures of honeycomb structure.

A higher-order correction can be made by including the interactions from second nearest neighbors. The Hamiltonian in Equation 3.25 can be written as the Hamiltonian with the interactions from first nearest neighbors ${\widehat{H}}_{1NN}$

, given by Equation 3.26, plus the Hamiltonian with the interactions from second nearest neighbors ${\widehat{H}}_{2NN}$ . From the unit cell in Figure 3.1, sublattice A is surrounded by six neighbors which are located at R_A + a₁, R_A + a₁ + a₂, R_A + a₂, R_A− a₁, R_A − a₁ − a₂, and R_A − a₂. Likewise, the second neighbors of sublattice B are located at R_B + a₁, R_B + a₁ + a₂, R_B + a₂, R_B − a₁, R_B − a₁ − a₂, and R_B − a₂. Therefore, the Hamiltonian up to the second neighbor interactions ${\widehat{H}}_{2NN}$ can be explicitly written as

${t^{\prime }}_A$

and ${t^{\prime }}_B$ are the hopping integrals between corresponding second neighbors and atoms at sites A and B, respectively. The second neighbors Hamiltonian is corresponding to the interactions between the same atomic species with a translation with lattice vector R. This correction will affect only on the diagonal elements, that is, H_AA and H_BB. The off-diagonal elements, that is, H_AB and H_BA, will be the same as the matrix in Equation 3.38. In this case, the first diagonal matrix element $H_{AA}=⟨k_A\left|\widehat{H}\right|k_A⟩$ is given by

Similarly, the second diagonal matrix element $H_{BB}=⟨k_B\left|\widehat{H}\right|k_B⟩$

is given by

With these forms of the diagonal elements, we can construct the Hamiltonian matrix, Equation 3.29, for each value of k,

and obtain the solutions by Equation 3.41. If the atom at sublattice A is identical to the atom at sublattice B(${t^{\prime }}_A={t^{\prime }}_B\equiv t^{\prime }$

), then the solution is

or

where ∊_1NN(k) is the dispersion relation from first neighbors interactions in Equation 3.47. f (k) is a function of k defined by Equation 3.45. According to Reich et al. (2002), the value of t′ for graphene is in the range of −0.2t ≤ t′ ≤ −0.02t depending on the parameterization of ∊ and t. This second neighbor interaction t′ ≠ 0 is responsible for asymmetric feature of the π and π* bands as shown in Figure 3.3.

The Hamiltonian with third nearest-neighbor interactions ${\widehat{H}}_{3NN}$

can be made by considering the unit cell in Figure 3.1. Sublattice A is surrounded by three neighbors which are located at R_B − a₂, R_B + a₂, and R_B − 2a₁ − a₂. Likewise, the third neighbors of sublattice B are located at R_A − a₂, R_A + a₂, and R_A + 2a₁ + a₂. Therefore, the Hamiltonian up to the third neighbors ${\widehat{H}}_{3NN}$ is given by

${t^{″}}_A$

and ${t^{″}}_B$ are the hopping integrals between corresponding third neighbors and atoms at sites A and B, respectively. The third neighbors Hamiltonian is corresponding to the interactions between different atomic species. This correction will affect the off-diagonal elements, that is, H_AB and H_BA, in Equation 3.52. The correction to the third nearest-neighbor interactions is given by the product of $⟨k_A\left|{\widehat{H}}_{3NN}\right|k_B⟩$ and $⟨k_B\left|{\widehat{H}}_{3NN}\right|k_A⟩$ , which are

and

With these corrections to the Hamiltonian matrix elements in Equation 3.52, all the elements in Equation 3.29 are given by

If the atom at sublattice A is identical to the atom at sublattice B, we can let ${t^{″}}_A={t^{″}}_B\equiv t^{″}$

. The eigenvalues of this matrix can be obtained by Equation 3.41 for each value of k. Usually, the second and third neighbor interactions are not significant and generally ignored. However, the change in these interactions can substantially affect the band structure. The value of t′ and t″ can change the number of Dirac points, their position, and their properties in graphene and graphene-related materials (Bena and Simon 2011). It is interesting to note that the simple first nearest-neighbors tight-binding results give a reasonable feature at the band edge near K point. Hence, for the larger system, first nearest-neighbor interactions will be our main focus.