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

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.

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]. Such techniques thus 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 such 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. One of the distinctive features of quantum system, however, 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 sparameterized 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, correspond- * f.shahandeh@uq.edu.au ing 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ˆ 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]. Such measurements, however, are not directly available for mechanical resonators. Instead, the common approach is to use a high-quality optical cavity in the resolved sideband regime, which would allow to transfer the mechanical state onto an optical field, and then use established techniques to reconstruct the optical state. However, this regime is experimentally very difficult to reach [14], and Wigner-function tomography of a mechanical system thus remains an outstanding challenge.
Here, we introduce an alternative approach which works outside the demanding resolved sideband regime with no, or only a low quality cavity. We show that, even in this regime, where the state transfer in the usual sense is not available, it is possible to perform a tomographic reconstruction of the s-parameterized quasiprobability distributions, including the Wigner function of the mechanical motion. Interestingly, our approach is independent of the single-photon coupling strength which makes it suitable for a wide range of optomechanical systems with current technology. Our approach works by imprinting mechanical quadrature distributions on a strong optical coherent pulse, from which it can then be extracted using standard optical methods. 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.
The paper is organised as follows. In Sec. I we provide a self-contained overview of phase-space quasiprobability dis-tributions, their tomographic reconstruction, and the optomechanical interaction. We then discuss the phase-space description of the optomechanical interaction in Sec. II. 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. In Sec. III, we provide a necessary condition for demonstrating mechanical nonclassicality using the smallest subset of data obtained from our scheme, namely, a single mechanical tomogram. Finally, in Sec. IV we describe an explicit experimental protocol and show that it can be implemented with current technology in a variety of optomechanical systems.

I. BACKGROUND
A. 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 operatorsT (α; s) = , for a complex-valued phase-space point α ∈ C and distribution parameter s ∈ [−1, 1] [8]. Here,D(ξ)= exp{ξa † − ξ * a} is the displacement operator by the value ξ ∈ C, and a (a † ) is the bosonic annihilation (creation) operator. Weyl-Wigner operators are normal, complete, and satisfy the orthogonality (duality) relation, TrT (α; s)T (β; −s)=πδ (2) respectively. As a result, we can expand any operatorΛ in this basis asΛ The distribution WΛ(α; s) = TrΛT (α; s) is known as the s-parameterized quasiprobability distribution of the operator Λ, which includes as special cases [8] the Husimi-Kano Qfunction [15,16] for s = −1, the Wigner [17] function for s = 0, and the Glauber-Sudarshan P-function [18,19] 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 or even be singular. If the operatorΛ is a quantum state, then such negative values are the primary signature of nonclassical behaviour [9,20], and can always be detected using direct observation via tomographic methods [12,[21][22][23] or nonclassicality criteria [24][25][26][27][28][29][30]. Using Eq. (4) and the duality relation of Eq. (3), the trace operation can be expressed in terms of phase-space distribu-tions as Finally, it is sometimes more natural to work with positionmomentum coordinates, in which case d 2 α = dReα dImα = dqdp 2 , and WΛ(α; s) = 2πWΛ(q, p; s), so that B. 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 φ-quadrature of a bosonic field is given bŷ 2 with eigenvalues and eigenvectors defined by the eigenvalue equationx(φ)|x(φ) = x(φ)|x(φ) . The set of all projections onto quadrature eigenstates {Π(x, φ) = |x(φ) x(φ)|} for x ∈ R and φ ∈ [0, 2π), known as tomographic operators, forms an informationally complete operator basis [12,31]. One can then define the set of tomograms {w(x, φ)} for a quantum stateˆ as Using Eq. (5), this tomographic relation can be expressed in terms of Wigner functions as, where Wˆ (α; 0) is the Wigner function of the quantum state, and is the Wigner function of the tomographic operatorΠ(x, φ). Equation (8), also known as a Radon transformation, can be interpreted as sampling the Wigner function of a quantum stateˆ with strings of zero width, see Fig. 1a. In turn, a sufficiently large set of tomograms can be used to reconstruct the Wigner function via an inverse Radon transformation, see Fig. 2. In short, The tomographic relation (7) can be generalized to arbitrary phase-space quasiprobability distributions with s 0 as [32] w(x, φ) = with the tomographic operator (12) In analogy with the Radon transform, Eq. (11) can be interpreted as sampling an s-parameterized phase-space quasiprobability distribution with a string with Gaussian cross section of finite width, see Fig. 1b. Notably, Eqs. (8) and (11) connect the same experimentally measured tomograms w(x, φ) to different quasiprobability representations of the underlying quantum state.
(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.

C. 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 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 ω o , the free evolution of the mechanics with frequency ω m , and the radiation pressure interaction. The latter couples the mechanical position operatorX ∝ b † + b to the optical photon number operatorN = a † a with the single photon coupling strength g 0 = Gx zpf , where G is the coupling constant, and x zpf = /2m eff ω 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 κ 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 the non-resonant terms, which is known as the rotating wave approximation, the optomechanical interaction is then dominated, either by a beam-splitterlike interaction for the red sideband, or a squeezing-like interaction for the blue sideband [33]. For tomographic purposes, 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.
Here, we are interested in the technically simpler scenario in which the optical field is a free-space mode 1 . 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. (13) must be considered. Suppose that the mechanical system is initially in a stateˆ m which then interacts with a displaced optical pulse,D(α)ˆ oD † (α) over a duration τ . At a time τ 0 after the interaction, the optical field is displaced by −α. The evolution of the state is thus given byD As we show in detail in Appendix A, under the condition that the optical field is strong (i.e. |α| 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Û oÛkÛom , whereÛ om ,Û k , andÛ o are unitaries representing the optomechanical interaction, an optical Kerr interaction, and a displaced optical free evolution, respectively, given bŷ U k = e iv[4a † a+2r(e iθ a † +e −iθ a)+(e 2iθ a †2 +e −2iθ a 2 )] , Here, we have defined the parameters as follows: α = re iθ with r ∈ R, µ = 1 − e iωmτ = ue iϕ with and It is important to note that, while the phase of the optical mode θ is being set by the input optical coherent state, the relevant mechanical phase ϕ is determined by the length of the input optical pulse τ .

A. Mode transformation
By defining a vector of input mode operatorsˆ A := (a † , b † , a, b), we aim to find the the output mode operators in the Heisenberg picture under the aforementioned unitary evo- . Since these unitary evolutions are single-mode and two-mode linear transformations, they correspond to elements of the symplectic groups Sp(2, R) and Sp(4, R), respectively. We can thus exploit the properties of symplectic transformations, in particular, that if [34,35]. Using Eqs. (14), (15), and (16), after some algebra and making reasonable restrictions on the choice of the optical pulse τ and the delay duration τ 0 , via the transformed optomechanical mode operators are found as Here,X(ϕ) = (b † e iϕ + be −iϕ )/ √ 2 andx = (a † + a)/ √ 2 are the ϕ-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. (20) and (21). Importantly, in Eq. (20) we require that 0 ω m τ 2π, both to obtain a nonzero coupling (u = 0) and to minimize the effect of mechanical losses. This can be achieved by assuming the pulse duration τ to be on the order of half a mechanical period, i.e., ω m τ ≈ π. It is evident from Eq. (21) 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 Wignerfunction 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, , the argument of which change according to the inverse of the corresponding symplectic transformation, i.e.,

B. Interaction in phase space
We are now ready to examine the optomechanical interactions in the phase-space representation. The Wigner-function of the joint initial optomechanical state, using Eq. (4), iŝ where W o and W m are the Wigner functions of light and mechanics, respectively. Allowing for a free mechanical evolution for a time τ d before the interaction with the optical field, the output state is given byˆ . Using the result of Lemma 1 and the inverse of the symplectic transformation corresponding to the overall evolutionÛ oÛkÛom given by Eq. (21), we obtain 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 in Eq. (1) to getˆ out;o = . Using the orthogonality relation of Weyl-Wigner operators in Eq. (3), the final measured Wigner-function of the optical field can be readily ob- . Switching to position and momentum coordinates (x , p ) of the optical output, we obtain showing that only the momentum quadrature of the optical output is modified by the interaction.

C. General case readout
In principle, an optical vacuum state would be sufficient as an optical probe pulse for our method. However, in the general case we now consider a momentum-squeezed vacuum state with Wigner function W o (x, p; 0) = π −1 exp{−(e −2ε x 2 + e 2ε p 2 )} where ε ∈ [0, ∞) 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. (22) and changing the variable p → 2χup , we obtain the Wigner function of the optical output as From Eq. (23) it is evident that the information about the mechanical (ϕ + ϕ d + π/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 In which w out;o (p, π 2 ) is the momentum tomogram of the optical output pulse. Interestingly, we can identify the Gaussian kernel in Eq. (24) with the s-parameterized quasiprobability of the tomographic operator in Eq. (12), obtaining with order parameter As we show in Appendix D, the order mismatch can be moved to the mechanical state, which leads to Here, we have defined the s-parameterized mechanical tomogram w m (x, ϕ; s) to be the marginal of the s-parameterized quasiprobability distribution of the mechanical quantum state. Importantly, in view of Eq. (9), the relation obtained in (27) represents a Radon transformation of the (−s )parameterized quasiprobability distribution W m (β; −s ).
Substituting Eq. (27) into Eq. (24) we find where s is given by Eq. (26). Equation (29) is the main result of this manuscript, implying that the (−s )-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 larger sensitivity in the optical homodyne measurement to access the optical tomograms over the required range of values. This can be avoided by setting 2χu = 1 and instead exploiting the momentum squeezing of the initial optical state to set the order parameter s . We will discuss in detail the experimental feasibility of all such assumptions in Sec. IV.

D. Strong interaction regime
Recalling Eq. (26), it turns out that the order parameter s goes to zero quadratically with the displacement amplitude r and exponentially with the squeezing parameter ε. Hence, for sufficiently large χue ε we can approximate which in turn, implies that the φ-tomogram of the mechanical Wigner function is imprinted in a rescaled form on the momentum quadrature of the optical output Wigner function,

III. NONCLASSICALITY VERIFICATION
The method introduced above allows for the full reconstruction of the 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, this 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 in Eq. (4). In this case the bosonic coherent state basis is given byT (α; −1) = |α α| and thus, 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, Pfunction negativity can, in principle, be directly observed using appropriate filtering of the homodyne data [21][22][23], or witnessed without full tomography, using only a small number of specific observables [25][26][27][28][29].
FIG. 4. Witnessing P-function non-classicality. Using the method outlined in the main text (c.f. Fig. 3) one can obtain a s-tomogram of the state of the mechanical resonator. Using this s-tomogram, one computes the fictitious s-parameterized quasiprobability distribution W f (x, p; s) = wm(x; s)w0(p; s), where w0(p; s) is an s-tomogram of the vacuum state. If the resulting function 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 [36], a single mechanical tomogram can be used to witness nonclassicality as follows. Given a Wigner tomogram w m (x) of a quantum state for some quadrature angle φ (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 de- } is a marginal of the vacuum state. If the fictitious Wigner function W f (x, p) fails to represent a legitimate quantum state, then the state of the mechanical resonator has been P-nonclassical.
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 → p in w m (x), to construct a ficticious Wigner functions as (w m (x), w m (x)) → W f (x, p) = w m (x)w m (p). Similarly, the failure of W f (x, p) to represent a legitimate quantum state implies the Pnonclassicality of the state of the mechanical resonator. Importantly, testing whether the result is a legitimate Wigner function can be relatively simple in many cases of interest [36]. 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 m (x; s), one can define the following two demarginalization maps, in which w 0 (x; s) is the s-parameterized tomogram of the vacuum state and W f (x, p; s) is the s-parameterized quasiprobability distribution of a fictitious quantum state. The failure of either of the distributions in Eq. (33) or (34) 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. [36] by using appropriate quasiprobability distributions to compute the matrix elements of the fictitious density operator using the trace relation of Eq. (5). This shows that our method can be used to verify P-function nonclassicality even in the regime of weak interactions.

IV. EXPERIMENTAL PROPOSAL
We now describe a practical protocol for performing tomography of a mechanical resonator using our technique. First, recall that the phase ϕ is set by the pulse duration τ via the condition of Eq. (20) and substituting into Eq. (17). For a given pulse duration and depending on the optomechanical system parameters, the best pulse amplitude is then |α| = r such that the energy-time uncertainty relation is satisfied. Although using these parameters fixes the readout quadrature to ϕ = ϕ 0 , Eq. (31) shows that arbitrary quadrature angles ϕ d can be probed by exploiting the free evolution of the mechanics in the time delay τ d between state preparation and readout, which can be precisely controlled in practice, see Fig 3 (a). 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. (31). 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. ε = 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 τ ∼ π/ω mec , with a relative stability on the order of 10 −4 , which, using Eqs. (17) and (20), implies u ∼ 2 and χ ∼ k/2 √ 3 for k ∈ 2Z + . 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 χ ∼ 3 to satisfy the condition of Eq. (30), up to the order of 10 −3 , which is below other typical noise sources in the experiment. Note that, Eq. (20) also requires a relative stability of the laser power on the order of 10 −3 .
The relevant parameters for a number of optomechanical systems, taken mainly from Ref. [1] are summarized in Tab. I. We note that although many of these systems are indeed in the resolved sideband regime, they could nonethe-less employ our technique by sufficiently decreasing the quality factor of the used optical cavities, as outlined in Tab. I. Let us now consider explicitly kHz-frequency resonators such as Si 3 N 4 membranes [37,38], which are promising candidates for ground-state cooling [39] 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. (18) and (20), the requirements on the optical power to reach χ ∼ 3 then depends on the optomechanical coupling strength g 0 for the respective system and the readout displacement amplitude. Note, that the power requirements of our method can be decreased significantly by using a low finesse optical cavity, such that the optical decay rate κ o is much larger than the mechanical frequency ω m . As an example, the system of Ref. [37] in Tab. I has a suitable cavity with a sideband resolution factor of κ o /ω mec ∼ 50, and is thus in the ideal regime for our method. Making use of the optical cavity with a finesse of ∼ 3 × 10 4 , the required optical power is on the order of 10 −13 W. Optical displacements [40] and homodyne tomography of the optical output field can then be performed using standard methods [41,42].
Finally, we would like to point out that the use of squeezed readout light, though not necessary, would be beneficial for our protocol through reducing the requirements on amplitude and precision of the optical displacement operations. Specifically, one could choose χ ∼ 1 which, in turn, requires a squeezing parameter ε ∼ ln(3) to reach the equivalent approximation as before. This is equivalent to a vacuum squeezing of ∼ 9.5dB, which is achievable with current technology. In this case, both the requirements of Eqs. (20) and (30) are satisfied easily without high precision in the control of the laser power. Importantly, the squeezing parameter does not enter into the interaction transformation of Eq. (21), implying that it need not even be precisely characterized, as long as Eq. (30) is satisfied.

V. CONCLUSION
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 resourceefficient way without requiring full quantum state reconstruc-  [48] 3.9 × 10 9 3.1 × 10 −16 3.9 × 10 3 9 × 10 5 5 × 10 8 0.13 ∼ 5.1 × 10 −10 ∼ 5.3 × 10 −11 Verhagen et al [49] 7.8 × 10 7 1.9 × 10 −12 3.4 × 10 3 3.4 × 10 3 7.1 × 10 6 0.09 ∼ 2.6 × 10 −8 ∼ 1.5 × 10 −9 Teufel et al [50] 1.1 × 10 7 4.8 × 10 −14 32 2 × 10 2 2 × 10 5 0.02 ∼ 1.8 × 10 −7 ∼ 8.5 × 10 −9 TABLE I. Experimental parameters for different optomechanical systems. For each system we provide the optimal pulse length τ for the probe pulse used by our scheme, together with the required pulse energy for Wigner-function tomography based on the parameter assumptions ε = 0, u = 2, and χ ∼ 3.
tion. 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.

ACKNOWLEDGMENTS
We thank Michael Vanner and Andrew White for helpful discussions. This work was supported in part by the Centre for Quantum Computation and Communication Technology (CE110001027) and the Engineering and Physical Sciences Research Council (grant number EP/N002962/1).

Appendix A: Disentangling the interaction
The joint state after all interactions can be written aŝ in whichD(α) = exp{αa † − α * a} is the optical displacement operator,τ is light pulse duration, and τ 0 is the time between the optomechanical interaction and the optical displacement. The r.h.s of Eq. (A1) follows from the fact that a † 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. [51], namely where , This giveŝ in whichD(α) = exp{αa † − α * a} is the optical displacement operator. In Eq. (A4), whenever needed we have used the fact that the optical and mechanical mode operators commute, and that a † 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 ω 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 relationŝ the second bracket is understood to be the free rotation of the optical field around a phase space point located at −α. 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. (A5), setting α = re iθ , and noting the fact that |α| = r 1, we can approximatê 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. (A4) as and linearise the exponent as where µ = 1 − e iωmτ = ue iϕ so that Curiously, the radiation pressure interaction term g0 ωm a † a(µb † + µ * b), has become negligible in Eq. (A8), since g0 ωm g0ru ωm and can thus be omitted. Note that this also requires that u be chosen large enough, which implies that the optimal pulse length τ is on the order of half the mechanical period. We can now define the new optical mode annihilation operatorã θ := ae −iθ + r 2 and the new mechanical annihilation operatorb ϕ = be −iϕ to get Note that, since 0 u 2 and we work in a parameter regime where χu = g0ru ωm ∼ 1 − 10, it follows that χ 2 χ and hence the optical Kerr effect in Eq. (A6) is not negligible compared to the interaction term of Eq. (A10).

(A11)
Appendix B: Transformation of the bosonic operators

Optomechanical interaction
As outlined in the main text, we can now use Eq. (A11) to compute the transformation of the vector of bosonic operatorŝ A := (a † , b † , a, b) under the optomechanical interaction given bŷ Exploiting the properties of symplectic transformations [34], we know that ifÛ = exp( 1 thenÛˆ AÛ † =ˆ AM . In terms of the new the new field operatorsã θ andb ϕ , which gives and thus, with det M = 1. Using the matrix of Eq. (B4), we find is the (ϕ + π 2 )-quadrature of the mechanics, andx(θ) = 1 √ 2 (a † e iθ + ae −iθ ) is the θ-quadrature of the optical field. Making use of Eq. (B1), and thus, by virtue of Eq. (B5), we find the relations of Eq. (21) in the main text, (B8)

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ˆ A := (a † , b † , a, b), we evaluate the output mode operatorŝ To find the result, we first define the new mode operators [35] a = ae −iθ + r 3 .
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ā, we haveÛ Using the same relation as in Sec. B 1, where v = χ 2 [ω m τ − sin(ω m τ )], and thus, in whichv = 2 √ 3v, and we have . Using the matrix of Eq. (B11), we find Equating this with Eq. (B12), we find The linear transformation corresponding to Eq. (B14) is thus given bŷ where the displacement vector D is defined as

Optical transformation
We now consider the effect of the last unitary evolution in Eq. (A11), namely,Û o =D(−α)e −i(τ0+τ )ωoa † aD (α). 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. (A5), while the action of the phase-rotation operator is well-known to be It is then straightforward to evaluate the following: This represents a rotation of the optical field in phase space by an angle of (τ 0 + τ )ω o around the phase-space point α.

The total evolution
To obtain the effect of the overall evolution operator as described in Eq. (A11), we simply combine the results of Eqs (B1), (B14), and (B18). First, combining the Kerr effect with the optomechanical interaction gives Next, after the optical displacement, the final transformation is obtained as (B20) In Eq. (B20), 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, θ, 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 θ = 0, that is α = r ∈ R. Consequently, Another simplification can be obtained by noting that one can arrange the delay for the displacement of the output optical field (τ 0 ) with high accuracy in such a way that (τ 0 + τ )ω o = m2π for m ∈ N. This gives Finally, if we choose the probe pulse duration τ in such a way that sin(v) = 0, we will obtain the simplest transformation equations. For this, we recall thatv = 2 In general, it is always possible to find a pair of parameters (τ, r), such that 0 ω m τ 2π, and Eq. (B23) is satisfied. However, as indicated earlier, the optimal pulse length to achieve u = 2 is ω m τ = π. In this case, the condition in Eq. (B23) can be satisfied by adjust the amplitude r of the optical coherent state, and thus the parameter χ. Chosing τ such thatv = kπ, we obtain where u and ϕ are given by Eq. (A9). We see that, for even k, the displacement term completely goes away, and we arrive at Eq. (21), that is, To summarize our parameter choices, consider the case where we choose the optimal pulse length τ such that ω m τ = π. This implies u = 2 and by Eq. (B23), for example for k = 32, χ = 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 χ ∼ 1 − 10, one can find a pulse length that satisfies the condition in Eq. (B23) with u = 0, for the purpose of Wigner function tomography. Importantly, for any pair of parameters (τ, r) satisfying Eq. (B23) such that 0 ω m τ 2π, a s-parametrized quasiprobability distribution of the mechanics is possible, as described in Sec. II C in the main text.
From Eq. (B25), after disregarding the mechanical constant displacement, we find the following symplectic matrices, det S = 1.
To check the legitimacy of the fictitious density operators resulting from these demarginalization maps, one may follow the standard arguments provided in Ref. [36], such as the Kastler-Loupias-Miracle-Sole test [52][53][54]. 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. (5) and the corresponding s-parametrized Fock basis projections given by [55] W |n m| (α; s) = ( 2 1 − s ) m+n+1 1 √ n!m! h n,m (α * , α| s 2 − 1 4 )e − 2|α| 2 1−s , in which h n,m (x, y| ) are the incomplete 2D Hermite polynomials defined as [56] h n,m (x, y| ) = min{m,n}