Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson
In this notebook we test the qutip stochastic Schrödinger equation solver (ssesolve) with a few examples. The same examples are solved using the stochastic master equation solver in the notebook development-smesolve-tests.
In [1]:
%pylab inline
In [2]:
from qutip import *
Stochastic Schrödinger equation for photocurrent detection, in Milburn's formulation, is
$$ d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + \left(\frac{a}{\sqrt{\langle a^\dagger a \rangle_\psi(t)}} - 1\right) |\psi(t)\rangle dN(t) - \frac{1}{2}\gamma\left(\langle a^\dagger a \rangle_\psi(t) - a^\dagger a \right) |\psi(t)\rangle dt $$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 Schrödinger equation on the form:
$$ d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + D_{1}[c] |\psi(t)\rangle dt + D_{2}[c] |\psi(t)\rangle dW $$where $c = \sqrt{\gamma} a$ is the collapse operator including the rate of the process as a coefficient in the operator. We can identify
$$ D_{1}[c] |\psi(t)\rangle = - \frac{1}{2}\left(\langle c^\dagger c \rangle_\psi(t) - c^\dagger c \right) $$$$ D_{2}[c] |\psi(t)\rangle = \frac{c}{\sqrt{\langle c^\dagger c \rangle_\psi(t)}} - 1 $$$$dW = dN(t)$$
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]:
psi0 = fock(N, 5)
In [7]:
sc_ops = [sqrt(gamma) * a]
In [8]:
e_ops = [a.dag() * a, a + a.dag(), (-1j)*(a - a.dag())]
In [9]:
result_ref = mesolve(H, psi0, times, sc_ops, e_ops)
In [10]:
plot_expectation_values(result_ref);
In [11]:
from qutip.cy.spmatfuncs import cy_expect_psi
In [12]:
# Function sigurature:
#
# def d(A, psi):
#
# psi = wave function at the current time step
#
# A[0] = c
# A[1] = c + c.dag()
# A[2] = c - c.dag()
# A[3] = c.dag() * c
#
# where c is a collapse operator. The combinations of c's stored in A are
# precomputed before the time-evolution is started to avoid repeated
# computations.
$\displaystyle D_{1}[a, |\psi\rangle] = \frac{1}{2} ( \langle \psi|a^\dagger a |\psi\rangle - a^\dagger a )|\psi\rangle $
In [13]:
def d1_psi_func(A, psi):
return 0.5 * (-A[3] * psi + cy_expect_psi(A[3], psi, 1) * psi)
$\displaystyle D_{2}[a, |\psi\rangle] = \frac{a|\psi\rangle}{\sqrt{\langle \psi| a^\dagger a|\psi\rangle}} - |\psi\rangle $
In [14]:
def d2_psi_func(A, psi):
psi_1 = A[0] * psi
n1 = sqrt(cy_expect_psi(A[3], psi, 1))
if n1 > 0.0:
return [psi_1 / n1 - psi]
else:
return [- psi]
In [15]:
result = ssesolve(H, psi0, times, sc_ops=sc_ops, e_ops=e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, d1=d1_psi_func, d2=d2_psi_func,
distribution='poisson', store_measurement=True, m_ops=[[None]])
In [16]:
plot_expectation_values([result, result_ref]);
In [17]:
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 [18]:
result = ssesolve(H, psi0, times, sc_ops=sc_ops, e_ops=e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, m_ops=[[None]],
d1=d1_psi_func, d2=d2_psi_func, distribution='poisson',
store_measurement=True, noise=result.noise)
In [19]:
plot_expectation_values([result, result_ref]);
In [20]:
fig, ax = subplots()
for m in result.measurement:
ax.step(times, dt * m.real)
In [21]:
result = ssesolve(H, psi0, times, sc_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
method='photocurrent', store_measurement=True)
In [22]:
plot_expectation_values([result, result_ref]);
In [23]:
result = ssesolve(H, psi0, times, sc_ops=sc_ops, e_ops=e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
method='photocurrent', store_measurement=True, noise=result.noise)
In [24]:
plot_expectation_values([result, result_ref]);
In [25]:
fig, ax = subplots()
for m in result.measurement:
ax.step(times, dt * m.real)
In [26]:
ntraj = 50
nsubsteps = 200
In [27]:
H = w0 * a.dag() * a + A * (a + a.dag())
In [28]:
result_ref = mesolve(H, psi0, times, sc_ops, e_ops)
Stochastic master equation for homodyne can be written as (see Gardiner and Zoller, Quantum Noise)
$$ d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + \frac{1}{2}\gamma \left( \langle a + a^\dagger \rangle_\psi a - a^\dagger a - \frac{1}{4}\langle a + a^\dagger \rangle_\psi^2 \right)|\psi(t)\rangle dt + \sqrt{\gamma}\left( a - \frac{1}{2} \langle a + a^\dagger \rangle_\psi \right)|\psi(t)\rangle dW $$where $dW(t)$ is a normal distributed increment with $E[dW(t)] = \sqrt{dt}$.
In QuTiP format we have:
$$ d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + D_{1}[c, |\psi(t)\rangle] dt + D_{2}[c, |\psi(t)\rangle] dW $$where $c = \sqrt{\gamma} a$, so we can identify
In [29]:
def d1_psi_func(A, psi):
e_x = cy_expect_psi(A[1], psi, 0)
return 0.5 * (e_x * (A[0] * psi) - (A[3] * psi) - 0.25 * e_x ** 2 * psi)
In [30]:
def d2_psi_func(A, psi):
e_x = cy_expect_psi(A[1], psi, 0)
return [(A[0] * psi) - 0.5 * e_x * psi]
In [31]:
result = ssesolve(H, psi0, times, sc_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
d1=d1_psi_func, d2=d2_psi_func, distribution='normal',
m_ops=[[(a + a.dag())]], dW_factors=[1/sqrt(gamma)],
store_measurement=True)
In [32]:
plot_expectation_values([result, result_ref]);
In [33]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real, 'b', alpha=0.1)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2)
ax.plot(times, result_ref.expect[1], 'k', lw=2)
ax.set_ylim(-15, 15);
In [34]:
result = ssesolve(H, psi0, times, sc_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
d1=d1_psi_func, d2=d2_psi_func, distribution='normal',
m_ops=[[(a + a.dag())]], dW_factors=[1/sqrt(gamma)],
store_measurement=True, noise=result.noise)
In [35]:
plot_expectation_values([result, result_ref]);
In [36]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real, 'b', alpha=0.1)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2)
ax.plot(times, result_ref.expect[1], 'k', lw=2);
ax.set_ylim(-20, 20);
In [37]:
result = ssesolve(H, psi0, times, sc_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
method='homodyne', store_measurement=True)
In [38]:
plot_expectation_values([result, result_ref]);
In [39]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real / sqrt(gamma), 'b', alpha=0.1)
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)
ax.set_ylim(-20, 20);
In [40]:
result = ssesolve(H, psi0, times, sc_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
method='homodyne', store_measurement=True, noise=result.noise)
In [41]:
plot_expectation_values([result, result_ref]);
In [42]:
fig, ax = subplots()
for m in result.measurement:
ax.plot(times, m[:, 0].real / sqrt(gamma), 'b', alpha=0.1)
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)
ax.set_ylim(-20, 20);
Stochastic Schrödinger equation for heterodyne detection can be written as
$$ d|\psi(t)\rangle = -iH|\psi(t)\rangle dt -\frac{1}{2}\gamma\left(a^\dagger a - \langle a^\dagger \rangle a + \frac{1}{2}\langle a \rangle\langle a^\dagger \rangle)\right)|\psi(t)\rangle dt + \sqrt{\gamma/2}\left( a - \frac{1}{2} \langle a + a^\dagger \rangle_\psi \right)|\psi(t)\rangle dW_1 -i \sqrt{\gamma/2}\left( a - \frac{1}{2} \langle a - a^\dagger \rangle_\psi \right)|\psi(t)\rangle dW_2 $$where $dW_i(t)$ is a normal distributed increment with $E[dW_i(t)] = \sqrt{dt}$.
In QuTiP format we have:
$$ d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + D_{1}[c, |\psi(t)\rangle] dt + \sum_n D_{2}^{(n)}[c, |\psi(t)\rangle] dW_n $$where $c = \sqrt{\gamma} a$, so we can identify
In [43]:
def d1_psi_func(A, psi):
e_c = cy_expect_psi(A[0], psi, 0)
e_cd = cy_expect_psi(A[0].T.conj(), psi, 0)
return (-0.5 * (A[3] * psi) + 0.5 * e_cd * (A[0] * psi) - 0.25 * e_c * e_cd * psi)
In [44]:
def d2_psi_func(A, psi):
x = 0.5 * cy_expect_psi(A[1], psi, 0)
y = 0.5 * cy_expect_psi(A[2], psi, 0)
d2_1 = np.sqrt(0.5) * ((A[0] * psi) - x * psi)
d2_2 = -1.0j * np.sqrt(0.5) * ((A[0] * psi) - y * psi)
return [d2_1, d2_2]
In [45]:
result = ssesolve(H, psi0, times, sc_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps,
d1=d1_psi_func, d2=d2_psi_func, d2_len=2, distribution='normal',
m_ops=[[(a + a.dag()), (-1j)*(a - a.dag())]],
dW_factors=[2/sqrt(gamma), 2/sqrt(gamma)],
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.plot(times, result_ref.expect[1], 'k', lw=2);
ax.plot(times, result_ref.expect[2], 'k', lw=2);
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);
In [48]:
result = ssesolve(H, psi0, times, sc_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, result_ref.expect[1], 'k', lw=2);
ax.plot(times, result_ref.expect[2], 'k', lw=2);
#ax.set_ylim(-15, 15)
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);
In [51]:
result = ssesolve(H, psi0, times, sc_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, result_ref.expect[1], 'k', lw=2);
ax.plot(times, result_ref.expect[2], 'k', lw=2);
ax.set_ylim(-25, 25)
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);
In [54]:
from qutip.ipynbtools import version_table
version_table()
Out[54]: