Stochastic Solver¶

Homodyne detection¶

Homodyne detection is an extension of the photocurrent method where the output is mixed with a strong external source allowing to get information about the phase of the system. With this method, the resulting detection rate depends is

\[ \begin{align}\begin{aligned}:label: jump_rate\\\tau = tr \left((\gamma^2 + \gamma (C+C^\dag) + C^\dag C)\rho \right)\end{aligned}\end{align} \]

With \(\gamma\), the strength of the external beam and \(C\) the collapse operator. When the beam is very strong \((\gamma >> C^\dag C)\), the rate becomes a constant term plus a term proportional to the quadrature of the system.

Closed system¶

In closed systems, the resulting stochastic differential equation is

(1)¶\[\delta \psi(t) = - i H \psi(t) \delta t - \sum_n \left( \frac{C_n^{+} C_n}{2} -\frac{e_n}{2} C_n + \frac{e_n^2}{8} \right) \psi \delta t + \sum_n \left( C_n - \frac{e_n}{2} \right) \psi \delta \omega\]

with

\[ \begin{align}\begin{aligned}:label: jump_rate\\e_n = \left<\psi(t)|C_n + C_n^{+}|\psi(t)\right>\end{aligned}\end{align} \]

Here \(\delta \omega\) is a Wiener increment.

In QuTiP, this is available with the function ssesolve.

In [1]: times = np.linspace(0.0, 10.0, 201)

In [2]: psi0 = tensor(fock(2, 0), fock(10, 5))

In [3]: a  = tensor(qeye(2), destroy(10))

In [4]: sm = tensor(destroy(2), qeye(10))

In [5]: H = 2 * np.pi * a.dag() * a + 2 * np.pi * sm.dag() * sm + 2 * np.pi * 0.25 * (sm * a.dag() + sm.dag() * a)

In [6]: data = ssesolve(H, psi0, times, sc_ops=[np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm], method="homodyne")
Total run time:   0.01s

In [7]: figure()
Out[7]: <Figure size 640x480 with 0 Axes>

In [8]: plot(times, data.expect[0], times, data.expect[1])
Out[8]: 
[<matplotlib.lines.Line2D at 0x1a25a75c18>,
 <matplotlib.lines.Line2D at 0x1a25a755c0>]

In [9]: title('Homodyne time evolution')
Out[9]: Text(0.5,1,'Homodyne time evolution')

In [10]: xlabel('Time')
Out[10]: Text(0.5,0,'Time')

In [11]: ylabel('Expectation values')
Out[11]: Text(0,0.5,'Expectation values')

In [12]: legend(("cavity photon number", "atom excitation probability"))
Out[12]: <matplotlib.legend.Legend at 0x1a263e9860>

In [13]: show()

Open system¶

In open systems, 2 types of collapse operators are considered, \(S_i\) represent the dissipation in the environment, \(C_i\) are monitored operators. The deterministic part of the evolution is the liouvillian with both types of collapses

(2)¶\[L(\rho(t)) = - i[H(t),\rho(t)] + \sum_n D(S_n, \rho) + \sum_i D(C_i, \rho),\]

with

(3)¶\[D(C, \rho) = \frac{1}{2} \left[2 C \rho(t) C^{+} - \rho(t) C^{+} C - C^{+} C \rho(t) \right].\]

The stochastic part is given by

(4)¶\[d_2 = \left(C \rho(t) + \rho(t) C^{+} - \rm{tr}\left(C \times \rho + \rho \times C^{+} \right)\rho(t) \right),\]

resulting in the stochastic differential equation

(5)¶\[\delta \rho(t) = L(\rho(t)) \delta t + d_2 \delta \omega\]

The function smesolve covert these cases in QuTiP.

Heterodyne detection¶

With heterodyne detection, two measurements are made in order to obtain information about 2 orthogonal quadratures at once.