Equation of motion for various operators in open quantum systems I am trying to reproduce the calculations presented on page 4 in arXiv:1511.03347. The Hamiltonian (Eq. 2.4) is given by
$H= \hbar \omega (a^{\dagger} a + \frac{1}{2}) - B \sqrt{\frac{\hbar g}{2 \omega}} (a^{\dagger} + a) $,
where $a^{\dagger}(a)$ denotes creation~(annihilation) operator, $B$ and $ g$ are real constants and $\omega$ is time-dependent. The time-evolution of the density matrix ($\rho$) is governed by the Lindblad master equation given by
$\frac{d \rho}{dt}= {\cal L}\rho=-\frac{i }{\hbar} [H, \rho] + \gamma (n_{\omega}+1)
\left[a \rho a^{\dagger} -\frac{1}{2} (a^{\dagger}a \rho + \rho a^{\dagger} a)  \right]
+ \gamma n_{\omega}
\left[a^{\dagger} \rho a -\frac{1}{2} (aa^{\dagger} \rho + \rho a a^{\dagger})  \right]
$,
where $\gamma$ is the coupling constant and the Bose-Einstein distribution at temperature $T$ is given by $n_{\omega}=1/(e^{\hbar \omega /T}-1)$.
In Eqs.~(2.7-2.9), equation of motion for three expectation values, namely $\langle aa\rangle $, $\langle a^{\dagger}a\rangle $ and $\langle a\rangle $, are presented
$
\partial_{t}\langle aa\rangle =-2 \left[ \frac{\gamma}{2} + i \omega \right] \langle aa\rangle  + i \sqrt{\frac{2g}{\hbar \omega}} B \langle a\rangle 
,\\
\partial_{t}\langle a^{\dagger}a\rangle =-\gamma \langle a^{\dagger} a\rangle +\gamma n_{\omega} 
+i \sqrt{\frac{g}{2 \hbar \omega}} B (\langle a^{\dagger} \rangle  -\langle a\rangle  )
,\\
\partial_{t}\langle a\rangle =-\left[ \frac{\gamma}{2} + i \omega \right] \langle a \rangle
+ i \sqrt{\frac{g}{2 \hbar \omega}} B.
$
Terms proportional to $\omega$ and $B$ in the above equations, originated from $[H, \rho]$, are easy to get. However, I am more concerned about dissipative terms proportional to $\gamma$.
I have tried to either work with ${\cal L}^{\dagger}$ or ${\cal L}$ as has been discussed in this and this threads. But in both cases, I ended up with higher-order correlation functions such as $\langle a \rho aaa^{\dagger}\rangle $ $\langle a^{\dagger} \rho aaa^{\dagger}\rangle $. As both terms with $\gamma$ in the master equation possess positive signs, even after using the commutation relations of operators, I couldn't simplify my equations to reproduce the above equations of motions. As I have encountered these relations in other papers, I assume it is a standard procedure. I would appreciate it if you could point me to a reference in which these calculations are presented in more detail or if you could answer this thread by deriving one of the above EOMs explicitly.
 A: First of all, you have a typo in the equations of motion: they are missing the derivative in their correspondent left-hand side. Second, the solution of the associated adjoint Lindblad master equation of the one you posted can be solved via the so-called "third quantization" (see, e.g., arXiv:0801.1257) or via Ansatz for $a_H(t)$ (see below).
Let $a_H(t) \equiv \exp(iHt/\hbar)a \exp(-iHt/\hbar)$ be the Heisenberg operator associated with $a$. Here, $H$ is the free Hamiltonian you provided. For a Lindblad equation of the form $\partial_t\rho(t) = \mathcal{L}[\rho(t)]$, where
$$ \mathcal{L}[\bullet] = -\frac{i}{\hbar}[H,\bullet] + \sum_{k}\gamma_k\left( L_k\bullet L_k^\dagger - \frac{1}{2}\{L_k^\dagger L_k, \bullet \}\right) $$
is the Lindbladian, one can find the associated adjoint Lindblad equation for a general Heisenberg operator $A_H(t)$ via $\text{Tr}\{A\mathcal{L}[\rho(t)]\} = \text{Tr}\{(\mathcal{L}^\dagger[A_H(t)])\rho(0)\}$:
$$ \partial_t A_H(t) = \mathcal{L}^\dagger[A_H(t)] =  \frac{i}{\hbar}[H,A_H(t)] + \sum_{k}\gamma_k\left( L_k^\dagger A_H(t) L_k - \frac{1}{2}\{L_k^\dagger L_k, A_H(t) \}\right). $$
Notes: (i) see what changed between $\mathcal{L}$ and $\mathcal{L}^\dagger$. (ii) we will denote $\mathcal{D}^\dagger[L_k][\bullet] = \gamma_k\left( L_k^\dagger \bullet L_k - \frac{1}{2}\{L_k^\dagger L_k, \bullet \}\right)$ as the adjoint dissipator of $L_k$.
In our case, $A_H(t)$ will be $ a_H(t), a_H(t)a_H(t),$ and $a^\dagger_H(t)a_H(t)$. Further, assume that $a_H(t) = f(t)a$, where $f(t)$ is some complex-valued function with the condition that $f(0) = 1$, i.e., the Schrödinger picture coincides with the Heisenberg one at $t = 0$. Let us replace $a_H(t) = f(t)a$ in the above adjoint Lindblad equation where we will use $[a,a^\dagger] = 1$: Term-by-term

*

*$\frac{i}{\hbar}[H,a_H(t)] = \frac{f(t)i}{\hbar}[\hbar \omega a^\dagger a - B\sqrt{\frac{\hbar g}{2\omega}}(a^\dagger + a), a] = -i\omega a_H(t) + iB\sqrt{\frac{g}{2\omega \hbar}}f(t)$.

*$ \mathcal{D}^\dagger[a][a_H(t)] = f(t)\gamma n_\omega\left[aaa^\dagger - \frac{1}{2}\{aa^\dagger, a \}   \right] = f(t)\gamma n_\omega\left[a^2a^\dagger - \frac{1}{2}(aa^\dagger a + a^2 a^\dagger) \right] =  \frac{f(t)\gamma n_\omega}{2}\left[a^2a^\dagger - a(a a^\dagger - 1) \right] =  \frac{\gamma n_\omega}{2} a_H(t).  $

*$ \mathcal{D}^\dagger[a^\dagger][a_H(t)] = f(t)\gamma (n_\omega+1)\left[a^\dagger aa - \frac{1}{2}\{aa^\dagger, a \}   \right] = -\frac{\gamma (n_\omega+1)}{2} a_H(t)$.

Hence, $$ \partial_t a_H(t) =  -i\omega a_H(t) + iB\sqrt{\frac{g}{2\omega \hbar}}f(t) - \frac{\gamma}{2}a_H(t). $$
Going back to the Schrödinger picture and taking the expectation value of the resultant equation gives the third equation you provided
$$ \partial_t \left\langle a \right\rangle(t) = -\left( \frac{\gamma}{2} +  i\omega  \right)\left\langle a \right\rangle(t) + iB\sqrt{\frac{g}{2\omega \hbar}}. $$
The other two are found in the same way.
