Tight binding Hamiltonian for a slab 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! 
 A: 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.


*

*Set a value for k and q

*Diagonalize hamiltonian in B-space

*Save E(k,q)

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