this initially prepared state with no external drive ($E=0$), the mechanical mode takes the exact form $\hat{b}(t)=e^{-\gamma_m t}\hat{b}+\sqrt{2\gamma_m}\int_0^t d\tau e^{-\gamma_m(t-\tau)}\hat{\xi}_{m}(\tau)$ from Eq. (3). It is the second noise drive term that maintains the invariant phonon number $\langle \hat{b}^\dagger\hat{b}\rangle=n_{th}$ under the thermal equilibrium, while the contribution from the first term lowers with time. The noise actions are thus essential to a quantum OMS.

In a general situation it is clearer to use the equation $$\begin{aligned}
\frac{d}{dt}\vec{\hat{c}}(t)=\hat{M}(t)\vec{\hat{c}}(t)+\vec{\lambda}(t)+\vec{\hat{\eta}}(t)
\label{d-equation}\end{aligned}$$ about the complete set $\vec{\hat{c}}(t)=(\hat{a}(t),\hat{a}^\dagger(t),\hat{b}(t),\hat{b}^\dagger(t))^T$ of the system modes, where $\vec{\lambda}=(\lambda_c,\lambda_c^\ast, \lambda_m, \lambda_m^\ast)^T$ with $\lambda_c(t)=i\kappa f(t)E$ and $\lambda_m(t)=ig|f(t)E|^2e^{i\omega_mt}$, $\vec{\hat{\eta}}=(\hat{\eta}_c, \hat{\eta}_c^\dagger, \hat{\eta}_m,\hat{\eta}_m^\dagger)^T$ with $\hat{\eta}_c(t)=\sqrt{2\kappa}\hat{\xi}_c(t)$ and $\hat{\eta}_m(t)=\sqrt{2\gamma_m}\hat{\xi}_m(t)$, and the detailed matrix $\hat{M}(t)$ can be found in [@supp]. The general solution to Eq. (\[d-equation\]) reads $$\begin{aligned}
\vec{\hat{c}}(t)&=&{\cal T}e^{\int_0^t d\tau \hat{M}(\tau)}\vec{\hat{c}}\nonumber\\
&+&\int_0^t d\tau~ {\cal T}e^{\int_\tau^t dt' \hat{M}(t')}\big(\vec{\lambda}(\tau)+\vec{\hat{\eta}}(\tau)\big),
\label{solution}\end{aligned}$$ with $\vec{\hat{c}}=(\hat{a},\hat{a}^\dagger,\hat{b},\hat{b}^\dagger)^T$. The time-ordered exponentials appear because we have $[\hat{M}(t),\hat{M}(t')]\neq 0$ for $t\neq t'$. In terms of the notation $d_{ij}(t,\tau)$ for the matrix elements $[{\cal T}e^{\int_\tau^t dt' \hat{M}(t')}]_{ij}$, the evolving thermal phonon number $n_m(t)=\langle \hat{b}^\dagger\hat{b}(t)\rangle-|\langle \hat{b}(t)\rangle|^2$ (the subtraction of the coherent drive $\vec{\lambda}(t)$’s contribution is to be discussed below) consists of two parts, $$\begin{aligned}
n^{(s)}_m(t)&=&|d_{41}(t,0)|^2+|d_{44}(t,0)|^2n_{th}\nonumber\\
&+&|d_{43}(t,0)|^2(n_{th}+1)
\label{s}\end{aligned}$$ from taking the expectation value of the first term on the right side 