# Optomechanical state reconstruction and nonclassicality verification beyond
the resolved-sideband regime

Farid Shahandeh<EMAIL_ADDRESS>Centre for Quantum Computation and
Communication Technology, School of Mathematics and Physics, University of
Queensland, St Lucia, Queensland 4072, Australia  Department of Physics,
Swansea University, Singleton Park, Swansea SA2 8PP, United Kingdom  Martin
Ringbauer  Centre for Engineered Quantum Systems, School of Mathematics and
Physics, University of Queensland, St Lucia, Queensland 4072, Australia
Institute for Experimental Physics, University of Innsbruck, Technikerstraße
25, 6020 Innsbruck, Austria

(August 27, 2024)

###### Abstract

Quantum optomechanics uses optical means to generate and manipulate quantum
states of motion of mechanical resonators. This provides an intriguing
platform for the study of fundamental physics and the development of novel
quantum devices. Yet, the challenge of reconstructing and verifying the
quantum state of mechanical systems has remained a major roadblock in the
field. Here, we present a novel approach that allows for tomographic
reconstruction of the quantum state of a mechanical system without the need
for extremely high quality optical cavities. We show that, without relying on
the usual state transfer presumption between light an mechanics, the full
optomechanical Hamiltonian can be exploited to imprint mechanical tomograms on
a strong optical coherent pulse, which can then be read out using well-
established techniques. Furthermore, with only a small number of measurements,
our method can be used to witness nonclassical features of mechanical systems
without requiring full tomography. By relaxing the experimental requirements,
our technique thus opens a feasible route towards verifying the quantum state
of mechanical resonators and their nonclassical behaviour in a wide range of
optomechanical systems.

##  1 Introduction

Optomechanics [1] where a mechanical oscillator interacts with an optical
field via radiation pressure, is a promising direction of research for
fundamental physics [2, 3, 4] and the development of novel weak-force sensors
[5]. Yet, preparing a large mechanical device in a quantum state of motion,
and then verifying and exploiting the quantum nature of such a system have
remained elusive goals. Most research to date has focused on the problem of
state preparation and there have been significant recent advances using
measurement-based techniques to overcome the limitations of weak coupling and
non-zero initial thermal occupation [6, 7]. These techniques bring the
preparation of non-classical states of motion within the realm of current
experimental capabilities. Here, we thus focus on the problem of observing
quantum states of motion.

Continuous-variable quantum systems, such as optical fields or mechanical
resonators, are best described using distribution functions over a quantum
phase-space spanned by two quadratures of interest, such as position and
momentum, or amplitude and phase. However, one of the distinctive features of
quantum systems is that any such distribution must allow for negative values
and can thus not be a bona fide probability measure. The class of viable
distributions form a single-parameter family, called the _$s$ -parameterized
quasiprobability distributions_ [8]. Of particular interest within this family
is the _Wigner function_ , corresponding to a value of $s=0$, as the only
distribution that faithfully represents the quantum state and, at the same
time, correctly reproduces the marginal quadrature distributions. Another
important case is the _P-function_ , corresponding to $s=1$, whose negativity
is one of the main signatures of nonclassical behaviour [9].

In principle, the Wigner function of an unknown quantum state $\hat{\varrho}$
can be reconstructed from quadrature measurements via an inverse Radon
transform. This technique has been demonstrated for optical fields [10, 11],
where quadrature measurements are readily available using homodyne detection
[12, 13]. For mechanical resonators, however, direct quadrature measurements
are not available. Instead, the common approach aims to use a high-quality
optical cavity in the _resolved sideband_ regime, which would allow to
transfer the mechanical state onto an optical field, which can be
reconstructed using established techniques. This regime, however, is
experimentally very difficult to reach [14]. Techniques that work outside the
demanding resolved sideband regime, where state transfer in the usual sense is
not available, are thus highly desirable. One such scheme, using a classical
optical readout field, has been proposed and experimentally implemented
recently [15, 6]. Under realistic conditions, however, limitations from
optical shot noise of the readout field ultimately preclude an experimental
reconstruction of the mechanical Wigner function using this approach [14].
Consequently, Wigner-function tomography of a mechanical system remains an
outstanding challenge.

Here, we introduce an alternative approach which works outside the resolved
sideband regime with no, or only a low quality cavity. Our technique works by
imprinting the mechanical quadrature distributions on a strong optical pulse
in a coherent or squeezed state, from which it can be extracted using standard
optical methods. We show that this makes it possible to overcome the noise
limitations of previous approaches and perform a tomographic reconstruction of
any $s$-parameterized quasiprobability distribution, including the Wigner
function of the mechanical motion. Furthermore, since our technique is
independent of the single-photon coupling strength it is suitable for a wide
range of optomechanical systems with current technology. Finally, we discuss
how our approach could be used to detect P-function nonclassicality of a
mechanical system with only a small number of measurements. We also prove the
robustness of our nonclassicality criteria against readout noise and detection
inefficiency.

##  2 Results

This work is organised as follows. We first provide an overview of the
optomechanical interaction without assuming resolved mechanical sidebands. We
then discuss the description of this interaction over the quantum mechanical
phase-space. There we also derive the required transformations, and present
our readout scheme for obtaining the mechanical phase-space distributions as
the main result of this paper. We then compare our scheme with the existing
techniques in the non-resolved sideband regime. Finally, we provide a
necessary condition for demonstrating mechanical nonclassicality using the
smallest subset of data obtained from our scheme, namely, a single mechanical
tomogram. We further show that this condition is robust against noise and
detection inefficiency, and describe an explicit experimental protocol which
we show can be implemented with current technology in a variety of
optomechanical systems.

###  2.1 Optomechanical interaction

We now consider a mechanical resonator coupled to an optical field via the
radiation pressure interaction. A single photon reflecting off a mechanical
resonator imparts a momentum kick proportional to the photon frequency to the
resonator, while obtaining a phase-shift proportional to the mechanical
position. The evolution of the joint system is then described by the
Hamiltonian

$\frac{1}{\hbar}\hat{H}=\omega_{\rm o}a^{\dagger}a+\omega_{\rm
m}b^{\dagger}b+g_{0}a^{\dagger}a(b^{\dagger}+b),$ (1)

where $a$ ($b$) is the optical (mechanical) bosonic annihilation operator. The
terms of the Hamiltonian correspond, respectively, to the free evolution of
the optical field with frequency $\omega_{\rm o}$, the free evolution of the
mechanics with frequency $\omega_{\rm m}$, and the radiation pressure
interaction. The latter couples the mechanical position operator
$\hat{X}\propto b^{\dagger}+b$ to the optical photon number operator
$\hat{N}=a^{\dagger}a$ with the single photon coupling strength $g_{0}=Gx_{\rm
zpf}$, where $G$ is the coupling constant, and $x_{\rm
zpf}=\sqrt{\hbar/2m_{\rm eff}\omega_{\rm m}}$ the mechanical zero-point
fluctuation. In this picture the mechanical response to a single photon is
captured by $g_{0}$, and can be linearly enhanced by a cavity that allows for
multiple interactions per photon [1].

In the cavity-enhanced case one can identify a number of interesting parameter
regimes. Much work to date has focused on the so-called _resolved sideband
regime_ , where the decay rate of the optical cavity $\kappa_{\text{o}}$ is
much smaller than the mechanical frequency. This makes it possible to excite,
to good approximation, only one of the two motional sidebands of the cavity
resonance. Ignoring non-resonant terms, which is known as the rotating wave
approximation, the optomechanical interaction is then dominated, either by a
beam-splitter-like interaction for the red sideband, or a squeezing-like
interaction for the blue sideband [16]. For tomography, however, this regime
is challenging due to the incompatible requirements of a high optical quality
factor for sideband resolution and a low optical quality factor to retain
control over the optical amplitude fluctuations inside the cavity.

Figure 1: Sketch of our tomography method. (a) Schematic experimental setup.
The initial (squeezed) vacuum state $\mathscr{W}_{\rm o}$ is first displaced
by $\alpha$ via interference with a strong laser pulse of length $\tau$
(dashed line) on a highly transmissive beamsplitter. The probe beam then
interacts with the mechanical system, is subsequently displaced by $-\alpha$,
and finally measured using homodyne detection. Both displacements are easily
implemented in an inherently stable manner using the same laser pulse, see
Secs. 2.4.1 and 2.5 for details. (b) In the phase-space representation (i) the
mechanical resonator is initially in an unknown state $\mathscr{W}_{\rm m}$
and the optical system is in the state $\mathscr{W}_{\rm o}$. (ii) The two
systems interact under the full optomechanical Hamiltonian, in the case of a
strong optical displacement given by Eqs. (2)–(4). (iii) After the
interaction, the $s$-parameterized $\phi$-tomogram of the mechanical system
will be imprinted on the momentum-quadrature of the optical Wigner-function,
which can be obtained from a standard homodyne measurement. Using a number of
such optical/mechanical tomograms one can use the inverse Radon transformation
in Fig. 3 to reconstruct the mechanical $s$-parameterized quasiprobability
distribution. For illustrative purposes only we chose an initial mechanical
state of the form
$\left(|0\rangle+2|1\rangle+|2\rangle+2|3\rangle\right)/\sqrt{10}$ in the
Fock-state basis. The displayed tomograms include points where we numerically
added $10\%$ noise to indicate measurement imprecision.

Here, we are interested in the technically simpler scenario in which the
optical field is a free-space mode, or weakly enhanced by a low-quanlity
factor optical cavity (with a cavity decay rate much faster than the
mechanical frequncy), see Fig. 1a. In this regime, the motional sidebands
cannot be resolved, implying that the rotating-wave approximation is not valid
and the full optomechanical Hamiltonian in Eq. (1) must be considered. Suppose
that the mechanical system is initially in a state $\hat{\varrho}_{\rm m}$
which then interacts with a displaced optical pulse,
$\hat{D}(\alpha)\hat{\varrho}_{\rm o}\hat{D}^{\dagger}(\alpha)$ over a
duration $\tau$. At a time $\tau_{0}$ after the interaction, the optical field
is displaced by $-\alpha$. The evolution of the state is thus given by
$\hat{D}(-\alpha)\exp{-i\tau_{0}(a^{\dagger}a+b^{\dagger}b)}\exp\\{-i\tau[\omega_{\rm
o}a^{\dagger}a+\omega_{\rm
m}b^{\dagger}b+g_{0}a^{\dagger}a(b^{\dagger}+b)]\\}\hat{D}(\alpha)$. As we
show in detail in Appendix A, under the condition that the optical field is
strong (i.e. $|\alpha|\gg 1$) and using the fact that we only perform
measurements on the optical field, the evolution of the system can be
disentangled and linearised to $\hat{\mathcal{U}}_{\rm o}\leavevmode\nobreak\
\hat{\mathcal{U}}_{\rm k}\leavevmode\nobreak\ \hat{\mathcal{U}}_{\rm om}$,
where $\hat{\mathcal{U}}_{\rm om}$, $\hat{\mathcal{U}}_{\rm k}$, and
$\hat{\mathcal{U}}_{\rm o}$ are unitaries representing the optomechanical
interaction, an optical Kerr interaction, and a displaced optical free
evolution, respectively, given by

$\displaystyle\hat{\mathcal{U}}_{\rm om}=e^{\chi
u(\tilde{a}^{\dagger}_{\theta}+\tilde{a}_{\theta})(\tilde{b}^{\dagger}_{\varphi}-\tilde{b}_{\varphi})},$
(2) $\displaystyle\hat{\mathcal{U}}_{\rm
k}=e^{iv[4a^{\dagger}a+2r(e^{i\theta}a^{\dagger}+e^{-i\theta}a)+(e^{2i\theta}a^{{\dagger}2}+e^{-2i\theta}a^{2})]},$
(3) $\displaystyle\hat{\mathcal{U}}_{\rm
o}=\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a}\hat{D}(\alpha).$ (4)

Here, we have defined the parameters as follows: $\alpha=re^{i\theta}$ with
$r\in\mathbb{R}^{+}$, $\mu=1-e^{i\omega_{\rm m}\tau}=ue^{i\varphi}$ with

$u{=}\sqrt{2(1{-}\cos{\omega_{\rm
m}\tau})},\quad\varphi{=}\tan^{-1}{\Big{(}\frac{\sin{\omega_{\rm
m}\tau}}{1-\cos{\omega_{\rm m}\tau}}\Big{)}},$ (5)

and

$\displaystyle\chi=\frac{g_{0}r}{\omega_{\rm m}},$ (6)
$\displaystyle\tilde{a}_{\theta}=ae^{-i\theta}+\frac{r}{2},\quad\tilde{b}_{\varphi}=be^{-i\varphi}.$
(7)

It is important to note that, while the phase of the optical mode $\theta$ is
being set by the input optical coherent state, the relevant mechanical phase
$\varphi$ is determined by the length of the input optical pulse $\tau$.

###  2.2 Optomechanical interaction in phase space

####  2.2.1 Mode transformation

By defining a vector of input mode operators
$\hat{\vec{A}}:=(a^{\dagger},b^{\dagger},a,b)$, we aim to find the the output
mode operators in the Heisenberg picture under the aforementioned unitary
evolution, i.e. $\hat{\vec{A}}^{\prime}=(\hat{\mathcal{U}}_{\rm
o}\hat{\mathcal{U}}_{\rm k}\hat{\mathcal{U}}_{\rm
om})\hat{\vec{A}}(\hat{\mathcal{U}}^{\dagger}_{\rm
om}\hat{\mathcal{U}}^{\dagger}_{\rm k}\hat{\mathcal{U}}^{\dagger}_{\rm o})$.
Since these unitary evolutions are single-mode and two-mode linear
transformations, they correspond to elements of the symplectic groups ${\rm
Sp}(2,\mathbb{R})$ and ${\rm Sp}(4,\mathbb{R})$, respectively. We can thus
exploit the properties of symplectic transformations, in particular, that if
$\hat{\mathcal{U}}=\exp\\{\frac{1}{2}\hat{\vec{A}}\ln(M)\Sigma\hat{\vec{A}}^{\mathsf{T}}\\}$,
with $\Sigma={\rm antidiag}\\{I,-I\\}$, then
$\hat{\mathcal{U}}\hat{\vec{A}}\hat{\mathcal{U}}^{\dagger}=\hat{\vec{A}}M$
[17, 18]. Using Eqs. (2), (3), and (4), after some algebra and making
reasonable restrictions on the choice of the optical pulse $\tau$ and the
delay duration $\tau_{0}$, via

$\begin{matrix}[l]\omega_{\rm m}\tau-\sin(\omega_{\rm
m}\tau)=k\pi/2\sqrt{3}\chi^{2},&\quad k\in 2\mathbb{Z}^{+},\\\
(\tau_{0}+\tau)\omega_{\rm o}=2\pi m,&\quad m\in\mathbb{N},\end{matrix}$ (8)

the transformed optomechanical mode operators are found as

$\begin{split}&a^{\prime}=a-i\sqrt{2}\chi u\hat{X}(\varphi+\frac{\pi}{2}),\\\
&b^{\prime}=b-\sqrt{2}\chi ue^{i\varphi}[\hat{x}+r].\end{split}$ (9)

Here, $\hat{X}(\varphi)=(b^{\dagger}e^{i\varphi}+be^{-i\varphi})/\sqrt{2}$ and
$\hat{x}=(a^{\dagger}+a)/\sqrt{2}$ are the $\varphi$-quadrature of the
mechanics and the $x$-quadrature of the optical field, respectively. We refer
the interested reader to Appendix B for the detailed derivation of the general
mode transformation, as well as the simplified Eqs. (8) and (9). Importantly,
in Eq. (8) we require that $0\ll\omega_{\rm m}\tau\ll 2\pi$, both to obtain a
nonzero coupling ($u\neq 0$) and to minimize the effect of mechanical losses.
This can be achieved by assuming the pulse duration $\tau$ to be on the order
of half a mechanical period, i.e., $\omega_{\rm m}\tau\approx\pi$. It is
evident from Eq. (9) that the output optical field now carries information
about the mechanical quadratures. At the same time, the mechanical mode is
affected by the back-action noise proportional to the input light quadrature.
Before we can proceed to the Wigner-function picture, we need the following
result, the proof of which is given in the Appendix C.

######  Lemma 1.

For the order parameter $s=0$, a product of Weyl-Wigner operators remains a
product under all linear transformations of the mode operators,
$\begin{pmatrix}a^{\prime{\dagger}}&b^{\prime{\dagger}}&a^{\prime}&b^{\prime}\end{pmatrix}=[\begin{pmatrix}a^{\dagger}&b^{\dagger}&a&b\end{pmatrix}S+\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm a}\end{pmatrix}]$, the argument of which
changes according to the inverse of the corresponding symplectic
transformation, i.e.,
$\hat{T}(\alpha;0)\otimes\hat{T}(\beta;0)=\hat{T}(\alpha^{\prime};0)\otimes\hat{T}(\beta^{\prime};0)$
where
$\begin{pmatrix}\alpha^{\prime*}&\beta^{\prime*}&\alpha^{\prime}&\beta^{\prime}\end{pmatrix}=[\begin{pmatrix}\alpha^{*}&\beta^{*}&\alpha&\beta\end{pmatrix}-\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm a}\end{pmatrix}]S^{-1}$.

####  2.2.2 Interaction in phase space

We are now ready to examine the optomechanical interactions in the phase-space
representation, see Fig. 1b. The Wigner decomposition (see Sec. 4.1) of the
joint initial optomechanical state is

$\hat{\varrho}_{\rm in}=\int\frac{d^{2}\alpha
d^{2}\beta}{\pi^{2}}\mathscr{W}_{\rm o}(\alpha;0)\mathscr{W}_{\rm
m}(\beta;0)\hat{T}_{\rm o}(\alpha;0)\otimes\hat{T}_{\rm m}(\beta;0),$

where $\mathscr{W}_{\rm o}$ and $\mathscr{W}_{\rm m}$ are the Wigner functions
of light and mechanics, respectively. Allowing for a free mechanical evolution
for a time $\tau_{\rm d}$ before the interaction with the optical field, the
output state is given by $\hat{\varrho}_{\rm out}=[\hat{\mathcal{U}}_{\rm
o}\hat{\mathcal{U}}_{\rm k}\hat{\mathcal{U}}_{\rm om}\hat{\mathcal{R}}_{\rm
m}(\varphi_{\rm d})]\hat{\varrho}_{\rm in}[\hat{\mathcal{R}}^{\dagger}_{\rm
m}(\varphi_{\rm d})\hat{\mathcal{U}}^{\dagger}_{\rm
om}\hat{\mathcal{U}}^{\dagger}_{\rm k}\hat{\mathcal{U}}^{\dagger}_{\rm o}]$
where $\hat{\mathcal{R}}_{\rm m}(\varphi_{\rm d})=\exp\\{-i\varphi_{\rm
d}b^{\dagger}b\\}$ with $\varphi_{\rm d}=\omega_{\rm m}\tau_{\rm d}$. This
results in $\hat{\varrho}_{\rm out}=\int\frac{d^{2}\alpha
d^{2}\beta}{\pi^{2}}\mathscr{W}_{\rm o}(\alpha;0)\mathscr{W}_{\rm
m}(\beta;0)\hat{T}_{\rm o}(\alpha^{\prime};0)\otimes\hat{T}_{\rm
m}(\beta^{\prime};0)$. Using Lemma 1 and the inverse of the symplectic
transformation corresponding to the overall evolution $\hat{\mathcal{U}}_{\rm
o}\hat{\mathcal{U}}_{\rm k}\hat{\mathcal{U}}_{\rm om}$ given by Eq. (9), we
obtain

$\begin{split}&\alpha^{\prime}=\alpha+i\sqrt{2}\chi uX(\varphi+\varphi_{\rm
d}+\frac{\pi}{2}),\\\ &\beta^{\prime}=\beta e^{-i\varphi_{\rm d}}+\sqrt{2}\chi
ue^{i\varphi}[x(\theta)+r].\end{split}$

Since we only perform measurements on the optical field, we can trace out the
mechanics and make use of the normality of the Weyl-Wigner operators to get
$\hat{\varrho}_{\rm out;o}=\int\frac{d^{2}\alpha
d^{2}\beta}{\pi^{2}}\mathscr{W}_{\rm o}(\alpha;0)\mathscr{W}_{\rm
m}(\beta;0)\hat{T}_{\rm o}(\alpha^{\prime};0)$. Using the orthogonality
relation of the Weyl-Wigner operators, the final measured Wigner-function of
the optical field can be readily obtained as $\mathscr{W}_{\rm
out;o}(\alpha^{\prime};0)=\int\frac{d^{2}\beta}{\pi}\mathscr{W}_{\rm
o}[\alpha^{\prime}-i\sqrt{2}\chi uX(\varphi+\varphi_{\rm
d}+\pi/2);0]\mathscr{W}_{\rm m}(\beta;0)$. Switching to position and momentum
coordinates $(x^{\prime},p^{\prime})$ of the optical output, we obtain

$\begin{split}&\mathscr{W}_{\rm out;o}(x^{\prime},p^{\prime};0)=\\\
&\int\frac{d^{2}\beta}{\pi}\mathscr{W}_{\rm o}[x^{\prime},p^{\prime}-2\chi
uX(\varphi+\varphi_{\rm d}+\frac{\pi}{2});0]\mathscr{W}_{\rm
m}(\beta;0),\end{split}$ (10)

showing that only the momentum quadrature of the optical output is modified by
the interaction.

###  2.3 Mechanical tomography

Having established the tools to investigate the information about the
mechanical quantum state imprinted in the refelcted probe ligh, we now
describe how to use these tools to achieve quantum tomography of the
mechanical state as depicted in Fig. 1b.

####  2.3.1 General case readout

In principle, an optical vacuum state would be sufficient as a probe pulse for
our method. However, in the general case we now consider a momentum-squeezed
vacuum state with Wigner function

$\mathscr{W}_{\rm
o}(x,p;0)=\frac{1}{\pi}\exp\\{-(e^{-2\varepsilon}x^{2}+e^{2\varepsilon}p^{2})\\},$
(11)

where $\varepsilon\in[0,\infty)$ is the squeezing parameter. While this is not
strictly necessary for our protocol, it adds an additional parameter that
could be convenient for some experimental implementations. Using this initial
state in Eq. (10) and changing the variable $p^{\prime}\rightarrow 2\chi
up^{\prime}$, we obtain the Wigner function of the optical output as

$\begin{split}&\mathscr{W}_{\rm out;o}(x^{\prime},2\chi
up^{\prime};0)=\frac{1}{\pi}e^{-e^{-2\varepsilon}x^{\prime 2}}\\\
&\quad\times\int\frac{d^{2}\beta}{\pi}e^{-4\chi^{2}u^{2}e^{2\varepsilon}[p^{\prime}-X(\varphi+\varphi_{\rm
d}+\frac{\pi}{2})]^{2}}\mathscr{W}_{\rm m}(\beta;0).\end{split}$ (12)

From Eq. (12) it is evident that the information about the mechanical
$(\varphi+\varphi_{\rm d}+\pi/2)$-quadrature is encoded within the momentum
quadrature of the optical field. Hence, we can integrate over the position of
the optical Wigner function and arrive at

$\begin{split}&w_{\rm out;o}(2\chi up^{\prime},\frac{\pi}{2})=\\\
&\qquad\frac{e^{\varepsilon}}{\sqrt{\pi}}\int\frac{d^{2}\beta}{\pi}e^{-4\chi^{2}u^{2}e^{2\varepsilon}[p^{\prime}-X(\varphi+\varphi_{\rm
d}+\frac{\pi}{2})]^{2}}\mathscr{W}_{\rm m}(\beta;0),\end{split}$ (13)

in which $w_{\rm out;o}(p,\frac{\pi}{2})$ is the momentum tomogram of the
optical output pulse. Interestingly, we can identify the Gaussian kernel in
Eq. (13) with the $s$-parameterized quasiprobability of the tomographic
operator (see Sec 4.2), obtaining

$\begin{split}&e^{-4\chi^{2}u^{2}e^{2\varepsilon}[p^{\prime}-X(\varphi+\varphi_{\rm
d}+\frac{\pi}{2})]^{2}}=\\\
&\qquad\qquad\qquad\sqrt{s^{\star}\pi}\leavevmode\nobreak\
\mathscr{W}_{\hat{\Pi}(p^{\prime},\varphi+\varphi_{\rm
d}+\frac{\pi}{2})}(\beta;-s^{\star}),\end{split}$ (14)

with order parameter

$s^{\star}=(2\chi ue^{\varepsilon})^{-2}.$ (15)

As we show in Appendix D, the order mismatch can be moved to the mechanical
state, which leads to

$\begin{split}&\int\frac{d^{2}\beta}{\pi}e^{-4\chi^{2}u^{2}e^{2\varepsilon}[p^{\prime}-X(\varphi+\varphi_{\rm
d}+\frac{\pi}{2})]^{2}}\mathscr{W}_{\rm m}(\beta;0)=\\\
&\sqrt{s^{\star}\pi}\int\frac{d^{2}\beta}{\pi}\mathscr{W}_{\hat{\Pi}(p^{\prime},\varphi+\varphi_{\rm
d}+\frac{\pi}{2})}(\beta;-s^{\star})\mathscr{W}_{\rm m}(\beta;0)=\\\
&\sqrt{s^{\star}\pi}\int\frac{d^{2}\beta}{\pi}\mathscr{W}_{\hat{\Pi}(p^{\prime},\varphi+\varphi_{\rm
d}+\frac{\pi}{2})}(\beta;0)\mathscr{W}_{\rm m}(\beta;-s^{\star})=\\\
&\qquad\sqrt{s^{\star}\pi}\leavevmode\nobreak\ w_{\rm
m}(p^{\prime},\varphi+\varphi_{\rm d}+\frac{\pi}{2};-s^{\star}).\end{split}$
(16)

Here, we have defined the $s$-parameterized mechanical tomogram $w_{\rm
m}(x,\varphi;s)$ to be the marginal of the $s$-parameterized quasiprobability
distribution of the mechanical quantum state. Importantly, the relation
obtained in (16) represents a Radon transformation (see Sec. 4.2) of the
$(-s^{\star})$-parameterized quasiprobability distribution $\mathscr{W}_{\rm
m}(\beta;-s^{\star})$.

$\begin{split}&\mathscr{W}_{\rm
m}(\beta;-s^{\star})\quad\xleftrightarrow[\text{Inv. Radon
transform}]{\text{Radon transform}}\quad\\{w_{\rm
m}(x,\phi;-s^{\star})\\},\end{split}$ (17)

Substituting Eq. (16) into Eq. (13) we find

$\begin{split}w_{\rm out;o}(2\chi up^{\prime},\frac{\pi}{2};0)=\frac{1}{2\chi
u}w_{\rm m}(p^{\prime},\varphi+\varphi_{\rm
d}+\frac{\pi}{2};-s^{\star}),\end{split}$ (18)

where $s^{\star}$ is given by Eq. (15).

Equation (18) is the main result of this manuscript, implying that the
$(-s^{\star})$-parameterized mechanical tomogram is encoded, up to a scaling
factor in the momentum quadrature of the optical Wigner function, which can be
easily read out using well-established techniques. Note, however, that due to
the scaling one might require much larger quadrature detection-range in the
optical homodyne measurement to access the optical tomograms over the required
range of values. This can be avoided by setting $2\chi u=1$ and instead
exploiting the momentum squeezing of the initial optical state to set the
order parameter $s^{\star}$. We will discuss in detail the experimental
feasibility of all such assumptions in our experimental proposal.

####  2.3.2 Strong interaction regime

Recalling Eq. (15), it turns out that the order parameter $s^{\star}$ goes to
zero quadratically with the displacement amplitude $r$ and exponentially with
the squeezing parameter $\varepsilon$. Hence, for sufficiently large $\chi
ue^{\varepsilon}$ we can approximate

$\begin{split}w_{\rm m}(p^{\prime},\varphi+\varphi_{\rm
d}+\frac{\pi}{2};-s^{\star})&\approx w_{\rm m}(p^{\prime},\varphi+\varphi_{\rm
d}+\frac{\pi}{2};0)\end{split}$ (19)

which in turn, implies that the $\phi$-tomogram of the mechanical _Wigner
function_ is imprinted in a rescaled form on the momentum quadrature of the
optical output Wigner function,

$\begin{split}w_{\rm out;o}(2\chi up^{\prime},\frac{\pi}{2};0)=\frac{1}{2\chi
u}w_{\rm m}(p^{\prime},\varphi+\varphi_{\rm d}+\frac{\pi}{2};0).\end{split}$
(20)

####  2.3.3 Advantages over standard pulsed tomography

The use of short pulses of coherent light for optomechanical state tomography
outside the resolved sideband regime was first proposed in Ref. [15]. Their
approach works by imprinting the mechanical quadrature distribution onto the
momentum quadrature of the probe pulse [14], which has recently been
demonstrated experimentally [6]. In practice, however, this simple technique
is limited by the fundamental quantum noise on the readout pulse. This
precludes experimental access to the Wigner function of the mechanical system
and indeed might not allow for the reconstruction of any quasiprobability
distribution with $s>-1$, as we discuss below.

The approach of Ref. [15, 6] uses classical readout pulses, which have a non-
negative quasiprobability distribution for $s=1$. The bosonic coherent state
basis in this case is given by
$\hat{T}(\alpha;-1)=|\alpha\rangle\langle\alpha|$, and thus

$\hat{\varrho}=\int\frac{d^{2}\alpha}{\pi}P(\alpha)|\alpha\rangle\langle\alpha|.$
(21)

The distribution corresponding to the optical field is then a Gaussian of the
form

$P(\alpha)=\frac{1}{2\pi\sigma_{x}\sigma_{p}}\
e^{-\frac{(x-\bar{x})^{2}}{2\sigma^{2}_{x}}-\frac{(p-\bar{p})^{2}}{2\sigma^{2}_{p}}},$
(22)

representing a (possibly squeezed) thermal state of light with
$x=\sqrt{2}\,{\rm Re}\alpha$ and $p=\sqrt{2}\,{\rm Im}\alpha$. The function
$P(\alpha)$ is non-negative and corresponds to a class of quantum states
called P-classical states. The method of Ref. [14] can now be recast in our
notation as

$\begin{split}w_{\rm
out;o}(p,\frac{\pi}{2};0)=\int\frac{dx}{\sqrt{\pi}}e^{-\frac{(p-\eta
x)^{2}}{1+2\sigma_{p}^{2}}}w_{\rm m}(x,\varphi;0)\end{split}$ (23)

or equivalently,

$w_{\rm out;o}(\eta
p,\frac{\pi}{2};0)=\int\frac{dx}{\sqrt{\pi}}e^{-\frac{\eta^{2}(p-x)^{2}}{1+2\sigma_{p}^{2}}}w_{\rm
m}(x,\varphi;0),$ (24)

where $\eta$ is the strength of the mechanical position measurement. Despite
similarities between Eq. (24) and our Eq. (13), there are important
differences providing an advantage to our method over that of Ref. [15, 14].
First, by assuming the input to be P-classical, one imposes the restriction
that the variance of the momentum quadrature of the readout is shot-noise
limited, i.e. $\sigma_{p}$ in Eq. (22) and hence in Eq. (24) is always non-
negative. Thus, the sharpest mechanical quadrature measurement possible within
that scheme is achieved in the limit $\sigma_{p}\to 0$, corresponding to an
optical coherent state readout. Consequently, nonclassical states, in
particular those with negative P-functions that possess ill-defined P-function
variances, e.g. squeezed coherent states, are not considered in this approach.
In contrast, our method makes no assumption about the classicality of the
readout optical field to arrive at Eqs. (13) and (18), allowing us to use a
nonclassical readout such as a squeezed coherent state of light. This makes it
possible to use an optical readout field with a momentum variance
significantly below the shot-noise limit.

To extract the mechanical Wigner tomograms $w_{\rm m;o}(x,\varphi;0)$ from the
output optical tomograms $w_{\rm out;o}(p,\frac{\pi}{2};0)$ using Ref. [15,
14], the Gaussian convolution kernel
$\exp\\{-\eta^{2}(p-x)^{2}/(1+2\sigma_{p}^{2})\\}$ on the r.h.s. of Eq. (24)
must be as narrow as possible to be approximated by a Dirac delta function.
Using the best classical readout with $\sigma_{p}=0$, there are two
possibilities to achieve this goal.

  1. (i)

Increasing the strength of the readout pulses such that $\eta\gg 1$ [6]. This
route, however, is problematic, since the scaling on the l.h.s. of Eq. (24)
implies that increasing $\eta$ also requires larger sensitivity in the optical
homodyne measurement to access the optical tomograms over the required range
of values. Hence, it is desirable to work in a regime where $\eta\approx 1$,
leaving us with the second possibility below.

  2. (ii)

Inverting the convolution on the r.h.s. of Eq. (24) by applying a
deconvolution map to the output optical tomogram $w_{\rm out;o}(\eta
p,\frac{\pi}{2};0)$ with the kernel $e^{\eta^{2}(p-x)^{2}}$. In theory this
allows one to compute the mechanical tomogram $w_{\rm m}(x,\varphi;0)$
exactly. In practice, however, the numerical process becomes very unstable
when applied to real experimental data [14], since the exponentially rising
kernel strongly amplifies any noise [19].

In contrast, our approach opens up a third possibility to overcome the
smoothing problem above, namely by using a _nonclassical_ readout pulse. A
careful comparison between Eq. (24) and our Eq. (13) shows that the extra free
parameter $\varepsilon$ introduced into our model via squeezing of the initial
optical probe neither acts as an scaling factor nor puts any strict lower
bounds on the variance of the momentum quadrature of the readout.
Consequently, using our approach one can reliably set the scaling parameter
$2\chi u\approx 1$, and manipulate the momentum squeezing to achieve, in
principle, an ordering parameter $s^{\star}$ arbitrarily close to zero. Our
approach thus allows for an enhanced mechanical quadrature measurement without
the previous difficulties, thereby overcoming one of the most challenging
problems in optomechanical tomography through the use of momentum-squeezed
probe light.

###  2.4 Certification of mechanical nonclassicality

The method introduced above allows for the full reconstruction of
$s$-parametrized quasiprobability distributions of a mechanical system in a
new experimental regime that is technically easier than previous approaches.
However, like every tomography method, it requires a large number of
tomograms. Below, we show in detail that, if the goal is not full mechanical
tomography, but merely demonstrating mechanical nonclassicality, this can be
done much more simply using a single $s$-parameterized mechanical tomogram.

In the following we focus on the nonclassicality signified by negativities of
the P-function, i.e. the quasiprobability distribution for $s=1$, as given in
Eq. (21). This is a favourable choice over the Wigner function, since
P-negativity is more resilient to noise and loss, and captures some non-
classical states, such as squeezed vacuum states, which have a positive Wigner
function. Furthermore, P-function negativity can, in principle, be directly
observed using appropriate filtering of the homodyne data [20, 21, 22], or
_witnessed_ without full tomography, using only a small number of observables
[23, 24, 25, 26, 27].

Figure 2: Witnessing P-function non-classicality. Using the method outlined in
the main text (c.f. Fig. 1) one can obtain a $s$-tomogram of the state of the
mechanical resonator. Using this $s$-tomogram, one computes (a) the fictitious
$s$-parameterized quasiprobability distribution $\mathscr{W}_{\rm
f}(x,p;s)=w_{\rm m}(x;s)w_{0}(p;s)$, where $w_{0}(p;s)$ is an $s$-tomogram of
the vacuum state, or (b) the fictitious distribution $\mathscr{W}_{\rm
f}(x,p;s)=w_{\rm m}(x;s)w_{\rm m}(p;s)$, using the same tomogram for both
quadratures. If either of the the resulting functions does not correspond to a
bona-fide density operator, then the state of the mechanical resonator has
been P-nonclassical.

Unfortunately, these methods do not translate directly to our optomechanical
scenario, where we only have access to a restricted set of observables,
namely, tomograms. Nonetheless, building on a criterion introduced by Park et
al [28], a single mechanical tomogram can be used to witness nonclassicality
as follows. Given a Wigner tomogram $w_{\rm m}(x)$ of a quantum state for some
quadrature angle $\phi$ (here, for brevity, we have dropped the phase from the
argument of the tomograms), one can construct a ficticious Wigner function (in
the position-momentum basis) using the so-called _first demarginalization_
maps $w_{\rm m}(x)\mapsto\mathscr{W}_{\rm f}(x,p)=w_{\rm m}(x)w_{0}(p)$, where
$w_{0}(p)=\sqrt{2}\exp\\{-p^{2}\\}$ is a marginal of the vacuum state. If the
fictitious Wigner function $\mathscr{W}_{\rm f}(x,p)$ is not a legitimate
quantum state, then the state of the mechanical resonator has been
P-nonclassical, see Fig 2a.

In case the first demarginalization cannot certify nonclassicality and the
results are inconclusive, one can consider the _second_ demarginalization map
which uses a duplication of a single tomogram by changing $x\rightarrow p$ in
$w_{\rm m}(x)$, to construct a ficticious Wigner functions as $(w_{\rm
m}(x),w_{\rm m}(x))\mapsto\mathscr{W}_{\rm f}(x,p)=w_{\rm m}(x)w_{\rm m}(p)$.
Similarly, the failure of $\mathscr{W}_{\rm f}(x,p)$ to represent a legitimate
quantum state implies the P-nonclassicality of the state of the mechanical
resonator, see Fig 2b. Importantly, testing whether the result is a legitimate
Wigner function can be relatively simple in many cases of interest [28]. These
simple, necessary but not sufficient, criteria are thereby able to detect the
nonclassicality of a large class of quantum states, including phonon added
thermal states, which are of particular current interest.

As we show in detail in the Appendix E, this procedure can be directly
extended to any $s$-parametrized quasiprobability distribution. Given an
$s$-parameterized tomogram of a quantum state, $w_{\rm m}(x;s)$, one can
define the following two demarginalization maps,

$\displaystyle(w_{\rm m}(x;s),w_{0}(x;s))\mapsto\mathscr{W}_{\rm
f}(x,p;s)=w_{\rm m}(x;s)w_{0}(p;s),$ (25) $\displaystyle(w_{\rm m}(x;s),w_{\rm
m}(x;s))\mapsto\mathscr{W}_{\rm f}(x,p;s)=w_{\rm m}(x;s)w_{\rm m}(p;s),$ (26)

in which $w_{0}(x;s)$ is the $s$-parameterized tomogram of the vacuum state
and $\mathscr{W}_{\rm f}(x,p;s)$ is the $s$-parameterized quasiprobability
distribution of a fictitious quantum state. The failure of either of the
distributions in Eq. (25) or (26) to represent a bona fide quantum state
implies the P-function nonclassicality of the corresponding quantum state. In
other words, the two maps presented above provide necessary nonclassicality
criteria.

To verify the legitimacy of the quantum state corresponding to the above
fictitious quasiprobability distributions, one can simply employ all the
arguments provided in Ref. [28] by using appropriate quasiprobability
distributions to compute the matrix elements of the fictitious density
operator using the phase-space trace relation. This shows that our method can
be used to verify P-function nonclassicality even in the regime of weak
interactions.

####  2.4.1 Effect of imperfections

An important question that remains is whether experimental noise and losses
could lead to a false positive result in our proposed nonclassicality-
verification procedure. Since any noise on the mechanical resonator can be
considered part of the mechanical quantum state, we only have to study the
effect of noise and losses on the optical readout signal.

Considering first the readout channel before the optomechanical interaction,
note that the optical input in Eq. (11) has a Gaussian Wigner function. Since
displacement operations as well as all experimental noise and losses are
represented by Gaussian channels [29], the state always remains Gaussian and
can easily be fully characterised; see for example [30]. One can then use this
experimentally characterized state in Eq. (10) and proceed as before.
Considering now the noise and losses after the interaction, we show in
Appendix F that these commute with the interaction and can thus be treated in
the same way.

Given that these imperfections are well characterized, they can then be taken
into account by correcting the parameter $s$ of the fictitious tomogram
$w_{0}(p;s)$ in Eq. (25). Curiously, this is not necessary for our second
nonclassicality criterion, Eq. (26), which is immune to noise and losses by
virtue of using the same measured tomogram for both quadratures. The latter
implies that both tomograms have matching order parameters. Nonetheless, since
the legitimacy tests in Appendix E, require using the $s$-parameterized phase-
space trace rule, a good estimate of the experimental $s$ is necessary for
both criteria. Notably, a benchmarking of the error is always possible
experimentally [28] by preparing classical states of the mechanics, e.g.
thermal states at different temperatures. As a result, the reliability of our
nonclassicality criteria is not affected by the reduction in the
experimentally achievable $s$ due to noise and loss effects. Moreover, since
the demarginalization maps of Eqs. (25) and (26) are informationally
equivalent for all values of $s$, the above mentioned reduction does not
decrease the power of our nonclassicality criteria.

We now comment briefly on the effect of displacement-amplitude fluctuations
that cause random fluctuations in the interaction strength $\chi$ given by Eq.
(6). First, we note that both displacements would in practice be performed
using the same laser and can thus be implemented using inherently stable
designs, such as Sagnac or polarization-based interferometers. Second, our
scheme properly measures only one mechanical quadrature at a time, and it is
well-known that a single quadrature of a harmonic oscillator is a continuous
quantum non-demolition (QND) observable that can be measured precisely and
continuously [31]. In our scenario, this is reflected in the fact that for any
pulse length $\tau$ the mechanical bosonic annihilation operator $b^{\prime}$
commutes with itself (see Eq. (69) in Appendix B), meaning that there is no
back-action from light onto the mechanics in our measurement scheme [32].
Third, in light of Eq. (12) it turns out that fluctuations in the amplitude
and consequently $\chi$ only affect the width of the tomographic kernel,
namely it corresponds to small fluctuations in the parameter $s^{\star}$ given
by Eq. (15). Therefore, as long as fluctuations give rise to small relative
errors in the ordering parameter, they can be safely ignored. Finally, we
emphasize here that amplitude fluctuations are not amplified when squeezing is
exploited as it takes place before the displacement.

Besides the above analysis, it is important to consider the possibility of
false positive (or false negative) results due to experimental noise and
finite statistics. For this we first use the following Lemma from Ref. [33]:

######  Lemma 2.

The set of P-classical states is nowhere dense in the whole bosonic state
space [33].

Lemma 2 together with the fact that the set of P-classical states is closed
implies that this set cannot have an accumulation point within the bosonic
state space that is not P-classical. The latter means that, for every
P-nonclassical state there exists a sufficiently small epsilon-ball (in state
space) which contains no P-classical states. In other words, given a
P-nonclassical state, it will remain P-nonclassical under sufficiently small
perturbations. As a consequence, our (or any) method for nonclassicality
detection is robust to fluctuations and finite-statistics effects for false-
positive results, as long as these are small enough.

In practice, the fluctuations might of course be larger, such that they could
in principle lead to false conclusions. Similar problems are well-known from
quantum state tomography, where statistical fluctuations might, for example,
make a state look entangled, although it is separable. Such false positives
can be identified using standard techniques such as Monte-Carlo resampling
according to the known experimental error model. Such techniques are, for
instance, well established in the context of discrete-variable quantum
tomography and produce reliable confidence regions for the reconstructed
quantum state. In these methods the finite-statistic problem is taken care of
by choosing a sufficiently large number of samples, typically determined by
the properties of the state of interest through Hoeffding’s inequality.

On the other hand, given an initial nonclassical state, noisy measurements can
wash out the nonclassical features and lead to a false-negative result. While
this cannot be avoided in general, it should be noted that the arguments above
imply that sufficiently small fluctuations will preserve the nonclassical
features. However, the exact magnitude of noise and finite statistics that can
be tolerated depends on the nonclassicality of the state under test and the
conditions of the experiment.

###  2.5 Experimental Proposal

We now describe a practical protocol for performing tomography of a mechanical
resonator using our technique. First, recall that the phase $\varphi$ is set
by the pulse duration $\tau$ via the condition of Eq. (8) and substituting
into Eq. (5). For a given pulse duration and depending on the optomechanical
system parameters, the best pulse amplitude is then $|\alpha|=r$ such that the
energy-time uncertainty relation is satisfied. Although using these parameters
fixes the readout quadrature to $\varphi=\varphi_{0}$, Eq. (20) shows that
arbitrary quadrature angles $\varphi_{\rm d}$ can be probed by exploiting the
free evolution of the mechanics in the time delay $\tau_{\rm d}$ between state
preparation and readout, which can be precisely controlled in practice, see
Fig 1a. For each quadrature angle, one then measures the momentum tomogram of
the Wigner function of the reflected light, from which the mechanical tomogram
can be extracted by simple rescaling according to Eq. (20). Repeating this
process for each required tomogram, on obtains enough information to fully
reconstruct the Wigner function of the mechanical resonator.

We will now study the parameters of a few state-of-the-art optomechanical
experiments to test the feasibility of our approach. For this analysis we will
consider no squeezing in the initial optical state, i.e. $\varepsilon=0$, as
squeezing remains a challenging task in practice. Furthermore, since the
optical pulses required for all these systems are relatively long, it incurs
no experimental complications to assume that the pulse length can be chosen
close to the optimal $\tau\sim\pi/\omega_{mec}$, with a relative stability on
the order of $10^{-4}$, which, using Eqs. (5) and (8), implies $u\sim 2$ and
$\chi\sim\sqrt{k/2\sqrt{3}}$ for $k\in 2\mathbb{Z}^{+}$. In practice, these
pulses would be generated by chopping continuous-wave lasers and can thus not
only be timed very precisely, but also inherit the bandwidth and coherence
properties of the parent laser enabling displacement-counter-displacement
operations (as shown in Fig. 1a) to high accuracy. We will study the more
difficult task of reconstructing the mechanical Wigner function, rather than
other $s$-parameterized quasiprobability distributions. We thus compute the
requirements for $\chi\sim 3$ to satisfy the condition of Eq. (19), up to the
order of $10^{-3}$, which is below other typical noise sources in the
experiment. Note that, Eq. (8) also requires a relative stability of the laser
power on the order of $10^{-3}$.

Reference |  $\omega_{\text{m}}$ [Hz] | mass [kg] |  $\Gamma_{\text{m}}$ [Hz] |  $g_{0}$ [Hz] |  $\kappa_{\text{o}}$ [Hz] | $\kappa_{\rm o}/\omega_{\text{m}}$ |  opt $\tau$ [s] | pulse Energy [J]
---|---|---|---|---|---|---|---|---
Kleckner et al [34] | $9.7\times 10^{3}$ | $1.1\times 10^{-10}$ | $1.3\times 10^{-2}$ | $22$ | $4.7\times 10^{5}$ | 55 | $\sim 2.1\times 10^{-4}$ | $\sim 5.5\times 10^{-13}$
Murch et al [40] | $4.2\times 10^{4}$ | $10^{-22}$ | $10^{3}$ | $6\times 10^{5}$ | $6.6\times 10^{5}$ | 15.7 | $\sim 4.8\times 10^{-5}$ | $\sim 1.4\times 10^{-20}$
Norte et al [35] | $1.5\times 10^{5}$ | $10^{-12}$ | $1.4\times 10^{-3}$ | $10^{2}$ | $10^{6}$ | 6 | $\sim 1.3\times 10^{-5}$ | $\sim 6.4\times 10^{-12}$
Thompson et al [41] | $1.3\times 10^{5}$ | $4\times 10^{-11}$ | $0.12$ | $50$ | $5\times 10^{5}$ | 3.7 | $\sim 1.5\times 10^{-5}$ | $\sim 1.9\times 10^{-11}$
Anguiano et al [42] | $20\times 10^{9}$ | $7.7\times 10^{-12}$ | $2\times 10^{6}$ | $4.8\times 10^{7}$ | $3.4\times 10^{10}$ | 1.72 | $\sim 10^{-10}$ | $\sim 4.9\times 10^{-13}$
Arcizet et al [43] | $8.1\times 10^{5}$ | $1.9\times 10^{-7}$ | $81$ | $1.2$ | $10^{6}$ | 1.3 | $\sim 2.5\times 10^{-6}$ | $\sim 1.3\times 10^{-6}$
Cuthbertson et al [44] | $10^{3}$ | $1.85$ | $2.5\times 10^{-6}$ | $1.2\times 10^{-3}$ | 275 | 0.9 | $\sim 2\times 10^{-3}$ | $\sim 2.0\times 10^{-6}$
Gröblacher et al [36] | $9.5\times 10^{5}$ | $1.4\times 10^{-10}$ | $1.4\times 10^{2}$ | $3.9$ | $2\times 10^{5}$ | 0.22 | $\sim 2.1\times 10^{-6}$ | $\sim 1.7\times 10^{-7}$
Chan et al [45] | $3.9\times 10^{9}$ | $3.1\times 10^{-16}$ | $3.9\times 10^{3}$ | $9\times 10^{5}$ | $5\times 10^{8}$ | 0.13 | $\sim 5.1\times 10^{-10}$ | $\sim 5.3\times 10^{-11}$
Verhagen et al [46] | $7.8\times 10^{7}$ | $1.9\times 10^{-12}$ | $3.4\times 10^{3}$ | $3.4\times 10^{3}$ | $7.1\times 10^{6}$ | 0.09 | $\sim 2.6\times 10^{-8}$ | $\sim 1.5\times 10^{-9}$
Teufel et al [47] | $1.1\times 10^{7}$ | $4.8\times 10^{-14}$ | $32$ | $2\times 10^{2}$ | $2\times 10^{5}$ | 0.02 | $\sim 1.8\times 10^{-7}$ | $\sim 8.5\times 10^{-9}$
Table 1: Experimental parameters for different optomechanical systems. For
each system we provide the optimal pulse length $\tau$ for the probe pulse
used by our scheme, together with the required pulse energy for Wigner-
function tomography based on the parameter assumptions $\varepsilon=0$, $u=2$,
and $\chi\sim 3$.

The relevant parameters for a number of optomechanical systems, taken mainly
from Ref. [1] are summarized in Tab. 1. We note that although many of these
systems are indeed in the resolved sideband regime, they could nonetheless
employ our technique by sufficiently decreasing the quality factor of the used
optical cavities, as outlined in Tab. 1. Let us now consider explicitly kHz-
frequency resonators such as $Si_{3}N_{4}$ membranes [34, 35], which are
promising candidates for ground-state cooling [36] and state preparation via
measurement-based methods [7]. In this case, the required pulse duration is on
the order of tens of microseconds, which can easily be achieved, even with
gain-switched CW-lasers. In view of Eqs. (6) and (8), the requirements on the
optical power to reach $\chi\sim 3$ then depends on the optomechanical
coupling strength $g_{0}$ for the respective system and the readout
displacement amplitude. In practice, optical powers on the order of mW or
below are sufficient for almost all systems presented in Tab. 1. Evidently, it
must be ensured that these powers are tolerated by the corresponding
mechanical resonators, which are typically operated at lower optical powers.
Should it be necessary to decrease the optical power impinging on the
resonator to prevent damaging it, one could consider using a train of readout
pulses with a repetition rate equal to the mechanical frequency and a total
duration much shorter than the mechanical decay time. This approach greatly
decreases the required optical power without incurring significant added
noise.

Note, that the input power requirements of our method can be decreased
significantly by using a low finesse optical cavity, such that the optical
decay rate $\kappa_{\rm o}$ is much larger than the mechanical frequency
$\omega_{\rm m}$. As an example, the system of Ref. [34] in Tab. 1 has a
suitable cavity with a sideband resolution factor of $\kappa_{\rm
o}/\omega_{\rm mec}\sim 50$, and is thus in the ideal regime for our method.
Making use of the optical cavity with a finesse of $\sim 3\times 10^{4}$, the
required optical power is on the order of $10^{-13}$ W. Optical displacements
[37] and homodyne tomography of the optical output field can then be performed
using standard methods [38, 39].

Finally, we would like to point out that the use of squeezed readout light,
though not necessary, would be highly beneficial for our protocol through
reducing the requirements on amplitude and precision of the optical
displacement operations. Specifically, one could choose $\chi\sim 1$ which, in
turn, requires a squeezing parameter $\varepsilon\sim\ln(3)$ to reach the same
approximation as before. This is equivalent to a vacuum squeezing of $\sim
9.5$dB, which is achievable with current technology. In this case, both the
requirements of Eqs. (8) and (19) are satisfied easily without high precision
in the control of the laser power. Importantly, we note that squeezing does
not affect the required stability in the laser power, as it takes place before
the displacement. Furthermore, the squeezing parameter does not enter into the
interaction transformation of Eq. (9), implying that it need not even be
precisely characterized, as long as Eq. (19) is satisfied.

##  3 Discussion

We have introduced a method that allows for the tomographic reconstruction of
the phase space quasiprobability distributions of a mechanical resonator in a
regime where the commonly used rotating-wave approximation does not hold. In
contrast to conventional methods in which motional sidebands are resolved with
the help of a high quality optical cavity, our approach fits within the bad or
no cavity regime. As a result, our technique is applicable to a wide range of
experiments. Taking into account the full optomechanical Hamiltonian, we
demonstrated that, even in this regime, the mechanical quadrature
distributions can be imprinted onto an optical field through the radiation
pressure interaction. Carefully controlling for back-action, optical
displacements, and the optical Kerr effect, we showed how this information can
be extracted using established methods and current technology. Finally, we
discussed how our approach could be used to witness nonclassicality of the
mechanical state in a much more resource-efficient way without requiring full
quantum state reconstruction. We anticipate that our tomographic technique
will enable the reconstruction of the quantum states of motion of mechanical
resonators, and verify their potential nonclassicality, in a wide range of
experimental platforms, without the challenging requirements of alternative
approaches.

We thank M. Vanner, A.G. White, W. Bowen, and A. Szorkovszky for helpful
discussions. This work was supported in part by the ARC Centres of Excellence
for Engineered Quantum Systems (CE110001013, CE170100009) and Quantum
Computation and Communication Technology (CE110001027), and the UK Engineering
and Physical Sciences Research Council (grant number EP/N002962/1). F.S.
developed the concepts. F.S. and M.R. worked out the calculations and
experimental parameters. All authors wrote the manuscript.

##  4 Methods

###  4.1 Phase-space representation

Bosonic continuous-variable quantum systems are best described using a quantum
phase-space representation with respect to two orthogonal quadratures, such as
position and momentum. For such a representation one requires an
informationally-complete set of basis operators, which are commonly taken to
be the Weyl-Wigner operators
$\hat{T}(\alpha;s)=\int\frac{d^{2}\xi}{\pi}\exp\\{\alpha\xi^{*}-\alpha^{*}\xi+\frac{s}{2}|\xi|^{2}\\}\hat{D}(\xi)$,
for a complex-valued phase-space point $\alpha\in\mathbb{C}$ and distribution
parameter $s\in[-1,1]$ [8]. Here, $\hat{D}(\xi){=}\exp\\{\xi
a^{\dagger}-\xi^{*}a\\}$ is the displacement operator by the value
$\xi\in\mathbb{C}$, and $a$ ($a^{\dagger}$) is the bosonic annihilation
(creation) operator. Weyl-Wigner operators are normal, complete, and satisfy
the orthogonality (duality) relation,

$\displaystyle{\rm Tr}\hat{T}(\alpha;s){=}1,$ (27)
$\displaystyle\int\frac{d^{2}\alpha}{\pi}\hat{T}(\alpha;s)=\hat{I},$ (28)
$\displaystyle{\rm
Tr}\hat{T}(\alpha;s)\hat{T}(\beta;-s){=}\pi\delta^{(2)}(\alpha{-}\beta),$ (29)

respectively. As a result, we can expand any operator $\hat{\Lambda}$ in this
basis as

$\hat{\Lambda}=\int\frac{d^{2}\alpha}{\pi}\mathscr{W}_{\hat{\Lambda}}(\alpha;s)\hat{T}(\alpha;-s).$
(30)

The distribution $\mathscr{W}_{\hat{\Lambda}}(\alpha;s)={\rm
Tr}\hat{\Lambda}\hat{T}(\alpha;s)$ is known as the $s$-parameterized
quasiprobability distribution of the operator $\hat{\Lambda}$, which includes
as special cases [8] the Husimi-Kano Q-function [48, 49] for $s=-1$, the
Wigner [50] function for $s=0$, and the Glauber-Sudarshan P-function [51, 52]
for $s=1$. Crucially, these distributions, which are easily generalized to the
multipartite case [9], are not, in general, legitimate probability
distributions, as they can take negative values. If the operator
$\hat{\Lambda}$ is a quantum state, then such negative values are the primary
signature of nonclassical behaviour [9, 53], and can always be detected using
direct observation via tomographic methods [12, 20, 21, 22] or nonclassicality
criteria [54, 23, 24, 25, 26, 27, 55].

Using Eq. (30) and the duality relation of Eq. (29), the trace operation can
be expressed in terms of phase-space distributions as

${\rm
Tr}\hat{\Lambda}\hat{\Upsilon}=\int\frac{d^{2}\alpha}{\pi}\mathscr{W}_{\hat{\Lambda}}(\alpha;s)\mathscr{W}_{\hat{\Upsilon}}(\alpha;-s).$
(31)

Finally, it is sometimes more natural to work with position-momentum
coordinates, in which case $d^{2}\alpha=d{\rm Re}\alpha\,d{\rm
Im}\alpha=\frac{dqdp}{2}$, and
$\mathscr{W}_{\hat{\Lambda}}(\alpha;s)=2\pi\mathscr{W}_{\hat{\Lambda}}(q,p;s)$,
so that

$\begin{split}\int\frac{d^{2}\alpha}{\pi}\mathscr{W}_{\hat{\Lambda}}(\alpha;s)=\int
dqdp\mathscr{W}_{\hat{\Lambda}}(q,p;s)=1.\end{split}$ (32)

###  4.2 Phase-space tomography

While the quasiprobability distribution of an unknown quantum state cannot be
measured directly, it can be tomographically reconstructed using a set of
experimentally accessible marginal quadrature distributions or _tomograms_.
Recall that the $\phi$-quadrature of a bosonic field is given by
$\hat{x}(\phi)=\left(ae^{-i\phi}+a^{\dagger}e^{i\phi}\right)/\sqrt{2}$ with
eigenvalues and eigenvectors defined by the eigenvalue equation
$\hat{x}(\phi)|x(\phi)\rangle=x(\phi)|x(\phi)\rangle$. The set of all
projections onto quadrature eigenstates
$\\{\hat{\Pi}(x,\phi)=|x(\phi)\rangle\langle x(\phi)|\\}$ for $x\in\mathbb{R}$
and $\phi\in[0,2\pi)$, known as _tomographic operators_ , forms an
informationally complete operator basis [12, 56]. One can then define the set
of tomograms $\\{w(x,\phi)\\}$ for a quantum state $\hat{\varrho}$ as

$w(x,\phi)={\rm Tr}\hat{\varrho}\hat{\Pi}(x,\phi)=\langle
x(\phi)|\hat{\varrho}|x(\phi)\rangle,$ (33)

Using Eq. (31), this tomographic relation can be expressed in terms of Wigner
functions as,

$w(x,\phi)=\int\frac{d^{2}\alpha}{\pi}\mathscr{W}_{\hat{\varrho}}(\alpha;0)\mathscr{W}_{\hat{\Pi}(x,\phi)}(\alpha;0),$
(34)

where $\mathscr{W}_{\hat{\varrho}}(\alpha;0)$ is the Wigner function of the
quantum state, and

$\mathscr{W}_{\hat{\Pi}(x,\phi)}(\alpha;0)=\delta(x-\frac{e^{i\phi}\alpha^{*}+e^{-i\phi}\alpha}{\sqrt{2}})$
(35)

is the Wigner function of the tomographic operator $\hat{\Pi}(x,\phi)$.
Equation (34), also known as a _Radon transformation_ , can be interpreted as
sampling the Wigner function of a quantum state $\hat{\varrho}$ with strings
of zero width, see Fig. 4a. In turn, a sufficiently large set of tomograms can
be used to reconstruct the Wigner function via an _inverse Radon
transformation_ , see Fig. 3.

Figure 3: Traditional Wigner-function tomography. (a) The continuous-variable
quantum system, in an initially unknown quantum state, is subject to
quadrature measurements, for example, through homodyne detection. (b) Using a
sufficiently large number of quadrature measurements (depending on the
complexity of the unknown state), an inverse Radon transformation (IRT) can be
used to reconstruct the Wigner function of the unknown state.

In short,

$\begin{split}&\mathscr{W}_{\hat{\varrho}}(\alpha;0)\quad\xmapsto[{\rm
Tr}\hat{\varrho}\hat{\Pi}(x,\phi)]{\text{Radon
transf.}}\quad\\{w(x,\phi)\\},\\\ &\\{w(x,\phi)\\}\quad\xmapsto{\text{Inverse
Radon transf.}}\quad\mathscr{W}_{\hat{\varrho}}(\alpha;0).\end{split}$ (36)

Figure 4: Radon transform over phase space. (a) We can think of the standard
Radon transform of Eq. (34) as sampling Wigner function of a quantum system
along different quadrature axes via a string of zero width (a Dirac delta
function cross section). (b) In the case of Eq. (37), one samples an
$s$-parametrized quasiprobability distribution with a string of finite width,
that is, a Gaussian cross section.

The tomographic relation (33) can be generalized to arbitrary phase-space
quasiprobability distributions with $s\geqslant 0$ as [57]

$w(x,\phi)=\int\frac{d^{2}\alpha}{\pi}\mathscr{W}_{\hat{\varrho}}(\alpha;s)\mathscr{W}_{\hat{\Pi}(x,\phi)}(\alpha;-s),$
(37)

with the tomographic operator

$\mathscr{W}_{\hat{\Pi}(x,\phi)}(\alpha;-s)=\frac{1}{\sqrt{s\pi}}\exp\\{-\frac{(x-\frac{e^{i\phi}\alpha^{*}+e^{-i\phi}\alpha}{\sqrt{2}})^{2}}{s}\\}.$
(38)

In analogy with the Radon transform, Eq. (37) can be interpreted as sampling
an $s$-parameterized phase-space quasiprobability distribution with a string
with Gaussian cross section of finite width, see Fig. 4b. Notably, Eqs. (34)
and (37) connect the _same_ experimentally measured tomograms $w(x,\phi)$ to
different quasiprobability representations of the underlying quantum state.

## Appendices

##  Appendix A Disentangling the interaction

The full unitary interaction can be written as

$\begin{split}&\hat{D}(-\alpha)e^{-i\tau_{0}\omega_{\rm
o}a^{\dagger}a-i\tau_{0}\omega_{\rm m}b^{\dagger}b}e^{-i\tau[\omega_{\rm
o}a^{\dagger}a+\omega_{\rm
m}b^{\dagger}b+g_{0}a^{\dagger}a(b^{\dagger}+b)]}\hat{D}(\alpha)=\\\
&\qquad\qquad\qquad\qquad\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a-i\tau_{0}\omega_{\rm m}b^{\dagger}b}e^{-i\tau[\omega_{\rm
m}b^{\dagger}b+g_{0}a^{\dagger}a(b^{\dagger}+b)]}\hat{D}(\alpha)\end{split}$
(39)

in which $\hat{D}(\alpha)=\exp\\{\alpha a^{\dagger}-\alpha^{*}a\\}$ is the
optical displacement operator, $\hat{\tau}$ is light pulse duration, and
$\tau_{0}$ is the time between the optomechanical interaction and the optical
displacement. The r.h.s of Eq. (39) follows from the fact that $a^{\dagger}a$
commutes with all other terms of the Hamiltonian. Now, to disentangle the
mechanical free evolution during the optomechanical interaction, we use the
disentangling relations (2.34) and (2.85)–(2.87) of Ref. [58], namely

$\begin{split}\exp\\{\theta[\alpha_{1}b^{\dagger}b+\alpha_{2}b+\alpha_{3}b^{\dagger}]\\}&=\exp\\{f_{3}b^{\dagger}\\}\exp\\{f_{1}b^{\dagger}b\\}\exp\\{f_{2}b\\}\exp\\{f_{\rm
x}(\theta)\\}\\\
&=\exp\\{f_{1}b^{\dagger}b\\}\exp\\{f_{4}b^{\dagger}\\}\exp\\{f_{2}b\\}\exp\\{f_{\rm
x}(\theta)\\}\\\
&=\exp\\{f_{1}b^{\dagger}b\\}\exp\\{f_{2}b+f_{4}b^{\dagger}\\}\exp\\{f_{\rm
x}(\theta)-\frac{1}{2}f_{2}f_{4}\\}\\\
&=\exp\\{f_{1}b^{\dagger}b\\}\exp\\{f_{2}b+f_{4}b^{\dagger}\\}\exp\\{f_{5}\\},\end{split}$
(40)

where

$\begin{split}&f_{1}=\alpha_{1}\theta,\qquad
f_{2}=\frac{\alpha_{2}}{\alpha_{1}}(e^{\theta\alpha_{1}}-1),\qquad
f_{3}=\frac{\alpha_{3}}{\alpha_{1}}(e^{\theta\alpha_{1}}-1),\qquad f_{\rm
x}=\frac{\alpha_{2}\alpha_{3}}{\alpha_{1}^{2}}(e^{\theta\alpha_{1}}-\theta\alpha_{1}-1),\\\
&f_{4}=f_{3}e^{-f_{1}}=\frac{\alpha_{3}}{\alpha_{1}}(1-e^{-\theta\alpha_{1}}),\qquad
f_{5}=\frac{\alpha_{2}\alpha_{3}}{\alpha_{1}^{2}}(\sinh(\theta\alpha_{1})-\theta\alpha_{1}).\end{split}$
(41)

This gives

$\begin{split}&\hat{D}(-\alpha)e^{-i\tau_{0}\omega_{\rm
o}a^{\dagger}a}e^{-i\tau_{0}\omega_{\rm m}b^{\dagger}b}e^{-i\tau[\omega_{\rm
o}a^{\dagger}a+\omega_{\rm
m}b^{\dagger}b+g_{0}a^{\dagger}a(b^{\dagger}+b)]}\hat{D}(\alpha)=\\\
&\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a}e^{-i(\tau+\tau_{0})\omega_{\rm
m}b^{\dagger}b}e^{\frac{g_{0}}{\omega_{\rm m}}a^{\dagger}a[(1-e^{i\omega_{\rm
m}\tau})b^{\dagger}+(e^{-i\omega_{\rm
m}\tau}-1)b]}e^{\frac{g_{0}^{2}}{\omega_{\rm
m}^{2}}(a^{\dagger}a)^{2}[i\omega_{\rm m}\tau-i\sin(\omega_{\rm
m}\tau)]}\hat{D}(\alpha)=\\\ &e^{-i(\tau+\tau_{0})\omega_{\rm
m}b^{\dagger}b}\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a}e^{\frac{g_{0}^{2}}{\omega_{\rm
mec}^{2}}(a^{\dagger}a)^{2}[i\omega_{\rm m}\tau-i\sin(\omega_{\rm
m}\tau)]}e^{\frac{g_{0}}{\omega_{\rm m}}a^{\dagger}a(\mu
b^{\dagger}-\mu^{*}b)}\hat{D}(\alpha)=\\\
&\underbrace{e^{-i(\tau+\tau_{0})\omega_{\rm
m}b^{\dagger}b}}\underbrace{\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a}\hat{D}(\alpha)}\underbrace{\hat{D}(-\alpha)e^{\frac{g_{0}^{2}}{\omega_{\rm
m}^{2}}(a^{\dagger}a)^{2}[i\omega_{\rm m}\tau-i\sin(\omega_{\rm
m}\tau)]}\hat{D}(\alpha)}\underbrace{\hat{D}(-\alpha)e^{\frac{g_{0}}{\omega_{\rm
m}}a^{\dagger}a(\mu b^{\dagger}-\mu^{*}b)}\hat{D}(\alpha)},\\\
&\mu=1-e^{i\omega_{\rm m}\tau},\end{split}$ (42)

in which $\hat{D}(\alpha)=\exp\\{\alpha a^{\dagger}-\alpha^{*}a\\}$ is the
optical displacement operator. In Eq. (42), whenever needed we have used the
fact that the optical and mechanical mode operators commute, and that
$a^{\dagger}a$ commutes with all other terms of the Hamiltonian. In the
resulting expression, we note that the first bracket is the free rotation of
the mechanics in phase space, which can be eliminated by moving to a rotating
frame with frequency $\omega_{\rm m}$.

It is important to note that the optical state that we will be reconstructing
contains all the evolutions described by the second, third, and fourth unitary
evolutions. Using the displacement relations

$\hat{D}(-\alpha)a\hat{D}(\alpha)=a+\alpha,\quad\hat{D}(-\alpha)a^{\dagger}\hat{D}(\alpha)=a^{\dagger}+\alpha^{*},$
(43)

the second bracket is understood to be the free rotation of the optical field
around a phase space point located at $-\alpha$. The inverse of this operation
can be easily applied to the reconstructed optical state and can thus be
neglected for now.

The third bracket represents a bare nonlinear Kerr effect on the optical mode.
The effect of this evolution can also be removed from the reconstructed
optical state using numerical techniques. However, in the linearized regime we
can take it into account as follows. Using Eq. (43), setting
$\alpha=re^{i\theta}$, and noting the fact that $|\alpha|=r\gg 1$, we can
approximate

$\begin{split}\hat{D}(-\alpha)e^{\frac{g_{0}^{2}}{\omega_{\rm
m}^{2}}(a^{\dagger}a)^{2}[i\omega_{\rm m}\tau-i\sin(\omega_{\rm
m}\tau)]}\hat{D}(\alpha)&\simeq e^{\frac{g_{0}^{2}r^{2}}{\omega_{\rm
m}^{2}}[4a^{\dagger}a+2r(e^{i\theta}a^{\dagger}+e^{-i\theta}a)+(e^{2i\theta}a^{{\dagger}2}+e^{-2i\theta}a^{2})+r^{2}+1][i\omega_{\rm
m}\tau-i\sin(\omega_{\rm m}\tau)]}\\\
&=e^{\chi^{2}[4a^{\dagger}a+2r(e^{i\theta}a^{\dagger}+e^{-i\theta}a)+(e^{2i\theta}a^{{\dagger}2}+e^{-2i\theta}a^{2})][i\omega_{\rm
m}\tau-i\sin(\omega_{\rm m}\tau)]},\\\ \chi=\frac{g_{0}r}{\omega_{\rm
m}},\end{split}$ (44)

where we have kept only the terms of order $r^{2}$ or higher, and have ignored
the overall phases corresponding to $r^{2}+1$. We can also rewrite the last
bracket of Eq. (42) as

$e^{\frac{g_{0}}{\omega_{\rm m}}(a^{\dagger}+\alpha^{*})(a+\alpha)(\mu
b^{\dagger}-\mu^{*}b)},$ (45)

and linearise the exponent as

$\begin{split}&\frac{g_{0}}{\omega_{\rm
m}}(a^{\dagger}+\alpha^{*})(a+\alpha)(\mu b^{\dagger}-\mu^{*}b)\\\
&\qquad\qquad=\frac{g_{0}}{\omega_{\rm m}}(\alpha^{*}\mu
ab^{\dagger}-\alpha\mu^{*}a^{\dagger}b)+\frac{g_{0}}{\omega_{\rm m}}(\alpha\mu
a^{\dagger}b^{\dagger}-\alpha^{*}\mu^{*}ab)+\frac{g_{0}}{\omega_{\rm
m}}a^{\dagger}a(\mu b^{\dagger}-\mu^{*}b)+\frac{g_{0}}{\omega_{\rm
m}}r^{2}(\mu b^{\dagger}-\mu^{*}b)\\\
&\qquad\qquad\simeq\frac{g_{0}ru}{\omega_{\rm
m}}(e^{i\theta}a^{\dagger}+e^{-i\theta}a+r)(e^{i\varphi}b^{\dagger}-e^{-i\varphi}b),\end{split}$
(46)

where $\mu=1-e^{i\omega_{\rm m}\tau}=ue^{i\varphi}$ so that

$\begin{split}&u=|\mu|=\sqrt{2(1-\cos{\omega_{\rm m}\tau})},\\\
&\varphi=\tan^{-1}{\left(\frac{\sin{\omega_{\rm m}\tau}}{1-\cos{\omega_{\rm
m}\tau}}\right)}=\tan^{-1}{\left(\cot{(\omega_{\rm
m}\tau/2)}\right)}=(m+1/2)\pi-\omega_{\rm m}\tau/2,\qquad
m\in\mathbb{Z}^{+}.\end{split}$ (47)

Curiously, the radiation pressure interaction term $\frac{g_{0}}{\omega_{\rm
m}}a^{\dagger}a(\mu b^{\dagger}+\mu^{*}b)$, has become negligible in Eq. (46),
since $\frac{g_{0}}{\omega_{\rm m}}\ll\frac{g_{0}ru}{\omega_{\rm m}}$ and can
thus be omitted. Note that this also requires that $u$ be chosen large enough,
which implies that the optimal pulse length $\tau$ is on the order of half the
mechanical period. We can now define the new optical mode annihilation
operator $\tilde{a}_{\theta}:=ae^{-i\theta}+\frac{r}{2}$ and the new
mechanical annihilation operator $\tilde{b}_{\varphi}=be^{-i\varphi}$ to get

$e^{\frac{g_{0}}{\omega_{\rm m}}(a^{\dagger}+\alpha^{*})(a+\alpha)(\mu
b^{\dagger}-\mu^{*}b)}\simeq e^{\chi
u(\tilde{a}^{\dagger}_{\theta}+\tilde{a}_{\theta})(\tilde{b}^{\dagger}_{\varphi}-\tilde{b}_{\varphi})}.$
(48)

Note that, since $0\leqslant u\leqslant 2$ and we work in a parameter regime
where $\chi u=\frac{g_{0}ru}{\omega_{\rm m}}\sim 1-10$, it follows that
$\chi^{2}\geqslant\chi$ and hence the optical Kerr effect in Eq. (44) is not
negligible compared to the interaction term of Eq. (48).

Using the linearised evolutions of Eqs. (44) and (48), we can now rewrite Eq.
(42) as

$\begin{split}&\hat{\mathcal{U}}_{\rm o}\leavevmode\nobreak\
\hat{\mathcal{U}}_{\rm k}\leavevmode\nobreak\ \hat{\mathcal{U}}_{\rm om},\\\
&\hat{\mathcal{U}}_{\rm om}=e^{\chi
u(\tilde{a}^{\dagger}_{\theta}+\tilde{a}_{\theta})(\tilde{b}^{\dagger}_{\varphi}-\tilde{b}_{\varphi})},\\\
&\hat{\mathcal{U}}_{\rm
k}=e^{\chi^{2}[4a^{\dagger}a+2r(e^{i\theta}a^{\dagger}+e^{-i\theta}a)+(e^{2i\theta}a^{{\dagger}2}+e^{-2i\theta}a^{2})][i\omega_{\rm
m}\tau-i\sin(\omega_{\rm m}\tau)]},\\\ &\hat{\mathcal{U}}_{\rm
o}=\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a}\hat{D}(\alpha).\end{split}$ (49)

##  Appendix B Transformation of the bosonic operators

###  B.1 Optomechanical interaction

As outlined in the main text, we can now use Eq. (49) to compute the
transformation of the vector of bosonic operators
$\hat{\vec{A}}:=(a^{\dagger},b^{\dagger},a,b)$ under the optomechanical
interaction given by

$\hat{\vec{A}}^{\prime}=\hat{\mathcal{U}}_{\rm
om}\hat{\vec{A}}\hat{\mathcal{U}}^{\dagger}_{\rm om}.$ (50)

Exploiting the properties of symplectic transformations [17], we know that if
$\hat{\mathcal{U}}=\exp(\frac{1}{2}\hat{\vec{A}}\ln
M\Sigma\hat{\vec{A}}^{\mathsf{T}})$, with $\Sigma=\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}$, then
$\hat{\mathcal{U}}\hat{\vec{A}}\hat{\mathcal{U}}^{\dagger}=\hat{\vec{A}}M$. In
terms of the new the new field operators $\tilde{a}_{\theta}$ and
$\tilde{b}_{\varphi}$,

$\begin{split}&\hat{\mathcal{U}}_{\rm om}=\exp\\{\frac{\chi
u}{2}\begin{pmatrix}\tilde{a}^{\dagger}_{\theta}&\tilde{b}^{\dagger}_{\varphi}&\tilde{a}_{\theta}&\tilde{b}_{\varphi}\end{pmatrix}\begin{pmatrix}0&-1&0&-1\\\
1&0&-1&0\\\ 0&-1&0&-1\\\
-1&0&1&0\end{pmatrix}\begin{pmatrix}\tilde{a}_{\theta}\\\
\tilde{b}_{\varphi}\\\ -\tilde{a}^{\dagger}_{\theta}\\\
-\tilde{b}^{\dagger}_{\varphi}\end{pmatrix}\\},\end{split}$ (51)

which gives

$\ln M=\chi u\begin{pmatrix}0&-1&0&-1\\\ 1&0&-1&0\\\ 0&-1&0&-1\\\
-1&0&1&0\end{pmatrix},$ (52)

and thus,

$\begin{split}&M=\begin{pmatrix}1&-\chi u&0&-\chi u\\\ \chi u&1&-\chi u&0\\\
0&-\chi u&1&-\chi u\\\ -\chi u&0&\chi u&1\end{pmatrix},\leavevmode\nobreak\
M^{-1}=\begin{pmatrix}1&\chi u&0&\chi u\\\ -\chi u&1&\chi u&0\\\ 0&\chi
u&1&\chi u\\\ \chi u&0&-\chi u&1\end{pmatrix},\end{split}$ (53)

with $\det M=1$. Using the matrix of Eq. (53), we find

$\begin{split}&\hat{\mathcal{U}}_{\rm
om}\tilde{a}^{\dagger}_{\theta}\hat{\mathcal{U}}_{\rm
om}^{\dagger}=\tilde{a}^{\dagger}_{\theta}+\chi
u(\tilde{b}^{\dagger}_{\varphi}-\tilde{b}_{\varphi})=a^{\dagger}e^{i\theta}+\frac{r}{2}-i\sqrt{2}\chi
u\hat{P}(\varphi),\\\ &\hat{\mathcal{U}}_{\rm
om}\tilde{b}^{\dagger}_{\varphi}\hat{\mathcal{U}}_{\rm
om}^{\dagger}=\tilde{b}^{\dagger}_{\varphi}-\chi
u(\tilde{a}^{\dagger}_{\theta}+\tilde{a}_{\theta})=b^{\dagger}e^{i\varphi}-\sqrt{2}\chi
u[\hat{x}(\theta)+r]\\\ &\hat{\mathcal{U}}_{\rm
om}\tilde{a}_{\theta}\hat{\mathcal{U}}_{\rm
om}^{\dagger}=\tilde{a}_{\theta}-\chi
u(\tilde{b}^{\dagger}_{\varphi}-\tilde{b}_{\varphi})=ae^{-i\theta}+\frac{r}{2}+i\sqrt{2}\chi
u\hat{P}(\varphi),\\\ &\hat{\mathcal{U}}_{\rm
om}\tilde{b}_{\varphi}\hat{\mathcal{U}}_{\rm
om}^{\dagger}=\tilde{b}_{\varphi}-\chi
u(\tilde{a}^{\dagger}_{\theta}+\tilde{a}_{\theta})=be^{-i\varphi}-\sqrt{2}\chi
u[\hat{x}(\theta)+r],\end{split}$ (54)

where
$\hat{P}(\varphi)=\frac{1}{i\sqrt{2}}(be^{-i\varphi}-b^{\dagger}e^{i\varphi})=\frac{1}{\sqrt{2}}(be^{-i(\varphi+\frac{\pi}{2})}+b^{\dagger}e^{i(\varphi+\frac{\pi}{2})})=\hat{X}(\varphi+\frac{\pi}{2})$
is the $(\varphi+\frac{\pi}{2})$-quadrature of the mechanics, and
$\hat{x}(\theta)=\frac{1}{\sqrt{2}}(a^{\dagger}e^{i\theta}+ae^{-i\theta})$ is
the $\theta$-quadrature of the optical field. Making use of Eq. (50),

$\begin{split}\hat{\mathcal{U}}_{\rm
om}\tilde{a}_{\theta}\hat{\mathcal{U}}^{\dagger}_{\rm
om}&=\hat{\mathcal{U}}_{\rm om}a\hat{\mathcal{U}}^{\dagger}_{\rm
om}e^{-i\theta}+\frac{r}{2}=a^{\prime}e^{-i\theta}+\frac{r}{2},\\\
\hat{\mathcal{U}}_{\rm om}\tilde{b}_{\varphi}\hat{\mathcal{U}}^{\dagger}_{\rm
om}&=\hat{\mathcal{U}}_{\rm om}b\hat{\mathcal{U}}^{\dagger}_{\rm
om}e^{-i\varphi}=b^{\prime}e^{-i\varphi},\end{split}$ (55)

and thus, by virtue of Eq. (54), we find the following relations similar to
that of Eq. (9) in the main text,

$\begin{split}&a^{\prime{\dagger}}=a^{\dagger}+i\sqrt{2}\chi
ue^{-i\theta}\hat{X}(\varphi+\frac{\pi}{2}),\\\
&b^{\prime{\dagger}}=b^{\dagger}-\sqrt{2}\chi
ue^{-i\varphi}[\hat{x}(\theta)+r],\\\ &a^{\prime}=a-i\sqrt{2}\chi
ue^{i\theta}\hat{X}(\varphi+\frac{\pi}{2}),\\\ &b^{\prime}=b-\sqrt{2}\chi
ue^{i\varphi}[\hat{x}(\theta)+r].\end{split}$ (56)

As a consequence of the properties of the symplectic group $M$ represents a
canonical (not necessarily unitary) transformation. This can be easily
verified by computing the bosonic commutation relations
$[a^{\prime{\dagger}},a^{\prime}]=[b^{\prime{\dagger}},b^{\prime}]=1$ and
$[a^{\prime},b^{\prime}]=[a^{\prime{\dagger}},b^{\prime{\dagger}}]=0$. Note
that the last term of the linearised Hamiltonian, Eq. (46), only causes a
constant displacement of the mechanics and could thus be absorbed into an
initial displacement $b\mapsto b-\sqrt{2}\chi ure^{i\varphi}$. Finally, the
transformation corresponding to Eq. (56) is given by

$\begin{split}S&=\begin{pmatrix}1&-\chi ue^{i\theta}e^{-i\varphi}&0&-\chi
ue^{i\theta}e^{i\varphi}\\\ i\chi
ue^{-i\theta}e^{i(\varphi+\frac{\pi}{2})}&1&-i\chi
ue^{i\theta}e^{i(\varphi+\frac{\pi}{2})}&0\\\ 0&-\chi
ue^{-i\theta}e^{-i\varphi}&1&-\chi ue^{-i\theta}e^{i\varphi}\\\ i\chi
ue^{-i\theta}e^{-i(\varphi+\frac{\pi}{2})}&0&-i\chi
ue^{i\theta}e^{-i(\varphi+\frac{\pi}{2})}&1\end{pmatrix}\\\
&=\begin{pmatrix}1&-\chi ue^{i\theta}e^{-i\varphi}&0&-\chi
ue^{i\theta}e^{i\varphi}\\\ -\chi ue^{-i\theta}e^{i\varphi}&1&\chi
ue^{i\theta}e^{i\varphi}&0\\\ 0&-\chi ue^{-i\theta}e^{-i\varphi}&1&-\chi
ue^{-i\theta}e^{i\varphi}\\\ \chi ue^{-i\theta}e^{-i\varphi}&0&-\chi
ue^{i\theta}e^{-i\varphi}&1\end{pmatrix},\\\ S^{-1}&=\begin{pmatrix}1&\chi
ue^{i\theta}e^{-i\varphi}&0&\chi ue^{i\theta}e^{i\varphi}\\\ \chi
ue^{-i\theta}e^{i\varphi}&1&-\chi ue^{i\theta}e^{i\varphi}&0\\\ 0&\chi
ue^{-i\theta}e^{-i\varphi}&1&\chi ue^{-i\theta}e^{i\varphi}\\\ -\chi
ue^{-i\theta}e^{-i\varphi}&0&\chi ue^{i\theta}e^{-i\varphi}&1\end{pmatrix},\\\
\det{S}&=1.\end{split}$ (57)

###  B.2 Kerr interaction

Along the same lines we can now evaluate the effect of the Kerr interaction on
the optical field. Again, given the vector of bosonic operators
$\hat{\vec{A}}:=(a^{\dagger},b^{\dagger},a,b)$, we evaluate the output mode
operators

$\hat{\vec{A}}^{\prime}=\hat{\mathcal{U}}_{\rm
k}\hat{\vec{A}}\hat{\mathcal{U}}^{\dagger}_{\rm k}.$ (58)

To find the result, we first define the new mode operators [18]

$\bar{a}=ae^{-i\theta}+\frac{r}{3}.$ (59)

Substituting Eq. (59) into the exponent of $\hat{\mathcal{U}}_{\rm k}$ in Eq.
(49), we obtain

$4a^{\dagger}a+2r(e^{i\theta}a^{\dagger}+e^{-i\theta}a)+(e^{2i\theta}a^{{\dagger}2}+e^{-2i\theta}a^{2})=4\bar{a}^{\dagger}\bar{a}+\bar{a}^{{\dagger}2}+\bar{a}^{2}-\frac{6}{9}r^{2}.$

The $r^{2}$ term is just an overall phase factor which can be changed
arbitrarily. Therefore, in terms of the new the new field operators $\bar{a}$,
we have

$\begin{split}&\hat{\mathcal{U}}_{\rm k}=\exp\\{\chi^{2}[i\omega_{\rm
m}\tau-i\sin(\omega_{\rm
m}\tau)]\begin{pmatrix}\bar{a}^{\dagger}&\bar{a}\end{pmatrix}\begin{pmatrix}2&-1\\\
1&-2\end{pmatrix}\begin{pmatrix}\bar{a}\\\
-\bar{a}^{\dagger}\end{pmatrix}\\},\end{split}$

Using the same relation as in Sec. B.1,

$\ln M=2iv\begin{pmatrix}2&-1\\\ 1&-2\end{pmatrix},$

where $v=\chi^{2}[\omega_{\rm m}\tau-\sin(\omega_{\rm m}\tau)]$, and thus,

$\begin{split}&M=\begin{pmatrix}\cos{\bar{v}}+\frac{2i}{\sqrt{3}}\sin{\bar{v}}&-\frac{i}{\sqrt{3}}\sin{\bar{v}}\\\
\frac{i}{\sqrt{3}}\sin{\bar{v}}&\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}}\end{pmatrix},\quad
M^{-1}=\begin{pmatrix}\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}}&\frac{i}{\sqrt{3}}\sin{\bar{v}}\\\
-\frac{i}{\sqrt{3}}\sin{\bar{v}}&\cos{\bar{\bar{v}}}+\frac{2i}{\sqrt{3}}\sin{\bar{v}}\end{pmatrix},\end{split}$
(60)

in which $\bar{v}=2\sqrt{3}v$, and we have . Using the matrix of Eq. (60), we
find

$\begin{split}\hat{\mathcal{U}}_{\rm k}\bar{a}\hat{\mathcal{U}}_{\rm
k}^{\dagger}&=-\frac{i}{\sqrt{3}}\bar{a}^{\dagger}\sin{\bar{v}}+\bar{a}[\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}}]\\\
&=-\frac{i}{\sqrt{3}}a^{\dagger}e^{i\theta}\sin{\bar{v}}+ae^{-i\theta}(\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}})+(\frac{1}{3}\cos{\bar{v}}-\frac{i}{\sqrt{3}}\sin{\bar{v}})r.\end{split}$
(61)

Next, we have,

$\begin{split}\hat{\mathcal{U}}_{\rm k}\bar{a}\hat{\mathcal{U}}^{\dagger}_{\rm
k}&=\hat{\mathcal{U}}_{\rm k}a\hat{\mathcal{U}}^{\dagger}_{\rm
k}e^{-i\theta}+\frac{r}{3}=a^{\prime}e^{-i\theta}+\frac{r}{3}.\end{split}$
(62)

Equating this with Eq. (61), we find

$a^{\prime}=-\frac{i}{\sqrt{3}}a^{\dagger}e^{2i\theta}\sin{\bar{v}}+a(\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}})+e^{i\theta}(\frac{1}{3}\cos{\bar{v}}-\frac{i}{\sqrt{3}}\sin{\bar{v}}-\frac{1}{3})r.$
(63)

The linear transformation corresponding to Eq. (63) is thus given by

$\hat{\vec{A}}^{\prime}=\hat{\vec{A}}\begin{pmatrix}\cos{\bar{v}}+\frac{2i}{\sqrt{3}}\sin{\bar{v}}&-\frac{i}{\sqrt{3}}e^{2i\theta}\sin{\bar{v}}\\\
\frac{i}{\sqrt{3}}e^{-2i\theta}\sin{\bar{v}}&\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}}\end{pmatrix}+\vec{D},$
(64)

where the displacement vector $\vec{D}$ is defined as

$\vec{D}=(D^{*}\quad
D)=\frac{r}{3}\left(e^{i\theta}(\cos{\bar{v}}+i\sqrt{3}\sin{\bar{v}}-1)\leavevmode\nobreak\
,\quad e^{-i\theta}(\cos{\bar{v}}-i\sqrt{3}\sin{\bar{v}}-1)\right).$ (65)

###  B.3 Optical transformation

We now consider the effect of the last unitary evolution in Eq. (49), namely,
$\hat{\mathcal{U}}_{\rm o}=\hat{D}(-\alpha)e^{-i(\tau_{0}+\tau)\omega_{\rm
o}a^{\dagger}a}\hat{D}(\alpha)$. This time, it is very easy to find the
corresponding transformation in the Wigner representation using Lemma 1 (see
the next section C for the proof this lemma). The action of the displacement
operator on the mode operators is given by Eq. (43), while the action of the
phase-rotation operator is well-known to be

$e^{-i\Theta a^{\dagger}a}ae^{i\Theta a^{\dagger}a}=ae^{i\Theta},\qquad
e^{-i\Theta a^{\dagger}a}a^{\dagger}e^{i\Theta
a^{\dagger}a}=a^{\dagger}e^{-i\Theta}.$ (66)

It is then straightforward to evaluate the following:

$a^{\prime}=\hat{\mathcal{U}}_{\rm o}a\hat{\mathcal{U}}_{\rm
o}^{\dagger}=ae^{i(\tau_{0}+\tau)\omega_{\rm
o}}+(e^{i(\tau_{0}+\tau)\omega_{\rm o}}-1)\alpha.$ (67)

This represents a rotation of the optical field in phase space by an angle of
$(\tau_{0}+\tau)\omega_{\rm o}$ around the phase-space point $\alpha$.

###  B.4 The total evolution

To obtain the effect of the overall evolution operator as described in Eq.
(49), we simply combine the results of Eqs (50), (63), and (67). First,
combining the Kerr effect with the optomechanical interaction gives

$\begin{split}&a^{\prime}=-\frac{i}{\sqrt{3}}a^{\dagger}e^{2i\theta}\sin{\bar{v}}+a(\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}})+(\frac{1}{3}\cos{\bar{v}}-\frac{i}{\sqrt{3}}\sin{\bar{v}}-\frac{1}{3})e^{i\theta}r-i\sqrt{2}\chi
ue^{i\theta}\hat{X}(\varphi+\frac{\pi}{2}),\\\ &b^{\prime}=b-\sqrt{2}\chi
ue^{i\varphi}[\hat{x}(\theta)+r].\end{split}$ (68)

Next, after the optical displacement, the final transformation is obtained as

$\begin{split}&a^{\prime}=-\frac{i}{\sqrt{3}}a^{\dagger}e^{-i[(\tau_{0}+\tau)\omega_{\rm
o}-2\theta]}\sin{\bar{v}}+ae^{i(\tau_{0}+\tau)\omega_{\rm
o}}(\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}})\\\
&\qquad+re^{i\theta}\left[\left(e^{i(\tau_{0}+\tau)\omega_{\rm
o}}-\frac{2}{3}\right)\cos{\bar{v}}-\frac{i}{\sqrt{3}}\left(e^{-i(\tau_{0}+\tau)\omega_{\rm
o}}+2e^{i(\tau_{0}+\tau)\omega_{\rm
o}}-2\right)\sin{\bar{v}}-\frac{1}{3}\right]\\\ &\qquad-i\sqrt{2}\chi
ue^{i\theta}\hat{X}(\varphi+\frac{\pi}{2}),\\\ &b^{\prime}=b-\sqrt{2}\chi
ue^{i\varphi}[\hat{x}(\theta)+r].\end{split}$ (69)

In Eq. (69), we have considered the most general situation to obtain the
transformations. However, these expressions can be simplified by some general
considerations. One instance, presented in the main text, is as follows. The
initial phase of the optical mode, $\theta$, is being set by the internal
clock of the probe laser which will be further used for the homodyne
measurement of the light. Therefore, we can set $\theta=0$, that is
$\alpha=r\in\mathbb{R}$. Consequently,

$\begin{split}&a^{\prime}=-\frac{i}{\sqrt{3}}a^{\dagger}e^{-i(\tau_{0}+\tau)\omega_{\rm
o}}\sin{\bar{v}}+ae^{i(\tau_{0}+\tau)\omega_{\rm
o}}(\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}})\\\
&\qquad+r\left[\left(e^{i(\tau_{0}+\tau)\omega_{\rm
o}}-\frac{2}{3}\right)\cos{\bar{v}}-\frac{i}{\sqrt{3}}\left(e^{-i(\tau_{0}+\tau)\omega_{\rm
o}}+2e^{i(\tau_{0}+\tau)\omega_{\rm
o}}-2\right)\sin{\bar{v}}-\frac{1}{3}\right]\\\ &\qquad-i\sqrt{2}\chi
u\hat{X}(\varphi+\frac{\pi}{2}),\\\ &b^{\prime}=b-\sqrt{2}\chi
ue^{i\varphi}[\hat{x}+r].\end{split}$ (70)

Another simplification can be obtained by noting that one can arrange the
delay for the displacement of the output optical field ($\tau_{0}$) with high
accuracy in such a way that $(\tau_{0}+\tau)\omega_{\rm o}=m2\pi$ for
$m\in\mathbb{N}$. This gives

$\begin{split}&a^{\prime}=-\frac{i}{\sqrt{3}}a^{\dagger}\sin{\bar{v}}+a(\cos{\bar{v}}-\frac{2i}{\sqrt{3}}\sin{\bar{v}})\\\
&\qquad+r\left(\frac{1}{3}\cos{\bar{v}}-\frac{i}{\sqrt{3}}\sin{\bar{v}}-\frac{1}{3}\right)\\\
&\qquad-i\sqrt{2}\chi u\hat{X}(\varphi+\frac{\pi}{2}),\\\
&b^{\prime}=b-\sqrt{2}\chi ue^{i\varphi}[\hat{x}+r].\end{split}$ (71)

Finally, if we choose the probe pulse duration $\tau$ in such a way that
$\sin(\bar{v})=0$, we will obtain the simplest transformation equations. For
this, we recall that $\bar{v}=2\sqrt{3}v=2\sqrt{3}\chi^{2}[\omega_{\rm
m}\tau-\sin(\omega_{\rm m}\tau)]$. As $\tau\neq 0\Rightarrow\bar{v}\neq 0$, we
must have

$\bar{v}=k\pi,\quad\text{for}\quad
k\in\mathbb{N}\qquad\Rightarrow\qquad\omega_{\rm m}\tau-\sin(\omega_{\rm
m}\tau)=\frac{k\pi}{2\sqrt{3}\chi^{2}}.$ (72)

In general, it is always possible to find a pair of parameters $(\tau,r)$,
such that $0\ll\omega_{\rm m}\tau\ll 2\pi$, and Eq. (72) is satisfied.
However, as indicated earlier, the optimal pulse length to achieve $u=2$ is
$\omega_{\rm m}\tau=\pi$. In this case, the condition in Eq. (72) can be
satisfied by adjust the amplitude $r$ of the optical coherent state, and thus
the parameter $\chi$. Chosing $\tau$ such that $\bar{v}=k\pi$, we obtain

$\begin{split}&\left\\{\begin{matrix}[l]a^{\prime}=-a-\frac{2}{3}r-i\sqrt{2}\chi
u\hat{X}(\varphi+\frac{\pi}{2})&\quad k\text{\leavevmode\nobreak\ odd}\\\
a^{\prime}=a-i\sqrt{2}\chi u\hat{X}(\varphi+\frac{\pi}{2})&\quad
k\text{\leavevmode\nobreak\ even}\end{matrix}\right.,\\\ &\quad
b^{\prime}=b-\sqrt{2}\chi ue^{i\varphi}[\hat{x}+r],\end{split}$ (73)

where $u$ and $\varphi$ are given by Eq. (47). We see that, for even $k$, the
displacement term completely goes away, and we arrive at Eq. (9), that is,

$\begin{split}&a^{\prime}=a-i\sqrt{2}\chi u\hat{X}(\varphi+\frac{\pi}{2}),\\\
&b^{\prime}=b-\sqrt{2}\chi ue^{i\varphi}[\hat{x}+r].\end{split}$ (74)

To summarize our parameter choices, consider the case where we choose the
optimal pulse length $\tau$ such that $\omega_{\rm m}\tau=\pi$. This implies
$u=2$ and by Eq. (72), for example for $k=32$,
$\chi=\sqrt{\frac{16}{\sqrt{3}}}\approx 3.04$, which is within the required
range for Wigner function tomography using our method, without exploiting
squeezing. Alternatively, if one does not have precise control over the
amplitude $r$ of the input optical coherent state, then for any value of
$\chi\sim 1-10$, one can find a pulse length that satisfies the condition in
Eq. (72) with $u\neq 0$, for the purpose of Wigner function tomography.
Importantly, for any pair of parameters $(\tau,r)$ satisfying Eq. (72) such
that $0\ll\omega_{\rm m}\tau\ll 2\pi$, a $s$-parametrized quasiprobability
distribution of the mechanics is possible, as described in the main text.

From Eq. (74), after disregarding the mechanical constant displacement, we
find the following symplectic matrices,

$\begin{split}&S=\begin{pmatrix}1&-\chi ue^{-i\varphi}&0&-\chi
ue^{i\varphi}\\\ -\chi ue^{i\varphi}&1&\chi ue^{i\varphi}&0\\\ 0&-\chi
ue^{-i\varphi}&1&-\chi ue^{i\varphi}\\\ \chi ue^{-i\varphi}&0&-\chi
ue^{-i\varphi}&1\end{pmatrix},\qquad S^{-1}=\begin{pmatrix}1&\chi
ue^{-i\varphi}&0&\chi ue^{i\varphi}\\\ \chi ue^{i\varphi}&1&-\chi
ue^{i\varphi}&0\\\ 0&\chi ue^{-i\varphi}&1&\chi ue^{i\varphi}\\\ -\chi
ue^{-i\varphi}&0&\chi ue^{-i\varphi}&1\end{pmatrix},\\\
&\det{S}=1.\end{split}$ (75)

##  Appendix C Proof of Lemma 1

Lemma 1. A product of Weyl-Wigner operators remains a product under all linear
transformations of the mode operators,
$\begin{pmatrix}a^{\prime{\dagger}}&b^{\prime{\dagger}}&a^{\prime}&b^{\prime}\end{pmatrix}=[\begin{pmatrix}a^{\dagger}&b^{\dagger}&a&b\end{pmatrix}S+\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm a}\end{pmatrix}]$, the argument of which
changes according to the inverse of the corresponding symplectic
transformation, i.e.,
$\hat{T}(\alpha)\otimes\hat{T}(\beta)=\hat{T}(\alpha^{\prime})\otimes\hat{T}(\beta^{\prime})$
where
$\begin{pmatrix}\alpha^{\prime*}&\beta^{\prime*}&\alpha^{\prime}&\beta^{\prime}\end{pmatrix}=[\begin{pmatrix}\alpha^{*}&\beta^{*}&\alpha&\beta\end{pmatrix}-\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm a}\end{pmatrix}]S^{-1}$.

###### Proof.

It is sufficient to prove the theorem for the transformations of Wigner
operators. A two-mode Wigner operator is given by [8]

$\begin{split}\hat{T}(\alpha)\otimes\hat{T}(\beta)&=\int\frac{d^{2}\xi
d^{2}\zeta}{\pi^{2}}e^{\alpha\xi^{*}-\alpha^{*}\xi+\beta\zeta^{*}-\beta^{*}\zeta}\hat{D}(\xi)\otimes\hat{D}(\zeta)\\\
&=\int\frac{d^{2}\xi
d^{2}\zeta}{\pi^{2}}\exp\\{\begin{pmatrix}\alpha^{*}&\beta^{*}&\alpha&\beta\end{pmatrix}\begin{pmatrix}-\xi\\\
-\zeta\\\ \xi^{*}\\\
\zeta^{*}\end{pmatrix}-\begin{pmatrix}a^{\dagger}&b^{\dagger}&a&b\end{pmatrix}\begin{pmatrix}-\xi\\\
-\zeta\\\ \xi^{*}\\\ \zeta^{*}\end{pmatrix}\\}.\end{split}$ (76)

Now, transforming the operators via a linear transformation
$\hat{\mathcal{U}}$ and making use of its correspondence with a symplectic
group element $S$ (with $|\det S|=1$), we have

$\begin{split}&\hat{\mathcal{U}}\hat{T}(\alpha)\otimes\hat{T}(\beta)\hat{\mathcal{U}}^{\dagger}\\\
&=\int\frac{d^{2}\xi
d^{2}\zeta}{\pi^{2}}e^{\alpha\xi^{*}-\alpha^{*}\xi+\beta\zeta^{*}-\beta^{*}\zeta}\hat{\mathcal{U}}\hat{D}(\xi)\otimes\hat{D}(\zeta)\hat{\mathcal{U}}^{\dagger}\\\
&=\int\frac{d^{2}\xi
d^{2}\zeta}{\pi^{2}}\exp\\{\begin{pmatrix}\alpha^{*}&\beta^{*}&\alpha&\beta\end{pmatrix}\begin{pmatrix}-\xi\\\
-\zeta\\\ \xi^{*}\\\
\zeta^{*}\end{pmatrix}-[\begin{pmatrix}a^{\dagger}&b^{\dagger}&a&b\end{pmatrix}S+\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm a}\end{pmatrix}]\begin{pmatrix}-\xi\\\
-\zeta\\\ \xi^{*}\\\ \zeta^{*}\end{pmatrix}\\}\\\ &=\int\frac{d^{2}\xi
d^{2}\zeta}{\pi^{2}}\exp\\{[\begin{pmatrix}\alpha^{*}&\beta^{*}&\alpha&\beta\end{pmatrix}-\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm
a}\end{pmatrix}]S^{-1}S\begin{pmatrix}-\xi\\\ -\zeta\\\ \xi^{*}\\\
\zeta^{*}\end{pmatrix}-\begin{pmatrix}a^{\dagger}&b^{\dagger}&a&b\end{pmatrix}S\begin{pmatrix}-\xi\\\
-\zeta\\\ \xi^{*}\\\ \zeta^{*}\end{pmatrix}\\}\\\ &=\frac{1}{|\det
S|}\int\frac{d^{2}\xi^{\prime}d^{2}\zeta^{\prime}}{\pi^{2}}e^{\alpha^{\prime}\xi^{\prime*}-\alpha^{\prime*}\xi^{\prime}+\beta^{\prime}\zeta^{\prime*}-\beta^{\prime*}\zeta^{\prime}}\hat{D}(\xi^{\prime})\otimes\hat{D}(\zeta^{\prime})\\\
&=\hat{T}(\alpha^{\prime})\otimes\hat{T}(\beta^{\prime}),\end{split}$ (77)

in which

$\begin{pmatrix}-\xi^{\prime}\\\ -\zeta^{\prime}\\\ \xi^{\prime*}\\\
\zeta^{\prime*}\end{pmatrix}=S\begin{pmatrix}-\xi\\\ -\zeta\\\ \xi^{*}\\\
\zeta^{*}\end{pmatrix},\quad\begin{pmatrix}\alpha^{\prime*}&\beta^{\prime*}&\alpha^{\prime}&\beta^{\prime}\end{pmatrix}=[\begin{pmatrix}\alpha^{*}&\beta^{*}&\alpha&\beta\end{pmatrix}-\begin{pmatrix}D_{\rm
a}^{*}&D_{\rm b}^{*}&D_{\rm a}&D_{\rm a}\end{pmatrix}]S^{-1}$ (78)

Clearly, this result also applies to single mode linear transformation. We
also bare in mind that $d^{2}\alpha=d{\rm Im}\alpha\times d{\rm
Re}\alpha=\frac{1}{2}d\alpha d\alpha^{*}=\frac{1}{2}dqdp$. ∎

##  Appendix D Order shift relation

Here we show that the overlap of two quasiprobability distributions with order
mismatch $\Delta$ is invariant under exchange of the mismatch between the
arguments of either distribution. We proceed as follows.

$\begin{split}A&=\int\frac{d^{2}\alpha}{\pi}\mathscr{W}_{\hat{\Lambda}}(\alpha;s+\Delta)\mathscr{W}_{\hat{\Upsilon}}(\alpha;t)\\\
&=\int\frac{d^{2}\alpha}{\pi}{\rm
Tr}\hat{\Lambda}\hat{T}(\alpha;s+\Delta)\cdot{\rm
Tr}\hat{\Upsilon}\hat{T}(\alpha;t)\\\
&=\int\frac{d^{2}\alpha}{\pi}\left[\frac{2}{-\Delta}\int\frac{d^{2}\xi}{\pi}e^{-\frac{2|\alpha-\xi|^{2}}{-\Delta}}{\rm
Tr}\hat{\Lambda}\hat{T}(\xi;s)\right]\cdot{\rm
Tr}\hat{\Upsilon}\hat{T}(\alpha;t)\\\ &=\int\frac{d^{2}\xi}{\pi}{\rm
Tr}\hat{\Lambda}\hat{T}(\xi;s)\cdot\left[\frac{2}{-\Delta}\int\frac{d^{2}\alpha}{\pi}e^{-\frac{2|\alpha-\xi|^{2}}{-\Delta}}{\rm
Tr}\hat{\Upsilon}\hat{T}(\alpha;t)\right]\\\ &=\int\frac{d^{2}\xi}{\pi}{\rm
Tr}\hat{\Lambda}\hat{T}(\xi;s)\cdot{\rm
Tr}\hat{\Upsilon}\hat{T}(\xi;t+\Delta).\end{split}$ (79)

In the third and fourth lines we have used the relation [8]

$\hat{T}(\alpha;s+\Delta)=\frac{2}{-\Delta}\int\frac{d^{2}\xi}{\pi}e^{-\frac{2|\alpha-\xi|^{2}}{-\Delta}}\hat{T}(\xi;s).$
(80)

##  Appendix E The case of _s_ -parameterized tomograms

In the main text we have shown how a single Wigner-function tomogram (i.e.
$s=0$) of the mechanical state—obtained using our method in the strong
interaction regime—can be used to detect the state’s P-function
nonclassicality. In the case of weaker interactions one generally obtains
smoothed $s$-parametrized tomograms of the mechanical state from which a
$s$-parametrized quasiprobability distribution can be reconstructed. We will
now generalize the criterion of Park et al [28] to $s$-parametrized tomograms
and show that even with weak interactions, one can verify the P-function
nonclassicality of mechanical states using our method.

Rewriting Eq. (21), in terms of $s$-parameterized quasiprobability
distributions using Eq. (30) we have

$W_{\hat{\varrho}}(\beta;s)=\int\frac{d^{2}\alpha}{\pi}P(\alpha)W_{\alpha}(\beta;s),$
(81)

where $W_{\alpha}(\beta;s)=2\exp\\{-2|\beta-\alpha|^{2}/(1-s)\\}/(1-s)$
represents the $s$-parameterized distribution corresponding to the coherent
state $|\alpha\rangle$. We now make the observation that, in the position-
momentum coordinates,

$\begin{split}W_{\alpha}(q,p;s)=w_{\alpha}(q,0;s)w_{\alpha}(p,\frac{\pi}{2};s)\equiv
W_{q_{0},p_{0}}(q,p;s),\end{split}$ (82)

where $w_{\alpha}(q,0;s)=e^{\frac{-(q-q_{0})^{2}}{1-s}}/\sqrt{\pi(1-s)}$ and
$w_{\alpha}(p,\frac{\pi}{2};s)=e^{\frac{-(p-p_{0})^{2}}{1-s}}/\sqrt{\pi(1-s)}$
are the two $s$-parameterized tomograms of the coherent state $|\alpha\rangle$
in which $q_{0}=\sqrt{2}{\rm Re}\alpha$ and $p_{0}=\sqrt{2}{\rm Im}\alpha$.
Note that, without loss of generality, we have assumed $\phi=0$. Substituting
Eq. (82) into Eq. (81), we obtain

$\begin{split}W_{\hat{\varrho}}(q,p;s)&=\int
dq^{\prime}dp^{\prime}P(q^{\prime},p^{\prime})W_{q^{\prime},p^{\prime}}(q,p;s)\\\
&=\frac{1}{\pi(1-s)}\int
dq^{\prime}dp^{\prime}P(q^{\prime},p^{\prime})e^{\frac{-(q-q^{\prime})^{2}}{1-s}}e^{\frac{-(q-q^{\prime})^{2}}{1-s}}.\end{split}$
(83)

Consequently, the $s$-parameterized tomogram of the state $\hat{\varrho}$,
using Eqs. (81) and (82), is given by

$\begin{split}w(q,0;s)&=\frac{1}{\sqrt{\pi(1-s)}}\int dq^{\prime}\left(\int
dp^{\prime}P(q^{\prime},p^{\prime})\right)e^{\frac{-(q-q^{\prime})^{2}}{1-s}}\\\
&=\frac{1}{\sqrt{\pi(1-s)}}\int
dq^{\prime}F(q^{\prime})e^{\frac{-(q-q^{\prime})^{2}}{1-s}},\end{split}$ (84)

with $F(q^{\prime})=\int dp^{\prime}P(q^{\prime},p^{\prime})$. Importantly, if
the state $\hat{\varrho}$ is a classical state, then
$P(q^{\prime},p^{\prime})\geqslant 0$ is a function with finite support and so
is $F(q^{\prime})$. Now, suppose that

$w_{\rm f}(p,\frac{\pi}{2};s)=\frac{1}{\sqrt{\pi(1-s)}}\int
dp^{\prime}G(p^{\prime})e^{\frac{-(p-p^{\prime})^{2}}{1-s}},$ (85)

where $G(p^{\prime})\geqslant 0$ is a positive function with a finite support.
We proceed by defining the fictitious $s$-parameterized distribution
$w(x,0;s)\mapsto W_{\rm f}(q,p;s)=w(x,0;s)w_{\rm f}(p,\frac{\pi}{2};s)$.
Hence,

$\begin{split}W_{\rm f}(q,p;s)=\frac{1}{\pi(1-s)}\int
dq^{\prime}dp^{\prime}F(q^{\prime})G(p^{\prime})e^{\frac{-[(q-q^{\prime})^{2}+(p-p^{\prime})^{2}]}{1-s}}.\end{split}$
(86)

A simple comparison between Eqs. (86) and (81) shows that $W_{\rm f}(q,p;s)$
corresponds to a bona fide P-classical quantum state, whenever the state
$\hat{\varrho}$ is classical. This is because, in this case, the new function
$P_{\rm f}(q^{\prime},p^{\prime})=F(q^{\prime})G(p^{\prime})$ is a positive
(normalized) function with finite support. Therefore, it represents a valid
P-function. Consequently, whenever $W_{\rm f}(q,p;s)$ fails to represent a
bona fide fictitious quantum state, we conclude that the initial state
$\hat{\varrho}$ was necessarily nonclassical.

Two choices for the fictitious tomogram $w_{\rm f}(p,\frac{\pi}{2};s)$, as
presented in the main text in Eqs. (25) and (26), are

  1. (i)

that of a vacuum state $w_{\rm
f}(p,\frac{\pi}{2};s)=w_{0}(p,\frac{\pi}{2};s)=e^{\frac{-p^{2}}{1-s}}/\sqrt{\pi(1-s)}$
corresponding to the map $w(x,0;s)\mapsto W_{\rm
f}(q,p;s)=w(x,0;s)w_{0}(p,\frac{\pi}{2};s)$; and,

  2. (ii)

the tomogram of the state itself, i.e., $w_{\rm
f}(p,\frac{\pi}{2};s)=w(p,0;s)$ corresponding to $(w(x,0;s),w(x,0;s))\mapsto
W_{\rm f}(q,p;s)=w(x,0;s)w(p,0;s)$.

To check the legitimacy of the fictitious density operators resulting from
these demarginalization maps, one may follow the standard arguments provided
in Ref. [28], such as the Kastler-Loupias-Miracle-Sole test [59, 60, 61]. In
the general formalism, to perform such tests one may also require evaluation
of various Fock basis elements of the fictitious density matrix. To this end,
one should use the trace relation of Eq. (31) and the corresponding
$s$-parametrized Fock basis projections given by [62]

$W_{|n\rangle\langle
m|}(\alpha;s)=(\frac{2}{1-s})^{m+n+1}\frac{1}{\sqrt{n!m!}}h_{n,m}(\alpha^{*},\alpha|\frac{s^{2}-1}{4})e^{-\frac{2|\alpha|^{2}}{1-s}},$

in which $h_{n,m}(x,y|\epsilon)$ are the incomplete 2D Hermite polynomials
defined as [63]

$h_{n,m}(x,y|\epsilon)=\sum_{i=0}^{\min\\{m,n\\}}\binom{n}{i}\binom{m}{i}i!\epsilon^{i}x^{n-i}y^{m-i}.$

##  Appendix F The effect of classical noise and losses on nonclassicality
certification

As discussed in the main text any noise and losses _before_ the optomechanical
interaction can be accounted for by calibrating the input optical field. We
now focus on the readout channel _after_ the interaction and show that all
noise and losses commute with the interaction and can thus be taken into
account as before. We start by rewriting the Gaussian optical kernel on the
r.h.s. of Eq. (12) as

$\mathscr{G}_{\rm o}({\bf x}^{\prime},{\bf
x}_{0};{\bf\Sigma})=\frac{1}{\pi}\exp\\{-\frac{1}{2}({\bf x^{\prime}}-{\bf
x}_{0}){\bf\Sigma}^{-1}({\bf x^{\prime}}-{\bf x}_{0})^{{}_{\rm T}}\\},$ (87)

where ${\bf x^{\prime}}=(x^{\prime},p^{\prime})$, ${\bf
x}_{0}=(0,X(\varphi+\varphi_{\rm d}+\frac{\pi}{2}))$, and
${\bf\Sigma}^{-1}={\rm
diag}(e^{-2\varepsilon},4\chi^{2}u^{2}e^{2\varepsilon})$ represents the
covariance matrix, so that

$\begin{split}\mathscr{W}_{\rm out;o}(x^{\prime},2\chi
up^{\prime};0)&=\frac{1}{\pi}e^{-e^{-2\varepsilon}x^{\prime
2}}\int\frac{d^{2}\beta}{\pi}e^{-4\chi^{2}u^{2}e^{2\varepsilon}[p^{\prime}-X(\varphi+\varphi_{\rm
d}+\frac{\pi}{2})]^{2}}\mathscr{W}_{\rm m}(\beta;0)\\\
&=\int\frac{d^{2}\beta}{\pi}\mathscr{G}_{\rm o}({\bf x}^{\prime},{\bf
x}_{0};{\bf\Sigma})\mathscr{W}_{\rm m}(\beta;0).\end{split}$ (88)

Since all classical noise and losses acting on the output optical field
correspond to Gaussian channels [29], they can be described as the action of a
Gaussian kernel $\mathscr{G}^{\rm noise}({\bf\tilde{x}},{\bf
x}^{\prime};\boldsymbol{\sigma})$ centred at $(x^{\prime},p^{\prime})$, with
covariance matrix $\boldsymbol{\sigma}$. Applying the this kernel to the
l.h.s. of Eq. (12) gives

$\begin{split}\mathscr{W}^{\rm noisy}_{\rm out;o}(\tilde{x},\tilde{p};0)&=\int
dx^{\prime}dp^{\prime}\mathscr{G}^{\rm noise}({\bf\tilde{x}},{\bf
x}^{\prime};\boldsymbol{\sigma})\mathscr{W}_{\rm out;o}(x^{\prime},2\chi
up^{\prime};0),\\\ &=\int dx^{\prime}dp^{\prime}\mathscr{G}^{\rm
noise}({\bf\tilde{x}},{\bf
x}^{\prime};\boldsymbol{\sigma})\Big{[}\int\frac{d^{2}\beta}{\pi}\mathscr{G}_{\rm
o}({\bf x}^{\prime},{\bf x}_{0};{\bf\Sigma})\mathscr{W}_{\rm
m}(\beta;0)\Big{]}.\end{split}$ (89)

It is now easy to see that noise and losses acting on the output optical field
commute with the optomechanical interaction and can thus be considered to act
on the input field. Formally, we change the order of integrations in Eq. (89)
to obtain

$\begin{split}\mathscr{W}^{\rm noisy}_{\rm
out;o}(\tilde{x},\tilde{p};0)&=\int\frac{d^{2}\beta}{\pi}\mathscr{W}_{\rm
m}(\beta;0)\Big{[}\int dx^{\prime}dp^{\prime}\mathscr{G}^{\rm
noise}({\bf\tilde{x}},{\bf x}^{\prime};\boldsymbol{\sigma})\mathscr{G}_{\rm
o}({\bf x}^{\prime},{\bf x}_{0};{\bf\Sigma})\Big{]}\\\
&=\int\frac{d^{2}\beta}{\pi}\mathscr{G}^{\rm noise}_{{}_{\rm
o}}({\bf\tilde{x}},{\bf
x}_{0};\boldsymbol{\sigma}+{\bf\Sigma})\mathscr{W}_{\rm
m}(\beta;0).\end{split}$ (90)

The latter equality follows from the fact that the expression in the square
brackets is a convolution of a Gaussian wave packet with a Gaussian kernel
which results in a Gaussian wave packet. It is now easy to see that by
appropriately rotating the coordinates $(\tilde{x},\tilde{p})$ to diagonalise
the covariance matrix $\boldsymbol{\sigma}+{\bf\Sigma}$, one obtains a
relation equivalent to Eq. (12) and can proceed as before. This shows that
classical noise and losses do not distort the picture in our general
formalism, assuming that a full characterization of the readout channel is
possible.

We would like to emphasize, however, that careful characterization of the
readout channel is crucial for _only one_ of our nonclassicality criteria,
namely the first demarginalization map in Eq. (25). Recall that any Gaussian
convolution of the form of Eq. (90) is equivalent to a degradation of a given
$s$-parameterized phase-space distribution to another distribution with
$s^{\prime}\leqslant s$. Consequently, since the criterion in Eq. (25),
requires that the $s$ parameter of the fictitious tomogram $w_{\rm
f}(p,\frac{\pi}{2};s)$ matches that of the measured tomogram $w_{\rm
f}(x,0;s)$, precise knowledge of the actual experimental value of $s^{\prime}$
is important. Relaxing this requirement somewhat, it is also possible to use a
fictitious tomogram with $s^{\prime}\leqslant s$ (i.e. an overestimation of
the actual noise), because, using Eq. (85), one possible fictitious tomogram
$w_{\rm f}(p,\frac{\pi}{2};s^{\prime})$ is given by

$w_{\rm f}(p,\frac{\pi}{2};s^{\prime})=\frac{1}{\sqrt{\pi(1-s^{\prime})}}\int
dp^{\prime}\mathscr{G}^{\rm
noise}(p^{\prime},0;\sigma)e^{\frac{-(p-p^{\prime})^{2}}{1-s^{\prime}}},$ (91)

where $\sigma=\frac{s-s^{\prime}}{(1-s^{\prime})(1-s)}$ and $\mathscr{G}^{\rm
noise}(p^{\prime},0;\sigma)$ is a Gaussian with finite support provided that
$s^{\prime}\leqslant s$. Consequently, in the presence of noise and loss, one
should be careful in using the nonclassicality criterion in Eq. (25).

Importantly, such a caution is not needed when using the second
demarginalization map, Eq. (26). In this case, the experimentally measured
tomogram is used for both quadratures and the ordering parameters are thus
always matched, regardless of the actual experimental value of $s$. As a
result, our second nonclassicality criteria is immune to void nonclassicality
detections due to noise and losses.

## References

  * Aspelmeyer  _et al._ [2014a] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt,  _Cavity optomechanics,_Rev. Mod. Phys. 86, 1391 (2014a).
  * Bose  _et al._ [1999] S. Bose, K. Jacobs, and P. L. Knight,  _Scheme to probe the decoherence of a macroscopic object,_Phys. Rev. A 59, 3204 (1999).
  * Marshall  _et al._ [2003] W. Marshall, C. Simon, R. Penrose, and D. Bouwmeester,  _Towards Quantum Superpositions of a Mirror,_Phys. Rev. Lett. 91, 130401 (2003).
  * Pikovski  _et al._ [2012] I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. S. Kim, and Č. Brukner,  _Probing Planck-scale physics with quantum optics,_Nat. Phys. 8, 393 (2012).
  * Rugar  _et al._ [2004] D. Rugar, R. Budakian, H. J. Mamin, and B. W. Chui,  _Single spin detection by magnetic resonance force microscopy,_Nature 430, 329 (2004).
  * Vanner  _et al._ [2013] M. R. Vanner, J. Hofer, G. D. Cole, and M. Aspelmeyer,  _Cooling-by-measurement and mechanical state tomography via pulsed optomechanics._Nat. Commun. 4, 2295 (2013).
  * Ringbauer  _et al._ [2016] M. Ringbauer, T. J. Weinhold, A. G. White, and M. R. Vanner,  _Generation of Mechanical Interference Fringes by Multi-Photon Quantum Measurement,_Preprint at arXiv:1602.05955 (2016).
  * Glauber [2007] R. J. Glauber,  _Quantum theory of optical coherence : selected papers and lectures_ (Wiley-VCH, 2007) p. 639.
  * Vogel and Welsch [2006] W. Vogel and D.-G. Welsch,  _Quantum Optics_ (Wiley-VCH, 2006).
  * Lvovsky  _et al._ [2001] A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller,  _Quantum State Reconstruction of the Single-Photon Fock State,_Phys. Rev. Lett. 87, 050402 (2001).
  * Zavatta  _et al._ [2004] A. Zavatta, S. Viciani, and M. Bellini,  _Tomographic reconstruction of the single-photon Fock state by high-frequency homodyne detection,_Phys. Rev. A 70, 053821 (2004).
  * Leonhardt [1997] U. Leonhardt,  _Cambridge Studies in Modern Optics_, edited by P. Knight and A. Miller (Cambridge University Press, 1997).
  * Hansen  _et al._ [2001] H. Hansen, T. Aichele, C. Hettich, P. Lodahl, A. I. Lvovsky, J. Mlynek, and S. Schiller,  _Ultrasensitive pulsed, balanced homodyne detector: application to time-domain quantum measurements._Opt. Lett. 26, 1714 (2001).
  * Vanner  _et al._ [2015] M. R. Vanner, I. Pikovski, and M. S. Kim,  _Towards optomechanical quantum state reconstruction of mechanical motion,_Ann. Phys. 527, 15 (2015).
  * Vanner  _et al._ [2011] M. R. Vanner, I. Pikovski, G. D. Cole, M. S. Kim, C. Brukner, K. Hammerer, G. J. Milburn, and M. Aspelmeyer,  _Pulsed quantum optomechanics,_Proc. Natl. Acad. Sci. 108, 16182 (2011).
  * Aspelmeyer  _et al._ [2014b] M. Aspelmeyer, T. J. Kippenberg, and F. Marquard,  _Cavity Optomechanics_, edited by M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Vol. 86 (Springer Berlin Heidelberg, 2014) pp. 1391–1452.
  * Wang  _et al._ [1994] X.-B. Wang, S.-X. Yu, and Y.-D. Zhang,  _Linear quantum transformation and normal product calculation of boson exponential quadratic operators,_J. Phys. A. 27, 6563 (1994).
  * Zhang and Tang [1994] Y.-d. Zhang and Z. Tang,  _General theory of linear quantum transformation of Bargmann-Fock space,_Nuovo Cim. B Ser. 11 109, 387 (1994).
  * Wallentowitz and Vogel [1996] S. Wallentowitz and W. Vogel,  _Unbalanced homodyning for quantum state measurements,_Phys. Rev. A 53, 4528 (1996).
  * Kiesel  _et al._ [2008] T. Kiesel, W. Vogel, V. Parigi, A. Zavatta, and M. Bellini,  _Experimental determination of a nonclassical Glauber-Sudarshan P-function,_Phys. Rev. A 78, 021804 (2008).
  * Kiesel  _et al._ [2011] T. Kiesel, W. Vogel, B. Hage, and R. Schnabel,  _Direct Sampling of Negative Quasiprobabilities of a Squeezed State,_Phys. Rev. Lett. 107, 113604 (2011).
  * Agudelo  _et al._ [2015] E. Agudelo, J. Sperling, W. Vogel, S. Köhnke, M. Mraz, and B. Hage,  _Continuous sampling of the squeezed-state nonclassicality,_Phys. Rev. A 92, 033837 (2015).
  * Shchukin  _et al._ [2005] E. Shchukin, T. Richter, and W. Vogel,  _Nonclassicality criteria in terms of moments,_Phys. Rev. A 71, 011802 (2005).
  * Shchukin and Vogel [2005] E. V. Shchukin and W. Vogel,  _Nonclassical moments and their measurement,_Phys. Rev. A 72, 043808 (2005).
  * Kiesel  _et al._ [2009] T. Kiesel, W. Vogel, B. Hage, J. DiGuglielmo, A. Samblowski, and R. Schnabel,  _Experimental test of nonclassicality criteria for phase-diffused squeezed states,_Phys. Rev. A 79, 022122 (2009).
  * Kiesel and Vogel [2012] T. Kiesel and W. Vogel,  _Complete nonclassicality test with a photon-number-resolving detector,_Phys. Rev. A 86, 032119 (2012).
  * Ryl  _et al._ [2015] S. Ryl, J. Sperling, E. Agudelo, M. Mraz, S. Köhnke, B. Hage, and W. Vogel,  _Unified nonclassicality criteria,_Phys. Rev. A 92, 011801 (2015).
  * Park  _et al._ [2017] J. Park, Y. Lu, J. Lee, Y. Shen, K. Zhang, S. Zhang, M. S. Zubairy, K. Kim, and H. Nha,  _Revealing nonclassicality beyond Gaussian states via a single marginal distribution,_Proc. Natl. Acad. Sci. 114, 891 (2017).
  * Eisert and Wolf [2007] J. Eisert and M. M. Wolf, in  _Quantum Inf. with Contin. Var. Atoms Light_ (Imperial College Press, London, 2007) pp. 23–42.
  * Parthasarathy and Sengupta [2015] K. R. Parthasarathy and R. Sengupta,  _From particle counting to gaussian tomography,_Infinite Dimensional Analysis, Quantum Probability and Related Topics 18, 1550023 (2015).
  * Thorne  _et al._ [1978] K. S. Thorne, R. W. P. Drever, C. M. Caves, M. Zimmermann, and V. D. Sandberg,  _Quantum nondemolition measurements of harmonic oscillators,_Phys. Rev. Lett. 40, 667 (1978).
  * Braginsky  _et al._ [1980] V. B. Braginsky, Y. I. Vorontsov, and K. S. Thorne,  _Quantum nondemolition measurements,_Science 209, 547 (1980).
  * Ferraro and Paris [2012] A. Ferraro and M. G. a. Paris,  _Nonclassicality Criteria from Phase-Space Representations and Information-Theoretical Constraints Are Maximally Inequivalent,_Phys. Rev. Lett. 108, 260403 (2012).
  * Kleckner  _et al._ [2011] D. Kleckner, B. Pepper, E. Jeffrey, P. Sonin, S. M. Thon, and D. Bouwmeester,  _Optomechanical trampoline resonators._Opt. Expr. 19, 19708 (2011).
  * Norte  _et al._ [2016] R. A. Norte, J. P. Moura, and S. Gröblacher,  _Mechanical Resonators for Quantum Optomechanics Experiments at Room Temperature,_Phys. Rev. Lett. 116, 147202 (2016).
  * Gröblacher  _et al._ [2009] S. Gröblacher, K. Hammerer, M. R. Vanner, and M. Aspelmeyer,  _Observation of strong coupling between a micromechanical resonator and an optical cavity field,_Nature 460, 724 (2009).
  * Becerra  _et al._ [2013] F. E. Becerra, J. Fan, and A. Migdall,  _Implementation of generalized quantum measurements for unambiguous discrimination of multiple non-orthogonal coherent states,_Nat. Commun. 4, 2028 (2013).
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer,  _Continuous-variable optical quantum-state tomography,_Rev. Mod. Phys. 81, 299 (2009).
  * Barbieri  _et al._ [2010] M. Barbieri, N. Spagnolo, M. G. Genoni, F. Ferreyrol, R. Blandino, M. G. A. Paris, P. Grangier, and R. Tualle-Brouri,  _Non-Gaussianity of quantum states: An experimental test on single-photon-added coherent states,_Phys. Rev. A 82, 063833 (2010).
  * Murch  _et al._ [2008] K. W. Murch, K. L. Moore, S. Gupta, and D. M. Stamper-Kurn,  _Observation of quantum-measurement backaction with an ultracold atomic gas,_Nat. Phys. 4, 561 (2008).
  * Thompson  _et al._ [2008] J. D. Thompson, B. M. Zwickl, A. M. Jayich, F. Marquardt, S. M. Girvin, and J. G. E. Harris,  _Strong dispersive coupling of a high-finesse cavity to a micromechanical membrane,_Nature 452, 72 (2008).
  * Anguiano  _et al._ [2016] S. Anguiano, A. E. Bruchhausen, B. Jusserand, I. Favero, F. R. Lamperti, L. Lanco, I. Sagnes, A. Lemaître, N. D. Lanzillotti-Kimura, P. Senellart, and A. Fainstein,  _Time-Resolved Cavity Nano-Optomechanics in the 20-100 GHz range,_Preprint at arXiv:1610.04179 (2016).
  * Arcizet  _et al._ [2006] O. Arcizet, P.-F. Cohadon, T. Briant, M. Pinard, and A. Heidmann,  _Radiation-pressure cooling and optomechanical instability of a micromirror,_Nature 444, 71 (2006).
  * Cuthbertson  _et al._ [1996] B. D. Cuthbertson, M. E. Tobar, E. N. Ivanov, and D. G. Blair,  _Parametric back-action effects in a high-Q cyrogenic sapphire transducer,_Rev. Sci. Instr. 67, 2435 (1996).
  * 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,  _Laser cooling of a nanomechanical oscillator into its quantum ground state,_Nature 478, 89 (2011).
  * Verhagen  _et al._ [2012] E. Verhagen, S. Deléglise, S. Weis, a. Schliesser, and T. J. Kippenberg,  _Quantum-coherent coupling of a mechanical oscillator to an optical cavity mode._Nature 482, 63 (2012).
  * Teufel  _et al._ [2011] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds,  _Sideband cooling of micromechanical motion to the quantum ground state,_Nature 475, 359 (2011).
  * Husimi [1940] K. Husimi,  _Some Formal Properties of the Density Matrix,_Proc. Phys. Math. Soc. Japan 3rd Series 22, 264 (1940).
  * Kano [1965] Y. Kano,  _A New Phase-Space Distribution Function in the Statistical Theory of the Electromagnetic Field,_J. Math. Phys. 6, 1913 (1965).
  * Wigner [1932] E. Wigner,  _On the Quantum Correction For Thermodynamic Equil,_ Phys. Rev. 40, 749 (1932).
  * Glauber [1963] R. J. Glauber,  _The Quantum Theory of Optical Coherence,_Phys. Rev. 130, 2529 (1963).
  * Sudarshan [1963] E. C. G. Sudarshan,  _Equivalence of Semiclassical and Quantum Mechanical Descriptions of Statistical Light Beams,_Phys. Rev. Lett. 10, 277 (1963).
  * Sperling [2016] J. Sperling,  _Characterizing maximally singular phase-space distributions,_Phys. Rev. A 94, 013814 (2016).
  * Vogel [2000] W. Vogel,  _Nonclassical States: An Observable Criterion,_Phys. Rev. Lett. 84, 1849 (2000).
  * Sperling  _et al._ [2017] J. Sperling, W. R. Clements, A. Eckstein, M. Moore, J. J. Renema, W. S. Kolthammer, S. W. Nam, A. Lita, T. Gerrits, W. Vogel, G. S. Agarwal, and I. A. Walmsley,  _Detector-independent verification of quantum light,_Phys. Rev. Lett. 118, 163602 (2017).
  * Schleich [2015] W. P. Schleich,  _Quantum optics in phase space_ (Wiley-VCH, 2015).
  * Bazrafkan and Nahvifard [2011] M. R. Bazrafkan and E. Nahvifard,  _Quantum-state tomogram from s-parameterized quasidistributions,_J. Russ. Laser Res. 32, 230 (2011).
  * Puri [2001] R. R. Puri,  _Mathematical Methods of Quantum Optics_, edited by W. T. R. Metz (Springer Berlin Heidelberg, 2001).
  * Kastler [1965] D. Kastler,  _The C*-algebras of a free boson field,_Commun. Math. Phys. 1, 14 (1965).
  * Loupias and Miracle-Sole [1966] G. Loupias and S. Miracle-Sole,  _C*-algèbres des systèmes canoniques. i,_Commun. Math. Phys. 2, 31 (1966).
  * Nha [2008] H. Nha,  _Complete conditions for legitimate wigner distributions,_Phys. Rev. A 78, 012103 (2008).
  * Shahandeh and Bazrafkan [2013] F. Shahandeh and M. R. Bazrafkan,  _The general boson ordering problem and its combinatorial roots,_Phys. Scr. T153, 014056 (2013).
  * Dattoli [2003] G. Dattoli,  _Incomplete 2D Hermite polynomials: properties and applications,_J. Math. Anal. Appl. 284, 447 (2003).