I want to study the surface states in a material using the tight binding model (the goal is to find surface states inside the bulk band gap). The material in question has rock salt crystal structure and I am interested in the (001)-surface.

Now, I have the Hamiltonian matrix, $H_0^{b}$, for the bulk written in a basis using the s, p and d orbitals for the respective elements in the crystal, so $H_0^b$ is a $36\times36$ matrix, when accounting for spin, and can be expressed as $$H_0^b = \begin{pmatrix}H_{11} & H_{12}\\ H_{21} & H_{22}\end{pmatrix},$$ where 1 and 2 label the different atoms in the unit cell.

Now I want to use "the slab method" to find the surface states. I have not found any detailed description about this anywhere, but the principle is that you model the material as a slab with finite thickness that is infinite in two dimensions. You then create a big unit cell and model the slab as two-dimensional. In the case of the rock-salt structure mentioned above with a (001)-surface, the Hamiltonian matrix for the slab would be a $36n\times36n$ matrix, where $n$ is the number of atomic layers in the slab. The structure of the Hamiltonian should be the following: $$H = \begin{pmatrix}H_0^s &H_I& 0& 0& ...\\ H_I^{\dagger} & H_0^s & H_I& 0 & 0& ... \\ 0 & H_I^{\dagger} & H_0^s &H_I&0&...\\ \vdots& 0&\ddots&\ddots&\ddots\end{pmatrix},$$ where $H_0^s$ is the Hamiltoninan for an isolated atomic layer of the slab, i.e. in this case a (001)-layer, and $H_I$ describes the interaction between the layers. I know that I should be able to construct H using the tight binding parameters for the bulk (at least to a good enough approximation), so the question is how to relate the matrices $H_0^s$ and $H_I$ to $H_0^b$. The problem with $H_0^b$ is that it contains terms that have elements containing $k_z$. Discussing this with some people I came to the conclusion that since the position of atom 2 relative to atom 1 in the unit cell is $\frac{a}{2}(1,0,0)$ and $\frac{a}{2}(1,0)$ in the 3D and 2D cases respectively, one should be able to get $H_0^s$ immediately from $H_0^b$ simply by removing all terms containing $k_z$. Then I also find that $$H_I = \begin{pmatrix}0 & H_{12}^I\\ H_{21}^I & 0\end{pmatrix},$$ where $H_{12}^I$ and $H_{21}^I$ are obtained from $H_{12}^b$ and $H_{21}^b$ by removing all terms that contain $k_x$ or $k_y$ (since we are only interested in the interaction between the single atomic layers in the z-direction).

When I implement this in MATLAB and calculate the eigenvalues of $H$ I do not get the expected results and I am wondering if I have misunderstood how to construct $H$. So any comments about my construction of $H$ or suggested references that describe tight binding for slabs in detail are welcome!

  • $\begingroup$ Been a while since I did this but I think the trick is to realize your unit cell is discretely transitionally invariant along the x and y directions. This allows you to apply Bloch's theorem and get a Kx and Ky quantum number and eigenstates and energies in terms of the unit cell Hamiltonian. Then, diagonalizing that Hamiltonian allows you to derive the band structure. I remember one small detail that can trip you up is relating the size of the unit cell to the size of the atomic spacing. Take care with that conversion. It might be why your matlab result is not correct. $\endgroup$ Commented Jun 14, 2017 at 21:44
  • $\begingroup$ Perhaps you could start with the 1D case with a two atom basis. I found a link on that problem that might be helpful. lampx.tugraz.at/~hadley/ss1/bands/tightbinding/tightbinding.php $\endgroup$ Commented Jun 14, 2017 at 21:46
  • $\begingroup$ @GregPetersen And with relating the size of the unit cell to the size of the atomic spacing you mean that the lattice constant is twice that of the atomic spacing? Because I'm quite sure I have taken that into account. $\endgroup$
    – ECA18
    Commented Jun 15, 2017 at 10:35
  • $\begingroup$ (Accidentally posted my last comment too early.) @GregPetersen Yes, but this only works when you have one atomic layer. When you have a slab, you need some kind of kz-dependence so that the layers can interact with each other. (But then you must put kz = 0 in the end since you are interested in the surface.) Thanks for the link! $\endgroup$
    – ECA18
    Commented Jun 15, 2017 at 10:41
  • $\begingroup$ In the case of the slab, (taking x and y to be the infinite directions and z to be the finite), your quantum number kz is determined from the diagonalization of your unit cell Hamiltonian. I'll see if I can recall some more specific details. $\endgroup$ Commented Jun 15, 2017 at 14:27

1 Answer 1


To solve this kind of issue, you need to apply a combination of Bloch's theorem and numerical diagonalization of the unit cell matrix. suppose you have a two-dimensional system with unit cell that forms a square lattice which repeats in the X and Y direction but not in the Z. The basis states can be represented by the states $$ |b: x, y\rangle $$ where b is the index of an atom in the unit cell and x and y are the positions of the unit cell in the XY plane. We can define an operator $\hat{B}$ where $$ $\hat{B}|b: x, y\rangle = -t|b': x + \Delta x, y + \Delta y\rangle $$ Note: t can be a function of $ \Delta x$, $\Delta y$, $b$ and $b'$, however, as to not overcrowd the notation we will omit them for the time being. We can also write B in terms of a creation/annialition operator pair as $$ \hat{B} = -t\hat{b'}_{x + \Delta x, y + \Delta y}^\dagger\hat{b}_{x,y} $$ These \hat{B} operators essentially make up the entire tight-binding Hamiltonian that we know and love. One can write any such Hamiltonian as a sum over these values $$ H = \sum_{x,y,b,b'} (-t\hat{b'}_{x + \Delta x, y + \Delta y}^\dagger\hat{b}_{x,y} + h.c.) $$ Up until this point we have simply defined the tight-binding Hamiltonian. Now let us begin to solve it. First, let us use our knowledge of the discrete transnational symmetry of the system in the X and Y directions. Let's start by taking the Fourier transform of the creation and annihilation operators in x and y. $$ \hat{b}_{x,y} = \sum_{k,q} \hat{b}_{k,q}e^{-ikx}e^{-iqy} $$ $$ \hat{b}^\dagger_{x,y} = \sum_{k',q'} \hat{b}^\dagger_{k',q'}e^{ik'x}e^{iq'y} $$ (Note: for non-square lattices you will need the dot product of the translation vector with the momentum vector but for this example we limit ourselves to a square lattice).

Inserting this equation back into the hamiltonian we have $$ H = \sum_{x,y,b,b'}-t\sum_{k,q} \hat{b}_{k,q}e^{-ikx}e^{iqy}\sum_{k',q'} \hat{b}^\dagger_{k',q'}e^{-ik'x + \Delta x}e^{iq'y + \Delta y} + h.c. $$

And rearranging $$ H = -\sum_{k,k',q,q',b,b'} t \hat{b}_{k,q}\hat{b}^\dagger_{k',q'}e^{ik' \Delta x}e^{iq' \Delta y} \sum_{x,y}e^{ix(k - k')}e^{iy(q - q)} + h.c. $$

We can then apply periodic boundary conditions to x and y which puts creates restriction

$$ k = 2\pi n / Na $$

$$ q = 2\pi m / Na $$

where as is the lattice constant, n and m are integers, and N are the number of repeated unit cells in the x and y direction (You will need to numerically compute the scaling laws later on if you want to simulate an infinite system. Check out Born-Von Karman boundary conditions for more information).

applying this leads us to a Dirac-delta function which gives k = k' and q = q' leading to

$$ H = -\sum_{k,q,b,b'} t \hat{b}_{k,q}\hat{b}^\dagger_{k,q}e^{ik \Delta x}e^{iq \Delta y} + h.c. $$

We now have a Hamiltonian that is diagonal in k and q, but not in b. The final step is to diagonalize our Hamiltonian in B space. However, there is now straightforward method to do this as we don't have any information on the symmetry of the system. Thus, most people resort to a numerical approach. The psuedo-code would be as follows.

  1. Set a value for k and q
  2. Diagonalize hamiltonian in B-space
  3. Save E(k,q)
  4. incremenent k,q

This will ultimately lead to the band structure.

Example: 2-atom basis in a 1D system.

Assuming an equal hopping energy, the Hamiltonian can be written as

$$ H = -t\sum_{x} (\hat{b_1}_{x + a}^\dagger\hat{b_2}_{x} + \hat{b_1}_{x}^\dagger\hat{b_2}_{x} + \hat{b_1}_{x}^\dagger\hat{b_2}_{x - a} + h.c.) $$

which can be written in Fourier space as $$ H = -t\sum_{x} (1 + e^{-iak} - e^{iak}) \hat{b_1}^\dagger\hat{b_2} + h.c. $$

You then just need to diagonalize a 2x2 Hamiltonian over all values of k.

  • $\begingroup$ Thank you for your detailed answer! It provided a lot of insight since it uses a different approach than I'm used to. Would have upvoted if I could (and I don't want to mark it as answered yet since there might be more to say about it). $\endgroup$
    – ECA18
    Commented Jul 5, 2017 at 10:02
  • $\begingroup$ Glad this helped. If a better answer exists I would be interested to learn it as well so please let us know what you come up with! $\endgroup$ Commented Jul 8, 2017 at 2:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.