Author: J. R. Johansson (robert@riken.jp), http://dml.riken.jp/~rob/
The latest version of this IPython notebook lecture is available at http://github.com/jrjohansson/qutip-lectures.
The other notebooks in this lecture series are indexed at http://jrjohansson.github.com.
In [1]:
# setup the matplotlib graphics library and configure it to show
# figures inline in the notebook
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
# make qutip available in the rest of the notebook
from qutip import *
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)
or with the rotating-wave approximation
where $\omega_c$ and $\omega_a$ are the frequencies of the cavity and atom, respectively, and $g$ is the interaction strength.
In [4]:
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 # avg number of thermal bath excitation
use_rwa = True
tlist = np.linspace(0,25,101)
In [5]:
# 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())
In [6]:
c_ops = []
# cavity relaxation
rate = kappa * (1 + n_th_a)
if rate > 0.0:
c_ops.append(sqrt(rate) * a)
# cavity excitation, if temperature > 0
rate = kappa * n_th_a
if rate > 0.0:
c_ops.append(sqrt(rate) * a.dag())
# qubit relaxation
rate = gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm)
In [7]:
output = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])
In [8]:
n_c = output.expect[0]
n_a = output.expect[1]
fig, axes = plt.subplots(1, 1, figsize=(10,6))
axes.plot(tlist, n_c, label="Cavity")
axes.plot(tlist, n_a, label="Atom excited state")
axes.legend(loc=0)
axes.set_xlabel('Time')
axes.set_ylabel('Occupation probability')
axes.set_title('Vacuum Rabi oscillations')
Out[8]:
In addition to the cavity's and atom's excitation probabilities, we may also be interested in for example the wigner function as a function of time. The Wigner function can give some valuable insight in the nature of the state of the resonators.
To calculate the Wigner function in QuTiP, we first recalculte the evolution without specifying any expectation value operators, which will result in that the solver return a list of density matrices for the system for the given time coordinates.
In [9]:
output = mesolve(H, psi0, tlist, c_ops, [])
Now, output.states contains a list of density matrices for the system for the time points specified in the list tlist:
In [10]:
output
Out[10]:
In [11]:
type(output.states)
Out[11]:
In [12]:
len(output.states)
Out[12]:
In [13]:
output.states[-1] # indexing the list with -1 results in the last element in the list
Out[13]:
Now let's look at the Wigner functions at the point in time when atom is in its ground state: $t = \\{5, 15, 25\\}$ (see the plot above).
For each of these points in time we need to:
In [14]:
# find the indices of the density matrices for the times we are interested in
t_idx = where([tlist == t for t in [0.0, 5.0, 15.0, 25.0]])[1]
tlist[t_idx]
Out[14]:
In [15]:
# get a list density matrices
rho_list = array(output.states)[t_idx]
In [17]:
# loop over the list of density matrices
xvec = np.linspace(-3,3,200)
fig, axes = plt.subplots(1,len(rho_list), sharex=True, figsize=(3*len(rho_list),3))
for idx, rho in enumerate(rho_list):
# trace out the atom from the density matrix, to obtain
# the reduced density matrix for the cavity
rho_cavity = ptrace(rho, 0)
# calculate its wigner function
W = wigner(rho_cavity, xvec, xvec)
# plot its wigner function
axes[idx].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-.25,.25), cmap=plt.get_cmap('RdBu'))
axes[idx].set_title(r"$t = %.1f$" % tlist[t_idx][idx], fontsize=16)
At $t =0$, the cavity is in it's ground state. At $t = 5, 15, 25$ it reaches it's maxium occupation in this Rabi-vacuum oscillation process. We can note that for $t=5$ and $t=15$ the Wigner function has negative values, indicating a truely quantum mechanical state. At $t=25$, however, the wigner function no longer has negative values and can therefore be considered a classical state.
In [22]:
t_idx = where([tlist == t for t in [0.0, 5.0, 10, 15, 20, 25]])[1]
rho_list = array(output.states)[t_idx]
fig_grid = (2, len(rho_list)*2)
fig = plt.figure(figsize=(2.5*len(rho_list),5))
for idx, rho in enumerate(rho_list):
rho_cavity = ptrace(rho, 0)
W = wigner(rho_cavity, xvec, xvec)
ax = plt.subplot2grid(fig_grid, (0, 2*idx), colspan=2)
ax.contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-.25,.25), cmap=plt.get_cmap('RdBu'))
ax.set_title(r"$t = %.1f$" % tlist[t_idx][idx], fontsize=16)
# plot the cavity occupation probability in the ground state
ax = plt.subplot2grid(fig_grid, (1, 1), colspan=(fig_grid[1]-2))
ax.plot(tlist, n_c, label="Cavity")
ax.plot(tlist, n_a, label="Atom excited state")
ax.legend()
ax.set_xlabel('Time')
ax.set_ylabel('Occupation probability');
In [23]:
from qutip.ipynbtools import version_table
version_table()
Out[23]: