# Solving linear Lindblad equations in the Heisenberg picture?

I'm interested in solving the explicit time dependence of operators in a simple open system described by a Lindblad equation. The concrete example I'm interested in is a harmonic oscillator with the usual thermalizing noise, described for example by the Lindblad equation for its density matrix (following the conventions in https://en.wikipedia.org/wiki/Lindbladian)

$$\dot{\rho} = -i [H,\rho] + \sum_{i=1,2} L_i^{\dagger} \rho L_i - \frac{1}{2} \left\{ L_i^{\dagger} L_i,\rho \right\}$$

with noise operators

$$L_1 = \sqrt{\bar{n} \gamma} a^{\dagger}, L_2 = \sqrt{(\bar{n}+1)\gamma} a$$

where $$\bar{n}$$ is the equilibrium mode occupation and $$\gamma$$ is a damping rate. For $$\rho_T = N e^{-\omega a^{\dagger} a/T}$$ the thermal state, it's easy to work out that $$\dot{\rho}_T = 0$$. What I'd like to do is compute some time-time correlation functions in this state, like

$$\langle x(t) x(0) \rangle$$.

I'm wondering how one would go about computing this in the Heisenberg picture. My naive approach to this would be to work with the adjoint Lindblad equation (or whatever it's really called, the thing labeled "Heisenberg picture" in the wikipedia link above):

$$\dot{O} = i [H,O] + \sum_{i} L_i O L_i^{\dagger} - \frac{1}{2} \left\{ L_i^{\dagger} L_i, O \right\}$$

for some operator $$O$$, as discussed in for example here: Lindblad equation for heisenberg operators?

What's got me stuck is: the resulting equations for $$\dot{a}, \dot{a}^{\dagger}$$ are some fairly gross-looking coupled cubic system. Concretely, these equations are

$$\dot{a} = -i \omega a + \frac{\gamma \bar{n}}{2} \left[ a,aa^{\dagger} \right] + \frac{\gamma (\bar{n}+1)}{2} \left[ a^{\dagger} a, a \right]$$

$$\dot{a}^{\dagger} = i \omega a^{\dagger} + \frac{\gamma \bar{n}}{2} \left[ a a^{\dagger},a^{\dagger} \right] + \frac{\gamma (\bar{n}+1)}{2} \left[ a^{\dagger}, a^{\dagger} a \right]$$

Of course, it would be lovely to just reduce this to a linear system using the canonical commutation relations, but those are not preserved since the evolution is not unitary :)

I imagine there's a smart way to get the solution, but it's eluding me. Normally I think of equations like this by appending an input noise field and using the input-output formalism in quantum optics, but I'm now curious how one would go about working directly with the Lindbladian form of this system. Any help appreciated!

The key point is the fact that the evolution is Markovian, driven by the linear Liouvillian superoperator $$\mathcal{L}$$. In particular, we have [1]: $$\rho(t)=\phi(t)[\rho(0)],$$ where $$\phi(t)=\exp{\mathcal{L}t}$$ is the Markovian quantum map and the Liouvillian $$\mathcal{L}$$ is the generator of the quantum dynamical semigroup (in the scenario with time-independent semigroup, as in your case). In Schrödinger picture, we have: $$\frac{d\rho(t)}{dt}=\mathcal{L}[\rho(t)],$$ which is the master equation that you have written in the first line. To switch to the Heisenberg picture, we need to introduce an adjoint quantum map $$\phi^\dagger(t)$$ and, correspondingly, the adjoint Liouvillian superoperator $$\mathcal{L}^\dagger$$ (see [1] for an extensive derivation). Accordingly, the evolution of any operator in Heisenberg picture is given by: $$\frac{dO(t)}{dt}=\mathcal{L}^\dagger[O(t)],\qquad O(t)=\phi^\dagger(t)=\exp{\mathcal{L}^\dagger t}[O(0)].$$ Therefore, we can find the solution of the dynamics by applying $$\mathcal{L}^\dagger$$ on $$O(0)$$ only, exactly as in the Schrödinger picture with $$\rho(0)$$. We will then diagonalize the adjoint Liouvillian and find the exponential dynamics at any time. In the case you are proposing everything is already diagonal, given that we get: $$\frac{da(0)}{dt}=\mathcal{L}^\dagger[a(0)]=(-i\omega-\gamma/2)a(0),$$ yielding $$a(t)=\exp(-i\omega t-\gamma t/2)a(0)$$, which is the correct solution. To calculate multi-time correlations functions, you should first compute the mean values of each operator as explained above (i.e. diagonalize the Liouvillian), and then rely on the quantum regression theorem [1].
• Thanks! I'm confused by the physical meaning of this answer, though. If we have $a(t) = e^{-i \omega t-\gamma t/2} a(0)$, then in the thermal state, I would compute in Heisenberg picture that the position variance $\langle x(t)^2 \rangle \sim e^{-\gamma t}$ will go to zero at late times (I know I was asking about $\langle x(t) x(0) \rangle$, so this is a slightly separate issue). But the thermal state is the stationary solution for these dynamics, so shouldn't we just have $\langle x(t)^2 \rangle = x^2_{thermal}$ the usual thermal uncertainty, a constant? Dec 14, 2020 at 17:58
• This is a tricky aspect of the GKLS Heisenberg picture. The point is that the Heisenberg picture is only providing us with the evolved operator to compute its precise mean value, i.e. if you have $a(t)$ then you can only calculate $\langle a(t)\rangle$. If you compose two Heisenberg operators $A(t)$ and $B(t)$, in general $\langle AB(t)\rangle\neq\langle A(t)B(t)\rangle$ (note that $A$ and $B$ are system operators only, different from $A\otimes\mathbb{I}_E$). This is due to the fact that, in general, $\mathcal{L}^\dagger[AB]\neq\mathcal{L}^\dagger[A]\mathcal{L}^\dagger[B]$. Dec 15, 2020 at 9:08
• Therefore, if you want to compute $\langle x(t)^2\rangle$, you will need to derive the Heisenberg expression for $x^2(t)$, following the same lines of my answer above. You will then recover the right expression for the stationary thermal state. The same applies for the computation of the mean value of the number of photons in the system, $\langle a^\dagger a(t)\rangle\neq\langle a^\dagger(t) a(t)\rangle$. Dec 15, 2020 at 9:12
• Ah, ok! Thanks much. So if I understand correctly, this is a consequence of the map being CPTP but not unitary. In a unitary setting I'd have $A(t) B(t) = U^{\dagger} A U U^{\dagger} B U = U^{\dagger} AB U$, but no equivalent thing can be done here, since $\phi^{\dagger} \neq \phi^{-1}$. Is that the right intuition? Dec 16, 2020 at 17:04