What are some ways to (approximately) symbolically diagonalize Hamiltonian operator? Specifically the Hamiltonian takes the form of
$$\hat H = \frac{\Delta }{2}{\hat \sigma _z} + {\omega _1}\hat a_1^\dagger {\hat a_1} + {\omega _2}\hat a_2^\dagger {\hat a_2} + {g_1}\left( {{{\hat a}_1}{{\hat \sigma }_ + } + \hat a_1^\dagger {{\hat \sigma }_ - }} \right) + {g_2}\left( {{{\hat a}_2}{{\hat \sigma }_ + } + \hat a_2^\dagger {{\hat \sigma }_ - }} \right),$$
a three body version of Jaynes-Cummings model.
I'm currently trying to diagonalize this Hamiltonian, a first step in our application of quantum Zeno effect to a three-body system. 
I guess this Hamiltonian simply has no close-form diagonalization, just like in classical physics there is no closed-form general solution for a three-body system. So my question is: what are several symbolic approximation techniques to diagonalize an Hermitian operator? Better if that techniques particularly suits this Hamiltonian. The values of $\Delta, \omega_1, \omega_2, g_1, g_2$ need not be general; they can be set, say, all equal in order to simplify calculation.
 A: At least for $\omega_1=\omega_2$ it is possible to solve the system exactly. Our Hamiltonian can be written
$$\hat H = \frac{\Delta }{2}{\hat \sigma _z} + {\omega}(\hat a_1^\dagger {\hat a_1} + \hat a_2^\dagger {\hat a_2}) + ({g_1}{\hat a}_1+{g_2}{\hat a}_2){{\hat \sigma }_ + } + ({g_1}{\hat a}_1^\dagger+{g_2}{\hat a}_2^\dagger){{\hat \sigma }_ - }$$
We apply a change of variables
$$a_1'=a_1 cos\theta-a_2sin\theta$$
$$a_2'=a_1 sin\theta+a_2cos\theta$$
This change of variables preserves the $\Delta$ and the $\omega$ terms in the Hamiltonian but rotates the $(g_1, g_2)$ vector. By choosing an appropriate $\theta$ we can achieve $g_2'=0$. This means the $(a_1', \sigma)$ system decouples from $a_2'$. The former is an ordinary Jaynes-Cummings model whereas the later is a Harmonic oscillator
A: This is a Jaynes-Cummings model with a two mode electromagnetic field. One can find some literature about (e.g. see here and refs therein) and some exact solutions are also known. But, if you content yourself with a perturbation soulution for the full parameter Hamiltonian, this can be accomplished in the following way.
Firstly, move to interaction picture. This will give the following Hamiltonian (I will remove the hats as there is no difficulty to tell operators):
$$H_I = {g_1}\left( {{{a}_1}{{\sigma }_ + }e^{i\Delta_1t} + a_1^\dagger {{\sigma }_ e^{-i\Delta_1t}- }} \right) + {g_2}\left( {{{a}_2}{{\sigma }_ + }e^{i\Delta_2t} + a_2^\dagger {{\sigma }_ - }e^{i\Delta_2t}} \right).$$
I obtained this after a unitary transformation on $H$, $U_0(t)=e^{-iH_0t}$, being $H_0=\frac{\Delta}{2}\sigma_z+\omega_1a_1^\dagger a_1+\omega_2a_2^\dagger a_2$ and having introduced the detunings $\Delta_1=\Delta-\omega_1$ and $\Delta_2=\Delta-\omega_2$. Interesting physics comes into play when these detunings are allowed to go to 0 but here we take them non null. Now, your problem is to solve the equation
$$H_I(t)U_I(t)=i\frac{\partial}{\partial t}U_I(t)$$
being $U_I(t)$ the time evolution operator in the interaction picture. Your solution will be obtained by computing $U(t)=e^{-iH_0t}U_I(t)$. The idea is to extract the diagonal contributions from $U_I(t)$ and evaluate from these the approximate eigenvalues for the Hamiltonian we started from. Now, the Schroedinger equation can be rewritte in integral form as
$$U_I(t)=I-i\int_0^tdt'H_I(t')U_I(t')$$
where I used the fact that $U_I(0)=I$ and I have put $t_0=0$ being this arbitrary. This is now an integral equation as you realize that the unknown is also under the integral sign. But this eqaution can be solved iteratively to give the so called Dyson series for time dependent perturbations starting with a solution $U_I(t)=I$ as a first iterate
$$U_I(t)=I-i\int_0^tdt'H_I(t')-\int_0^tdt'H_I(t')\int_0^{t'}dt''H_I(t'')+\ldots$$
Now, we insert our Hamiltonian in the interaction picture into this series and we get (remember that $(\sigma_+)^2=(\sigma_-)^2=0$)
$$U_I(t)=I-i\int_0^tdt'\left[{g_1}\left( {{{a}_1}{{\sigma }_ + }e^{i\Delta_1t'} + a_1^\dagger {{\sigma }_ e^{-i\Delta_1t'}- }} \right) + {g_2}\left( {{{a}_2}{{\sigma }_ + }e^{i\Delta_2t'} + a_2^\dagger {{\sigma }_ - }e^{i\Delta_2t'}} \right)\right]-$$
$$-\int_0^tdt'\int_0^{t'}dt''\left[g_1^2\left(a_1a_1^\dagger\sigma_+\sigma_-e^{i\Delta_1(t'-t'')}+a_1^\dagger a_1\sigma_-\sigma_+e^{-i\Delta_1(t'-t'')}\right)\right.$$
$$+g_2^2\left(a_2a_2^\dagger\sigma_+\sigma_-e^{i\Delta_2(t'-t'')}+a_2^\dagger a_2\sigma_-\sigma_+e^{-i\Delta_2(t'-t'')}\right)$$
$$\left.2g_1g_2\left(a_1a_2^\dagger\sigma_+\sigma_-e^{i(\Delta_2 t'-\Delta_1t'')}+a_1^\dagger a_2\sigma_-\sigma_+e^{-i(\Delta_2 t'-\Delta_1 t'')}\right)\right]+\ldots$$
The terms we are interested in are those having a "secular" behavior as they increses lienarly with time and these are just the first temrs of the expansion of a time evolution operator in the form $e^{-i\delta H_0 t}\approx I-i\delta H_0 t+\ldots$ being $\delta H_0$ the correction to the unperturbed Hamiltonian. These come out from the second order terms in the Dyson series producimg
$$U_I(t)=I-\ldots-i\frac{g_1^2}{\Delta_1}\left(a_1a_1^\dagger\sigma_+\sigma_-+a_1^\dagger a_1\sigma_-\sigma_+\right)t$$
$$-i\frac{g_2^2}{\Delta_2}\left(a_2a_2^\dagger\sigma_+\sigma_-+a_2^\dagger a_2\sigma_-\sigma_+\right)t+\ldots$$
At this point we will use the fact that $[a_i,a_i^\dagger]=1$ obtaining
$$U_I(t)=I-\ldots-i\frac{g_1^2}{\Delta_1}a_1^\dagger a_1t-i\frac{g_1^2}{\Delta_1}\sigma_+\sigma_-t$$
$$-i\frac{g_2^2}{\Delta_2}a_2^\dagger a_2t-i\frac{g_2^2}{\Delta_2}\sigma_+\sigma_-t+\ldots$$
Turning back to $U_0(t)=I-i\frac{\Delta}{2}\sigma_zt-i\omega_1a_1^\dagger a_1t-i\omega_2a_2^\dagger a_2t+\ldots$ and collecting all together we recognize a new unperturbed Hamiltonian given by
$$H_0^{'}=\frac{\Delta}{2}\sigma_z+\frac{g_1^2}{\Delta_1}\frac{1}{2}(1+\sigma_z)+\frac{g_2^2}{\Delta_2}\frac{1}{2}(1+\sigma_z)+ \left( \omega_1+\frac{g_1^2}{\Delta_1} \right) a_1^\dagger a_1$$
$$+\left(\omega_2+\frac{g_2^2}{\Delta_2}\right)a_2^\dagger a_2$$
from which you can read off the new eigenvalues. You will notice that the excited state goes higher. This is exactly what one should expect from Rayleigh-Schroedinger stationary perturbation method. From the other terms of the series you will able to recover also the eigenvectors (see here). You can also extend this approach to the case when one of the two detunings is zero as one of the two contribution is just the well-knwon Jaynes-Cummings Hamiltonian for a single mode.
