Diagonalizing Hamiltonian in Second Quantization I have a fermionic system with states 1,2. They are coupled by a harmonic oscillator. The Hamiltonian of the system should then be
$$
H=\left[\gamma(a^\dagger+a)-\delta\right]\left( c^\dagger_1 c_2+c^\dagger_2 c_1 \right) + \xi \left ( c^\dagger_1 c_1-c^\dagger_2 c_2 \right) + \hbar\omega\left(a^\dagger a + \frac{1}{2} \right)
$$
where $\gamma$ is the coupling constant, $\delta$ some initial overlap integral and $\xi$ the energy difference between the uncoupled states.
In order to get the eigenvalues of this system I would like to diagonalize $H$. However I fail to see how to approach this because of the terms with three operators in them such as $a^\dagger c^\dagger_1 c_2$, which make me wonder how to even display $H$ as a matrix.
 A: 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.
A: The hamiltonian operator we have is that of three harmonic oscillator (two fermionic and one bosonic) in the presence of a interaction. If we assume of treating this term as a perturbation over the free theory. We can find the spectrum of the full hamiltonian by using time-independent perturbation theory. If ${\rm O}\left(H_{int}\right) $ where $\lambda$ is the perturbative parameter, the eigenvalue problem reads
\begin{gather}
\left( \hat{H}_0 + \hat{H}_{int} \right)\left\vert n\right\rangle = E_n \left\vert n \right\rangle \\
{\rm where}\\
H_0 = \xi\left( \hat{n}_{f1} - \hat{n}_{f2} \right) + \hbar\omega\left( \hat{n}_b + \frac{1}{2}\right) \\
H_{int} = \left[ \gamma\left(\hat{a}^\dagger + \hat{a} \right) - \delta \right] \Big( \hat{c}_1^\dagger \hat{c}_2 + \hat{c}_2^\dagger \hat{c}_1 \Big) 
\end{gather}
here $n_i$ is a shorthand for the oscillator's occupation number, while $n$ is the label of the eigenstate. 
If we schematically indicate the eigenstates of the free theory with $\vert n^{(0)} \rangle=\left\vert n_b,n_{f1},n_{f2} \right\rangle$, that is
\begin{align}
\hat{H}_0 \vert n^{(0)} \rangle  = E^{(0)}_n \vert n^{(0)} \rangle
\end{align}
the full theory eigenstates and eigenvalues, in time-independent perturbation theory, are given in terms of the exact free-theory solutions as eq.(1.44) of Sakurai,Napolitano - Modern quantum mechanics
\begin{align}
\vert n \rangle &= \vert n^{(0)} \rangle +  \lambda \sum_{k\neq n} \vert k^{(0)} \rangle 
\frac{\langle k^{(0)} \vert H^{(int)} \vert n^{(0)} \rangle}{E^{(0)}_n - E^{(0)}_k} + \lambda^2...
\end{align}
Before you apply to the problem you should notice that the fermionic creation operator gives zero if the corresponding state is occupied. There are few possible states


*

*$ \vert n_b,0,0 \rangle$

*$ \vert n_b,0,1 \rangle$

*$ \vert n_b,1,0 \rangle$

*$ \vert n_b,1,1 \rangle$
And $H_{int}$ conserves the sum of occupation numbers of fermionic oscillators, that is $n_{f1} + n_{f2}$ of the ket must be the same of the bra here $\langle k^{(0)} \vert H^{(int)} \vert n^{(0)} \rangle = 0$ otherwise it gives zero. The energy spectrum you're looking for is eq. (1.53b) of Sakurai, Napolitano.
Therefore, I think that the problem may further reduced to finding only the state $\vert\alpha\rangle$ since the hamiltonian eigenstates can be written as the tensor product $ \vert \alpha \rangle \times \Big( \vert 0,1 \rangle + \vert 1,0 \rangle \Big)$
