# Reservoir engineering of a mechanical resonator: generating a macroscopic
superposition state and monitoring its decoherence

Muhammad Asjad  School of Science and Technology, Physics Division, University
of Camerino, Camerino (MC), Italy  David Vitali  School of Science and
Technology, Physics Division, University of Camerino, Camerino (MC), and INFN,
Sezione di Perugia, Italy

###### Abstract

A deterministic scheme for generating a macroscopic superposition state of a
nanomechanical resonator is proposed. The nonclassical state is generated
through a suitably engineered dissipative dynamics exploiting the
optomechanical quadratic interaction with a bichromatically driven optical
cavity mode. The resulting driven dissipative dynamics can be employed for
monitoring and testing the decoherence processes affecting the nanomechanical
resonator under controlled conditions.

##  I Introduction

Quantum reservoir engineering generally labels a strategy which exploits the
non-unitary evolution of a system in order to generate robust quantum coherent
states and dynamics Diehl  _et al._ (2008). The idea is in some respect
challenging the intuitive expectation that in order to obtain quantum coherent
dynamics one should guarantee that the evolution is unitary at all stages. Due
to the noisy and irreversible nature of the processes which generate the
target dynamics, strategies based on quantum reservoir engineering are in
general more robust against variations of the parameters than protocols solely
based on unitary evolution Diehl  _et al._ (2008); Verstraete  _et al._
(2009). A prominent example of quantum reservoir engineering is laser cooling,
achieving preparation of atoms and molecules at ultralow temperatures by means
of an optical excitation followed by radiative decay Wineland  _et al._
(1978). The idea of quantum reservoir engineering has been formulated in Ref.
Poyatos  _et al._ (1996), and further pursued in Ref. Carvalho  _et al._
(2001). Proposals for quantum reservoir engineering of many-body systems have
been then discussed in the literature Diehl  _et al._ (2008); Verstraete  _et
al._ (2009) and first experimental realizations have been reported Syassen
_et al._ (2008); Krauter  _et al._ (2011).

In particular reservoir engineering has been proposed and already used Krauter
_et al._ (2011) for the generation of steady state nonclassical states, such
as linear superposition (Schrödinger cat) states Poyatos  _et al._ (1996);
Carvalho  _et al._ (2001) or entangled states in microwave cavities Pielawa
_et al._ (2007, 2010). In this case one has the advantage that the desired
target state is largely independent of the specific initial states, and at the
same time is robust with respect to a large class of decoherence processes.
These ideas have been recently extended also to the field of cavity
optomechanics for the generation of entangled states of two cavity modes Wang
and Clerk (2013), and of two mechanical modes Tan  _et al._ (2013).

Here we apply reservoir engineering for the deterministic generation of robust
macroscopic superpositions of coherent states of a mechanical resonator (MR).
First proposals for the generation of superposition states exploited the
intrinsic nonlinearity of radiation pressure interaction Bose  _et al._
(1999); Marshall  _et al._ (2003) but are hard to realize due to the extremely
weak nonlinear coupling. More recent proposals focused on the conditional
generation of those linear superposition states Paternostro (2011); Vanner
_et al._ (2013), exploiting for example the effective measurement of the MR
position _squared_ in order to generate a superposition of two spatially
separated states Romero-Isart  _et al._ (2011, 2011); Jacobs  _et al._ (2009,
2011). These latter schemes do not suffer from weak radiation pressure
nonlinearities, but are probabilistic and strongly dependent upon the
efficiency of the conditional measurement on the optical mode. As underlined
above, the generation of a linear superposition state through reservoir
engineering is instead deterministic and extremely robust, because the state
is reached asymptotically as a result of a dissipative irreversible evolution,
and is less sensitive to the details of the preparation process.

Here we propose to generate a superposition of coherent states of a MR by
exploiting the nonlinearity associated with the _quadratic_ interaction of the
MR with an optical cavity mode, appropriately driven by a bichromatic field
(see also Ref. Tan  _et al._ (2013)). We study the resulting dynamics,
determined by the joint action of the engineered reservoir realized by the
driven cavity mode and of the standard (and unavoidable) thermal reservoir of
the MR. We show that a high-quality superposition state can be generated in a
transient time interval, which then decoheres at longer times due to the
action of thermal reservoir. The present scheme is particularly useful for
monitoring decoherence processes affecting nanomechanical resonators,
similarly to what has been done for cavities Deléglise  _et al._ (2008) or
trapped ions Myatt  _et al._ (2000), and could also be useful for testing
alternative decoherence models (see Ref. Romero-Isart  _et al._ (2011) and
references therein).

In Sec. II we describe the properties of the required engineered reservoir. In
Sec. III we show how to engineer such a reservoir by tailoring the
optomechanical interaction with a bichromatically driven cavity mode. Sec. IV
describes the resulting dynamics under realistic scenarios, and we verify that
a superposition state can be efficiently generated in a transient time
interval and that its decoherence can be monitored. Sec. V is for concluding
remarks.

##  II The desired dissipative evolution

Let $\rho$ be the reduced density matrix of the MR, and
$\rho_{\infty}=|\psi_{\infty}\rangle\langle\psi_{\infty}|$ the target linear
superposition state we want to generate in the steady state of the MR, the so-
called even Schrödinger cat state

$|\psi_{\infty}\rangle=(|\beta\rangle+|-\beta\rangle)/\mathcal{N},$ (1)

where $|\beta\rangle$ denotes a coherent state of the MR with complex
amplitude $\beta$ and $\mathcal{N}=\sqrt{2[1+\exp(-2|\beta|^{2})]}$ is the
normalization constant. Reservoir engineering means in the present case
tailoring the interaction with the optical cavity mode in order to have an
effective reduced dynamics of the MR described by the master equation

$\frac{\partial}{\partial t}\rho=\mathcal{L}\rho\,,$ (2)

for which $\rho_{\infty}$ is a fixed point, namely,

$\mathcal{L}\rho_{\infty}=0\,.$ (3)

A simple solution is to take the Lindbladian ${\mathcal{L}}$

$\displaystyle\mathcal{L}\rho$ $\displaystyle=$
$\displaystyle\Gamma{\mathcal{D}}(C)\rho$ (4)
$\displaystyle{\mathcal{D}}(C)\rho$ $\displaystyle=$
$\displaystyle\left(2C\rho C^{\dagger}-C^{\dagger}C\rho-\rho
C^{\dagger}C\right),$ (5)

with $\Gamma$ a model-dependent rate, and with the operator $C$ such that
$|\psi_{\infty}\rangle$ is eigenstate of $C$ with zero eigenvalue. Such a
condition is satisfied by

$C=(b^{2}-\beta^{2}),$ (6)

where $b$ is the annihilation operator of the MR, and $\beta$ is just the
complex amplitude of the target linear superposition state.

Notice that state $\rho_{\infty}$ is not the unique solution of Eq. (3)
because any coherent or incoherent superposition of $|\beta\rangle$ and
$|-\beta\rangle$ solves Eq. (3). However for our purposes it is sufficient
that, at least for a physically realizable class of initial states of the MR,
the dissipative evolution asymptotically drives it only to $\rho_{\infty}$ and
not to other states of the convex set of states ${\mathcal{C}_{L}}$ defined by
Eq. (3). In this respect one can profit from an additional symmetry of the
Lindbladian of Eq. (4), i.e., the fact that it commutes with the parity
operator ${\mathcal{P}}=(-1)^{b^{\dagger}b}$, and therefore ${\mathcal{P}}$ is
conserved as long as the dynamics is driven by $\mathcal{L}$ only or at least
parity non-conserving terms are negligible in the time evolution generator. In
such a case, since $|\psi_{\infty}\rangle$ is the unique pure state of
${\mathcal{C}_{L}}$ which is even, that is, eigenstate of ${\mathcal{P}}$ with
eigenvalue $+1$, the asymptotic steady state of the MR will also have parity
$+1$. In particular it is possible to see that if the initial state is pure
and even, the asymptotic state will be $|\psi_{\infty}\rangle$. A natural case
of this kind is provided by a vacuum initial state $|0\rangle\langle 0|$,
which is obtained if the MR is initially cooled to its ground state.

Therefore our goal is to generate the effective reduced dynamics of the MR
driven by the above Lindbladian of Eqs. (4)-(6) when the cavity mode is
adiabatically eliminated. In practice however, the MR dynamics will be
affected not only by the cavity mode “engineered reservoir” but also by the
standard thermal reservoir. Therefore we have to establish the effect of these
undesired latter terms, and to determine if and when they are negligible.

##  III Engineering the dissipative process

Our starting point is the Hamiltonian of an optomechanical system formed by a
driven cavity mode interacting _quadratically_ with a MR. Such a quadratic
interaction is achieved for example in a membrane-in-the-middle (MIM) setup,
when the membrane is placed at a node, or exactly at an avoided crossing point
within the cavity Thompson  _et al._ (2008); Sankey  _et al._ (2010); Karuza
_et al._ (2013). Alternatively, such a quadratic coupling can be obtained by
trapping levitating nanoparticles around an intensity maximum of a cavity mode
Barker (2010); Li  _et al._ (2011); Gieseler  _et al._ (2012); Kiesel  _et
al._ (2013). We assume that the cavity is bichromatically driven, that is

$\displaystyle H=\hbar\omega_{m}b^{\dagger}b+\hbar\omega_{c}a^{\dagger}a+\hbar
g_{2}a^{\dagger}a(b+b^{\dagger})^{2}$ (7)
$\displaystyle+\mathrm{i}\hbar\left[(E_{0}e^{-\mathrm{i}\omega_{\mathrm{L}}t}+E_{1}e^{-\mathrm{i}(\omega_{\mathrm{L}}+\Omega)t})a^{\dagger}\right.$
$\displaystyle\left.-(E_{0}e^{\mathrm{i}\omega_{\mathrm{L}}t}+E_{1}e^{\mathrm{i}(\omega_{\mathrm{L}}+\Omega)t})a\right],$

where $\omega_{m}$ is the resonance frequency of the MR, $\omega_{c}$ the
cavity mode frequency,
$E_{0}=\sqrt{2P_{0}\kappa_{0}/\hbar\,\omega_{\mathrm{L}}}$,
$E_{1}=\sqrt{2P_{1}\kappa_{0}/\hbar\,(\omega_{\mathrm{L}}+\Omega)}$, with
$\kappa_{0}$ the cavity decay rate through the input mirror, and $P_{0}$ and
$P_{1}$ (with $P_{0}\gg P_{1}$) the respective input power at the two driving
frequencies. $g_{2}$ is the quadratic optomechanical coupling rate, which is
equal to $g_{2}=\Theta(\partial^{2}\omega_{c}/\partial
z_{0}^{2})(\hbar/2m\omega_{m})$ in the MIM case, with $\Theta$ the transverse
overlap between the mechanical and optical mode at the membrane, and $m$ the
membrane mode effective mass Karuza  _et al._ (2013).

We then move to the frame rotating at the main laser frequency
$\omega_{\mathrm{L}}$, where the system Hamiltonian becomes

$\displaystyle H=\hbar\omega_{m}b^{\dagger}b+\hbar\Delta_{0}a^{\dagger}a+\hbar
g_{2}a^{\dagger}a(b+b^{\dagger})^{2}$ (8)
$\displaystyle+\mathrm{i}\hbar\left[(E_{0}+E_{1}e^{-\mathrm{i}\Omega
t})a^{\dagger}-(E_{0}+E_{1}e^{\mathrm{i}\Omega t})a\right],$

where $\Delta_{0}=\omega_{c}-\omega_{\mathrm{L}}$ is the cavity mode detuning.
The dynamics is however also driven by fluctuation-dissipation processes
associated with the coupling of the cavity mode with the optical vacuum field
outside the cavity, and of the MR with its thermal reservoir characterized by
a temperature $T$ and a mean thermal phonon number
$\bar{n}=\left[\exp(\hbar\omega_{m}/k_{B}T)-1\right]^{-1}$. In the usual
Markovian approximation Gardiner and Zoller (2000), we have that optical
dissipation is described by

$\kappa_{\rm T}{\mathcal{D}}(a)\rho_{om},$ (9)

where $\rho_{om}$ is the density matrix of the total optomechanical system,
and $\kappa_{\rm T}$ is the total cavity decay rate, while mechanical
fluctuation-dissipation effects are described by the following terms in the
master equation Gardiner and Zoller (2000)

$\frac{\gamma_{\rm
m}}{2}(\bar{n}+1){\mathcal{D}}(b)\rho_{om}+\frac{\gamma_{\rm
m}}{2}\bar{n}{\mathcal{D}}(b^{\dagger})\rho_{om},$ (10)

where $\gamma_{m}=\omega_{m}/Q_{m}$ is the mechanical damping, and $Q_{m}$ is
the mechanical quality factor. Therefore the time evolution of the system is
described by the following general master equation

$\displaystyle\frac{\partial}{\partial
t}\rho_{om}=-\frac{\mathrm{i}}{\hbar}\left[H,\rho_{om}\right]+\kappa_{\rm
T}{\mathcal{D}}(a)\rho_{om}$ (11) $\displaystyle+\frac{\gamma_{\rm
m}}{2}(\bar{n}+1){\mathcal{D}}(b)\rho_{om}+\frac{\gamma_{\rm
m}}{2}\bar{n}{\mathcal{D}}(b^{\dagger})\rho_{om},$

with $H$ given by Eq. (8).

The intense driving associated with the laser field at the carrier frequency
$\omega_{\mathrm{L}}$ generates a stationary intracavity state of the cavity
mode with large coherent amplitude

$\alpha_{s}=\frac{E_{0}}{\kappa_{\rm T}+\mathrm{i}\Delta_{0}},$ (12)

and it is convenient to look at the dynamics of the quantum fluctuations of
the cavity mode, performing the displacement $a=\alpha_{s}+\delta a$. After
some algebra and using Eq. (12), the master equation of Eq. (11) becomes

$\displaystyle\frac{\partial}{\partial
t}\rho_{om}=-\frac{\mathrm{i}}{\hbar}\left[H_{\delta},\rho_{om}\right]+\kappa_{\rm
T}{\mathcal{D}}(\delta a)\rho_{om}$ (13) $\displaystyle+\frac{\gamma_{\rm
m}}{2}(\bar{n}+1){\mathcal{D}}(b)\rho_{om}+\frac{\gamma_{\rm
m}}{2}\bar{n}{\mathcal{D}}(b^{\dagger})\rho_{om},$

with the modified Hamiltonian

$\displaystyle H_{\delta}=\hbar\tilde{\omega}_{m}b^{\dagger}b+\hbar
g_{2}|\alpha_{s}|^{2}\left(b^{2}+b^{\dagger\,2}\right)+\hbar\Delta_{0}\delta
a^{\dagger}\delta a$
$\displaystyle+\mathrm{i}\hbar\left[E_{1}e^{-\mathrm{i}\Omega t}\delta
a^{\dagger}-E_{1}e^{\mathrm{i}\Omega t}\delta a\right]$ (14)
$\displaystyle+\hbar g_{2}\left(\alpha_{s}^{*}\delta a+\alpha_{s}\delta
a^{\dagger}\right)(b+b^{\dagger})^{2}+\hbar g_{2}\delta a^{\dagger}\delta
a(b+b^{\dagger})^{2},$

where $\tilde{\omega}_{m}=\omega_{m}+2g_{2}|\alpha_{s}|^{2}$ is the
renormalized mechanical frequency.

We now take $\Omega=\Delta_{0}$, i.e., we assume that the second, less intense
beam is exactly resonant with the cavity mode, and move to the interaction
picture with respect to

$H_{0}=\hbar\tilde{\omega}_{m}b^{\dagger}b+\hbar\Delta_{0}\delta
a^{\dagger}\delta a.$ (15)

Within such a picture, the dissipative terms in the master equation of Eq.
(13) remain unchanged, while the Hamiltonian becomes

$\displaystyle H_{\delta}^{int}=\hbar
g_{2}|\alpha_{s}|^{2}\left(b^{2}e^{-2\mathrm{i}\tilde{\omega}_{m}t}+b^{\dagger\,2}e^{2\mathrm{i}\tilde{\omega}_{m}t}\right)+\mathrm{i}\hbar\left[E_{1}\delta
a^{\dagger}-E_{1}\delta a\right]$ $\displaystyle+\hbar
g_{2}\left(\alpha_{s}^{*}\delta ae^{-\mathrm{i}\Delta_{0}t}+\alpha_{s}\delta
a^{\dagger}e^{\mathrm{i}\Delta_{0}t}\right)(be^{-\mathrm{i}\tilde{\omega}_{m}t}+b^{\dagger}e^{\mathrm{i}\tilde{\omega}_{m}t})^{2}$
$\displaystyle+\hbar g_{2}\delta a^{\dagger}\delta
a(be^{-\mathrm{i}\tilde{\omega}_{m}t}+b^{\dagger}e^{\mathrm{i}\tilde{\omega}_{m}t})^{2}.$
(16)

We have made no approximation up to now. We now take the following _resonance
condition_ , $\Delta_{0}=\Omega=2\tilde{\omega}_{m}$, which means that the
second driving beam is resonant not only with the cavity, but also with the
second order sideband of the carrier beam at $\omega_{\mathrm{L}}$, and make
the two following approximations: i) we neglect the last, higher order
interaction term $\hbar g_{2}\delta a^{\dagger}\delta
a(be^{-\mathrm{i}\tilde{\omega}_{m}t}+b^{\dagger}e^{\mathrm{i}\tilde{\omega}_{m}t})^{2}$,
which is justified whenever $|\delta a|\ll|\alpha_{s}|$; ii) we make the
rotating wave approximation (RWA) and neglect all the fast-oscillating terms
at $\tilde{\omega}_{m}$ and $2\tilde{\omega}_{m}$, which is justified in the
weak coupling limit $g_{2}|\alpha_{s}|\ll\tilde{\omega}_{m}$. The effective
interaction picture Hamiltonian of Eq. (16) reduces to

$H_{\rm eff}=\hbar g_{2}\alpha_{s}^{*}\delta
a\left(b^{\dagger\,2}-\mathrm{i}E_{1}/g_{2}\alpha_{s}^{*}\right)+{\rm H.C.},$
(17)

and therefore under the conditions specified above, the dynamics of the
optomechanical system is described by Eq. (13) with $H_{\delta}$ replaced by
$H_{\rm eff}$. An analogous effective optomechanical dynamics has been
considered in Ref. Tan  _et al._ (2013), where however the renormalization of
the mechanical frequency $\omega_{m}\to\tilde{\omega}_{m}$ has been neglected.

###  III.1 Reduced dynamics of the mechanical resonator

In the bad cavity limit, i.e., when $\kappa_{\rm T}\gg
g_{2}|\alpha_{s}|,\gamma_{m}\bar{n}$, the optical mode fluctuations $\delta a$
can be adiabatically eliminated because they quickly decay and their state
always remains close to the vacuum state (see for example Ref. Gardiner and
Zoller (2000), pag. 147 and Ref. Wiseman and Milburn (1993)). One gets the
following final effective master equation for the reduced density matrix of
the MR, $\rho$,

$\frac{\partial}{\partial t}\rho=\Gamma{\mathcal{D}}(C)\rho+\frac{\gamma_{\rm
m}}{2}(\bar{n}+1){\mathcal{D}}(b)\rho+\frac{\gamma_{\rm
m}}{2}\bar{n}{\mathcal{D}}(b^{\dagger})\rho,$ (18)

where $\Gamma=g_{2}^{2}|\alpha_{s}|^{2}/\kappa_{\rm T}$ and $C$ is just given
by Eq. (6), with the superposition state amplitude
$\beta^{2}=E_{1}/{\mathrm{i}}g_{2}\alpha_{s}$. The first term is just the
desired term, i.e., the engineered dissipative evolution able to drive the MR
asymptotically to the target superposition state $|\psi_{\infty}\rangle$.
However, the time evolution of the MR state is also driven by the second and
third terms which are due to the coupling with the thermal reservoir at
temperature $T$. The latter are “undesired” terms, because they drive the MR
to a thermal state rather than the desired even Schrödinger cat state, and
also because they do not conserve the parity. Due to the joint action of these
two dissipative evolutions, the asymptotic state achieved by the MR at long
times will be different from the desired even superposition state
$|\psi_{\infty}\rangle$. However, if $\Gamma\gg\gamma_{m}\bar{n}$ so that the
effect of the thermal reservoir is negligible, we expect that the target state
can be generated at least for a reasonable transient time interval around
$\bar{t}\sim 1/\Gamma$. This condition, together with the conditions
$|\alpha_{s}|\gg 1$, and $\kappa_{\rm T},\tilde{\omega}_{m}\gg
g_{2}|\alpha_{s}|$ which are needed for deriving Eq. (18), represent the
parameter conditions for realizing the robust generation of a superposition
state of the MR.

##  IV Results

Let us now we verify if and when the proposal is implementable in a state-of-
the-art optomechanical setup and a nanomechanical resonator can be prepared
with high fidelity, at least for a long-lived transient, in the macroscopic
superposition state $|\psi_{\infty}\rangle$. We consider parameter values
achievable in state-of-the-art MIM setups Thompson  _et al._ (2008); Sankey
_et al._ (2010); Karuza  _et al._ (2013); Wilson  _et al._ (2009); Purdy  _et
al._ (2012); Flowers-Jacobs  _et al._ (2012); Karuza  _et al._ (2012); Purdy
_et al._ (2013). For the mechanical resonator we take $\omega_{m}=10$ MHz,
$\gamma_{m}=0.1$ Hz (implying $Q_{m}=10^{8}$), $m=1$ ng and we can take
$\partial^{2}\omega_{c}/\partial z_{0}^{2}=2\pi\times 20$ GHz/nm2 Flowers-
Jacobs  _et al._ (2012), yielding $g_{2}\simeq 5$ Hz. We then take a laser
with frequency $\omega_{\rm L}=1.77\times 10^{15}$ Hz (corresponding to a
wavelength $\lambda=1064$ nm) and input power $P_{0}=40$ mW. We also choose a
cavity with total decay rate $\kappa_{\rm T}=10^{5}$ Hz and with decay rate
through the input mirror $\kappa_{0}\sim\kappa_{\rm T}/2$, yielding $E_{0}\sim
1.5\times 10^{11}$ Hz. The corresponding value of the intracavity amplitude
from Eq. (12) is $|\alpha_{s}|\sim 3.45\times 10^{3}$. As a consequence
$g_{2}|\alpha_{s}|\sim 1.7\times 10^{4}$ Hz, which therefore agrees with the
assumptions made. Moreover we have an effective decay rate
$\Gamma=g_{2}^{2}|\alpha_{s}|^{2}/\kappa_{\rm T}\sim 2.98$ kHz which is
reasonably larger than the thermal decay rate $\gamma_{m}\bar{n}$ as long as
$\bar{n}\lesssim 100$. Even though nontrivial, this latter condition is
achievable in current optomechanical experiments because cryogenic
environments at temperatures $T\simeq 10$ mK are feasible and, with the chosen
value $\omega_{m}=10$ MHz, this corresponds just to $\bar{n}\simeq 100$.

The amplitude $\beta$ of the target state is determined by $E_{1}$ and
therefore by the input power $P_{1}$. Assuming $P_{1}\sim 1$ pW, one gets
$E_{1}\sim 10^{6}$ Hz and therefore $|\beta|\sim 23.6$, which corresponds to a
quite macroscopic superposition state; here however, in order to verify
numerically the proposal in a not too large operational Hilbert space, we have
taken $P_{1}\sim 0.01$ pW, yielding $E_{1}\sim 10^{5}$ Hz and therefore
$|\beta|\sim 2.36$.

###  IV.1 Cat state generation starting from the mechanical ground state

As discussed above, we expect to generate a long-lived transient even
Schrödinger cat state of the MR when $\Gamma\gg\gamma_{m}\bar{n}$ and if we
start from the mechanical ground state, which is pure and even. Since in the
considered scenario it is very hard to go below $\bar{n}\sim 100$ with
cryogenic techniques only, this initial state could be achieved, at least in
principle, by first laser cooling the MR to its ground state, i.e., by first
considering a _linear_ optomechanical interaction with a cavity mode and
driving it on its first red sideband Marquardt  _et al._ (2007); Wilson-Rae
_et al._ (2007); Genes  _et al._ (2008); Chan  _et al._ (2011); Verhagen  _et
al._ (2012). Then, one should switch to the quadratic optomechanical
interaction (either by displacing the membrane or by driving a different
appropriate cavity mode) soon after ground state cooling is attained.

We have numerically solved the time evolution of the optomechanical system
density matrix $\rho_{om}$ as described by Eq. (13) with the Hamiltonian of
Eq. (17), starting from the mechanical ground state and the vacuum state for
the cavity mode fluctuations. Plots of the Wigner representation of the
reduced state $\rho$ of the MR at different times are shown in Fig. 1, which
refers to the set of parameters described above, and $\bar{n}=100$. These
plots confirm our expectations and that the state $|\psi_{\infty}\rangle$ with
$|\beta|\sim 2.36$ is generated in the transient regime $t\sim 1/\Gamma$ due
to the appropriate bichromatic driving and the quadratic optomechanical
interaction. The superposition state then decoheres on a time scale governed
by $\gamma_{m}\left(2\bar{n}+1\right)$. These results are consistent with
those of Ref. Tan  _et al._ (2013) which also studies the generation of a cat
state of a MR starting from the ground state in a bichromatically driven
quadratic optomechanical system by means of the Wigner function of the reduced
MR state.

Figure 1: Time evolution of the Wigner function of the reduced state of the MR
starting from an initial factorized state in which both the optical cavity
fluctuations and the mechanical mode are in their ground state. The set of
parameters is given in the text, and $\bar{n}=100$. (a) Wigner function of the
initial state of the MR; (b) Wigner function of the MR state at at time
$t=0.71/\Gamma=2.39\times 10^{-4}$ s; (c) Wigner function of the MR state at
time $t=100/\Gamma=0.03355$ s. The even superposition state is successfully
generated in a short time of the order of $1/\Gamma$, and it slowly loses its
nonclassical interference fringes at a longer timescale, of the order of
$\left[\gamma_{m}\left(2\bar{n}+1\right)\right]^{-1}$.

This qualitative analysis based on the Wigner function is confirmed by a
quantitative analysis based on the time evolution of the fidelity of the state
with respect to the target state $|\psi_{\infty}\rangle$. Rather than the more
common Uhlmann fidelity Uhlmann (1976); Jozsa (1994), in order to simplify the
numerical calculation, here we use the Hilbert-Schmidt fidelity introduced in
Ref. Wang  _et al._ (2008)

${\mathcal{F}}(\rho_{0},\rho_{1})=\frac{\left|{\rm
Tr}\left\\{\rho_{0}\rho_{1}\right\\}\right|}{\sqrt{{\rm
Tr}\left\\{\rho_{0}^{2}\right\\}{\rm Tr}\left\\{\rho_{1}^{2}\right\\}}}.$ (19)

When $\rho_{0}$ is pure, this fidelity coincides with the probability of
finding the state $\rho_{0}$ being in $\rho_{1}$, divided by the square root
of the purity $\sqrt{{\rm Tr}\left\\{\rho_{1}^{2}\right\\}}$. In Fig. 2 we
plot the time evolution of
${\mathcal{F}}(t)={\mathcal{F}}(\rho_{\infty},\rho(t))$ corresponding to the
same parameter condition of Fig. 1. The fidelity reaches a maximum
${\mathcal{F}}\simeq 0.9992$ at $t\simeq 1/\Gamma$ when an almost perfect cat
state is generated, which then decays so that ${\mathcal{F}}\simeq
1/\sqrt{2}\simeq 0.7$.

Figure 2: Plot of fidelity
${\mathcal{F}}(t)={\mathcal{F}}(\rho_{\infty},\rho(t))$ as a function of time.
The inset shows the behavior at short times. Parameters are those given in the
text and coinciding with those of Fig. 1.

###  IV.2 Cat state generation after two-phonon cooling

Fast switching from the linear optomechanical interaction needed for cooling
to the mechanical ground state to the quadratic optomechanical interaction
necessary for generating the even cat state is quite challenging in practical
experimental situations. However, one could exploit the quadratic interaction
also for pre-cooling the MR and avoid using a different cavity mode and
different driving field. To be more specific one could use the same
interaction Hamiltonian and parameter conditions described in the previous
section and consider the special case $E_{1}=\beta=0$, i.e., with the weak
resonant field turned off. In this case, the engineered interaction with the
cavity mode induces a two-phonon cooling process driving the MR to its ground
state. The joint dynamics in the presence of nonlinear two-phonon damping and
standard decay to the thermal equilibrium with $\bar{n}$ thermal phonons has
been already studied in Ref. Nunnenkamp  _et al._ (2010), where it is shown
that in the limit $\Gamma\gg\gamma_{m}\bar{n}$ we are considering, cooling is
good even though not perfect, being the MR steady state a mixture of the zero
and one phonon state, with probabilities $\rho_{11}(\infty)=n_{\rm
eff}=(4+1/\bar{n})^{-1}$ and $\rho_{00}(\infty)=1-\rho_{11}(\infty)$.
Therefore a feasible cat state generation protocol is to first cool the MR
with the two-phonon cooling process with $E_{1}=0$, and then switch on the
weak resonant field with $E_{1}\neq 0$ for generating the even cat state as
discussed above. We now see that despite the initial approximate $25\%$
probability of being in the odd one phonon state, the cat state generation
process is still quite efficient, showing that such a robust macroscopic
superposition can be generated in achievable quadratic optomechanical setups.

We have in fact numerically solved the master equation for the optomechanical
system density matrix $\rho_{om}$ of Eq. (13) with the Hamiltonian of Eq.
(17), now taking as initial state the vacuum state for the cavity mode
fluctuations and the above mixture of the zero and one phonon state for the
MR, using the same set of parameters of the previous subsection (we have
verified that with this set of parameters one actually cools the MR to this
mixture of states). Plots of the Wigner representation of the reduced state
$\rho$ of the MR at different times are shown in Fig. 3, which refer to
$\bar{n}=10$ and in Fig. 4, which refers to $\bar{n}=100$. In both cases the
target cat state is generated with high fidelity at $t\sim 1/\Gamma$, despite
the residual excitation in the one-phonon state. This is confirmed by time
evolution of the fidelity of the state with respect to the target state
$|\psi_{\infty}\rangle$, which is shown in Fig. 5 for $\bar{n}=10$ (a) and
$\bar{n}=100$ (b). The fidelity reaches a maximum ${\mathcal{F}}\simeq 0.94$
at $t\simeq 1/\Gamma$ which does not depend upon $\bar{n}$ and then decays to
${\mathcal{F}}\simeq 1/\sqrt{2}\simeq 0.7$. The superposition state decoheres
to a mixture of two Gaussian states on a time scale governed by the thermal
decoherence rate given by $\gamma_{\rm
dec}=2\gamma_{m}|\beta|^{2}\left(2\bar{n}+1\right)$ Kennedy and Walls (1988);
Kim and Buzěk (1992); Brune  _et al._ (1996); Deléglise  _et al._ (2008).

Figure 3: Time evolution of the Wigner function of the reduced state of the MR
with mean phonon number $\bar{n}(0)=10$. (a) Wigner function of the initial
state of the MR at time t=0; (b) Wigner function of the MR state at time
$t=1/\Gamma=3.3547\times 10^{-4}$ s; (c) Wigner function of the MR state at
time $t=1000/\Gamma=0.3355$ s. The other parameters are given in the text and
coincide with those of Fig. 2.

Figure 4: Time evolution of the Wigner function of the reduced state of the MR
with mean phonon number $\bar{n}(0)=100$. (a) Wigner function of the initial
state of the MR at time t=0; (b) Wigner function of the MR state at time
$t=1/\Gamma=3.3547\times 10^{-4}$ s; (c) Wigner function of the MR state at
time $t=100/\Gamma=0.03355$ s. The other parameters are given in the text and
coincide with those of Fig. 2.

Figure 5: Plot of fidelity
${\mathcal{F}}(t)={\mathcal{F}}(\rho_{\infty},\rho(t))$ as a function of time
for (a) $\bar{n}=10$ and (b) $\bar{n}=100$. The other parameters are given in
the text and coincide with those of Fig. 2. The insets show the behavior at
short times.

###  IV.3 Approximate description of the progressive decoherence of the
generated superposition state

The above analysis shows that the combined action of the engineered reservoir
term with rate $\Gamma$ and the thermal reservoir terms with rate
$\gamma_{m}\bar{n}$, when $\Gamma\gg\gamma_{m}\bar{n}$, generates a
superposition state at time $t\simeq 1/\Gamma$ which then decoheres with
decoherence rate $2\gamma_{m}|\beta|^{2}\left(2\bar{n}+1\right)$. In
particular, Figs. 1-4 suggest that the MR decoheres to an asymptotic state
given by the mixture of the two coherent states
$|\pm\beta\rangle\langle\pm\beta|$, with $\beta=\sqrt{E_{1}/ig_{2}\alpha_{s}}$
just the amplitude of the target superposition state. To state it in other
words, the combined action of the engineered and “natural” reservoirs tends to
stabilize such a mixture of coherent states emerging after the decoherence
process. Taking into account the well-established theory of decoherence of
superposition of two coherent state in the presence of a thermal reservoir
Kennedy and Walls (1988); Kim and Buzěk (1992); Brune  _et al._ (1996);
Deléglise  _et al._ (2008), one is led to approximate the time evolution of
the reduced MR state after a transient time $t\geq t_{0}\simeq 1/\Gamma$ with
the following expression

$\displaystyle\rho_{\rm app}(t>t_{0})={\cal
N}(t-t_{0})^{-1}\left\\{|\beta\rangle\langle\beta|+|-\beta\rangle\langle-\beta|\right.$
(20)
$\displaystyle\left.+e^{-\left(1+2\bar{n}\right)\gamma_{m}\left(t-t_{0}\right)}\left[|\beta\rangle\langle-\beta|+|-\beta\rangle\langle\beta|\right]\right\\},$

with ${\cal
N}(t)=2\left[1+e^{-2|\beta|^{2}}e^{-\left(1+2\bar{n}\right)\gamma_{m}t}\right]$,
describing a decohering cat state, which decoheres to its corresponding
mixture just at the rate $2\gamma_{m}|\beta|^{2}\left(2\bar{n}+1\right)$.

We can check the validity of this approximate description at $t>t_{0}$ by
using again the Hilbert-Schmidt fidelity of Eq. (19) for measuring the overlap
between the actual reduced MR state $\rho(t)$ given by the solution of the
master equation of Eq. (13) and the approximate solution $\rho_{\rm app}$ of
Eq. (20). In Fig. 6 we plot the “distance” between the two states,
$D(t)=1-{\mathcal{F}}\left(\rho(t),\rho_{\rm app}(t)\right)$ for the same set
of parameters of Fig. 2, and we find a very good agreement for the proposed
solution. Therefore the generated Schrödinger cat state can be used for
verifying experimentally the decoherence processes affecting the
nanomechanical resonator and eventually testing alternative decoherence
models, as suggested in Ref. Romero-Isart  _et al._ (2011).

Figure 6: Plot of the distance between the actual solution of the master
equation $\rho(t)$ and the approximate MR state of Eq. (20) as a function of
time for (a)$\bar{n}=100$ and initial ground state for the MR; (b)
$\bar{n}=10$ and the mixture of zero and one phonon state as initial state of
the MR; (c) $\bar{n}=100$ and the mixture of zero and one phonon state as
initial state of the MR. The other parameters are given in the text and
coincide with those of Fig. 2.

###  IV.4 Cat state decoherence as decay of non-Gaussianity

The decoherence process affecting the MR state can also be described as a
dynamical “Gaussification” process in which the non-Gaussian even cat state
generated at short times by the engineered two-phonon reservoir becomes at
long times a convex mixture of Gaussian state, i.e., the equal-weight
incoherent superposition of the two coherent states $|\pm\beta\rangle$. This
suggests an alternative quantitative description of the above loss of quantum
coherence caused by the interplay between the engineered and natural reservoir
in terms of a measure of _quantum non-Gaussianity_ recently proposed in Refs.
Genoni  _et al._ (2013); Palma  _et al._ (2013).

A state is quantum non-Gaussian if it cannot be written as a convex sum of
Gaussian states, and a simple sufficient condition for non-Gaussianity can be
given in terms of the value of the Wigner function of the state at the phase
space origin $W[\rho](0)$ Genoni  _et al._ (2013): $\rho$ is quantum non-
Gaussian if $W[\rho](0)<(2/\pi)\exp[-2\langle n\rangle(\langle n\rangle+1)]$,
where $\langle n\rangle={\rm Tr}\left\\{\rho b^{\dagger}b\right\\}$ is the
mean number of excitations. However this condition does not detect many
quantum non-Gaussian states (for example even cat states) and a more efficient
condition for detecting quantum non-Gaussian states has been derived in Ref.
Palma  _et al._ (2013): $\rho$ is quantum non-Gaussian if there is a Gaussian
map ${\mathcal{E}}$ such that

$NG=W[{\mathcal{E}}(\rho)](0)-\frac{2}{\pi}\exp[-2\langle
n_{\mathcal{E}}\rangle(\langle n_{\mathcal{E}}\rangle+1)<0,$ (21)

where ${\mathcal{E}}(\rho)$ is the state transformed by the Gaussian map and
$n_{\mathcal{E}}$ is the mean excitation number of the transformed state.

Figure 7: Plot of the non-Gaussianity $NG$ of Eq. 21 versus time (soon after
the cat state generation) for (a) $\bar{n}=100$ starting from the mechanical
ground state, (b) $\bar{n}=10$ starting from the two-phonon cooling initial
state; (c) $\bar{n}=100$ starting from the two-phonon cooling initial state.
The other parameters are given in the text and coincide with those of Fig. 2.

We have calculated the quantity $NG$ quantifying non-Gaussianity by
restricting to Gaussian unitary maps formed by a composition of the phase
space displacement operator $D(\alpha)=\exp\left[\alpha
b^{\dagger}-\alpha^{*}b\right]$ and of the squeezing operator
$S(s)=\exp\left[(s/2)(b^{\dagger})^{2}-(s^{*}/2)b^{2}\right]$, and minimizing
$NG$ over $\alpha$ and $s$. The values $\alpha=0.35i$ and $s=0.01$ work very
well at all time instants after the cat state generation, either when starting
from the mechanical ground state and when starting from the mixture of the
vacuum and one phonon state obtained with two-phonon cooling. Plot of the time
evolution of $NG$ soon after the cat state generation, in the three cases
studied above, i.e., starting from the ground state and $\bar{n}=100$ (a),
starting from two-phonon cooling and $\bar{n}=10$ (b), and $\bar{n}=100$ (c),
are shown in Fig. 7. In all cases we see an exponential-like “decay” of non-
Gaussianity to the Gaussian limit $NG=0$, as expected, which is faster in the
cases when $\bar{n}=100$; the non-Gaussianity decay rate is in good agreement
with the usual decoherence rate
$2\gamma_{m}|\beta|^{2}\left(2\bar{n}+1\right)$. Therefore the measure of non-
Gaussianity of Eq. (21) proposed in Ref. Palma  _et al._ (2013) detects very
well the non-Gaussian property, and for the present even cat state the
dynamics of non-Gaussianity provides a satisfactory description of the
decoherence process.

##  V Conclusions

We have proposed a scheme for the deterministic generation of a linear
superposition of two coherent states of a MR based on the implementation of an
engineered reservoir realized by a bad cavity mode, bichromatically driven and
coupled _quadratically_ with the MR. The proposal extends in various aspects
the proposal of Ref. Tan  _et al._ (2013) and is feasible adopting either MIM
optomechanical setups or levitated nanospheres trapped around an intensity
maximum of the optical cavity mode. The interplay between the engineered
reservoir and the natural thermal reservoir of the MR allows the efficient
generation of the linear superposition state in a transient regime if the rate
of the engineered reservoir $\Gamma$ is larger than $\gamma_{m}\bar{n}$, which
is experimentally achievable in cryogenic environments at about $T\sim 10$ mK.
The generation of an even superposition of two coherent states of opposite
phases is almost ideal when starting from the MR ground state. This initial
condition could be obtained by laser pre-cooling the MR through a linear
optomechanical interaction, which however must be then suddenly switched to a
quadratic interaction, by shifting for example the membrane to a node of the
cavity mode. However the cat state generation is very efficient also when
precooling is realized by exploiting only the two-phonon relaxation processes
associated with the quadratic interaction Nunnenkamp  _et al._ (2010), which
is much easier to implement since it is based on the same configuration
allowing the cat state generation.

At longer times, the thermal reservoir is responsible for the progressive
decoherence of the generated superposition state, which asymptotically tends
to a steady state given by the incoherent mixture of the two coherent states
of the superposition, and which can be satisfactorily approximated by a simple
analytical expression. For this reason, the present protocol is ideal for
testing decoherence models acting on nanomechanical resonators.

An important issue is also the development of an efficient detection of the
generated MR state. A satisfactory detection could be obtained by realizing a
homodyne tomography D Ariano  _et al._ (1994) of the Wigner function of the
generated state. Homodyne tomography of the MR state could be obtained by
first transferring such state to an auxiliary cavity mode, weakly linearly
coupled to the MR, as suggested in Ref. Vitali  _et al._ (2007) or adopting
the pulsed homodyne measurement scheme of Ref. Vanner  _et al._ (2011). When
the auxiliary cavity mode is driven on its first red sideband and can be
adiabatically eliminated, its output field $a_{2}^{\rm out}$ is proportional
to the MR annihilation operator $b$ plus additional noise Vitali  _et al._
(2007), and therefore a calibrated homodyne detection of this output field at
various phases could be exploited for a tomographic reconstruction of the MR
Wigner function. The presence of the driven, weakly linearly coupled detection
cavity mode affects the two-phonon processes creating the engineered
reservoir, and therefore the detection process should be turned on only after
the cat state generation has been completed.

##  VI Acknowledgments

This work has been supported by the European Commission (ITN-Marie Curie
project cQOM), and by MIUR (PRIN 2011).

## References

  * Diehl  _et al._ (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. B chler, and P. Zoller, Nat. Phys., 4, 878 (2008).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. I. Cirac, Nat. Phys., 5, 633 (2009).
  * Wineland  _et al._ (1978) D. Wineland, R. Drullinger, and F. Walls, Phys. Rev. Lett., 40, 1639 (1978).
  * Poyatos  _et al._ (1996) J. Poyatos, J. Cirac, and P. Zoller, Phys. Rev. Lett., 77, 4728 (1996).
  * Carvalho  _et al._ (2001) A. R. R. Carvalho, P. Milman, R. L. de Matos Filho, and L. Davidovich, Phys. Rev. Lett., 86, 4988 (2001).
  * Syassen  _et al._ (2008) N. Syassen, D. Bauer, M. Lettner, T. Volz, D. Dietze, J. Garcia-Ripoll, J. Cirac, G. Rempe, and S. Dürr, Science, 320, 1329 (2008).
  * Krauter  _et al._ (2011) H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M. Petersen, J. I. Cirac, and E. S. Polzik, Phys. Rev. Lett., 107, 080503 (2011).
  * Pielawa  _et al._ (2007) S. Pielawa, G. Morigi, D. Vitali, and L. Davidovich, Phys. Rev. Lett., 98, 240401 (2007).
  * Pielawa  _et al._ (2010) S. Pielawa, L. Davidovich, D. Vitali, and G. Morigi, Phys. Rev. A, 81, 043802 (2010).
  * Wang and Clerk (2013) Y.-D. Wang and A. A. Clerk, Phys. Rev. Lett., 110, 253601 (2013).
  * Tan  _et al._ (2013) H. Tan, G. Li, and P. Meystre, Phys. Rev. A, 87, 033829 (2013a).
  * Bose  _et al._ (1999) S. Bose, K. Jacobs, and P. L. Knight, Phys. Rev. A, 59, 3204 (1999).
  * Marshall  _et al._ (2003) W. Marshall, C. Simon, R. Penrose, and D. Bouwmeester, Phys. Rev. Lett., 91, 130401 (2003).
  * Paternostro (2011) M. Paternostro, Phys. Rev. Lett., 106, 183601 (2011).
  * Vanner  _et al._ (2013) M. R. Vanner, M. Aspelmeyer, and M. S. Kim, Phys. Rev. Lett., 110, 010504 (2013).
  * Romero-Isart  _et al._ (2011) O. Romero-Isart, A. C. Pflanzer, F. Blaser, R. Kaltenbaek, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, Phys. Rev. Lett., 107, 020405 (2011a).
  * Romero-Isart  _et al._ (2011) O. Romero-Isart, A. C. Pflanzer, M. L. Juan, R. Quidant, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, Phys. Rev. A, 83, 013803 (2011b).
  * Jacobs  _et al._ (2009) K. Jacobs, L. Tian, and J. Finn, Phys. Rev. Lett., 102, 057208 (2009).
  * Jacobs  _et al._ (2011) K. Jacobs, J. Finn, and S. Vinjanampathy, Phys. rev. Lett., 83, 041801(R) (2011).
  * Tan  _et al._ (2013) H. Tan, F. Bariani, G. Li, and P. Meystre, Phys. Rev. A, 88, 023817 (2013b).
  * Deléglise  _et al._ (2008) S. Deléglise, I. Dotsenko, C. Sayrin, J. Bernu, M. Brune, J.-M. Raimond, and S. Haroche, Nature (London), 455, 510 (2008).
  * Myatt  _et al._ (2000) C. J. Myatt, B. E. King, Q. A. Turchette, C. A. Sackett, D. Kielpinski, W. M. Itano, C. Monroe, and D. J. Wineland, Nature (London), 403, 269 (2000).
  * Thompson  _et al._ (2008) J. D. Thompson, B. M. Zwickl, A. M. Jayich, F. Marquardt, S. M. Girvin, and J. G. E. Harris, Nature (London), 452, 72 (2008).
  * Sankey  _et al._ (2010) J. C. Sankey, C. Yang, B. M. Zwickl, A. M. Jayich, and J. G. E. Harris, Nat. Phys., 6, 707 (2010).
  * Karuza  _et al._ (2013) M. Karuza, M. Galassi, C. Biancofiore, C. Molinelli, R. Natali, P. Tombesi, G. Di Giuseppe, and D. Vitali, J. Opt., 15, 025704 (2013).
  * Barker (2010) P. F. Barker, Phys. Rev. Lett., 105, 073002 (2010).
  * Li  _et al._ (2011) T. Li, S. Kheifets, and M. G. Raizen, Nat. Phys., 7, 527 (2011).
  * Gieseler  _et al._ (2012) J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny, Phys. Rev. Lett., 109, 103603 (2012).
  * Kiesel  _et al._ (2013) N. Kiesel, F. Blaser, U. Delić, D. Grass, R. Kaltenbaek, and M. Aspelmeyer, arXiv:1304.6679v1 [quant-ph] (2013).
  * Gardiner and Zoller (2000) C. W. Gardiner and P. Zoller,  _Quantum Noise_ (Springer, 2000).
  * Wiseman and Milburn (1993) H. Wiseman and G. J. Milburn, Phys. Rev. A, 47, 642 (1993).
  * Wilson  _et al._ (2009) D. J. Wilson, C. A. Regal, S. B. Papp, and H. J. Kimble, Phys. Rev. Lett., 103, 207204 (2009).
  * Purdy  _et al._ (2012) T. P. Purdy, R. W. Peterson, P.-L. Yu, and C. A. Regal, New J. Phys., 14, 115021 (2012).
  * Flowers-Jacobs  _et al._ (2012) N. E. Flowers-Jacobs, S. W. Hoch, J. C. Sankey, A. Kashkanova, A. M. Jayich, C. Deutsch, J. Reichel, and J. G. E. Harris, Appl. Phys. Lett., 101, 221109 (2012).
  * Karuza  _et al._ (2012) M. Karuza, C. Molinelli, M. Galassi, C. Biancofiore, R. Natali, P. Tombesi, G. Di Giuseppe, and D. Vitali, New J. Phys., 14, 095015 (2012).
  * Purdy  _et al._ (2013) T. P. Purdy, R. W. Peterson, and C. A. Regal, Science, 339, 801 (2013).
  * Marquardt  _et al._ (2007) F. Marquardt, J. P. Chen, A. A. Clerk, and S. M. Girvin, Phys. Rev. Lett., 99, 093902 (2007).
  * Wilson-Rae  _et al._ (2007) I. Wilson-Rae, N. Nooshi, W. Zwerger, and T. Kippenberg, Phys. Rev. Lett., 99, 093901 (2007).
  * Genes  _et al._ (2008) C. Genes, D. Vitali, P. Tombesi, S. Gigan, and M. Aspelmeyer, Phys. Rev. A, 77, 033804 (2008).
  * Chan  _et al._ (2011) J. Chan, T. P. M. Alegre, A. H. Safavi-Naeini, J. T. Hill, A. Krause, S. Gröblacher, M. Aspelmeyer, and O. Painter, Nature (London), 478, 89 (2011).
  * Verhagen  _et al._ (2012) E. Verhagen, S. Deléglise, S. Weis, A. Schliesser, and T. J. Kippenberg, Nature (London), 482, 63 (2012).
  * Uhlmann (1976) A. Uhlmann, Rep. Math. Phys, 9, 273 (1976).
  * Jozsa (1994) R. Jozsa, J. Mod. Opt., 41, 2315 (1994).
  * Wang  _et al._ (2008) X. Wang, C.-S. Yu, and X. Yi, Phys. Lett. A, 373, 58 (2008).
  * Nunnenkamp  _et al._ (2010) A. Nunnenkamp, K. Børkje, J. G. E. Harris, and S. M. Girvin, Phys. Rev. A, 82, 021806(R) (2010).
  * Kennedy and Walls (1988) T. Kennedy and D. Walls, Phys. Rev. A, 37, 152 (1988).
  * Kim and Buzěk (1992) M. Kim and V. Buzěk, Phys. Rev. A, 46, 4239 (1992).
  * Brune  _et al._ (1996) M. Brune, E. Hagley, J. Dreyer, X. Maître, A. Maali, C. Wunderlich, J. M. Raimond, and S. Haroche, Phys. Rev. Lett., 77, 4887 (1996).
  * Genoni  _et al._ (2013) M. G. Genoni, M. Palma, T. Tufarelli, S. Olivares, M. Kim, and M. Paris, Phys. Rev. A, 87, 062104 (2013).
  * Palma  _et al._ (2013) M. L. Palma, J. Stammers, M. G. Genoni, T. Tufarelli, S. Olivares, M. S. Kim, and M. G. A. Paris, arXiv:1309.4221v1 (2013).
  * D Ariano  _et al._ (1994) G. D Ariano, C. Macchiavello, and M. Paris, Phys. Rev. A, 50, 4298 4302 (1994).
  * Vitali  _et al._ (2007) D. Vitali, S. Gigan, A. Ferreira, H. R. Bohm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Phys. Rev. Lett., 98, 030405 (2007).
  * Vanner  _et al._ (2011) M. R. Vanner, I. Pikovski, G. D. Cole, M. S. Kim, Č Brukner, K. Hammerer, G. J. Milburn, and M. Aspelmeyer, Proc. Nat. Acad. Sci, 108, 16182 (2011).