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

(1)\[\tau = \tr \left((\gamma^2 + \gamma (C+C^\dag) + C^\dag C)\rho \right)\]

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

(2)\[\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\]


(3)\[e_n = \left<\psi(t)|C_n + C_n^{+}|\psi(t)\right>\]

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

In QuTiP, this is available with the function ssesolve.

times = np.linspace(0.0, 10.0, 201)
psi0 = tensor(fock(2, 0), fock(10, 5))
a  = tensor(qeye(2), destroy(10))
sm = tensor(destroy(2), qeye(10))

H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a)
data = ssesolve(H, psi0, times, sc_ops=[np.sqrt(0.1) * a], e_ops=[a.dag()*a, sm.dag()*sm], method="homodyne")

plt.plot(times, data.expect[0], times, data.expect[1])
plt.title('Homodyne time evolution')
plt.ylabel('Expectation values')
plt.legend(("cavity photon number", "atom excitation probability"))

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

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


(5)\[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

(6)\[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

(7)\[\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.