Second quantization and Hamiltonian diagonalization So I want to diagonalize my Hamiltonian (it is bosonic hamiltonian) which is: 
$H=(E+\Delta)a^{\dagger}a + 1/2\Delta(a^{\dagger}a^{\dagger} + aa)$
My class didn't cover this material so I don't really know how to procede. I would be grateful for any literature which covers this topics and a problem book with solutions would be great too. 
What I tried to do was writing my Hamiltonian in matrix form which would be: $\begin{pmatrix}
  1/2 \Delta & 1/2(E+\Delta)  \\
  1/2(E+\Delta) & 1/2 \Delta   \\
 \end{pmatrix}$
And then diagonalize it, find eigenstates etc. Is this the correct way?
 A: Diagonalizing the Hamiltonian means you want to bring it into the form $H=\omega b^\dagger b$, and it is pretty obvious that $b$ should be a linear combination of $a$ and $a^\dagger$, and $b$ should satisfy the canonical commutation of annihilation operators, namely $[b,b^\dagger]=1, [b,b]=0$. 
Now let's write $b=ua+va^\dagger$ (this is called the Bogoliubov transformation, by the way). The condition $[b,b^\dagger]=1$ leads to $|u|^2-|v|^2=1$. Let us expand out $b^\dagger b$:
$$
b^\dagger b= |u|^2 a^\dagger a+ |v|^2 a a^\dagger + u^*v a^\dagger a^\dagger + uv^* aa.
$$
Therefore
$$
\omega(|u|^2+|v|^2)=E+\Delta, \omega u^*v = \frac{1}{2}\Delta.
$$
Together with $|u|^2-|v|^2=1$, we have three equations for three variables ($u, v, \omega$). In fact, in this case one can safely assume $u$ and $v$ are both real. The rest is just algebra.
A: Diagonalising an operator means finding its eigenstates.
Without loss of generality your Hamiltonian can be written as 
$$
H = c_1 a^{\dagger}a + c_2 a^{\dagger}a^{\dagger} + c_3 a a
$$
with $a^{\dagger},a$ being operators of the type $a^{\dagger}\colon \mathcal{H}_n\mapsto \mathcal{H}_{n+1}$ (and conversely for $a$), where $\mathcal{H}_n$ is the $n$-particle Hilbert space contributing to the Fock space $\mathcal{F}= \oplus^{\infty}_n\mathcal{H}_n$.
There must be a few errors in your equation if you really mean that in a second quantisation procedure. First of all there is no general $a^{\dagger},a$ operator, rather you have one for each momentum $k$, that is $a^{\dagger}_k,a_k$ create and destroy (in quotation marks) particles with momentum $k$; there is no $k$ in your initial Hamiltonian, whereas the general form must be $\sum_k c_k\,a^{\dagger}_ka_k$.
Second of all: according to whether your particles are fermions or bosons the corresponding operators behave in a different way: for instance $a^{\dagger}_ka^{\dagger}_k=0$ for fermions.
If the Hamiltonian acts on a subspace of the Fock space with a certain number of particles $\mathcal{H}_n$, then the last two terms in your equation would bring the action onto $\mathcal{H}_{n\pm2}$, therefore the rhs will live in $\mathcal{H}_n +\mathcal{H}_{n+2} +\mathcal{H}_{n-2}$, which does not really make any sense since no prescription on how to sum elements in different Hilbert spaces is given (the last two pieces).
Either you assign a precise prescription to achieve the above, or there must be errors elsewhere in the formula, as pointed out; try giving more context so that one can work out what you mean. This said, suggested literature on how to write any Hamiltonian in second quantisation and find the corresponding solutions is, for example:


*

*http://www.springer.com/us/book/9783540850618

*http://www.amazon.com/Quantum-Field-Theory-Franz-Mandl/dp/0471496847
A: How about just using Matrix representation.
https://en.wikipedia.org/wiki/Creation_and_annihilation_operators#Matrix_representation
You can have any number of Bosons from 0 to infinity. That will be your basis, and your wave function is represented as a vector in which element 0 gives probability amplitude of having 0 quanta, element 1 gives amplitude for 1 quantum, 2 for 2 quanta etc. in the system.
Calculating with matrices is easy:
N = 1000;
a = zeros(N);
for i=1:N-1
a(i,i+1) = sqrt(i);
end
H = 10*a*a' + 5 / 2 * (a*a+a'*a');
eig(N)

Disclaimer: I have worked with Fermions almost always, except for some quantum course 8 years ago.
