The band structure of graphene

Graphene has a hexagonal lattice that can be described using the lattice vectors ${𝐫}_0=a(1,0,0)$, ${𝐫}_1=a(-1/2,\sqrt{3}/2,0)$, and ${𝐫}_2=a(0,0,1)$. We set the lattice constant to a = 2.5 Å, which only is meant to reflect the correct order of magnitude. Strictly speaking, only two two-dimensional vectors are required to describe the lattice. However, we have here made the vectors three-dimensional and included a third lattice vector ${𝐫}_2$ that is perpendicular to the lattice to simplify the calculation of the reciprocal lattice vectors. The resulting lattice is shown in the image below.

Indicated in the image is also the unit cell and the A and B sublattice sites within the unit cell. We further define

which can be verified to be the vector from A to its three nearest neighbor B sites.

The reciprocal lattice vectors can now be calculated using

Here ${𝐤}_2$ just like ${𝐫}_2$ is related to the $z$-direction and will be dropped from here on. The two remaining reciprocal lattice vectors ${𝐤}_0$ and ${𝐤}_1$ defines the reciprocal lattice of interest to us. Below the reciprocal lattice is shown.

In the the image above we have also indicated the first Brillouin zone in red, and outlined the path $\mathrm{\Gamma }\to M\to K\to \mathrm{\Gamma }$ along which the band structure will be calculated. The high symmetry points can be verified to be given by

We will use a simple tight-binding model of graphene of the form

Here t is a hopping amplitude connecting nearest neighbor sites which we take to be t = 3 eV (again, only correct up to order of magnitude). The index $𝐢$ runs over the unit cells, while $𝜹$ takes on the three values that makes $𝐢+𝜹,B$ the nearest neighbor of $𝐢,A$.

Because of translation invariance, it is possible to transform the expression above to a simple expression in reciprocal space. To do so we write

Here ${𝐑}_{𝐢}$ is the position of the $A$ atom in unit cell $𝐢$, while ${𝐫}_{𝜹}\in \{{𝐫}_{AB}^{(0)},{𝐫}_{AB}^{(1)},{𝐫}_{AB}^{(2)}\}$. ¹¹1 Note that the definition of $c_{𝐢+𝜹,B}$ is threefold redundant since for each site $B$ in unit cell $𝐢$ there are three ways to label the same site. For each of the three vectors ${𝐫}_{AB}^{(0)}$, ${𝐫}_{AB}^{(1)}$, and ${𝐫}_{AB}^{(2)}$ there is a different unit cell from which the $B$ site in cell $𝐢$ can be reached. However, each of these three redundant definitions can be verified to be equivalent and the redundant expression simplifies the mathematical notation.

Inserting this into the expression above we get

In the second line we have summed over $𝐢$ and used the delta function property of ${\sum }_{𝐢}{e}^{i(𝐤-{𝐤}^{\prime })\cdot {𝐫}_{𝐢}}/N$ to eliminat ${𝐤}^{\prime }$, while in the last line we have made the sum over $𝜹$ explicit.

We are now ready to implement the calculation using TBTK. To do so, we begin by specifying the parameters to be used in the calculation

In lines 2-5, we specify what units that dimensionful parameters should be understood in terms of. For our purposes here, this means that energies and distances are measured in electron Volts (eV) and Ångström (Ao), respectively.

We next specify the value for t and a as given in the physical description above. The value ‘BRILLOUIN_ZONE_RESOLUTION’ determines the number of subintervals each of the reciprocal lattice vectors is divided into when generating a mesh for the Brillouin zone. ‘numMeshPoints’ combine this information into a vector. Similarly, ‘K_POINTS_PER_PATH’ determines the number of $k$-points to use along each of the paths $\mathrm{\Gamma }\to M$, $M\to K$, and $K\to \mathrm{\Gamma }$ when calculating the band structure. Finally, the energy window that is used to calculate the DOS is specified in the three last lines.

The lattice vectors ${𝐫}_0$, ${𝐫}_1$, and ${𝐫}_2$ are now defined as follows.

The nearest neighbor vectors ${𝐫}_{AB}^{(0)}$, ${𝐫}_{AB}^{(1)}$, and ${𝐫}_{AB}^{(2)}$ are similarly defined using

Next, the reciprocal lattice vectors are calculated using the expressions at the beginning of the Brillouin zone section.

Note that here ‘*’ indicates the cross product, while the scalar product between u and v is calculated using Vector3d::dotProduct(u, v).

Having defined the lattice vectors, we can now construct a Brillouin zone.

Here the first argument is a list of lattice vectors, while the second argument specifies the type of mesh that is going to be associated with this Brillouin zone. SpacePartition::MeshType::Nodal means that when the Brillouin zone is divided into small parallelograms, the mesh points will lie on the corners of these parallelograms. If SpacePartition::MeshType::Interior is used instead, the mesh points will lie in the middle of the parallelograms.

Having constructed a Brillouin zone, we finally generate a mesh by passing the number of segments to divide each reciprocal lattice vector to the following function.

The model is set up using the following code.

Here the loop runs over each lattice point in the mesh. To understand line 4-7, we first note that the mesh contains $k$-points of the form $(k_x,k_y)$ with real numbers $k_x$ and $k_y$. However, to specify a model, we need discrete indices with integer-valued subindices to label the points. Think of this as each $k$-value being of the form ${𝐤}_{𝐢}$, where ${𝐤}_{𝐢}$ is a real-valued vector while $𝐢$ is an integer-valued vector that indexes the different $k$-points. The BrillouinZone solves this problem for us by automatically providing a mapping from continuous indices to discrete indices. Given a real index in the first Brillouin zone and information about the number of subdivisions of the lattice vectors, brillouinZone.getMinorCellIndex() finds the closest discrete point and returns its Index representation. This index can then be used to refer to the k-point whenever an Index object is required.

In line 10, the mesh point is converted to a Vector3d object that is then used in lines 11-30 to calculate the matrix element we derived in the section about the Hamiltonian. The matrix element, together with its Hermitian conjugate, is then fed to the model in lines 33-37. Here we construct the full index for the model by combining the integer kIndex[0] and kIndex[1] values, together with the sublattice index, into an index structure of the form kx, ky, sublattice. A zero in the sublattice index refers to an A site, while one refers to a B site. Finally, the model is constructed in line 39.

We are now ready to solve the model and will use diagonalization to do so. In this case, we have a problem that is block-diagonal in the k-index. For this reason, we will use the solver Solver::BlockDiagonalizer that can take advantage of this block structure. To set up and run the diagonalization we write

We then set up the corresponding property extractor and configure the energy window using

Extracting, smoothing, and plotting the DOS is now done as follows.

We are now ready to calculate the band structure. To do so, we begin by defining the high symmetry points

These are then packed in pairs to define the three paths Γ→M, M →K, and K →Γ

We now loop over each of the three paths and calculate the band structure at each k-point along these paths.

Here we first set up the array ‘bandStructure’ that is to contain the band structure with two eigenvalues per k-point. The variable ‘interpolator’ is then defined, which is an array of values from 0 to 1 that will be used to interpolate between the paths start and endpoints. In lines 9 and 10 we then select the start and endpoints for the current path, and the loop over this path begins in line 13.

In lines 20-23, we calculate an interpolated k-point value between the first and last points in the given path. Similarly, as when setting up the model, we then request the Index that corresponds to the given k-point. Having obtained this Index, we extract the lowest eigenvalue for the corresponding block in lines 34-39 and store it in the array ‘bandStructure’. Similarly, the highest eigenvalue is obtained for the same block in lines 40-45.

Next, we calculate the maximum and minimum value of the band structure along the calculated paths. We do this to allow for drawing vertical lines of the correct height at the M and K points when plotting the data.

Finally, we plot the band structure.

Here the ‘_a_’ flag in the call to bandStructure.getSlice(0, _a_, “color”, “black”) on line 5 (and 6) is a wildcard. It indicates that we want to get a new array that contains all the elements of the original array that has a zero in the first index and an arbitrary index in the second. ‘a’ stands for ‘all’.