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]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
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 * np.pi
times = np.linspace(0, 15, 150)
dt = times[1] - times[0]
gamma = 0.25
A = 2.5
ntraj = 50
nsubsteps = 50
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 = [np.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);
$\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]$
$\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 [11]:
result = photocurrent_mesolve(H, rho0, times, c_ops=[], sc_ops=c_ops, e_ops=e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
store_measurement=True, noise=1234)
In [12]:
plot_expectation_values([result, result_ref]);
In [13]:
for m in result.measurement:
plt.step(times, dt * m.real)
In [14]:
result = photocurrent_mesolve(H, rho0, times, c_ops=[], sc_ops=c_ops, e_ops=e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, store_measurement=True, noise=1234)
In [15]:
plot_expectation_values([result, result_ref]);
In [16]:
for m in result.measurement:
plt.step(times, dt * m.real)
In [17]:
H = w0 * a.dag() * a + A * (a + a.dag())
In [18]:
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 [19]:
L = liouvillian(H, c_ops=c_ops).data
def d1_rho_func(t, rho_vec):
return cy.spmv(L, 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 [20]:
n_sum = spre(c_ops[0]) + spost(c_ops[0].dag())
n_sum_data = n_sum.data
def d2_rho_func(t, rho_vec):
e1 = cy.cy_expect_rho_vec(n_sum_data, rho_vec, False)
out = np.zeros((1,len(rho_vec)),dtype=complex)
out += cy.spmv(n_sum_data, rho_vec) - e1 * rho_vec
return out
In [21]:
result = general_stochastic(ket2dm(rho0), times, d1=d1_rho_func, d2=d2_rho_func,
e_ops=[spre(op) for op in e_ops], ntraj=ntraj, solver="platen",
m_ops=[spre(a + a.dag())], dW_factors=[1/np.sqrt(gamma)],
nsubsteps=nsubsteps, store_measurement=True, map_func=parallel_map)
In [22]:
plot_expectation_values([result, result_ref]);
In [23]:
for m in result.measurement:
plt.plot(times, m[:, 0].real, 'b', alpha=0.025)
plt.plot(times, result_ref.expect[1], 'k', lw=2);
plt.ylim(-15, 15)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2);
In [24]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps, solver="pc-euler",
method='homodyne', store_measurement=True)
In [25]:
plot_expectation_values([result, result_ref]);
In [26]:
for m in result.measurement:
plt.plot(times, m[:, 0].real / np.sqrt(gamma), 'b', alpha=0.025)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0].real / np.sqrt(gamma), 'r', lw=2);
plt.plot(times, result_ref.expect[1], 'k', lw=2)
Out[26]:
In [27]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps, solver="pc-euler",
method='homodyne', store_measurement=True, noise=result.noise)
In [28]:
plot_expectation_values([result, result_ref]);
In [29]:
for m in result.measurement:
plt.plot(times, m[:, 0].real / np.sqrt(gamma), 'b', alpha=0.025)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0].real / np.sqrt(gamma), 'r', lw=2);
plt.plot(times, result_ref.expect[1], 'k', lw=2)
Out[29]:
In [30]:
e_ops = [a.dag() * a, a + a.dag(), -1j * (a - a.dag())]
In [31]:
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 [32]:
#def d1_rho_func(A, rho_vec):
# return A[7] * rho_vec
L = liouvillian(H, c_ops=c_ops).data
def d1_rho_func(t, rho_vec):
return cy.spmv(L, 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 [33]:
n_sump = spre(c_ops[0]) + spost(c_ops[0].dag())
n_sump_data = n_sump.data/np.sqrt(2)
n_summ = spre(c_ops[0]) - spost(c_ops[0].dag())
n_summ_data = -1.0j*n_summ.data/np.sqrt(2)
def d2_rho_func(A, rho_vec):
out = np.zeros((2,len(rho_vec)),dtype=complex)
e1 = cy.cy_expect_rho_vec(n_sump_data, rho_vec, False)
out[0,:] += cy.spmv(n_sump_data, rho_vec) - e1 * rho_vec
e1 = cy.cy_expect_rho_vec(n_summ_data, rho_vec, False)
out[1,:] += cy.spmv(n_summ_data, rho_vec) - e1 * rho_vec
return out
#def d2_rho_func(t, rho_vec):
# e1 = cy.cy_expect_rho_vec(n_sum_data, rho_vec, False)
# out = np.zeros((1,len(rho_vec)),dtype=complex)
# out += cy.spmv(n_sum_data, rho_vec) - e1 * rho_vec
# return out
In [34]:
result = general_stochastic(ket2dm(rho0), times, d1=d1_rho_func, d2=d2_rho_func,
e_ops=[spre(op) for op in e_ops], solver="platen", # order=1
ntraj=ntraj, nsubsteps=nsubsteps, len_d2=2,
m_ops=[spre(a + a.dag()), (-1j)*spre(a - a.dag())],
dW_factors=[2/np.sqrt(gamma), 2/np.sqrt(gamma)],
store_measurement=True, map_func=parallel_map)
In [35]:
plot_expectation_values([result, result_ref])
Out[35]:
In [36]:
for m in result.measurement:
plt.plot(times, m[:, 0].real, 'r', alpha=0.025)
plt.plot(times, m[:, 1].real, 'b', alpha=0.025)
plt.ylim(-20, 20)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2);
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,1].real, 'b', lw=2);
plt.plot(times, result_ref.expect[1], 'k', lw=2);
plt.plot(times, result_ref.expect[2], 'k', lw=2);
In [37]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps, solver="milstein", # order=1
method='heterodyne', store_measurement=True)
In [38]:
plot_expectation_values([result, result_ref]);
In [39]:
for m in result.measurement:
plt.plot(times, m[:, 0, 0].real / np.sqrt(gamma), 'r', alpha=0.025)
plt.plot(times, m[:, 0, 1].real / np.sqrt(gamma), 'b', alpha=0.025)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,0].real / np.sqrt(gamma), 'r', lw=2);
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,1].real / np.sqrt(gamma), 'b', lw=2);
plt.plot(times, result_ref.expect[1], 'k', lw=2);
plt.plot(times, result_ref.expect[2], 'k', lw=2);
In [40]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps, solver="milstein", # order=1
method='heterodyne', store_measurement=True, noise=result.noise)
In [41]:
plot_expectation_values([result, result_ref]);
In [42]:
for m in result.measurement:
plt.plot(times, m[:, 0, 0].real / np.sqrt(gamma), 'r', alpha=0.025)
plt.plot(times, m[:, 0, 1].real / np.sqrt(gamma), 'b', alpha=0.025)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,0].real / np.sqrt(gamma), 'r', lw=2);
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,1].real / np.sqrt(gamma), 'b', lw=2);
plt.plot(times, result_ref.expect[1], 'k', lw=2);
plt.plot(times, result_ref.expect[2], 'k', lw=2);
plt.axis('tight')
plt.ylim(-25, 25);
In [43]:
from qutip.ipynbtools import version_table
version_table()
Out[43]:
In [ ]: