J.R. Johansson and P.D. Nation
This ipython notebook demonstrates how to simulate the quantum vacuum rabi oscillations in the Jaynes-Cumming model, using QuTiP: The Quantum Toolbox in Python.
For more information about QuTiP see project web page: http://code.google.com/p/qutip/
# import the required python packages
%pylab inline
from qutip import *
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
The Jaynes-Cumming model is the simplest possible model of quantum mechanical light-matter interaction, describing a single two-level atom interacting with a single electromagnetic cavity mode. The Hamiltonian for this system is (in dipole interaction form)
$H = \hbar \omega_c a^\dagger a + \frac{1}{2}\hbar\omega_a\sigma_z + \hbar g(a^\dagger + a)(\sigma_- + \sigma_+)$
or with the rotating-wave approximation
$H_{\rm RWA} = \hbar \omega_c a^\dagger a + \frac{1}{2}\hbar\omega_a\sigma_z + \hbar g(a^\dagger\sigma_- + a\sigma_+)$
where $\omega_c$ and $\omega_a$ are the frequencies of the cavity and atom, respectively, and $g$ is the interaction strength.
Here we use units where $\hbar = 1$:
wc = 1.0 * 2 * pi # cavity frequency
wa = 1.0 * 2 * pi # atom frequency
g = 0.05 * 2 * pi # coupling strength
kappa = 0.005 # cavity dissipation rate
gamma = 0.05 # atom dissipation rate
N = 15 # number of cavity fock states
n_th_a = 0.0 # temperature in frequency units
use_rwa = True
tlist = linspace(0,25,100)
# intial state
psi0 = tensor(basis(N,0), basis(2,1)) # start with an excited atom
# operators
a = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))
# Hamiltonian
if use_rwa:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
else:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * (sm + sm.dag())
c_op_list = []
rate = kappa * (1 + n_th_a)
if rate > 0.0:
c_op_list.append(sqrt(rate) * a)
rate = kappa * n_th_a
if rate > 0.0:
c_op_list.append(sqrt(rate) * a.dag())
rate = gamma
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm)
Here we evolve the system with the Lindblad master equation solver, and we request that the expectation values of the operators $a^\dagger a$ and $\sigma_+\sigma_-$ are returned by the solver by passing the list [a.dag()*a, sm.dag()*sm]
as the fifth argument to the solver.
output = mesolve(H, psi0, tlist, c_op_list, [a.dag() * a, sm.dag() * sm])
Here we plot the excitation probabilities of the cavity and the atom (these expectation values were calculated by the mesolve
above). We can clearly see how energy is being coherently transferred back and forth between the cavity and the atom.
figure(figsize=(8,5))
plot(tlist, output.expect[0], label="Cavity")
plot(tlist, output.expect[1], label="Atom excited state")
legend()
xlabel('Time')
ylabel('Occupation probability')
title('Vacuum Rabi oscillations')
show()