Remarks on the tight-binding model of graphene

We address a simple but fundamental issue arising in the study of graphene, as well as of other systems that have a crystalline structure with more than one atom per unit cell. For these systems, the choice of the tight-binding basis is not unique. For monolayer graphene two bases are widely used in the literature. While the expectation values of operators describing physical quantities should be independent of basis, the form of the operators may depend on the basis, especially in the presence of disorder or of an applied magnetic field. Using an inappropriate form of certain operators may lead to erroneous physical predictions. We discuss the two bases used to describe monolayer graphene, as well as the form of the most commonly used operators in the two bases. We repeat our analysis for the case of bilayer graphene.


Introduction
A peculiar characteristic of graphene is the presence of two atoms per unit cell.The solidstate theory for such systems necessitates the introduction of multi-dimensional tight-binding bases, whose choice is not unique.The expectation values of physically measurable quantities are of course independent of basis; however, in practice this is oftentimes not straightforward to see.In particular, if the expectation values of certain operators are to be independent of basis, their form must be basis-dependent.
There appears to exist a rather bit of confusion in the literature about the form of various operators in the two tight-binding bases most commonly used to describe graphene.The operators that are most commonly misidentified are the k-space Hamiltonian, the density, the density of states, and the single-impurity potential.Some of these operators are used to describe the effects of impurity scattering in graphene [1,2,3,4,5,6].Using the correct form of these operators is essential in correctly computing the density of states in the presence of impurities, which is measured in STM experiments [7,8,9].
Our purpose is to clarify the subtleties associated with the correct form of these operators.We present carefully the two bases, and write down the tight-binding Hamiltonian and its low energy expansion in first-quantized language.We also describe the corresponding secondquantized formalism, and show that the choice of basis is equivalent to choosing the manner of taking the Fourier transform of the second-quantized operators.This allows us to write down the form of various operators in the two languages.
For monolayer graphene, one can choose a basis [10] in which only one point per unit cell is used as the origin for the Bloch wave-functions.This basis consists of two p z orbital wavefunctions centered on the two carbon atoms of the unit cell; these wavefunctions have the same phase factor, determined by the position of the "origin" of the unit cell.Alternatively, one can use a second basis, in which the positions of the two atoms in the unit cell are used as "centers" for Bloch's theorem; hence the second basis also consists of two p z orbital wavefunctions centered at the two carbon atoms, but their phase factors (determined by the position of the corresponding atom) are different [11,12].
Bilayer graphene on the other hand has four atoms per unit cell.Consequently, there are at least two choices of tight-binding basis.We present the canonical form, which is widely used in the literature [13,14], and in which all four p z orbital wavefunctions have different phases (given by the positions of the four atoms in the unit cell).We also discuss an alternative basis, in which the four wavefunctions have the same phase factor.
In section 2 we present the two tight-binding bases and the tight-binding Hamiltonian for monolayer graphene and its low energy expansion using a first-quantized formalism and Bloch's theorem.In section 3 we present the second-quantized formalism.In section 4 and section 5 we present the density operator, and the impurity potential respectively.In section 6 we discuss the case of bilayer graphene and we conclude in section 7. Given the honeycomb hexagonal lattice of graphene with two atoms per unit cell, one can use Bloch's theorem to write down the eigenstates of the lattice Hamiltonian.In the tight-binding approximation, one searches for eigenfunctions of the Hamiltonian as linear combinations Ψ k ( r) of atomic wave functions.A common representation of this combination is

Lattice considerations
where N is the number of elementary cells, and the functions φ( r) are the wave-functions of the p z orbitals of the carbon atoms.As described below, the coefficients c A/B I are chosen such that Ψ k ( r) is an eigenstate of the tight-binding Hamiltonian.The vectors R j = n a 1 + m a 2 with j = (n, m) specify the position of one graphene unit cell, with a 1 = a √ 3x/2 + 3aŷ/2, and a 2 = −a √ 3x/2 + 3aŷ/2, where a is the distance between two nearest neighbors.Also, R A/B j are the positions of the A and B atoms respectively.For simplicity we took the positions of the unit cells to be given by the positions of the A atoms, In our choice of the coordinate system, the B atoms are located at R B j = R j + δ 3 , where the vector δ 3 ≡ δ AB is one of the three vectors connecting an atom A with its three nearest neighbors: aŷ/2 and δ 3 = −aŷ, as depicted in Fig. 1.Note that the choice of the origin, as well as of the axes of the coordinate system is arbitrary, but once the choice has been made it has to be used consistently in later analysis.
In this representation of the tight-binding Hamiltonian eigenstates, one first constructs a combination of the atomic wave functions within the unit cell, then attaches a phase factor to each cell to construct a Bloch function.This is the "textbook procedure" (see for example Ashcroft and Mermin Eq. 10.26 [10]).
In the second representation one writes the Hamiltonian eigenstates as linear combinations of two Bloch functions corresponding respectively to the A and B atoms, but with a different phase factor attached to each atom A and B.
This second representation is used for example in the paper by Wallace on the band structure of graphite [11], and in many recent papers on graphene [12].Note that in each representation we have chosen a tight-binding basis {Ψ Ak ν ( r), Ψ Bk ν ( r)} where ν = I/II, and Ψ ).We can see that the two bases differ by relative phase factors between their components.The eigenstates of the tight-binding Hamiltonian are linear combinations of each basis wavefunctions.We will show that, while the coefficients of the linear combinations are basis-dependent, the eigenfunctions of the tight-binding Hamiltonian are the same in both bases.Also, the expectation value of any physical quantity is independent of the basis chosen.

Tight-binding Hamiltonian
The tight-binding Hamiltonian used to describe graphene allows for hopping between nearest neighbors (j, A) and (i, B), such that electrons on an atom of the type A/B can hop on the three nearest B/A atoms respectively.Thus we can write where |φ A/B j is the standard notation for wavefunctions φ ).The eigenequations for the coefficients c A ( k) and c B ( k) in Eqs.(1,3) are straightforwardly obtained from evaluating φ A/B j |H|Ψ k ν using Eq.( 1), where in the first basis, or in the second basis.Defining the Hamiltonian density is written as (ν = I or II) with the eigenvalues We should note that the eigenvalues of the Hamiltonian (which give the energy dispersion of the two bands of graphene) are the same in both bases, as expected.This is because in the two representations, the two functions f I and f II differ simply by a phase factor: where δ AB ≡ δ 3 = −aŷ is the vector connecting the A and B atoms in a unit cell.
Given that f , the Hamiltonian density in the first representation can be also rewritten as The k dependence of this phase is shown in Figure (2).One can see clearly the two inequivalent Brillouin zone corners K and K ′ .Each of the two points is equivalent to all the points that can be be obtained by translations with the reciprocal lattice vectors.
The k dependence the phase θ I ( k) is represented by small segments in the two-dimensional k space.One sees clearly the two inequivalent BZ corners K and K ′ having different topologies, and being characterized by opposite Berry phases [15].
In the second representation, the Hamiltonian carries an inconvenient phase : We can go back and rewrite the eigenfunctions of the tight-binding Hamiltonian in the two bases: where the ± signs correspond to the eigenfunctions describing the conduction band, and the valence band respectively.In the second representation,

Considering that θ
, one checks easily that the two representations lead to the same expression of the eigenfunctions Ψ k II ( r) = Ψ k I ( r).

Low energy expansions
As well known, the energy vanishes at the Dirac points, which are at Here ξ = ± is the valley index (there are two such points for each elementary cell of the reciprocal space).Each point is equivalent to all the points in the reciprocal space that have the same ξ but different (m, n) and that can be be obtained by translations with the reciprocal lattice vectors.The ξ = ± pair that is chosen most often contains two corners of the first Brillouin zone, K ≡ K + 00 and K ′ ≡ K − 00 , as described in Figs.1,2,3.Note that Thus we can expand the Hamiltonian in the first basis around the Dirac points to find where k = K ξ mn + q and ξ = ±1 is the valley index.We see that in this basis the expansion does not depend on the choice of (m, n).The six corners of the BZ appear thus equivalent to either K or K ′ and can be recovered by a translation of K and K ′ by various reciprocal lattice vectors.Given the above choice of vectors a 1 and a 2 , one obtains where v = 3t/(2a).We can now write the low-energy Hamiltonian density in the 4 × 4 space defined by ( KA, KB, K ′ B, K ′ A) as: which can be expressed in the compact form: where σ and τ are the usual Pauli spin matrices.Alternatively we can write the low-energy Hamiltonian density in Eq. ( 17) as: where θ( q) = arctan(q y /q x ).
In the second basis, the expression of the Hamiltonian is less convenient because it contains the phase factor e −i K ξ mn • δ AB = e i K ξ mn •ŷa where is independent of the valley index ξ, but depends on the index (m, n).This makes the six corners of the BZ appear inequivalent.Thus in basis II, in the 4 × 4 space defined by ( K + mn A, K + mn B, K − mn B, K − mn A), the Hamiltonian density is where the phase factor z mn = e 2iπ(m+n)/3 depends on the choice of the vector K ± mn in the reciprocal space.In the standard choice for the two valley-points ( K = K + 00 and K ′ = K − 00 ) we have m = n = 0 and the low-energy expansion of the Hamiltonian is the same in both bases.

Second quantization
In the second quantized formalism we can define the operators a † j , b † j that correspond to creating electrons on the sublattices A and B, at sites R A j and R B j respectively.From Eq.(1) we see that the Fourier transform (FT) of the a j and b j operators should depend on the basis, such that and where the sum is taken over all lattice unit cells.The inverse Fourier transform of these operators will be: where we define k∈BZ ≡ BZ d 2 k S BZ , and S BZ = 8π 2 /3 √ 3. Given the choice for the origin of the unit cell, R j = R A j , we can see easily that a Thus the change of basis described in the previous section introduces a different momentum-dependent phase factor in the definition of the k-space Fourier-transformed operators.
In the second quantized formalism we can write the tight-binding Hamiltonian as: where t is the nearest neighbor hoping amplitude, and ij denotes summing over the nearest neighbors.In momentum space the tight-binding Hamiltonian becomes: where ν = I/II, and the f functions are defined in Eqs.(7,8) in the previous section.We can see that, exactly like in the first-quantized formalism, the form of the Hamiltonian is unique in real space, but it depends on the basis in momentum space.

The density and density of operators
It is quite interesting to keep track consistently of the correct form of a few other operators in the two bases.We first focus on the local charge-density operator: In the absence of disorder, the density will be independent of position.However, if impurities are present the density will fluctuate, and it is useful to define its Fourier transform: whose expectation value can be related to the results of FTSTS measurements [7,8,9].In basis I, the FT of the charge density becomes (see the Appendix for the complete derivation): Similarly we can redo the analysis in the basis II: Note that for systems that conserve momentum (translationally invariant), as in the absence of disorder and magnetic fields, only the q = 0 term is non-zero.Furthermore, if one is interested in average quantities, only the q = 0 term is relevant.Hence, in these cases the operators have the same form in the two bases.
Another operator of interest is the local density of states (LDOS), given by the number of electrons of energy ω at a given position.Its integral over ω gives the total density described above.The previous formulas can be trivially extended to the local density of states, by taking all operators at a specific energy ω.
The expectation values of the density of states operator at various positions on the two sublattices and at energy ω are given by: and S BZ is performed over the first BZ, and q ≡ d 2 q 4π 2 is performed over the entire reciprocal space.One can straightforwardly show that if r = R A j only the a † a terms contribute, and the b † b terms vanish; conversely, if r = R B j only the b † b terms contribute, and the a † a terms vanish.
Evaluating the density of A and B electrons in the unit cell is also different in the two bases.Since in basis I both the A and the B operators are defined at the origin of the unit cell, both densities have to be evaluated at this position (which we chose to be the position of the A atom, R A j ).In the basis II, the density of states is evaluated for each atom at its corresponding position ( R A j or R B j ).

Impurity potential
We can also write down the form of a delta-function impurity potential.For an impurity located on sublattice A, we have while for an impurity on the sublattice B In the case of a single impurity, it is most convenient to choose the origin of the coordinate system such that R imp j = 0. Thus, in basis I the impurity potential will be independent of momentum, regardless of whether the impurity is on the A or on the B site.
In basis II, R A imp j = 0, and no phase factors will appear when the impurity is on sublattice A (at the origin of the coordinate system).However, when the impurity is on sublattice B, R B imp j = δ AB , and a momentum-dependent phase factor will appear in the form of the impurity potential.We should note that this phase factor comes from choosing the origin of the coordinate system on an A atom.Indeed, if one performs a FT of the Friedel oscillations generated by an impurity at a B atom while using a coordinate system with the origin at a neighboring A atom, one generates a momentum-dependent phase factor in the FT.This can be eliminated by changing the origin of the coordinate system from the A atom to the B atom, and by carefully tracking the change in the form of the other operators.

Bilayer graphene
We can generalize the formalism presented in the previous sections to systems with arbitrary numbers of electrons per unit cell.Bilayer graphene is made of two coupled graphene monolayers (see Fig. 4), and there are four atoms per unit cell, two for each layer.In real space the δ δ δ a a The operators a † i , b † i denote the creation of particles at sites A and B in layer 1, while ã † i , b † i denote the creation of particles at sites Ã and B in layer 2. The sites B in the first layer lie on top of the sites Ã in the second layer, and there is a non-zero t p hopping of electrons between them.Also ij 1,2 denotes summing over the nearest neighbors in layers 1 and 2 respectively; jc denotes summing only over the sites B in the first layer which are on top of sites Ã in the second layer.
As for monolayer graphene, we can define two types of Fourier transform, consistent with two different tight-binding bases.In the first basis one first constructs a combination of the atomic wave functions within the unit cell, then attaches a phase factor to each cell to construct a Bloch function.The corresponding Fourier transformed operators in secondquantized formalism are given by: The vectors R j = n a 1 + m a 2 , with j = (n, m), specify the position of the unit cell of the top layer, which we also take to be the origin of the four-atom unit cell of bilayer graphene (see Fig. 4).Here we chose the origin of the coordinate system on an A atom in layer 1.
In the second basis the positions of the four atoms in the unit cell are used as "centers" for Bloch's theorem, and the corresponding Fourier transformed operators are: where R A j = j and R B j = R j + δ AB are the positions of the A and B atoms in layer 1, while R Ã j = R j + δ AB and R B j = R j + 2 δ AB are the positions of the Ã and B atoms in layer 2. This is the basis that is most often used in the literature to describe bilayer graphene.
In momentum space the tight-binding Hamiltonian becomes: where ν = I/II and the f 's are the same as the ones defined in Eqs.(7,8) for monolayer graphene.
The density operator is given by and its Fourier transform is: where ν = I/II, β I = αI = e i q• δ AB , and βI = e 2i q• δ AB , while in basis II there are no relative phase factors, β II = αII = βII = 1.We can also evaluate the density and the density of states at various positions: where as before k∈BZ ≡ BZ d 2 k S BZ is performed over the first BZ, while q ≡ d 2 q 4π 2 is performed over the entire reciprocal space.
Note that (like in the case of monolayer graphene) when working in basis I, all the four densities of states are evaluated at the position of the unit cell vector R j (which we chose to be the position of the A atom, R A j ).In II however, the density of states is evaluated for each atom at its corresponding position: R A/B/ Ã/ B j .

Conclusions
We analyzed the two tight-binding bases used to describe monolayer graphene.We showed that, while the eigenstates of the tight-binding Hamiltonian, as well as the expectation values of physical quantities are independent of basis, the form of certain operators depends on the basis.We also showed that the choice of basis is equivalent to the choice of the manner of performing the Fourier transforms of second-quantized operators.We wrote down in the two languages the Hamiltonian, the density, and the local density of states (LDOS), as well as the impurity potential.We also analyzed the case of bilayer graphene and presented two possible choices of tight-binding basis, and the form of the aforementioned operators in these bases.
For the case of the Fourier-transformed density operator, it is important to note that due to the (arbitrary) choice of the coordinate system, its expectation value can only be related to the FT of the experimental data if this FT is taken using the same coordinate system (axes and origin as depicted in Fig. 1).If the FT is taken using a different coordinate system, a momentum-dependent phase factor is introduced and needs to be accounted for before comparing theory and experiment.This inadvertence may lead in some cases to a simple rotation of the data, but other more complicated phase factors can also be introduced by a mismatch of the origins of the coordinate systems.For the case of a single impurity, it is most convenient to use a coordinate system with the origin at the impurity site.However, if multiple impurities are present, one needs to keep track of the relative phase factors introduced by their spatial distribution.
We should also comment that our careful tracking of the phase factors generated by the change of basis is in general not relevant if one is only interested in uniform properties or in spatial averages.However, a careful analysis of the phase factors is crucial if one studies systems with disorder, in the presence of an applied magnetic field, and more generally with broken translational invariance.

Appendix
In basis I the FT of the charge density is written as: The sum over the lattice unit cells j can be performed to obtain RL δ( k− k ′ + q+ Q RL ), where Q RL is any vector of the reciprocal lattice.However, as the integral over k ′ is constrained to the BZ, not all the terms of the sum contribute to the result, but only those for which k ′ = k + q + Q RL is in the Brillouin zone.For each q and k there is an unique Q RL that satisfies this condition.Consequently we have: However, given the FT definitions in Eq.( 19), a I ( k + Q RL ) = a I ( k), and b I ( k + Q RL ) = b I ( k), so we have: Similarly we can redo the analysis for the second basis: The sum over the sites R A j can be performed to obtain: RL δ( k − k ′ + q + Q RL ), where Q RL is again any vector of the reciprocal lattice.However, the sum over the sites R B j gives RL exp(−i Q RL • δ AB )δ( k − k ′ + q + Q RL ).Consequently we get: From the definitions in Eq.(20) we see that a II ( k + Q RL ) = a II ( k) and b II ( k + Q RL ) = b II ( k)e i Q RL • δ AB , and thus

Figure 1 :
Figure 1: Hexagonal honeycomb lattice of graphene (a), and its band structure (b).In b) the equal energy contours are drawn, and the Brillouin zone is indicated by dashed lines.The Dirac points K and K ′ are marked by arrows, and the reciprocal lattice vectors a * 1,2 are also drawn.

Figure 3 :
Figure 3: The k dependence of the phase θ II ( k) carries an inconvenient addition, so that all the six K points of the first Brillouin zone appear different.