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.
Best Answer
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
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.