---
abstract: 'We study the evolution of the hybrid entangled states in a bipartite (ultra) strongly coupled qubit-oscillator system. Using the generalized rotating wave approximation the reduced density matrices of the qubit and the oscillator are obtained. The reduced density matrix of the oscillator yields the phase space quasi probability distributions such as the diagonal $P$-representation, the Wigner $W$-distribution and the Husimi $Q$-function. In the strong coupling regime the $Q$-function evolves to uniformly separated macroscopically distinct Gaussian peaks representing ‘kitten’ states at certain specified times that depend on *multiple* time scales present in the interacting system. For the ultra-strong coupling realm a large number of interaction-generated modes arise with a complete randomization of their phases. A stochastic averaging of the dynamical quantities sets in while leading to the decoherence of the system. The delocalization in the phase space of the oscillator is studied by using the Wehrl entropy. The negativity of the $W$-distribution, while registering its departure from the classical states, allows us to compare the information-theoretic measures such as the Wehrl entropy with the Wigner entropy. Other features of nonclassicality such as the existence of the squeezed states and appearance of negative values of the Mandel parameter are realized during the course of evolution of the bipartite system. In the parametric regime studied here these properties do not survive after a time-averaging process.'
---

\
R. Chakrabarti$^{\dagger}$ and V. Yogesh$^{\ddagger}$\

*$^{\dagger}$Chennai Mathematical Institute, H1 SIPCOT IT Park,\
Siruseri, Kelambakkam 603 103, India\
$^{\ddagger}$ Department of Theoretical Physics, University of Madras,\
Maraimalai Campus, Guindy, Chennai 600 025, India\
*

Introduction
============

The quantum mechanical model that describes a two-level system (qubit) interacting with an electromagnetic field mode (oscillator) in a cavity possesses a rich theoretical structure. Its dynamical features, particularly in the regime of a weak qubit-oscillator coupling and also for a small detuning between the qubit and the oscillator frequencies, are described by the exactly solvable Jaynes-Cummings model \[[@JC1963]\] that employs the rotating wave approximation. Recently, however, the stronger coupling domain, where the rotating wave approximation no longer holds, has been experimentally investigated. Various novel realizations such as a nanomechanical resonator capacitively coupled to a Cooper-pair box \[[@ABS2002]\], a quantum semiconductor microcavity undergoing excitonic transitions \[[@ALT2009]\], a flux-biased superconducting quantum circuit that uses large nonlinear inductance of a Josephson junction to produce ultra-strong coupling with a coplanar waveguide resonator \[[@ND2010]\] fall in this category. The superconducting qubits and circuits, in specific, offer much flexibility in selecting the control parameters. Consequently, they are regarded as suitable candidates for the quantum simulators \[[@BAN2011]-[@GAN2014]\]. Moreover, integrated hybrid quantum circuits involving atoms, spins, cavity photons, and superconducting qubits with nanomechanical resonators have triggered much interest \[[@ZAYN2013]\].

To study the strongly interacting qubit-oscillator system where the Hamiltonian incorporates terms that do not preserve the total excitation number, the authors of Refs. \[[@IGMS2005], [@AN2010]\] have advanced an adiabatic approximation scheme that holds in the parametric realm where the oscillator frequency dominates the qubit frequency. The resultant separation of the time scales facilitates decoupling of the full Hamiltonian into sectors related to each scale, and allows approximate evaluation \[[@IGMS2005]\] of the eigenstates of the system. To augment the parametric domain so that it includes the exact resonance as well as the large detuning regime, a generalization of the rotating wave approximation has been proposed \[[@I2007]\]. This generalization utilizes the basis states obtained via the adiabatic limit, and subsequently, the excitation number conserving argument *à la* the rotating wave approximation is applied \[[@I2007]\] to the Hamiltonian in the new basis. The resultant block diagonalized Hamiltonian now furnishes \[[@I2007]\] energy eigenvalues which are approximately valid for the (ultra) high coupling strength and also for a wide range of the detuning parameter.

The hybrid entangled states of the bipartite coupled qubit-oscillator system provide instances of the entanglement of the microscopic atomic states and, say, the photonic Schrödinger cat states that may be regarded as macroscopic for reasonably large values of the coherent state amplitude. These states play crucial roles in many areas. For instance, they facilitate the non-destructive measurement \[[@Haroche2007]\] of the photon number in a field stored in a cavity. Bell inequality tests involving these qubit-field entangled states have been proposed \[[@PSSAJ2012]\]. The hybrid entanglement involving a discrete and a continuous quantum variable may offer advantages such as achieving near-deterministic quantum teleportation \[[@LJ2013]\]. Moreover, the quantum bus \[[@VL2008]\] framework utilizes the hybrid entanglement where direct qubit-qubit interactions are eliminated while providing a commonly coupled continuous mode that mediates among qubits. Realization of the micro-macro entangled states via a controllable interaction of a single-mode microwave cavity field with a superconducting charge qubit, and the subsequent creation of the superposition of macroscopically distinguishable field modes by virtue of measuring the charge states of the qubit have been proposed \[[@LWN2005]\]. Much current experimental activity \[[@A2014]\] is focused towards generation of such hybrid entanglement. Recently the optical hybrid entanglement has been observed \[[@Jeong2014]\] by the quantum superposition of non-Gaussian operations on distinct modes. Following a measurement based procedure, the hybrid entanglement between two remote nodes residing in the Hilbert spaces of different dimensionality has been established \[[@Morin2014]\]. The superconducting circuits have also been used for controllable and deterministic generation of complex superposition of states \[[@LWN2004]-[@Hofheinz2009]\].

Our objective in the present work is as follows. Within the generalized rotating wave approximation scheme we analytically study the evolution of an initial hybrid entangled state in a coupled qubit-oscillator system. Our study includes both the strong and the ultra-strong coupling domains. Tracing over the complementary degree of freedom we obtain the time-evolution of the reduced density matrices of the qubit and the oscillator, respectively. The qubit reduced density matrix provides the von Neumann entropy of the system that measures the entanglement and the mixedness of the state. The oscillator reduced density matrix, in turn, yields the phase space quasi probability distributions \[[@Schleich2001]\] such as the diagonal $P$-representation, the Wigner $W$-distribution, and the Husimi $Q$-function. The Wehrl entropy \[[@W1978]\] constructed via the $Q$-function measures the delocalization of the oscillator in the phase space. In the strong coupling regime where a quadratic approximation to the effective interaction described by the Laguerre polynomials holds, a long-range quasi periodic time dependence is visible for the Wehrl entropy. This is similar to its behavior in the Kerr-like nonlinear self-interacting photonic models, where it has been observed \[[@MTK1990]-[@MBWI2001]\] that an initial coherent state therein evolves, at rational fractions of the time period, to the superposition of certain macroscopic coherent states popularly known as ‘kitten’ states. While paralleling this phenomenon the current qubit-oscillator interacting model has, however, an important distinction that stems from the existence of interaction-generated multiple time scales. Quantum fluctuations at shorter time scale that signifies the energy exchange between the qubit and the oscillator cause, at rational fractions of the long-range time period, a doubling of the number of coherent state peaks (‘kittens’) while the Wehrl entropy evolves from a local minimum to a local maximum of the short-range undulations. A frequency modulation of the short time scale oscillation is also observed. In the ultra-strong coupling limit the randomization of the phases of a large number of incommensurate quantum modes enforces a statistical equilibrium, and a stabilization of the occupation on the phase space sets in. Another object of study is the Wigner entropy \[[@SKD2012]\] that probes the quantumness of the mixed states linked with the negativity of the $W$-distribution. Our comparison between the Wigner entropy and the Wehrl entropy reveals the close connection of their interrelation with the negativity of the $W$-distribution. Quantum fluctuations in the strong coupling domain induce formations of *almost pure* squeezed states. Using the Mandel parameter \[[@M1979]\] the nonclassicality of the photon statistics is also observed.

The Hamiltonian and the reduced density matrices
================================================

In natural units ($\hbar=1$) the Hamiltonian \[[@IGMS2005]-[@I2007]\] of the coupled qubit-oscillator system reads H = a\^ a + \_x + \_z (a\^ + a), \[qOH\] where the field mode of frequency $\omega$ is described by the annihilation and creation operators ($a, a^{\dagger}| \hat{n} \equiv a^{\dagger} a$), and the two-level atom having a transition frequency $\Delta$ is expressed via the spin variables $(\sigma_{x}, \sigma_{z})$. The qubit-oscillator coupling strength equals $\lambda$. The Fock states $\{\hat{n} |n\rangle = n |n\rangle,\,n = 0, 1,\ldots;\;a \,|n\rangle = \sqrt{n}\,|n - 1\rangle, 
a^{\dagger}\, |n\rangle = \sqrt{n + 1}\,|n + 1\rangle\}$ provide the basis for the oscillator, whereas the eigenstates $\sigma_z |\pm 1\rangle = \pm \,|\pm 1\rangle$ span the space of the qubit. As the Hamiltonian (\[qOH\]) is not known to be exactly solvable, numerous approximation schemes adapted to various ranges of parameters have been advanced. For a small detuning between the oscillator and the qubit frequencies, and also for a weak qubit-oscillator coupling, the dynamical behavior of the interacting system is accurately described \[[@JC1963]\] by the rotating wave approximation. In another approach, the adiabatic approximation \[[@IGMS2005], [@AN2010]\] scheme is found to be appropriate in the large detuning limit ($\Delta \ll \omega$) as it utilizes the difference between the time scales of the slow-moving qubit and that of the fast-moving oscillator. Combining the virtues of these two approximation schemes a new procedure has been proposed \[[@I2007]\] that maintains a wide range of validity in the regime of large values of both the coupling strength and the detuning of frequencies $(\lambda \sim O(\omega), \Delta \lesssim \omega )$. This generalization of the rotating wave approximation employs \[[@I2007]\], as a suitable intermediate step, a change of basis to that introduced by the adiabatic approximation \[[@IGMS2005], [@AN2010]\]. Mimicking an inherent feature of the ordinary rotating wave approximation, the generalization advanced in \[[@I2007]\] retains, in the transformed basis, the ‘energy-conserving’ one-particle transition elements in the Hamiltonian matrix which now assumes a direct sum of $2 \times 2$ diagonal blocks apart from an uncoupled ground state.

For the purpose of making our notations and subsequent applications clear, we, following \[[@I2007]\], now briefly review the diagonalization of the Hamiltonian in the generalized rotating wave approximation scheme, and explicitly list its eigenstates. The adiabatic approximation physically signifies that the rapidly moving oscillator adjusts itself quickly to the slow changing qubit observable $\sigma_{z}$, which may be considered to reside in one of its eigenstates $|\pm1\rangle$. Setting $\Delta = 0$, the bipartite Hamiltonian (\[qOH\]) is readily diagonalized in the tensored basis $\ket{\pm 1,n_{\pm}}$, where the displaced number states read: $\ket{n_{\pm}}= \mathrm{D}^{\dagger}\left(\pm \frac{\lambda}{\omega}\right) \ket{n},
 \, \mathrm{D}\left(\alpha \right) = \exp\left(\alpha a^{\dagger}- \alpha^{*}a \right),\, \alpha \in \mathbb{C}$. The eigenvalues of both the displaced basis states $\ket{\pm 1,n_{\pm}}$ are degenerate: $E_{n} = 
 \omega\big(n - \frac{x}{4}\big), \, x=4\lambda^{2}/\omega^2$. The matrix elements for the projections of the displaced number states \[[@IGMS2005]\] are listed below: \_[m,n]{} =

(-1)\^[m-n]{} \^ (-) L\_n\^[(m-n)]{}(x), & m n\
[x]{}\^ (-) L\_m\^[(n-m)]{}(x) & m &lt; n, \[Mmn\]

where the associated Laguerre polynomial reads $L_{n}^{(j)}(x) = \sum_{k = 0}^{n} 
\,(-1)^{k}\, \binom{n + j}{n - k}\,\frac{x^{k}}{k!}$. The reality of the matrix element (\[Mmn\]) leads to the identity \_[m,n]{} = (-1)\^[m+n]{} \_[n,m]{} \[IdenM\] that will be used later to simplify the density matrix elements. In the adiabatic approximation scheme, the degeneracy of the states are lifted by the qubit Hamiltonian that causes the mixing of the equi-displaced bipartite states $\ket{\pm 1,n_{\pm}}$. Diagonalizing the resulting $2 \times 2$ blocks associated with the same photon number the energy eigenvalues and the corresponding eigenstates are obtained \[[@IGMS2005]\]: E\^[()]{}\_[n(0)]{} = ( n- ) L\_[n]{}(x), = (-x/2), = ( ). \[E\_adiabatic\]

The basis obtained in (\[E\_adiabatic\]) is now utilized \[[@I2007]\] for presenting the generalized rotating wave approximation scheme. The Hamiltonian (\[qOH\]), when expressed in the adiabatic eigenstates (\[E\_adiabatic\]), shows characteristics similar to the standard rotating wave approximation. Proceeding parallely, the tridiagonal matrix elements that contain the ‘energy conserving’ one particle transitions are retained \[[@I2007]\], and the remote off-diagonal matrix elements that manifest higher order transitions are neglected within the generalized rotating wave approximation scheme. The ground state is now decoupled from all other states, whereas the remaining matrix elements are organized in a tower of $2 \times 2$ diagonal blocks \[[@I2007]\] as

E\^[(+)]{}\_[n-1]{} & \_[n]{}\
\_[n]{} & E\^[(-)]{}\_[n]{}

,\_[n]{} = \_[n-1,n]{} = L\^[(1)]{}\_[n-1]{}(x), n 1. \[2X2GRWA\] The energies of the singlet ground state, and the infinite tower of the doublets read \[[@I2007]\] \_[0]{} & & E\^[(-)]{}\_[0]{} =--,\
\^[()]{}\_[n (1)]{}& = &(n - - ) + (L\_[n-1]{} (x) - L\_[n]{}(x))\
& &. \[energy\_GRWA\] The corresponding eigenstates are explicitly given by |[E]{}\_0 = ( - ), |[E]{}\_[n (1)]{}\^[()]{} = \^[()]{}\_[n]{} | E\^[(+)]{}\_[[n-1]{}]{} \^[()]{}\_[n]{}|E\^[(-)]{}\_[n]{} , \[Hn\_eigenstate\] where we abbreviate: $\chi_n=\sqrt{\zeta_n^2+\varepsilon_{n}^2},\,\varepsilon_{n} 
=\frac{E^{(+)}_{n-1}-E^{(-)}_{n}}{2},\, \mu^{(\pm)}_{n} = 
\sqrt{\frac{\chi_n \pm \varepsilon_{n}} {2\chi_n}}$. The completeness requirement of the bipartite basis states (\[Hn\_eigenstate\]) now reads: |[E]{}\_[0]{} \_[0]{}| + \_[n = 1]{}\^ (|[E]{}\_[n]{}\^[(+)]{} \_[n]{}\^[(+)]{}| + |[E]{}\_[n]{}\^[(-)]{} \_[n]{}\^[(-)]{}|) = \_[n = 0]{}\^ (| E\_[n]{}\^[(+)]{} E\_[n]{}\^[(+)]{}| + |E\_[n]{}\^[(-)]{} E\_[n]{}\^[(-)]{}|) = \[completeness\]

With the above construction of the energy eigenstates via the generalized rotating wave approximation procedure in hand, we now study the time evolution of the entanglement of the bipartite system. The initial qubit-oscillator hybrid entangled states read \_[()]{}= ( + ) ( - ) , = (). \[t0quasiBell\] As the Schrödinger cat states $\sim\, \ket{\alpha} \pm \ket{-\alpha}$ are mutually orthogonal, the initial states (\[t0quasiBell\]) may be regarded as examples of the hybrid Bell states. Even though our subsequent analysis may be developed with both the initial states (\[t0quasiBell\]), we, for the purpose of notational simplicity, quote the results for the choice $\ket{\Psi(0)}_{(-)}$ and omit the subscript hereafter. The complete basis states (\[Hn\_eigenstate\]) immediately implement the mode expansion = \_[0]{} + \_[n=1]{}\^ \_[n]{}\^[()]{} . \[InExpansion\] The coefficients of the vector space expansion (\[InExpansion\]) are listed below: \_[0]{} &=& - + i\_ ,\_ = (),\
\_[n]{}\^[()]{} &=& - -i\_ P\^[(+)]{}\_[n-1]{} \^[()]{}\_[n]{} \^[()]{}\_[n]{} .\
& & . - - + i\_ P\^[(-)]{}\_[n-1]{} \^[()]{}\_[n]{} \^[()]{}\_[n]{} , \[C\_def\] where the following definitions are used: $\alpha_{\pm} = \alpha \pm \frac{\sqrt{x}}{2}, P_{n}^{(\pm)} = \frac{1 \pm (-1)^{n}}{2}$. The basis of approximate eigenstates (\[Hn\_eigenstate\]) of the Hamiltonian (\[qOH\]) readily yields the time evolution of the initial state (\[InExpansion\]): = \_[0]{}(t) + \_[n=1]{}\^ \_[n]{}\^[()]{}(t) , \[t\_BipState\] where the coefficients read: $\mathcal{C}_{0}(t) = 
\mathcal{C}_{0} \exp\left(-i\mathcal{E}_{0} t\right),\,
\mathcal{C}_{n}^{(\pm)}(t)=\mathcal{C}_{n}^{(\pm)} \exp\Big(-i\mathcal{E}^{(\pm)}_{n} t\Big)$. The time dependent density matrix of the bipartite pure state has the tensored form (t) = . \[t\_DenMa\]

The reduced density matrices of the individual subsystems are extracted by partial tracing over the complementary subspaces in the full Hilbert space. For instance, the tracing over the oscillator states yields the reduced density matrix of the qubit: \_[Q]{} \_ =

+ &\
\^[\*]{} & -

, \[Qubit\_DenMa\] where the matrix elements read &=& \^[\*]{} + \_[n=1]{}\^ [\_[n+1]{}(t)]{}\^[\*]{} \_[n]{}(t) ,\
&=& - |\_[0]{}(t)|\^[2]{} \_[0,0]{}- \_[n=1]{}\^ ( [\_[0]{}(t)]{}\^[\*]{} \_[n]{}(t) + (-1)\^[n]{} \_[0]{}(t) [\_[n]{}(t)]{}\^[\*]{} ) \_[0,n-1]{}\
& &+ ( \^[\*]{} \_[n]{}(t) + (-1)\^[n]{} \_[0]{}(t) [\_[n]{}(t)]{}\^[\*]{} ) \_[0,n]{} + \_[n,m=1]{}\^ \^[\*]{} \_[m]{}(t)\_[n-1,m-1]{}\
& & - [\_[n]{}(t)]{}\^[\*]{} \_[m]{}(t) \_[n,m]{} + ([\_[n]{}(t)]{}\^[\*]{} \_[m]{}(t) + (-1)\^[m+n]{} \_[n]{}(t) [\_[m]{}(t)]{}\^[\*]{}) \_[n-1,m]{} ,\
\_[n]{}(t) &=& \_[n]{}\^[(+)]{} \_[n]{}\^[(+)]{}(t) + \_[n]{}\^[(-)]{} \_[n]{}\^[(-)]{}(t),\_[n]{}(t) = \_[n]{}\^[(-)]{} \_[n]{}\^[(+)]{}(t) - \_[n]{}\^[(+)]{} \_[n]{}\^[(-)]{}(t) . \[2X2\_DenMaElts\] In the above expression we have used the identity (\[IdenM\]). The expectation value of the qubit spin variable $\langle\sigma_{z}\rangle$ measures the statistical average of the population inversion. For the density matrix (\[Qubit\_DenMa\]) it is given by \_[z]{}(\_[z]{}\_) = 2 . \[PopInv\] In general, the pair of eigenvalues of the qubit density matrix (\[Qubit\_DenMa\]) ,= \[rho\_eigenvalue\] allows us to compute its von Neumann entropy $S(\rho_{\cal Q}) \equiv - \hbox {Tr} (\rho_{\cal Q} \, \log \rho_{\cal Q})$ as S\_[Q]{} = - ( + ) ( + ) - ( - ) ( - ). \[entropy\_defn\]

The bipartite density matrix (\[t\_DenMa\]) also similarly produces the reduced density matrix of the oscillator via partial tracing on the qubit Hilbert space: $\rho_{\cal{O}} \equiv \hbox{Tr}_{\cal Q} \, \rho$. Its explicit construction is listed below: \_(t) &=& |\_[0]{}(t)|\^[2]{} P\_[0,0]{}\^[(+)]{} + \_[n=1]{}\^ ( \_[0]{}(t) [\_[n]{}(t)]{}\^[\*]{} P\_[0,n-1]{}\^[(-)]{} + \^[\*]{} \_[n]{}(t) P\_[n-1,0]{}\^[(-)]{} + \_[0]{}(t) [\_[n]{}(t)]{}\^[\*]{} P\_[0,n]{}\^[(+)]{}\
& &+ [\_[0]{}(t)]{}\^[\*]{} \_[n]{}(t) P\_[n,0]{}\^[(+)]{} ) + \_[n,m=1]{}\^ ( \_[n]{}(t) [\_[m]{}(t)]{}\^[\*]{} P\_[n-1,m-1]{}\^[(+)]{} + \_[n]{}(t) [\_[m]{}(t)]{}\^[\*]{} P\_[n,m]{}\^[(+)]{}\
& &+ \_[n]{}(t) [\_[m]{}(t)]{}\^[\*]{} P\_[n,m-1]{}\^[(-)]{} + \_[n]{}(t) [\_[m]{}(t)]{}\^[\*]{} P\_[n-1,m]{}\^[(-)]{} ), \[density\_oscillator\] where the symmetrized projection operators read $P_{n,m}^{(\pm)}= \frac{1}{2} \left( \ket{n_{+}}\bra{m_{+}} \pm \ket{n_{-}}\bra{m_{-}} \right),\,(n, m = 0, 1, \ldots)$. The density matrix (\[density\_oscillator\]) obeys the normalization condition: $\mathrm{Tr}\,\rho_{\mathcal{O}}(t) = 1$. For later use we quote the initial oscillator density matrix: $\rho_{\mathcal{O}}(0) =\frac{1}{2}\,(\ket{\alpha}\bra{\alpha} + \ket{-\alpha}\bra{-\alpha})$.

To obtain measures of the entanglement and the mixedness of the bipartite system the von Neumann entropy of the reduced density matrix of the qubit may be considered. It is well-known \[[@AL1970]\] that if a composite system, comprising of two subsystems, resides in a pure state, the entropies of both subsystems are equal. In the present example it holds for the oscillator with an infinite dimensional Hilbert space, and the two-level qubit interacting with it: $S_{\cal Q} = S_{\cal O} \equiv S$.

Quasi-probability distributions on the oscillator phase space
=============================================================

The diagonal Sudarshan-Glauber $P$-representation {#Sudarshan}
-------------------------------------------------

It is well-known that the coherent state representation is overcomplete \[[@Schleich2001]\]. Employing the overcompleteness property, it has been established \[[@S1963], [@G1963]\] that an arbitrary oscillator density matrix $\rho_{\cal O}$ may be realized via diagonal coherent state projectors: \_[O]{} = P(,\^[\*]{}) \^[2]{}. \[P\_def\] The normalization of the density matrix $\rho_{\cal O}$ ensures the property P(, \^[\*]{}) \^[2]{}= 1, \[P\_norm\] whereas its Hermiticity implies $P(\beta, \beta^{*})$ to be real. The $P$-representation, however, is a general distribution of indeterminate sign as it possibly includes highly singular derivatives of the $\delta$-function \[[@S1963]\]. For the coherent state, which corresponds as closely as possible to a classical harmonic oscillator state with a complex amplitude, the $P$-representation is a positive definite $\delta$-function \[[@Schleich2001]\] which is a legitimate measure of the classical probability density. For this reason the coherent state is considered to be classical. On the other hand, if the diagonal kernel $P(\beta,\beta^{*})$ can not be considered a valid probability measure due to its negativity, the oscillator state given by the density matrix $\rho_{\cal O}$ exhibits nonclassical features. The diagonal $P$-representation is known to be unique, even though it may be realized in different equivalent forms \[[@M2009]\]. Towards inverting the diagonal representation (\[P\_def\]) the Fourier transform is employed \[[@Schleich2001]\]: P(,\^[\*]{})= \_ (||\^[2]{}) (\^[\*]{} - \^[\*]{} ) \^[2]{}. \[P\_evaluation\] The invertibility of the relations (\[P\_def\], \[P\_evaluation\]) suggests that both the density matrix $\rho_{\cal{O}}$ and the diagonal representation $P(\beta,\beta^{*})$ encapsulate equivalent informations. We now explicitly evaluate the pseudo probability distribution corresponding to the oscillator density matrix (\[density\_oscillator\]): P(,\^[\*]{}) &=& |\_[0]{}(t)|\^[2]{} D\^[(+)]{}\_[0,0]{} + \_[0]{}(t)\^[\*]{} \_[n=1]{}\^ ( \_[n]{}(t) D\^[(-)]{}\_[0,n-1]{} + \_[n]{}(t) D\^[(+)]{}\_[0,n]{} )\
& &+ \_[n,m=1]{}\^ \_[n]{}(t)\^[\*]{} \_[m]{}(t) D\^[(-)]{}\_[n,m-1]{} + \_[n,m=1]{}\^ ( \_[n]{}(t)\^[\*]{} \_[m]{}(t) D\^[(+)]{}\_[n-1,m-1]{}\
& &+ \_[n]{}(t)\^[\*]{} \_[m]{}(t) D\^[(+)]{}\_[n,m]{} ), \[P\_function\] where the arbitrarily singular derivatives of the $\delta$-functions read D\_[n,m]{}\^[()]{} &=& (|\_[+]{}|\^[2]{}) ( -)\^[n]{} ( -)\^[m]{} \^[(2)]{}(\_[+]{})\
& &(|\_[-]{}|\^[2]{}) ( -)\^[n]{} ( -)\^[m]{} \^[(2)]{}(\_[-]{}) , \_=. \[D\_singular\] The diagonal $P$-representation (\[P\_function\]) satisfies the normalization condition (\[P\_norm\]), and is obviously real. At $t=0$, the initial hybrid entangled state (\[t0quasiBell\]) gives rise to a positive definite $P$-representation: $P(\beta,\beta^{*})|_{{t=0}} = 2^{-1}\, \exp(-(|\alpha|^{2}-|\beta|^{2}))\,
(\delta^{(2)}(\alpha-\beta) + \delta^{(2)}(\alpha+\beta))$. At an arbitrary time, however, a large number of interaction-dependent modes develop, and produce quantum interferences among themselves. The quantum properties in (\[P\_function\]) are evident due to the appearance of rapidly oscillating singular derivatives of $\delta$-function in the $P$-representation, which at time $t > 0$ is clearly not a nonnegative definite quantity. The existence of negativity of the diagonal $P$-representation (\[P\_function\]) points towards the manifestations of specific quantum effects such as anti-bunching \[[@M1982]\] and quadrature squeezing \[[@Yurke1985]\]. Other phase space pseudo probability functions such as Wigner $W$-distribution and the Husimi $Q$-function may be obtained from the diagonal $P$-representation after its suitable smoothing by the Gaussian kernels.

The Wigner $W$-distribution {#Wigner}
---------------------------

The Wigner phase space distribution for a quantum state given by the oscillator density matrix (\[density\_oscillator\]) is described \[[@Schleich2001]\] as W(, \^[\*]{}) = (()\_) ( \^[\*]{} - \^[\*]{} ) \^[2]{}, \[W\_def\] which is a real-valued function maintaining the normalization property: W(,\^[\*]{}) \^[2]{}= 1. \[W\_normalization\] A characteristic of the Wigner distribution that distinguishes it from a classical probability density, is that its integral over a given subregion of the phase space may be negative or greater than one. As an efficient computational algorithm the $W$-distribution has been expressed \[[@MK1993]\] as a sum of an infinite series of the expectation values of the oscillator density matrix $\rho^{\phantom{A}}_{\cal{O}}$ in the displaced number basis states: W (, \^[\*]{}) = \_[k=0]{}\^ (-1)\^[k]{} \^\_ , = () . \[W\_series\] On the other hand, the convolution of a smoothing Gaussian function of variance $1/2$ on the phase space with the diagonal $P$-representation leads to the $W$-distribution that is free from the divergences brought forward by the derivatives of the $\delta$-function contained in the expression (\[P\_function\]) of the $P$-representation: W(,\^[\*]{}) = P(, \^[\*]{}) (-2|-|\^[2]{}) \^[2]{}. \[P\_int\_W\] Employing the explicit evaluation (\[P\_function\]) of the $P$-representation and the convolution integral (\[P\_int\_W\]) we now provide a direct derivation of the $W$-distribution. Towards this we employ the identity $$\begin{aligned}
\int &\exp \Big( -2|\beta - \gamma|^{2}+|\gamma |^{2}\Big) 
\Big( -\dfrac{\partial}{\partial \gamma^{*}} \Big)^{n} 
\Big( -\dfrac{\partial}{\partial \gamma} \Big)^{m} \delta^{(2)}(\gamma)\, \mathrm{d}^{2}\gamma\nn\\
&= 2^{n+m} \beta^{n} \beta^{*m}  \exp(-2|\beta|^{2})\,{}_2F_0\Big( \! \! -n,-m;\phantom{}_{-} ;
-\frac{1}{4|\beta|^{2}}\Big),
\label{Int_hyper}\end{aligned}$$ where the hypergeometric sum reads $ \displaystyle{ {}_2F_0(\mathsf{x},\mathsf{y};\phantom{}_{-}; \tau)
= \sum_{k=0}^{\infty} (\mathsf {x})_{k} (\mathsf{y})_{k}\, \frac{\tau^{k}}{k!},\,
(\mathsf {x})_{k} = \prod_{\ell = 0}^{k - 1}(\mathsf {x} + \ell) }$. The above hypergeometric function with the negative integral numerator may be expressed \[[@M1939]\] via the Charlier polynomial: $\mathrm{c}_{k} (\ell; \tau) = {}_{2}F_0\left(-k,-\ell;\phantom{}_{-} ;
-\frac{1}{\tau}\right)\,\forall \tau > 0$. To emphasize the equal footing of the integers $(k, \ell)$ in (\[Int\_hyper\]) we, however, do not use this notation here. The integral representation (\[P\_int\_W\]) aided by the identity (\[Int\_hyper\]) now generate the $W$-distribution: W(,\^[\*]{}) &=& |\_[0]{}(t)|\^[2]{} \_[0,0]{}\^[(+)]{} + 2 \_[0]{}(t)\^[\*]{} \_[n=1]{}\^ ( \_[n]{}(t) \_[0,n-1]{}\^[(-)]{}(,\^[\*]{}) + \_[n]{}(t) \_[0,n]{}\^[(+)]{}(,\^[\*]{}) )\
& +& \_[n,m=1]{}\^ \_[n]{}(t)\^[\*]{} \_[m]{}(t) \_[n,m-1]{}\^[(-)]{} (,\^[\*]{}) + \_[n,m=1]{}\^ ( \_[n]{}(t)\^[\*]{} \_[m]{}(t) \_[n-1,m-1]{}\^[(+)]{} (,\^[\*]{})\
&+& \_[n]{}(t)\^[\*]{} \_[m]{}(t) \_[n,m]{}\^[(+)]{} (,\^[\*]{}) ), \[wigner\_function\] where the weights involving the Gaussian functions are given by \^[()]{}\_[n,m]{}(, \^[\*]{}) &=& \_[+]{}\^[n]{} \_[+]{}\^[\*\^m]{} (-2|\_[+]{}|\^[2]{}) \_2F\_0( -n,-m;\_[-]{} ; -)\
& &\_[-]{}\^[n]{} \_[-]{}\^[\*\^m]{} (-2|\_[-]{}|\^[2]{}) \_2F\_0( -n,-m;\_[-]{} ; -) . \[HW\]

An alternate derivation of the Wigner $W$-distribution follows from sum rule (\[W\_series\]). Towards this we list the matrix elements = (-1)\^[k]{} (- i \_) \_2F\_0( -n,-k;\_[-]{}; -) \[NAlphaK\] that facilitate the evaluation of the series (\[W\_series\]). Another necessary tool for the present derivation is the identity $$\begin{aligned}
\sum_{k=0}^{\infty} \frac{(-1)^{k} \tau^{k}}{k!}\, {}_2F_0\Big( -n,-k;\phantom{}_{-} ; -\frac{1}{\tau}\Big) &
\,{}_2F_0\Big( -k,-m;\phantom{}_{-} ; -\dfrac{1}{\tau}\Big) = \nn \\
& 2^{n+m} \exp(-\tau) \; {}_2F_0\Big( -n,-m;\phantom{}_{-} ; -\frac{1}{4\tau}\Big)
\label{IdentityHyper}\end{aligned}$$ that readily follows \[[@M1939]\] from the bilinear generating function of the Charlier polynomials. The series sum (\[W\_series\]), in conjunction with the results (\[NAlphaK\], \[IdentityHyper\]), now reproduce the $W$-distribution obtained in (\[wigner\_function\]). This serves as a consistency check on our derivations (\[P\_function\], \[wigner\_function\]) of the phase space distributions in the present model.

The negativity of the Wigner function is an indication towards the existence of nonclassical properties of the states. At $t = 0$ the $W$-distribution for our initial hybrid entangled state (\[t0quasiBell\]) is positive definite: $W(\beta,\beta^{*})|_{t=0} = \pi^{-1}\,
(\exp(-2|\alpha - \beta|^{2}) + \exp(-2|\alpha + \beta|^{2}))$. At a later time, however, interactions produce a large number of modes that interfere with each other. Interference fringes gives rise to appearance of phase space domains representing negative values of the $W$-distribution. A fruitful measure of negativity relate to the volume of the negative part of the Wigner function on the phase space \[[@KZ2004]\]: \_[W]{} = |W()| -1. \[negativity\] We will return to the topic in Sec. \[DelocalQuantum\].

The Husimi $Q$-function {#Husimi}
-----------------------

The Husimi $Q$-function \[[@Schleich2001]\] is a quasi probability distribution defined as expectation value of the oscillator density matrix in an arbitrary coherent state. It assumes nonnegative values on the phase space in contrast to the other phase space quasi probabilities. Being easily computable it is extensively employed \[[@SA2002], [@IWAH2002]\] in the study of the occupation on the phase space. For our reduced density matrix of the oscillator (\[density\_oscillator\]) the $Q$-function Q(,\^[\*]{}) = \_[O]{}\[Q\_defn\] maintains the normalization restriction: $\int Q(\beta, \beta^{*}) \mathrm{d}^{2}\beta = 1$ and the bounds: $0 \leq Q(\beta,\beta^{*}) \leq \frac{1}{\pi}$. Our construction of the oscillator density matrix (\[density\_oscillator\]) now yields the time-evolution of the $Q$-function: Q(,\^[\*]{}) &=& |\_[0]{}(t)|\^[2]{} H\^[(+)]{}\_[0,0]{}(,\^[\*]{}) + \_[0]{}(t)\^[\*]{} \_[n=1]{}\^ ( \_[n]{}(t) H\^[(-)]{}\_[0,n-1]{}(,\^[\*]{})\
&+& \_[n]{}(t) H\^[(+)]{}\_[0,n]{}(,\^[\*]{}) ) + \_[n,m=1]{}\^ \_[n]{}(t)\^[\*]{} \_[m]{}(t) H\^[(-)]{}\_[n-1,m]{}(,\^[\*]{})\
&+& \_[n,m=1]{}\^ ( \_[n]{}(t)\^[\*]{} \_[m]{}(t) H\^[(+)]{}\_[n-1,m-1]{}(,\^[\*]{}) +\_[n]{}(t)\^[\*]{} \_[m]{}(t) H\^[(+)]{}\_[n,m]{}(,\^[\*]{}) ), \[Q\_function\] where the weight functions on the phase space read H\^[()]{}\_[n,m]{} (, \^[\^\*]{})= ( \_[+]{}\^[n]{} \_[+]{}\^[\*\^m]{} (-|\_[+]{}|\^[2]{}) \_[-]{}\^[n]{} \_[-]{}\^[\*\^m]{} (-|\_[-]{}|\^[2]{}) ). \[Q\_weight\] A convolution integral of the diagonal $P$-representation with a Gaussian weight of unit variance produce \[[@Schleich2001]\] the well-behaved nonnegative Husimi $Q$-function as follows: Q(,\^[\*]{}) = P(,\^[\*]{}) (-|-|\^[2]{}) \^[2]{}. \[P\_int\_Q\] The explicit constructions of the $P$-representation (\[P\_function\]) and that of the $Q$-function (\[Q\_function\]) obey the consistency check (\[P\_int\_Q\]). Another convolution property that express the $Q$-function via the $W$-distribution reads \[[@Schleich2001]\]: Q(,\^[\*]{})=W(,\^[\*]{}) (-2|-|\^[2]{}) \^[2]{}. \[Q\_int\_W\] Employing the integral representation (\[Int\_hyper\]) it follows that the phase space quasi probability densities (\[wigner\_function\], \[Q\_function\]) maintain the integral sum rule (\[Q\_int\_W\]). The convolution relations (\[P\_int\_W\], \[Q\_int\_W\]) suggest an interesting feature. The diagonal $P$-representation containing singular derivatives of $\delta$-functions is smoothed at two stages by a Gaussian weight of variance $1/2$. After the first smoothing operation the $W$-distribution, which is nonsingular but not non-negative, appears. Subsequent smoothing operation produces a nonsingular *and* nonnegative Husimi $Q$-function. If a single smoothing operation (\[P\_int\_Q\]) is performed on $P$-representation via the Gaussian weight of variance $1$, the $Q$-function is directly reproduced. We symbolically summarize this as $[[P]_{1/2}]_{1/2} = [P]_{1}$, where the notation reads $[\mathcal{X}]_{s} =  
\frac{1}{\pi\, s}\int \mathcal{X}(\gamma,\gamma^{*})  
\exp\left(-\frac{|\beta-\gamma|^{2}}{s}\right)\,\mathrm{d}^{2}\gamma$.

One of the utilities of the $Q$-function is that it provides a convenient evaluation of the expectation values of the operators expressed in their antinormal ordered form \[[@Schleich2001]\]. The first and second moments of the quadrature variable defined as X\_= (a (-i)+ a\^ (i)) \[angle\_hermitian\_operator\] may be expressed \[[@Schleich2001]\] via the following phase space integrals over the complex plane: X\_ && (X\_ \_[O]{}(t)) = (e\^[-i]{}+ \^\* e\^[i]{}) Q(, \^[\*]{}) d\^2 ,\
X\_\^2 && (X\_\^2 \_[O]{}(t)) = ( (e\^[-i]{}+ \^\* e\^[i]{})\^2-1 ) Q(, \^[\*]{}) d\^2 . \[moment\_def\] The $Q$-function (\[Q\_function\]) now allows explicit determination of the above quadrature moments: X\_ &=& (\_[1]{}(t)(-i)) - , \[1stMoment\]\
X\_\^2 &=& (\_[2]{}(t)(-i2) - \_[1]{}(t) (1+(-i2)))\
& & + ( \_[1]{}(t) - ) + \^[2]{}\[2ndMoment\] where the time-dependent coefficients $\{\mathbf{G}_{k}(t), \mathbf{N}_{k}(t), \mathbf{F}_{k}(t)|k\in (1, 2)\}$ expressed as sums over Fourier modes read: \_[k]{}(t) &=& \_[0]{}(t)\^[\*]{} \_[k]{}(t) + \_[n=1]{}\^ \_[n]{}(t)\^[\*]{} \_[n+k]{}(t)\
& & + \_[n]{}(t)\^[\*]{} \_[n+k]{}(t) ,\
\_[k]{}(t) &=& |\_[0]{}(t)|\^[2]{}+ \_[n=1]{}\^ n\^[k]{} \_[n]{}(t)\^[\*]{} \_[n]{}(t) + (n+1)\^[k]{} \_[n]{}(t)\^[\*]{} \_[n]{}(t) ,\
\_[k]{}(t) &=& k \_[0]{}(t)\^[\*]{} \_[2]{}(t) + \_[n=1]{}\^ (n+1)\_[k-1]{} \_[n]{}(t)\^[\*]{} \_[n]{}(t)\
& & + (n+2)\_[k-1]{} \_[n]{}(t)\^[\*]{} \_[n+2]{}(t) . \[time\_dependent\_coeff\] The variance of the quadrature variable is given by $V_{\theta}= \langle X_{\theta}^{2} \rangle - \langle X_{\theta} \rangle ^2,$ which we will later employ for studying emergence of the squeezed states during the time evolution. For later utilization we also quote here the mean photon number $\langle\hat{n}\rangle
 = \hbox{Tr} (a^{\dagger} a\,\rho_{\cal O}(t))$ and its variance $\langle (\Delta\hat{n})^{2}\rangle = 
 \hbox{Tr} (\hat{n}^{2}\,\rho_{\cal O}(t)) - \langle\hat{n}\rangle^{2}$: &=& \_[1]{}(t) + - \_[1]{}(t) -1 \[avg\_photon\_number\]\
&=& \_[2]{}(t) + \_[1]{}(t) + (\_[2]{}(t)) - + (\_[1]{}(t)) - 2 (\_[2]{}(t))\
&&- \_[1]{}(t) - (\_[1]{}(t)) \^[2]{}. \[variance\_num\]

Another dynamical quantity that is useful in the study of ‘kitten’-like states is the polar phase density of the Husimi $Q$-function \[[@TMG1993]\] obtained via its radial integration on the phase space: ()=\_[0]{}\^ Q(, \^[\*]{}) ||||, =|| (i), \[QPD\] which is a convenient tool for describing the splitting of the $Q$-function. Towards obtaining a series expansion of the phase density $\mathcal{Q}(\theta)$ we define the weight factors depending on the polar angle: \^[()]{}\_[n,m]{} ()= \_[0]{}\^ H\^[()]{}\_[n,m]{} (, \^[\^\*]{}) |||| \[H\_theta\] that admit, via (\[Q\_weight\]), explicit evaluation as follows: \^[()]{}\_[n,m]{} () &=& \_[= 0]{}\^[n]{}\_[=0]{}\^[m]{} (i (- ) ) (-) ()\^\
&&P\^[()]{}\_[n+m--]{} ( + 1) \_1F\_1 (- ; ; -)\
&&-P\^[()]{}\_[n+m--]{} () \_1F\_1 (-; ; -), \[H\_nm\] where we have used the integration given by \_[0]{}\^ r\^p (-(r )\^2) dr &=& ( ) \_1F\_1 (- ; ; - \^[2]{})\
&& ( + 1) \_1F\_1 ( ; ; -\^[2]{}). \[r\_shift\_int\] The above compendium of expressions (\[QPD\], \[Q\_function\], \[H\_nm\]) now provides the evolution of the polar phase density $\mathcal{Q}(\theta)$ for our initial state (\[t0quasiBell\]) () &=& |\_[0]{}(t)|\^[2]{} \^[(+)]{}\_[0,0]{}() + \_[0]{}(t)\^[\*]{} \_[n=1]{}\^ ( \_[n]{}(t) \^[(-)]{}\_[0,n-1]{}()\
&&+ \_[n]{}(t) \^[(+)]{}\_[0,n]{}() ) + \_[n,m=1]{}\^ \_[n]{}(t)\^[\*]{} \_[m]{}(t) \^[(-)]{}\_[n-1,m]{}()\
&& + \_[n,m=1]{}\^ ( \_[n]{}(t)\^[\*]{} \_[m]{}(t) \^[(+)]{}\_[n-1,m-1]{}() +\_[n]{}(t)\^[\*]{} \_[m]{}(t) \^[(+)]{}\_[n,m]{}() ). \[Q\_PD\]

Delocalization on the phase space {#DelocalQuantum}
=================================

Wehrl entropy {#Wehrl}
-------------

The Wehrl entropy defined as \[[@W1978]\] S\_[Q]{}=-Q(, \^[\*]{}) Q(, \^[\*]{}) \^[2]{}\[WehrlDef\] is an information-theoretic measure estimating the delocalization of the system in the oscillator phase space. It is considered \[[@BKK1995]\] as a count of an equivalent number of widely separated coherent states necessary for covering the existing phase space occupation of the coupled oscillator. Being subject to the restriction originating from the Heisenberg uncertainty principle, the Wehrl entropy (\[WehrlDef\]) is a positive definite quantity \[[@W1978]\]. In the present case we employ the definition (\[WehrlDef\]) and the evolution (\[Q\_function\]) of the $Q$-function to numerically study the long-range time dependence of $S_{Q}$ for various values of the coupling strength. We note that here and hereafter the time is measured in the natural unit: $\omega^{-1}$. We observe the following properties: () In the long time limit $t \gtrsim (x^{2} \widetilde{\Delta})^{-1}$ the quasi periodicity of the Wehrl entropy is manifest in the coupling strength regime $\lambda \lesssim 0.05 \,\omega$, where the Laguerre polynomials $L_{n}^{(j)}(x)$ are well-approximated by their quadratic components $O(x^{2})$. Frequency modes $O(x^{2} \widetilde{\Delta})$ and their harmonics produced via the interaction now give rise to the quasi periodicity of $S_{Q}$, where the long range time period maintains the property: $T_{\mathrm{long}} \propto  
\lambda^{-4}\, \exp (x/2)$ (Figs. \[Wehrl\_Time\] (a)-(c)). For a smaller value of the qubit frequency $\Delta$ the quasi periodic behavior persists for a comparatively higher coupling strength $\lambda$. This follows from the requirement that for the quasi periodicity to hold, the phase change caused by the higher order fluctuations $\{O(x^{n})| n > 2\}$ during the time span $T_{\mathrm{long}}$ is to be negligibly small: $T_{\mathrm{long}}\, x^{3} \widetilde{\Delta} \ll 1$. The time period observed in Figs. \[Wehrl\_Time\](a)-(c) are noted below: $T_{\mathrm{long}} = 6.9914 \times 10^{6}$ (for $\lambda = 0.008\,\omega$), $T_{\mathrm{long}} = 2.9568 \times 10^{6}$ (for $\lambda = 0.010\, \omega$), $T_{\mathrm{long}} = 
 1.4743 \times 10^{6}$ (for $\lambda = 0.012\, \omega$), respectively. The near equality of the product $T_{\mathrm{long}}\, \lambda^{4}\, \exp (-x/2)$ in the respective cases ($0.02863, 0.02956, 0.03056)$ validates our argument that the quantum fluctuations $O(x^{2})$ produce the observed long range time period. In the instance of nonlinear Kerr-like medium similar behavior in the time evolution of $S_{Q}$ was previously noticed \[[@JO1994]\], where its local minima corresponded with the formations of finite superposition of coherent states. In the present bipartite interacting model, however, the qubit-oscillator interaction superimposes short time span fluctuations of the frequency $O(x \widetilde{\Delta})$ on the long range oscillations that now act as an envelope of the total time evolution of the Wehrl entropy. This introduces important distinctions to the present model. We will discuss this in the Subsec. \[kitten\]. () In the ultra-strong coupling regime $\lambda \gtrsim  0.1\, \omega$ all frequency modes $\{O(x^{n} \widetilde{\Delta})| n= 0, 1, \ldots\}$ and their harmonics arise. Random phase differences between a large number of incommensurate modes cause the resultant interference to average out, while ensuring an effective stabilization of the occupation of the phase space after an initial build up (Fig. \[Wehrl\_Time\] (d)). The rapid high frequency $O(\omega)$ fluctuations are of small amplitude: $|\Delta S_{Q}|/ S_{Q} \ll 1$, and may be removed by a suitable coarse graining process \[[@JJ2007]\]. In this regime it is observed that for a fixed initial state parameter $\alpha$ the time-averaged value of the Wehrl entropy gradually increases with increasing coupling strength. Following (\[avg\_photon\_number\]) it is evident that higher coupling strength leads to an enhancement of the photon expectation value $\braket{\hat{n}}$ that causes a wider spread of the $Q$-function resulting in an increment in $S_{Q}$.

() Sufficiently localized states are transiently observed at large macroscopic values of the coherent state amplitude $\alpha$ even in the ultra-strong coupling domain. These states are realized (Fig.\[single\_peak\]) at the local minimum of the Wehrl entropy $S_{Q}$. Physically, the energy exchange between the qubit and the oscillator degrees of freedom induces the revival and collapse of the qubit density matrix elements. A revival of the qubit matrix element indicates that the oscillator has, reciprocally, less energy available to it. This acts as a constraint on its delocalization on the phase space. Moreover, a large value of the amplitude $\alpha$ signifies comparatively higher magnitude of energy residing in the interference pattern (Fig.\[single\_peak\] (b)) of the quantum oscillations. Qualitatively, when the energy associated with the stochastic randomization of the modes is less than the energy of the coherent interference pattern, a localization on the phase space takes place and, consequently, a relative decrease in $S_{Q}$ develops. The quantum interference pattern (Fig.\[single\_peak\] (b)) has the shape of a localized peak with ‘twisted arms’, where oscillations develop perpendicular to these arms causing a transport of energy necessary for the localization. The phase density diagram (Fig. \[single\_peak\] (d)) shows that the initial state (\[t0quasiBell\]) consisting of two almost maximally mixed Gaussian peaks coalesce at the local minimum of $S_{Q}$ to produce a partially pure transient state of the oscillator with its von Neumann entropy given by $S = 0.5883$. () The fast initial rise of $S_{Q}$ in the ultra-strong coupling domain $\lambda \sim \omega$ may be understood as follows. Increased qubit-oscillator coupling leads to the generation of all high-frequency quantum fluctuation modes. The phase randomization of the modes of incommensurate frequencies results in the rapid initial spreading on the phase space, and consequent fast production of $S_{Q}$. Once the modes are statistically populated, the Wehrl entropy $S_{Q}$, except for high frequency $O(\omega)$ quantum fluctuations of relatively small amplitude, maintains an almost stationary value. The initial production time of the Wehrl entropy $T_{\hbox{\tiny{ent. prod.}}}$ follows from the asymptotic behavior of associated Laguerre polynomials $L_n^{(j)}(x)$ at large $n \gg 1$, fixed $j$, and $x > 0$ \[[@S1975]\]: L\_n\^[(j)]{}(x) = (2 - (j +) ) + O(n\^[-]{}). \[Lagurre\_asymptotic\] The asymptotic limit of the energy eigenvalues (\[Hn\_eigenstate\]) is now readily obtained: \^[()]{}\_[n (1)]{} = (n 1 - - ) (n x)\^[-]{} (2 - )+ O(n\^[-]{}), \[E\_asymptotic\] where the leading interaction-generated part $O\Big(n^{-\frac{1}{4}}\Big)$ on the rhs provides the effective statistical stabilization of the occupation in the phase space. Recognizing this, the coupling strength dependence of the typical time scale for the generation of $S_{Q}$ is now given by $T_{\mathrm{ent. prod.}}\, \frac{\Delta}{\sqrt{\pi}}\,(n x)^{-\frac{1}{4}} \sim 1 \Rightarrow 
 T_{\mathrm{ent. prod.}} \propto \sqrt{\lambda}$. The rapid initial increase of the Wehrl entropy is described in Fig. \[Wehrl\_initial rise\]. For the data presented in Fig. \[Wehrl\_initial rise\] the proportionality constant $T_{\mathrm{ent. prod.}} /\sqrt{\lambda}$ read $30.99, 29.66, 33.11$ for the coupling strengths $\lambda=0.9\, \omega, 1.1\, \omega, 1.3\, \omega$, respectively. The discrepancy ($\sim 6\%$) in the observed data occurs since the local fluctuations in the time evolution of $S_{Q}$ play an important role in determining $T_{\mathrm{ent. prod.}}$. A suitable coarse-graining process \[[@JJ2007]\] to smooth the high frequency fluctuations may be adopted for fuller agreement.

![The initial rise of the Wehrl entropy $S_{Q}$ for the parameters $\Delta=0.5\, \omega , \alpha =2.5$, and the coupling strength $\lambda = 0.9 \,\omega\hbox{(blue)}, 1.1\, \omega\hbox{(green)}, 1.3 \, \omega \hbox{(red)}$. The horizontal dotted lines (for the respective colors) measure the corresponding time-averaged values of $S_{Q}$ that read $4.0420, 4.1280, 4.1890$ in the said order. The averaging is done for a time interval $\omega t=4000$ in each case. The vertical dotted lines signify the first crossing of the averaged $S_{Q}$ by the corresponding time evolution graph. This is regarded as the measure of $  T_{\mathrm{ent. prod.}}$. The arbitrariness in this determination can be improved by a local smoothing operation. For the respective coupling strengths these estimates read $29.40, 31.11, 37.76$ in the said order.[]{data-label="Wehrl_initial rise"}](Final-entropy-initial-rise-vertical-line-ws-del-0-5-alpha-2-5-tevolve.pdf){width="6cm" height="4.5cm"}

Kitten states and multiple time scales {#kitten}
--------------------------------------

For the Kerr-type nonlinear self-interacting photonic models, the local minima in the time evolution of $S_{Q}$ of an initial coherent state are associated \[[@JO1994]\] with transient formations of the superposition of a finite number of coherent states \[[@MTK1990]-[@MBWI2001]\] maintaining a uniform angular separation on the complex plane. These superpositions are realized at rational submultiples of the time period of $S_{Q}$. Recently such nonclassical superposition of multiple coherent states in a Kerr medium has been experimentally achieved \[[@Kirchmair2013]\]. Formation of cat-like states in a finite dimensional bosonic system that admits applying a displacement operator on its ground state has also been studied \[[@MPPN2014]\] in a Kerr medium.

In our model we study the emergence of these transitory ‘kitten’ states using the $W$-distribution and its smoothed analog the $Q$-function. In the strong coupling limit $\lambda \lesssim 0.05\, \omega$, and at specific times given by the rational submultiples of $T_{\mathrm{long}}: \{T_{p, q} = (p/q)\,T_{\mathrm{long}}|(p,q) = 1,
p \le q\}$ density matrices comprising of a finite number of macroscopic coherent states with uniform angular separation on the phase space are observed (Fig. \[Kitten\_rational\]). Starting with the initial hybrid Bell state (\[t0quasiBell\]) of the composite system, the evolution of $S_{Q}$ in the long range quasi periodic regime shows (Fig. \[Wehrl\_Time\] (b)) the existence of the local minima at rational submultiples of $T_{\mathrm{long}}$. The presence of *numerous* time scales due to the qubit-oscillator interaction in the present model, however, introduces another novel interference related feature. In particular, the interaction-dependent *linear* mode with frequency $O(x \widetilde{\Delta})$ causes an energy transfer, in a *short* time scale, between the qubit and the oscillator degrees of freedom. In the vicinity of the said times $T_{p, q}$ the oscillations of the period $T_{\mathrm{long}}$ produce a locally minimum occupation on the phase space, whereas the short time period fluctuations engineer the spread of the occupation by splitting of the Gaussian peaks. This manifests as an ordered bifurcation (*à la* Figs. \[Kitten\_rational\](a$_{2}$, a$_{4}$) and (a$_{3}$, a$_{5}$), say) of the $q$ quasi-probability peaks to $2 q$ peaks representing mixed state oscillator density matrices, while evolving from the local minima to the maxima of the *short* time period oscillations of the Wehrl entropy. These *local* variations of $S_{Q}$ in the neighborhood of $(p/q)\,T_{\mathrm{long}}$ are given in Figs. \[Kitten\_rational\] (a$_{1}$)-(c$_{1}$), where we fix $p =1; q=1, 2, 3$, respectively. Moreover, the short range time period ($T_{\mathrm{short}}$) associated with the splitting and subsequent rejoining of the peaks at a particular rational fraction $(p/q)\,T_{\mathrm{long}}$ scales inversely with $q$. For instance, from the Figs. \[Kitten\_rational\] (a$_{1}$)-(a$_{3}$), respectively, we observe that $T_{\mathrm{short}} (T_{\mathrm{long}}) \approx 6750, 
 T_{\mathrm{short}} (T_{\mathrm{long}}/2) \approx 3375, T_{\mathrm{short}} (T_{\mathrm{long}}/3) \approx 2250$. This suggests the scaling relation $T_{\mathrm{short}} (T_{\mathrm{long}}/q) \approx 
 (1/q) \,T_{\mathrm{short}} (T_{\mathrm{long}})$. In other words, due to the complex nature of the qubit-oscillator interaction in the strong coupling regime $\lambda \sim 0.05\, \omega$ the exchange of energy is realized between multiple interaction-dependent modes, and effectively the high frequency quantum oscillation $O(x \widetilde{\Delta})$ is *frequency modulated* by the low frequency component $O(x^{2} \widetilde{\Delta})$. It is worth mentioning that the oscillator at the local minima of the Wehrl entropy (Figs. \[Kitten\_rational\](a$_{2}$, a$_{4}$),(b$_{2}$, b$_{4}$), (c$_{2}$, c$_{4}$)) is close to pure states. Their respective von Neumann entropy $S$ read $0.26616, 0.17983, 0.53918$, whereas the corresponding maxima (Figs. \[Kitten\_rational\](a$_{3}$, a$_{5}$),(b$_{3}$, b$_{5}$), (c$_{3}$, c$_{5}$)) describe almost maximally mixed states.

Lastly, we note the geometry of the domain on the phase space that supports the $W$-distribution (Figs. \[Kitten\_rational\] (a$_{2}$)-(c$_{2}$), (a$_{3}$)-(c$_{3}$)). The interference pattern realized between two Gaussian peaks occurs in the intermediate phase space giving rise to oscillations in a direction perpendicular to the line joining the peaks. Alternate lines with the positive and the negative values of the $W$-distribution appear with relative phase differences of $\pi$. As the number of peaks increase, the interference pattern becomes more complex while being restricted within a regular polygon with peaks lying at its corners. The Gaussian peaks of the $Q$-functions (Figs. \[Kitten\_rational\] (a$_{4}$)-(c$_{4}$), (a$_{5}$)-(c$_{5}$)) appear as smoothed versions of the $W$-distributions.

\
\

Wigner entropy and negativity {#Wigner_entropy}
-----------------------------

It is also of interest to study the quantum entropy based on the modulus of the Wigner distribution $|W(\beta, \beta^{*})|$ \[[@SKD2012]\] which is a nonnegative quantity: S\_[W]{}=-|W(, \^[\*]{})| |W(, \^[\*]{})| d\^2 . \[S\_W\] As the $W$-distribution contains more information on the phase space structure of a quantum state than its smoothed analog the $Q$-function, a comparative study of the Wigner entropy $S_{W}$ (\[S\_W\]) and the Wehrl entropy $S_{Q}$ (\[WehrlDef\]) is expected to throw a light on the nonclassicality of the state. It is evident from the Figs. \[Wigner\_fig\] (a), (b) that the time evolution of the Wigner entropy (\[S\_W\]) and the negativity parameter (\[negativity\]) have close kinship with each other. This was observed for certain oscillator density functions in \[[@SKD2012]\]. In this sense the Wigner entropy $S_W$ reveals the extent of nonclassicality of a quantum density matrix. We distinguish between two possible scenarios depending upon the qubit-oscillator coupling strength. [(i)]{} In the strong coupling regime $(\lambda/\omega \lesssim 0.1)$ we observe (Fig. \[Wigner\_fig\] (a)) a periodic structure that may be identified with the revival and collapse of the qubit density matrix elements reflecting the energy exchanges between the qubit and the oscillator degrees of freedom. The collapse, say, of the qubit density matrix elements coincides with a wider spread of the phase space distributions of the oscillator. The order of the time period of the revival and collapse of the qubit matrix elements, obtained via retaining up to the linear terms in the Laguerre polynomials, is $O(2 \pi / x \widetilde{\Delta})$. Therefore the variables such as the Wigner entropy $S_{W}$, the negativity $\delta_{W}$, and the Wehrl entropy $S_Q$ display similar periodic patterns in the said time scale. For a dominant value of $\delta_{W}$ the quantum interference effects are overwhelming, and we, expectedly, find $S_{W} > S_{Q}$, as an increased negativity necessitates an increment in the magnitude $|W(\beta, \beta^{*})|$ for maintaining the normalization property (\[W\_normalization\]). This, in turn, leads to increased value of the entropy $S_{W}$. On the other hand, for a low negativity domain $\delta_{W} \ll 1$ the inequality is reversed: $S_{W} < S_{Q}$. The underlying reason is that the $Q$-function is obtained from the $W$-distribution (\[Q\_int\_W\]) after suitable smearing with a positive definite Gaussian kernel on the phase space, and therefore it incorporates less information on the quantum state than the latter \[[@MF2000]\]. [(ii)]{} In the ultra-strong coupling regime $(\lambda/\omega \gg 0.1)$ all modes for the qubit-oscillator interaction with incommensurate frequencies are excited and a fully randomized interference pattern evolves very quickly. Smearing the clear periodic structures observed earlier these large number of interaction-generated modes lead to quasi stationary values of the phase space observables (Fig. \[Wigner\_fig\] (b)). However, despite the statistical stabilization of the occupation on the phase space the nonclassicality of the state *remains prominent* due to the interferences occurring between multiple modes. These interferences necessarily develop (Fig. \[Wigner\_fig\] (c)) significant domains on the phase space with negative values of the $W$-distribution leading to dominant values of the negativity parameter $\delta_{W}$. The average value of $\delta_{W}$ increases with that of the coupling strength as more interfering modes come into existence. As mentioned before, this results in an increment of the entropy $S_{W}$. In the quasi stationary state the negativity $\delta_{W}$ is stochastically preserved. The stochastic stabilization of $\delta_{W}$ occurs after a suitable decoherence time. Therefore in Fig. \[Wigner\_fig\] (b) we observe that except for a brief initial period the Wigner entropy $S_W$ is consistently more than the Wehrl entropy $S_{Q}$, even though the nonnegative $Q$-function may be viewed as a smeared form of the $W$-distribution. Our results (Figs. \[Wigner\_fig\] (a), (b)) suggest that when the negativity $\delta_{W}$ assumes more than a threshold value $\delta_{W} \gtrsim 0.4$ the quantum fluctuations ensure the entropy relation: $S_{W} > S_{Q}$.

Another feature revealed in Figs. \[Wigner\_fig\] (a) and (b) is that the relative fluctuations of the Wigner entropy $S_{W}$ overwhelms that of the Wehrl entropy $S_Q: |\delta S_Q/S_Q| \ll |\delta S_W/S_W|$. It signifies that the quantum interferences in the evolution of the $W$-distribution induce rapid reversal of its sign, whereas the fluctuations of different modes are significantly attenuated in the smoothing introduced towards obtaining the $Q$-function.

Evolution to squeezed states {#evolution_squeezing}
----------------------------

A feature of nonclassicality such as squeezing is observed during the time evolution of the state due to the quantum interferences between various modes. At relatively small value of the phase space separation $(\alpha \lesssim 1)$ and for the strong coupling regime $(\lambda \sim 0.1\, \omega)$, the squeezing is noticed both for the large detuning (Fig. \[squeezing\]) and the resonant frequencies. The plot of the $W$-distribution (Fig. \[squeezing\] (c)) makes the quadrature squeezing evident. The signature of the squeezing is observed (Fig. \[squeezing\] (d)) when the variance $V_{\theta}$ of the quadrature variable, say at $\theta=0$, is rendered less than its classical value $1/4$. It follows from the polar plot (Fig. \[squeezing\] (d)) that at the scaled time $\omega t=28.80$ the quadrature variance $V_{\theta}$ reaches a minimum value $0.15967$ at an angle $\theta=179.55^{\circ}$. The polar angle at which the minimum quadrature variance is realized varies with time and depends on the dynamical state of interference of the quantum modes. A qualitative understanding of the squeezing of the state may be described as follows. Passing to the interaction picture for the Hamiltonian (\[qOH\]), an effective Hamiltonian may be obtained *à la* \[[@JJ2007]\] in an order by order perturbation theory. This effective Hamiltonian contains two photon terms $(a^{2}, a^{\dagger\,2})$ at the order $O(x)$ of the coupling strength. These two photon terms give rise to the squeezing of the state. An increase in the coupling $\lambda$ first enhances the squeezing as it augments the strength of the two photon terms in the effective Hamiltonian. Multiple photon terms, however, soon appear in the effective Hamiltonian \[[@JJ2007]\] with increased value of the coupling strength. The resulting randomness of the phase relationships of the higher order terms eliminates the squeezing property in the ultra-strong coupling limit. It is obvious in Fig. \[squeezing\] (a) that the quadrature squeezing, even though it may be present during part of the dynamical evolution of the oscillator state, is not realized throughout the oscillatory cycle. The quadrature variance remains above the threshold value: $V_{\theta = 0} > 1/4$ during the major part of the evolution of the state. The instances of squeezing decreases with increasing $\lambda$. To illustrate the feature we plot (Fig. \[squeezing\] (e)) the time-averaged quadrature variance $\langle V_{\theta = 0}\rangle_{\mathrm{time\, av.}}$ w.r.t. the coupling strength $\lambda$. We observe that, for the parametric regime studied here $\omega \gtrsim \Delta$, the time-averaged quadrature variance has no squeezing property: $\langle V_{\theta = 0}\rangle_{\mathrm{time\, av.}} > 1/4$, and it smoothly increases with the rising $\lambda$. The effective Hamiltonian approach \[[@JJ2007]\], however, suggests that the two photon terms $(a^{2}, a^{\dagger\,2})$ underlying the squeezing of the state survive the smearing of the fast oscillatory terms only in the limit of the large qubit frequency: $\omega \ll \Delta$. We will return to this topic somewhere else.

It is interesting to note that the squeezed state described in Fig. \[squeezing\] represents almost pure state of the oscillator at the given time as its von Neumann entropy $S\lvert_{\omega t=28.80}\,=\,0.15690$ is much less than its maximal value. The almost pure state of the oscillator is reciprocated by the corresponding nearly pure state of the qubit. For the sake of completeness we include the qubit density matix (\[Qubit\_DenMa\]) at the given time:

0.90760 & -0.10668 + i 0.19300\
-0.10668 - i 0.19300 & 0.09240

. \[Qubit\_Pure\] Its eigenvalues and corresponding eigenvectors read: $\lambda_{1}=0.96342$, $\lambda_{2}=0.03658, \; \ket{\lambda_{1}}=(0.46897- 
i \; 0.84843)\ket{1}-0.24543\ket{-1}$ and $\ket{\lambda_{2}}=(-0.11873+ i\; 0.21480)\ket{1} 
 - 0.96941 \ket{-1}$. The magnitude of the largest eigenvalue of the qubit density matix (\[Qubit\_DenMa\]) may be regarded as measure of the purity of the state. For smaller values of the coherent state amplitude $\alpha \lesssim 0.1$ the purity of the state may exceed $99\%$. Appearance of a squeezed state requires emergence of appropriate phase relations between various quantum modes that is facilitated by the proximity to a pure state. A statistical mixture of a large number of pure states is likely to destroy the phase relationships and, consequently, the squeezing property is eliminated.

Mandel parameter and nonclassicality {#Mandel}
------------------------------------

In experiments allowing a direct detection of photons the Mandel parameter \[[@M1979]\] Q\_[M]{} = ()\^[2]{} / -1 \[Q\_M\] is a convenient tool for studying the classical-quantum boundary. For the coherent state the number operator follows the Poissonian statistics with its signature property $Q_{M} = 0$ that is regarded as the threshold for the classical characteristics. When the quantum correlations in the system suppress the fluctuations in the photon number, it assumes negative values: $Q_{M} < 0$ and captures the sub-Poissonian behavior of the photon statistics that acts as a measure of nonclassicality. The mean value (\[avg\_photon\_number\]) and the variance (\[variance\_num\]) of the photon number operator obtained in our model now provide a direct construction of the Mandel parameter $Q_{M}$. For the coupling strength $\lambda \lesssim 0.1\, \omega$ transient sub-Poissonian photon statistics is manifest (Fig. \[Mandel\_Evolution\]) during parts of the evolution of the system. The time evolution of the Mandel parameter $Q_{M}$ in Fig. \[Mandel\_Evolution\] (a), however, passes alternately between the classical and the nonclassical regimes. The said sub-Poissonian behavior within short sample intervals is offset by the super-Poissonian photon statistics $Q_{M} > 0$ during most of the time evolution of the oscillator so that the averaged value over a time scale $t_{\mathrm{avg}} \gtrsim 2 \pi/(x \widetilde {\Delta})$ characterising the qubit-oscillator energy fluctuations always remains nonnegative. The classicality of the photon statistics, therefore, emerge via the time-averaging procedure even though at a shorter time scale the underlying quantum nature of the photon emission process is evident. At higher coupling strengths Fig. \[Mandel\_Evolution\] (b), (c) the randomization of the phases generated by a large number of incommensurate modes with the frequencies $\{O(x^{n} \widetilde{\Delta})| n = 1,2,\ldots\}$ partially erase the quantum correlations, while making the negativity of $Q_{M}$ less common. The effective Hamiltonian approach \[[@JJ2007]\], as noted in Subsec. \[evolution\_squeezing\], produce multiple photon terms that give rise to the negativity of $Q_{M}$. Such terms, however, contribute after the smoothing of the rapidly fluctuating components only in the limit $\Delta \gg \omega$. The time-averaged properties the Mandel parameter $Q_{M}$ is summarized in Fig. \[Mandel\_Evolution\] (d), where, in the coupling regime $\lambda \lesssim 0.1\, \omega$ the time-averaged value $\langle Q_{M} \rangle_{\mathrm {time\; av.}} \sim 0$ is maintained so that the overall photon statistics effectively remains to be Poissonian. After a brief transition zone the increasing coupling strength ($ \lambda \gg 0.1\, \omega$) triggers a full randomization of the phase relations among the interfering modes, which now constitute an effective bath leading to progressive emergence of classical properties. The resultant stochasticity introduces a scaling behavior for the averaged value of the Mandel parameter: $\langle Q_{M} \rangle_{\mathrm {time\; av.}} \propto (\lambda - \lambda_{0})$ with a suitable $\lambda_{0}$ that can be determined from the Fig. \[Mandel\_Evolution\] (d).

Conclusion
==========

Employing the generalized rotating wave approximation we have studied an interacting qubit-oscillator bipartite system for both the strong and the ultra-strong coupling domains for the choice of an initial hybrid Bell state. The evolution of the reduced density matrices of the qubit and the oscillator are obtained via the partial tracing of the complementary part of the respective degrees of freedom. On the oscillator phase space its density matrix furnishes the diagonal $P$-representation that is highly singular due to the presence of rapidly oscillating derivatives of the $\delta$-function. Two successive smoothing performed by the Gaussian kernels on the singular $P$-representation produce first the Wigner $W$-distribution, and then the Husimi $Q$-function. The quasi probability $W$-distribution admits negative values due to the quantum interferences. Its negativity measure $\delta_{W}$ marks the departure of the state from a positive definite distribution on the phase space. The nonnegative $Q$-function provides the Wehrl entropy $S_{Q}$ that acts as measure of delocalization on the oscillator phase space. Even though the $W$-distribution may be thought of encapsulating more informations than the $Q$-function, the Wigner entropy $S_{W}$ defined via the magnitude $|W|$ overwhelms the Wehrl entropy $S_{Q}$ whenever the negativity measure $\delta_{W}$ is dominant. The presence of multiple time scales induced by the interaction introduces a novel feature in the generation of ‘kitten’ states. In the coupling range $\lambda \lesssim 0.05 \,\omega$ the long-term quasi periodicity is observed to follow via the terms $O(x^{2})$ in the interaction Hamiltonian. In this domain the $Q$-function evolves at rational fractions of $T_{\mathrm{long}}$ to a collection of uniformly separated Gaussian peaks representing the ‘kitten’ states. A shorter time scale $O((x \widetilde{\Delta})^{-1})$ now causes further bifurcation of the Gaussian peaks coinciding with the collapse of the qubit density matrix elements. In the chaotic $\lambda \sim \omega$ regime all interaction-dependent modes of frequencies $\{O(x^{n}\widetilde{\Delta})| n = 0, 1, \ldots\}$ with a randomization of their phases. The decoherence time may be estimated as the transient production time of $S_{Q}$ as it approaches its stochastic stabilization. Using the asymptotic behavior of the associated Laguerre polynomials $L_{n}^{(j)}(x)$ the decoherence time is estimated as proportional to $\sqrt{\lambda}$. Nonclassical features such as squeezing and negativity of the Mandel parameter arise due to appearance of multiple photon terms induced by the interaction. In the parametric regime studied here $\omega \gtrsim \Delta$, such effects do not survive a suitable coarse graining process that smooths the rapidly oscillatory components in the fluctuations.

\[summary\]

Acknowledgement {#acknowledgement .unnumbered}
===============

One of us (VY) acknowledges the support from DST (India) under the INSPIRE Fellowship scheme.

[99]{}

E.T. Jaynes, F.W. Cummings, Proc. IEEE [**51**]{}, 89 (1963).

A.D.Armour, M.P. Blencowe, K.C. Schwab, Phys. Rev. Lett. [**88**]{}, 148301 (2002).

A.A. Anappara, S.D. Liberato, A. Tredicucci, C. Ciuti, G. Biasiol, L. Sorba, F. Beltram, Phys. Rev. B [**79**]{}, 201303 (2009).

T. Niemczyk, F. Deppe, H. Huebl, E.P. Menzel, F. Hocke, M.J. Schwarz, J.J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, R. Gross, Nature Physics [**6**]{}, 772 (2010).

I. Buluta, S. Ashhab, F. Nori, Rep. Progr. Phys. 74 (2011) 104401.

I.M. Georgescu, S. Ashhab, F. Nori, Rev. Modern Phys. 86 (2014) 153.

Z.L. Ziang, S. Ashab, J.Q. You, F. Nori, Rev. Modern Phys. 85 (2013) 623.

E.K. Irish, J. Gea-Banacloche, J. Martin, K.C. Schwab, Phys. Rev. B 72, 195410 (2005).

S. Ashhab, F. Nori, Phys. Rev. A 81, 042311 (2010).

E.K. Irish, Phys. Rev. Lett. [**99**]{}, 173601 (2007).

C. Guerlin, J. Bernu, S. Deléglise, C. Sayrin, S. Gleyzes, S. Kuhr, M. Brune, J.M. Raimond, S. Haroche, Nature (London) [**448**]{}, 889 (2007).

J. Park, H. Saunders, Y. Shin, K. An, H. Jeong, Phys. Rev. A[**85**]{}, 022120 (2012).

S.W. Lee, H. Jeong, Phys. Rev. A[**87**]{}, 022326 (2013).

P. van Loock, W.J. Munro, K. Nemoto, T.P. Spiller, T.D. Ladd, S.L. Braunstein, G.J. Milburn, Phys. Rev. A[**78**]{}, 022303 (2008).

Y.X. Liu, L.F. Wei, F. Nori, Phys. Rev. A [**71**]{}, 063820 (2005).

U.L. Anderson, J.S. Neergard-Nielsen, P. van Loock, A. Furusawa, *Hybrid quantum information processing*, arXiv:1409.37 \[quant-ph\] (2014).

H. Jeong, S. Zavatta, M. Kang, S.W. Lee, Nature Photonics [**8**]{}, 564 (2014).

O. Morin, K. Huang, J. Liu, H. Le Jeannie, C. Fabre, J. Laurat, Nature Photonics [**8**]{}, 570 (2014).

Y.X. Liu, L.F. Wei, F. Nori, Europhys. Lett. [**67**]{}, 941 (2004).

L.F. Wei, Y.X. Liu, M.J. Storcz, F. Nori, Phys. Rev. A [**73**]{}, 052307 (2006).

M. Hofheinz, E.M. Weig, M. Ansmann, R.C. Bialczak, E. Lucero, M. Neeley, A.D. O’Connell, H. Wang, J.M. Martinis, A.N. Cleland, Nature [**454**]{}, 310 (2008).

M. Hofheinz, H. Wang, M. Ansmann, R.C. Bialczak, E. Lucero, M. Neeley, A.D. O’Connell, D. Sank, J. Wenner J.M. Martinis, A.N. Cleland, Nature [**459**]{}, 546 (2009).

W.P. Schleich, *Quantum Optics in Phase Space*, Wiley-VCH, Berlin(2001).

A. Wehrl, Rev. Mod. Phys. [**50**]{}, 221 (1978).

A. Miranowicz, R. Tanas, S. Kielich, Quant. Opt. [**2**]{}, 253 (1990).

I. Jex, A. Orlowski, J. Mod. Phys. [**41**]{}, 2301 (1994).

R. Tanas, A. Miranowicz, T. Gantsog, Phys. Scripta [**T48**]{}, 53 (1993).

A. Miranowicz, J. Bajer, M.R.B. Wahiddin, N. Imoto, J. Phys. A [**34**]{}, 3887 (2001).

P. Sadeghi, S. Khademi, A.H. Daroone Phys. Rev. A [**86**]{}, 012119 (2012).

L. Mandel, Opt. Lett. [**4**]{}, 205 (1979).

H. Araki, E.H. Lieb, Comm. Math. Phys. [**18**]{}, 160 (1970).

E.C.G. Sudarshan, Phys. Rev. Lett. [**10**]{}, 277 (1963).

R.J. Glauber, Phys. Rev. [**131**]{}, 2766 (1963).

C.L. Mehta, J. Phys.: Conf. Ser. [**196**]{}, 012014 (2009).

L. Mandel, Phys. Rev. Lett. [**49**]{}, 136 (1982).

R.E. Slusher, L.W. Hollberg, B. Yurke, J.C. Mertz, J.F. Valley, Phys. Rev. Lett. [**55**]{}, 2409 (1985).

H. Moya-Cessa and P.L. Knight, Phys. Rev. A [**48**]{}, 2479 (1993).

J. Meixner, Math. Z. [**44**]{}, 531 (1939).

A. Kenfack, K. Życzkowski, J. Opt. B: Quantum Semiclass. Opt. [**6**]{}, 396 (2004).

A. Sugita, H. Aiba, Phys. Rev. E [**65**]{}, 036205 (2002).

G.-L. Ingold, A. Wobst, C. Aulbach, P. Hänggi, Lecture Notes in Physics [**630**]{}, 85, eds. T. Brandes, S. Ketterman, Springer-Verlag (Berlin) (2003).

V. Buzek, C.H. Keitel, P.L. Knight, Phys. Rev. A [**51**]{}, 2594 (1995).

D.F.V. James, J. Jerke, Can. J. Phys. 85 (2007) 625.

S. Szegö, *Orthogonal Polynomials*, Amer. Math. Soc., Providence (1975).

G. Kirchmair, B. Vlastakis, Z. Legtas, S.E. Nigg, H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S.M. Girvin, R.J. Schoelkopf, Nature [**495**]{}, 205 (2013).

A. Miranowicz, M. Paprzycka, A. Pathak, F. Nori, Phys. Rev. A [**89**]{}, 033812 (2014).

G. Manfredi and M.R. Feix, Phys. Rev. E [**62**]{}, 4665 (2000).
