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.
%pylab inline
Populating the interactive namespace from numpy and matplotlib
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.
N = 15
w0 = 0.5 * 2 * pi
times = linspace(0, 15, 150)
gamma = 0.15
a = destroy(N)
H = w0 * a.dag() * a
rho0 = fock(N, 5)
c_ops = [sqrt(gamma) * a]
e_ops = [a.dag() * a, a + a.dag()]
result_ref = mesolve(H, rho0, times, c_ops, e_ops)
plot_expectation_values(result_ref);
# 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]$
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$
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]
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)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 6.71s. Est. remaining time: 00:00:00:26. Completed: 40.0%. Elapsed time: 13.41s. Est. remaining time: 00:00:00:20. Completed: 60.0%. Elapsed time: 20.11s. Est. remaining time: 00:00:00:13. Completed: 80.0%. Elapsed time: 26.81s. Est. remaining time: 00:00:00:06. Elapsed time: 33.52s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.step(times, m)
/usr/local/lib/python3.3/dist-packages/numpy/core/numeric.py:460: ComplexWarning: Casting complex values to real discards the imaginary part return array(a, dtype, copy=False, order=order)
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)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 6.01s. Est. remaining time: 00:00:00:24. Completed: 40.0%. Elapsed time: 12.03s. Est. remaining time: 00:00:00:18. Completed: 60.0%. Elapsed time: 18.04s. Est. remaining time: 00:00:00:12. Completed: 80.0%. Elapsed time: 24.05s. Est. remaining time: 00:00:00:06. Elapsed time: 30.06s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.step(times, m)
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='photocurrent', store_measurement=True)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 5.63s. Est. remaining time: 00:00:00:22. Completed: 40.0%. Elapsed time: 11.25s. Est. remaining time: 00:00:00:16. Completed: 60.0%. Elapsed time: 16.85s. Est. remaining time: 00:00:00:11. Completed: 80.0%. Elapsed time: 22.45s. Est. remaining time: 00:00:00:05. Elapsed time: 28.07s
plot_expectation_values([result, result_ref]);
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)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 4.92s. Est. remaining time: 00:00:00:19. Completed: 40.0%. Elapsed time: 9.84s. Est. remaining time: 00:00:00:14. Completed: 60.0%. Elapsed time: 14.76s. Est. remaining time: 00:00:00:09. Completed: 80.0%. Elapsed time: 19.69s. Est. remaining time: 00:00:00:04. Elapsed time: 24.62s
plot_expectation_values([result, result_ref]);
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)$
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$
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]
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)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 4.45s. Est. remaining time: 00:00:00:17. Completed: 40.0%. Elapsed time: 8.90s. Est. remaining time: 00:00:00:13. Completed: 60.0%. Elapsed time: 13.35s. Est. remaining time: 00:00:00:08. Completed: 80.0%. Elapsed time: 17.81s. Est. remaining time: 00:00:00:04. Elapsed time: 22.26s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
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)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 4.46s. Est. remaining time: 00:00:00:17. Completed: 40.0%. Elapsed time: 8.92s. Est. remaining time: 00:00:00:13. Completed: 60.0%. Elapsed time: 13.37s. Est. remaining time: 00:00:00:08. Completed: 80.0%. Elapsed time: 17.84s. Est. remaining time: 00:00:00:04. Elapsed time: 22.30s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='homodyne', store_measurement=True)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 4.41s. Est. remaining time: 00:00:00:17. Completed: 40.0%. Elapsed time: 8.83s. Est. remaining time: 00:00:00:13. Completed: 60.0%. Elapsed time: 13.24s. Est. remaining time: 00:00:00:08. Completed: 80.0%. Elapsed time: 17.66s. Est. remaining time: 00:00:00:04. Elapsed time: 22.08s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='homodyne', store_measurement=True, noise=result.noise)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 4.42s. Est. remaining time: 00:00:00:17. Completed: 40.0%. Elapsed time: 8.83s. Est. remaining time: 00:00:00:13. Completed: 60.0%. Elapsed time: 13.25s. Est. remaining time: 00:00:00:08. Completed: 80.0%. Elapsed time: 17.65s. Est. remaining time: 00:00:00:04. Elapsed time: 22.08s
plot_expectation_values([result, result_ref]);
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$
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\}$
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]
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)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 8.24s. Est. remaining time: 00:00:00:32. Completed: 40.0%. Elapsed time: 16.47s. Est. remaining time: 00:00:00:24. Completed: 60.0%. Elapsed time: 24.68s. Est. remaining time: 00:00:00:16. Completed: 80.0%. Elapsed time: 32.91s. Est. remaining time: 00:00:00:08. Elapsed time: 41.13s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='heterodyne', store_measurement=True)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 8.13s. Est. remaining time: 00:00:00:32. Completed: 40.0%. Elapsed time: 16.27s. Est. remaining time: 00:00:00:24. Completed: 60.0%. Elapsed time: 24.42s. Est. remaining time: 00:00:00:16. Completed: 80.0%. Elapsed time: 32.57s. Est. remaining time: 00:00:00:08. Elapsed time: 40.70s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=5, nsubsteps=100,
method='heterodyne', store_measurement=True, noise=result.noise)
smesolve smesolve_generic Completed: 0.0%. Elapsed time: 0.00s. Est. remaining time: 00:00:00:00. Completed: 20.0%. Elapsed time: 8.09s. Est. remaining time: 00:00:00:32. Completed: 40.0%. Elapsed time: 16.20s. Est. remaining time: 00:00:00:24. Completed: 60.0%. Elapsed time: 24.31s. Est. remaining time: 00:00:00:16. Completed: 80.0%. Elapsed time: 32.43s. Est. remaining time: 00:00:00:08. Elapsed time: 40.54s
plot_expectation_values([result, result_ref]);
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m)
ax.plot(times, array(result.measurement).mean(axis=0), 'k', lw=5);
from qutip.ipynbtools import version_table
version_table()
Software | Version |
---|---|
Python | 3.3.1 (default, Apr 17 2013, 22:30:32) [GCC 4.7.3] |
QuTiP | 2.3.0.dev-6fe28dc |
OS | posix [linux] |
SciPy | 0.13.0.dev-7c6d92e |
IPython | 1.0.0-dev |
matplotlib | 1.4.x |
Numpy | 1.8.0.dev-b307a8a |
Cython | 0.19.1 |
Wed Jul 31 13:31:13 2013 JST |