Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson
In [1]:
%pylab inline
In [2]:
from qutip import *
In [3]:
from matplotlib import rcParams
rcParams['font.family'] = 'STIXGeneral'
rcParams['mathtext.fontset'] = 'stix'
rcParams['font.size'] = '14'
Here we follow an example from Wiseman and Milburn, Quantum measurement and control, section. 4.8.1.
Consider cavity that leaks photons with a rate $\kappa$. The dissipated photons are detected with an inefficient photon detector, with photon-detection efficiency $\eta$. The master equation describing this scenario, where a separate dissipation channel has been added for detections and missed detections, is
$\dot\rho = -i[H, \rho] + \mathcal{D}[\sqrt{1-\eta} \sqrt{\kappa} a] + \mathcal{D}[\sqrt{\eta} \sqrt{\kappa}a]$
To describe the photon measurement stochastically, we can unravelling only the dissipation term that corresponds to detections, and leaving the missed detections as a deterministic dissipation term, we obtain [Eq. (4.235) in W&M]
$d\rho = \mathcal{H}[-iH -\eta\frac{1}{2}a^\dagger a] \rho dt + \mathcal{D}[\sqrt{1-\eta} a] \rho dt + \mathcal{G}[\sqrt{\eta}a] \rho dN(t)$
or
$d\rho = -i[H, \rho] dt + \mathcal{D}[\sqrt{1-\eta} a] \rho dt -\mathcal{H}[\eta\frac{1}{2}a^\dagger a] \rho dt + \mathcal{G}[\sqrt{\eta}a] \rho dN(t)$
where
$\displaystyle \mathcal{G}[A] \rho = \frac{A\rho A^\dagger}{\mathrm{Tr}[A\rho A^\dagger]} - \rho$
$\displaystyle \mathcal{H}[A] \rho = A\rho + \rho A^\dagger - \mathrm{Tr}[A\rho + \rho A^\dagger] \rho $
and $dN(t)$ is a Poisson distributed increment with $E[dN(t)] = \eta \langle a^\dagger a\rangle (t)$.
In QuTiP we write the stochastic master equation on the form:
$\displaystyle d\rho(t) = -i[H, \rho] dt + \mathcal{D}[B] \rho dt + D_{1}[A]\rho(t) dt + D_{2}[A]\rho(t) d\xi$
where the first two term gives the deterministic master equation (Lindblad form with collapse operator $B$). Here $A = \sqrt{\eta\gamma} a$ and $B = \sqrt{(1-\eta)\gamma} $a.
We can identify
$\displaystyle D_{1}[A]\rho(t) = - \frac{1}{2}\eta\gamma \mathcal{H}[a^\dagger a] \rho(t) = -\frac{1}{2}\eta\gamma \left( a^\dagger a\rho + \rho a^\dagger a - \mathrm{Tr}[a^\dagger a\rho + \rho a^\dagger a] \rho \right) = -\frac{1}{2}\left( A^\dagger A\rho + \rho A^\dagger A - \mathrm{Tr}[A^\dagger A\rho + \rho A^\dagger A] \rho \right)$
$\displaystyle D_{2}[A]\rho(t) = \mathcal{G}[\sqrt{\eta\gamma}a] \rho = \frac{A\rho A^\dagger}{\mathrm{Tr}[A\rho A^\dagger]} - \rho$
and $d\xi = dN(t)$ with a Poisson distribution.
In [4]:
N = 15
w0 = 0.5 * 2 * pi
times = linspace(0, 15, 150)
gamma = 0.1
In [5]:
a = destroy(N)
In [6]:
H = w0 * a.dag() * a
In [7]:
rho0 = fock(N, 5)
In [8]:
e_ops = [a.dag() * a]
In [9]:
# The argument A in the d1 and d2 callback functions is a list with the following
# precomputed superoperators, where c is the stochastic collapse operator given
# to the solver (called once for each operator, if more than one is given)
#
# A[0] = spre(c)
# A[1] = spost(c)
# A[2] = spre(c.dag())
# A[3] = spost(c.dag())
# A[4] = spre(n)
# A[5] = spost(n)
# A[6] = (spre(c) * spost(c.dag())
# A[7] = lindblad_dissipator(c)
$\displaystyle D_{1}[A]\rho(t) = -\frac{1}{2}\left( A^\dagger A\rho + \rho A^\dagger A - \mathrm{Tr}[A^\dagger A\rho + \rho A^\dagger A] \rho \right)$
In [10]:
def d1_rho_func(A, rho_vec):
n_sum = A[4] + A[5]
return 0.5 * (- n_sum * rho_vec + expect_rho_vec(n_sum, rho_vec) * rho_vec)
$\displaystyle D_{2}[A]\rho(t) = \frac{A\rho A^\dagger}{\mathrm{Tr}[A\rho A^\dagger]} - \rho$
In [11]:
def d2_rho_func(A, rho_vec):
e1 = expect_rho_vec(A[6], rho_vec) + 1e-16 # add a small number to avoid division by zero
return [(A[6] * rho_vec) / e1 - rho_vec]
In [12]:
eta = 0.7
c_ops = [sqrt(1-eta) * sqrt(gamma) * a] # collapse operator B
sc_ops = [sqrt(eta) * sqrt(gamma) * a] # stochastic collapse operator A
In [13]:
result_ref = mesolve(H, rho0, times, c_ops+sc_ops, e_ops)
In [14]:
result1 = smesolve(H, rho0, times, c_ops=c_ops, sc_ops=sc_ops, e_ops=e_ops,
ntraj=1, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='poisson', store_measurement=True)
In [15]:
result2 = smesolve(H, rho0, times, c_ops=c_ops, sc_ops=sc_ops, e_ops=e_ops,
ntraj=10, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='poisson', store_measurement=True)
In [16]:
fig, axes = subplots(2,2, figsize=(12,8), sharex=True)
axes[0,0].plot(times, result1.expect[0], label=r'Stochastic ME (ntraj = 1)', lw=2)
axes[0,0].plot(times, result_ref.expect[0], label=r'Lindblad ME', lw=2)
axes[0,0].set_title("Cavity photon number (ntraj = 1)")
axes[0,0].legend()
axes[0,1].plot(times, result2.expect[0], label=r'Stochatic ME (ntraj = 10)', lw=2)
axes[0,1].plot(times, result_ref.expect[0], label=r'Lindblad ME', lw=2)
axes[0,1].set_title("Cavity photon number (ntraj = 10)")
axes[0,1].legend()
axes[1,0].step(times, np.cumsum(result1.measurement[0]), lw=2)
axes[1,0].set_title("Cummulative photon detections (ntraj = 1)")
axes[1,1].step(times, np.cumsum(array(result2.measurement).sum(axis=0)) / 10, lw=2)
axes[1,1].set_title("Cummulative avg. photon detections (ntraj = 10)")
fig.tight_layout()
In [17]:
eta = 0.1
c_ops = [sqrt(1-eta) * sqrt(gamma) * a] # collapse operator B
sc_ops = [sqrt(eta) * sqrt(gamma) * a] # stochastic collapse operator A
In [18]:
result_ref = mesolve(H, rho0, times, c_ops+sc_ops, e_ops)
In [19]:
result1 = smesolve(H, rho0, times, c_ops=c_ops, sc_ops=sc_ops, e_ops=e_ops,
ntraj=1, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='poisson', store_measurement=True)
In [20]:
result2 = smesolve(H, rho0, times, c_ops=c_ops, sc_ops=sc_ops, e_ops=e_ops,
ntraj=10, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='poisson', store_measurement=True)
In [21]:
fig, axes = subplots(2,2, figsize=(12,8), sharex=True)
axes[0,0].plot(times, result1.expect[0], label=r'Stochastic ME (ntraj = 1)', lw=2)
axes[0,0].plot(times, result_ref.expect[0], label=r'Lindblad ME', lw=2)
axes[0,0].set_title("Cavity photon number (ntraj = 1)")
axes[0,0].legend()
axes[0,1].plot(times, result2.expect[0], label=r'Stochatic ME (ntraj = 10)', lw=2)
axes[0,1].plot(times, result_ref.expect[0], label=r'Lindblad ME', lw=2)
axes[0,1].set_title("Cavity photon number (ntraj = 10)")
axes[0,1].legend()
axes[1,0].step(times, np.cumsum(result1.measurement[0]), lw=2)
axes[1,0].set_title("Cummulative photon detections (ntraj = 1)")
axes[1,1].step(times, np.cumsum(array(result2.measurement).sum(axis=0)) / 10, lw=2)
axes[1,1].set_title("Cummulative avg. photon detections (ntraj = 10)")
fig.tight_layout()
In [22]:
from qutip.ipynbtools import version_table
version_table()
Out[22]: