term $\vec{\hat{\lambda}}(t)$. The numerical solution to this differential equation can be found by calculating the time-ordered exponentials (in Eq. (5) of the main text) as the products of the matrices $(I+\hat{M}(t_i)\delta t)$ at the different $t_i$ [@s-laser], where $\delta t$ is a small iteration step. The corresponding quantum master equation for the reduced system state $\rho(t)$ reads $$\begin{aligned}
\dot{\rho}(t)&=&-i[H_{couple}(t)-H_{displace}(t),\rho(t)]-\kappa\big(\hat{o}_a^\dagger\hat{o}_a\rho(t)+\rho(t)\hat{o}_a^\dagger\hat{o}_a-2\hat{o}_a\rho(t)\hat{o}_a^\dagger\big)\nonumber\\
&-&\gamma_m(n_{th}+1)\big(\hat{b}^\dagger\hat{b}\rho(t)+\rho(t)\hat{b}^\dagger\hat{b}-2\hat{b}\rho(t)\hat{b}^\dagger\big)-\gamma_m n_{th}\big(\hat{b}\hat{b}^\dagger\rho(t)+\rho(t)\hat{b}\hat{b}^\dagger-2\hat{b}^\dagger\rho(t)\hat{b}\big)\end{aligned}$$ with the notations in Eq. (\[Hi\]), where $\hat{o}_a=\hat{a}-i\frac{e^{i\Delta t}-1}{\Delta}E$.

As a comparison, in the previous approach of linearization with the decomposition $\hat{a}=\alpha+\delta \hat{a}$ or the displacement $\hat{a}\rightarrow  \hat{a}+\alpha$ with a steady-state value $\alpha$ of the cavity mode $\hat{a}$, the linearized Hamiltonian excluding the coherent and noise drive terms takes the form [@supp2] $$\begin{aligned}
H_{linear}=G(\delta \hat{a} e^{-i\Delta t}+\delta \hat{a}^\dagger e^{i\Delta t})(\hat{b}e^{-i\omega_m t}+\hat{b}^\dagger e^{i\omega_m t}),
\label{lnh}\end{aligned}$$ with $|G|=g|\alpha|$ being a constant. A variation form of this linearized Hamiltonian has been used to study real-time cooling processes [@c6]. In reality, however, the constant optomechanical coupling intensity $G$, which equals to $g$ multiplied by the steady value $\alpha$ appearing after a certain period of time, can not exist before that time, for instance, shortly after turning on the cooling laser beam at $t=0$. The cavity field will certainly take time to evolve to 