1.1 Graphene Band Structure

Graphene is a two-dimensional structure of carbon atoms. Single-layer graphene consists of only one layer of carbon atoms. A carbon atom has six electrons occupying the 1s², 2s², 2p_x, and 2p_y atomic orbitals. When carbon atoms are brought together, one electron from the 2s orbital is promoted to the 2p_z orbital for the formation of hybrid orbitals. In diamonds, the 2s, 2p_x, 2p_y, and 2p_z orbitals are mixed to form four sp³ hybrid orbitals for each carbon atom; therefore, each carbon atom is joined with four neighbors by overlapping their sp³ hybrid orbitals. In graphite, instead of four sp³ hybrid orbitals, three sp² hybrid orbitals are formed through mixing of the 2s, 2p_x, and 2p_y orbitals, while the fourth orbital remains as 2p_z. Overlapping sp² hybrid orbitals from two adjacent atoms creates a strong σ covalent bond (C‒C bond); these in-plane σ bonds connect each carbon atom to three neighbors. The remaining 2p_z orbitals of these carbon atoms form π bonds, which are responsible for binding carbon layers together in graphite. Because π bonds are much weaker than σ bonds, graphite has a low shear strength so that its carbon layers can be easily detached. For monolayer graphene, these nearly free π electrons are responsible for most of its experimentally observed electronic and optical properties. Because the Pauli exclusion principle requires that π electrons from different carbon atoms do not occupy the same state, the large number of closely packed carbon atoms in graphene causes degenerate energy levels to split into continuously distributed nondegenerate levels of allowed energy states, forming energy bands.

The real-space two-dimensional honeycomb lattice of graphene is shown in Figure 1.1(a). The distance between two neighboring carbon atoms in graphene is

$a\approx 0.142\text{ nm}.$(1.1)

The primitive lattice vectors are

${\mathbf{a}}_1=\sqrt{3}a\left(\frac{\sqrt{3}}{2}\widehat{x}-\frac{1}{2}\widehat{y}\right)\text{and}{\mathbf{a}}_2=\sqrt{3}a\left(\frac{\sqrt{3}}{2}\widehat{x}+\frac{1}{2}\widehat{y}\right).$(1.2)

Therefore, the lattice constant is

$a_0=\left|a_1\right|=\left|a_2\right|=\sqrt{3}a\approx 0.246\text{ nm}.$ (1.3)

The area A_p of a primitive unit cell can be obtained as

$A_{\text{p}}=\left|a_1\times a_2\cdot a_3\right|=\frac{3\sqrt{3}}{2}{a}^2=\frac{\sqrt{3}}{2}{a}_0^2,$(1.4)

where $a_3=a_1\times a_2/A_{\text{p}}=\widehat{z}$ is the unit vector that points in the direction perpendicular to the plane of the two-dimensional graphene lattice. We also find that the vectors connecting adjacent carbon atoms are

${\mathbf{\delta }}_1=a\left(\frac{1}{2}\widehat{x}-\frac{\sqrt{3}}{2}\widehat{y}\right),{\mathbf{\delta }}_2=a\left(\frac{1}{2}\widehat{x}+\frac{\sqrt{3}}{2}\widehat{y}\right)\text{, }\text{and}{\mathbf{\delta }}_3=-a\widehat{x}.$(1.5)

as shown in Figure 1.1(a).

Figure 1.1 (a) Real-space honeycomb graphene lattice. The lattice consists of two overlapping Bravais sublattices, A (gray dots) and B (black dots). The primitive unit cell is drawn as a shaded area. $a_1$ and $a_2$ are the primitive lattice vectors. ${\mathbf{\delta }}_1$ , ${\mathbf{\delta }}_2$ , and ${\mathbf{\delta }}_3$ are vectors pointing from a B atom to its nearest A atoms. (b) Brillouin zone of graphene drawn as a shaded area. $b_1$ and $b_2$ are the primitive vectors of the reciprocal lattice. Dirac points $\text{Κ}$ (black dots) and ${\rm K}\prime$ (gray dots) are marked.

The corresponding Brillouin zone of the lattice in Figure 1.1(a) is shown in Figure 1.1(b). The primitive vectors of the reciprocal lattice are

$b_1=2\pi \frac{a_2\times a_3}{a_1\times a_2\cdot a_3}=\frac{4\pi }{3a}\left(\frac{1}{2}\widehat{x}-\frac{\sqrt{3}}{2}\widehat{y}\right),b_2=2\pi \frac{a_3\times a_1}{a_1\times a_2\cdot a_3}=\frac{4\pi }{3a}\left(\frac{1}{2}\widehat{x}+\frac{\sqrt{3}}{2}\widehat{y}\right).$(1.6)

We can also find the vector $b_3=2\pi a_1\times a_2/\left(a_1\times a_2\cdot a_3\right)=2\pi \widehat{z}$ that is perpendicular to the plane of the two-dimensional reciprocal lattice. Therefore,

$a_i\cdot b_j=2\pi {\delta }_{ij},$(1.7)

where ${\delta }_{ij}$ is the Kronecker delta function: ${\delta }_{ij}=1$ if $i=j$ , and ${\delta }_{ij}=0$ if $i\ne j$ .

The honeycomb lattice of graphene is not a Bravais lattice because the locations of all carbon atoms cannot be generated by the translation

$R=m{a}_1+n{a}_2,$(1.8)

where m and n are integers. The lattice structure of graphene consists of two overlapping Bravais sublattices A and B, and thus a primitive unit cell contains two carbon atoms, as shown in Figure 1.1(a). All of the carbon atoms on the same sublattice, but not those on the two different sublattices, are connected by the vector $R$ .

The electronic band structure of graphene can be obtained by using the tight-binding model. We start from a tight-binding Hamiltonian $\widehat{H}$ for the Schrödinger equation:

$\widehat{H}\psi \left(r\right)=E\psi \left(r\right).$(1.9)

The wave function $\psi \left(r\right)$ is given by the linear superposition of orbital functions ${\varphi }_l\left(r\right)$ :

$\psi \left(r\right)=\frac{1}{\sqrt{N}}{\displaystyle \sum _l{c}_l{\text{e}}^{\text{i}k\cdot R_l}{\varphi }_l\left(r\right)}=\frac{1}{\sqrt{N}}{\displaystyle \sum _l{c}_l{\text{e}}^{\text{i}k\cdot R_l}\varphi \left(r-R_l\right)},$(1.10)

where N is the number of unit cells, ${\varphi }_l\left(r\right)=\varphi \left(r-R_l\right)$ is the ground state of the $2p_z$ electron of an isolated carbon atom that is located at $R_l$ , and the index l runs over all carbon atom points on the graphene lattice. Equation (1.10) is of the Bloch form because $\psi \left(r+a_i\right)=\text{exp}\left(\text{i}k\cdot a_i\right)\psi \left(r\right)$ , where $i=1$ or 2. From the symmetry point of view, all atoms on sublattice A are geometrically identical; in other words, the surrounding is the same when viewing the graphene lattice from any atom on sublattice A. The same can be said for all atoms on sublattice B. No difference can be seen when viewing the graphene lattice from different atoms on the same sublattice. However, the atoms on sublattice A are not geometrically identical to those on sublattice B. By comparing the scenery in a certain direction, it is possible to tell the difference between viewing the graphene lattice from one atom on sublattice A and viewing it from one on sublattice B. Therefore, the wave function in (1.10) has only two independent coefficients for $c_l$ : $c_{\text{A}}$ and $c_{\text{B}}$ , which respectively represent the amplitudes of the wave functions at carbon sites on sublattices A and B. Thus, (1.10) can be rewritten as a linear superposition of two Bloch functions ${\psi }_{\text{A}}$ and ${\psi }_{\text{B}}$ that are respectively wave functions for sublattices A and B:

$\psi \left(r\right)=c_{\text{A}}{\psi }_{\text{A}}\left(r\right)+c_{\text{B}}{\psi }_{\text{B}}\left(r\right)=\frac{c_{\text{A}}}{\sqrt{N}}{\displaystyle \sum _{l_{\text{A}}}{\text{e}}^{\text{i}k\cdot R_{l_{\text{A}}}}{\varphi }_{l_{\text{A}}}\left(\mathbf{\text{r}}\right)}+\frac{c_{\text{B}}}{\sqrt{N}}{\displaystyle \sum _{l_{\text{B}}}{\text{e}}^{\text{i}k\cdot R_{l_{\text{B}}}}{\varphi }_{l_{\text{B}}}\left(\mathbf{\text{r}}\right)},$(1.11)

where $l_{\text{A}}$ and $l_{\text{B}}$ run over all of the carbon atoms on sublattices A and B, respectively. The wave functions ${\psi }_{\text{A}}$ and ${\psi }_{\text{B}}$ are normalized, and the overlap of ${\psi }_{\text{A}}$ and ${\psi }_{\text{B}}$ is negligible [Reference Reich, Maultzsch, Thomsen and Ordejón1], so that $〈{\psi }_i|{\psi }_j〉={\delta }_{ij}$ . Therefore, $|{\psi }_{\text{A}}〉$ and $|{\psi }_{\text{B}}〉$ form the basis for the graphene wave functions.

By multiplying both sides of (1.9) by ${\psi }_{\text{A}}^{\ast }\left(r\right)$ , we obtain

$c_{\text{A}}〈{\psi }_{\text{A}}\left|\widehat{H}\right|{\psi }_{\text{A}}〉+c_{\text{B}}〈{\psi }_{\text{A}}\left|\widehat{H}\right|{\psi }_{\text{B}}〉=c_{\text{A}}E.$(1.12)

In (1.12), all of the interactions among A atoms and those between A atoms and B atoms are considered. However, it is very difficult to solve for the eigenenergy from (1.12) where all interactions are accounted for. To avoid such difficulty, (1.12) is simplified by taking the approximation of considering only the self-interaction of each carbon atom at a lattice point and the interactions of the atom with its three nearest neighbors; all long-range interactions are assumed to be much weaker than these interactions and are thus ignored. Within this approximation, (1.12) takes the simple form:

$c_{\text{A}}〈{\varphi }_{l_{\text{A}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{A}}}〉+c_{\text{B}}〈{\varphi }_{l_{\text{A}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{B}}}〉\left({\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_1}+{\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_2}+{\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_3}\right)=c_{\text{A}}E.$(1.13)

By multiplying both sides of (1.11) by ${\psi }_{\text{B}}^{\ast }\left(r\right)$ and taking the same approximation, we have

$c_{\text{B}}〈{\varphi }_{l_{\text{B}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{B}}}〉+c_{\text{A}}〈{\varphi }_{l_{\text{B}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{A}}}〉\left({\text{e}}^{\text{i}k\cdot {\mathbf{\delta }}_1}+{\text{e}}^{\text{i}k\cdot {\mathbf{\delta }}_2}+{\text{e}}^{\text{i}k\cdot {\mathbf{\delta }}_3}\right)=c_{\text{B}}E.$(1.14)

Equations (1.13) and (1.14) can be expressed in a matrix form for the eigenenergy equation of the system:

$\left[\begin{array}{cc}\epsilon & -{\gamma }_0\left({\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_1}+{\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_2}+{\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_3}\right)\\ -{\gamma }_0{\left({\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_1}+{\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_2}+{\text{e}}^{-\text{i}k\cdot {\mathbf{\delta }}_3}\right)}^{\ast }& \epsilon \end{array}\right]\left[\begin{array}{c}{c}_{\text{A}}\\ c_{\text{B}}\end{array}\right]=E\left[\begin{array}{c}{c}_{\text{A}}\\ c_{\text{B}}\end{array}\right],$(1.15)

where

$\epsilon =〈{\varphi }_{l_{\text{A}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{A}}}〉=〈{\varphi }_{l_{\text{B}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{B}}}〉$(1.16)

is the self-interaction energy, and

${\gamma }_0=-〈{\varphi }_{l_{\text{A}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{B}}}〉=-〈{\varphi }_{l_{\text{B}}}\left|\widehat{H}\right|{\varphi }_{l_{\text{A}}}〉$(1.17)

is the nearest-neighbor hopping energy in graphene. The nearest-neighbor hopping energy ${\gamma }_0$ is experimentally measured to have a value of ${\gamma }_0=3.16\text{ eV}$ [Reference Peres2–Reference Abergel, Apalkov, Berashevich, Ziegler and Chakraborty4], which is the same as that found for graphite.

The nontrivial eigenvalues of (1.15) yield two eigenenergies:

$E\left(k\right)=\epsilon \pm {\gamma }_0\left|{\text{e}}^{\text{i}k\cdot {\mathbf{\delta }}_1}+{\text{e}}^{\text{i}k\cdot {\mathbf{\delta }}_2}+{\text{e}}^{\text{i}k\cdot {\mathbf{\delta }}_3}\right|,$(1.18)

where the plus sign is for the conduction band, which has higher energy, and the minus sign is for the valence band, which has lower energy. The conduction band and the valence band described by (1.18) are also called the ${\mathit{\pi }}^{\ast }$ band and the $\pi$ band, respectively. Equation (1.18) can be written in another form that is frequently seen in the literature [Reference Wallace5,Reference Castro Neto, Guinea, Peres, Novoselov and Geim6]:

$E\left(k\right)=\epsilon \pm {\gamma }_0\sqrt{3+2\text{cos}\left(\sqrt{3}{k}_{y}a\right)+4\text{cos}\left(\frac{\sqrt{3}}{2}{k}_{y}a\right)\text{cos}\left(\frac{3}{2}{k}_{x}a\right)}.$(1.19)

The energy band structure described by (1.19) is plotted in Figure 1.2(a). At the Г point ($k=0$ ), we find the eigenenergies $E=\epsilon \pm 3{\gamma }_0$ from (1.19), which are associated with the eigenstates $\left[c_{\text{A}},c_{\text{B}}\right]=\frac{1}{\sqrt{2}}\left[1,\mp 1\right]$ , respectively. The eigenenergy $E=\epsilon -3{\gamma }_0$ for the valence band is associated with symmetric bonding orbitals such that the amplitudes $c_{\text{A}}$ and $c_{\text{B}}$ of the wave functions for sublattices A and B have the same sign. By contrast, the eigenenergy $E=\epsilon +3{\gamma }_0$ for the conduction band is associated with antibonding orbitals such that the amplitudes $c_{\text{A}}$ and $c_{\text{B}}$ of the wave functions for sublattices A and B have opposite signs. An eigenstate of graphene is not necessarily bonding or antibonding, however; most eigenstates are a mixture of both bonding and antibonding orbitals. As we shall see in the following, only along a specific set of directions can the eigenstates be described by purely bonding or antibonding orbitals. The conduction band (${\text{π}}^{\ast }$ band) and the valence band ($\text{π}$ band) are symmetric with respect to the $E=0$ plane, as shown in Figure 1.2(a). Note that the long-range interactions of carbon atoms beyond the nearest neighbors are ignored in the above analysis. When these long-range interactions are accounted for, we find that the ${\text{π}}^{\ast }$ and $\text{π}$ bands are actually not symmetric to each other, as seen in Figure 1.2(c), where we also show the ${\sigma }^{\ast }$ and $\sigma$ bands contributed by the $\sigma$ bonds.

Figure 1.2 (a) Band structure of graphene plotted using (1.19) while setting the self-interaction energy to be $\epsilon =0$ . The top surface of a higher energy is the conduction band (${\text{π}}^{\ast }$ band), and the bottom surface of a lower energy is the valence band ($\text{π}$ band). The conduction band is projected onto a two-dimensional plane on the top. The circled region near the Dirac point ${\rm K}\prime$ is enlarged in (b). (c) Full band structure of graphene. The thick curves are the ${\text{π}}^{\ast }$ and $\text{π}$ bands, and the thin curves are the ${\sigma }^{\ast }$ and $\sigma$ bands [Reference Sedelnikova, Bulusheva and Okotrub7].

As can be seen in Figure 1.2(a), the conduction band touches the valence band at $E=\epsilon$ . The self-interaction energy $\epsilon$ can be taken as the reference level by setting it to be 0 for simplicity so that an electron in a state on the conduction band has a positive energy and that in a state on the valence band has a negative energy. Then, (1.19) can be simplified as

$E\left(k\right)=\pm {\gamma }_0\sqrt{3+2\text{cos}\left(\sqrt{3}{k}_{y}a\right)+4\text{cos}\left(\frac{\sqrt{3}}{2}{k}_{y}a\right)\text{cos}\left(\frac{3}{2}{k}_{x}a\right)},$(1.20)

where the plus sign represents the conduction band and the minus sign represents the valence band. The points at $E\left(k\right)=0$ are called the Dirac points; the vicinities of these points are referred to as the valleys, where the graphene band structure can be described by the relativistic Dirac equation, which is further discussed in the following section. Corresponding to the two distinct sublattices A and B in the real space, there are two distinct groups of Dirac points $\text{K}$ and ${\rm K}\prime$ in the k space. Inside the Brillouin zone drawn in Figure 1.1(b), two Dirac points can be found at the k-space locations represented by the reciprocal space vectors:

$K=\frac{4\pi }{3\sqrt{3}a}\left(\frac{\sqrt{3}}{2}\widehat{x}+\frac{1}{2}\widehat{y}\right),K\prime =\frac{4\pi }{3\sqrt{3}a}\left(\frac{\sqrt{3}}{2}\widehat{x}-\frac{1}{2}\widehat{y}\right).$(1.21)

Another popular choice of the k-space locations of $\text{K}$ and ${\rm K}\prime$ in the literature is to select the $\text{K}$ and ${\rm K}\prime$ points on the k_y axis so that they are mirror symmetric with respect to the k_x axis; the physics and the results obtained in the following are unchanged by this alternative selection. It can be easily shown using (1.20) that $E\left(K\right)=E\left(K\prime \right)=0$ at these Dirac points.

For perfectly intrinsic graphene, the chemical potential, i.e., the Fermi level, is located at the energy level of the Dirac points, $E=0$ , so that the valence band is completely filled and the conduction band is completely empty. As most physics and carrier transitions happen at energy levels around the chemical potential, we shall rewrite (1.20) by shifting the origin of the k space to one of the Dirac points and assume a small k to simplify (1.20). Substituting $k$ in (1.20) with $K+k$ or $K\prime +k$ , we obtain, for $\left|k\right|a\ll 1$ in the vicinity of a Dirac point,

$E=\pm {\gamma }_0\sqrt{\frac{9}{4}{k}_y^2{a}^2+\frac{9}{4}{k}_x^2{a}^2}=\pm \frac{3}{2}a{\gamma }_0\left|k\right|=\pm \hslash {\upsilon }_{\text{F}}k,$(1.22)

where $\hslash$ is the reduced Planck’s constant and

${\upsilon }_{\text{F}}=\frac{3a{\gamma }_0}{2\hslash }\approx 1.0\times {10}^6\text{m}{\text{s}}^{-1}$(1.23)

is the Fermi velocity in graphene, the physical meaning of which will become clear in Section 1.5. Equation (1.22) indicates that there is a region in the k space around each Dirac point where the carrier energy is linearly proportional to the wave number measured with respect to the given Dirac point. This region is called the Dirac cone, as shown in the insert in Figure 1.2(b). Because (1.22) can also be derived from the massless Dirac equation, electrons on the Dirac cone are also called Dirac electrons.

The Hamiltonian around the $\text{K}$ point is obtained in the same way as (1.22) is obtained. By substituting $k$ with $K+k$ in (1.15) to set the origin of the k space at the $\text{K}$ point, the off-diagonal matrix elements for $\left|k\right|a\ll 1$ become

$〈{\psi }_{\text{A}}\left|{\widehat{H}}_{\text{K}}\right|{\psi }_{\text{B}}〉=-{\gamma }_0\left({\text{e}}^{-\text{i}\left(K+k\right)\cdot {\mathbf{\delta }}_1}+{\text{e}}^{-\text{i}\left(K+k\right)\cdot {\mathbf{\delta }}_2}+{\text{e}}^{-\text{i}\left(K+k\right)\cdot {\mathbf{\delta }}_3}\right)=\frac{\hslash {\upsilon }_{\text{F}}}{2}\left[\sqrt{3}{k}_x+k_y+\text{i}\left(k_x-\sqrt{3}{k}_y\right)\right]$(1.24)

and

$〈{\psi }_{\text{B}}\left|{\widehat{H}}_{\text{K}}\right|{\psi }_{\text{A}}〉={〈{\psi }_{\text{A}}\left|{\widehat{H}}_{\text{K}}\right|{\psi }_{\text{B}}〉}^{\ast }.$(1.25)

The coefficients of $k_x$ and $k_y$ can be simplified by rotating both the $k_x$ and $k_y$ coordinate axes counterclockwise by 30 degrees through the unitary transformation

$\left[\begin{array}{l}k{\prime }_x\\ k{\prime }_y\end{array}\right]=\left[\begin{array}{cc}\text{cos}\frac{\pi }{6}& \text{sin}\frac{\pi }{6}\\ -\text{sin}\frac{\pi }{6}& \text{cos}\frac{\pi }{6}\end{array}\right]\left[\begin{array}{l}{k}_x\\ k_y\end{array}\right]=\frac{1}{2}\left[\begin{array}{l}\sqrt{3}{k}_x+k_y\\ -k_x+\sqrt{3}{k}_y\end{array}\right].$(1.26)

Using (1.26), we obtain the relation $〈{\psi }_{\text{A}}\left|{\widehat{H}}_{\text{K}}\right|{\psi }_{\text{B}}〉=\hslash {\upsilon }_{\text{F}}\left(k{\prime }_x-\text{i}k{\prime }_y\right)={〈{\psi }_{\text{B}}\left|\widehat{H}\text{K}\right|{\psi }_{\text{A}}〉}^{\ast }$ from (1.24) and (1.25). By dropping the prime symbols, $k{\prime }_x\to k_x$ and $k{\prime }_y\to k_y$ , to simplify the notation, the Hamiltonian near the $\text{K}$ point for $\left|k\right|a\ll 1$ can be expressed as

${\widehat{H}}_{\text{K}}=\left[\begin{array}{cc}0& \hslash {\upsilon }_{\text{F}}\left(k_x-\text{i}{k}_y\right)\\ \hslash {\upsilon }_{\text{F}}\left(k_x+\text{i}{k}_y\right)& 0\end{array}\right].$(1.27)

The corresponding eigenfunctions are

${\psi }_{k,\pm }^{\text{K}}\left(r\right)=\frac{{\text{e}}^{\text{i}k\cdot r}}{\sqrt{2}}\left[\begin{array}{c}1\\ \pm \frac{k_x+\text{i}{k}_y}{k}\end{array}\right]=\frac{{\text{e}}^{\text{i}k\cdot r}}{\sqrt{2}}\left[\begin{array}{c}1\\ \pm {\text{e}}^{\text{i}{\theta }_k}\end{array}\right],$(1.28)

where the $\pm$ signs correspond to the signs of the eigenenergies given in (1.22), and ${\theta }_k$ is the polar angle measured between $k$ and the rotated $k_x$ axis as shown in Figure 1.3(a) for the $\text{K}$ point.

Figure 1.3 Wave vector $k$ near (a) the $\text{K}$ point, and (b) the ${\rm K}\prime$ point. The Hamiltonian near the $\text{K}$ or ${\rm K}\prime$ point is obtained by rotating the coordinate from the gray axes to the black axes. The angle ${\theta }_k$ is measured with respect to the rotated $k_x$ axis.

Similarly, the Hamiltonian near the ${\rm K}\prime$ point is obtained by substituting $k$ with $K\prime +k$ in (1.15) and taking the unitary transformation by rotating the axes 30 degrees clockwise; the result for $\left|k\right|a\ll 1$ is

${\widehat{H}}_{{\text{K}}^{\prime }}=\left[\begin{array}{cc}0& \hslash {\upsilon }_{\text{F}}\left(k_x+\text{i}{k}_y\right)\\ \hslash {\upsilon }_{\text{F}}\left(k_x-\text{i}{k}_y\right)& 0\end{array}\right]={\widehat{H}}_{\text{K}}^{\text{T}}.$(1.29)

The eigenfunctions are

${\psi }_{k,\pm }^{{\text{K}}^{\prime }}\left(r\right)=\frac{{\text{e}}^{\text{i}k\cdot r}}{\sqrt{2}}\left[\begin{array}{c}\pm {\text{e}}^{\text{i}{\theta }_k}\\ 1\end{array}\right],$(1.30)

where the $\pm$ signs correspond to the signs of the eigenenergies given in (1.22), and ${\theta }_k$ is the polar angle measured between $k$ and the rotated $k_x$ axis at the $\text{K}\prime$ point, as shown in Figure 1.3(b). From (1.28) and (1.30), it is clear that for the conduction band, the eigenstate is antibonding at ${\theta }_k=\pi$ and bonding at ${\theta }_k=0$ , as shown in Figure 1.4(a). By contrast, the eigenstate is bonding at ${\theta }_k=\pi$ and antibonding at ${\theta }_k=0$ for the valence band, as shown in Figure 1.4(b).

Equations (1.27) and (1.29) are often written as

${\widehat{H}}_{\text{K}}={\upsilon }_{\text{F}}\mathbf{\sigma }\mathbf{\cdot }p$(1.31)

for an electron near a $\text{K}$ point and as

${\widehat{H}}_{{\text{K}}^{\prime }}={\upsilon }_{\text{F}}{\mathbf{\sigma }}^{\ast }\cdot p$(1.32)

for an electron near a $\text{K}\prime$ point, where $p=\hslash k$ and $\mathbf{\sigma }$ is the Pauli vector:

$\mathbf{\sigma }=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\widehat{x}+\left[\begin{array}{cc}0& -\text{i}\\ \text{i}& 0\end{array}\right]\widehat{y}+\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\widehat{z}.$(1.33)

Note that $p=\hslash k=\hslash \left(k_x\widehat{x}+k_y\widehat{y}\right)$ with $k_z=0$ for the $\text{π}$ electrons in monolayer graphene because they only move through the graphene lattice on the xy plane of a monolayer graphene sheet.

Figure 1.4 Antibonding (solid thin lines) and bonding (dotted lines) states in (a) the conduction band and (b) the valence band.

In fact, the Hamiltonian of graphene near a $\text{K}$ or ${\rm K}\prime$ point has the form identical to that of a spin one-half particle in a magnetic field; its eigenfunctions are two-component spinors, as seen in (1.28) and (1.30). In the case of graphene, $\mathbf{\sigma }$ does not really represent an electronic spin, but a pseudospin that indicates the sublattice on which the electron is located. We can define a helicity operator $\widehat{h}=\mathbf{\sigma }\cdot p/p$ to find the projection of the momentum along the pseudospin direction; then it can be shown that $\widehat{h}{\psi }_{k,\pm }^{\text{K}}=\pm {\psi }_{k,\pm }^{\text{K}}$ . The positive sign gives positive helicity, meaning that the momentum of an electron on the conduction band is parallel to the pseudospin, whereas the negative sign gives negative helicity, meaning that the momentum of an electron on the valence band is antiparallel to the pseudospin. The relation among the helicity, the carrier momentum, and the pseudospin are shown in Figure 1.5. Because of the conservation of pseudospin, a right-moving electron on the conduction band cannot be scattered into a left-moving electron on the conduction band or into a right-moving electron on the valence band. Such scattering requires the pseudospin to be flipped from $+1$ to $-1$ for an electron near the K point, or flipped from $-1$ to $+1$ for an electron near the ${\rm K}\prime$ point, neither of which is allowed. Similarly, scattering of a left-moving electron on the conduction band into a right-moving electron on the conduction band or into a left-moving electron on the valence band is also forbidden. Nevertheless, such scattering events can happen when there exists a short-range potential that acts differently on sublattices A and B, thus breaking the symmetry between the two sublattices.

Figure 1.5 Helicity h, pseudospin $\mathbf{\sigma }$ , and momentum $p$ of an electron near (a) the $\text{K}$ point, and (b) the ${\rm K}\prime$ point. The directions of $p$ and $\mathbf{\sigma }$ are shown. In the case of a hole, the direction of momentum is reversed. Due to the conservation of pseudospin, an electron or a hole cannot be scattered into an electron or a hole of a different pseudospin.

1.2 Density of States and Carrier Concentration

With the band structure near the Dirac point described by (1.22), the density of electron states of graphene in the energy range between $E$ and $E+\text{d}E$ for $E>0$ near the conduction band edge is

${\tilde{\rho }}_{\text{c}}\left(E\right)\text{d}E=\frac{1}{A}{\displaystyle \sum _{k,E_k>0}\delta \left(E_k-E\right)\text{d}E}=\frac{g}{4{\pi }^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}\delta \left(\hslash {\upsilon }_{\text{F}}k-E\right)2\pi k\text{d}k}\text{d}E=\frac{2E}{\pi {\hslash }^2{\upsilon }_{\text{F}}^2}\text{d}E,$(1.34)

and that for $E<0$ near the valence band edge is

${\tilde{\rho }}_{\text{v}}\left(E\right)\text{d}E=\frac{1}{A}{\displaystyle \sum _{k\text{, }{E}_k<0}\delta \left(E_k-E\right)\text{d}E}=\frac{g}{4{\pi }^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}\delta \left(-\hslash {\upsilon }_{\text{F}}k-E\right)2\pi k\text{d}k}\text{d}E=-\frac{2E}{\pi {\hslash }^2{\upsilon }_{\text{F}}^2}\text{d}E,$(1.35)

where $A$ is the area of the graphene sheet and $g=4$ is the total degeneracy due to the spin degeneracy and the valley degeneracy (two valleys $\text{K}$ and ${\rm K}\prime$ in one Brillouin zone). Because of the symmetry between the conduction band and the valence band with respect to the Dirac point, ${\tilde{\rho }}_{\text{c}}\left(E\right)\text{d}E={\tilde{\rho }}_{\text{v}}\left(-E\right)\text{d}E$ for the same absolute value $\left|E\right|$ of energy.

From (1.34) and (1.35), the surface concentrations of electrons and holes of a graphene sheet are, respectively,

${\tilde{n}}_0={\displaystyle \underset{0}{\overset{\infty }{\int }}{\tilde{\rho }}_{\text{c}}\left(E\right)f_0\left(E\right)\text{d}E},{\tilde{p}}_0={\displaystyle \underset{-\infty }{\overset{0}{\int }}{\tilde{\rho }}_{\text{v}}\left(E\right)[1-f_0\left(E\right)]\text{d}E},$(1.36)

where

$f_0\left(E\right)=\frac{1}{{\text{e}}^{\left(E-\mu \right)/k_{\text{B}}T}+1}$(1.37)

is the equilibrium Fermi–Dirac distribution function for a chemical potential $\mu$ at a temperature T, and $k_{\text{B}}$ is the Boltzmann constant. The probability of finding an electron in a state at the energy level E is $f_0\left(E\right)$ , and the probability of finding a hole in a state at the energy level E is $1-f_0\left(E\right)$ . When the condition $\left|\mu \right|\gg k_{\text{B}}T$ is satisfied, $f_0\left(E\right)$ is approximately a Heaviside step function $H\left(\mu -E\right)$ , which has a value of 1 for $E\le \mu$ and 0 for $E>\mu$ . Therefore, for $\left|\mu \right|\gg k_{\text{B}}T$ ,

${\tilde{n}}_0+{\tilde{p}}_0\approx {\displaystyle \underset{0}{\overset{\infty }{\int }}{\tilde{\rho }}_{\text{c}}\left(E\right)H\left(\mu -E\right)\text{d}E}+{\displaystyle \underset{-\infty }{\overset{0}{\int }}{\tilde{\rho }}_{\text{v}}\left(E\right)[1-H\left(\mu -E\right)]\text{d}E}=\frac{{\mu }^2}{\pi {\hslash }^2{\upsilon }_{\text{F}}^2}.$(1.38)

In this book, we use a tilde for a 2D quantity to clearly distinguish it from its 3D counterpart in order to avoid confusion. Therefore, the surface electron density $\tilde{n}$ , the surface hole density $\tilde{p}$ , and the densities of states ${\tilde{\rho }}_{\text{c}}\left(E\right)\text{d}E$ and ${\tilde{\rho }}_{\text{v}}\left(E\right)\text{d}E$ of graphene are all numbers per unit area in the unit of ${\text{m}}^{-2}$ , or ${\text{cm}}^{-2}$ ; whereas their 3D counterparts, $n$ , $p$ , ${\rho }_{\text{c}}\left(E\right)\text{d}E$ , and ${\rho }_{\text{v}}\left(E\right)\text{d}E$ , are all numbers per unit volume in the unit of ${\text{m}}^{-3}$ , or ${\text{cm}}^{-3}$ .

1.3 Fermi Energy, Chemical Potential, and Fermi Level

Here we clarify the terminologies and notations that are sometimes swapped or not clear in some books and literature: the Fermi energy $E_{\text{F}}$ , the chemical potential $\mu$ , and the Fermi level. The Fermi energy $E_{\text{F}}$ is defined as the energy difference between the highest energy level and the lowest energy level that are occupied by electrons at the absolute zero temperature, $T=0$ . Therefore, the Fermi energy is defined only at $T=0$ ; it is theoretically defined but experimentally cannot be exactly measured because the absolute temperature can never truly reach zero. For a metal, the Fermi energy is usually taken to be the kinetic energy of the electron that occupies the highest energy level at $T=0$ ; therefore, it is always positive. For graphene, we take the Fermi energy to be the total energy of the electron that occupies the highest energy level at $T=0$ , which can be either positive or negative, or zero, if the reference point of zero energy is defined at the Dirac point, as seen below. By contrast, the chemical potential is meaningful at any temperature because it is a statistical reference for the total energy, including the potential energy and the kinetic energy, of a particle in a system in thermal equilibrium. In solid-state physics, the chemical potential $\mu$ usually refers to the change of free energy when an additional electron is added to a system. The effect of $\mu$ is best observed in the Fermi–Dirac distribution $f_0\left(E\right)$ , where one can see that at zero temperature, the minimum energy required to add an additional electron is $\mu$ because below it electronic states are all occupied. In the electronic devices context, the term Fermi level usually refers to the chemical potential $\mu$ . However, in semiconductor physics, the Fermi level is usually denoted by the symbol $E_{\text{F}}$ though it refers to the chemical potential $\mu$ . To avoid confusion with the Fermi energy, in this book we use the symbol $E_{\text{F}}$ to represent only the Fermi energy while using the symbol $\mu$ to represent the chemical potential and the Fermi level.

From the above discussions, we see that $\mu$ is a function of temperature: $\mu \left(T\right)$ . In the case when $\mu \left(0\right)$ lies within the conduction band, $E_{\text{F}}$ is positive and can be found from the definition $E_{\text{F}}\equiv \mu \left(0\right)$ at zero temperature. In the case when $\mu \left(0\right)$ lies within the valence band, we define $E_{\text{F}}$ to be negative using the definition $E_{\text{F}}\equiv \mu \left(0\right)$ at zero temperature. The temperature dependence of $\mu \left(T\right)$ of graphene can be found by recognizing the fact that each additional thermally excited electron always comes with the creation of a hole to result in an electron–hole pair. Therefore, for a given piece of graphene, the difference between the electron concentration and the hole concentration, ${\tilde{n}}_0-{\tilde{p}}_0$ , is a conserved quantity that is independent of temperature. In the case when $\mu \left(0\right)=E_{\text{F}}>0$ at zero temperature so that the chemical potential lies within the conduction band, we have ${\tilde{p}}_0=0$ and ${\tilde{n}}_0-{\tilde{p}}_0={\tilde{n}}_0+{\tilde{p}}_0=E_{\text{F}}^2/\pi {\hslash }^2{\upsilon }_{\text{F}}^2$ from (1.38). Therefore we can relate $\mu$ and $E_{\text{F}}$ by combining (1.36) and (1.38) to write

$\frac{E_{\text{F}}^2}{\pi {\hslash }^2{\upsilon }_{\text{F}}^2}={\displaystyle \underset{0}{\overset{\infty }{\int }}{\tilde{\rho }}_{\text{c}}\left(E\right)f_0\left(E\right)\text{d}E}-{\displaystyle \underset{-\infty }{\overset{0}{\int }}{\tilde{\rho }}_{\text{v}}\left(E\right)[1-f_0\left(E\right)]\text{d}E},$(1.39)

which can be written in a compact form:

$\frac{E_{\text{F}}^2}{2}={\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{E\text{sinh}\left(\mu /k_{\text{B}}T\right)}{\text{cosh}\left(E/k_{\text{B}}T\right)+\text{cosh}\left(\mu /k_{\text{B}}T\right)}\text{d}E},$(1.40)

from which $\mu$ can be numerically obtained. In the case when $\mu \left(0\right)$ lies within the valence band so that $E_{\text{F}}\equiv \mu \left(0\right)<0$ , we have ${\tilde{n}}_0=0$ and $-{\tilde{n}}_0+{\tilde{p}}_0={\tilde{n}}_0+{\tilde{p}}_0=E_{\text{F}}^2/\pi {\hslash }^2{\upsilon }_{\text{F}}^2$ from (1.38). Then, following a procedure similar to the above, we obtain

$\frac{E_{\text{F}}^2}{2}={\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{E\text{sinh}\left(-\mu /k_{\text{B}}T\right)}{\text{cosh}\left(E/k_{\text{B}}T\right)+\text{cosh}\left(\mu /k_{\text{B}}T\right)}\text{d}E}.$ (1.41)

When $\left|\mu \right|\gg k_{\text{B}}T$ , the integrals for both (1.40) and (1.41) approach ${\mu }^2/2$ , and we have $\mu \approx E_{\text{F}}$ .

1.4 Temperature Dependence of Carrier Concentration

The dependence of the free carrier concentration ${\tilde{n}}_0$ of graphene on temperature is shown in Figure 1.6 for $E_{\text{F}}=50\text{meV}$ . The characteristics obtained under two different approximations are compared to the exact temperature dependence of the carrier concentration obtained without approximations. The dotted line is the low-temperature approximation given by (1.38), and the dashed curve is obtained from (1.36) by assuming $\mu =E_{\text{F}}$ independent of temperature. The solid curve is the exact result calculated from (1.36), with $\mu$ varying with temperature obtained from (1.40). As can be seen, the dotted line is accurate only when the temperature is low, and the dashed curve overestimates the carrier concentration. In general, a higher temperature leads to a higher carrier concentration because of thermally excited electron–hole pairs.

Figure 1.6 Carrier concentration for $E_{\text{F}}=50\text{meV}$ as a function of temperature calculated under different assumptions. The dotted line and the dashed curve are obtained under the low-temperature approximation and the constant $\mu$ approximation, respectively, while the solid curve is the exact result.

1.5 Carrier Velocity and Effective Mass

In the preceding sections, we regard a charge carrier in graphene as a wave that has a momentum of $\hslash k$ and an energy of $E$ . In this section, we treat the wave packet as a particle that has a velocity equal to the group velocity of its wave packet. For an isotropic two-dimensional material such as graphene, the phase velocity ${\upsilon }_{\text{p}}$ and the group velocity ${\upsilon }_{\text{g}}$ of an electron that has a momentum of $\hslash k$ and an energy of $E$ are, respectively,

${\upsilon }_{\text{p}}=\frac{E}{\hslash k}\text{and}{\upsilon }_{\text{g}}=\frac{\partial E}{\hslash \partial k}.$ (1.42)

For an electron in graphene near a Dirac point, the carrier energy is linearly proportional to the momentum: $E=\hslash {\upsilon }_{\text{F}}k$ ; thus, we find that ${\upsilon }_{\text{p}}={\upsilon }_{\text{g}}={\upsilon }_{\text{F}}\approx 1\times {10}^6{\text{ m}\text{s}}^{-1}$ , independent of the electron energy. The same conclusion applies to a hole in graphene near a Dirac point. Therefore, the parameter ${\upsilon }_{\text{F}}$ represents both the group velocity and the phase velocity of a carrier in graphene. The parameter ${\upsilon }_{\text{F}}$ also represents the Fermi velocity because a carrier that has a kinetic energy equal to the Fermi energy also moves at a group velocity of ${\upsilon }_{\text{g}}={\upsilon }_{\text{F}}$ in the vicinity of a Dirac point. This result is distinctively different from the energy-dependent group velocity ${\upsilon }_{\text{g}}=\sqrt{2E/m^{\ast }}$ of a carrier in an ordinary material, where $m^{\ast }$ is the effective mass of an electron in the conduction band of the material or that of a hole in the valence band.

The linear energy–momentum relation of the carriers in graphene resembles the relativistic Dirac equation. Unlike the Schrödinger equation, the Dirac equation accounts for special relativity through the relativistic energy–momentum relation:

$E^2=p^2{c}^2+m_0^2{c}^4,$(1.43)

where $p=\hslash k$ is the particle momentum, $m_0$ is the rest mass of the particle, and c is the speed of light. If we compare the energy–momentum relation given in (1.22) for a carrier in graphene with the relativistic energy–momentum relation given in (1.43), we find that a carrier in graphene has an effective speed of $c^{\ast }={\upsilon }_{\text{F}}$ and a zero effective rest mass of $m_0=0$ . Therefore, the carriers in graphene are often called massless Dirac fermions.

In contrast to the effective rest mass, an effective mass associated with the acceleration of a carrier in graphene under the influence of an external force can be defined. For example, when a carrier in graphene is accelerated by an external magnetic field, the carrier follows a circular orbit dictated by Maxwell’s equations. If the area enclosed by the orbit is $A\left(E\right)$ , the cyclotron effective mass is defined as [Reference Castro Neto, Guinea, Peres, Novoselov and Geim6,Reference Miller, Kubista and Rutter8,Reference Novoselov, Geim and Morozov9]

$m^{\ast }=\frac{{\hslash }^2}{2\pi }\frac{\partial A\left(E\right)}{\partial E}.$(1.44)

With $E=\hslash {\upsilon }_{\text{F}}k$ and $A\left(E\right)=\pi k^2=\pi E^2/{\left(\hslash {\upsilon }_{\text{F}}\right)}^2$ , we obtain

$m^{\ast }=\frac{E}{{\upsilon }_{\text{F}}^2}=\frac{\hslash k}{{\upsilon }_{\text{F}}}.$(1.45)

Therefore, the cyclotron effective mass is energy dependent, and (1.45) manifests the mass–energy equivalence from special relativity, namely, $E=m^{\ast }{c}^{\ast }{}^2$ with $c^{\ast }={\upsilon }_{\text{F}}$ .

In solid-state physics, the effective mass of a charge carrier is often calculated as $m^{\ast }={\hslash }^2{\left({\partial }^{2}E/\partial k^2\right)}^{-1}$ because the energy of a carrier is quadratically dependent on its momentum near a conduction band minimum or a valence band maximum of an ordinary three-dimensional solid-state material. However, in the case of graphene, this relation is invalid because it gives a mathematically divergent result due to the fact that ${\partial }^{2}E/\partial k^2=0$ for a carrier in graphene. The reason is that the energy of a carrier in graphene near a Dirac point varies with its momentum not quadratically but linearly: $E=\pm \hslash {\upsilon }_{\text{F}}k$ from (1.22). Thus the assumption of a constant effective mass given by $m^{\ast }={\hslash }^2{\left({\partial }^{2}E/\partial k^2\right)}^{-1}$ is not valid for a carrier in graphene. Instead, $m^{\ast }$ for a carrier in graphene can be calculated using the momentum relation $m^{\ast }{\upsilon }_{\text{g}}=m^{\ast }{\upsilon }_{\text{F}}=\hslash k$ so that

$m^{\ast }=\frac{\hslash k}{{\upsilon }_{\text{g}}}=\frac{\hslash k}{{\upsilon }_{\text{F}}}=\frac{E}{{\upsilon }_{\text{F}}^2},$(1.46)

which gives a result that is consistent with (1.45). This relation yields $m^{\ast }=0$ at the Dirac point where $k=0$ ; this result is also conceptually consistent with the zero effective rest mass for a carrier in graphene discussed above.

The comparison of parabolic 2D/3D systems and monolayer graphene for various physical quantities is presented in Table 1.1.

Table 1.1 Physical quantities for parabolic 2D/3D systems and monolayer graphene.^a

Quantity	Parabolic 3D	Parabolic 2D	Graphene	Equation
$E_{\text{F}}=\mu \left(0\right)\ge 0$	$\frac{{\hslash }^2}{2m^{\ast }}{\left(\frac{6{\pi }^{2}n}{g}\right)}^{2/3}$	$\frac{2{\hslash }^2\pi \tilde{n}}{m^{\ast }g}$	$\hslash {\upsilon }_{\text{F}}{\left(\frac{4\pi \tilde{n}}{g}\right)}^{1/2}$	(1.38)
$k_{\text{F}}$	${\left(\frac{6{\pi }^{2}n}{g}\right)}^{1/3}$	${\left(\frac{4\pi \tilde{n}}{g}\right)}^{1/2}$	${\left(\frac{4\pi \tilde{n}}{g}\right)}^{1/2}$	–
${\rho }_{\text{c}}\left(E\right),{\tilde{\rho }}_{\text{c}}\left(E\right)$	$\frac{g{\left(2m^{\ast }\right)}^{3/2}}{4{\pi }^2{\hslash }^3}{E}^{1/2}$	$\frac{g{m}^{\ast }}{2\pi {\hslash }^2}$	$\frac{gE}{2\pi {\hslash }^2{\upsilon }_{\text{F}}^2}$	(1.34)
${\rho }_{\text{c}}\left(E_{\text{F}}\right),{\tilde{\rho }}_{\text{c}}\left(E_{\text{F}}\right)$	$\frac{m^{\ast }}{\pi {\hslash }^2}{\left(\frac{3g^{2}n}{4\pi }\right)}^{1/3}$	$\frac{g{m}^{\ast }}{2\pi {\hslash }^2}$	$\frac{1}{\hslash {\upsilon }_{\text{F}}}{\left(\frac{g\tilde{n}}{\pi }\right)}^{1/2}$	–
$m^{\ast }$	Energy independent		$\frac{E}{{\upsilon }_{\text{F}}^2}$	(1.45)
${\upsilon }_{\text{F}}$	${\left(\frac{2E_{\text{F}}}{m^{\ast }}\right)}^{1/2}$		Energy independent	–
${\upsilon }_{\text{g}}$	${\left(\frac{2E}{m^{\ast }}\right)}^{1/2}$		${\upsilon }_{\text{F}}$	–

a: In this table, the parameters of an electron in a conduction band are listed; those of a hole in a valence band can be similarly expressed. $E$ is the energy of the electron measured with respect to the conduction band edge, by setting $E_{\text{c}}=0$ for the conduction band minimum of a parabolic system and $E_{\text{Dirac}}=0$ for graphene. $n$ is the electron density of a 3D system and $\tilde{n}$ is the electron density of a 2D system including graphene, in the limit that $T\to 0$ . $g$ is the total degeneracy including spin degeneracy.

1.6 Band Structure of Multilayer Graphene

When multiple layers of graphene are brought together to form a sheet of multilayer graphene, the interlayer interaction can fundamentally change the band structure. In this section, the band structures of AA and AB (Bernal) stacking orders are derived. Other possible stacking arrangements are discussed at the end.

1.6.1 AB Stacking Order

In AB-stacked multilayer graphene, the arrangement of layers follows the αβαβ … order. Figure 1.7 illustrates the AB stacking by showing one set of α and β layers. Nearest-neighbor intralayer and interlayer interactions are considered, as shown in Figure 1.7(a); high-order interactions beyond nearest neighbors are neglected for simplicity. The intralayer interaction energy is ${\gamma }_0$ , defined in (1.17), whereas the interlayer interaction energy is characterized by a different off-diagonal parameter ${\gamma }_1$ in the Hamiltonian. The top layer (β layer) is shifted relative to the bottom layer (α layer) by one C–C distance, as shown in Figure 1.7(b).

Figure 1.7 (a) Intralayer (${\gamma }_0$ ) and interlayer (${\gamma }_1$ ) interactions of AB-stacked multilayer graphene showing one set of α and β layers. (b) Top view of (a), showing the overlap between α (bottom) and β (top) layers. (c) Side view of (a). A chain of alternating A and B carbon atoms connected through interlayer interactions of energy ${\gamma }_1$ along the z direction is shown.

The band structure of N-layer AB-stacked multilayer graphene can be obtained in a manner similar to the derivation of the band structure of monolayer graphene. As shown in (1.27), after proper transformation the Hamiltonian at the $\text{K}$ point for intralayer interactions can be expressed in terms of the energy parameters ${\gamma }_{\pm }=\hslash {\upsilon }_{\text{F}}\left(k_x\pm \text{i}{k}_y\right)=3a{\gamma }_0\left(k_x\pm \text{i}{k}_y\right)/2$ with $\hslash {\upsilon }_{\text{F}}=3a{\gamma }_0/2$ . Therefore, the tight-binding Hamiltonian at the $\text{K}$ point for AB-stacked multilayer graphene can be expressed in a form similar to (1.27) as [Reference Min and MacDonald10]

${\widehat{H}}_{\text{K}}^{\text{AB}}=\left[\begin{array}{c}\begin{array}{cccccc}0& {\gamma }_{-}& 0& 0& 0& 0\\ {\gamma }_{+}& 0& {\gamma }_1& 0& 0& 0\\ 0& {\gamma }_1& 0& {\gamma }_{-}& 0& {\gamma }_1\\ 0& 0& {\gamma }_{+}& 0& 0& 0\\ 0& 0& 0& 0& 0& {\gamma }_{-}\\ 0& 0& {\gamma }_1& 0& {\gamma }_{+}& 0\\ & & & ⋮\end{array}\cdots \end{array}\right].$(1.47)

In (1.47), the basis is taken to be $\left\{|{\psi }_1^{\text{A}}〉,|{\psi }_1^{\text{B}}〉,|{\psi }_2^{\text{A}}〉,|{\psi }_2^{\text{B}}〉,\cdots ,|{\psi }_N^{\text{A}}〉,|{\psi }_N^{\text{B}}〉\right\}$ , with the eigenvector $|\psi 〉={\left[c_1^{\text{A}},c_1^{\text{B}},c_2^{\text{A}},c_2^{\text{B}},\dots ,c_N^{\text{A}},c_N^{\text{B}}\right]}^{\text{T}}$ , where $|{\psi }_n^{\text{A,B}}〉$ and $c_n^{\text{A,B}}$ correspond to the sublattice A or B of layer n. By applying the eigenvalue equation ${\widehat{H}}_{\text{K}}^{\text{AB}}|\psi 〉=E|\psi 〉$ , we find from (1.47) the following coupled difference equations,

$\{\begin{array}{l}E{c}_n^{\text{A}}={\gamma }_{-}{c}_n^{\text{B}}\\ E{c}_n^{\text{B}}={\gamma }_1\left(c_{n-1}^{\text{A}}+c_{n+1}^{\text{A}}\right)+{\gamma }_{+}{c}_n^{\text{A}}\end{array}\text{for}\text{odd}n,$(1.48)

$\{\begin{array}{l}E{c}_n^{\text{A}}={\gamma }_1\left(c_{n-1}^{\text{B}}+c_{n+1}^{\text{B}}\right)+{\gamma }_{-}{c}_n^{\text{B}}\\ E{c}_n^{\text{B}}={\gamma }_{+}{c}_n^{\text{A}}\end{array}\text{ for}\text{even}n,$(1.49)

where $n=1,2,\dots ,N$ and $c_0^{\text{A}}=c_0^{\text{B}}=c_{N+1}^{\text{A}}=c_{N+1}^{\text{B}}=0$ .

We first check if $E=0$ is a nontrivial solution for the coupled equations of (1.48) and (1.49). For $E=0$ , we have ${\gamma }_{-}{c}_n^{\text{B}}=0$ for odd n and ${\gamma }_{+}{c}_n^{\text{A}}=0$ for even n. If we assume nonzero ${\gamma }_{-}$ and ${\gamma }_{+}$ , we obtain from (1.48) and (1.49) a trivial solution with all c coefficients being zero. However, if we assume ${\gamma }_{-}={\gamma }_{+}=0$ , we have rows in (1.47) that are identically zero; thus, $c_n^{\text{B}}$ for even n and $c_n^{\text{A}}$ for odd n are decoupled and can assume any value. Other coefficients can be either zero or nonzero, determined by (1.48) and (1.49) and the number of layers. Therefore,

$E={\gamma }_{+}={\gamma }_{-}=0$(1.50)

represents a nontrivial solution for the coupled equations of (1.48) and (1.49).

For $E\ne 0$ , we can simplify (1.48) and (1.49) as

$\left(E-{\gamma }_{+}{\gamma }_{-}/E\right)c_n^{\text{B}}={\gamma }_1\left(c_{n-1}^{\text{A}}+c_{n+1}^{\text{A}}\right)\text{ for odd}n,$(1.51)

$\left(E-{\gamma }_{+}{\gamma }_{-}/E\right)c_n^{\text{A}}={\gamma }_1\left(c_{n-1}^{\text{B}}+c_{n+1}^{\text{B}}\right)\text{ for even}n.$(1.52)

Then, by defining $c_n=c_n^{\text{B}}$ for odd n and $c_n=c_n^{\text{A}}$ for even n, the coupled equations of (1.51) and (1.52) can be expressed in a single simplified form:

$E_0{c}_n={\gamma }_1\left(c_{n-1}+c_{n+1}\right),$(1.53)

where $E_0=E-{\gamma }_{+}{\gamma }_{-}/E$ . Thus we have transformed the original problem expressed in the coupled equations of (1.48) and (1.49) into a chain of N atoms with eigenenergy $E_0$ while considering only the nearest-neighbor interactions characterized by the interlayer interaction energy ${\gamma }_1$ , as shown in Figure 1.7(c). With the constraint that $c_0=c_{N+1}=0$ , we find the solutions of (1.53) as

$c_n=\text{sin}\left(\pm r\frac{\pi n}{N+1}\right),$(1.54)

where $r=1,2,\dots ,N$ . Therefore, we find N eigenvectors for N-layer graphene. By plugging (1.54) back into (1.53), we find the eigenenergies:

$E_0=2{\gamma }_1\text{cos}\left(\pm r\frac{\pi }{N+1}\right).$(1.55)

To gain more physical insight for (1.55), we can regard N-layer graphene as a Fabry–Pérot cavity in which an electromagnetic wave is confined in the direction perpendicular to the multilayer graphene surface, which is taken to be the z direction [Reference Mak, Sfeir, Misewich and Heinz11]. The boundary condition $c_0=c_{N+1}=0$ implies that the wave function vanishes at $n=0$ and $n=N+1$ beyond the first and $N\text{th}$ layers. The distance between the locations of $n=0$ and $n=N+1$ layers is $\left(N+1\right)d$ ; therefore, the effective length of the Fabry–Pérot cavity is $d_{\text{eff}}=\left(N+1\right)d$ , where $d=0.335\mathrm{nm}$ is the distance between two neighboring graphene layers in AB stacking. A standing wave can form within the multilayer graphene structure with two opposite wave vectors of equal magnitude in the z direction:

$k_z=\pm r\frac{\pi }{d_{\text{eff}}}=\pm r\frac{\pi }{\left(N+1\right)d}.$(1.56)

Therefore, as the thickness of the multilayer graphene increases with the number N of layers, the number of the allowed values of quantized $k_z$ also increases. By using (1.56), we can rewrite (1.55) as

$E_0=2{\gamma }_1\text{cos}{k}_{z}d,$(1.57)

which gives the eigenenergy of N interacting carbon atoms with $k_z$ being the wave number along the direction normal to the graphene surface. Indeed, if the in-plane wave numbers have zero values so that $k_x=k_y=0$ and ${\gamma }_{+}={\gamma }_{-}=0$ , the original problem represented by (1.48) is tantamount to the chain of carbon atoms along the z direction shown in Figure 1.7(c). The eigenenergy is simply $E=E_0$ because ${\gamma }_{+}={\gamma }_{-}=0$ in the relation $E_0=E-{\gamma }_{+}{\gamma }_{-}/E$ .

When $k_x$ and $k_y$ are not both zero so that ${\gamma }_{+}\ne 0$ and ${\gamma }_{-}\ne 0$ , the eigenenergy E can be found from the relation $E_0=E-{\gamma }_{+}{\gamma }_{-}/E$ , which gives

$E=\frac{1}{2}\left(E_0\pm \sqrt{4{\left(\hslash {\upsilon }_{\text{F}}{k}_{\parallel }\right)}^2+E_0^2}\right),$(1.58)

where $k_{\parallel }={\left(k_x^2+k_y^2\right)}^{1/2}$ is the in-plane wave number, in contrast to the out-of-plane wave number $k_z$ . In the limit that $\left|E_0\right|\gg \hslash {\upsilon }_{\text{F}}{k}_{\parallel }$ , (1.58) has two small-$k_{\parallel }$ solutions:

$E=E_0+\frac{{\left(\hslash {\upsilon }_{\text{F}}{k}_{\parallel }\right)}^2}{E_0}$(1.59)

and

$E=-\frac{{\left(\hslash {\upsilon }_{\text{F}}{k}_{\parallel }\right)}^2}{E_0}.$(1.60)

Note that in both (1.59) and (1.60), $E$ varies quadratically, not linearly, with $k_{\parallel }$ . By defining the effective mass for these quadratic bands as $m^{\ast }=\left|E_0\right|/2{\upsilon }_{\text{F}}^2$ , the first solution (1.59) has the form of $E=E_0\pm {\left(\hslash k_{\parallel }\right)}^2/2m^{\ast }$ , and the second solution (1.60) has the form of $E=\mp {\left(\hslash k_{\parallel }\right)}^2/2m^{\ast }$ , where the upper signs are taken when $E_0>0$ and the lower signs are taken when $E_0<0$ . Accordingly, the energy bands given by (1.59) and (1.60) are called massive bands, in contrast to the massless Dirac bands of monolayer graphene.

For even N, there are N different values of $E_0$ given by (1.55). Therefore, (1.59) and (1.60) respectively give N energy bands and thus together a total of 2N massive energy bands. Because half of $E_0$ values are positive and the other half are negative, there are N conduction bands and N valence bands. As an example, the band structure of AB-stacked bilayer graphene ($N=2$ ) is shown in Figure 1.8(a), where the solid curves are obtained from (1.58) without approximation and the dashed curves are obtained from (1.59) and (1.60) under the small-$k_{\parallel }$ approximation, using ${\gamma }_0=3.16\mathrm{eV}$ and ${\gamma }_1=0.37\mathrm{eV}$ obtained from experiments [Reference Abergel, Apalkov, Berashevich, Ziegler and Chakraborty4]. As can be seen, there are two bands for each value of r (marked next to the curves), and there are a total of $2N=4$ bands for the bilayer system, all of which are massive bands. The effective mass in the small-$k_{\parallel }$ limit is $m^{\ast }=0.032m_{\text{e}}$ , where $m_{\text{e}}=9.11\times {10}^{-31}\text{kg}$ is the rest mass of a free electron. The eigenvectors can be obtained accordingly from (1.54) for $c_1^{\text{B}}=c_1$ and $c_2^{\text{A}}=c_2$ , which in turn give $c_1^{\text{A}}$ and $c_2^{\text{B}}$ from (1.48): $|\psi 〉={\left[\frac{{\gamma }_{-}{c}_1}{E},c_1,c_2,\frac{{\gamma }_{+}{c}_2}{E}\right]}^{\text{T}}$ , where the eigenenergy E is given in (1.58).

For odd N, we can always find a value of r such that $E_0=0$ given by (1.57), for which (1.58) gives the linear bands $E=\pm \hslash {\upsilon }_{\text{F}}{k}_{\parallel }$ . Therefore, there are $2\left(N-1\right)$ massive conduction and valence bands described by (1.59) and (1.60), and two linear bands defining one massless Dirac cone given by $E=\pm \hslash {\upsilon }_{\text{F}}{k}_{\parallel }$ . As an example, the band structure of AB-stacked trilayer of $N=3$ is shown in Figure 1.8(b), where six bands can be seen; four are massive bands and two are massless linear Dirac bands.

Figure 1.8 Band structures of AB-stacked (a) bilayer and (b) trilayer graphene. The solid curves in (a) for both $r=1$ and $2$ and those in (b) for $r=1$ and $3$ are obtained from (1.58), and the dashed curves are the small-$k_{\parallel }$ approximation obtained from (1.59). The linear bands for $r=2$ in (b) are the Dirac bands.

From the above discussions, N-layer AB-stacked graphene has $2N$ bands. For even N, all $2N$ bands are massive bands. For odd N, there are $2\left(N-1\right)$ massive bands plus two massless linear Dirac bands.

The band structure of natural graphite follows the AB stacking order; it can be obtained from (1.56) and (1.57) by taking the limit that $N\to \infty$ in (1.54). In this limit, $k_z$ is no longer discrete but forms a continuous spectrum. In a reduced Brillouin zone scheme, $k_z$ spans from $-\pi /2d$ to $\pi /2d$ , where the H point is located, as shown in Figure 1.9 [Reference Slonczewski and Weiss12]. At $k_{\parallel }=0$ , there are three energy bands: one energy band for $E=0$ given by (1.50), and two energy bands for $E=E_0$ given by (1.57) with $E_0>0$ for a conduction band and $E_0<0$ for a valence band, as shown in Figure 1.9. By using (1.58), one can map out the band structure of graphite as a function of $k_{\parallel }$ at different values of $k_z$ . The band structures of bilayer and trilayer graphene are shown in Figure 1.8. For bilayer graphene, two bands at $k_z=\pm \pi /3d$ , respectively, in the reduced zone scheme are shown in Figure 1.9(a). For trilayer graphene, the massive bands at $k_z=\pm \pi /4d$ and the massless Dirac bands at $k_z=\pm \pi /2d$ are also plotted, as shown in Figure 1.9(b). Note that the mirrored band structure does not represent double degeneracy because the positive and negative values of $k_z$ do not represent two independent eigenstates but represent two counter-propagating waves that form the standing wave in multilayer graphene. As can be seen from (1.54) and (1.56), for any solution of $k_z$ , the values of $c_n$ obtained using the positive and negative values of $k_z$ only differ by a negative sign. Therefore, both sets of the coefficients $c_n$ obtained from the positive and negative values of $k_z$ give the same eigenstate $|\psi 〉$ . The truly independent states are the standing waves of different wave numbers represented by different values of $\left|k_z\right|$ .

Figure 1.9 Band structures of AB-stacked (a) bilayer graphene and (b) trilayer graphene, obtained by cutting the band structure of graphite at certain $k_z$ values.