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
where \(dW_n\) is a Wiener increment, which has the expectation values \(E[dW] = 0\) and \(E[dW^2] = dt\).
Stochastic Schrodinger Equation
The stochastic Schrodinger equation is given by (see section 4.4, [Wis09])
where \(H\) is the Hamiltonian, \(S_n\) are the stochastic collapse operators, and \(e_n\) is
In QuTiP, this equation can be solved using the function ssesolve
,
which is implemented by defining \(d_1\) and \(d_{2,n}\) from Equation (1) as
and
The solver ssesolve
will construct the operators
\(d_1\) and \(d_{2,n}\) once the user passes the Hamiltonian (H
) and
the stochastic operator list (sc_ops
). As with the 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 ssesolve
.
Stochastic Master Equation
When the initial state of the system is a density matrix \(\rho\), the stochastic master equation solver qutip.stochastic.smesolve
must be used.
The stochastic master equation is given by (see section 4.4, [Wis09])
where
and
In QuTiP, solutions for the stochastic master equation are obtained using the solver
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_1\) in Equation
(1), takes into account all operators \(C_i\) and \(S_n\):
The stochastic part, \(d_{2,n}\), is given solely by the operators \(S_n\)
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 \(\kappa\) and the
\(\Delta\) 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
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
Changing the Number of Trajectories, Reproducibility and Running trajectories in parallel for details.