Diagonalising the Hamilton operator, why does this magic work? Let the Hamilton operator $H= \omega_1 a_1^\dagger a_1 + \omega_2 a_2^\dagger a_2 + \frac{J}{2} (a_1^\dagger a_2 + a_1 a_2^\dagger)$ be given, of course $a_j$ and $a_j^\dagger$ are the creation and annihilation operators, respectively. This operator can be rewritten as
$$
H=
\begin{pmatrix} a_1^\dagger & a_2^\dagger \end{pmatrix}
\begin{pmatrix}
\omega_1 & J/2\\ J/2& \omega_2
\end{pmatrix}
\begin{pmatrix} a_1 & a_2 \end{pmatrix}
$$
The Eigenvalues of the matrix in the middle, which may be rewritten as
$$
\begin{pmatrix}
\omega_0 - \frac{\Delta}{2\omega_0} & J/2\\ J/2& \omega_0 + \frac{\Delta}{2\omega_0}
\end{pmatrix}
$$
where $\omega_0 = \frac{\omega_1 + \omega_2}{2}$ and $\Delta = \omega_2 - \omega_1$, are
$$
\lambda_{\pm} = \omega_0 \pm \frac{1}{2}\sqrt{J^2 + \Delta^2}.
$$
Question: Could these be the Eigenvalues of that thing $H$ that maps vectors from an infinite dimensional space to others ? Why is that ? How could the diagonalisation of the little $2\times2$-matrix have anything to do with diagonalising $H$, whose true matrix is infinite dimensional ?
Please mention the subject that rigorously explains this and some references about that.
 A: Using the Schwinger realization, this $H$ can be rewritten as
$$
H=A\hat N+B\hat J_z+C\hat J_y+D\hat J_x
$$
where
\begin{align}
\hat N&=\frac{1}{2}\left(\hat a^\dagger_1\hat a_1+
\hat a^\dagger_2\hat a_2\right)\, ,\\
\hat J_z&=\frac{1}{2}\left(\hat a^\dagger_1\hat a_1-
\hat a^\dagger_2\hat a_2\right)\, ,\\
\hat J_+&=\hat a_1^\dagger \hat a_2\, , \\
\hat J_-&=\hat J_+^\dagger\, .
\end{align}
You can verify that $\hat N$ commutes with everything else so that basically inside each subspace of fixed eigenvalue $n$ of $\hat N$ your $H$ is just proportional to a rotation $R$ of the $\hat J_z$ operator.  The change of basis that will give you this rotated $\hat J'_z=R J_z R^{-1}$ is just given by the eigenvectors of your matrix.
My favourite source for this is the (old but nice) textbook by Gordon Baym.
A: I will provide an answer using a technique which (suitably generalized) is very often used in condensed matter physics (e.g. in superconductivity, magnetism, superfluidity). This is what @Mauricio was referring to by the Bogoliubov transformation (although again, that is really a slightly fancier version of the following technique.)
Assume $[a_1,a_2]=[a_1,a_2^\dagger] = [a_1^\dagger,a_2]= [a_1^\dagger,a_2^\dagger] = 0$. (It is not too much trouble to handle the case where they anticommute either.) Let the Hamiltonian be given by
\begin{equation}
H = \mathbf{a}^\dagger h \mathbf{a},
\end{equation}
where $\mathbf{a} = (a_1,a_2)^T$ is the vector of annihilation operators, and $h$ is the $2\times2$ matrix you wrote. We can diagonalize $h$ by writing $h = U^{^\dagger} D U$. Here, $U$ is a unitary change of basis matrix and $D$ is the matrix composed of the eigenvalues $\lambda_1$ and $\lambda_2$ along the diagonal. With this, we then have
\begin{align}
H 
&= \mathbf{a}^\dagger h \mathbf{a}\\
&= \mathbf{a}^\dagger U^{-1} D U \mathbf{a}\\
&= \mathbf{b}^\dagger D \mathbf{b},
\end{align}
where we have defined the new operators $\mathbf{b}= U \mathbf{a}$. Writing out this matrix product to be explicit, we have
\begin{align}
b_1 &= U_{11}a_1 + U_{12}a_2\\
b_2 &= U_{21}a_1 + U_{22}a_2\\
\end{align}
If you expand this product, because $D$ is diagonal, $H$ is simply
\begin{equation}
H = \lambda_1 b_1^\dagger b_1+\lambda_2 b_2^\dagger b_2.
\end{equation}
From this expression, there is a natural interpretation. $b_1^\dagger$, $b_1$, $b_2^\dagger$ and $b_2$ are just the creation and annihilation operators of new excitations (this will be discussed more at the end of the post). $b_1^\dagger b_1$ counts the number of excitations with energy $\lambda_1$, and $b_2^\dagger b_2$ counts the number of excitations with energy $\lambda_2$. From this, it should be clear that the energy of the system is
$$E_{n_1,n_2} = \lambda_1 n_1 + \lambda_2 n_2,$$
where $n_1$ counts the number of $b_1$ excitations and same for $b_2$. Here, the eigenstates of the system are given by $|n_1,n_2\rangle$, which is defined by $b_1^\dagger b_1|n_1,n_2\rangle = n_1 |n_1,n_2\rangle$, and $b_2^\dagger b_2|n_1,n_2\rangle = n_2 |n_1,n_2\rangle$.
At this point, we're basically done! By changing to the eigenbasis of $h$, we have essentially found the eigenstates of $H$.
Now, although we have interpreted the $b^\dagger$ and $b$ operators as creation and annihilation operators, we haven't actually checked whether they satisfy the correct commutation relations: $[b_1,b_1^\dagger] = 1$, $[b_2,b_2^\dagger] = 1$, $[b_1,b_2] = [b_1,b_2^\dagger] =[b_1^\dagger,b_2] =[b_1^\dagger,b_2^\dagger] = 0$. This actually follows from the unitarity of $U$ (and the original commutation relations for the original $a$ operators.) I may edit this post to add the proof of this later, though please feel free to verify this yourself.
Edit: Here, I give a more explicit picture of the eigenstates. Let's talk about the ground state specifically, first. The ground state of this Hamiltonian ($| 0,0\rangle$) is defined as the state with zero $b_1$ or $b_2$ excitations: i.e. $b_1 |\mathrm{gs}\rangle = b_2|\mathrm{gs}\rangle= 0$. Since $b_1$ and $b_2$ are linear combinations of $a_1$ and $a_2$, the ground state with zero $b$ excitations is also the state with zero $a$ excitations. (Note: in the more fancy application of the Bogoliubov transformations, the ground state will be more complicated to think of in terms of the original $a$ operators, although it can still be defined as the state with no $b$ excitations.)
For excited states,
\begin{equation}
| n_1, n_2 \rangle = \frac{(b_1^\dagger)^{n_1}(b_2^\dagger)^{n_2}}{\sqrt{n_1! n_2!}}|\mathrm{gs}\rangle.
\end{equation}
One can write out $b_1^\dagger$ and $b_2^\dagger$ in terms of the $a^\dagger$'s, if we want to rewrite everything in terms of the original operators.
A: We are dealing here essentially with a one-particle Hamiltonian - more precisely, a Hamiltonian for many non-interacting particles. The Hamiltonian in question is simply the second quantized version of
$$
H=\omega_1 |1\rangle\langle 1|+ \omega_2 |2\rangle\langle 2| +
\frac{J}{2}\left(|2\rangle\langle 1| + |1\rangle\langle 2|\right),
$$
where the second quantized version is obtained via usual prescription
$$
\hat{\Psi}=a_1|1\rangle + a_2|2\rangle,\\
\hat{\Psi}^\dagger=a_1^\dagger\langle 1| + a_2^\dagger\langle 2|,\\
\mathcal{H}=\Psi^\dagger H\Psi
$$
Thus, we could have diagonalized the original single particle Hamiltonian (which is really a 2-by-2 matrix) and then introduce the second quantization OR we can first perform the second quantization and then diagonalize the Hamiltonian in "the space of creation and annihilation operators" as is done in the OP.
