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):
HI=g1(a1σ+eiΔ1t+a†1σ−iΔ1te−)+g2(a2σ+eiΔ2t+a†2σ−eiΔ2t).
I obtained this after a unitary transformation on H, U0(t)=e−iH0t, being H0=Δ2σz+ω1a†1a1+ω2a†2a2 and having introduced the detunings Δ1=Δ−ω1 and Δ2=Δ−ω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
HI(t)UI(t)=i∂∂tUI(t)
being UI(t) the time evolution operator in the interaction picture. Your solution will be obtained by computing U(t)=e−iH0tUI(t). The idea is to extract the diagonal contributions from UI(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
UI(t)=I−i∫t0dt′HI(t′)UI(t′)
where I used the fact that UI(0)=I and I have put t0=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 UI(t)=I as a first iterate
UI(t)=I−i∫t0dt′HI(t′)−∫t0dt′HI(t′)∫t′0dt″HI(t″)+…
Now, we insert our Hamiltonian in the interaction picture into this series and we get (remember that (σ+)2=(σ−)2=0)
UI(t)=I−i∫t0dt′[g1(a1σ+eiΔ1t′+a†1σ−iΔ1t′e−)+g2(a2σ+eiΔ2t′+a†2σ−eiΔ2t′)]−
−∫t0dt′∫t′0dt″[g21(a1a†1σ+σ−eiΔ1(t′−t″)+a†1a1σ−σ+e−iΔ1(t′−t″))
+g22(a2a†2σ+σ−eiΔ2(t′−t″)+a†2a2σ−σ+e−iΔ2(t′−t″))
2g1g2(a1a†2σ+σ−ei(Δ2t′−Δ1t″)+a†1a2σ−σ+e−i(Δ2t′−Δ1t″))]+…
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δH0t≈I−iδH0t+… being δH0 the correction to the unperturbed Hamiltonian. These come out from the second order terms in the Dyson series producimg
UI(t)=I−…−ig21Δ1(a1a†1σ+σ−+a†1a1σ−σ+)t
−ig22Δ2(a2a†2σ+σ−+a†2a2σ−σ+)t+…
At this point we will use the fact that [ai,a†i]=1 obtaining
UI(t)=I−…−ig21Δ1a†1a1t−ig21Δ1σ+σ−t
−ig22Δ2a†2a2t−ig22Δ2σ+σ−t+…
Turning back to U0(t)=I−iΔ2σzt−iω1a†1a1t−iω2a†2a2t+… 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
+(ω2+g22Δ2)a†2a2
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.
This post has been migrated from (A51.SE)