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.
This post has been migrated from (A51.SE)
Q&A (4871)
