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}[\frac{1}{2}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 = \frac{1}{2}(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)dt$.
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) = - \frac{1}{2}\gamma \mathcal{H}[a^\dagger a] \rho(t) = -\gamma \frac{1}{2}\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}[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 = 10
w0 = 0.5 * 2 * pi
times = linspace(0, 15, 150)
dt = times[1] - times[0]
gamma = 0.25
A = 2.5
ntraj = 50
nsubsteps = 100
In [4]:
a = destroy(N)
x = a + a.dag()
In [5]:
H = w0 * a.dag() * a
In [6]:
#rho0 = coherent(N, 5)
rho0 = fock(N, 5)
In [7]:
c_ops = [sqrt(gamma) * a]
In [8]:
e_ops = [a.dag() * a, x]
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 \frac{1}{2}\left( a^\dagger a\rho + \rho a^\dagger a - \mathrm{Tr}[a^\dagger a\rho + \rho a^\dagger a] \right) \rightarrow - \frac{1}{2}(\{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 0.5 * (expect_rho_vec(n_sum, rho_vec, False) * rho_vec - n_sum * 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, False) + 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=ntraj, nsubsteps=nsubsteps, d1=d1_rho_func, d2=d2_rho_func, m_ops=[[None]],
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, dt * m.real)
Note that we multiply the measurement records with dt to show the count events instead of the count rate.
In [17]:
result = smesolve(H, rho0, times, c_ops=[], sc_ops=c_ops, e_ops=e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
d1=d1_rho_func, d2=d2_rho_func, distribution='poisson', m_ops=[[None]],
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, dt * m.real)
In [20]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
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=ntraj, nsubsteps=nsubsteps,
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, dt * m.real)
In [25]:
H = w0 * a.dag() * a + A * (a + a.dag())
In [26]:
result_ref = mesolve(H, rho0, times, c_ops, e_ops)
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 [27]:
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 [28]:
def d2_rho_func(A, rho_vec):
n_sum = A[0] + A[3]
e1 = expect_rho_vec(n_sum, rho_vec, False)
return [n_sum * rho_vec - e1 * rho_vec]
In [29]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
d1=d1_rho_func, d2=d2_rho_func, distribution='normal',
m_ops=[[(a + a.dag())]], dW_factors=[1/sqrt(gamma)],
store_measurement=True)
In [30]:
plot_expectation_values([result, result_ref]);
In [31]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real, 'b', alpha=0.025)
ax.plot(times, result_ref.expect[1], 'k', lw=2);
ax.set_ylim(-15, 15)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2);
In [32]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, d1=d1_rho_func, d2=d2_rho_func,
m_ops=[[(a + a.dag())]], dW_factors=[1/sqrt(gamma)],
distribution='normal', store_measurement=True, noise=result.noise)
In [33]:
plot_expectation_values([result, result_ref]);
In [34]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real, 'b', alpha=0.025)
ax.plot(times, result_ref.expect[1], 'k', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2);
In [35]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=100,
method='homodyne', store_measurement=True)
In [36]:
plot_expectation_values([result, result_ref]);
In [37]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real / sqrt(gamma), 'b', alpha=0.025)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real / sqrt(gamma), 'r', lw=2);
ax.plot(times, result_ref.expect[1], 'k', lw=2);
In [38]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
method='homodyne', store_measurement=True, noise=result.noise)
In [39]:
plot_expectation_values([result, result_ref]);
In [40]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real / sqrt(gamma), 'b', alpha=0.025)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real / sqrt(gamma), 'r', lw=2)
ax.plot(times, result_ref.expect[1], 'k', lw=2);
In [41]:
e_ops = [a.dag() * a, a + a.dag(), -1j * (a - a.dag())]
In [42]:
result_ref = mesolve(H, rho0, times, c_ops, e_ops)
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 [43]:
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 [44]:
def d2_rho_func(A, rho_vec):
n_sum = A[0] + A[3]
e1 = expect_rho_vec(n_sum, rho_vec, False)
drho1 = n_sum * rho_vec - e1 * rho_vec
n_sum = A[0] - A[3]
e1 = expect_rho_vec(n_sum, rho_vec, False)
drho2 = n_sum * rho_vec - e1 * rho_vec
return [1.0/sqrt(2) * drho1, -1.0j/sqrt(2) * drho2]
In [45]:
result = smesolve(H, rho0, times, [], c_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, d1=d1_rho_func, d2=d2_rho_func, d2_len=2,
m_ops=[[(a + a.dag()), (-1j)*(a - a.dag())]],
dW_factors=[2/sqrt(gamma), 2/sqrt(gamma)],
distribution='normal', store_measurement=True)
In [46]:
plot_expectation_values([result, result_ref]);
In [47]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0, 0].real, 'r', alpha=0.025)
ax.plot(times, m[:, 0, 1].real, 'b', alpha=0.025)
ax.set_ylim(-20, 20)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,0].real, 'r', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,1].real, 'b', lw=2);
ax.plot(times, result_ref.expect[1], 'k', lw=2);
ax.plot(times, result_ref.expect[2], 'k', lw=2);
In [48]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
method='heterodyne', store_measurement=True)
In [49]:
plot_expectation_values([result, result_ref]);
In [50]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0, 0].real / sqrt(gamma), 'r', alpha=0.025)
ax.plot(times, m[:, 0, 1].real / sqrt(gamma), 'b', alpha=0.025)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,0].real / sqrt(gamma), 'r', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,1].real / sqrt(gamma), 'b', lw=2);
ax.plot(times, result_ref.expect[1], 'k', lw=2);
ax.plot(times, result_ref.expect[2], 'k', lw=2);
In [51]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
method='heterodyne', store_measurement=True, noise=result.noise)
In [52]:
plot_expectation_values([result, result_ref]);
In [53]:
fig, ax = subplots(figsize=(12,6))
for m in result.measurement:
ax.plot(times, m[:, 0, 0].real / sqrt(gamma), 'r', alpha=0.025)
ax.plot(times, m[:, 0, 1].real / sqrt(gamma), 'b', alpha=0.025)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,0].real / sqrt(gamma), 'r', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,1].real / sqrt(gamma), 'b', lw=2);
ax.plot(times, result_ref.expect[1], 'k', lw=2);
ax.plot(times, result_ref.expect[2], 'k', lw=2);
ax.axis('tight')
ax.set_ylim(-25, 25);
In [54]:
from qutip.ipynbtools import version_table
version_table()
Out[54]: