Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson
In this notebook we test the qutip stochastic master equation solver (smesolve) with a few textbook examples taken from the book Quantum Optics, by Walls and Milburn, section 6.7.
In [1]:
%pylab inline
In [2]:
from qutip import *
Stochastic master equation in Milburn's formulation
$\displaystyle d\rho(t) = dN(t) \mathcal{G}[a] \rho(t) - dt \gamma \mathcal{H}[a^\dagger a] \rho(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)] = \gamma \langle a^\dagger a\rangle (t)$.
In QuTiP we write the stochastic master equation on the form (in the interaction picture, with no deterministic dissipation):
$\displaystyle d\rho(t) = D_{1}[A]\rho(t) dt + D_{2}[A]\rho(t) dW$
where $A = \sqrt{\gamma} a$, so we can identify
$\displaystyle D_{1}[A]\rho(t) = - \gamma \mathcal{H}[a^\dagger a] \rho(t) = -\gamma \left( a^\dagger a\rho + \rho a^\dagger a - \mathrm{Tr}[a^\dagger a\rho + \rho a^\dagger a] \rho \right) = -\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}[a] \rho = \frac{A\rho A^\dagger}{\mathrm{Tr}[A\rho A^\dagger]} - \rho$
and
$dW = dN(t)$
and $A = \sqrt{\gamma} a$ is the collapse operator including the rate of the process as a coefficient in the operator.
In [3]:
N = 15
w0 = 0.5 * 2 * pi
times = linspace(0, 15, 150)
gamma = 0.15
In [4]:
a = destroy(N)
In [5]:
H = w0 * a.dag() * a
In [6]:
rho0 = fock(N, 5)
In [7]:
c_ops = [sqrt(gamma) * a]
In [8]:
e_ops = [a.dag() * a, a + a.dag()]
In [9]:
result_ref = mesolve(H, rho0, times, c_ops, e_ops)
In [10]:
plot_expectation_values(result_ref);
In [11]:
# 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] = -\gamma \left( a^\dagger a\rho + \rho a^\dagger a - \mathrm{Tr}[a^\dagger a\rho + \rho a^\dagger a] \right) \rightarrow - (\{A^\dagger A\}_L + \{A^\dagger A\}_R)\rho_v + \mathrm{E}[(\{A^\dagger A\}_L + \{A^\dagger A\}_R)\rho_v]$
In [12]:
def d1_rho_func(A, rho_vec):
n_sum = A[4] + A[5]
return - 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 \rightarrow \frac{A_LA^\dagger_R \rho_v}{\mathrm{E}[A_LA^\dagger_R \rho_v]} - \rho_v$
In [13]:
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 [14]:
result = smesolve(H, rho0, times, c_ops=[], sc_ops=c_ops, e_ops=e_ops,
ntraj=5, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='poisson', store_measurement=True)
In [15]:
plot_expectation_values([result, result_ref]);
In [16]:
fig, ax = subplots()
for m in result.measurement:
ax.step(times, m)
In [17]:
result = smesolve(H, rho0, times, c_ops=[], sc_ops=c_ops, e_ops=e_ops,
ntraj=5, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='poisson', store_measurement=True, noise=result.noise)
In [18]:
plot_expectation_values([result, result_ref]);
In [19]:
fig, ax = subplots()
for m in result.measurement:
ax.step(times, m)
In [20]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='photocurrent', store_measurement=True)
In [21]:
plot_expectation_values([result, result_ref]);
In [22]:
result = smesolve(H, rho0, times, c_ops=[], sc_ops=c_ops, e_ops=e_ops, ntraj=5, nsubsteps=100,
method='photocurrent', store_measurement=True, noise=result.noise)
In [23]:
plot_expectation_values([result, result_ref]);
In [24]:
fig, ax = subplots()
for m in result.measurement:
ax.step(times, m)
Stochastic master equation for homodyne in Milburn's formulation
$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + \gamma\mathcal{D}[a]\rho(t) dt + dW(t) \sqrt{\gamma} \mathcal{H}[a] \rho(t)$
where $\mathcal{D}$ is the standard Lindblad dissipator superoperator, and $\mathcal{H}$ is defined as above, and $dW(t)$ is a normal distributed increment with $E[dW(t)] = \sqrt{dt}$.
In QuTiP format we have:
$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + D_{1}[A]\rho(t) dt + D_{2}[A]\rho(t) dW$
where $A = \sqrt{\gamma} a$, so we can identify
$\displaystyle D_{1}[A]\rho(t) = \gamma \mathcal{D}[a]\rho(t) = \mathcal{D}[A]\rho(t)$
In [25]:
def d1_rho_func(A, rho_vec):
return A[7] * rho_vec
$\displaystyle D_{2}[A]\rho(t) = \sqrt{\gamma} \mathcal{H}[a]\rho(t) = A\rho + \rho A^\dagger - \mathrm{Tr}[A\rho + \rho A^\dagger] \rho \rightarrow (A_L + A_R^\dagger)\rho_v - \mathrm{Tr}[(A_L + A_R^\dagger)\rho_v] \rho_v$
In [26]:
def d2_rho_func(A, rho_vec):
n_sum = A[0] + A[3]
e1 = expect_rho_vec(n_sum, rho_vec)
return [n_sum * rho_vec - e1 * rho_vec]
In [27]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=5, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='normal', store_measurement=True)
In [28]:
plot_expectation_values([result, result_ref]);
In [29]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
In [30]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=5, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func,
distribution='normal', store_measurement=True, noise=result.noise)
In [31]:
plot_expectation_values([result, result_ref]);
In [32]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
In [33]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='homodyne', store_measurement=True)
In [34]:
plot_expectation_values([result, result_ref]);
In [35]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
In [36]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='homodyne', store_measurement=True, noise=result.noise)
In [37]:
plot_expectation_values([result, result_ref]);
In [38]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
Stochastic master equation for heterodyne in Milburn's formulation
$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + \gamma\mathcal{D}[a]\rho(t) dt + \frac{1}{\sqrt{2}} dW_1(t) \sqrt{\gamma} \mathcal{H}[a] \rho(t) + \frac{1}{\sqrt{2}} dW_2(t) \sqrt{\gamma} \mathcal{H}[-ia] \rho(t)$
where $\mathcal{D}$ is the standard Lindblad dissipator superoperator, and $\mathcal{H}$ is defined as above, and $dW_i(t)$ is a normal distributed increment with $E[dW_i(t)] = \sqrt{dt}$.
In QuTiP format we have:
$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + D_{1}[A]\rho(t) dt + D_{2}^{(1)}[A]\rho(t) dW_1 + D_{2}^{(2)}[A]\rho(t) dW_2$
where $A = \sqrt{\gamma} a$, so we can identify
$\displaystyle D_{1}[A]\rho = \gamma \mathcal{D}[a]\rho = \mathcal{D}[A]\rho$
In [39]:
def d1_rho_func(A, rho_vec):
return A[7] * rho_vec
$D_{2}^{(1)}[A]\rho = \frac{1}{\sqrt{2}} \sqrt{\gamma} \mathcal{H}[a] \rho = \frac{1}{\sqrt{2}} \mathcal{H}[A] \rho = \frac{1}{\sqrt{2}}(A\rho + \rho A^\dagger - \mathrm{Tr}[A\rho + \rho A^\dagger] \rho) \rightarrow \frac{1}{\sqrt{2}} \left\{(A_L + A_R^\dagger)\rho_v - \mathrm{Tr}[(A_L + A_R^\dagger)\rho_v] \rho_v\right\}$
$D_{2}^{(2)}[A]\rho = \frac{1}{\sqrt{2}} \sqrt{\gamma} \mathcal{H}[-ia] \rho = \frac{1}{\sqrt{2}} \mathcal{H}[-iA] \rho = \frac{-i}{\sqrt{2}}(A\rho - \rho A^\dagger - \mathrm{Tr}[A\rho - \rho A^\dagger] \rho) \rightarrow \frac{-i}{\sqrt{2}} \left\{(A_L - A_R^\dagger)\rho_v - \mathrm{Tr}[(A_L - A_R^\dagger)\rho_v] \rho_v\right\}$
In [40]:
def d2_rho_func(A, rho_vec):
n_sum = A[0] + A[3]
e1 = expect_rho_vec(n_sum, rho_vec)
drho1 = n_sum * rho_vec - e1 * rho_vec
n_sum = A[0] - A[3]
e1 = expect_rho_vec(n_sum, rho_vec)
drho2 = n_sum * rho_vec - e1 * rho_vec
return [1.0/sqrt(2) * drho1, -1.0j/sqrt(2) * drho2]
In [41]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=5, nsubsteps=100, d1=d1_rho_func, d2=d2_rho_func, d2_len=2,
distribution='normal', store_measurement=True)
In [42]:
plot_expectation_values([result, result_ref]);
In [43]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
In [44]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='heterodyne', store_measurement=True)
In [45]:
plot_expectation_values([result, result_ref]);
In [46]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
In [47]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='heterodyne', store_measurement=True, noise=result.noise)
In [48]:
plot_expectation_values([result, result_ref]);
In [49]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
In [50]:
from qutip.ipynbtools import version_table
version_table()
Out[50]: