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!