I think I have a conceptual doubt on the way you calculate the matrix elements in momentum space for the tight binding Hamiltonian of graphene. I will break down my question into sections to make it as organized and self contained as possible.
The general background:
The electronic Hamiltonian operator for the system in first quantization is
$\hat{H}=-t\sum_{i,j} (|i\rangle \langle j|+|j\rangle \langle i|)$
where $i$ and $j$ are nearest neighbors, $t$ hopping parameter. As usual we can rewrite this as:
$\hat{H}=-t\sum_{i \in A,\delta} (|i\rangle \langle i+\delta|+|i+\delta\rangle \langle i|)$
where the sum over $\delta$ is carried out over the nearest-neighbor vectors $\delta_{1}$, $\delta_{2}$,$\delta_{3}$ that connect an atom of the $A$ sub lattice to the closest 3 atoms of the $B$ sub lattice in graphene.
Calculation of H matrix
Now let's say I want to calculate the matrix for this hamiltonian operator in momentum space. Since the system is translational invariant, the hamiltonian in momentum space will be block diagonal, with each block corresponding to a value of $k$. So I just need to calculate:
$h(k)=\begin{pmatrix} \langle k_{A}|\hat{H}|k_{A}\rangle &\langle k_{A}|\hat{H}|k_{B}\rangle\\ \langle k_{B}|\hat{H}|k_{A}\rangle & \langle k_{B}|\hat{H}|k_{B}\rangle\end{pmatrix}$
where by $|k_{A}\rangle$ I indicate a state where the electron has momentum k and it is completely delocalized on sub lattice $A$ (so no amplitude on $B$).
Diagonal matrix elements
Now let's look at those matrix elements. What people usually write is:
$\langle k_{A}|\hat{H}|k_{A}\rangle=\langle k_{B}|\hat{H}|k_{B}\rangle=0$
the reason is because $i \in A$ and $i+\delta \in B$, so, for example, in the first case I get:
$\langle k_{A}|\hat{H}|k_{A}\rangle=-t\sum_{i \in A,\delta}(\langle k_{A}|i\rangle\langle i+\delta|k_{A}\rangle + \langle k_{A}|i+\delta \rangle\langle i|k_{A}\rangle)$
And $\langle k_{A}|i+\delta\rangle$ is zero because, as I said, $| k_{A}\rangle$ has no amplitude on $B$ .
So far so good!
Off diagonal matrix elements
But now let's look at the off diagonal matrix elements. The usual calculation that I found in many resources goes as:
$\langle k_{A}|\hat{H}|k_{B}\rangle=-t\sum_{i \in A,\delta}(\langle k_{A}|i\rangle\langle i+\delta|k_{B}\rangle + \langle k_{A}|i+\delta \rangle\langle i|k_{B}\rangle)$
By the same argument as before the second term is zero so I just have:
$\langle k_{A}|\hat{H}|k_{B}\rangle=-t\sum_{i \in A,\delta}\langle k_{A}|i\rangle\langle i+\delta|k_{B}\rangle$
Now $\langle k_{A}|i\rangle$ is nothing but the conjugate of the wave function $\langle i|k_{A}\rangle$, which is the wave function, written in real space, of a state with momentum $k$, delocalized on the sub lattice A. So $\langle k_{A}|i\rangle=\frac{1}{\sqrt{N/2}}e^{-ik r_{i}}$. The same goes for the other factor and at the end I get:
$\langle k_{A}|\hat{H}|k_{B}\rangle=-t\sum_{i \in A,\delta} \frac{1}{N/2}(e^{-ik r_{i}}
e^{ik(r_{i}+\delta)} )=-t\sum_{i \in A,\delta}\frac{1}{N/2} e^{i k \delta} =-t \sum_{\delta} e^{ik \delta}$.
Similarly I get the conjugate of that for the other off diagonal element.
Final matrix
So at the end you get:
$h(k)=\begin{pmatrix} 0 & \sum_{\delta} e^{i k \delta} \\ \sum_{\delta} e^{-i k \delta} & 0 \end{pmatrix}$
which is the usual tight binding hamiltonian in momentum space that you find in the literature for graphene.
Here comes my question:
I find it very strange that the product $\langle k_{A}|i\rangle\langle i+\delta|k_{B}\rangle$ (hence the off diagonal elements of the H matrix) gives a non zero result! Because what we are doing is to multiply two wave functions, where one has an amplitude only on $A$ sub lattice (i.e. $\langle k_{A}|i\rangle$), whereas the other only has amplitude on $B$ sub lattice (i.e. $\langle i+\delta|k_{B}\rangle$). So where one is non-zero the other is zero and viceversa. So why don't I get zero?
I think this problem emerges clearly if you choose to represent your wave function as a two component spinor where the first component corresponds to sub lattice $A$ and the second to $B$:
$\psi(x)=\begin{pmatrix}\psi_{a}(x)\\\psi_{b}(x) \end{pmatrix}$.
So my basis would be a collection of spinors containing delta functions centered on each lattice site:
$\phi_{i,a}=\begin{pmatrix}\delta(x_{i})\\0 \end{pmatrix}$
and
$\phi_{j,b}=\begin{pmatrix}0\\\delta(x_{j}) \end{pmatrix}$
where $i$ goes over all sites of the $A$ sub lattice and $j$ same for $B$.
If you do that and now try to calculate the above off diagonal matrix element you would get:
$\langle k_{A}|i\rangle=\begin{pmatrix} \frac{e^{-ik r_{i}}}{\sqrt{N/2}} \\0 \end{pmatrix}^{T}$,
$\langle i+\delta|k_{B}\rangle=\begin{pmatrix} 0 \\\frac{e^{ik (r_{i}+\delta)}} {\sqrt{N/2}} \end{pmatrix}$
because now the wave functions are spinors. And then
$\langle k_{A}|\hat{H}|k_{B}\rangle=-t\sum_{i \in A,\delta}\langle k_{A}|i\rangle\langle i+\delta|k_{B}\rangle=-t\sum_{i \in A,\delta} \frac{1}{N/2} \begin{pmatrix} e^{-ik r_{i}} & 0 \end{pmatrix} \cdot \begin{pmatrix} 0 \\e^{ik (r_{i}+\delta)} \end{pmatrix}=0$
where I just did the inner product between spinors.
To conclude
Why the value of those off diagonal matrix elements is $\sum_{\delta} e^{\pm ik \delta}$ and not zero? What is wrong in my argument that they are zero?