Author: J.R. Johansson, robert@riken.jp
Latest version of this ipython notebook lecture is available at: http://github.com/jrjohansson/qutip-lectures
The example in this lecture is based on an example by P.D. Nation.
# setup the matplotlib graphics library and configure it to show figures inline in the notebook
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
# make qutip available in the rest of the notebook
from qutip import *
#qutip.settings.qutip_graphics=False
#qutip.settings.qutip_gui=False
The Quantum Monte-Carlo trajectory method is an equation of motion for a single realization of the state vector $\left|\psi(t)\right>$ for a quantum system that interacts with its environment. The dynamics of the wave function is given by the Schrodinger equation,
where the Hamiltonian is an effective Hamiltonian that, in addition to the system Hamiltonian $H(t)$, alos contains a non-Hermitian contribution due to the interaction with the environment:
... incomplete ...
$\delta p = \delta t \sum_n \left<\psi(t)|c_n^\dagger c_n|\psi(t)\right>$
$\left|\psi(t+\delta t)\right> = c_n \left|\psi(t)\right>/\left<\psi(t)|c_n^\dagger c_n|\psi(t)\right>^{1/2}$
This is a Monte-Carlo simulation showing the decay of a cavity Fock state $\left|1\right>$ in a thermal environment with an average occupation number of $n=0.063$ .
Here, the coupling strength is given by the inverse of the cavity ring-down time $T_c = 0.129$ .
The parameters chosen here correspond to those from S. Gleyzes, et al., Nature 446, 297 (2007), and we will carry out a simulation that corresponds to these experimental results from that paper:
N = 4 # number of basis states to consider
kappa = 1.0/0.129 # coupling to heat bath
nth = 0.063 # temperature with <n>=0.063
tlist = linspace(0,0.6,100)
Here we create QuTiP Qobj
representations of the operators and state that are involved in this problem.
a = destroy(N) # cavity destruction operator
H = a.dag() * a # harmonic oscillator Hamiltonian
psi0 = basis(N,1) # initial Fock state with one photon: |1>
# collapse operator list
c_op_list = []
# decay operator
c_op_list.append(sqrt(kappa * (1 + nth)) * a)
# excitation operator
c_op_list.append(sqrt(kappa * nth) * a.dag())
Here we start the Monte-Carlo simulation, and we request expectation values of photon number operators with 1, 5, 15, and 904 trajectories (compare with experimental results above).
ntraj = [1, 5, 15, 904] # list of number of trajectories to avg. over
mc = mcsolve(H, psi0, tlist, c_op_list, [a.dag()*a], ntraj)
# get expectation values from mc data (need extra index since ntraj is list)
ex1 = mc.expect[0][0] # for ntraj=1
ex5 = mc.expect[1][0] # for ntraj=5
ex15 = mc.expect[2][0] # for ntraj=15
ex904 = mc.expect[3][0] # for ntraj=904
For comparison with the averages of single quantum trajectories provided by the Monte-Carlo solver we here also calculate the dynamics of the Lindblad master equation, which should agree with the Monte-Carlo simultions for infinite number of trajectories.
# run master equation to get ensemble average expectation values
me = mesolve(H, psi0, tlist, c_op_list, [a.dag()*a])
# calulate final state using steadystate solver
final_state = steadystate(H, c_op_list) # find steady-state
fexpt = expect(a.dag()*a, final_state) # find expectation value for particle number
import matplotlib.font_manager
leg_prop = matplotlib.font_manager.FontProperties(size=10)
fig, axes = subplots(4, 1, sharex=True, figsize=(8,12))
fig.subplots_adjust(hspace=0.1) # reduce space between plots
for idx, n in enumerate(ntraj):
axes[idx].plot(tlist, mc.expect[idx][0], 'b', lw=2)
axes[idx].plot(tlist, me.expect[0], 'r--', lw=1.5)
axes[idx].axhline(y=fexpt, color='k', lw=1.5)
axes[idx].set_yticks(linspace(0, 2, 5))
axes[idx].set_ylim([0, 1.5])
axes[idx].set_ylabel(r'$\left<N\right>$', fontsize=14)
if idx == 0:
axes[idx].set_title("Ensemble Averaging of Monte Carlo Trajectories")
axes[idx].legend(('Single trajectory', 'master equation', 'steady state'), prop=leg_prop)
else:
axes[idx].legend(('%d trajectories' % n, 'master equation', 'steady state'), prop=leg_prop)
axes[3].xaxis.set_major_locator(MaxNLocator(4))
axes[3].set_xlabel('Time (sec)',fontsize=14)
<matplotlib.text.Text at 0xb46df50>