\betagin{document}

\title[Visualizing the quantum interaction picture in phase space]{Visualizing the quantum interaction picture in phase space}
\author{Bahar Mehmani$^1$, Andrea Aiello$^{1,2}$}
\address{$^1$ Max Planck Institute for the Science of Light, G\"{u}nter-Scharowsky-Stra$\ss$e 1/Bau 24, D-91058 Erlangen, Germany}
\address{$^2$ Institute for Optics, Information and Photonics, University Erlangen-N\"{u}rnberg, Staudtstra$\ss$e 7/B2, D-91058 Erlangen, Germany}
\ead{bahar.mehmani@mpg.mpl.de} 

\betagin{abstract}
We illustrate the correspondence between the quantum \textnormal{Interaction}\ Picture-evolution of the state of a quantum system in Hilbert space and a combination of local and global transformations of its Wigner function in phase space. To this aim, we consider the time-evolution of a quantized harmonic oscillator driven by both a linear and a quadratic (in terms of bosonic creation and annihilation operators) potentials and employ the Magnus series to derive the exact form of the time-evolution operator. In this case, the Interaction Picture corresponds to a local transformation of phase space-reference frame into the one that is co-moving with the Wigner function. 
\end{abstract}
\hat at pacs{03.65.-w, 03.65.Ca}
\submitto{\NJP}
\maketitle
 \section{Introduction}\lambdabel{intro}
One of the most important lessons that one learns in quantum mechanics courses is to choose a proper picture  for describing the dynamics of a system. An appropriate picture is the one in which the physical properties of the system can be easily revealed and the mathematical calculations involved are relatively simple. In general, there are three pictures for this aim: the Schr\"{o}dinger, the \textnormal{Heisenberg}, and the \textnormal{Interaction}\ picture~\cite{Parisi}. The Schr\"{o}dinger\  picture (SP) is more suitable for studying closed and conservative systems. In this picture the only time-dependent quantity of the system is its quantum state. Correspondingly, observables of the system are constant in time. However, the number of situations in which a quantum system can be considered as a closed one is limited and, furthermore, most of the interesting phenomenon in the quantum world such as decoherence and optical cooling, to just name two, happen when a quantum system couples to its surrounding environment, to an external force or to both of them. Hence, the Hamiltonian describing the system is more complicated due to the presence of the environment and/or the external force. In many of such situations it is impossible to find an exact solution to the Schr\"{o}dinger\ equation. However, sometimes, adopting the \textnormal{Heisenberg}\ picture (HP) makes it possible to find the time-evolution of the expectation values of the observables of the system. In HP the state does not vary with time while the observables are time-dependent and their evolution is described by the \textnormal{Heisenberg}\ equation of motion. 
Finally, the last scheme, namely the \textnormal{Interaction}\ picture (IP) one, is the most suitable picture when the Hamiltonian of a system can be written as a sum of two parts: a time-independent term of which the eigenstates and eigenenergies are known, and a usually time-dependent term which influences the dynamics of the system. Employing this picture enables one to set the dynamics arising from the time-independent Hamiltonian aside and focus on the influence of its time-dependent part on the evolution of the system. That is why this picture is frequently employed in quantum optics where  matter interacts with the radiation field. 

All these three pictures are described in terms of unitary transformations in Hilbert space. The aim of this work is to map such transformations in phase space and illustrate the simplicity of the IP time-evolution of a quantum system in comparison to its SP counterpart. Although SP and HP have been mapped into classical active and passive transformations~\cite{Parisi}, to the best of our knowledge, there is no such direct comparison for the interaction picture in the physics literature.
Here we fill this gap and show that by taking active local and global transformations one can transform the IP to the SP or HP. In this way, it is easier to understand the concept of performing different unitary transformations in Hilbert space and its consequences on the evolution of the state of the system of interest.

To this end, we take a single harmonic oscillator as our quantum system, initially prepared in an ideal squeezed state~\cite{Mandel} and first study its free time-evolution in SP and HP. Then by employing the Wigner-Weyl description we calculate the time-evolution of its corresponding Wigner function. The basic advantage of the Wigner-Weyl representation~\cite{Case08,Hillery,Agarwal,Schleich} is that the operators $\hat at x$ and $\hat at p$ turn into $c$-numbers  making the correspondence with canonical transformations possible. The Wigner function of the coherent and squeezed states of a quantized harmonic oscillator is a Gaussian function of these $c-$number variables and its corresponding phase space distribution can be pictured by circles and ellipses, respectively~\cite{Kim90,Kim88,Kim89}. This is the subject of sections~\ref{Sec1:HO} and~\ref{Sec:Passive}.
In section~\ref{Sec:Interaction} we introduce two types of a time-dependent potential to the system: a linear and a quadratic term (in terms of the bosonic creation and annihilation operators). Then we take the IP to study the dynamics of the initially squeezed state under the driving potentials.  We explicitly show that the IP- time-evolution of the state in Hilbert space corresponds to a transformation to a local frame that is co-rotating with the Wigner function in phase space. Furthermore, we show switching from one picture to another is equivalent to performing different active and/or passive transformations in phase space.

\section{Basics of the quantized harmonic oscillator}\lambdabel{Sec1:HO}

A single mode harmonic oscillator with unit mass and frequency $\omega_0$ is described by the sum of its kinetic and potential energy
\betagin{eqnarray}\lambdabel{Hconv}
\hat H_0 = \frac{ \omega_0}{2} (\hat at p^2 + \hat at x^2).
\end{eqnarray}
Here $\hat at x$ and $\hat at p$ represent quadratures of the quantized harmonic oscillator and are defined in terms of the annihilation operator $\hat a$ and  the creation operator $\hat ad$ as
\numparts
\betagin{eqnarray}\lambdabel{quad-def}
\hat at x & = \frac{1}{\sqrt{2}rt{2}} ( \hat at a + \hat at a^\dagger ),\\
\hat at p & = \frac{1}{\rmi \sqrt{2}rt{2}}  ( \hat at a - \hat at a^\dagger ),
\end{eqnarray}
\endnumparts
where $\hbar$ is set to one. 
Alternatively, this system may be described in terms of $\hat a$ and $\hat ad$ as
\betagin{eqnarray}\lambdabel{HDirac}
\hat at H_0 =\omega_0 \left(\hat at a^\dagger \hat at a + \frac{1}{2}\right),
\end{eqnarray}
where the commutator of the creation and the annihilation operator is $[\hat a, \hat ad]= 1$. 
$\hat H_0$ has the energy eigenvalues $E_n$ defined as
\betagin{eqnarray}\lambdabel{En}
E_n =\omega_0 \left(n+ \frac{1}{2}\right), \thetaxtnormal{with} \qquad n = 0, 1,2, \cdots.
\end{eqnarray}
The energy eigenstates (number states) $\vert n \rangle$ such that:
\betagin{eqnarray}\lambdabel{n}
\hat H_0 \vert n \rangle = E_n \vert n \rangle.
\end{eqnarray}
Given $\hat H_0$, the time-evolution of quadrature operators $\hat at x$ and $\hat at p$ is given by the Heisenberg equation of motion as
\numparts
\betagin{eqnarray}\lambdabel{x,p:t}
\hat at x(t) & =\hat at x(0) \cos (\omega_0 t)  + \hat at p(0) \sin (\omega_0 t) ,\\
\hat at p(t) & = - \hat at x(0)  \sin (\omega_0 t) + \hat at p(0) \cos (\omega_0 t) .\lambdabel{p,x:t}
\end{eqnarray}
\endnumparts

Any superposition of the number states is a solution to the Schr\"{o}dinger equation for a quantized harmonic oscillator described by the Hamiltonian $\hat H_0$. There are two classes of such superpositions that are of great interest: coherent and squeezed states.
A coherent state $\vert \alpha \rangle$ is described as the displaced vacuum state $\vert 0 \rangle$ and is generated by a canonical transformation of the vacuum state $\left \vert 0 \right \rangle$ in phase space
\betagin{eqnarray}\lambdabel{cohstatevac}
\vert \alpha \rangle = \hat D(\alpha) \vert 0 \rangle,
\end{eqnarray}
where $\hat D(\alpha)$ is the unitary displacement operator defined as~\cite{Glauber,Scully}
\betagin{eqnarray}
\lambdabel{D_def}
\hat D(\alpha)= e^{( \alpha \hat at a^{\dagger} - \alpha^* \hat a)}.
\end{eqnarray}
The displacement parameter $\alpha$, is a complex number and may be decomposed to
\betagin{eqnarray}\lambdabel{alpha-def}
\alpha \equiv \frac{1}{\sqrt{2}} ( a + \rmi b),
\end{eqnarray}

Squeezed states are of interest because they provide reduced fluctuations in one quadrature, when compared with the coherent state. A squeezed state is described by the unitary squeezing operator~\cite{Mandel,Yuen76,Caves81,Walls83,Tapia93} $\hat S(\zeta)$ in a similar fashion. However, the difference lies in the fact that $\hat S(\zeta)$ is not linear in $\hat a$ and $\hat ad$ and is defined as 
\betagin{eqnarray}\lambdabel{SqueeOp-def}
\hat S(\zeta) = \exp \left[ \left(  \zeta^* \hat a^2- \zeta  {\hat a}^{\dagger 2}\right)/2 \right],\qquad \zeta = s e^{\rmi \theta},
\end{eqnarray}
where $s = \vert \zeta \vert$ is the strength of the squeezing and $\theta = \arg(\zeta)$ determines the direction along which the squeezing is performed. In phase space the squeezed states correspond to Gaussian distributions with unequal widths as opposed to the symmetric distribution of the coherent state.~\cite{Kim88}

Throughout this paper, for the sake of simplicity and without the loss of generality, we take both the displacement parameter $\alpha$ and the squeezing parameter $\zeta$ as real numbers.  Moreover, we assume the squeezed state to be an ideal squeezed state~\cite{Mandel}. 

Since we want to find a classical correspondence for the unitary transformations in the Interaction picture, we need to map the time-evolution of the ideal squeezed state mentioned above into classical phase space and study its time-evolution. This is best described with the help of the Wigner function and Weyl transform~\cite{Case08,Hillery,Agarwal}.
Knowing $\hat{\rho}(t)$ the state of the quantum system at any time, we can build up the symmetric characteristic function~\cite{Barnett} $\Tr[\hat{\rho}(t) \hat D(\chi)]$ and from there the corresponding SP-Wigner function $W(x,p;t)$ as
\betagin{eqnarray}\lambdabel{Wignerdef}
W(x,p;t) = \frac{1}{(2 \hat at pi)^2} \int \rmd\hat at xi \int \rmd\eta\,\rme^{-\rmi (\hat at xi x + \eta p)}\,\Tr[\hat{\rho}(t) \hat D(\chi)],
\end{eqnarray}
where $\hat D(\chi)$ is the displacement operator with the parameter $\chi$ defined as
\betagin{eqnarray}\lambdabel{chi-def}
\chi \equiv \frac{1}{\sqrt{2}} (- \eta +\rmi \hat at xi).
\end{eqnarray}
Both integrations are from $-\infty$ to $\infty$. 

\section{Schr\"{o}dinger and Heisenberg pictures vs. passive and active transformations}\lambdabel{Sec:Passive}

Before embarking on the classical correspondence of the IP, let us first explicitly show for a single-mode harmonic oscillator how SP and HP map to active and passive transformations in phase space, respectively.

For a harmonic oscillator with the Hamiltonian given by (\ref{HDirac}), the unitary time-evolution operator $\hat U_0$, which governs the evolution of the state, reads
\betagin{eqnarray}\lambdabel{U0}
\hat U_0(t) = \rme^{- \rmi \hat H_0 t} = \rme^{- \rmi \omega_0 t (\hat ad \hat a + \hat alf)}.
\end{eqnarray}
If the system is initially in an ideal squeezed state, 
\betagin{eqnarray}\lambdabel{int_ro_ini}
\hat{\rho}(0) = \hat D(\m) \hat S \left(s\right) \left \vert 0 \right \rangle \lambdac \hat S \left(s\right)m \hat D(- \m),
\end{eqnarray}
with $\mu_x$ representing the displacement along the $x$-axis in phase space and $s$ denoting the real squeezing parameter, then the time-evolution of the state is given by
\betagin{equation}
\lambdabel{rhot_sq-free}
\hat{\rho}(t) = \hat at U_0(t) \hat{\rho}(0) \hat at U^\dagger_0(t).
\end{equation}
\betagin{figure}
\subfloat[$\omega_0t=0$]{\lambdabel{fig:12}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig12.eps}}
\subfloat[$\omega_0t=\frac{\hat at pi}{4}$]{\lambdabel{fig:13}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig13}}
\subfloat[$\omega_0t=\frac{\hat at pi}{2}$]{\lambdabel{fig:14}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig14}}\\
\subfloat[$\omega_0t=\frac{3 \hat at pi}{4}$]{\lambdabel{fig:15}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig15}}
\subfloat[$\omega_0t=\hat at pi$]{\lambdabel{fig:16}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig16}}
\subfloat[$\omega_0t=\frac{5\hat at pi}{4}$]{\lambdabel{fig:17}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig17}}\\
\subfloat[$\omega_0t=\frac{3 \hat at pi}{2}$]{\lambdabel{fig:18}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig18.eps}}
\subfloat[$\omega_0t=\frac{7\hat at pi}{4}$]{\lambdabel{fig:19}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig19.eps}}
\subfloat[$\omega_0t=2 \hat at pi$]{\lambdabel{fig:20}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig20}}
\caption{Free time-evolution of the $\frac{1}{\rme}$-contour of the SP-Wigner function of an ideal squeezed state given by (\ref{Wshro_final}) in units of $\omega_0 t$ ($\sigma_x = 1, \sigma_p = \hat alf,\mu_x = 4$, in arbitrary units). The evolution is a combination of a local counterclockwise rotation of the ellipse about its centroid through angle $\omega_0 t$, and the global clockwise rotation of its center about the origin of the phase space on the dashed-line circle through the same angle.}
\lambdabel{fig.1}
\end{figure}
Knowing $\hat{\rho}(t)$, the Wigner function at time $t$ can be derived by first calculating the symmetric characteristic function as the expectation value of the displacement operator $\hat D(\chi)$ as
\betagin{eqnarray}\lambdabel{trfin_sq-free}
\fl \Tr\left[\hat{\rho}(t) \hat D(\chi)\right]&= \exp{\left[\rmi \hat at xi  \mu_x \cos (\omega_0 t) - \rmi \eta \mu_x \sin (\omega_0 t) \right]}
\exp{\left[- \frac{\hat at xi^2}{4}  \left( \frac{\cos^2(\omega_0 t)}{2 \sigma^2_p} + \frac{\sin^2 (\omega_0 t)}{2 \sigma^2_x}\right)\right]}\nonumber \\
 & \times\exp{\left[- \frac{\eta^2}{4}  \left(\frac{\cos^2(\omega_0 t)}{2 \sigma^2_x} + \frac{\sin^2 (\omega_0 t)}{2 \sigma^2_p} \right)\right]} \exp{\left[ - \frac{\hat at xi \eta}{8} \sin (2 \omega_0 t) \left( \frac{1}{\sigma^2_x}- \frac{1}{\sigma^2_p}  \right)\right] }.\nonumber \\
\end{eqnarray}
where we substitute $e^{- 2 s} = 2\sigma^2_x$ and $e^{ 2 s} =2 \sigma^2_p$ in the coefficients of $\hat at xi^2$,  $\eta^2$ and $\hat at xi \eta$. 
Substituting (\ref{trfin_sq-free}) into the definition of the Wigner function given by (\ref{Wignerdef}) and performing two Gaussian integrals over $\hat at xi$ and $\eta$ then yields to the SP-Wigner function as
\betagin{eqnarray}
\lambdabel{Wshro_final}
\fl W^{\textnormal{SP}}_0(x,p;t) = \frac{1}{\hat at pi} \exp{\left\{ - \hat alf \left[\bi r - R(-\omega_0\,t) \bmu \right]^T \cdot \Gamma^{- 1}(t) \cdot \left[\bi r - R(-\omega_0\,t) \bmu \right]\right\}},
\end{eqnarray}
where the superscripts ``$T$'' and ``${-1}$'' denote matrix transposition and inversion operations, respectively. The suffix $\thetaxtnormal{SP}$ denotes the Schr\"{o}dinger picture and the index $0$ in $W_0(x,p;t)$ represents the fact that the system evolves freely in time.
The vectors $\bmu$ and $\bi r$ represents the initial displacement of the Wigner function and the coordinates of phase space, respectively, as
\betagin{eqnarray}
\bmu  = \left(\betagin{array}{c}\mu_{x} \\0\end{array}\right),\lambdabel{mu_def} \qquad
\bi r= \left(\betagin{array}{c}x \\hat at p\end{array}\right).\lambdabel{vec_r_def}
\end{eqnarray}
The matrix $\Gamma^{- 1}(t)$ is defined as the clockwise rotation of the squeezing matrix through the angle $\omega_0 t$
\betagin{eqnarray}\lambdabel{Gamma_def}
\Gamma^{-1}(t) \equiv R(- \omega_0 t)\cdot \Xi^{-1} \cdot R( \omega_0 t),
\end{eqnarray}
where the diagonal $2\times2$ squeezing matrix $\Xi^{-1}$ is given by 
\betagin{eqnarray}\lambdabel{covmatrix_int}
\Xi^{-1}\equiv \left(\betagin{array}{cc}\frac{1}{\sigma^2_x} & 0 \\0 &\frac{1}{\sigma^2_p}\end{array}\right).
\end{eqnarray}
Notice that the matrix $R(\omega_0 t)$ defines a counterclockwise rotation by $\omega_0 t$ as
\betagin{eqnarray}\lambdabel{Rotation_def}
R(\omega_0 t) = \left(\betagin{array}{cc}\cos\omega_0 t &  - \sin\omega_0 t \\ \sin\omega_0 t & \cos\omega_0 t\end{array}\right).
\end{eqnarray} 

The free-time evolution of the $\frac{1}{\rme}$-contour of the SP-Wigner function of an ideal squeezed state is represented in figure~\ref{fig.1}. As it is illustrated in figures~\ref{fig:12}-\ref{fig:20}, the elliptical distribution performs two types of rotations: ${\it i)}$ a global clockwise rotation (active transformation) of the distribution about the origin of the phase space on a circular path defined by the initial amount of displacement; ${\it ii)}$ a local counterclockwise rotation  of the the squeezing direction through the same angle $\omega_0 t$ about its centroid. From the mathematical point of view, the global rotation is produced by the clockwise rotation matrix $R(-\omega_0t)$ acting on the vector $\bmu$, whereas the action of the  rotation matrix $R(\omega_0 t)$ on the squeezing matrix $\Xi^{-1}$ amounts to its local counterclockwise rotations.

In (\ref{Wshro_final}), if we take $R(- \omega_0 t)$ from left in to the first bracket and $R(\omega_0 t)$ from right in to the second bracket, the resulting Wigner function can be written in terms of a new phase space coordinate which is described by
\betagin{equation}\lambdabel{r-r'}
\bi{r'} = R(\omega_0 t)\, \bi{r}.
\end{equation}
In this new local frame, attached to the centroid of the Wigner function, the local counterclockwise rotation of the squeezing direction of the Wigner function about its centroid cancels out its global clockwise rotations about the origin
\betagin{eqnarray}\lambdabel{W':HP}
W^{\textnormal{HP}}_0(x',p') = \frac{1}{\hat at pi} \exp{\left[- \hat alf \left(\bi{r'} -  \bmu \right)^T\, \cdot \Xi^{-1}\cdot \, \left(\bi{r} - \bmu \right)\right]}.
\end{eqnarray}
The passive transformation (\ref{r-r'}) corresponds to changing the reference frame to the one in which $x(t)$ and $p(t)$ are rotated back to their initial values $x(0)$ and $p(0)$. The suffix \mbox{HP} represents that the Wigner function is described in the Heisenberg picture, in which the Wigner function is  fixed in time while the operators $\hat at x$ and $\hat at p$ evolve according to (\ref{x,p:t}) and (\ref{p,x:t}). 
Hence we see that the passive transformation (\ref{r-r'}) in phase space corresponds to the Heisenberg picture.

The main message of this section is that the action of the unitary transformation $\hat U_0$ on an ideal squeezed state in the Hilbert space corresponds to two types of rotations of its Wigner function in the phase space: ${\it i)}$ a global clockwise rotation about the origin through angle $\omega_0 t$, ${\it ii)}$ a local counterclockwise rotation of the Wigner function through angle $\omega_0 t$ about its centroid. Further, the action of $\hat U_0$ on the system in Hilbert space is equivalent to a counterclockwise rotation of the coordinate system with angle $\omega_0 t$, or equivalently, performing a passive transformation in phase space.

\section{Interaction Picture and mixed rotations}\lambdabel{Sec:Interaction}

The interaction picture is the most frequently employed picture in quantum dynamics. Specifically, it is the preferred picture in quantum optics where a system is driven by the time-dependent electric field of light. 

Considering the single mode harmonic oscillator described in section~\ref{Sec1:HO}, we now assume it undergoes a time-dependent potential $\hat V(t)$. Thus the total Hamiltonian is given by 
\betagin{equation}\lambdabel{interaction}
\hat H = \hat H_0 + \hat V(t),
\end{equation}
with $\hat H_0$ being defined in (\ref{HDirac}). $\hat V(t)$ influences the temporal behavior of the system and is a linear driving potential described as
\betagin{equation}\lambdabel{H1_def}
\hat V(t) =  g (\rme^{- \rmi \omega_1 t} \alpha \hat a + \rme^{\rmi \omega_1 t} \alpha^* \hat ad),
\end{equation}
where $g$ is a real coupling constant. Physically, this corresponds to a forced harmonic oscillator with the driving frequency $\omega_1$.

In the IP both the state and observables are time-dependent. The observables evolve in time according to the Heisenberg equation of motion as
\betagin{equation}
\lambdabel{Sch-Int-op}
\hat at {\Or}_{\thetaxtnormal I}(t) = \hat at U^\dagger_0(t) \, \hat at{\Or} \, \hat at U_0(t),
\end{equation}
where $\hat at U_0(t)$ is given by (\ref{U0}). 
Accordingly, the time-dependent potential in the interaction picture, $\hat V_{\thetaxtnormal I}$, reads as
\betagin{equation}\lambdabel{VI}
\hat V_{\thetaxtnormal I}(t) = g (e^{- \rmi \Omega t} \alpha \hat a + e^{\rmi \Omega t} \alpha^* \hat ad),
\end{equation}
where $\Omegaega$ is defined as the sum of the natural frequency $\omega_0$ and the driving frequency $\omega_1$:
\betagin{equation}\lambdabel{Omega}
\Omega \equiv \omega_0 + \omega_1.
\end{equation}
When $\Omega = 0\ (\omega_1 = - \omega_0)$, the interaction picture potential becomes time-independent. We shall study this case in more details in section~\ref{sec:Omzero}. For the moment we assume $\Omega \neq 0$.

The time evolution of the state is described by the interaction picture unitary transformation $\hat U_{\textnormal I}(t)$ as
\betagin{equation}\lambdabel{int-w-sol}
\hat{\rho}_{\thetaxtnormal I}(t) = \hat U_{\textnormal I}(t) \hat{\rho}(0) \hat at U^\dagger_{\thetaxtnormal I}(t),
\end{equation}
where $\hat U_{\textnormal I}(t)$ is given by the Magnus series~\cite{Magnus}
\betagin{equation}\lambdabel{UMagnus}
\hat U_{\textnormal I}(t) = \exp{\left\{ - \rmi \left[\hat at A_1(t) + \hat at A_2(t) + \hat at A_3(t) + \cdots \right]\right\}}.
\end{equation}
The first three terms in the exponential are given by
\numparts
\betagin{eqnarray}\lambdabel{Magnus}
\fl\hat at A_1(t) =  \int_0^t \rmd t_1\, \hat at V_{\thetaxtnormal I}(t_1),\lambdabel{A1Mag}\\
\fl\hat at A_2(t) = \frac{1}{2! \rmi } \int_0^t \rmd t_2 \int_0^{t_2} \rmd t_1\, \left[ \hat at V_{\thetaxtnormal I}(t_1), \hat at V_{\thetaxtnormal I}(t_2)\right],\\\lambdabel{A2Mag}
\fl\hat at A_3(t) = \frac{1}{3! (\rmi^2 )} \int_0^t \rmd t_3 \int_0^{t_3} \rmd t_2 \int_0^{t_2} \left\{ \left[\hat at V_{\thetaxtnormal I}(t_1),\left[ \hat at V_{\thetaxtnormal I}(t_2), \hat at V_{\thetaxtnormal I}(t_3\right] \right] + \left[\hat at V_{\thetaxtnormal I}(t_3), \left[ \hat at V_{\thetaxtnormal I}(t_2), \hat at V_{\thetaxtnormal I}(t_1)\right] \right]\right\}\, \rmd t_1.\nonumber \\
\lambdabel{A3Mag}
\end{eqnarray}
\endnumparts
where $\hat V_{\textnormal I}(t)$ is given by (\ref{VI}).
Thus the time evolution of the state of the system can be derived by calculating $\hat at A_1(t), \hat at A_2(t), \dots$, and subsequently constructing
 $\hat U_{\textnormal I}(t)$, which ultimately enables us to calculate $\rho_{\thetaxtnormal I}(t)$. Then it is straight forward to calculate the IP-Wigner function and present a classical correspondence to $\hat U_{\textnormal I}(t)$. We will ultimately transform the IP-density matrix back to the Schr\"{o}dinger\ picture $\hat{\rho}_{\thetaxtnormal{S}}(t)$ with the help of $\hat at U(t)$ as
\betagin{equation}\lambdabel{Int-Sch-state}
\hat{\rho}_{\thetaxtnormal{S}}(t) = \hat at U_0(t)  \hat{\rho}_{\thetaxtnormal I}(t) \hat at U^\dagger_0(t),
\end{equation}
and study the corresponding SP-Wigner function.

For the first two term $\hat A_1(t)$, and $\hat A_2(t)$ we have
\betagin{eqnarray}\lambdabel{A1}
\hat at A_1(t) = \rmi \left[ \nu(t) \hat ad - \nu^*(t) \hat a +\gO\left(\alpha^* \hat ad - \alpha \hat a \right)\right],
\end{eqnarray}
where $\nu(t)$ is defined as
\betagin{equation}\lambdabel{nu}
\nu(t) \equiv - \frac{g}{\Omega} \alpha^* e^{\rmi \Omega t},
\end{equation}
with $\alpha^{*}$ being the complex conjugate of the displacement parameter defined in (\ref{alpha-def}).

Since the commutator $[\hat a,\hat ad] = 1$, the commutator of $\hat V_{\textnormal I}(t)$ at two different times turns into a purely imaginary function of time. Hence the higher order terms in Magnus series given by (\ref{Magnus}) vanish and we are able to write down the exact solution for the equations of motion. This conclusion justifies our choice of $\hat H_1(t)$ and enables us to illustrate the classical correspondence of the IP in phase space.
For $A_2(t)$ we get
\betagin{eqnarray}\lambdabel{A2}
A_2(t) = \frac{g^2 \vert \alpha \vert^2}{\Omega^2} \left(\Omega t - \sin\Omega t \right).
\end{eqnarray}
As a result $\hat U_{\textnormal I}(t)$ is given by
\betagin{eqnarray}\lambdabel{UI}
\hat U_{\textnormal I}(t) = \rme^{- \rmi \frac{g^2\vert \alpha \vert^2}{\Omega^2}\left( \Omega t - 2 \sin\Omega t \right)} \hat D \left(\nu(t)\right) \hat D (\gO \alpha^*).
\end{eqnarray}

\betagin{figure}
\subfloat[$\Omega t = 0$]{\lambdabel{fig:22}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig22.eps}}
\subfloat[$\Omega t=\frac{\hat at pi}{4}$]{\lambdabel{fig:23}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig23.eps}}
\subfloat[$\Omega t=\frac{\hat at pi}{2}$]{\lambdabel{fig:24}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig24.eps}}\\
\subfloat[$\Omega t=\frac{3 \hat at pi}{4}$]{\lambdabel{fig:25}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig25.eps}}
\subfloat[$\Omega t=\hat at pi$]{\lambdabel{fig:26}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig26.eps}}
\subfloat[$\Omega t=\frac{5\hat at pi}{4}$]{\lambdabel{fig:27}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig27.eps}}\\
\subfloat[$\Omega t=\frac{3 \hat at pi}{2}$]{\lambdabel{fig:28}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig28.eps}}
\subfloat[$\Omega t=\frac{7\hat at pi}{4}$]{\lambdabel{fig:29i}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig29.eps}}
\subfloat[$\Omega t=2 \hat at pi$]{\lambdabel{fig:29j}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig29f.eps}}
\caption{Time-evolution of the $\frac{1}{\rme}$-contour of the IP-Wigner function  ($\sigma_x = 1, \sigma_p = \hat alf, \mu_x = -2, a=1, b= -1, \Omega =3, g = 5$) in one period. The evolution corresponds to its global counterclockwise parallel transformation about the origin on the dotted circle described by equation (\ref{breath}).}
\lambdabel{fig.2}
\end{figure}

For a system that is initially in an ideal squeezed state described by (\ref{int_ro_ini}), the state at later time $t$ in the interaction picture reads
\betagin{equation}\lambdabel{int_ro_t}
\hat{\rho}_{\thetaxtnormal I}(t) = \hat D \left(\nu(t)\right) \hat D (\frac{g}{\Omega} \alpha^*)\hat{\rho}(0)\hat D (- \frac{g}{\Omega} \alpha^*)\hat D \left(- \nu(t)\right), 
\end{equation}
where we substituted Eq. (\ref{UI}) for $\hat U_{\textnormal I}(t)$ in (\ref{int-w-sol}).
Knowing the state, we can construct the IP-Wigner function following the same line of calculations  as in section~\ref{Sec:Passive} and get

\betagin{eqnarray}
W^{\thetaxtnormal{I}}_{1}(x,p;t) = \frac{1}{\hat at pi}\exp{\left\{- \hat alf \left[ \bi r - \bi {\bmu_{\thetaxtnormal I}}(t) \right]^T \cdot \Xi^{- 1} \cdot \left[ \bi r - \bi{\bmu}_{\textnormal I}(t) \right]\right\}},
\end{eqnarray}
where the suffix $\thetaxtnormal I$ represents the IP, and the index $1$ represents the evolution of the Wigner function the takes place in the presence of the interaction Hamiltonian $\hat H_1(t)$ given by (\ref{H1_def}). The vector $\bi{\bmu}_{\textnormal I}(t)$ is defined as
\betagin{equation}\lambdabel{r_int_mean}
\bi{\bmu}_{\textnormal I}(t)=\left(\betagin{array}{c}\mu_{x} \\0\end{array}\right) + \gO \left(\betagin{array}{c}a \\- b\end{array}\right) -\gO R(\Omega t) \left(\betagin{array}{c}a \\- b\end{array}\right).
\end{equation}
In order to see how the Wigner function evolves in the IP, we look at the equation of motion of its centroid given by
\betagin{equation}\lambdabel{breath}
\left(\lambdangle x_{\thetaxtnormal I}(t) \rangle - \mu_x\right)^2 + \left( \lambdangle p_{\thetaxtnormal I}(t) \rangle \right)^2 = 2 \left(\gO \right)^2 (a^2 + b^2) (1 - \cos \Omega t).
\end{equation}
As it is shown in figure~\ref{fig.2}, the time-evolution of the $\frac{1}{\rme}$-contour of the IP-Wigner function corresponds to a displacement of the centroid of its elliptical distribution in phase space together with a counterclockwise parallel transformation of the distribution on a circle with a time-dependent radius described by equation (\ref{breath}). There are two important points that can be inferred from figure~\ref{fig.2}: ${\it i)}$ the squeezing direction remains invariant during the evolution in the interaction picture, which in our example was in the $x-$direction. ${\it ii)}$ The evolution of the IP-Wigner function corresponds to the displacement of its centroid around a circle as it was the case for the free Hamiltonian evolution of the free evolution of SP-Wigner function of an ideal squeezed state shown in figure \ref{fig.1}. The reason is that the IP-unitary transformation $\hat U_{\textnormal I}(t)$ is defined as two successive displacement operators apart from a phase factor rather than a rotation. It is exactly  because of this property that for the Hamiltonian (\ref{interaction}) choosing the IP makes the study of the time-evolution of the state much simpler. We will see in the next part that both of the above mentioned properties get spoiled when we transform back to the SP. 

\subsection{Back to the Schr\"{o}dinger picture}\lambdabel{Sec:back}

Transforming back to the SP amounts to performing the unitary transformation  
\betagin{equation}\lambdabel{rho_int_schro}
\hat{\rho}(t) = \hat U_0\hat{\rho}_{\thetaxtnormal I}(t) \hat U_0^{\dagger},
\end{equation}
which as we showed in section~\ref{Sec:Passive} corresponds to a global and a local  rotations of the corresponding SP-Wigner function through angle $\omega_0 t$.
Then the symmetric characteristic function is straightforward to establish. As a result the SP-Wigner function of the system in the presence of the interaction Hamiltonian $\hat H_1(t)$ is given by

\betagin{eqnarray}\lambdabel{backschrowig}
W^{\thetaxtnormal{SP}}_1(x,p;t) = \frac{1}{\hat at pi} \exp{\left\{- \hat alf \left[ \bi r - \bi{\bmu}_{\thetaxtnormal S}(t)\right]^{T}\cdot \Gamma^{-1}\cdot \left[ \bi r - \bi{\bmu}_{\thetaxtnormal S}(t) \right] \right\}}.
\end{eqnarray}
The $2 \times 2$ matirx $\Gamma^{-1}$ is given by (\ref{Gamma_def}) and the mean value vactor $\bi{\bmu}_{\thetaxtnormal S}(t)$ is defined as the clockwise rotation of the initial displacement about the origin through angle $\omega_0 t$ plus two types of transformations of the interaction-induced displacement about the origin: a clockwise rotation through angle $\omega_0t$ and a counterclockwise rotation through angle $\omega_1 t$:
\betagin{equation}\lambdabel{meanr_def}
\bi{\bmu}_{\thetaxtnormal S}(t) = R(- \omega_{0}t)\left(\betagin{array}{c}\mu_x\\0\end{array}\right) + \gO\left[ R(- \omega_0 t) \left(\betagin{array}{c}a \\- b\end{array}\right) -  R(\omega_1 t) \left(\betagin{array}{c}a \\- b\end{array}\right) \right].
\end{equation}

\betagin{figure}
  \centering
  \subfloat[$\omega_1 = \frac{1}{3}$]{\lambdabel{fig:o1_3}\includegraphics[width=0.3\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig48.eps}}                
  \subfloat[$\omega_1 =  \frac{2}{3}$]{\lambdabel{fig:o2_3}\includegraphics[width=0.3\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig49.eps}}
  \subfloat[$\omega_1 =  \frac{3}{5}$]{\lambdabel{fig:o3_5}\includegraphics[width=0.31\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig50.eps}}\\
  \subfloat[$\omega_1 = 3$]{\lambdabel{fig:o3}\includegraphics[width=0.3\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig51.eps}}                
  \subfloat[$\omega_1 = 4$]{\lambdabel{fig:o4}\includegraphics[width=0.3\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig52.eps}}
  \subfloat[$\omega_1 = 4 \frac{5}{7}$]{\lambdabel{fig:o33_5}\includegraphics[width=0.313\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig53.eps}}\\
  \subfloat[$\omega_1 = 5$]{\lambdabel{fig:o5}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig54.eps}}                
  \subfloat[$\omega_1 =  6 \frac{4}{5}$]{\lambdabel{fig:o34_5}\includegraphics[width=0.3\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig55.eps}}
  \subfloat[$\omega_1 = 9$]{\lambdabel{fig:o9}\includegraphics[width=0.313\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig56.eps}}
  \caption{The trajectory of the center of SP-Wigner function in the presence of the interaction for different values of $\omega_1$. $\omega_0$ is fixed to $1$ in arbitrary units. The parameters are fixed at $\sigma_x = 1, \sigma_p = \hat alf, \mu_x = -2, a=1, b= -1, g=5$.}
  \lambdabel{fig.3}
\end{figure}
The first point we notice by looking at the Wigner function of (\ref{backschrowig}) is that when we go back to the SP the local and global rotations of the Wigner function appear again and this is due to the presence of the unitary transformation $\hat U_0$ . Furthermore, the global rotation of the centroid of the distribution occurs not on a circle, like in the free evolution case illustrated in figure~\ref{fig.2}, but on a more complex trajectory given by
\betagin{equation}\lambdabel{breathS}
\left(\lambdangle x_{\thetaxtnormal S}(t) \rangle - \mu_x \cos \omega_0 t\right)^2 + \left( \lambdangle p_{\thetaxtnormal S}(t) \rangle + \mu_x \sin \omega_0 t \right)^2 = 2 \left(\gO \right)^2 (a^2 + b^2) (1 - \cos \Omega t).
\end{equation}
Equation (\ref{breathS}) resembles a glissette~\cite{Roulette}. A glissette is defined as the locus of a generator point which is moving along a given curve. The generator point in (\ref{breathS}) is $(\mu_x \cos \omega_0 t, - \mu_x \sin \omega_0 t)$ described as the center of the circle in the right hand side of (\ref{breathS}). This point is rolling  along the time-dependent curve given by the left hand side of (\ref{breathS}). The resulting glissette is shown in figure~\ref{fig.3} for different values of $\Omega$.
\betagin{figure}
\subfloat[$\omega_1 t = 0$]{\lambdabel{fig:39}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig39.eps}}
\subfloat[$\omega_1 t = \frac{\hat at pi}{4}$]{\lambdabel{fig:40}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig40.eps}}
\subfloat[$\omega_1 t = \frac{\hat at pi}{2}$]{\lambdabel{fig:41}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig41.eps}}\\
\subfloat[$\omega_1 t= \frac{3 \hat at pi}{4}$]{\lambdabel{fig:42}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig42.eps}}
\subfloat[$\omega_1 t = \hat at pi$]{\lambdabel{fig:43}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig43.eps}}
\subfloat[$\omega_1 t = \frac{5 \hat at pi}{4}$]{\lambdabel{fig:44}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig44.eps}}\\
\subfloat[$\omega_1 t = \frac{3 \hat at pi}{2}$]{\lambdabel{fig:45}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig45.eps}}
\subfloat[$\omega_1 t = \frac{7 \hat at pi}{4}$]{\lambdabel{fig:46}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig46.eps}}
\subfloat[$\omega_1 t = 2 \hat at pi$]{\lambdabel{fig:47}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig47.eps}}
\caption{Time-evolution of the $\frac{1}{\rme}$-contour of $W^{\thetaxtnormal SP}_{1}(x,p;t)$  given by (\ref{backschrowig}) ($\sigma_x = 1, \sigma_p = \hat alf, \mu_x = -2, a=1, b= -1, \omega_1 = 2, \omega_0 =1, g = 5$). The evolution is a combination of local rotations of the center of the ellipse through two different angles $\omega_{0}t$ and $\omega_{1}t$ and a global rotation of the ellipse through angle $\omega_{0}t$. This leads to  the ``dancing'' of $W^{\thetaxtnormal{SP}}_{1}(x,p;t)$. The dotted curve describes the motion of the centroid given by (\ref{breathS}).}
\lambdabel{fig.4}
\end{figure}

On top of the rotations induced by $\hat U_0$, there are two more rotations of the driven displacement given by $\alpha^{*}$ due to the presence of the interaction. ${\it i)}$ Its clockwise rotation through angle $\omega_{0}t$; ${\it ii)}$ its counterclockwise rotation through angle $\omega_1t$. As a result, the SP-Wigner function in the presence of the interaction ``dances" in a complicated manner, in comparison with its parallel transform in IP, in phase space. This is illustrated  in figure~\ref{fig.4} for a fixed value of $\Omega$. A comparison between the figure \ref{fig.2} and figure \ref{fig.4} makes it clear why in studying the evolution of the state of our system of interest the interaction picture is preferred. 

Choosing the counterclockwise rotation of the $xp$-plane given by (\ref{r-r'}) transforms the SP-Wigner function into the IP-Wigner function in the local $x'p'$ frame. This corresponds to performing the unitary transformation $\hat U_0$ on the observables in Hilbert space. Correspondingly, in phase space replacing $\bi r$ in (\ref{backschrowig}) with $R(- \omega_0\,t)\,\bi{r'}$ yields the same result. Then the vector $\bi{r} - \bi{\bmu}_{\thetaxtnormal S}(t)$ in the local reference frame reads 
\betagin{eqnarray}\lambdabel{rminmus_new}
\fl\bi{r} - \bi{\bmu}_{\thetaxtnormal S}(t) \rightarrow R(- \omega_0\, t) \left(\bi{r'} - \bi{\bmu}\right) - \gO R(- \omega_0\, t) \left(\betagin{array}{c}a \\- b\end{array}\right) + \gO R(\omega_1 t) \left(\betagin{array}{c}a \\- b\end{array}\right).
\end{eqnarray}
Substituting (\ref{rminmus_new}) into the expression for the SP-Wigner function given by (\ref{backschrowig}) and taking into account that $\Gamma^{-1}$ is defined as $R(-\omega_0\,t) \Xi^{-1} R(\omega_0\,t)$, we get
\betagin{eqnarray}\lambdabel{backshro_int_new}
\fl W^{\thetaxtnormal{SP}}_{1}(x,p;t) \rightarrow W^{\thetaxtnormal{I}}_{1}(x',p';t) = \frac{1}{\hat at pi} \exp{\left[- \hat alf \left( \bi{r'} - \bi{\bmu}_{\thetaxtnormal I}(t)\right)^{T}\cdot \Xi^{-1}\cdot \left( \bi{r'} - \bi{\bmu}_{\thetaxtnormal I}(t) \right) \right]}.
\end{eqnarray}
In fact, by choosing the counterclockwise rotation of the coordinates, we let the observables of the system evolve under the free Hamiltonian $\hat H_0$. Thus the evolution of the state in the new local frame is only governed by the interaction Hamiltonian. In other words, we separate the time-evolution of the state from that of the observables. But this is exactly the definition of the Schr\"{o}dinger interaction picture~\cite{SHIP}.

We can also transform the Wigner function to the Heisenberg interaction picture~\cite{SHIP}. Such transformation may be performed by ${\it i)}$ a counterclockwise rotation of $xp$ frame, i.e. a passive transformation of $xp$-frame to $x'p'$-frame; ${\it ii)}$ displacing the origin of $x'p'$-frame by $\gO \tiny{\left(\betagin{array}{c}a \\ - b\end{array}\right)}$; ${\it iii)}$ counterclockwise rotation of the new origin of the $x'p'$-frame with the frequency $\Omega$. Thus the transformation from $xp$-plane to the new $x''p''$-plane is defined as
\betagin{eqnarray}\lambdabel{r-r''}
\bi r \rightarrow \bi{r''} =  R(\omega_0\,t) \bi r - \gO \left(\betagin{array}{c}a \\- b\end{array}\right) + \gO R(\Omega\,t) \left(\betagin{array}{c}a \\- b\end{array}\right).
\end{eqnarray}
The last two transformations are precisely compensating the unitary transformation arising from the interaction Hamiltonian. 
It is then straightforward to describe the vector $\bi r - \bi{\bmu}_{\thetaxtnormal S}(t)$  in the $x''p''$-frame as the clockwise rotations of $(\bi{r''} - \bmu)$ through angle $\omega_0 t$:
\betagin{equation}\lambdabel{mufree}
\bi r -\bi{\bmu}_{\thetaxtnormal S}(t) \rightarrow R(- \omega_0 t) (\bi{r''} - \bmu).
\end{equation}
Inserting (\ref{mufree}) into (\ref{backschrowig}) gives us the Wigner function in the local $x''p''$-plane as
\betagin{equation}\lambdabel{W''}
W^{\thetaxtnormal{SP}}_{1}(x,p;t) \rightarrow W^{\thetaxtnormal{SP}}_0(x'',p'') = \frac{1}{\hat at pi} \exp{\left[ - \hat alf \left ( \bi{r''} - \bmu \right)^T \cdot \Xi^{-1} \cdot \left( \bi{r''}- \bi{\bmu}\right) \right]}.
\end{equation}
We notice that in the local $x''p''-$frame, only $x''$ and $p''$ are evolving in time according to $\hat H(t)$ and the state does not explicitly depend on time.

\subsection{The case $\Omega = 0$}\lambdabel{sec:Omzero}

It is interesting to see what happens when the driving frequency of the interaction Hamiltonian given in equation (\ref{H1_def}) is equal to that of $\hat H_{0}$ with opposite sing. The first thing we notice is that the interaction picture potential given by equation (\ref{VI}) is independent of time. As a result the Wigner function is described by
\betagin{eqnarray}\lambdabel{Wigne_int_Om0}
W^{\thetaxtnormal{I}}_{\Omega = 0}(x,p;t) = \frac{1}{\hat at pi} \exp{\left[- \hat alf \left( \bi r - \bi{\bmu}_{\textnormal I}(t)\vert_{\Omega=0}\right)^{T}\cdot\Xi \cdot \left( \bi r - \bi{\bmu}_{\textnormal I}(t)\vert_{\Omega=0}\right) \right]},
\end{eqnarray}
where the components of $\bi{\bmu}_{\textnormal I}(t)\vert_{\Omega=0}$ are given by
\numparts
\betagin{eqnarray}\lambdabel{xmean_int_Om0}
\lambdangle x_{\thetaxtnormal I}(t) \rangle\vert_{\Omega = 0}& = \mu_{x} - g b t,\\
\lambdabel{pmean_int_Om0}
\lambdangle p_{\thetaxtnormal I}(t) \rangle\vert_{\Omega = 0} &= - g a t.
\end{eqnarray}
\endnumparts
\betagin{figure}
\centering
 \includegraphics[width=0.3\thetaxtwidth,clip=]{Mehmani_quantuminteractionpictureFig57.eps}
 \caption{A snapshot of the interaction picture Wigner function $W^{\thetaxtnormal{I}}_{1}(x,p;t)\vert_{\Omega = 0}$ at some time $t$ with the parameters $\sigma_x = 1, \sigma_p = \hat alf,\mu_x = -2, a=1, b= -1, \omega_1 = - \omega_0 , g = 5$. The dotted ellipse represent the Winger function at time $t = 0$. The dashed line is the trajectory of the center of the ellipse.}
\lambdabel{fig.5}
\end{figure}
Thus the equation of motion for the centroid of the ellipse will be a line instead of a circle and the Wigner function will only be parallel displaced along this line in time. The displacement takes place on a line with slope $a/b$, as it is shown in figure~\ref{fig.5}, and is given by the following equation
\betagin{equation}\lambdabel{meanx,meanp_int_Om0}
\lambdangle p_{\thetaxtnormal I}(t) \rangle\vert_{\Omega = 0} = \frac{a}{b} \left( \lambdangle x_{\thetaxtnormal I}(t) \rangle\vert_{\Omega = 0} - \mu_{x} \right).
\end{equation}

\subsection{Quadratic Hamiltonian}\lambdabel{sec:qua}

\betagin{figure}
\subfloat[$t = 0$]{\lambdabel{fig:70}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig70.eps}}
\subfloat[$ t = \frac{\hat at pi}{4}$]{\lambdabel{fig:71}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig71.eps}}
\subfloat[$t = \frac{\hat at pi}{2}$]{\lambdabel{fig:72}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig72.eps}}\\
\subfloat[$t= \frac{3 \hat at pi}{4}$]{\lambdabel{fig:73}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig73.eps}}
\subfloat[$t = \hat at pi$]{\lambdabel{fig:74}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig74.eps}}
\subfloat[$ t = \frac{5 \hat at pi}{4}$]{\lambdabel{fig:75}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig75.eps}}\\
\subfloat[$t = \frac{3 \hat at pi}{2}$]{\lambdabel{fig:76}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig76.eps}}
\subfloat[$t = \frac{7 \hat at pi}{4}$]{\lambdabel{fig:77}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig77.eps}}
\subfloat[$t = 2 \hat at pi$]{\lambdabel{fig:78}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig78.eps}}
\caption{Time-evolution of the $\frac{1}{\rme}$-contour of  $W^{\thetaxtnormal {IP}}_{2}(x,p;t)$  given by (\ref{Wqu:Int}) in a time period of $2 \hat at pi$ ($\sigma_x = 1, \sigma_p = \hat alf, \omega_0 =1, \kappa= 0.1$). The amount of the squeezing changes in time. We see that in the course of time the initially squeezed in the $p$-direction distribution function becomes squeezed in the $x$-direction due to the presence of the quadratic interaction term. The dashed circle represents the vacuum coherent  state.}
\lambdabel{fig.6}
\end{figure}

So far we have considered an interaction Hamiltonian that is linear in the creation and annihilation operators. From figures \ref{fig.1}--\ref{fig.5} it is apparent that such a Hamiltonian preserves the covariance matrix of the Gaussian distribution of the Wigner function, i.e., an initial coherent/ squeezed state remains coherent/squeezed under time evolution, although it might perform different types of rotations. Mathematically, this means the covariance matrix given by (\ref{covmatrix_int}) is invariant under a linear time-dependent interaction term. In this section we consider a quadratic interaction term in the Hamiltonian to illustrate the deformation of the covariance matrix. We show that in this case the contour of the Wigner function initially squeezed in the $p$-direction expands into a coherent state and then becomes squeezed in the $x$-direction. We choose the quadratic Hamiltonian to be
\betagin{eqnarray}
\lambdabel{Hqu}
\hat H(t) = \omega_0 \hat ad \hat a + \rmi \kappa \left[ e^{ 2 \rmi \omega_0 t} \hat a^2-  e^{- 2 \rmi \omega_0 t} (\hat ad)^2  \right],
\end{eqnarray}
where $\kappa$ is a real coupling constant. A physical realization of the above Hamiltonian may be the interaction of a coherent light beam with a non-linear optical medium. In that case, $\kappa$ contains the non-linear susceptibility and the pump field amplitude. The group-theoretical approach to the time-evolution of the Wigner function with a quadratic potential has been immensely studied~\cite{Wod,Gerry,Dattoli2,Dattoli1,Han,Orlowski,Vourdas,Zoubi}.
In what follows we use the Lie algebra of $\mbox{SU(1,1)}$\cite{Wod,Gerry} realized in terms of $\hat a$ and $\hat ad$,
\betagin{eqnarray}\lambdabel{K0pm:def}
\hat{K}_0 = \hat alf( \hat ad \hat a + \hat alf),\qquad \hat{K}_+ = \hat alf (\hat ad)^2,\qquad \hat{K}_- = \hat alf \hat a^2,
\end{eqnarray}
such that 
\betagin{eqnarray}\lambdabel{K0pm:comm}
\left[\hat{K}_0, \hat{K}_{\hat at pm} \right] = \hat at pm \hat{K}_{\hat at pm},\qquad \left[\hat{K}_-,\hat{K}_+ \right] = 2 \hat{K}_0.
\end{eqnarray}
Thus in terms of these operators we have a Hamiltonian with linear interaction term
\betagin{eqnarray}
\lambdabel{Hqu:K}
\hat H(t) = 2 \omega_0 \hat{K}_0  + \rmi \kappa\left(e^{ 2 \rmi \omega_0 t} \hat{K}_- - e^{- 2 \rmi \omega_0 t} \hat{K}_+\right).
\end{eqnarray}
It is straight forward to calculate the unitary time-evolution operator in the interaction picture using (\ref{UMagnus}). In this case the interaction picture potential $\hat V_{\textnormal I}(t)$ becomes independent of time and as a result the Magnus series contains only the first term $\hat A_1(t)$ given by (\ref{Magnus}). Thus $\hat U_{\textnormal I}(t)$ reads
\betagin{eqnarray}\lambdabel{Uqu:int}
\hat U_{\textnormal I}(t) = \rme^{\kappa t (\hat{K}_- - \hat{K}_+)}.
\end{eqnarray}
$\hat U_{\textnormal I}(t)$ resembles a squeezing operator with a time-dependent parameter, as a generalization of the time-dependent displacement operator in (\ref{UI}) for the linear interaction term in terms of $\hat a$ and $\hat ad$ calculated in the previous section.

\betagin{figure}
\subfloat[$\omega_0 t = 0$]{\lambdabel{fig:80}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig80.eps}}
\subfloat[$\omega_0 t = \frac{\hat at pi}{4}$]{\lambdabel{fig:81}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig81.eps}}
\subfloat[$\omega_0 t = \frac{\hat at pi}{2}$]{\lambdabel{fig:82}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig82.eps}}\\
\subfloat[$\omega_0 t= \frac{3 \hat at pi}{4}$]{\lambdabel{fig:83}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig83.eps}}
\subfloat[$\omega_0 t = \hat at pi$]{\lambdabel{fig:84}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig84.eps}}
\subfloat[$\omega_0 t = \frac{5 \hat at pi}{4}$]{\lambdabel{fig:85}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig85.eps}}\\
\subfloat[$\omega_0 t = \frac{3 \hat at pi}{2}$]{\lambdabel{fig:86}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig86.eps}}
\subfloat[$\omega_0 t = \frac{7 \hat at pi}{4}$]{\lambdabel{fig:87}\includegraphics[width=0.3\thetaxtwidth]{Mehmani_quantuminteractionpictureFig87.eps}}
\subfloat[$\omega_0 t = 2 \hat at pi$]{\lambdabel{fig:88}\includegraphics[width=0.313\thetaxtwidth]{Mehmani_quantuminteractionpictureFig88.eps}}
\caption{Time-evolution of the contour of $\frac{1}{\rme}$ value of $W^{\thetaxtnormal {SP}}_{2}(x,p;t)$  given by (\ref{Wqu:Sch}) in one period of $ \omega_0 t$ ($\sigma_x = 1, \sigma_p = \hat alf, \kappa= 0.1, \omega_0 = 1$). The distribution function starts to expand in the $x$-direction and then squeeze in the $p$-direction. On top of that it also performs a clockwise rotation around the origin. The dashed circle represents the vacuum coherent state.}
\lambdabel{fig.7}
\end{figure}

Now we can calculate the state of the system at later time $t$ needed for the expectation values of the displacement operator. For the sake of simplicity, we assume the system starts its evolution in the IP from a vacuum squeezed state. This way, the centroid of the distribution coincides with the origin of the phase space. Then the symmetric characteristic function reads
\betagin{eqnarray}
\lambdabel{chiqu:Int}
\Tr\left[\hat \rho_{\textnormal I} (t) \hat D(\chi)\right] = \exp{\left\{- \frac{1}{8} \left[\frac{\hat at xi^2}{\left(\sigma_p \rme^{\kappa t}\right)^2} + \frac{\eta^2}{\left(\sigma_x \rme^{- \kappa t}\right)^2}\right]\right\}}.
\end{eqnarray}
We immediately realize that in the presence of a quadratic interaction the elements of the covariance matrix given by (\ref{covmatrix_int}) undergo a scaling (conformal) transformation. Furthermore, at a certain time the two elements of the covariance matrix become identical thus a squeezed states transforms into a coherent state.
Substituting (\ref{chiqu:Int}) into the definition of the Wigner function given by (\ref{Wignerdef}) and performing two Gaussian integrals over $\hat at xi$ and $\eta$ gives the IP-Wigner function for the quadratic Hamiltonian as
\betagin{eqnarray}
\lambdabel{Wqu:Int}
W^{\textnormal{IP}}_2(x,p;t) = \frac{1}{\hat at pi} \exp{\left[-\hat alf \bi r^T\cdot \Xi ^{-1}_2(t)\cdot \bi r\right]},
\end{eqnarray}
where the index $2$ represents the quadratic interaction term added to the free Hamiltonian. $\Xi^{- 1}_2(t)$ is the ``breathing'' covariance matrix defined as
\betagin{eqnarray}\lambdabel{covmatrixqu:int}
\Xi^{-1}_2(t) = \left(\betagin{array}{cc}\frac{\rme^{2 \kappa t}}{\sigma^2_x} & 0 \\0 &\frac{\rme^{-2 \kappa t}}{\sigma^2_p}\end{array}\right).
\end{eqnarray}
This is illustrated in figure~\ref{fig.6}.

As it is showed in section~\ref{Sec:back} the transformation of the time-dependent density matrix back to the SP, amounts to a global clockwise rotation of the contour of the SP-Wigner function in phase space. In the case of quadratic interaction term this is straight forward to prove and yields to
\betagin{eqnarray}\lambdabel{Wqu:Sch}
W^{\textnormal{SP}}_2(x,p;t) = \frac{1}{\hat at pi} \exp{-\hat alf \left[\bi r^T\cdot\Gamma^{-1}_2(t) \cdot \bi r\right]},
\end{eqnarray}
where $\Gamma^{-1}_2(t)$ is defined as the clockwise rotated breathing covariance matrix,
\betagin{eqnarray}\lambdabel{Gammaqu}
\Gamma^{-1}_2(t)\equiv R (- \omega_0 t)\cdot \Xi ^{-1}_2(t) \cdot R(\omega_0 t).
\end{eqnarray}
This is shown in figure~\ref{fig.7}.
It is worth noticing that since the centroid of the Wigner function coincides with the origin of the phase space, the local counterclockwise rotation of the elliptical distribution around the centroid is canceled out by the global clockwise rotations around the origin as illustrated in figure \ref{fig.7}. However, this is not the case for an ideal squeezed state as it is illustrated in figure \ref{fig.1}.

\section{Summary}
We presented a pictorial representation of time evolution of a harmonic oscillator driven linearly and quadratically in Schr\"{o}dinger, \textnormal{Heisenberg}, and \textnormal{Interaction}\ picture (both Schr\"{o}dinger\ interaction picture, in which the evolution of the state is due to the interaction potential, and \textnormal{Heisenberg}\ interaction picture, where the the free evolution is contained in the state ). This has been done by employing the Wigner-Weyl representation to map the density matrices and observables in Hilbert space on distribution functions and variables, correspondingly, in phase space. As an example, we took the initial state of the system an ideal squeezed one. We showed that in the presence of a linear interaction term the time evolution of the corresponding Wigner function in the Schr\"{o}dinger\ picture amounts to different types of local and global rotations. The interaction picture time-evolution is mapped to a parallel transformation of the distribution function on a circle around the origin of the phase space. In transforming back to the Schr\"{o}dinger\ picture, we showed that the evolution of the Wigner function becomes more complicated due to the fact that it performs different types of rotations with difference frequencies. However, one can transform from one picture to another by performing passive and local transformations of the reference frame. Adding a quadratic interaction term to the Hamiltonian amounts to squeezing of the state with a time-dependent squeezing parameter. Classical picture of such interaction corresponds to the breathing of the Wigner function contours.

\section*{References}
\betagin{thebibliography}{22}

\bibitem{Parisi} Auletta G and Parisi G 2001 {\it Foundations and Interpretation of Quantum Mechanics} (James Bennett Pty Ltd) p~ 48

\bibitem{Mandel} Mandel L and Wolf  E 1995 {\it Optical Coherence and Quantum Optics} (Cambridge University Press) 1042--45

\bibitem{Case08} Case W B 2008 ``Wigner functions and Weyl transforms for pedestrians'' {\it Am. J. Phys.} {\bf 76}, 937--46

\bibitem{Hillery} Hillery M, O'Connnell R F, Scully M O and Wigner E P  1984``Distribution functions in Physics: Fundamentals''  {\it Phys. Rep.}  {\bf 106}(3), 121--67

\bibitem{Agarwal} Agarwal G S 1987 ``Wigner-function description of quantum noise in interferometers'' {\it J. Mod. Opt.} {\bf 34}(6/7) 909--21

\bibitem{Schleich} Schleich W P 2001 {\it Quantum Optics in Phase Space} 1st. ed. (Wiley-VCH Verlag Berlin GmbH) 321--42

\bibitem{Kim90} Kim Y S and Wigner E P 1990 ``Canonical transformation in quantum mechanics'' {\it Am. J. Phys.} {\bf 58}  439--48

\bibitem{Kim88}Han D, Kim Y S and Noz M E 1988 ``Linear canonical transformations of coherent and squeezed states in the Wigner phase space'' {\it Phys. Rev.} A {\bf 37} 807--14

\bibitem{Kim89}Han D, Kim Y S and Noz M E 1989 ``Linear canonical transformations of coherent and squeezed states in the Wigner phase space. II. Quantitative analysis'' {\it Phys. Rev.} A {\bf 40}, 902--12

\bibitem{Glauber} Glauber R J 1963 ``Coherent and incoherent states of the radiation field'' {\it Phys. Rev.} {\bf 131} 2766--88

\bibitem{Barnett} Barnett S M and Radmore P M 1997 {\it Methods in Theoretical Quantum Optics} (Clarendon Press,Oxford) p~ 106

\bibitem{Scully} Scully M O and Zubairy M S 2002 {\it Quantum Optics} (Cambridge University Press) p~ 51

\bibitem{Yuen76}Yuen H P  1976 ``Two-photon coherent states of the radiation field'' {\it Phys. Rev.} A {\bf 13} 2226--43

\bibitem{Caves81} Caves C M 1981 ``Quantum-mechanical noise in an interferometer'' {\it Phys. Rev.} D {\bf 23} 1693--708

\bibitem{Walls83} Walls D F 1983  ``Squeezed states of light'' {\it Nature} {\bf 306} 141--6

\bibitem{Tapia93} Mu\~{n}oz-Tapia R  1993``Quantum mechanical squeezed state'' {\it Am. J. Phys.} {\bf 61} 1005--8

\bibitem{Magnus} Magnus W 1954 ``On the exponential solution of differential equations for a linear operator''  {\it Communications on Pure and Applied Mathematics} {\bf 7} 649--73. For a general review on the Magnus series both from the Mathematical and the Physical point of view see Blanes S, Casas F, Oteo J A and Ros J  2009 ``The Magnus expansion and some of its applications'' {\it Phys. Rep} {\bf 470} 151--238

\bibitem{Roulette} Walker G 1939 The theory of roulettes and glissettes
{\it National Mathematics Magazine} {\bf 13} No. 5, 223--29

\bibitem{SHIP} Barnett S M and Radmore P M 1997 {\it Methods in Theoretical Quantum Optics} (Clarendon Press,Oxford) p~ 14

\bibitem{Wod} W\`{o}dkiewicz K and Eberly J H  1985 ``Coherent states, squeezed fluctuations, and the SU(2) and SU(1,1) groups in quantum-optics applications'' {\it J. Opt. Soc. Am. B} Mar 1;{\bf2(3)} 458--66. 

\bibitem{Gerry} Gerry C C 1985 ``Dynamics of SU(1,1) coherent states'' {\it Phys. Rev. A} {\bf31} 2721--23.

\bibitem{Dattoli2} Dattoli G and Torre A. 1988 ``Algebraic view to the quantum anharmonic oscillator'' {\it Phys. Rev. A} {\bf 37(5)} pp.1571--75.

\bibitem{Dattoli1} Dattoli G, Solimeno S, Torre A. 1987  `` Algebraic view of the optical propagation in a nonhomogeneous medium''. {\it Phys. Rev. A.} {\bf 35} 1668--72. 

\bibitem{Han} Han D, Kim YS, Noz ME. 1989 ``Linear canonical transformations of coherent and squeezed states in the Wigner phase space. II. Quantitative analysis''  {\it Phys. Rev. A.} {\bf 40(2)} 902--12. 

\bibitem{Orlowski} Orlowski A, W\'{o}dkiewicz K  1990 ``On the SU(1, 1) Phase-space Description of Reduced and Squeezed Quantum Fluctuations''  {\it Journal of Modern Optics}  {\bf 37(3)} 295--301. 

\bibitem{Vourdas} Vourdas A 1990  ``SU(2) and SU(1,1) phase states'' {\it Phys. Rev. A.}  {\bf 41}1653--61. 

\bibitem{Zoubi} Zoubi H and Ben-Aryeh Y. 1998 ``The evolution of harmonic oscillator Wigner functions described by the use of group representation'' {\it Quantum Semiclass. Opt.} {\bf 10} 447--58.

\end{thebibliography}

\end{document}