Two implementations of heterodyne: direct heterodyne and as two homodyne measurements

Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson

Warning: WIP: Requires latest development version of QuTiP.


In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
from qutip import *

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 [3]:
N = 15
w0 = 1.0 * 2 * pi
A = 0.1 * 2 * pi
times = linspace(0, 15, 150)
gamma = 0.25

ntraj = 200
nsubsteps = 250

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, sqrt(5.0), method='analytic')
c_ops = [sqrt(gamma) * a]
e_ops = [a.dag() * a, x, y]

In [4]:
result_ref = mesolve(H, rho0, times, c_ops, e_ops)

In [5]:
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 [6]:
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 [7]:
def d2_rho_func(A, rho_vec):

    B = A[0] + A[3]
    e1 = expect_rho_vec(B, rho_vec)
    drho1 = B * rho_vec - e1 * rho_vec

    B = A[0] - A[3]
    e1 = expect_rho_vec(B, rho_vec)
    drho2 = B * rho_vec - e1 * rho_vec

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

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

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

$J_y(t) = \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), sqrt(2)].


In [8]:
result = smesolve(H, rho0, times, [], c_ops, e_ops, 
                  ntraj=ntraj, nsubsteps=nsubsteps,
                  d1=d1_rho_func, d2=d2_rho_func, d2_len=2, dW_factors=[sqrt(2), sqrt(2)], m_ops=[[x, y]],
                  distribution='normal', store_measurement=True)


Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time: 461.21s. Est. remaining time: 00:01:09:10.
Completed: 20.0%. Elapsed time: 927.10s. Est. remaining time: 00:01:01:48.
Completed: 30.0%. Elapsed time: 1381.08s. Est. remaining time: 00:00:53:42.
Completed: 40.0%. Elapsed time: 1836.64s. Est. remaining time: 00:00:45:54.
Completed: 50.0%. Elapsed time: 2287.43s. Est. remaining time: 00:00:38:07.
Completed: 60.0%. Elapsed time: 2738.17s. Est. remaining time: 00:00:30:25.
Completed: 70.0%. Elapsed time: 3188.93s. Est. remaining time: 00:00:22:46.
Completed: 80.0%. Elapsed time: 3639.65s. Est. remaining time: 00:00:15:09.
Completed: 90.0%. Elapsed time: 4090.24s. Est. remaining time: 00:00:07:34.
Elapsed time: 4541.58s

In [9]:
plot_expectation_values([result, result_ref]);



In [10]:
fig, ax = subplots(figsize=(8,4))

for m in result.measurement:
    ax.plot(times, m[:, 0, 0].real, 'b', alpha=0.05)
    ax.plot(times, m[:, 0, 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_xlim(0, times.max())
ax.set_xlabel('time', fontsize=12)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,0].real, 'k', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:,0,1].real, 'k', lw=2);


Heterodyne implementation #2: using two homodyne measurements

We can also write the heterodyne equation as

$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + \frac{1}{2}\gamma\mathcal{D}[a]\rho(t) dt + \frac{1}{\sqrt{2}} dW_1(t) \sqrt{\gamma} \mathcal{H}[a] \rho(t) + \frac{1}{2}\gamma\mathcal{D}[a]\rho(t) dt + \frac{1}{\sqrt{2}} dW_2(t) \sqrt{\gamma} \mathcal{H}[-ia] \rho(t)$

And using the QuTiP format for two stochastic collapse operators, we have:

$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + D_{1}[A_1]\rho(t) dt + D_{2}[A_1]\rho(t) dW_1 + D_{1}[A_2]\rho(t) dt + D_{2}[A_2]\rho(t) dW_2$

so we can also identify

$\displaystyle D_{1}[A_1]\rho = \frac{1}{2}\gamma \mathcal{D}[a]\rho = \mathcal{D}[\sqrt{\gamma}a/\sqrt{2}]\rho = \mathcal{D}[A_1]\rho$

$\displaystyle D_{1}[A_2]\rho = \frac{1}{2}\gamma \mathcal{D}[a]\rho = \mathcal{D}[-i\sqrt{\gamma}a/\sqrt{2}]\rho = \mathcal{D}[A_2]\rho$

$D_{2}[A_1]\rho = \frac{1}{\sqrt{2}} \sqrt{\gamma} \mathcal{H}[a] \rho = \mathcal{H}[A_1] \rho$

$D_{2}[A_2]\rho = \frac{1}{\sqrt{2}} \sqrt{\gamma} \mathcal{H}[-ia] \rho = \mathcal{H}[A_2] \rho $

where $A_1 = \sqrt{\gamma} a / \sqrt{2}$ and $A_2 = -i \sqrt{\gamma} a / \sqrt{2}$.

In summary we have

$\displaystyle d\rho(t) = -i[H, \rho(t)]dt + \sum_i\left\{\mathcal{D}[A_i]\rho(t) dt + \mathcal{H}[A_i]\rho(t) dW_i\right\}$

which is a simultaneous homodyne detection with $A_1 = \sqrt{\gamma}a/\sqrt{2}$ and $A_2 = -i\sqrt{\gamma}a/\sqrt{2}$

Here the two heterodyne currents for the $x$ and $y$ quadratures are

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

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

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

In qutip we can use the predefined homodyne solver for solving this problem.


In [11]:
result = smesolve(H, rho0, times, [], [sqrt(gamma/2) * a, -1.0j * sqrt(gamma/2) * a],
                  e_ops, ntraj=ntraj, nsubsteps=nsubsteps, m_ops=[[x], [y]], dW_factors=[sqrt(2)],
                  method='homodyne', store_measurement=True)


Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time: 455.67s. Est. remaining time: 00:01:08:21.
Completed: 20.0%. Elapsed time: 908.73s. Est. remaining time: 00:01:00:34.
Completed: 30.0%. Elapsed time: 1362.95s. Est. remaining time: 00:00:53:00.
Completed: 40.0%. Elapsed time: 1816.76s. Est. remaining time: 00:00:45:25.
Completed: 50.0%. Elapsed time: 2269.20s. Est. remaining time: 00:00:37:49.
Completed: 60.0%. Elapsed time: 2722.45s. Est. remaining time: 00:00:30:14.
Completed: 70.0%. Elapsed time: 3176.23s. Est. remaining time: 00:00:22:41.
Completed: 80.0%. Elapsed time: 3630.57s. Est. remaining time: 00:00:15:07.
Completed: 90.0%. Elapsed time: 4084.32s. Est. remaining time: 00:00:07:33.
Elapsed time: 4537.40s

In [12]:
plot_expectation_values([result, result_ref]);



In [13]:
fig, ax = 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_xlim(0, times.max())
ax.set_ylim(-15, 15)
ax.set_xlabel('time', fontsize=12)
ax.plot(times, array(result.measurement).mean(axis=0)[:,0].real, 'k', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:,1].real, 'k', lw=2);


Implementation #3: builtin function for heterodyne


In [14]:
result = smesolve(H, rho0, times, [], [sqrt(gamma) * a],
                  e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
                  method='heterodyne', store_measurement=True)


Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time: 448.84s. Est. remaining time: 00:01:07:19.
Completed: 20.0%. Elapsed time: 896.30s. Est. remaining time: 00:00:59:45.
Completed: 30.0%. Elapsed time: 1346.22s. Est. remaining time: 00:00:52:21.
Completed: 40.0%. Elapsed time: 1796.25s. Est. remaining time: 00:00:44:54.
Completed: 50.0%. Elapsed time: 2243.91s. Est. remaining time: 00:00:37:23.
Completed: 60.0%. Elapsed time: 2691.39s. Est. remaining time: 00:00:29:54.
Completed: 70.0%. Elapsed time: 3138.75s. Est. remaining time: 00:00:22:25.
Completed: 80.0%. Elapsed time: 3583.92s. Est. remaining time: 00:00:14:55.
Completed: 90.0%. Elapsed time: 4027.13s. Est. remaining time: 00:00:07:27.
Elapsed time: 4473.34s

In [15]:
plot_expectation_values([result, result_ref]);



In [16]:
fig, ax = subplots(figsize=(8,4))

for m in result.measurement:
    ax.plot(times, m[:, 0, 0].real, 'b', alpha=0.05)
    ax.plot(times, m[:, 0, 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_xlim(0, times.max())
ax.set_ylim(-15, 15)
ax.set_xlabel('time', fontsize=12)
ax.plot(times, array(result.measurement).mean(axis=0)[:, 0, 0].real, 'k', lw=2);
ax.plot(times, array(result.measurement).mean(axis=0)[:, 0, 1].real, 'k', lw=2);


Versions


In [17]:
from qutip.ipynbtools import version_table

version_table()


Out[17]:
SoftwareVersion
Python3.3.1 (default, Apr 17 2013, 22:30:32) [GCC 4.7.3]
OSposix [linux]
SciPy0.13.0.dev-7c6d92e
QuTiP2.3.0.dev-2283478
matplotlib1.4.x
Cython0.19.1
IPython2.0.0-dev
Numpy1.8.0.dev-b307a8a
Thu Sep 26 21:20:52 2013 JST