The band structure of graphene

Most recent TBTK release at the time of writing: v1.1.1
Updated to work with: v2.0.0

In condensed matter physics, the electronic band structure is one of the most commonly used tools for understanding the electronic properties of a material. Here we take a look at how to set up a tight-binding model of graphene and calculate the band structure along paths between certain high symmetry points in the Brillouin zone. We also calculate the density of states (DOS). A basic understanding of band theory is assumed and if further background is required a good starting point is to look at the nearly free electron model, Bloch waves, Bloch’s theorem, and the Brillouin zone.

Physical description

Lattice

Graphene has a hexagonal lattice that can be described using the lattice vectors $\mathbf{r}_0 = a(1, 0, 0)$ , $\mathbf{r}_1 = a(-1/2, \sqrt{3}/2, 0)$ , and $\mathbf{r}_2 = a(0, 0, 1)$ . We set the lattice constant to $a = 2.5 \textrm{\AA}$ , 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 $\mathbf{r}_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

$\begin{aligned} \mathbf{r}_{AB}^{(0)} =& \frac{\mathbf{r}_0 + 2\mathbf{r}_1}{3},\\ \mathbf{r}_{AB}^{(1)} =& -\mathbf{r}_1 + \mathbf{r}_{AB}^{(0)},\\ \mathbf{r}_{AB}^{(2)} =& -\mathbf{r}_0 - \mathbf{r}_1 + \mathbf{r}_{AB}^{(0)}, \end{aligned}$

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

Brillouin zone

The reciprocal lattice vectors can now be calculated using

$\begin{aligned} \mathbf{k}_0 =& \frac{2\pi \mathbf{r}_1\times\mathbf{r}_2}{\mathbf{r}_0\cdot\mathbf{r}_1\times\mathbf{r}_2},\\ \mathbf{k}_1 =& \frac{2\pi \mathbf{r}_2\times\mathbf{r}_0}{\mathbf{r}_1\cdot\mathbf{r}_2\times\mathbf{r}_0},\\ \mathbf{k}_2 =& \frac{2\pi\mathbf{r}_0\times\mathbf{r}_1}{\mathbf{r}_2\cdot\mathbf{r}_0\times\mathbf{r}_1}. \end{aligned}$

Here $\mathbf{k}_2$ just like $\mathbf{r}_2$ is related to the z-direction and will be dropped from here on. The two remaining reciprocal lattice vectors $\mathbf{k}_0$ and $\mathbf{k}_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 $\Gamma \rightarrow M \rightarrow K \rightarrow \Gamma$ along which the band structure will be calculated. The high symmetry points can be verified to be given by

$\begin{aligned} \mathbf{\Gamma} =& (0, 0, 0),\\ \mathbf{M} =& (\frac{\pi}{a}, \frac{-\pi}{\sqrt{3}a}, 0),\\ \mathbf{K} =& (\frac{4\pi}{3a}, 0, 0). \end{aligned}$

Hamiltonian

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

$H = -t\sum_{\mathbf{i}\delta}c_{\mathbf{i},A}^{\dagger}c_{\mathbf{i}+\boldsymbol{\delta},B} + H.c.$

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 $\mathbf{i}$ runs over the unit cells, while $\boldsymbol{\delta}$ takes on the three values that makes $\mathbf{i}+\boldsymbol{\delta},B$ the nearest neighbor of $\mathbf{i},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

$\begin{aligned} c_{\mathbf{i},A} =& \frac{1}{\sqrt{N}}\sum_{\mathbf{k}}c_{\mathbf{k},A}e^{ik\cdot \mathbf{R}_{\mathbf{i}}},\\ c_{\mathbf{i}+\boldsymbol\delta},B} =& \frac{1}{\sqrt{N}}\sum_{\mathbf{k}}c_{\mathbf{k},B}e^{i\mathbf{k}\cdot \left(\mathbf{R}_{\mathbf{i}} + \mathbf{r}_{\boldsymbol\delta}\right)}. \end{aligned}$

Here $\mathbf{R}_{\mathbf{i}}$ is the position of the A atom in unit cell $\mathbf{i}$ , while $\mathbf{r}_{\boldsymbol{\delta}} \in \{\mathbf{r}_{AB}^{(0)}, \mathbf{r}_{AB}^{(1)}, \mathbf{r}_{AB}^{(2)}\}$ .¹

Inserting this into the expression above we get

$\begin{aligned} H =& \frac{-t}{N}\sum_{i\mathbf{k}\mathbf{k}'\boldsymbol{\delta}}c_{\mathbf{k},A}^{\dagger}c_{\mathbf{k}',B}e^{-i\mathbf{k}\cdot \mathbf{R}_{\mathbf{i}}}e^{i\mathbf{k}'\cdot\left(\mathbf{R}_{\mathbf{i}} + \mathbf{r}_{\boldsymbol{\delta}}\right)} + H.c\\ =& -t\sum_{\mathbf{k}\boldsymbol{\delta}}c_{\mathbf{k},A}^{\dagger}c_{\mathbf{k},B}e^{i\mathbf{k}\cdot \mathbf{r}_{\boldsymbol{\delta}}} + H.c\\ =& -t\sum_{\mathbf{k}}c_{\mathbf{k},A}^{\dagger}c_{\mathbf{k},B}\left(e^{i\mathbf{k}\cdot\mathbf{r}_{AB}^{(0)}} + e^{i\mathbf{k}\cdot\mathbf{r}_{AB}^{(1)}} + e^{i\mathbf{k}\cdot\mathbf{r}_{AB}^{(2)}}\right) + H.c. \end{aligned}$

In the second line we have summed over $\mathbf{i}$ and used the delta function property of $\sum_{\mathbf{i}}e^{i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}_{\mathbf{i}}}/N$ to eliminat $\mathbf{k}'$ , while in the last line we have made the sum over $\boldsymbol{\delta}$ explicit.

Implementation

Parameters

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

//Set the natural units for this calculation.
UnitHandler::setScales({
        "1 rad", "1 C", "1 pcs", "1 eV",
        "1 Ao", "1 K", "1 s"
});

//Define parameters.
double t = 3;   //eV
double a = 2.5; //Ångström
unsigned int BRILLOUIN_ZONE_RESOLUTION = 1000;
vector<unsigned int> numMeshPoints = {
    BRILLOUIN_ZONE_RESOLUTION,
    BRILLOUIN_ZONE_RESOLUTION
};
const int K_POINTS_PER_PATH = 100;
const double ENERGY_LOWER_BOUND = -10;
const double ENERGY_UPPER_BOUND = 10;
const int ENERGY_RESOLUTION = 1000;

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 $\Gamma \rightarrow M$ , $M \rightarrow K$ , and $K \rightarrow \Gamma$ when calculating the band structure. Finally, the energy window that is used to calculate the DOS is specified in the three last lines.

Real and reciprocal lattice

The lattice vectors $\mathbf{r}_{0}$ , $\mathbf{r}_{1}$ , and $\mathbf{r}_{2}$ are now defined as follows.

Vector3d r[3];
r[0] = Vector3d({a,     0,              0});
r[1] = Vector3d({-a/2,  a*sqrt(3)/2,    0});
r[2] = Vector3d({0,     0,              a});

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

Vector3d r_AB[3];
r_AB[0] = (r[0] + 2*r[1])/3.;
r_AB[1] = -r[1] + r_AB[0];
r_AB[2] = -r[0] - r[1] + r_AB[0];

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

Vector3d k[3];
for(unsigned int n = 0; n < 3; n++){
    k[n] = 2*M_PI*r[(n+1)%3]*r[(n+2)%3]/(
        Vector3d::dotProduct(
            r[n],
            r[(n+1)%3]*r[(n+2)%3]
        )
    );
}

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.

BrillouinZone brillouinZone(
    {
        {k[0].x, k[0].y},
        {k[1].x, k[1].y}
    },
    SpacePartition::MeshType::Nodal
);

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.

vector<vector<double>> mesh
    = brillouinZone.getMinorMesh(
        numMeshPoints
    );

Set up the model

The model is set up using the following code.

Model model;
for(unsigned int n = 0; n < mesh.size(); n++){
    //Get Index representation of the current
    //k-point.
    Index kIndex = brillouinZone.getMinorCellIndex(
        mesh[n],
        numMeshPoints
    );

    //Calculate the matrix element.
    Vector3d k({mesh[n][0], mesh[n][1], 0});
    complex<double> h_01 = -t*(
        exp(
            -i*Vector3d::dotProduct(
                k,
                r_AB[0]
            )
        )
        + exp(
            -i*Vector3d::dotProduct(
                k,
                r_AB[1]
            )
        )
        + exp(
            -i*Vector3d::dotProduct(
                k,
                r_AB[2]
            )
        )
    );

    //Add the matrix element to the model.
    model << HoppingAmplitude(
        h_01,
        {kIndex[0], kIndex[1], 0},
        {kIndex[0], kIndex[1], 1}
    ) + HC;
}
model.construct();

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 $\mathbf{k}_{\mathbf{i}}$ , where $\mathbf{k}_{\mathbf{i}}$ is a real-valued vector while $\mathbf{i}$ 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.

Selecting solver

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

Solver::BlockDiagonalizer solver;
solver.setModel(model);
solver.run();

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

PropertyExtractor::BlockDiagonalizer
    propertyExtractor(solver);
propertyExtractor.setEnergyWindow(
    ENERGY_LOWER_BOUND,
    ENERGY_UPPER_BOUND,
    ENERGY_RESOLUTION
);

Extracting the DOS

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

//Calculate the density of states (DOS).
Property::DOS dos
    = propertyExtractor.calculateDOS();

//Smooth the DOS.
const double SMOOTHING_SIGMA = 0.03;
const unsigned int SMOOTHING_WINDOW = 51;
dos = Smooth::gaussian(
    dos,
    SMOOTHING_SIGMA,
    SMOOTHING_WINDOW
);

//Plot the DOS.
Plotter plotter;
plotter.plot(dos);
plotter.save("figures/DOS.png");

Extract the band structure

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

Vector3d Gamma({0,        0,                 0});
Vector3d M({M_PI/a,       -M_PI/(sqrt(3)*a), 0});
Vector3d K({4*M_PI/(3*a), 0,                 0});

These are then packed in pairs to define the three paths $\Gamma \rightarrow M$ , $M \rightarrow K$ , and $K \rightarrow \Gamma$

vector<vector<Vector3d>> paths = {
    {Gamma, M},
    {M, K},
    {K, Gamma}
};

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

Array<double> bandStructure(
    {2, 3*K_POINTS_PER_PATH},
    0
);
Range interpolator(0, 1, K_POINTS_PER_PATH);
for(unsigned int p = 0; p < 3; p++){
    //Select the start and endpoints for the
    //current path.
    Vector3d startPoint = paths[p][0];
    Vector3d endPoint = paths[p][1];

    //Loop over a single path.
    for(
        unsigned int n = 0;
        n < K_POINTS_PER_PATH;
        n++
    ){
        //Interpolate between the paths start
        //and end point.
        Vector3d k = (
            interpolator[n]*endPoint
            + (1 - interpolator[n])*startPoint
        );

        //Get the Index representation of the
        //current k-point.
        Index kIndex = brillouinZone.getMinorCellIndex(
            {k.x, k.y},
            numMeshPoints
        );

        //Extract the eigenvalues for the
        //current k-point.
        bandStructure[
            {0, n+p*K_POINTS_PER_PATH}
        ] = propertyExtractor.getEigenValue(
            kIndex,
            0
        );
        bandStructure[
            {1, n+p*K_POINTS_PER_PATH}
        ] = propertyExtractor.getEigenValue(
            kIndex,
            1
        );
    }
}

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.

double min = bandStructure[{0, 0}];
double max = bandStructure[{1, 0}];
for(
    unsigned int n = 0;
    n < 3*K_POINTS_PER_PATH; n++ ){ if(min > bandStructure[{0, n}])
        min = bandStructure[{0, n}];
    if(max < bandStructure[{1, n}])
        max = bandStructure[{1, n}];
}

Finally, we plot the band structure.

plotter.clear();
plotter.setLabelX("k");
plotter.setLabelY("Energy");
plotter.plot(
    bandStructure.getSlice({0, _a_}),
    {{"color", "black"}}
);
plotter.plot(
    bandStructure.getSlice({1, _a_}),
    {{"color", "black"}}
);
plotter.plot(
    {K_POINTS_PER_PATH, K_POINTS_PER_PATH},
    {min, max},
    {{"color", "black"}}
);
plotter.plot(
    {2*K_POINTS_PER_PATH, 2*K_POINTS_PER_PATH},
    {min, max}
    {{"color", "black"}}
);
plotter.save("figures/BandStructure.png ");

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’.

Results

Below the results are shown. The DOS is characteristically v-shaped around $E=0$ and is peaked 1/3 of the way to the band edge. In the band structure, we can see the Dirac cone at the K-point (the second vertical line). The low density of states at $E=0$ can be understood as a consequence of the small number of k-points around the K-point for which the energy is close to zero. The cone has a finite slope, and therefore only a single k-point is at $E=0$ (one more point is at the $K'$ -point, which we have not considered here but which is at another corner of the Brillouin zone). Compare this with the flat region at the M-point (first vertical line), which results in a significant amount of k-points with roughly the same energy. This flat region is responsible for the strong peaks in the DOS.

DOS

Band structure

Full code

The full code is available in src/main.cpp in the project 2019_07_05 of the Second Tech code package. See the README for instructions on how to build and run.

Note that the definition of $c_{\mathbf{i}+\boldsymbol{\delta},B}$ is threefold redundant since for each site B in unit cell $\mathbf{i}$ there are three ways to label the same site. For each of the three vectors $\mathbf{r}_{AB}^{(0)}$ , $\mathbf{r}_{AB}^{(1)}$ , and $\mathbf{r}_{AB}^{(2)}$ there is a different unit cell from which the B site in cell $\mathbf{i}$ can be reached. However, each of these three redundant definitions can be verified to be equivalent and the redundant expression simplifies the mathematical notation.