in addition to those commonly used ones to obtain the system-reservoir coupling Hamiltonian in Eq. (\[couple1\]). The neglecting of the cubic term $H_{nl}(t)$ has a physically different meaning from that of omitting the higher-order fluctuation terms in the previously used linearization around the steady states of the averaged Eq. (\[e2\]), because it is not necessary for our completely quantum treatment to require that the cavity mean field should overwhelm the corresponding cavity field fluctuation.

By using the associated Ito’s rules for the stochastic part [@supp1], one will obtain the following dynamical equation $$\begin{aligned}
\frac{d}{dt}\left( 
\begin{array}{c}
\hat{a}\\ 
\hat{a}^\dagger\\
\hat{b}\\
\hat{b}^\dagger
\end{array}
\right) 
&=&
\underbrace{\left( 
\begin{array}{cccc}
-\kappa & 0 & gE f(t)e^{-i\omega_m t} & gE f(t)e^{i\omega_m t}\\ 
0 &-\kappa & gE f^\ast(t)e^{-i\omega_m t}& gE f^\ast(t)e^{i\omega_m t}\\
-gE f^\ast (t)e^{i\omega_m t}&gE f(t)e^{i\omega_m t}&-\gamma_m & 0\\
gE f^\ast(t)e^{-i\omega_m t} &-gE f(t)e^{-i\omega_m t} &0 &-\gamma_m\\
\end{array}
\right)}_{\hat{M}(t)}\left( 
\begin{array}{c}
\hat{a}\\ 
\hat{a}^\dagger\\
\hat{b}\\
\hat{b}^\dagger
\end{array}
\right) \nonumber\\
 &+&
\underbrace{\left( \begin{array}{c}
i\kappa f(t)E \\ 
-i\kappa f^\ast(t)E\\
ig|f(t)E|^2e^{i\omega_mt}\\
-ig|f(t)E|^2e^{-i\omega_mt}\\
\end{array}
\right)}_{\vec{\lambda}(t)}+\underbrace{\left( \begin{array}{c}
\sqrt{2\kappa}\hat{\xi}_c(t) \\ 
\sqrt{2\kappa}\hat{\xi}^\dagger_c(t)\\
\sqrt{2\gamma_m}\hat{\xi}_m(t)\\
\sqrt{2\gamma_m}\hat{\xi}^\dagger_m(t)\\
\end{array}
\right)}_{\vec{\hat{\eta}}(t)}
\label{dynamics}\end{aligned}$$ with the effective Hamiltonian $H^{in}(t)+H_{nl}(t)$. This extended form of Eqs. (2)-(3) in the main text contains the factor $f(t)=\frac{e^{i\Delta t}-1}{\Delta}$ in the dynamical matrix $\hat{M}(t)$ and the coherent drive 