# Room-temperature Mechanical Resonator with a Single Added or Subtracted Phonon

Rishi N. Patel, Timothy P. McKenna, Zhaoyou Wang, Jeremy D. Witmer, Wentao Jiang, Raphaël Van Laer, Christopher J. Sarabalis, Amir H. Safavi-Naeini

RRoom-temperature Mechanical Resonator with a Single Added or Subtracted Phonon

Rishi N. Patel, ∗ Timothy P. McKenna, Zhaoyou Wang, Jeremy D. Witmer, WentaoJiang, Rapha¨el Van Laer, Christopher J. Sarabalis, and Amir H. Safavi-Naeini † Ginzton Laboratory, Stanford University Department of Applied Physics

A room-temperature mechanical oscillator undergoes thermal Brownian motion with an amplitudemuch larger than the amplitude associated with a single phonon of excitation. This motion can beread out and manipulated using laser light using a cavity-optomechanical approach. By performinga strong quantum measurement, i.e. , counting single photons in the sidebands imparted on a laser,we herald the addition and subtraction of single phonons on the 300 K thermal motional stateof a 4 GHz mechanical oscillator. To understand the resulting mechanical state, we implement atomography scheme and observe highly non-Gaussian phase-space distributions. Using a maximumlikelihood method, we infer the density matrix of the oscillator and conﬁrm the counter-intuitivedoubling of the mean phonon number resulting from phonon addition and subtraction.

A mechanical oscillator at room temperature will be-have in nearly perfect accordance with the laws of clas-sical statistical physics. Nonetheless, interaction withoptical-frequency photons can lead to behaviour that isnonclassical. For the mechanical motion of trapped ions,laser cooling and coupling of motion to the internal elec-tronic degrees of freedom has long been pursued as a pathto realizing a scalable quantum computer [1]. For solid-state mechanical devices, the signatures of quantum noiseand back-action have been observed at room tempera-ture and proposed as a means to realize quantum sen-sors [2, 3]. A key feature of such optomechanical devicesis their ability to eﬃciently generate correlations betweenmotion and light [4]. A strong quantum measurement ofthe resulting light ﬁeld, e.g, by using a single-photon de-tector can consequently alter the state of the mechanicalsystem [5–10].Experiments in nonlinear optics have shown thatadding or subtracting photons fundamentally alters thestate of an optical degree of freedom. For example, inphoton addition, parametric down conversion followedby post-selection of an idler photon allows the prepara-tion of a single-photon state from vacuum [11]. Photonsubtraction, when operating on squeezed light is used togenerate and enlarge Schr¨odinger’s cat states [12, 13].Both photon addition and subtraction are essential for anumber of tasks in continuous variable quantum informa-tion processing, in which non-Gaussian states are oftenrequired [14]. Highly thermal states of light, where sig-niﬁcant classical noise would be expected to wash awayany quantum eﬀects, have also been converted to non-classical states by the addition of one photon [15–18].Such states, combined with photon subtraction, havealso allowed a direct test of quantum commutation re-lations [19, 20].Mechanical oscillators have emerged as an importantavenue for realizing quantum technologies. Cryogeni-cally cooled oscillators are used more widely due to sig- ∗ [email protected] † [email protected] niﬁcantly reduced dissipation and reduction of thermalnoise from the environment, which usually masks quan-tum features. Motivated by the long intrinsic relax-ation times possible in cryogenically cooled mechanicaloscillators [21], and success in strongly coupling themto superconducting qubits [22–25], proposals for realiz-ing quantum machines that leverage their coherence andsmall size [26–28] have emerged recently. Separately,cavity-optomechanical approaches for quantum sensingand transduction are being pursued by groups aroundthe world. Nonetheless, there is signiﬁcant interest inoperating quantum systems at higher temperatures sincefor many applications operation at ambient conditions isessential.In this work, we perform experiments on a mechan-ical oscillator in a highly thermal state at room tem-perature, and use its interaction with optical photonsto perform quantum-limited read-out and state control.First, we use single photon counting to prepare phonon-added and subtracted mechanical states in a regimewhere kT (cid:29) (cid:126) ω m . Second, building on tomography tech-niques in quantum optics [11, 14, 29, 30] and optome-chanics [31, 32], we combine heralded single-phonon addi-tion with the continuous measurement of mechanical am-plitude and phase ﬂuctuations. We reconstruct the den-sity matrix of a mechanical oscillator in a phonon-addedor phonon-subtracted thermal state with initial meanphonon occupancy of approximately kT / (cid:126) ω m = 1580.With future improvements, our technique may be ex-tended to prepare and characterize more complex statesof mechanical motion.The process of phonon-addition and subtraction arisesfrom the inelastic scattering of light from a laser dueto mechanical motion in a cavity. By energy conserva-tion, photons scattered by mechanical motion are shiftedto a lower (higher) frequency corresponding to addition(subtraction) of a phonon in the mechanical resonator(Fig. 1a). Selection of the desired process (addition orsubtraction) is achieved by using the optical cavity res-onance as a ﬁlter and tuning the laser to its blue or redside. This can be understood by considering the Hamilto-nian that describes the optical and mechanical systems. a r X i v : . [ qu a n t - ph ] F e b For blue-detuning it is given by ˆ H = (cid:126) G (ˆ a † ˆ b † + h.c. )where G is the linearized optomechanical coupling rate,and the annihilation operators for optical and mechanicaloscillators are given by ˆ a and ˆ b respectively. Similarly,for red-detuning we have ˆ H = (cid:126) G (ˆ a † ˆ b + h.c. ) whereby aphonon is annihilated while a photon is generated. In ourexperiments, the optical cavity has a decay rate κ thatis much faster than G , and its photons are sent to a de-tector. Instead of coherent dynamics between the opticalﬁeld and motion, detection leads to a quantum opera-tion with corresponding jump operators proportional toˆ b † and ˆ b for the two Hamiltonians, respectively. There-fore, the detection of an individual photon at the cavityresonance frequency heralds the addition or subtractionof a phonon to the mechanical mode characterized bythe operation of the respective jump operator. The stateof the mechanical oscillator will be thermal before therecord of detections is taken into account. This reﬂectsour lack of knowledge about the motional state and thatthe mechanical oscillator is in equilibrium with its envi-ronment. Starting with the thermal density matrix ˆ ρ th ,detection of a photon corresponds to updating the statewith the jump operator corresponding to the correct de-tuning. For the blue-detuned case, we use the jump op-erator proporational to ˆ b † to obtain the phonon-addedstate: ˆ ρ th → ˆ b † ˆ ρ th ˆ b . In the red-detuned case, we havephonon-subtraction represented by: ˆ ρ th → ˆ b ˆ ρ th ˆ b † .Our experimental setup is shown in simpliﬁed form inFig.1b. First, we send light from a laser into an op-tomechanical crystal cavity, either red-detuned or blue-detuned from the cavity resonance by the mechanical fre-quency. The ﬁber to chip coupling eﬃciency is η f ≈

76 %.We split oﬀ some of the light before interacting withthe device, to use as a local oscillator. A delay is ap-plied using about 100 meters of ﬁber to approximatelycompensate for the signal path. A continuous wave, fre-quency upshifted (40 MHz), probe tone is generated usingan acousto-optic-modulator. The reﬂected light from thecavity is then split into two paths. One path, for sin-gle photon counting, contains two cascaded high ﬁnesseﬁber fabry-perot cavity ﬁlters with free spectral rangeof 15 GHz and ﬁnesse of 300 (Micron Optics FFP-I) tosuppress the pump light and pass through only the pho-tons due to scattering from the mechanical resonator.Another path, for heterodyne detection, contains a bal-anced heterodyne receiver (Thorlabs 75 MHz PDB425C-AC) whose output is sent to a 12 bit 500 MS s − digitizer(Alazartech ATS 9350). The local oscillator and signalare combined on a variable optical coupler before the bal-anced detector. The DC output level of the two photo-diodes that comprise the receiver is monitored, and thecoupler splitting ratio is adjusted until the voltages areapproximately equal. We note that the phase of our lo-cal oscillator is left unlocked. This simpliﬁes the exper-iment but means that we extract no information aboutthe phase of the system, causing our inferred states tohave rotationally symmetric quasi-probability distribu-tions. Our justiﬁcation is that the initial state is a ro- tationally symmetric thermal state, and that the pho-ton addition and subtraction processes occur at randomtimes and have no associated phase.In the photon counting path, we send the light reﬂectedfrom the room temperature device to a superconductingnanowire single photon detector (SNSPD) which resideson the still plate of a dilution refrigerator (Bluefors) at ≈ µ A.Our optomechanical crystal device, used to read outand control mechanical motion, is similar to those pre-sented in previous work [33]. The optical cavity modehas a center wavelength of λ c = 1551 nm, total decayrate κ/ π = 822 MHz, and external coupling rate of κ e / π = 190 MHz. The probability that a cavity pho-ton leaks into the detected waveguide channel is given bythe cavity eﬃciency, η c = κ e /κ ≈

23 %. The mechani-cal mode frequency is ω m / π ≈ .

96 GHz. From a sweepof laser power, we determine the intrinsic, backaction-free, mechanical linewidth of γ i / π = (2 . ± .

01) MHzand single-photon optomechanical coupling rate g / π =(1 . ± .

02) MHz. The measured value of g deviates byabout 5 % from the simulated g / π = 961 kHz, a diﬀer-ence which we can attribute to systematic errors in powercalibration as well as uncertainties in the material’s pho-toelastic parameters. In terms of the system parametersabove, the single-photon generation rate per phonon inthe mechanical resonator, is proportional to the optome-chanical measurement rate: γ OM = 4 g n cav /κ , where n cav is the number of optical intracavity photons (on or-der 10 in this experiment).The ﬁrst experiment we perform is phonon-addition.We tune the laser on the blue side of the optical cavityresonance. The experimental data is collected by trig-gering the digitizer on a single photon click, collecting ≈ µ s of heterodyne data, and estimating the in-phaseand in-quadrature components of the down-convertedmechanical signal. The resulting complex voltage sam-ples, v = √ G ( X + iP ) comprise the dataset with whichwe perform tomography. Here, G is the overall detectiongain and X and P are the in-phase and in-quadraturecomponents respectively. We collect the quadrature sam-ples in two interleaved phases, one in which data collec-tion is triggered by a single photon click, the other inwhich the clicks are ignored. When the measurement istriggered by a single photon, a phonon-added thermalstate is heralded. The histogram of quadrature samples,obtained by binning the raw data for a mechanical ther-mal state, is shown in Fig. 2a. These data were binnedinto 101 bins in both the X and P directions. The tomog-raphy of the thermal state shows a Gaussian distributionof quadrature amplitudes (Fig. 2a). By contrast, we ob-serve a clear non-Gaussian rotationally symmetric distri-bution in the phonon-added thermal state in Fig. 2b.The Husimi Q function for the post-selected phonon- FPC

AOM

300 K λ ~ FP OMC

Single-Photon-Detector(SNSPD) -

20 MHz

40 MHz

Digitizer BHDMatchedFilterIQ Samples Trigger IQ (a)(b) n m Phonon Addition Phonon Subtraction

FIG. 1.

Concept and principal components of the ex-perimental setup (a) Concept showing on left (right) how aRaman scattered photon on cavity-resonance at frequency ω c heralds the addition (subtraction) of a phonon at the mechan-ical frequency ω m . Insets from left to right show simulatedmechanical and optical mode proﬁles respectively. (b) Dia-gram showing the main components of the setup. A tunablelaser is stabilized to a Fabry-Perot cavity (FP) and is fre-quency upshifted using an acousto-optic modulator (AOM).The resulting signal polarization is set using a ﬁber polar-ization controller (FPC) and sent into the optomechanicalcrystal device (OMC) whose temperature is stabilized near300 K. The reﬂected light is incident on a beamsplitter andis split into two paths. One (going to the bottom) containspump-rejection ﬁltering with a pair of FP ﬁlters, followed bysingle photon detection using a superconducting nanowire sin-gle photon detector (SNSPD). The other path (going to theleft) performs balanced heterodyne detection (BHD) wherethe resulting signal is further digitally downconverted from20 MHz, and ﬁltered using a digitizer and GPU. The digitizercaptures ≈ µ s of data every time it is triggered by a singlephoton pulse from the SNSPD. A schematic of the decayingcavity ﬁeld is shown inset in the ﬁgure. τ denotes the time-lagbetween when the click occurs and when the heterodyne datacollection starts, as described in the text. added thermal state has the form: Q post ( α ) ∝ (cid:104) α | ˆ b † ˆ ρ th ˆ b | α (cid:105) ∝ | α | e −| α | / (¯ n th +1) , (1)where ¯ n th is the mean thermal phonon occupancy and α = X + iP . Q post ( α ) describes the measurement statis-tics in phase space, in the absence of technical noise.Since the Q function is a probability distribution, the ef-fect of added Gaussian noise can be represented by itsconvolution with a Gaussian distribution N , with zero mean and a variance of n added : Q ( α ) = (cid:16) Q post ∗ N (cid:17) ( α ) . (2)(See Appendix for an analytic expression for the mea-sured Q ( α ) which depends only on ¯ n th and n added ). Torelate to experimental observation, an additional scal-ing parameter G is needed. This distribution is binned,and ﬁt to the data via a maximum-likelihood method.Since the bath temperature is known ( T ≈

300 K), weﬁx ¯ n th = 1578 in Eq. 1 and ﬁt Eq. 2 via two freeparameters: detection gain G , and number of addednoise phonons n added . The parameter G is used to re-scale the data as v = √ Gα . The result of the ﬁt with( G, n added ) = (8 . × − V , n th = 1597 . ± . n post = 3158 . ± . p vac,th = (6 . ± . × − to p vac,post = (5 . ± . × − . In all quantities es-timated from the state reconstruction, the quoted errorsreﬂect statistical uncertainty obtained using a bootstrap-ping method of the entire dataset. This method works byre-sampling the entire dataset of 1 . × samples withreplacement 50 times, and reconstructing the density ma-trix for each trial. Lastly, we perform the reconstructionof the phonon-subtracted state, shown in (Fig. 3c). -100-50050100 P X -100-50050100 P -100 -50 0 50 100 X(a) (b)(c) (d)Experiment ExperimentTheory Theory -100 -50 0 50 10002468101214 C oun t s (e) X10 Pre-detection Post-detection

FIG. 2.

Measurement of quadrature histograms forthermal and phonon-added thermal state (a) Experi-mental result showing a histogram of 1 . × quadraturesamples. The data are binned into 101 bins along the X and P axes. Color bars show counts per bin. The statistics areGaussian, corresponding to a mechanical thermal state. (b)Post-selected results, triggered by single-photon clicks. Theobserved statistical distribution now corresponds to a phonon-added state and is non-Gaussian. (c,d) Corresponding theoryplots to (a,b), obtained by ﬁtting the analytically determinedHusimi Q functions via two parameters. The plots shownare the binned theory Q functions for both the thermal andpostselected case. The dotted black circles, shown for ref-erence in the plots, have radii r = √ ¯ n th , the characteristiclength scale for thermal ﬂuctuations. ¯ n th ≈ .

96 GHz oscilla-tor at room temperature. (e) Blue (red) bars correspond toa linecut of the 2D experiment histograms at P = 0 for thethermal (post-selected) datasets. The theory ﬁts for the ther-mal (post-selected) results are shown in solid (dotted) blacklines. Next, we use our reconstruction results to investi-gate how the expected number of phonons changes afterpost-selection. We compute the ratio of the two meanphonon numbers as ¯ n post / ¯ n th = 1 . ± . or subtracting a phonon doubles the mean number ofphonons in a resonator, ¯ n post ≈ n th , is best understoodby considering the information gained about the mechan-ical system from the optical single photon measurement.Before the measurement of a photon occurs, the a-priori probability distribution over each phonon energy level, n , is given by the familiar exponentially decaying Boltz-mann factor: P ( n ) ∝ e − βn where β is the inverse tem-perature (Fig. 3b, blue curve). Once a click has occurredhowever, the observer gains information about the state,and we must update these probabilities via Bayes’ rule.Letting the number of resonant cavity photons in a smalltime interval be N , we have the following update rule: P ( n | N = 1) ∝ P ( N = 1 | n ) P ( n ), a rescaling of the a-priori distribution. Now, the probability of a photonscattering event itself depends on the phonon number: P ( N = 1 | n ) ∝ n . Thus we see that the a-posteriori prob-ability is the prior distribution re-scaled by n , leading tothe suppression of probability for small phonon numbers(Fig. 3b). This causes the average phonon number to beincreased.One can verify this intuitive argument by direct calcu-lation of the phonon-added state. Writing the thermalstate as a sum (for n ≥

0) using the prior, thermal,probability distribution: ˆ ρ th = (cid:80) n P ( n ) | n (cid:105)(cid:104) n | , we cal-culate the post-selected phonon-added state as ˆ ρ post =(ˆ b † ˆ ρ th ˆ b ) / Tr (cid:16) ˆ b ˆ b † ˆ ρ th (cid:17) . Recalling that ˆ b † | n (cid:105) = √ n + 1 | n (cid:105) and simplifying, we get: ˆ ρ post = (1 /n th ) (cid:80) n nP ( n ) | n (cid:105)(cid:104) n | .Notice that the priors have been updated: P ( n ) → nP ( n ) /n th . This analysis, while illustrated for phonon-addition, applies as well to phonon-subtraction in thelarge thermal occupation limit (see Eq. A4).As noted above we observe a near, though inexact,doubling in mean phonon number after post-selection.To understand this discrepancy with theory, we indepen-dently measure the dark count rate in our measurement.Dark counts introduce a loss in ﬁdelity of the heraldedstate. More precisely, the heralding ﬁdelity is deﬁned bythe quantity χ = Γ sig / (Γ sig + Γ dark ). Here, Γ sig denotesthe count rate of thermal signal phonons, while Γ dark de-notes the total dark count rate, caused by a sum of intrin-sic SNSPD dark counts and pump feedthrough from im-perfect pump rejection ﬁltering. We have independentlyestimated these quantities, using the measured photoncount rate on-resonance, and oﬀ-resonance. We measureΓ sig ≈ (278 ±

5) kHz and Γ dark ≈ (5 . ± .

1) kHz. Fromthis we estimate the ﬁdelity, χ ≈ . ± .

02. Written interms of ﬁdelity, the theoretically expected ratio of meanoccupancy is: ¯ n post / ¯ n th ≈ χ . The measured ﬁdelity χ thus closely explains the observed ratio ¯ n post / ¯ n th .Finally, we map out the time-evolution phonon-addedmechanical thermal state. We sweep the delay time τ ,which controls the temporal mode matching between thesignal at the heterodyne detector and the single photoncounter (Fig. 1b). As a function of delay, we compute (a) ThermalPhonon-addedn -4 nn Radius P r obab ili t y -2 (b)(c) Phonon-subtracted FIG. 3.

Reconstructed density matrix elements (a) Ra-dial histogram of measurement results for the thermal state(ﬁlled blue points), the phonon-added thermal state (open redcircles) and the phonon-subtracted thermal state (green tri-angles). The respective dotted lines show the theory ﬁts. (b)Diagonal density matrix elements reconstructed from experi-mental data. Estimates for gain and added noise are providedby the ﬁts to the quadrature histogram data. The black dot-ted lines show the result of theory for the thermal and phonon-added states with thermal bath occupancy set to ¯ n th = 1578phonons ( ≈

300 K mode temperature). (c) Experimentally re-constructed density matrix elements for a phonon-subtractedthermal state. the total variance of the quadrature histograms. An ex-ponential decay is observed, with a decay time on the or-der of the mechanical lifetime ( τ ≈

100 ns). We observethat the noise distribution near zero delay is distinctlynon-Gaussian, but tends towards a Gaussian thermal dis-tribution at large delay. The results, normalized to thevariance of the mechanical thermal noise Gaussian, areshown in Fig. 4.In this work we have experimentally demonstratedphonon addition and subtraction followed by state to-mography in an optomechanical system. These capabil-ities open up several directions for future studies. Firstof all, the single-phonon-added thermal states demon-strated here have theoretically been shown to be non-classical at all temperatures. This is ensured by the nega-tivity of their Glauber-Sudarshan P functions [16, 19, 36].Although this negativity is diﬃcult to detect experi-mentally at room temperature, we hope that our workwill motivate further studies into whether weakly non- -400 -300 -200 -100 0 100 200

Delay Time ( n s ) N o r m a li z ed V a r i an c e FIG. 4.

Time-evolution of phonon-added state

The topinsets from left to right show how the raw quadrature his-tograms change as the delay time τ is swept. The noise vari-ance, normalized to the variance of the thermal state, is plot-ted below as a function of τ . Red (blue) points denote theresults for post-selected (thermal) distributions respectively.The temporal mode-matching between the single photon pulseand the heterodyne signal is well deﬁned for τ ≥

0. As suchwe shade out the negative delay region, in which the matchedﬁlter of the detection chain only partially overlaps with thedecaying mechanical signal. classical, but non-Gaussian, states may prove useful asa quantum resource [37, 38]. Finally, we point out thatour experiment demonstrates the ability to add a nodeto the Q function of a single-mode oscillator. Since itis known that any pure state whose Q function con-tains nodes is non-classical [39], one could prepare non-classical states beyond the phase insensitive ones we havereconstructed. In particular, we expect with near-termimprovements, including the implementation of phase-sensitive detection, our technique will allow the prepara-tion of phonon-added coherent states of motion [40, 41].Likewise, phonon-subtraction, performed on a squeezedmechanical steady state, gives a route to cat state gener-ation in the mechanical domain [12, 42].During the preparation of this manuscript we becameaware of related work demonstrating variance doubling inphonon-added and subtracted mechanical thermal states[8].

ACKNOWLEDGMENTS

This work was funded by the U.S. governmentthrough the Department of Energy through Grant No.DE-SC0019174, and the U.S. Army Research Oﬃce(ARO)/Laboratory for Physical Sciences (LPS) Cross-Quantum Systems Science & Technology (CQTS) pro-gram (Grant No. W911NF-18-1-0103). The authors wishto thank NTT Research for their ﬁnancial and technicalsupport. We acknowledge the David and Lucille PackardFellowship, and the Stanford University Terman Fellow-ship. We thank Pieter-Jan C. Stas, Alex Wollack, HubertStokowski and Marek Pechal for experimental support.Device fabrication was performed at the Stanford NanoShared Facilities (SNSF) and the Stanford Nanofabri-cation Facility (SNF). The SNSF is supported by theNational Science Foundation under Grant No. ECCS-2026822. RNP was partly supported by the NationalScience Foundation Graduate Research Fellowship underGrant No. DGE-1656518.

APPENDIXAppendix A: Phonon-added and subtracted states

In this section we state results for the density matrixelements of phonon added and subtracted states. We alsogive analytic expressions of the Husimi Q function, fun-damental to analyzing the output of heterodyne tomog-raphy experiments, both for the noiseless (zero technicalnoise) and added-noise cases.

1. Density Matrix

In the formulas below, we use the shorthand notation β = 1 /kT and (cid:126) ω = 1. The mechanical thermal state isgiven by: ˆ ρ th = 1¯ n + 1 ∞ (cid:88) n =0 exp( − βn ) | n (cid:105)(cid:104) n | (A1)where Tr[ˆ n ˆ ρ th ] = ¯ n = 1 / (exp( β ) −

1) is the mean phononoccupancy. Following post-selection, the phonon-addedstate is:ˆ ρ a = ˆ b † ˆ ρ th ˆ b Tr (cid:104) ˆ b ˆ b † ˆ ρ th (cid:105) = 1¯ n (¯ n + 1) ∞ (cid:88) n =0 n exp( − βn ) | n (cid:105)(cid:104) n | . (A2)In this state, the mean phonon number is Tr[ˆ n ˆ ρ a ] = 2¯ n +1. The phonon-subtracted state is:ˆ ρ s = ˆ b ˆ ρ th ˆ b † Tr (cid:104) ˆ b † ˆ b ˆ ρ th (cid:105) = 1(¯ n + 1) ∞ (cid:88) n =0 ( n + 1) exp( − βn ) | n (cid:105)(cid:104) n | (A3)or: ˆ ρ s = 1¯ n + 1 (ˆ ρ th + ¯ n ˆ ρ a ) (A4)Note that for ¯ n (cid:29)

1, as in this experiment, ˆ ρ s ≈ ˆ ρ a

2. Husimi Q functions for phonon-added andsubtracted states

The Husimi Q function is deﬁned as: Q ( α ) = 1 π (cid:104) α | ˆ ρ | α (cid:105) . (A5)Substituting the three states of interest into this equationgive the following results.For a thermal state with mean phonon number ¯ n theQ function is given by: Q th ( α ) = 1 π (¯ n + 1) exp (cid:18) −| α | ¯ n + 1 (cid:19) . (A6)The phonon-added state Q function is: Q a ( α ) = 1 π (¯ n + 1) | α | exp (cid:18) −| α | ¯ n + 1 (cid:19) . (A7)Note in particular that this distribution is no longerGaussian, and goes to 0 at α = 0 for all temperature(all values of ¯ n ).The phonon subtracted Q function can be expressedas a weighted sum of the previous two functions. Statedexplicitly: Q s ( α ) = 1 π (¯ n + 1) (cid:0) n | α | (cid:1) exp (cid:18) −| α | ¯ n + 1 (cid:19) . (A8)This equation implies that the non-Gaussian characterof the phonon-subtracted state increases with increasingtemperature.

3. Husimi Q functions in presence of technicalnoise

The detected Q function for the post-selected state, inthe presence of technical noise is: Q ( r ) = 14 π (cid:18) σ σ (cid:19) (cid:18) σ + ˜ σ σ r (cid:19) exp (cid:18) − r σ + σ ) (cid:19) , (A9)where σ = (¯ n + 1) / σ ≈ n added / σ = σ σ σ + σ (A10) Appendix B: Derivation of Time-Independent Fieldsafter Filtering

In this section, we present an input-output theory anal-ysis of the detected noise, and show the dependence ofdetected noise on the matched ﬁlter bandwidth and delaytime. We show the exponential decay of the signal vari-ance with time, which we have detected in experiment.

1. Output ﬁelds from input-output theory

The output cavity ﬁeld is sampled continuously in thisexperiment. The output of the detector is a time inde-pendent quantity, obtained by integrating the ﬁeld witha ﬁlter function. Our treatment of this problem is similarto [30], although we work in a regime in which thermalnoise from the mechanical mode plays a signiﬁcant role.In the following we derive expressions for the integratedoutput cavity ﬁeld:ˆ A out = (cid:90) f ( t )ˆ a out ( t ) dt (B1)And the resulting measured photocurrent:ˆ I = α LO ( ˆ A out + ˆ A † out ) (B2)where α LO denotes the local oscillator strength (whosephase is random in our experiment, but which we taketo be real for simplicity), and the ﬁlter function is givenby: f ( t ) = Θ( t − τ ) (cid:112) β exp (cid:18) − β t − τ ) (cid:19) . (B3)In the above equation β denotes the matched ﬁlter en-ergy decay rate, and the function Θ( t − τ ) is the Heavisidestep function starting at time τ > (cid:90) ∞−∞ (cid:107) f ( t ) (cid:107) dt = 1 . (B4)To proceed, we write the Heisenberg-Langevin equa-tions describing the dynamics of the system operators ina frame rotating at the mechanical frequency: ddt ˆ a ( t ) = − κ a ( t ) − iG ˆ b ( t ) −√ κ e ˆ a in ( t ) −√ κ i ˆ a in,i ( t ) , (B5) ddt ˆ b ( t ) = − γ b ( t ) − √ γ i ˆ b in ( t ) + ˆ F BA ( t ) (B6)along with the input-output boundary condition:ˆ a out ( t ) = ˆ a in ( t ) + √ κ e ˆ a ( t ) . (B7)Note that in Eq. B5 the ﬁnal term arises from vacuumnoise due to undetected channels [43]. In Eq. B6 we havetaken into account the eﬀect of the optical mode into themechanical damping rate by writing γ = γ i ± γ OM , where γ OM = G κ is the optomechanical measurement rate, andthe minus sign is chosen for blue laser-cavity detuning.Assuming κ (cid:29) G, γ, β throughout, we set the left handside of Eq. B5 to 0 and substitute ˆ b ( t ) with the formalsolution of Eq. B6:ˆ b ( t ) = ˆ b (0) exp (cid:16) − γ t (cid:17) −√ γ i (cid:90) t exp (cid:16) − γ t − t (cid:48) ) (cid:17) ˆ b in ( t (cid:48) ) dt (cid:48) (B8)to obtain an equation for ˆ a ( t ) and ˆ a out ( t ):ˆ a ( t ) = − iGκ (cid:20) ˆ b (0) exp (cid:16) − γ t (cid:17) − √ γ i (cid:90) t exp (cid:16) − γ t − t (cid:48) ) (cid:17) ˆ b in ( t (cid:48) ) dt (cid:48) (cid:21) − √ κ e κ ˆ a in ( t ) − √ κ i κ ˆ a in,i ( t ) (B9)ˆ a out ( t ) = − i √ η c γ OM (cid:20) ˆ b (0) exp (cid:16) − γ t (cid:17) − √ γ i (cid:90) t exp (cid:16) − γ t − t (cid:48) ) (cid:17) ˆ b in ( t (cid:48) ) dt (cid:48) (cid:21) + (1 − η c )ˆ a in ( t ) − √ κ i κ e κ ˆ a in,i ( t ) (B10)where in Eq. B10 we substituted the deﬁnition of cav-ity eﬃciency η c = κ e /κ and γ OM . We have ignored herethe contribution due to the optical read-out ﬁeld on themechanical motion ˆ F BA ( t ), as this back-action noise ismuch smaller than the noise we will be measuring [43, 44].The ﬁrst two terms of this equation represent contribu- tions from the decaying mechanical mode and the ther-mal noise respectively. The last two terms represent con-tributions from optical vacuum noise arising from bothdriven and undetected channels.To obtain the time independent ﬁeld via Eq. B1 weevaluate:ˆ A out = − i (cid:112) βη c γ OM exp( βτ / (cid:34) ˆ b (0) (cid:90) ∞ Θ( t − τ ) exp (cid:18) − γ + β t (cid:19) dt −√ γ i (cid:90) ∞ (cid:90) t Θ( t − τ ) exp (cid:18) − γ + β t (cid:19) exp (cid:16) γ t (cid:48) (cid:17) ˆ b in ( t (cid:48) ) dt (cid:48) dt (cid:35) +(1 − η c ) exp( βτ / (cid:90) ∞ Θ( t − τ ) exp (cid:18) − β t (cid:19) ˆ a in ( t ) dt − √ κ i κ e κ exp( βτ / (cid:90) ∞ Θ( t − τ ) exp (cid:18) − β t (cid:19) ˆ a in,i ( t ) dt (B11)which reduces to:ˆ A out = − i √ βη c γ OM γ + β (cid:34) ˆ b (0) exp (cid:16) − γ τ (cid:17) − √ γ i exp( βτ / (cid:90) ∞ τ exp (cid:18) − β t (cid:19) ˆ b in ( t ) dt − √ γ i exp( − γτ / (cid:90) τ exp (cid:16) γ t (cid:17) ˆ b in ( t ) dt (cid:35) + (cid:112) β (1 − η c ) exp( βτ / (cid:90) ∞ τ exp (cid:18) − β t (cid:19) ˆ a in ( t ) dt − √ βκ i κ e κ exp( βτ / (cid:90) ∞ τ exp (cid:18) − β t (cid:19) ˆ a in,i ( t ) dt. (B12)This equation allows us to compute the ﬂuctuations ofthe photocurrent as a function of matched ﬁlter linewidthand start time: (cid:104) ˆ I † ˆ I (cid:105) ( τ, β ) = α (cid:104) ( ˆ A out + ˆ A † out ) (cid:105) . (B13)Neglecting cross-correlations that appear in Eq. B13, weevaluate all expectation values in angular brackets givena state ρ = ˆ ρ a (cid:78) ˆ ρ thermal , where ˆ ρ a is the post-selectedstate after heralding, given by Eq. A2, and where all in-put operators act on the ˆ ρ thermal state. In our derivationwe note that the expectation value of the system opera-tors in the post-selected state is: (cid:104) ˆ b † (0)ˆ b (0) (cid:105) ≈ n γ i γ .(Itdiﬀers from 2¯ n due to backaction from the laser drive).The ﬁnal result for τ ≥ n (cid:29) (cid:104) ˆ I † ˆ I (cid:105) ( τ, β ) = 8 α η c γ OM ( γ + β ) ¯ nγ i (cid:104) βγ (1 + e − γτ ) + 1 (cid:105) + α . (B14)The ﬁrst set of terms in the brackets represent contri-butions from the state and the input thermal noise. Thelast term denotes a constant noise ﬂoor due to the opticalvacuum noise (shot noise). Notice that the ﬂuctuationsat 0 delay roughly double compared to their large-delayvalue in the limit where β (cid:29) γ . In our experiment, thematched ﬁlter bandwidth is limited by preceding ﬁlters.If β is made too large, then the shot-noise term (andother added technical noise contributions) begin to dom-inate the measurement signal. Appendix C: Detector and Device Characterization

In this section we present details related to our hetero-dyne down-conversion scheme, and device characteriza- tion of the mechanical mode in our experiment.We perform down-conversion of the mechanical sig-nal by shifting the frequency of our laser using single-sideband suppressed-carrier modulation. Such a mod-ulation scheme avoids added shot-noise from sidebandsthat do not contribute to signal gain. We implementthe single-sideband modulation using a quadrature phaseshift keying (QPSK) modulator (Optilab QPSK-OM-23).Figure 5 describes the relevant tones in frequency do-main. In this diagram, the vertical arrows denote theprincipal laser frequencies used in the experiment, where α LO denotes the local oscillator tone generated frommodulating the carrier, α carrier denotes the carrier toneprior to any modulation, and α probe denotes the probewhich is sent to the optomechanical crystal. The probeis up-shifted from the carrier by ω AOM / π ≈

40 MHz us-ing an acousto-optic modulator. In Fig. 5a, we show theschematic for phonon-addition, where the probe is blue-detuned from the cavity frequency ω c by a mechanicalfrequency ω m . In order to down-convert the mechanicalsignal to a chosen intermediate frequency ∆ IF , we gener-ate a down-shifted RF drive at the frequency: ω RF = ω m − ∆ IF − ω AOM (C1)Similarly, for phonon-subtraction and red-side drivingof the optical cavity, inspection of Fig. 5b gives the re-quired RF frequency for upconversion: ω RF = ω m − ∆ IF + ω AOM . (C2)We choose ∆ IF / π ≈

20 MHz throughout our experi-ment. A measurement of the mechanical noise spectrum (a)(b)

FIG. 5. Heterodyne scheme: generating single-sideband lo-cal oscillator for detecting mechanical motion. (a) Schematicshowing the down-conversion scheme with probe (vertical redarrow) tuned on the blue side of the optical cavity (grey dot-ted curve). The local oscillator tone, α LO , beats against themechanical signal generating a signal at ∆ IF . (b) Schematicfor the case where the probe tone is red-detuned from cavityresonance. down-shifted to ∆ IF is shown in the red curve of Fig. 6.We use the residual beating tone between the probe andcarrier, at 40 MHz, to balance the optical paths in ourdetection, and to optimize polarization of the local oscil-lator.To determine the optomechanical coupling rate wemeasure mechanical linewidth narrowing, due to back-action, versus laser power. The result is shown in Fig. 7. Appendix D: Phonon-subtraction Results

In the main text we described the results of phonon-subtraction using radially binned histograms. In Figure 8we show the two dimensional quadrature histogram forthe post-selected experimental data along with a theoryﬁt.

Appendix E: Radial MaxLik Tomography1. Deﬁnition of Radial POVMs

The POVM (positive-operator-value-measure) opera-tors that describe heterodyne measurement results aregiven by: ˆΠ( α ) = 1 π | α (cid:105)(cid:104) α | (E1) Frequency (MHz) -60-50-40-30-20-100 P o w e r ( d B m ) Electronic noise floorSignal ONSignal OFF-10 -9 -8 -7 -6 -5 -4

LO power (dBm) -55-54-53-52-51-50-49-48 O p t i c a l N o i s e P o w e r ( d B m ) (a)(b) FIG. 6. (a) Noise power spectra obtained at the output ofthe balanced photodetector. Total power in a 25 kHz band-width is shown vs. frequency. The blue curve shows the noiselevel when the optical local oscillator is oﬀ. With the localoscillator on, but the signal port blocked, we observe an in-crease in the noise ﬂoor due to shot-noise (yellow curve). Thered curve shows the heterodyned signal. Mechanical thermalnoise from the device is visible as a peak at 20 MHz, wherethe black line denotes a Lorentzian ﬁt. The tone at 40 MHzis due to the probe signal beating against the residual opticalcarrier. (b) Scaling of the the noise ﬂoor, with the electronicnoise subtracted, vs local oscillator power. The power law ﬁtgives p ≈ . for all complex α . However, for a measurement apparatussuch as ours which is insensitive to phase, we deﬁne a setof rotationally symmetric POVMs via phase averagingwith ˆΠ RS ( r ) r = (cid:90) π ˆΠ( α ) dφ (E2)where we take α = re iφ .The resulting POVM operators corresponding to mea-surements | α | = r are diagonal with the elements: (cid:104) n | ˆΠ RS ( r ) | n (cid:105) = 2 e − r r n +1 n ! . (E3)Note that we have deﬁned Π RS ( r ) such that the followingcondition required by POVMs is satisﬁed: (cid:90) ∞ ˆΠ RS ( r ) dr = ˆ I (E4)where ˆ I is the identity matrix.In a numerical implementation, the POVM must beevaluated on a discrete vector of radius points, and0 Intracavity Photon Number M e c han i c a l L i ne w i d t h ( M H z ) FIG. 7. Mechancial linewidth versus intracavity photon num-ber obtained from a laser power sweep. The cavity is pumpedon the blue side. From the ﬁt we extract the optomechanicalcoupling rate: g / π = (1 . ± .

02) MHz. -100 0 100-150-100-50050100150 -100 0 100 05001000150020002500

XXP Experiment Theory

FIG. 8. Experimental data and theory ﬁt for quadrature his-tograms obtained after phonon-subtraction. then numerically integrated to generate a probability foreach radial bin. This must be done before applying amaximum-likelihood tomography technique, as describedin the following sections. Full matrix multiplication ismade unnecessary by the choice of radial POVMs, andboth the POVMs and the density matrix can be repre-sented by vectors in a numerical implementation.Using these POVMs we can write a likelihood functionfor generating a given dataset in the experiment. Givena set of measurement results { p i } where each p i denotesthe experimentally measured probability for a result lyingin the i th radial bin. More precisely, given a set of IQdatapoints { s } , we construct the radially binned datasetby simply computing: p i = P (cid:16) r i < (cid:112) Re[ s ] + Im[ s ] < r i +1 (cid:17) (E5)where the smallest radial bin edge r = 0 and the largest, r N bin − = R max is chosen to capture the largest magni-tudes in the data.The cost function, C ,to be minimized via the densitymatrix elements (the optimization parameters) of ρ is thenegative log of the likelihood function L constructed fromthe POVMs and the experimentally observed data: C = − log( L ) = − N bin − (cid:88) i =0 p i log Tr (cid:104) ˆΠ i ρ (cid:105) (E6)where we have used the shorthand:ˆΠ i = (cid:90) r i +1 r i ˆΠ RS ( r ) dr (E7)In other words, the set { Π i } are the binned POVMs.Minimizing Eq. E6 via the set of parameters { ˆ ρ n,n } can in principle be done using a variety of methods. Inthis work, we use the iterative MaxLik algorithm that hasfound widespread use in continuous variable tomographyexperiments [30, 34]. We form the R matrices using thequantitites deﬁned above as: R = (cid:88) i p i Tr (cid:104) ˆΠ i ρ (cid:105) ˆΠ i (E8)The binned POVMs described by Eq. E7 satisfy therequirement: (cid:88) i ˆΠ i = I (E9)Of course, an inﬁnite fock space is required for this tobe true numerically, but this is not an issue in the radialestimation scenario here because we have a-priori knowl-edge on the size of the states, and can choose N fock (cid:29) ¯ n .Starting from a uniform initial guess: ρ = N fock I we update ρ on each step via: ρ → ˆ Rρ ˆ R Tr (cid:104) ˆ Rρ ˆ R (cid:105) (E10)In practice, a learning rate (step-size) is chosen to di-lute R and improve convergence [45]. In this case, R isreplaced: R → I + (cid:15)R where I is the identity matrix,and (cid:15) is a learning rate. The exact value of the rate isunimportant, but a value that is too large results in os-cillations in the log-likelihood value vs. iteration. We set (cid:15) = 10 − .Finally, we stop the sovler when the trace distance be-tween successive ρ falls below a threhsold, ε : (cid:107) ˆ ρ k − ˆ ρ k − (cid:107) < ε (E11)We obtain good results with ε = 10 − , veriﬁed using aknown simulated thermal state as input to the solver.1 λ ~ QPSK Path MatchingQuadrature Hybrid4 GHz IQ BiasController90%1%90% -

20 MHzDigitizer BHDMatchedFilterIQ Samples

FM in FFP

TEC

EOM MS VidiaITLA λ -meter EOMSNA OUT12 GHz AOM VOA40 MHz MS PMEL SNSPDELEL

TECTEC

FPC Amplifier/TTL conversionTrigger A (single photon)Trigger B (uncorrelated) PDEDFAVOA

VOC

OMC 300 K

SNA INLaser Stabilization Source SelectionSingle Sideband LOGeneration High Speed Detection Low SpeedDetectionSideband Filtering IsolatorServo

FIG. 9. Detailed diagram of optomechanical tomography setup. EOM: Electro-optic modulator. FFP: Fiber fabry-perot ﬁlter(see main text for details). EL: inline eigenlight power tap. MS: MEMS optical switch (the BAR state is solid green; theCROSS state is dotted black). FPC: Fiber polarization controller. SNA: Scalar network analyzer. VOA: Variable opticalattenuator. VOC: Variable optical coupler. QPSK: Quadrature phase shift keying modulator. TEC: Thermo-electric coolercontroller. PD: Photodiode. PM: Power meter. EDFA: Erbium doped ﬁber ampliﬁer. SNSPD: Superconducting nanowiresingle photon detector. BHD: Balanced heterodyne detection.

2. POVMs for states with added noise

The POVMs given above allow reconstruction of quan-tum states assuming no additional uncorrelated techni-cal noise has been added to the signal. Following thetreatment in [30, 35, 46], added technical noise can becaptured by modifying Eq. E1 accordingly:ˆΠ( α ) = 1 π ˆ D ( α )ˆ ρ th ˆ D † ( α ) (E12) Here ˆ D ( α ) is the displacement operator, and ˆ ρ th is de-ﬁned as in Eq. A1 where the bath temperature ¯ n is re-placed by the added eﬀective noise power, in units ofphonons. i.e. ¯ n = n added . As described in the main text,in our setup this is estimated by blocking the signal port,and measuring the noise that results. Note in particularthat setting n added = 0 gives ˆ ρ th = | (cid:105)(cid:104) | and recoversEq. E1 from Eq. E12. One can obtain the radially sym-metric POVMs in Eq. E2 by phase averaging over thePOVMs in Eq. E12.2In our estimation problem, we require N fock (cid:29) ¯ n . Fora 4 GHz oscillator at room temperature, ¯ n ≈ · there-fore requiring about N fock ≈ · to reliably reconstructthe state. Thus, given this enormous Hilbert space di-mension, directly evaluating the displacement operatorsin Eq. E12 is prohibitive numerically. Since the radialPOVMs are diagonal matrices, we do not require the fullcomputation of the displacement operators. In the nextsection we provide a formula for the on-diagonal elementsof Eq. E12.

3. Photon number statistics of a displaced thermalstate