Actually, after reading my reply I understood that all of this is by no means necessary. The solution is extremely simple.
\begin{align}
\left( \hat{H}_0 + \hat{H}_{int} \right) \vert \Psi_k \rangle = E_k \vert \Psi_k \rangle
\end{align}
The most general state $\vert \Psi_k \rangle$ is in the tensor space therefore it is of the form
\begin{align}
\vert \Psi_k \rangle = \sum_{n_1,f_1,f_2} c_{1,2,3} \;\vert n_1 \rangle \times \vert f_1\rangle \times \vert f_2 \rangle
\end{align}
The fermion product $\vert f_1 \rangle \times \vert f_2 \rangle $ is one of the 4 possible candidates
1) $\vert0,0\rangle$
2) $\vert0,1\rangle$
3) $\vert1,0\rangle$
4) $\vert1,1\rangle$
The free hamiltonian has eigenvectors $\vert n_b, f_1, f_2 \rangle$ with eigenvalues $E(n_b,f_1,f_2) = \hbar\left(n_b + \frac{1}{2}\right) + \xi(f_1-f_2)$. Here the term proportional to $\xi$ contributes only if $\vert f_1 \neq f_2 \vert$. However, the states 1-4 are not eigenstates of $H_{int}$. A linear combination of these eigenstates is still an eigenstate of $H_0$. Therefore I seek for eigenstates of $H_{int}$ in the form of such linear combination.
Let us analyse the interaction term $\hat{H}_{int}$. It is the sum of two term, the first in which there are three operators ($a c^\dagger c$ etc.) and the last with two operator.
\begin{align}
\hat{H}_{int} &= \hat{H}_{int,1} + \hat{H}_{int,2}\\
&=\gamma \Big( \hat{a}^\dagger + \hat{a} \Big)\left(\hat{c}_1^\dagger \hat{c}_2 + \hat{c}_2^\dagger \hat{c}_1\right) - \delta \left(\hat{c}_1^\dagger \hat{c}_2 + \hat{c}_2^\dagger \hat{c}_1\right)
\end{align} Let us start from the latter - which is easier. The eigenvectors of $H_{int2}$ are the rays $c_1\Big( \vert \alpha ,1,0\rangle + \vert \alpha ,0,1 \rangle\Big)$ for any $\vert\alpha \rangle \equiv \sum_{n_b=0}^{\infty} c_{n_b} \vert n_b\rangle$. Also, we have the linear combinations $c_1 \vert \alpha ,0,0\rangle + c_2 \vert \alpha ,1,1\rangle $.
Finally, let us turn to $H_{int1}$. The fermionic part is analogue, however this time we have the bosonic creation/annihilation operators. Therefore any linear combination
$c_1 \vert \alpha ,0,0\rangle + c_2 \vert \alpha ,1,1\rangle $ is eigenstate (with zero eigenvalue). The other case is where I use the state $c_1\Big( \vert \alpha ,1,0\rangle + \vert \alpha ,0,1 \rangle\Big)$. Here $\vert \alpha\rangle $ is a generic linear combination of number-eigenstates of the bosonic operator $n_b$. The result is
\begin{align}
\hat{H}_{int1} \vert \alpha\rangle \times \vert(...)\rangle &=\gamma \Big(\hat{a}^\dagger + \hat{a} \Big) \sum_{n_b=0}^{\infty} c_{n_b} \vert n_b\rangle \times \vert(...)\rangle \\
&=\gamma \Big(\sum_{n_b=1}^{\infty} \sqrt{n_b}\,c_{n_b-1} \vert n_b\rangle + \sum_{n_b=0}^{\infty} \sqrt{n_b+1}\,c_{n_b+1} \vert n_b\rangle \Big)\times \vert(...) \\
&=\gamma \Big(\sum_{n_b=1}^{\infty} (\sqrt{n_b}\,c_{n_b-1}+ \sqrt{n_b+1}\,c_{n_b+1}) \vert n_b\rangle\Big)\times \vert(...) + c_1 \vert0\rangle\times \vert(...)
\end{align}
Comparing the lhs with rhs I've $c_1=c_0$ and $c_{n_b} = (\sqrt{n_b}\,c_{n_b-1}+ \sqrt{n_b+1}\,c_{n_b+1})$ therefore $c_{n_b + 1} = \frac{c_{n_b} - \sqrt{n_b}\,c_{n_b-1}}{\sqrt{n_b+1}}$
Modulo few errors I've made in ladder operators factors etc., you can find the state $\vert \alpha \rangle$ and the eigenvectors of the full $\hat{H}$ taking the intersection of the three sets.