\downarrow }  =  
    - i [ \hat H_0  ,\hat \rho_{ \uparrow  \downarrow } ]  - i \{ \hat H_{\mathrm{int}}(t) ,\hat \rho_{ \uparrow  \downarrow } \}+ \gamma ( {\bar n}_\mathrm{th} + 1 ) \mathcal{D} [ \hat c ]\hat \rho_{ \uparrow  \downarrow }  + \gamma {\bar n}_\mathrm{th} \mathcal{D} [\hat c^\dag ]\hat \rho_{ \uparrow  \downarrow },
\label{SEq:meqcoh}$$ which is a direct extension, with a time-modulation $F(t)$ in interaction $\hat H_{\mathrm{int}} \left( t\right)=  {\lambda} F(t) \hat n \hat \sigma_z /2$, of the technique used in Ref. [@Utami2007]. Here we use a constant coefficient $\lambda$ to keep track of orders in expansion on the coupling; by the end of the calculation, one can always set $\lambda=1$. We stress that if we replace the time-independent coupling $\lambda$ with a time-dependent one, the relevant derivations in Ref. [@Utami2007] still hold rigorously, and we refer interested readers to this paper for more detail.

Without loss of generality, the system initial state can be chosen as a product state between the qubit and the cavity, with the cavity in thermal equilibrium. Thus, Wigner function $W(x,p;t)$ of the coherence operator $\hat \rho _{ \uparrow  \downarrow } (t) $ is Gaussian throughout time evolution. Moreover, for the Fourier transform of $W(x,p;t)$, 