Diagonalizing a given Hamiltonian The following Hamiltonian, which has to be diagonalized, is given:

$H = \epsilon(f^{\dagger}_1f_1 + f_2^{\dagger}f_2)+\lambda(f_1^{\dagger}f_2^{\dagger}+f_1f_2)$

$f_i^{\dagger}$ and $f_i$ represent fermionic creation and annihiliation operators.
Right now I am not sure how to approach this problem. My idea is to use some kind of Bogoliubov transformation. I would be thankful for ideas on how to approach this problem.
 A: Since we have only two fermions creation operators, we are dealing with a finite dimensional system. In that case, I find that it is often easier to write out matrices and do algebra on those.
In the basis $(|00\rangle,|01\rangle,|10\rangle,|11\rangle)$ we have :
\begin{align}
f_1 &= \begin{pmatrix} 0 & 0 & 1 & 0\\
0&0&0&1\\
0&0&0&0\\
0&0&0&0\\
\end{pmatrix}\\
f_2 &= \begin{pmatrix} 0&1&0&0\\
0&0&0&0 \\
0&0&0&1\\
0&0&0&0
\end{pmatrix}
\end{align}
Therefore :
$$H = \begin{pmatrix} 
0 & 0&0 & \lambda\\
0 & \epsilon & 0& 0 \\
0 & 0 & \epsilon &0 \\
\lambda & 0& 0 & 2\epsilon
\end{pmatrix}$$
The states $|01\rangle$ and $|10\rangle$ are eigenstates with eigenvalue $\epsilon$.
On the orthogonal subspace, generated by $|00\rangle$ and $|11\rangle$, the induced Hamiltonian is $\begin{pmatrix} 0 & \lambda \\ \lambda & 2\epsilon \end{pmatrix} = \epsilon \mathbb I_2 + \lambda \sigma_x - \epsilon \sigma_z$. Therefore the eigenvalues are $\epsilon \pm \sqrt{\lambda^2 + \epsilon^2}$ and the eigenvectors are :
$$|+\rangle = \frac{1}{\sqrt{2(\lambda^2+\epsilon^2 -\epsilon \sqrt{\epsilon^2 +\lambda^2}) }}\begin{pmatrix}-\epsilon + \sqrt{\epsilon^2 + \lambda^2} \\ \lambda\end{pmatrix} $$
and
$$|-\rangle = \frac{1}{\sqrt{2(\lambda^2+\epsilon^2 +\epsilon \sqrt{\epsilon^2 +\lambda^2}) }}\begin{pmatrix}-\epsilon - \sqrt{\epsilon^2 + \lambda^2} \\ \lambda\end{pmatrix} $$
In the basis $|01\rangle,|10\rangle,|+\rangle,|-\rangle$, we have :
$$H = \begin{pmatrix} \epsilon & 0 & 0 & 0\\ 0 & \epsilon &0 & 0\\
0& 0 & \epsilon+\sqrt{\lambda^2+\epsilon^2} & 0 \\ 0 & 0  & 0& \epsilon -\sqrt{\lambda^2 + \epsilon^2}\end{pmatrix}$$
To write this as a Bogoliubov transform, we remark that we can write $H = E_0 + E_1 ( c_1^\dagger c_1 + c_2^\dagger c_2)$, with $c_1$ and $c_2$ independent fermion annihilation operators, when :
\begin{align}
E_0 &= \epsilon \\
E_1 &= \sqrt{\lambda^2 + \epsilon^2}
\end{align}
and the eigenstates of $c_1^\dagger c_1$ and $c_2^\dagger c_2$ are :
\begin{align}
|00\rangle' &= |-\rangle \\
|01\rangle' &= |01\rangle\\
|10\rangle' &= |10\rangle\\
|11\rangle' &= |+\rangle
\end{align}
Solving for $c_1,c_2$ as linear combination of $f_1,f_2,f_1^\dagger,f_2^\dagger$ (as in a Bogoliubov transform), we get :
\begin{align}
c_1 &=\frac{1}{\sqrt{2(\lambda^2+\epsilon^2 +\epsilon \sqrt{\epsilon^2 +\lambda^2}) }} \left( (-\epsilon -\sqrt{\epsilon^2+ \lambda^2})f_1 + \lambda f_2^\dagger\right) \\ &\qquad+\frac{1}{\sqrt{2(\epsilon^2 + \lambda^2 - \epsilon\sqrt{\epsilon^2+ \lambda^2})}} \left( (-\epsilon +\sqrt{\epsilon^2+ \lambda^2})f_2 + \lambda f_1^\dagger\right) \\
c_2 &=\frac{1}{\sqrt{2(\lambda^2+\epsilon^2 -\epsilon \sqrt{\epsilon^2 +\lambda^2}) }} \left( (-\epsilon + \sqrt{\epsilon^2+ \lambda^2})f_1 + \lambda f_2^\dagger\right) \\ &\qquad+\frac{1}{\sqrt{2(\epsilon^2 + \lambda^2 + \epsilon\sqrt{\epsilon^2+ \lambda^2})}} \left( (-\epsilon -\sqrt{\epsilon^2+ \lambda^2})f_2 + \lambda f_1^\dagger\right) \\
\end{align}
There might be smarter/more efficient ways do perform the calculations, but this does the job.
