Questions about Creation/Annhiliation Operators of electrons I am trying to solve a question that asks me to find the two energy bands of a NaCl crystal, where "c" is the distance between Na and Cl where $a_{\vec{R}}^{\dagger}$ creates and electron centered on Cl at location $\vec{R}$, and $b_{\vec{R}+\vec{\delta_i}}^{\dagger}$ creates and electron centered on Cl at location $\vec{R}+\vec{\delta_i}$ ($\delta_i$ are the six nearest neighbors). 
It is given that the Hamiltonian is $H=-t\Sigma_{\vec{R}+\vec{\delta_i}}(a_{\vec{R}}^{\dagger}b_{\vec{R}+\vec{\delta_i}}+b_{\vec{R}+\vec{\delta_i}}^{\dagger}a_{\vec{R}})$. 
To solve this, I Fourier transformed the Hamiltonian by using 
$a_{\vec{R}}=\frac{1}{\sqrt{N}}\Sigma_{k}e^{-i\vec{k}\cdot\vec{r}}a_{\vec{k}}$, and $b_{\vec{R}+\vec{\delta_i}}=\frac{1}{\sqrt{N}}\Sigma_{k}e^{-i\vec{k}\cdot(\vec{r}+\vec{\delta_i})}b_{\vec{k}}$.
I got some Hamiltonian expressed in terms of a,b,a dagger and b dagger. The professor told me to restructure the Hamiltonian in the form, $H=\Sigma_{\vec{k}}(a_{\vec{k}}^{\dagger}~~ b_{\vec{k}}^{\dagger})M\begin{pmatrix}a_{\vec{k}}\\a_{\vec{k}}   \end{pmatrix}$
Then said that if we take the eigenvalues of M, it would be the energies.
However, there are few things that I am not clear of :


*

*Why did we change the hamiltonian to the form $H=\Sigma_{\vec{k}}(a_{\vec{k}}^{\dagger}~~ b_{\vec{k}}^{\dagger})M\begin{pmatrix}a_{\vec{k}}\\a_{\vec{k}}   \end{pmatrix}$? H was initially a 'scalar' of a, b, and so on. Why did we try to 'matrixize' it? What does it mean to go from H to M? (i.e. from Scalar to Matrix form?)

*Why is the Eigenvalue of M equal to the eigenvalue of H? Is it because the Eigen values of 'Fourier transformed' H is different from the eigenvalues of actual H, and the actual eigenvalues of H is the eigenvalues of M? What happens to the eigenvalues of a matrix when it is Fourier Transformed? Does it stay the same?
Thank you for reading these questions:)
 A: I assume that in the matrix formulation, the vector to the right of the matrix is $\begin{pmatrix} a_{\vec{k}} \\ b_{\vec{k}} \end{pmatrix}$ and not with $a_{\vec{k}}$ in the lower entry as you wrote.
The thing about writing the Hamiltonian in that way is that it is very convenient, and makes it straight-forward, to understand its structure and solutions. Similarly to Fourier Transforming, by changing the basis in which we represent the modes will allow us to derive its solution and find the eigenmodes. The matrix representation helps us in doing so.
Let's say we have an operator $$A = \sum_n (E^{a}_n a^{\dagger}_n a_n + E^{b}_n b^{\dagger}_n b_n) + \sum_n (\lambda_n a^{\dagger}_n b_n + \rm{h.c.})$$ and I want to find its eigenmodes, where $a_n$ and $b_n$ are cannonical creation and annihilation operators (either bosonic or fermionic). Since this is a quadratic oeprator, the eigenmodes will be linear in the creation and annihilation operators, and we already see that the operator doesn't mix between different $n$'s. So we can write the solution generally as $c_n = \alpha_n a_n + \beta_n b_n$. Now we want to solve the eigenmode equation $[c_n , A] = \epsilon_n c_n$. Writing it explicitly we get two sets of equations - one equating the $a_n$ terms on the left-hand-side and the right-hand-side, and one equating the $b_n$ terms. As this is a set of two linear equations in $\alpha_n$ and $\beta_n$ you can already sense that solving it using matrix is much easier, and indeed you can put this in matrix form and solve.
Things, however, can be made even nicer and more general! The fact that $c_n$ is an eigenmode of $A$ means that after we write everything in terms of $c_n$ the operator will be of the form $A = \sum_n \epsilon_n c^{\dagger}_n c_n$. This is just a basis change to the basis of the eigenmodes. And $\alpha_n$ and $\beta_n$ represent the basis change coefficients. We demand that $\{ c^{\dagger}_n, c_n \} = 1$, and if we take the basis change to be done by a unitary matrix $U_n$ is immediately maintained! So we see that now we have reduced the problem to finding a way in which we can transform 
$$ A = \sum_n \begin{pmatrix} a^{\dagger}_{n} & b^{\dagger}_{n} \end{pmatrix} M \begin{pmatrix} a_{n} \\ b_{n} \end{pmatrix}$$
which is a "scalar" (in the space of $a_n$ and $b_n$) into diagonal form, by pushing an appropriate $U_n$ that will transform its basis. This is just a $2\times 2$ diagonalization problem of the matrix $A$, with the eigenvalues giving us $\epsilon_n$ and the rotation matrix $U_n$ giving us $c_n$ as a linear superposition of $a_n$ and $b_n$. We are guaranteed that when we do that, we will get $A$ in a diagonal form $A = \sum_n \epsilon_n c^{\dagger}_n c_n$ as we wants it to be, with $\{c^{\dagger}_m, c_n\}=\delta_{m,n}$ the new basis of fermionic creation and annihilation operators.
