Stochastic Solver: Heterodyne Detection¶

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from qutip import (Options, about, coherent, destroy, general_stochastic,
                   ket2dm, lindblad_dissipator, liouvillian, mesolve,
                   parallel_map, plot_expectation_values, smesolve, spost,
                   spre, stochastic_solvers)
from qutip.expect import expect_rho_vec

%matplotlib inline

Introduction¶

Homodyne and hetrodyne detection are techniques for measuring the quadratures of a field using photocounters. Homodyne detection (on-resonant) measures one quadrature and with heterodyne detection (off-resonant) both quadratures can be detected simulateously.

The evolution of a quantum system that is coupled to a field that is monitored with homodyne and heterodyne detector can be described with stochastic master equations. This notebook compares two different ways to implement the heterodyne detection stochastic master equation in QuTiP.

Deterministic reference¶

In [2]:

N = 15
w0 = 1.0 * 2 * np.pi
A = 0.1 * 2 * np.pi
times = np.linspace(0, 15, 201)
gamma = 0.25

ntraj = 50
nsubsteps = 50

a = destroy(N)
x = a + a.dag()
y = -1.0j * (a - a.dag())

H = w0 * a.dag() * a + A * (a + a.dag())

rho0 = coherent(N, np.sqrt(5.0), method="analytic")
c_ops = [np.sqrt(gamma) * a]
e_ops = [a.dag() * a, x, y]

In [3]:

result_ref = mesolve(H, rho0, times, c_ops, e_ops)

In [4]:

plot_expectation_values(result_ref);

Heterodyne implementation #1¶

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 [5]:

L = liouvillian(H)
D = lindblad_dissipator(c_ops[0])
d1_operator = L + D


def d1_rho_func(t, rho_vec):
    return d1_operator * 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 [6]:

B1 = spre(c_ops[0]) + spost(c_ops[0].dag())
B2 = spre(c_ops[0]) + spost(c_ops[0].dag())


def d2_rho_func(t, rho_vec):
    e1 = expect_rho_vec(B1.data, rho_vec, False)
    drho1 = B1 * rho_vec - e1 * rho_vec

    e1 = expect_rho_vec(B2.data, rho_vec, False)
    drho2 = B2 * rho_vec - e1 * rho_vec

    return np.vstack([1.0 / np.sqrt(2) * drho1, -1.0j / np.sqrt(2) * drho2])

The heterodyne currents for the $x$ and $y$ quadratures are

$J_x(t) = \sqrt{\gamma}\left<x\right> + \sqrt{2} \xi(t)$

$J_y(t) = \sqrt{\gamma}\left<y\right> + \sqrt{2} \xi(t)$

where $\xi(t) = \frac{dW}{dt}$.

In qutip we define these measurement operators using the m_ops = [[x, y]] and the coefficients to the noise terms dW_factor = [sqrt(2/gamma), sqrt(2/gamma)].

In [7]:

result = general_stochastic(
    ket2dm(rho0),
    times,
    d1_rho_func,
    d2_rho_func,
    e_ops=[spre(op) for op in e_ops],
    len_d2=2,
    ntraj=ntraj,
    nsubsteps=nsubsteps,
    solver="platen",
    dW_factors=[np.sqrt(2 / gamma), np.sqrt(2 / gamma)],
    m_ops=[spre(x), spre(y)],
    store_measurement=True,
    map_func=parallel_map,
)

10.0%. Run time:   7.87s. Est. time left: 00:00:01:10
20.0%. Run time:  14.13s. Est. time left: 00:00:00:56
30.0%. Run time:  20.97s. Est. time left: 00:00:00:48
40.0%. Run time:  27.39s. Est. time left: 00:00:00:41
50.0%. Run time:  33.87s. Est. time left: 00:00:00:33
60.0%. Run time:  40.45s. Est. time left: 00:00:00:26
70.0%. Run time:  47.03s. Est. time left: 00:00:00:20
80.0%. Run time:  53.59s. Est. time left: 00:00:00:13
90.0%. Run time:  60.35s. Est. time left: 00:00:00:06
100.0%. Run time:  66.84s. Est. time left: 00:00:00:00
Total run time:  66.85s

In [8]:

plot_expectation_values([result, result_ref]);

In [9]:

fig, ax = plt.subplots(figsize=(8, 4))

for m in result.measurement:
    ax.plot(times, m[:, 0].real, "b", alpha=0.05)
    ax.plot(times, m[:, 1].real, "r", alpha=0.05)

ax.plot(times, result_ref.expect[1], "b", lw=2)
ax.plot(times, result_ref.expect[2], "r", lw=2)

ax.set_ylim(-10, 10)
ax.set_xlim(0, times.max())
ax.set_xlabel("time", fontsize=12)
ax.plot(times, np.array(result.measurement).mean(axis=0)[:, 0].real, "k", lw=2)
ax.plot(times, np.array(result.measurement).mean(axis=0)[:, 1].real, "k", lw=2);