dynamics-stochastic.rst

Stochastic Solver

When a quantum system is subjected to continuous measurement, through homodyne detection for example, it is possible to simulate the conditional quantum state using stochastic Schrodinger and master equations. The solution of these stochastic equations are quantum trajectories, which represent the conditioned evolution of the system given a specific measurement record.

In general, the stochastic evolution of a quantum state is calculated in QuTiP by solving the general equation

dρ(t) = d₁ρ dt + ∑_nd_2, nρ dW_n,

where dW_n is a Wiener increment, which has the expectation values E[dW] = 0 and E[dW²] = dt.

Stochastic Schrodinger Equation

The stochastic Schrodinger equation is given by (see section 4.4, [Wis09])

$$d \psi(t) = - i H \psi(t) dt - \sum_n \left( \frac{S_n^\dagger S_n}{2} -\frac{e_n}{2} S_n + \frac{e_n^2}{8} \right) \psi(t) dt + \sum_n \left( S_n - \frac{e_n}{2} \right) \psi(t) dW_n,$$

where H is the Hamiltonian, S_n are the stochastic collapse operators, and e_n is

$$e_n = \left<\psi(t)|S_n + S_n^\dagger|\psi(t)\right>$$

In QuTiP, this equation can be solved using the function ~qutip.solver.stochastic.ssesolve, which is implemented by defining d₁ and d_2, n from Equation general_form as

$$d_1 = -iH - \frac{1}{2} \sum_n \left(S_n^\dagger S_n - e_n S_n + \frac{e_i^2}{4} \right),$$

and

$$d_{2, n} = S_n - \frac{e_n}{2}.$$

The solver ~qutip.solver.stochastic.ssesolve will construct the operators d₁ and d_2, n once the user passes the Hamiltonian (H) and the stochastic operator list (sc_ops). As with the ~qutip.solver.mcsolve.mcsolve, the number of trajectories and the seed for the noise realisation can be fixed using the arguments: ntraj and seeds, respectively. If the user also requires the measurement output, the options entry {"store_measurement": True} should be included.

Per default, homodyne is used. Heterodyne detections can be easily simulated by passing the arguments 'heterodyne=True' to ~qutip.solver.stochastic.ssesolve.

Stochastic Master Equation

When the initial state of the system is a density matrix ρ, the stochastic master equation solver qutip.stochastic.smesolve must be used. The stochastic master equation is given by (see section 4.4, [Wis09])

dρ(t) =  − i[H, ρ(t)]dt + D[A]ρ(t)dt + ℋ[A]ρdW(t)

where

$$D[A] \rho = \frac{1}{2} \left[2 A \rho A^\dagger - \rho A^\dagger A - A^\dagger A \rho \right],$$

and

ℋ[A]ρ = Aρ(t) + ρ(t)A^† − tr[Aρ(t) + ρ(t)A^†].

In QuTiP, solutions for the stochastic master equation are obtained using the solver ~qutip.solver.stochastic.smesolve. The implementation takes into account 2 types of collapse operators. C_i (c_ops) represent the dissipation in the environment, while S_n (sc_ops) are monitored operators. The deterministic part of the evolution, described by the d₁ in Equation general_form, takes into account all operators C_i and S_n:

d₁ =  − i[H(t), ρ(t)] + ∑_iD[C_i]ρ + ∑_nD[S_n]ρ,

The stochastic part, d_2, n, is given solely by the operators S_n

d_2, n = S_nρ(t) + ρ(t)S_n^† − tr(S_nρ(t)+ρ(t)S_n^†) ρ(t).

As in the stochastic Schrodinger equation, heterodyne detection can be chosen by passing heterodyne=True.

Example

Below, we solve the dynamics for an optical cavity at 0K whose output is monitored using homodyne detection. The cavity decay rate is given by κ and the Δ is the cavity detuning with respect to the driving field. The measurement operators can be passed using the option m_ops. The homodyne current J_x is calculated using

J_x = ⟨x⟩ + dW/dt,

where x is the operator passed using m_ops. The results are available in result.measurements.

# parameters DIM = 20 # Hilbert space dimension DELTA = 5 * 2 * np.pi # cavity detuning KAPPA = 2 # cavity decay rate INTENSITY = 4 # intensity of initial state NUMBER_OF_TRAJECTORIES = 500

# operators a = destroy(DIM) x = a + a.dag() H = DELTA * a.dag() * a

rho_0 = coherent(DIM, np.sqrt(INTENSITY)) times = np.arange(0, 1, 0.0025)

stoc_solution = smesolve(

H, rho_0, times, c_ops=[], sc_ops=[np.sqrt(KAPPA) * a], e_ops=[x], ntraj=NUMBER_OF_TRAJECTORIES, options={"dt": 0.00125, "store_measurement": True,}

)

fig, ax = plt.subplots() ax.set_title('Stochastic Master Equation - Homodyne Detection') ax.plot(times[1:], np.array(stoc_solution.measurement).mean(axis=0)[0, :].real, 'r', lw=2, label=r'$J_x$') ax.plot(times, stoc_solution.expect[0], 'k', lw=2, label=r'$langle x rangle$') ax.set_xlabel('Time') ax.legend()

The stochastic solvers share many features with .mcsolve, such as end conditions, seed control and running in parallel. See the sections monte-ntraj, monte-seeds and monte-parallel for details.