<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Symbolic derivation qubit-cavity Hamiltonian
<style>
p {
    font-family
Step1:  The Jaynes-Cummings model
The Jaynes-Cummings model is one of the most elementary quantum mechanical models light-matter interaction. It describes a single two-level atom that interacts with a single harmonic-oscillator mode of a electromagnetic cavity.
The Hamiltonian for a two-level system in its eigenbasis can be written as
$$
H = \frac{1}{2}\omega_q \sigma_z
$$
and the Hamiltonian of a quantum harmonic oscillator (cavity or nanomechanical resonator, NR) is
$$
H = \hbar\omega_r (a^\dagger a + 1/2)
$$
$$
H = \hbar\omega_{NR} (b^\dagger b + 1/2)
$$
The atom interacts with the electromagnetic field produced by the cavity (NR) mode $a + a^\dagger$ ($b + b^\dagger$) through its dipole moment. The dipole-transition operators is $\sigma_x$ (which cause a transition from the two dipole states of the atom). The combined atom-cavity Hamiltonian can therefore be written in the form
$$
H = 
\hbar \omega_r (a^\dagger a + 1/2)
+
\hbar \omega_{NR} (b^\dagger b + 1/2)
+ \frac{1}{2}\hbar\Omega\sigma_z 
+
\hbar
\sigma_x \left( g(a + a^\dagger)
+
\lambda(b + b^\dagger)\right)
$$
Although the Jaynes-Cumming Hamiltonian allow us to evolve the given initial state according to the Schrödinger Equation, in an experiment we would like to predicte the response of the coupled cavity-NR-qubit system under the influence of driving fields for the cavity and qubit, and also account for the effects of dissipation and dephasing (not treated here)
The external coherent-state input may be incorporated in the Jaynes-Cummings Hamiltonian by addition of terms involving the amplitude of the driving field $\vec{E_d} \left(\vec{E_s}\right)$ and it's frequency $\omega_d\left(\omega_s\right)$
$$
H_{cavity} = 
E_d \left(e^{i\omega_dt}a +e^{-i\omega_dt}a^\dagger\right)
$$
$$
H_{qubit} = 
E_s \left(e^{i\omega_st}\sigma_- +e^{-i\omega_st}\sigma_+\right)
$$
$$
H_{NR} = 
E_p \left(e^{i\omega_st}b +e^{-i\omega_st}b^\dagger\right)
$$
To obtain the Jaynes-Cumming Hamiltonian 
$$
?
H = 
\hbar\omega_r (a^\dagger a + 1/2)
%-\frac{1}{2}\Delta\sigma_x 
+ \frac{1}{2}\hbar\Omega\sigma_z 
+
\hbar
\sigma_x \left( g(a + a^\dagger)
+
\lambda(b + b^\dagger)\right)
$$
we also need to perform a rotating-wave approximation which simplifies the interaction part of the Hamiltonian. In the following we will begin with looking at how these two Hamiltonians are related.
To represent the atom-cavity Hamiltonian in SymPy we creates an instances of the operator classes BosonOp and SigmaX, SigmaY, and SigmaZ, and use these to construct the Hamiltonian (we work in units where $\hbar = 1$).
Step2:  Unitary transformation to interaction picture
CPW (a)
Step3:  Nanoresonator - b
Step4:  Qubit - Sigma_z
Step5:  Now, in the rotating-wave approximation we can drop the fast oscillating terms containing the factors
$e^{\pm i(\omega_q + \omega_r)t} $ and $e^{\pm i(\omega_q + \omega_{NR})t}  $
Step6:  This is the interaction term of in the Jaynes-Cumming model in the interaction picture. If we transform back to the Schrödinger picture we have
Step7:  Linearized interaction
First we apply the unitary transformation $U = e^{i \omega_d a^\dagger a t}$
Step8:  We can now perform a rotating-wave approximation (RWA) by eliminating all terms that rotate with frequencies $2\omega_d$
Step9:  Introduce the detuning $\Delta_d = \omega_r - \omega_d$
Step10:  Second we apply the unitary transformation $U = e^{i \omega_p b^\dagger b t}$
Step11:  Substiture $\sigma_z = 2\sigma_+\sigma_-$
Step12:  Second we apply the unitary transformation $U = e^{i \omega_s \sigma_+ \sigma_- t}$
Step13:  Substitute $ \sigma_+\sigma_-=\dfrac{\sigma_z}{2}$
Step14:  Introduce the detuning $\Delta_{SP} = \omega_S - \omega_D$ and $\Delta_{SD} = \omega_S - \omega_D$
Step15:  This is the Jaynes-Cumming model give above, and we have now seen that it is obtained to the dipole interaction Hamiltonian through the rotating wave approximation.
Dispersive regime
In the dispersive regime, where the two-level system is detuned from the cavity by much more than the interaction strength, $\Delta \gg g, \lambda$, an effective Hamiltonian can be dervied which describes the Stark shift of the two-level system (which depends on the number of photons in the cavity) and the frequency shift of the cavity (which depend on the state of the two-level system).
This effective Hamiltonian, which is correct up to second order in the small paramter $g/\Delta$, is obtained by performing the unitary transformation
$$
U = e^{\frac{g}{\Delta}(a \sigma_- - a^\dagger \sigma_+)}
$$
$$
U = e^{\frac{\lambda}{\Delta}(b \sigma_- - b^\dagger \sigma_+)}
$$
Step16:  This is the Hamiltonian of the Jaynes-Cummings model in the the dispersive regime. It can be interpreted as the resonator having a qubit-state-dependent frequency shift, or alternatively that the qubit is feeling a resonator-photon-number dependent Stark-shift.
<ASSISTANT_TASK:>
Python Code:
from sympy import *
init_printing()
from sympsi import *
from sympsi.boson import *
from sympsi.pauli import *
Explanation: Symbolic derivation qubit-cavity Hamiltonian
<style>
p {
    font-family: "Liberation Serif", serif;
    font-size: 12pt;
}
</style>
Based on: J. R. Johansson (robert@riken.jp), http://jrjohansson.github.io, and Eunjong Kim.
Setup modules
End of explanation
# CPW, qubit and NR energies
omega_r, omega_q, omega_nr = symbols("omega_r, omega_q, omega_{NR}") 
# Coupling CPW-qubit, NR_qubit
g, L, chi, eps = symbols("g, lambda, chi, epsilon")
# Detuning 
# Drives and detunnings
Delta_d, Delta_s, Delta_p  = symbols(" Delta_d, Delta_s, Delta_p ")
A, B, C = symbols("A,B,C") # Electric field amplitude
omega_d, omega_s, omega_p = symbols("omega_d, omega_s, omega_p") # drive frequencies
# Detunning CPW-qubit, NR-qubit
Delta_CPW, Delta_NR = symbols("Delta_{CPW},Delta_{NR}")
# auxilary variables
y, x, t, Hsym = symbols("x, y, t, H ") 
Delta_SP, Delta_SD = symbols("Delta_{SP},Delta_{SD}")
# omega_r, omega_q, g, Delta_d, Delta_s, t, x, chi, Hsym = symbols("omega_r, omega_q, g, Delta_d, Delta_s, t, x, chi, H")
# A, B, C = symbols("A,B,C") # Electric field amplitude
# omega_d, omega_s = symbols("omega_d, omega_s") # 
# omega_nr, L = symbols("omega_{NR},lambda")
# Delta,Delta_t = symbols("Delta,Delta_t")
# y, omega_t = symbols("y, omega_t")
sx, sy, sz, sm, sp = SigmaX(), SigmaY(), SigmaZ(), SigmaMinus(), SigmaPlus()
a = BosonOp("a")
b = BosonOp("b")
H = omega_r * Dagger(a) * a + omega_q/2 * sz + omega_nr * Dagger(b) * b 
H_int = sx * (g  * (a + Dagger(a)) + L * (b + Dagger(b)))
H_drive_r = A * (exp(I*omega_d*t)*a + exp(-I*omega_d*t)*Dagger(a))
H_drive_q = B * (exp(I*omega_s*t)*sm + exp(-I*omega_s*t)*sp)
H_drive_NR = C * (exp(I*omega_p*t)*b + exp(-I*omega_p*t)*Dagger(b))
H_total =  H+  H_int + H_drive_r + H_drive_q + H_drive_NR
Eq(Hsym,H_total)
Explanation: The Jaynes-Cummings model
The Jaynes-Cummings model is one of the most elementary quantum mechanical models light-matter interaction. It describes a single two-level atom that interacts with a single harmonic-oscillator mode of a electromagnetic cavity.
The Hamiltonian for a two-level system in its eigenbasis can be written as
$$
H = \frac{1}{2}\omega_q \sigma_z
$$
and the Hamiltonian of a quantum harmonic oscillator (cavity or nanomechanical resonator, NR) is
$$
H = \hbar\omega_r (a^\dagger a + 1/2)
$$
$$
H = \hbar\omega_{NR} (b^\dagger b + 1/2)
$$
The atom interacts with the electromagnetic field produced by the cavity (NR) mode $a + a^\dagger$ ($b + b^\dagger$) through its dipole moment. The dipole-transition operators is $\sigma_x$ (which cause a transition from the two dipole states of the atom). The combined atom-cavity Hamiltonian can therefore be written in the form
$$
H = 
\hbar \omega_r (a^\dagger a + 1/2)
+
\hbar \omega_{NR} (b^\dagger b + 1/2)
+ \frac{1}{2}\hbar\Omega\sigma_z 
+
\hbar
\sigma_x \left( g(a + a^\dagger)
+
\lambda(b + b^\dagger)\right)
$$
Although the Jaynes-Cumming Hamiltonian allow us to evolve the given initial state according to the Schrödinger Equation, in an experiment we would like to predicte the response of the coupled cavity-NR-qubit system under the influence of driving fields for the cavity and qubit, and also account for the effects of dissipation and dephasing (not treated here)
The external coherent-state input may be incorporated in the Jaynes-Cummings Hamiltonian by addition of terms involving the amplitude of the driving field $\vec{E_d} \left(\vec{E_s}\right)$ and it's frequency $\omega_d\left(\omega_s\right)$
$$
H_{cavity} = 
E_d \left(e^{i\omega_dt}a +e^{-i\omega_dt}a^\dagger\right)
$$
$$
H_{qubit} = 
E_s \left(e^{i\omega_st}\sigma_- +e^{-i\omega_st}\sigma_+\right)
$$
$$
H_{NR} = 
E_p \left(e^{i\omega_st}b +e^{-i\omega_st}b^\dagger\right)
$$
To obtain the Jaynes-Cumming Hamiltonian 
$$
?
H = 
\hbar\omega_r (a^\dagger a + 1/2)
%-\frac{1}{2}\Delta\sigma_x 
+ \frac{1}{2}\hbar\Omega\sigma_z 
+
\hbar
\sigma_x \left( g(a + a^\dagger)
+
\lambda(b + b^\dagger)\right)
$$
we also need to perform a rotating-wave approximation which simplifies the interaction part of the Hamiltonian. In the following we will begin with looking at how these two Hamiltonians are related.
To represent the atom-cavity Hamiltonian in SymPy we creates an instances of the operator classes BosonOp and SigmaX, SigmaY, and SigmaZ, and use these to construct the Hamiltonian (we work in units where $\hbar = 1$).
End of explanation
U = exp(I  * omega_r * t * Dagger(a)*a)
# U
H1 = hamiltonian_transformation(U, H_total.expand(), independent=True)
# H1
Explanation: Unitary transformation to interaction picture
CPW (a)
End of explanation
U = exp(I  * omega_nr * t * Dagger(b)*b)
# U
H1a = hamiltonian_transformation(U, H1.expand(), independent=True)
# H1a
Explanation: Nanoresonator - b
End of explanation
U = exp(I * omega_q * t * sp * sm)
# U
H2 = hamiltonian_transformation(U, H1a.expand())
H2 = H2.subs(sx, sm + sp).expand()
H2 = powsimp(H2)
# H2
# H2
# trick to simplify exponents
def simplify_exp(e):
    if isinstance(e, exp):
        return exp(simplify(e.exp.expand()))
    if isinstance(e, (Add, Mul)):
        return type(e)(*(simplify_exp(arg) for arg in e.args)) 
    return e
H3a = simplify_exp(H2).subs(-omega_r + omega_q, Delta_CPW)
H3 = simplify_exp(H3a).subs(-omega_nr + omega_q, Delta_NR)
# H3
Explanation: Qubit - Sigma_z
End of explanation
H4 = drop_terms_containing(H3, [exp( I * (omega_q + omega_r) * t),
                                exp(-I * (omega_q + omega_r) * t),
                                exp(I * (omega_q + omega_nr) * t),
                                exp(-I * (omega_q + omega_nr) * t)],
                          )
H4 = drop_c_number_terms(H4.expand())
Eq(Hsym, H4)
Explanation: Now, in the rotating-wave approximation we can drop the fast oscillating terms containing the factors
$e^{\pm i(\omega_q + \omega_r)t} $ and $e^{\pm i(\omega_q + \omega_{NR})t}  $
End of explanation
U = exp(-I * omega_r * t * Dagger(a) * a)
H5 = hamiltonian_transformation(U, H4.expand(), independent=True)
# H5
U = exp(-I * omega_nr * t * Dagger(b) * b)
H5a = hamiltonian_transformation(U, H5.expand(), independent=True)
# H5a
U = exp(-I * omega_q * t * sp * sm)
H6 = hamiltonian_transformation(U, H5a.expand())
# H6
H7a = simplify_exp(H6).subs(Delta_CPW, omega_q - omega_r)
H7 = simplify_exp(H7a).subs(Delta_NR, omega_q - omega_nr)
H7 = simplify_exp(powsimp(H7)).expand()
H7 = drop_c_number_terms(H7)
H = collect(H7, [A,B, C,g,L])
Eq(Hsym, H)
Explanation: This is the interaction term of in the Jaynes-Cumming model in the interaction picture. If we transform back to the Schrödinger picture we have:
End of explanation
U = exp(I * Dagger(a) * a * omega_d * t)
# U
H1 = hamiltonian_transformation(U, H, independent=True)
# H1
Explanation: Linearized interaction
First we apply the unitary transformation $U = e^{i \omega_d a^\dagger a t}$:
End of explanation
H2 = drop_terms_containing(H1.expand(), [exp(-2*I*omega_d*t), exp(2*I*omega_d*t)])
H2 = H2.collect([A,B,Dagger(a)*a,g])
# H2
# Eq(Symbol("H_{rwa}"), H2)
Explanation: We can now perform a rotating-wave approximation (RWA) by eliminating all terms that rotate with frequencies $2\omega_d$:
End of explanation
H3 = H2.subs(-omega_d+omega_r, Delta_d).expand()
H3 = H3.collect([A,B,C,Dagger(a)*a,Dagger(b)*b,g,L])
# H3
H3
Explanation: Introduce the detuning $\Delta_d = \omega_r - \omega_d$:
End of explanation
U = exp(I * Dagger(b) * b * omega_p * t)
H3a = hamiltonian_transformation(U, H3, independent=True)
H3b = drop_terms_containing(H3a.expand(), [exp(-2*I*omega_p*t), exp(2*I*omega_p*t)])
H3b = H3b.collect([A,B,C,L,Dagger(a)*a,Dagger(b)*b,g])
H3c = H3b.subs(-omega_p+omega_nr, Delta_p).expand()
H3c = H3c.collect([A,B,C,L,Dagger(a)*a,Dagger(b)*b,g])
H3c
Explanation: Second we apply the unitary transformation $U = e^{i \omega_p b^\dagger b t}$:
End of explanation
H3 = H3c.subs(sz,2*sp*sm).expand()
# H3
H3 = H3.collect([A,B,C,L,Dagger(a)*a,Dagger(b)*b,g])
# H3
# H3
Explanation: Substiture $\sigma_z = 2\sigma_+\sigma_-$
End of explanation
U = exp(I * sp*sm* omega_s * t)
# U
H4 = hamiltonian_transformation(U, H3, independent=True)
H4 =H4.expand()
H4
Explanation: Second we apply the unitary transformation $U = e^{i \omega_s \sigma_+ \sigma_- t}$:
End of explanation
H4 = H4.subs(omega_q*sp*sm,omega_q*sz/2)
# H4
H4 = H4.subs(-omega_s*sp*sm,-omega_s*sz/2)
# H4 = H4.subs(-omega_s*sp*sm,-omega_s*sz/2)
H4= H4.collect([A,B,C,L,Dagger(a)*a,Dagger(b)*b,g,sz])
H4
H4 = H4.subs(omega_q/2 - omega_s/2,Delta_s/2).expand()
H4= H4.collect([A,B,C,L,Dagger(a)*a,Dagger(b)*b,g])
H4
Explanation: Substitute $ \sigma_+\sigma_-=\dfrac{\sigma_z}{2}$
End of explanation
H4 = simplify_exp(powsimp(H4))
H4 = H4.subs(omega_p - omega_s,-Delta_SP).expand()
H4 = simplify_exp(powsimp(H4))
H4 = H4.subs(omega_d - omega_s,-Delta_SD).expand()
H5 = H4.collect([A,B,C,L,Dagger(a)*a,Dagger(b)*b,g])
# H5
H5 = H5.expand()
# H5
H6 = simplify_exp(powsimp(H5))
# H6 = simplify_exp(H6).subs(omega_t - omega_s, Delta_t)
H = H6.collect([A,B,C,Dagger(a)*a,g,sz,L])
Eq(Hsym, H)
Explanation: Introduce the detuning $\Delta_{SP} = \omega_S - \omega_D$ and $\Delta_{SD} = \omega_S - \omega_D$:
End of explanation
U = exp(x * (a * sp).expand())
# U
#H1 = unitary_transformation(U, H, allinone=True, expansion_search=False, N=3).expand()
#H1 = qsimplify(H1)
#H1
H1 = hamiltonian_transformation(U, H, expansion_search=False, N=3).expand()
H1a = qsimplify(H1)
# H1a
U = exp(y * (b * sp).expand())
# U
H1a = hamiltonian_transformation(U, H1a, expansion_search=False, N=3).expand()
H1 = qsimplify(H1a)
# H1
U = exp(-x * (Dagger(a) * sm).expand())
# U
H2a = hamiltonian_transformation(U, H1, expansion_search=False, N=3).expand()
H2a = qsimplify(H2a)
# H2a
U = exp(-y * (Dagger(b) * sm).expand())
# U
H2 = hamiltonian_transformation(U, H2a, expansion_search=False, N=3).expand()
H2 = qsimplify(H2)
# H2
# H3 = drop_terms_containing(H2.expand(), [x**2,x**3, x**4,x**5,x**6,x**7, y**2,y**3, y**4,y**5,y**6,y**7])
H3 = drop_terms_containing(H2.expand(), [x**3, x**4,x**5,x**6,x**7, y**3, y**4,y**5,y**6,y**7])
# H3
H4 = H3.subs(x, g/Delta_CPW)
H4a = H4.subs(y, L/Delta_NR)
# H4
H4b = qsimplify(H4a)
H5 = drop_c_number_terms(H4b)
# H6a = drop_terms_containing(H5.expand(), [exp(I*Delta_SP*t), exp(-I*Delta_SP*t)])
# H6 = drop_terms_containing(H6a.expand(), [exp(I*Delta_SD*t), exp(-I*Delta_SD*t)])
# H6
H6 = collect(H5, [A,B,C,Dagger(a)*a,Dagger(b)*b,sz,g,L])
# H6 
U = exp(I * omega_r * t * Dagger(a) * a)
H7 = hamiltonian_transformation(U, H6.expand()).expand()
H7a = qsimplify(H7)
U = exp(I * omega_nr * t * Dagger(b) * b)
H7b = hamiltonian_transformation(U, H7a).expand()
H7b = qsimplify(H7b)
U = exp(I * omega_q * t * sp * sm)
H7c = hamiltonian_transformation(U, H7b, expansion_search=False, N=3).expand()
H7c = qsimplify(H7c)
H8 = drop_terms_containing(H7, [exp(I * omega_r * t), exp(-I * omega_r * t),
                                exp(I * omega_nr * t), exp(-I * omega_nr * t),
                                exp(I * omega_q * t), exp(-I * omega_q * t)])
U = exp(- I * omega_r * t * Dagger(a) * a)
H9 = hamiltonian_transformation(U, H8.expand(), expansion_search=False, N=3).expand()
H9 = qsimplify(H9)
U = exp(-I * omega_nr * t * Dagger(b) * b)
H10 = hamiltonian_transformation(U, H9, expansion_search=False, N=3).expand()
H10 = qsimplify(H10)
U = exp(-I * omega_q * t * sp * sm)
H11 = hamiltonian_transformation(U, H10, expansion_search=False, N=3).expand()
H11 = qsimplify(H11)
# H8 = qsimplify(H7.expand())
# H9 = collect(H8, [A,B,C,Dagger(a)*a,g,sz,L])
# # H9
# H10 = qsimplify(H9)
H12 = collect(H11, [A,B,C,Dagger(a)*a,g,sz,L])
# H10
Explanation: This is the Jaynes-Cumming model give above, and we have now seen that it is obtained to the dipole interaction Hamiltonian through the rotating wave approximation.
Dispersive regime
In the dispersive regime, where the two-level system is detuned from the cavity by much more than the interaction strength, $\Delta \gg g, \lambda$, an effective Hamiltonian can be dervied which describes the Stark shift of the two-level system (which depends on the number of photons in the cavity) and the frequency shift of the cavity (which depend on the state of the two-level system).
This effective Hamiltonian, which is correct up to second order in the small paramter $g/\Delta$, is obtained by performing the unitary transformation
$$
U = e^{\frac{g}{\Delta}(a \sigma_- - a^\dagger \sigma_+)}
$$
$$
U = e^{\frac{\lambda}{\Delta}(b \sigma_- - b^\dagger \sigma_+)}
$$
End of explanation
H12
Explanation: This is the Hamiltonian of the Jaynes-Cummings model in the the dispersive regime. It can be interpreted as the resonator having a qubit-state-dependent frequency shift, or alternatively that the qubit is feeling a resonator-photon-number dependent Stark-shift.
End of explanation
<END_TASK>