of Eq. (\[solution\]) with respect to the initial cavity vacuum state and mechanical thermal state, and $$\begin{aligned}
 n^{(n)}_m(t)&=&2\kappa\int_0^{t} d\tau |d_{41}(t,\tau)|^2+2\gamma_m \int_0^{t} d\tau |d_{44}(t,\tau)|^2n_{th}\nonumber\\
&+& 2\gamma_m \int_0^{t} d\tau |d_{43}(t,\tau)|^2(n_{th}+1)
\label{n}\end{aligned}$$ from averaging the following noise drive term over the reservoir states by means of Eq. (\[noise\]). A meaningful scenario beyond thermal equilibrium is cooling—the noise contribution in Eq. (\[n\]) finally stabilizes to a thermal phonon number $n_{m,f}$ less than the initial occupation $n_{th}$, while the contribution $n_m^{(s)}(t)$ from the evolved system operators gradually tends to zero with time.

Going back to the specific terms in Eqs. (2)-(3), one finds that an increased magnitude of the third term on their right sides (from a squeezing type coupling) can enhance the phonon numbers in both Eq. (\[s\]) and Eq. (\[n\]). The cooling action, on the other hand, manifests as the second term of a BS coupling in these equations. To let the BS action dominate, one could set the detuning $\Delta$ to be the mechanical frequency $\omega_m$, so that the factor $f(t)e^{-i\omega_mt}$ in Eq. (2) will become $\Gamma_b(t)/\omega_m=(1-e^{-i\omega_m t})/\omega_m$ containing a non-oscillating term. Meanwhile, the factor $f(t)e^{i\omega_mt}$ in the squeezing coupling term will be $\Gamma_s(t)/\omega_m=(e^{i2\omega_m t}-e^{i\omega_m t})/\omega_m$. A cooling can be performed with 