The Hamiltonian
$$H=\begin{pmatrix}
\alpha^2+a^\ast a&\alpha a+\beta a^\ast\\
\alpha a^\ast+\beta a&\beta^2+a^\ast a
\end{pmatrix}
$$
is known in the physics literature as the anisotropic Rabi Hamiltonian. (In the most general case there is an additional term $\Delta\sigma_z$.) I give some pointers to the literature in this Physics SE posting. The eigenvalues can be computed from a recursive scheme, but closed-form expressions for the spectrum only exist for either $\alpha=\beta$ or $\alpha\beta=0$.
$\bullet$ Consider first the case $\alpha=\beta$. A unitary transformation $H'=UHU^\ast$ with $U=e^{i\pi\sigma_y/4}$ brings the Hamiltonian to the diagonal form
$$H'=\begin{pmatrix}
b_+^\ast b_+&0\\
0&b_-^\ast b_-
\end{pmatrix},\;\; b_\pm=a\pm\alpha.
$$
The eigenvalues are the integers $N=0,1,2,\ldots$, each twofold degenerate. The corresponding eigenstates $|N,\pm\rangle$ are obtained from the eigenstates $|N\rangle$ of the harmonic oscillator by acting on these with the displacement operator,
$$| N,\pm\rangle=e^{\pm\alpha(a-a^\ast)}|N\rangle.$$
$\bullet$ At the other extreme, we can take one of the two parameters $\alpha, \beta$ much smaller than the other. Let me set $\beta=0$, $\alpha\neq 0$. (The spectrum is the same for $\alpha=0$, $\beta\neq 0$.) The Hamiltonian
$$H=\begin{pmatrix}
\alpha^2+a^\ast a&\alpha a\\
\alpha a^\ast&a^\ast a
\end{pmatrix}=(\tfrac{1}{2}\alpha^2+a^\ast a)I+\alpha(\sigma_+a+\sigma_-a^\ast)+\tfrac{1}{2}\alpha^2\sigma_z$$
is known in physics as the Jaynes-Cummings Hamiltonian.
The eigenvalues can be computed exactly, because the Hamiltonian decomposes into an infinite direct product $H_n$ of $2\times 2$-matrix Hamiltonians in the basis $|g,n+1\rangle$ (two-level system in the lower state, $n+1$ quanta excited in the oscillator) and $|e,n\rangle$ (two-level system in upper state, $n$ quanta excited):
$$H_n\begin{pmatrix}
|g,n+1\rangle\\
|e,n\rangle
\end{pmatrix}=
\begin{pmatrix}n+1&\alpha\sqrt{n+1}\\
\alpha\sqrt{n+1}&n+\alpha^2
\end{pmatrix}
\begin{pmatrix}
|g,n+1\rangle\\
|e,n\rangle
\end{pmatrix}.$$
The eigenvalues then follow directly,
$$\Omega_{n,\pm}=\tfrac{1}{2}(\alpha^2+1)+n\pm\tfrac{1}{2}\sqrt{(\alpha^2+1)^2+4\alpha^2 n}.$$
This sequence of eigenvalue pairs $\Omega_{n,+},\Omega_{n,-}$ exists for $n\in\{0,1,2,\ldots\}$. In addition, there is an unpaired ground state $|g,0\rangle$ without any excitations, which is annihilated by $H$ so it has eigenvalue $\Omega_{0}=0$ independent of $\alpha$.
$\bullet$ The general case of arbitrary $\alpha,\beta$ does not have a simple closed form expression for the eigenvalues. For $\alpha,\beta\ll 1$ a perturbative solution is given in Appendix C of arXiv:1008.1317.
The ground state does have a simple exact result: the lowest eigenvalue equals zero for any $\alpha,\beta\in\mathbb{R}$. The zero-mode – the eigenstate which is annihilated by $H$ – is a coherent state $|\xi\rangle=e^{\xi(a^\ast-a)}|0\rangle$ for the harmonic oscillator (such that $a|\xi\rangle=\xi|\xi\rangle$). Substitution of $H{p\choose q}\otimes|\xi\rangle=0$ gives the two zero-modes
$$\Psi_\pm={{\sqrt\beta}\choose{\mp\sqrt\alpha}}\otimes|\!\pm\!\!\sqrt{\alpha\beta}\rangle.$$
(Thanks to Michal Pacholski for helping me with this.)
Notation: $I$ is the $2\times 2$ unit matrix, $\sigma_x$, $\sigma_y$, $\sigma_z$ are the three Pauli matrices, and $\sigma_\pm=\tfrac{1}{2}(\sigma_x\pm i\sigma_y)$. Also please note that the eigenvalues in the OP are shifted by one unit relative to those given here, because I reordered $aa^\ast\mapsto a^\ast a$.
Since there was an issue with the appearance of spurious "parabolas" in the comparison with the numerics by yarchik (in a now deleted answer), I show here a numerical check of the spectrum as a function of $\alpha$ for $\beta=0$, using a slightly modified version of yarchik's code.
Solid curves are the numerical result, dashed and dotted curve are $\Omega_{n,\pm}$. Note that the ground state is at zero energy for any $\alpha$ when $\beta=0$. The plot below, for $\beta=2\alpha$ shows this also holds for nonzero $\beta$.
