Simulate dual Lindblad master equations in the Heisenberg picture in QuTiP

Question

In QuTiP, it is possible to solve Lindblad master equations describing the time evolution of an open quantum system $\rho$:

$$ \dot{\rho}(t)=-\frac{i}{\hbar}[H(t), \rho(t)]+\sum_n \frac{1}{2}\left[2 C_n \rho(t) C_n^{\dagger}-\rho(t) C_n^{\dagger} C_n-C_n^{\dagger} C_n \rho(t)\right] $$

where $H$ is the Hamiltonian of the system and $(C_n)_n$ are collapse operators. The simulation is done using the mesolve function.

Is it possible to simulate the dual of this evolution in the Heisenberg picture ? I would like to simulate the evolution of an operator $E$ via:

$$ \dot{E}(t)= +\frac{i}{\hbar}[H(t), E(t)]+\sum_n \frac{1}{2}\left[2 C_n^{\dagger} E(t) C_n-E(t) C_n^{\dagger} C_n-C_n^{\dagger} C_n E(t)\right] $$

Answer 1

The function mesolve (see here) takes in an array of collapse operators (c_ops) which, if given as superoperators, will be treated as direct contributions to the total Liouvillian. In this way, if you pass in $E(t=0)$ as your initial "density matrix", you should vectorize the Heisenberg contributions $C_n^\dagger (\cdot) C_n, \ \mathbb{I} (\cdot) C_n^\dagger C_n, \ C_n^\dagger C_n (\cdot) \mathbb{I}$ and use them as your c_ops, which should produce the correct evolution for $\dot{E}(t)$. You can also pass in an array of expectation values (e_ops) to evaluate; the initial state $\rho_0$ of the Schrodinger evolution can be used as the "observable" to measure in the Heisenberg picture, so you can compare the final expectation values between the two pictures to make sure everything agrees.