QuTiP lecture: Quantum Monte-Carlo Trajectories

Author: J.R. Johansson, [email protected]

http://dml.riken.jp/~rob/

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.

In [1]:
# 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)'.
In [2]:
# make qutip available in the rest of the notebook
from qutip import *

#qutip.settings.qutip_graphics=False
#qutip.settings.qutip_gui=False

Introduction to the Quantum Monte-Carlo trajectory method

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,

$\frac{d}{dt}\left|\psi(t)\right> = - \frac{i}{\hbar} H_{\rm eff} \left|\psi(t)\right>$

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:

$H_{\rm eff}(t) = H(t) - \frac{i\hbar}{2}\sum_n c_n^\dagger c_n$

... 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}$

Decay of a single-photon Fock state in a cavity

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:

Problem parameters

In [3]:
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)

Create operators, Hamiltonian and initial state

Here we create QuTiP Qobj representations of the operators and state that are involved in this problem.

In [4]:
a = destroy(N)      # cavity destruction operator
H = a.dag() * a     # harmonic oscillator Hamiltonian
psi0 = basis(N,1)   # initial Fock state with one photon: |1>

Create a list of collapse operators that describe the dissipation

In [5]:
# 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())

Run Monte-Carlo simulation

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).

In [31]:
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

Lindblad master-equation simulation and steady state

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.

In [32]:
# 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

Plot the results

In [36]:
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)
Out[36]:
<matplotlib.text.Text at 0xb46df50>