# Supplemental Material:
Ponderomotive squeezing of light by a levitated nanoparticle in free space

Andrei Militaru  These authors contributed equally to this work.  Photonics
Laboratory, ETH Zürich, CH-8093 Zürich, Switzerland  Massimiliano Rossi  These
authors contributed equally to this work.  Photonics Laboratory, ETH Zürich,
CH-8093 Zürich, Switzerland  Felix Tebbenjohanns  present address: Department
of Physics, Humboldt-Universität zu Berlin, 10099 Berlin, Germany  Photonics
Laboratory, ETH Zürich, CH-8093 Zürich, Switzerland  Oriol Romero-Isart
Institute for Quantum Optics and Quantum Information of the Austrian Academy
of Sciences, A-6020, Innsbruck, Austria  Institute for Theoretical Physics,
University of Innsbruck, A-6020 Innsbruck, Austria  Martin Frimmer  Photonics
Laboratory, ETH Zürich, CH-8093 Zürich, Switzerland  Lukas Novotny  Photonics
Laboratory, ETH Zürich, CH-8093 Zürich, Switzerland  Quantum Center, ETH
Zürich, CH-8093 Zürich, Switzerland

###### Contents

  1. I Theory
    1. I.1 System of interest
    2. I.2 Optomechanical Hamiltonian
      1. I.2.1 Tweezer Hamiltonian: optical potential
      2. I.2.2 Recoil Hamiltonian: quantum optical scattering
      3. I.2.3 Interacting modes
      4. I.2.4 Information patterns
      5. I.2.5 Total Hamiltonian
    3. I.3 Quantum Langevin equations
      1. I.3.1 Heisenberg measurement-disturbance relation
      2. I.3.2 Coupling to a thermal bath
  2. II Experimental details
    1. II.1 Systematic error and calibration
    2. II.2 Fitting the spectra
    3. II.3 Wigner quasidistribution reconstruction

##  I Theory

In this section, we give a complete derivation of the theory used to analyse
our experimental data. This theory is based on a quantum optics theoretical
framework which assumes a linearized dipolar light-matter interaction. The
derivation is based on Refs. Chang  _et al._ (2010); Romero-Isart  _et al._
(2011); Rodenburg  _et al._ (2016); Gonzalez-Ballestero  _et al._ (2019);
Magrini  _et al._ (2021).

###  I.1 System of interest

We consider a lossless and isotropic dielectric nanoparticle of mass $m$ and
real scalar polarizability $\alpha$. The center-of-mass (COM) motion of the
nanoparticle is described by its momentum operator
$\bm{p}=(p_{x},p_{y},p_{z})^{T}$ and by its corresponding position operator
$\mathbf{r}=(x,y,z)^{T}$. In the limit of a particle much smaller than any
wavelength of interest, the interaction between the nanoparticle and the
surrounding electromagnetic field is given by the dipole interaction
$-\bm{d}(\mathbf{r})\cdot\bm{E}(\mathbf{r})/2$, with
$\bm{d}(\mathbf{r})=\alpha\bm{E}(\mathbf{r})$ the induced dipole moment and
$\bm{E}(\mathbf{r})$ the electric field operator evaluated at the
nanoparticle’s position. The total Hamiltonian can then be described by three
terms,

$H=H_{\mathrm{np}}+H_{\mathrm{f}}+H_{\mathrm{I}},$ (S1a) with the nanoparticle
Hamiltonian $H_{\mathrm{np}}=\frac{\bm{p}^{2}}{2m},$ (S1b) the free field
Hamiltonian
$H_{\mathrm{f}}=\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\hbar\omega a_{\varepsilon}^{\dagger}(\mathbf{k})a_{\varepsilon}(\mathbf{k}),$
(S1c) and the interaction Hamiltonian
$H_{\mathrm{I}}=-\frac{1}{2}\alpha\bm{E}^{2}(\mathbf{r}).$ (S1d)

In Eq. (S1c), $\hbar$ is the reduced Planck constant, $\varepsilon$ is the
polarization index, $\omega=c|\mathbf{k}|$ represents the angular frequency of
each free field mode, and $a_{\varepsilon}(\mathbf{k})$ and
$a^{\dagger}_{\varepsilon}(\mathbf{k})$ are the corresponding ladder operators
that obey standard bosonic commutation relations:
$[a_{\varepsilon}(\mathbf{k}),a_{\varepsilon^{\prime}}(\mathbf{k}^{\prime})]=[a^{\dagger}_{\varepsilon}(\mathbf{k}),a^{\dagger}_{\varepsilon^{\prime}}(\mathbf{k}^{\prime})]=0$
and
$[a_{\varepsilon}(\mathbf{k}),a^{\dagger}_{\varepsilon^{\prime}}(\mathbf{k}^{\prime})]=\delta(\mathbf{k}-\mathbf{k}^{\prime})\delta_{\varepsilon\varepsilon^{\prime}}$.
The electric field is evaluated at the nanoparticle’s position operator
$\mathbf{r}$, such that Eq. (S1d) represents indeed an interaction between
$\bm{E}$ and COM motion. At this point, we deploy the basis of plane waves,
such that the annihilation (creation) operators $a_{\varepsilon}(\mathbf{k})$
($a^{\dagger}_{\varepsilon}(\mathbf{k})$) in Eq. (S1c) refer to a plane wave
mode with polarization along the unit vector
$\bm{u}^{(\mathbf{k})}_{\varepsilon}$ and wavevector $\mathbf{k}$.

In many levitodynamics experiments in free space, the nanoparticle is trapped
in the optical potential of a focused laser beam Maurer  _et al._ (2021);
Ricci  _et al._ (2017); Jain  _et al._ (2016); Tebbenjohanns  _et al._ (2019a,
2020); Hebestreit  _et al._ (2018); Hoang  _et al._ (2018). Within this
framework, we can account for the presence of the trapping beam, referred to
in the following as tweezer field, by splitting the electric field operator
$\bm{E}$ into a strong (classical and time dependent) component
$\bm{E}_{\mathrm{tw}}(\mathbf{r},t)$ displaced by a coherent amplitude and the
quantum fluctuations $\bm{E}_{\mathrm{f}}(\mathbf{r})$ (taken in the
Schrödinger picture and hence time independent) Gonzalez-Ballestero  _et al._
(2019):

$\bm{E}(\mathbf{r},t)=\bm{E}_{\mathrm{tw}}(\mathbf{r},t)+\bm{E}_{\mathrm{f}}(\mathbf{r}).$
(S2a) The free field in Eq. (S2a) takes the standard form
$\bm{E}_{\mathrm{f}}(\mathbf{r})=\text{\rm
i}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\sqrt{\frac{\hbar\omega}{16\pi^{3}\varepsilon_{0}}}\
\bm{u}^{(\mathbf{k})}_{\varepsilon}\left(a_{\varepsilon}(\mathbf{k})\text{\rm
e}^{\text{\rm
i}\mathbf{k}\cdot\mathbf{r}}-a^{\dagger}_{\varepsilon}(\mathbf{k})\text{\rm
e}^{-\text{\rm i}\mathbf{k}\cdot\mathbf{r}}\right),$ (S2b) with
$\varepsilon_{0}$ being the vacuum permittivity. In this work, the tweezer
field is chosen to be an $x$-polarised Gaussian beam whose focus lies at the
origin of our frame of reference. In the neighbourhood of the origin, the
tweezer field can be written as Novotny (2017)
$\bm{E}_{\mathrm{tw}}(\mathbf{r},t)=\frac{E_{0}\
\exp\left(-\frac{x^{2}}{w_{x}^{2}}-\frac{y^{2}}{w_{y}^{2}}\right)}{\sqrt{1+\left(\frac{z}{z_{\mathrm{R}}}\right)^{2}}}\
\cos\big{(}k_{0}z-\mathrm{atan}(z/z_{\mathrm{R}})-\omega_{0}t\big{)}\bm{u}_{x}.$
(S2c)

In Eq. (S2c), $E_{0}$ is the field’s amplitude at the origin, $w_{x}$ and
$w_{y}$ are the beam radii along the axes $x$ and $y$,
$k_{0}=\omega_{0}/c=2\pi/\lambda_{0}$ is the wave number of the beam,
$z_{\mathrm{R}}$ is an effective Rayleigh range, and $\bm{u}_{x}$ is the unit
vector pointing along $x$.

The interaction Hamiltonian in Eq. (S1d) contains a term quadratic in the
tweezer field $\bm{E}_{\mathrm{tw}}$, a mixed term and a term quadratic in the
fluctuation field $\bm{E}_{\mathrm{f}}$. We can simplify the interaction
hamiltonian by retaining only the terms that depend on the tweezer field,
which is assumed to be much stronger than the free field
$\bm{E}_{\mathrm{f}}$. By neglecting the term quadratic in the free field, we
are not taking into account particle-mediated interactions between different
electromagnetic modes. This approximation is justified because we consider the
limit of an infinitely small nanoparticle Pflanzer  _et al._ (2012). For
finite sized objects, instead, cross terms should be included in Eq. (S1c)
Pflanzer  _et al._ (2012), or a more suitable scattering mode basis should be
developed (as done in Ref. Maurer  _et al._ (2021)). The simplified
interaction Hamiltonian is now

$H_{\mathrm{I}}\approx-\frac{1}{2}\alpha\bm{E}_{\mathrm{tw}}\cdot\bm{E}_{\mathrm{tw}}-\alpha\bm{E}_{\mathrm{tw}}\cdot\bm{E}_{\mathrm{f}}=:H_{\mathrm{tw}}+H_{\mathrm{rc}}.$
(S3)

We have defined the first and second terms of Eq. (S3) respectively as tweezer
Hamiltonian $H_{\mathrm{tw}}$ and recoil Hamiltonian $H_{\mathrm{rc}}$. The
tweezer Hamiltonian is governed by the classical component of the
electromagnetic field, while the quantum properties of light are carried by
the recoil Hamiltonian. In Secs. I.2.1 and I.2.2 we will show that the tweezer
Hamiltonian is responsible for generating the optical trapping potential while
the recoil Hamiltonian accounts for the light scattered by the particle and
for the related photon recoil.

###  I.2 Optomechanical Hamiltonian

In this section we manipulate the separate terms of Eqs. (S1) and derive a
form of the Hamiltonian that is simple and that can be of direct use for many
experiments where the small particle limit is valid.

####  I.2.1 Tweezer Hamiltonian: optical potential

The tweezer Hamiltonian is straightforward to analyse. When inserting Eq.
(S2c) into $H_{\mathrm{tw}}$ from Eq. (S3) we get

$H_{\mathrm{tw}}=-\frac{\frac{1}{2}\alpha
E_{0}^{2}}{1+\left(\frac{z}{z_{\mathrm{R}}}\right)^{2}}\cos^{2}\big{(}k_{0}z-\mathrm{atan}(z/z_{\mathrm{R}})-\omega_{0}t\big{)}\exp\left(-2\left(\frac{x^{2}}{w_{x}^{2}}+\frac{y^{2}}{w_{y}^{2}}\right)\right).$
(S4)

In most experiments, feedback cooling is applied on the COM motion, such that
the particle’s displacement from the tweezer’s focus is very small and we can
expand Eq. (S4) to second order. Once we expand Eq. (S4), we can neglect the
terms oscillating at twice the tweezer frequency, in other words we apply the
rotating wave approximation. We are then left with

$H_{\mathrm{tw}}\approx\frac{1}{2}m\left(\Omega_{x}^{2}x^{2}+\Omega_{y}^{2}y^{2}+\Omega_{z}^{2}z^{2}\right),$
(S5a) where we defined the three resonance angular frequencies
$\Omega_{x}^{2}=\frac{\alpha
E_{0}^{2}}{mw_{x}^{2}},\qquad\Omega_{y}^{2}=\frac{\alpha
E_{0}^{2}}{mw_{y}^{2}},\qquad\Omega_{z}^{2}=\frac{\alpha
E_{0}^{2}}{2mz_{\mathrm{R}}^{2}}.$ (S5b)

Equation (S5a) represents the harmonic optical potential generated by the
trapping laser. Equation (S5a) tells us that the COM motion of the
nanoparticle represents a set of three independent harmonic oscillators, one
for each axis of motion, whose respective angular frequencies are proportional
to the input beam’s amplitude $E_{0}$ and inversely proportional to the beams
focal spot size, Eq. (S5b).

####  I.2.2 Recoil Hamiltonian: quantum optical scattering

The recoil Hamiltonian
$H_{\mathrm{rc}}=-\alpha\bm{E}_{\mathrm{tw}}\cdot\bm{E}_{\mathrm{f}}$ in Eq.
(S3) represents the interaction between the trapping beam
$\bm{E}_{\mathrm{tw}}$ and the free field $\bm{E}_{\mathrm{f}}$, which is
mediated by the nanoparticle. This interaction lies at the heart of most
levitodynamics experiments: on the one hand, every photon that is rerouted by
the particle from the tweezer field to the free field contributes to disturb
the particle’s motion through photon recoil. On the other hand, as we will see
in Sec. I.2.3, measuring the field scattered by the particle allows us to
reconstruct its position. The precision of the position measurement is then
linked to the photon recoil through the Heisenberg measurement-disturbance
relation.

We consider the common case in which the displacements of the nanoparticle are
very small compared to the trapping beam’s wavelength. We can then simplify
both the tweezer and the free fields in Eqs. (S2b) and (S2c) by expanding the
complex exponentials to first order:

$\displaystyle\bm{E}_{\mathrm{f}}(\mathbf{r})$
$\displaystyle\approx\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\sqrt{\frac{\hbar\omega}{16\pi^{3}\varepsilon_{0}}}\ \left[\text{\rm
i}\left(a_{\varepsilon}(\mathbf{k})-a^{\dagger}_{\varepsilon}(\mathbf{k})\right)-\mathbf{k}\cdot\mathbf{r}\left(a_{\varepsilon}(\mathbf{k})+a^{\dagger}_{\varepsilon}(\mathbf{k})\right)\right]\bm{u}^{(\mathbf{k})}_{\varepsilon},$
(S6a) $\displaystyle\bm{E}_{\mathrm{tw}}(\mathbf{r},t)$ $\displaystyle\approx
E_{0}\left(\cos\big{(}\omega_{0}t\big{)}+Ak_{0}z\sin\big{(}\omega_{0}t\big{)}\right)\bm{u}_{x}.$
(S6b)

Here, we introduced the geometric factor $A=1-1/(k_{0}z_{\mathrm{R}})$ to
account for the Gouy phase in the neighbourhood of the origin. We then
substitute Eqs. (S6) into the recoil Hamiltonian $H_{\mathrm{rc}}$ in Eq. (S3)
and retain only terms linear in the position. The result reads

$\displaystyle H_{\mathrm{rc}}$ $\displaystyle\approx-\text{\rm
i}\cos\big{(}\omega_{0}t\big{)}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
G(\mathbf{k},\varepsilon)\left(a_{\varepsilon}(\mathbf{k})-a^{\dagger}_{\varepsilon}(\mathbf{k})\right)$
$\displaystyle\qquad+\cos\big{(}\omega_{0}t\big{)}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
G(\mathbf{k},\varepsilon)\
\left(a_{\varepsilon}(\mathbf{k})+a^{\dagger}_{\varepsilon}(\mathbf{k})\right)\mathbf{k}\cdot\mathbf{r}$
$\displaystyle\qquad-\text{\rm
i}\sin\big{(}\omega_{0}t\big{)}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
G(\mathbf{k},\varepsilon)\
\left(a_{\varepsilon}(\mathbf{k})-a^{\dagger}_{\varepsilon}(\mathbf{k})\right)Ak_{0}z.$
(S7)

In Eq. (I.2.2), the function $G(\mathbf{k},\varepsilon)$ is defined in analogy
to Ref. Gonzalez-Ballestero  _et al._ (2019) as

$G(\mathbf{k},\varepsilon)=\alpha
E_{0}\sqrt{\frac{\hbar\omega}{16\pi^{3}\varepsilon_{0}}}(\bm{u}^{(\mathbf{k})}_{\varepsilon}\cdot\bm{u}_{x}),$
(S8)

and expresses the strength of the coupling between the tweezer and each plane
wave mode as a function of the wavevector $\mathbf{k}$ and of the polarization
$\varepsilon$. The time dependence of Eq. (I.2.2) can be eliminated by moving
to a frame rotating at the tweezer frequency, i.e., we substitute
$a_{\varepsilon}(\mathbf{k})\to
a_{\varepsilon}(\mathbf{k})\exp\left(-\text{\rm i}\omega_{0}t\right)$ and
$a^{\dagger}_{\varepsilon}(\mathbf{k})\to
a^{\dagger}_{\varepsilon}(\mathbf{k})\exp\left(\text{\rm
i}\omega_{0}t\right)$, and neglecting the fast rotating terms:

$\displaystyle H_{\mathrm{rc}}$ $\displaystyle\approx-\text{\rm
i}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\frac{G(\mathbf{k},\varepsilon)}{2}\left(a_{\varepsilon}(\mathbf{k})-a^{\dagger}_{\varepsilon}(\mathbf{k})\right)$
$\displaystyle\quad-\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\frac{G(\mathbf{k},\varepsilon)}{2}\
\left(a_{\varepsilon}(\mathbf{k})+a^{\dagger}_{\varepsilon}(\mathbf{k})\right)\left(\mathbf{k}\cdot\mathbf{r}-Ak_{0}z\right).$
(S9)

Equation (I.2.2) is the main result of this subsection. It consists of two
terms. The first one, position independent, displaces the free field modes to
a coherent state whose strength as a function of orientation is proportional
the coupling function $G(\mathbf{k},\varepsilon)$. This first term represents
the quantum description of the Rayleigh scattering of the nanoparticle. The
second term describes instead a position dependent displacement of the free
field, whose radiation pattern is different from the Rayleigh scattering due
to the presence of the term $\mathbf{k}$. In the next section, we explicitly
solve the summation over all modes and show that only three free space modes
are sufficient to describe the optomechanical interaction present in Eq.
(I.2.2).

####  I.2.3 Interacting modes

So far in this work, we have used the plane-wave basis for the free space
modes. Within our small-particle approximation, this plane-wave basis provides
a convenient mathematical tool that has allowed us to derive the simplified
recoil Hamiltonian of Eq. (I.2.2). In Eq. (I.2.2), every COM degree of freedom
is coupled to every mode $a_{\varepsilon}(\mathbf{k})$ with a coupling
strength that depends on the orientation of the plane wave and on its
frequency. In this subsection, we effectively solve the integrals over plane
wave orientation in Eq. (I.2.2) and introduce a new mode basis, consisting of
what we define as interacting modes. We use the interacting modes to derive
the main theoretical result of this work.

When considering Eq. (I.2.2), we can see that both terms involve a summation
over all plane-wave modes with different weights. For the $x$, $y$ and $z$
dependent terms, we find weights given by $G(\mathbf{k},\varepsilon)k_{x}$,
$G(\mathbf{k},\varepsilon)k_{y}$, and
$G(\mathbf{k},\varepsilon)(k_{z}-k_{0}A)$ respectively. It is then appealing
to define interacting modes as

$a_{j}(\omega)\propto\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
G(\mathbf{k},\varepsilon)\delta\left(k-\frac{\omega}{c}\right)\tilde{k}_{j}a_{\varepsilon}(\mathbf{k}),$
(S10)

where the Dirac $\delta$ selects the modes with angular frequency $\omega$ and
$\tilde{k}_{j}=k_{j}-A\delta_{zj}$. By choosing suitable numerical prefactors
to guarantee

$\left[a_{j}(\omega),a^{\dagger}_{j^{\prime}}(\omega^{\prime})\right]=\delta_{jj^{\prime}}\delta(\omega-\omega^{\prime}),$
(S11)

and using Eq. (S8) we reach the expressions

$\displaystyle a_{x}(\omega)$
$\displaystyle=-\sqrt{\frac{15c}{8\pi\omega^{2}}}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\delta\left(k-\frac{\omega}{c}\right)\frac{k_{x}}{k}\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)a_{\varepsilon}(\mathbf{k}),$
(S12a) $\displaystyle a_{y}(\omega)$
$\displaystyle=-\sqrt{\frac{15c}{16\pi\omega^{2}}}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\delta\left(k-\frac{\omega}{c}\right)\frac{k_{y}}{k}\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)a_{\varepsilon}(\mathbf{k}),$
(S12b) $\displaystyle a_{z}(\omega)$
$\displaystyle=-\sqrt{\frac{15c}{8\pi\omega^{2}(2+5A^{2})}}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\delta\left(k-\frac{\omega}{c}\right)\left(\frac{k_{z}-k_{0}A}{k}\right)\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)a_{\varepsilon}(\mathbf{k}),$
(S12c)

where the negative signs have been chosen for later notational convenience.
With Eqs. (S12), the commutator condition in Eq. (S11) is indeed fulfilled. We
show this explicitly for the $x$ interacting mode (the verification of the
other modes proceeds analogously). Using
$[a_{\varepsilon}(\mathbf{k}),a_{\varepsilon^{\prime}}(\mathbf{k}^{\prime})]=\delta_{\varepsilon\varepsilon^{\prime}}\delta(\mathbf{k}-\mathbf{k}^{\prime})$,
we get

$\displaystyle[a_{x}(\omega),a_{x}^{\dagger}(\omega^{\prime})]$
$\displaystyle=\frac{15c}{8\pi\omega^{2}}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\delta\left(k-\frac{\omega}{c}\right)\delta\left(k-\frac{\omega^{\prime}}{c}\right)\frac{k_{x}^{2}}{k^{2}}\sum_{\varepsilon}\left(\bm{u}_{x}\cdot\bm{u}_{\varepsilon}^{(\mathbf{k})}\right)$
$\displaystyle=\delta(\omega-\omega^{\prime})\frac{15c^{2}}{8\pi\omega^{2}}\int_{4\pi}\mathrm{d}\Omega\,k_{x}^{2}\sum_{\varepsilon}\left(\bm{u}_{x}\cdot\bm{u}_{\varepsilon}^{(\mathbf{k})}\right),$
(S13)

where $\mathrm{d}\Omega$ represents the infinitesimal solid angle. We can
solve the sum over the polarization degrees of freedom using the closure
relation Gonzalez-Ballestero  _et al._ (2019)

$\sum_{\varepsilon}\left(\bm{u}^{(\mathbf{k})}_{\varepsilon}\cdot\bm{v}\right)^{2}=\bm{v}\cdot\bm{v}-\frac{(\bm{v}\cdot\mathbf{k})^{2}}{k^{2}},$
(S14)

with the help of which we write Eq. (I.2.3) as

$\displaystyle[a_{x}(\omega),a_{x}^{\dagger}(\omega^{\prime})]$
$\displaystyle=\delta(\omega-\omega^{\prime})\frac{15}{8\pi}\int_{4\pi}\mathrm{d}\Omega\,\frac{k_{x}^{2}}{k^{2}}\left(1-\frac{k_{x}^{2}}{k^{2}}\right)$
$\displaystyle=\delta(\omega-\omega^{\prime})\frac{15}{8\pi}\int_{0}^{2\pi}\mathrm{d}\varphi\int_{0}^{\pi}\mathrm{d}\theta\,\sin^{3}(\theta)\cos^{2}(\varphi)\left(1-\cos^{2}(\varphi)\sin^{2}(\theta)\right)$
$\displaystyle=\delta(\omega-\omega^{\prime}),$ (S15)

where we made use of the polar and azimuthal spherical coordinates,
respectively $\theta$ and $\varphi$.

In addition to the interacting modes of Eqs. (S12), the first line of Eq.
(I.2.2) excites the $x$-oriented dipole mode

$a_{0}(\omega)=\sqrt{\frac{3c}{8\pi\omega^{2}}}\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
\delta\left(k-\frac{\omega}{c}\right)\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)a_{\varepsilon}(\mathbf{k}),$
(S16)

which describes Rayleigh scattering, as we will show in Sec. I.2.4. While the
dipole mode in Eq. (S16) does not affect the COM degrees of freedom directly
[the first line of Eq. (I.2.2) contains no position-dependent factor], it does
impact the dynamics of the $z$ motion of the nanoparticle. Its effect can be
understood by noticing that 111We define
$a^{\dagger}_{j}(\omega)=\left[a_{j}(\omega)\right]^{\dagger}$

$\left[a_{0}(\omega),a^{\dagger}_{z}(\omega^{\prime})\right]=\sqrt{\frac{5A^{2}}{2+5A^{2}}}\
\delta(\omega-\omega^{\prime}),$ (S17)

in other words the $z$-interacting mode and the Rayleigh scattering mode are
not orthogonal. In order to limit our analysis to the COM motion and to the
interacting modes only, we can decompose the Rayleigh scattering mode in two
components,

$a_{0}(\omega)=\sqrt{\frac{5A^{2}}{2+5A^{2}}}\
a_{z}(\omega)+g_{\mathrm{n}}a_{\mathrm{n}}(\omega),$ (S18)

where the first one accounts for Eq. (S17) and the second one—whose exact
expression is irrelevant for our purpose—is orthogonal to all interacting
modes and can be neglected because it does not interact with the COM motion.

Using the interacting modes in Eqs. (S12) and the coupling functions

$\displaystyle g_{x}(\omega)$ $\displaystyle=\frac{\alpha
E_{0}k}{\hbar}\sqrt{\frac{\pi\hbar\omega^{3}}{15\varepsilon_{0}(2\pi
c)^{3}}},$ (S19a) $\displaystyle g_{y}(\omega)$ $\displaystyle=\frac{\alpha
E_{0}k}{\hbar}\sqrt{\frac{2\pi\hbar\omega^{3}}{15\varepsilon_{0}(2\pi
c)^{3}}},$ (S19b) $\displaystyle g_{z}(\omega)$ $\displaystyle=\frac{\alpha
E_{0}k}{\hbar}\sqrt{\frac{(2+5A^{2})\pi\hbar\omega^{3}}{15\varepsilon_{0}(2\pi
c)^{3}}},$ (S19c) $\displaystyle g_{\mathrm{rp}}(\omega)$
$\displaystyle=\frac{\alpha
E_{0}k}{\hbar}\sqrt{\frac{\pi\hbar\omega^{3}}{3\varepsilon_{0}(2\pi
c)^{3}}}\sqrt{\frac{5A^{2}}{2+5A^{2}}}.$ (S19d)

we can write then the recoil Hamiltonian simply as

$\displaystyle H_{\mathrm{rc}}=$
$\displaystyle\sum_{j\in\\{x,y,z\\}}\hbar\int_{\mathbb{R}}\mathrm{d}\omega\
\left[-g_{j}(\omega)r_{j}\left(a_{j}(\omega)+a_{j}^{\dagger}(\omega)\right)+\text{\rm
i}\delta_{zj}\frac{g_{\mathrm{rp}}}{k}(\omega)\left(a_{z}(\omega)-a_{z}^{\dagger}(\omega)\right)\right],$
(S20)

where $r_{j}=\mathbf{r}\cdot\bm{u}_{j}$ represents an alternative notation for
the position operator along the $j$ axis. Note that we extended the lower
bound of the integrals over the angular frequency $\omega$ to $-\infty$ rather
then $-\omega_{0}$. This approximation is convenient and valid as long as the
particle’s dynamics is much slower than the tweezer frequency $\omega_{0}$, in
other words as long as the Markovian approximation is valid Gardiner and
Collett (1985).

The resulting recoil Hamiltonian Eq. (S20) is known in literature as
optomechanical interaction Clerk  _et al._ (2010); Aspelmeyer  _et al._
(2014). Every position degree of freedom $r_{j}$ interacts with infinitely
many modes $a_{j}(\omega)$ (one for every angular frequency $\omega$)
belonging to only one spatial mode, which is the reason why we named these
spatial modes interacting modes. The interaction is such that each interacting
modes is coherently displaced by an amount proportional to the nanoparticle’s
COM displacement, and the proportionality factors are detailed in Eqs. (S19).
Moreover, the $z$ interacting mode presents a position-independent populating
term [second part of Eq. (S20)], which in turn affects $z$ through their
optomechanical interaction. This position-independent term is responsible for
the radiation pressure that the tweezer beam exerts on the nanoparticle, as we
will explain in Sec. I.3.

####  I.2.4 Information patterns

We have seen that the interaction between the mechanical degrees of the
nanoparticle and the free electromagnetic field becomes considerably simpler
when described in terms of the interacting modes. In this section, we further
investigate their meaning and physical form.

To this end, let us conduct a thought experiment. We use the spherical
coordinates $\theta$ and $\varphi$, given by the polar angle with respect to
the $z$ axis and the azimuthal angle, respectively. We can then define the
photon number operator $N$ as

$N=\sum_{\varepsilon}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}k\
a^{\dagger}_{\varepsilon}(\mathbf{k})a_{\varepsilon}(\mathbf{k})=\int_{4\pi}\mathrm{d}\Omega\
\underbrace{\sum_{\varepsilon}\int_{\mathbb{R}^{+}}\mathrm{d}k\
k^{2}a^{\dagger}_{\varepsilon}(\mathbf{k})a_{\varepsilon}(\mathbf{k})}_{\mathcal{N}(\theta,\varphi)},$
(S21)

where $\Omega$ refers to the solid angle along $(\theta,\varphi)$. In Eq.
(S21), $N$ quantifies the number of photons present in a state independent
from their frequency, orientation or polarization. As shown in Eq. (S21), it
is possible to define the density of photons $\mathcal{N}$ along a direction
$(\theta,\varphi)$. Integrating $\mathcal{N}$ over the full solid angle yields
then the total number of photons 222The integral over the total solid angle
takes the standard definition
$\int_{4\pi}\mathrm{d}\Omega=\int_{0}^{\pi}\mathrm{d}\theta\
\int_{0}^{2\pi}\mathrm{d}\varphi$. Equipped with the concept of density of
photons $\mathcal{N}$, we can introduce the radiation pattern $\rho$
corresponding to a state $|\psi\rangle$ of the optical field as

$\rho(\theta,\varphi):=\langle\mathcal{N}(\theta,\varphi)\rangle=\langle\psi|\mathcal{N}(\theta,\varphi)|\psi\rangle.$
(S22)

Equation (S22) is our key to understand the interacting modes. We consider the
optical vacuum state $|0\rangle$ and excite one of the interacting modes with
a coherent state. Rather that focusing on a specific frequency of the
interacting modes, $a_{j}(\omega)$, it is beneficial to consider the multimode
operator

$a_{j}^{\dagger}=\int_{\mathbb{R}^{+}}\mathrm{d}\omega\
f(\omega)a_{j}^{\dagger}(\omega),$ (S23)

whose commutator with $a_{\varepsilon}(\mathbf{k})$ is

$[a_{\varepsilon}(\mathbf{k}),a^{\dagger}_{j}]=\frac{c}{N_{j}}\frac{\tilde{k}_{j}}{k}f(kc)\left(\bm{u}_{x}\cdot\bm{u}_{\varepsilon}^{(\mathbf{k})}\right)=:\Lambda_{j}(\varepsilon,\mathbf{k}),$
(S24)

where

$N_{j}^{2}=\frac{8\pi\omega^{2}}{15c}\begin{pmatrix}1\\\ 2\\\
2+5A^{2}\end{pmatrix},$ (S25)

and $\tilde{k}_{j}=k_{j}-\delta_{zj}A$. The weight function $f$, normalized
such that $\int_{\mathbb{R}^{+}}\mathrm{d}\omega\ |f|^{2}=1$, can be chosen
arbitrarily narrow or wide around a frequency of interest. Using the multimode
operator in Eq. (S23), rather than Eqs. (S12), allows us to bypass the
singular nature of the commutator of continuous mode operators. We can now
excite the interacting mode with a coherent state of amplitude $\beta$. Our
excited interacting mode takes then the form

$|\beta\rangle_{j}=\mathcal{D}_{j}(\beta)|0\rangle=\text{\rm e}^{\beta
a_{j}^{\dagger}-\beta^{*}a_{j}}|0\rangle=\text{\rm
e}^{-|\beta|^{2}/2}\text{\rm e}^{\beta a^{\dagger}_{j}}|0\rangle.$ (S26)

We can see from Eq. (S24) and Eq. (S26) that 333Use
$[a_{\varepsilon},a_{j}^{\dagger
n}]=\Lambda_{j}(\varepsilon,\mathbf{k})na_{j}^{\dagger(n-1)}$ and the Taylor
expansion of $\exp(\beta a_{j}^{\dagger})$.

$a_{\varepsilon}(\mathbf{k})|\beta\rangle_{j}=\beta\Lambda_{j}(\varepsilon,\mathbf{k})|\beta\rangle_{j},$
(S27)

which we substitute in Eq. (S22) to obtain

$\displaystyle\rho_{x}$
$\displaystyle=|\beta|^{2}\frac{15}{8\pi}\sum_{\varepsilon}\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)^{2}\frac{k_{x}^{2}}{k^{2}},$
(S28a) $\displaystyle\rho_{y}$
$\displaystyle=|\beta|^{2}\frac{15}{16\pi}\sum_{\varepsilon}\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)^{2}\frac{k_{y}^{2}}{k^{2}},$
(S28b) $\displaystyle\rho_{z}$
$\displaystyle=|\beta|^{2}\frac{15}{8\pi(2+5A^{2})}\sum_{\varepsilon}\left(\bm{u}_{x}\cdot\bm{u}^{(\mathbf{k})}_{\varepsilon}\right)^{2}\left(\frac{k_{z}}{k}-A\right)^{2},$
(S28c)

with the angle dependence implicitly contained in $\mathbf{k}$ and
$\bm{u}_{\varepsilon}^{(\mathbf{k})}$. Using the polarization closure relation
Eq. (S14) we get

$\displaystyle\rho_{x}$
$\displaystyle=|\beta|^{2}\frac{15}{8\pi}\big{(}1-\cos^{2}\varphi\sin^{2}\theta\big{)}\cos^{2}\varphi\sin^{2}\theta,$
(S29a) $\displaystyle\rho_{y}$
$\displaystyle=|\beta|^{2}\frac{15}{16\pi}\big{(}1-\cos^{2}\varphi\sin^{2}\theta\big{)}\sin^{2}\varphi\sin^{2}\theta,$
(S29b) $\displaystyle\rho_{z}$
$\displaystyle=|\beta|^{2}\frac{15}{8\pi(2+5A^{2})}\big{(}1-\cos^{2}\varphi\sin^{2}\theta\big{)}\left(\cos\theta-A\right)^{2}.$
(S29c)

When looking at Eqs. (S29) we see that the radiation patterns of the
interacting modes are not new to the community. In fact, these patterns have
been already derived in a semiclassical calculation Tebbenjohanns  _et al._
(2019b). In the context of Ref. Tebbenjohanns  _et al._ (2019b),
$\rho_{j}(\theta,\varphi)$ represented the information density about the
position $r_{j}$ of a dipolar scatterer contained along the direction
$(\theta,\varphi)$. In this work, these information patterns play the role of
orthonormal electromagnetic modes that interact respectively with the position
$r_{j}$ only. Within the current framework, it is natural that an optimal
measurement of $r_{j}$ must involve the radiation pattern $\rho_{j}$.

A second way to understand Eqs. (S29), which is mathematically simpler but
further away from a classical intuition, is to populate the interacting modes
with one photon, in other words to use $|\psi\rangle=a_{j}^{\dagger}|0\rangle$
in Eq. (S22). The result coincides once again with Eqs. (S29) for $\beta=1$.
The information patterns can be then interpreted as the probability to find a
photon along $(\theta,\varphi)$ when the interacting mode contains one photon.

Finally, we can repeat the same procedure on the dipole mode, Eq. (S16), and
derive its radiation pattern

$\rho_{0}=|\beta|^{2}\frac{3}{8\pi}\big{(}1-\cos^{2}\varphi\sin^{2}\theta\big{)},$
(S30)

which coincides indeed with the dipole radiation from classical
electrodynamics Novotny (2017).

####  I.2.5 Total Hamiltonian

We conclude this section by putting together the results achieved so far. In
Sec. I.2.1 we have seen how the tweezer Hamiltonian generates a harmonic
potential for the nanoparticle. In addition to this, the nanoparticle mediates
a coupling between the positions $r_{j}$ and the free space electromagnetic
mode. Importantly, each position $r_{j}$ interacts with only one spatial
electromagnetic mode when the basis of interacting modes is chosen 444Note
that we use the expression interacting mode to refer to the spatial shape of
the optical modes described by Eqs. (S12). Each interacting mode consists then
of infinitely many frequency modes.. The total optomechanical Hamiltonian that
describes the problem consists thus of three mutually commuting terms:

$H=H_{x}+H_{y}+H_{z},$ (S31a) where each of the Hamiltonians $H_{j}$ takes the
form $\displaystyle H_{j}$
$\displaystyle=\frac{p_{j}^{2}}{2m}+\frac{1}{2}m\Omega_{j}^{2}r_{j}^{2}+\int_{\mathbb{R}}\mathrm{d}\omega\
\hbar\omega a_{j}^{\dagger}(\omega)a_{j}(\omega)-\hbar
g_{j}r_{j}\int_{\mathbb{R}}\mathrm{d}\omega\
\left(a_{j}(\omega)+a^{\dagger}_{j}(\omega)\right)$
$\displaystyle\qquad+\text{\rm
i}\delta_{jz}\hbar\frac{g_{\mathrm{rp}}}{k}\int_{\mathbb{R}}\mathrm{d}\omega\
\left(a_{j}(\omega)-a^{\dagger}_{j}(\omega)\right),$ (S31b)

and where the Markovian approximation has been applied, i.e. the coupling
constants are evaluated at the tweezer frequency: $g_{j}(\omega)\approx
g_{j}(\omega_{0})=g_{j}$(and analogously for $g_{\mathrm{rp}}$). This
approximation is valid as long as the dynamics of the particle is much slower
than the tweezer frequency.

Equation (S31) contains many terms, nevertheless they can be readily
understood. The first three terms contain the harmonic dynamics of the
nanoparticle and the free evolution of the $j$-interacting mode. This
$j$-interacting mode interacts with the position $r_{j}$ through a bilinear
Hamiltonian given by the fourth term. Finally, the last line of Eq. (S31) is
nonzero only for $H_{z}$ and represents a coherent population of the
$z$-interacting mode. Because of the fourth term of Eq. (S31), the coherent
drive translates into a constant force pushing the particle along the $z$
axis, in other words this term describes the radiation pressure (see Sec.
I.3). Within this framework, the radiation pressure can be understood as a
consequence of the fact that the Rayleigh scattering mode has a finite mode
overlap with the $z$-interacting mode.

Equation (S31) is cast in a form that is ideally suited to derive the quantum
Langevin description of the system, which we address in the following section.

###  I.3 Quantum Langevin equations

In this section, we focus the attention to the dynamics related to the motion
along the $z$ axis only. The equations of motion for the transverse modes $x$
and $y$ can be derived in an analogous fashion. From $H_{z}$ in Eq. (S31), we
obtain the following Heisenberg equations of motion

$\displaystyle\dot{z}$ $\displaystyle=p_{z}/m,$ (S32a)
$\displaystyle\dot{p}_{z}$ $\displaystyle=-m\Omega_{z}^{2}z+\hbar
g\int_{\mathbb{R}}\mathrm{d}\omega\
\left(a(\omega)+a^{\dagger}(\omega)\right),$ (S32b)
$\displaystyle\dot{a}(\omega,t)$ $\displaystyle=-\text{\rm i}\omega
a(\omega,t)+\tilde{g}(t),$ (S32c)

where $\tilde{g}(t)=g_{\text{rp}}/k+\text{\rm i}gz(t)$ is the drive of the
field and where we have neglected the $z$ subscripts for conciseness. Assuming
the optical mode to have an initial condition $a(\omega,t_{-})$ in the past
($t_{-}<t$), the solution of Eq. (S32c) is Gardiner and Collett (1985)

$\displaystyle a(\omega,t)=a(\omega,t_{-})\text{\rm e}^{\text{\rm
i}\omega(t_{-}-t)}+\int_{t_{-}}^{t}\mathrm{d}s\ \tilde{g}(s)\text{\rm
e}^{\text{\rm i}\omega(s-t)}.$ (S33)

Equivalently, by choosing the boundary conditions $a(\omega,t_{+})$ in the
future ($t_{+}>t$) we have

$\displaystyle a(\omega,t)=a(\omega,t_{+})\text{\rm e}^{\text{\rm
i}\omega(t_{+}-t)}-\int_{t}^{t_{+}}\mathrm{d}s\ \tilde{g}(s)\text{\rm
e}^{\text{\rm i}\omega(s-t)}.$ (S34)

It is useful at this point to define the input (output) field as

$\displaystyle a_{\mathrm{in}}(t)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{d}\omega\
a(\omega,t_{-})\text{\rm e}^{\text{\rm i}\omega(t_{-}-t)},$ (S35a)
$\displaystyle a_{\mathrm{out}}(t)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{d}\omega\
a(\omega,t_{+})\text{\rm e}^{\text{\rm i}\omega(t_{+}-t)},$ (S35b)

which has the commutation relations
$[a_{i}(t),a^{\dagger}_{j}(t^{\prime})]=\delta_{ij}\delta(t-t^{\prime})$,
where $j=\\{\text{in, out}\\}$. The input (output) field represents the free
evolution of the optical field before (after) the interaction with the
mechanical degrees of freedom occurred and propagated until the present time
$t$ assuming no interaction. For later convenience, we describe the field in
terms of Hermitian quadrature amplitudes (later simply termed quadratures)
according to $a_{\text{out}}=\left(X_{\text{out}}+\text{\rm
i}Y_{\text{out}}\right)/\sqrt{2}$.

By integrating both hand sides of Eqs. (S33) and (S34) and using Eqs. (S35) we
can obtain the following input-output relations for the quadratures Gardiner
and Collett (1985):

$\displaystyle X_{\text{out}}(t)$
$\displaystyle=X_{\text{in}}(t)+\sqrt{2\Gamma_{\mathrm{R}}},$ (S36a)
$\displaystyle Y_{\text{out}}(t)$
$\displaystyle=Y_{\text{in}}(t)+\sqrt{4\Gamma_{\text{qba}}}q(t),$ (S36b)

where we have introduced the scattering rate $\Gamma_{\mathrm{R}}=2\pi
g^{2}_{\text{rp}}/k^{2}$ and the quantum backaction decoherence rate
$\Gamma_{\text{qba}}=2\pi g^{2}_{z}z_{\text{zpf}}^{2}$, with the zero point
fluctuations $z_{\text{zpf}}^{2}=\hbar/(2m\Omega_{z})$ and the dimensionless
position $q=z/(\sqrt{2}z_{\text{zpf}})$.

By substituting Eq. (S33) in Eqs. (S32a) and (S32b), the equations of motion
for the mechanical position and momentum read

$\displaystyle\dot{q}$ $\displaystyle=\Omega_{z}p,$ (S37a)
$\displaystyle\dot{p}$
$\displaystyle=-\Omega_{z}q+\sqrt{2\Gamma_{\text{qba}}\Gamma_{\mathrm{R}}}+\sqrt{4\Gamma_{\text{qba}}}X_{\text{in}}(t),$
(S37b)

where $p=p_{z}/(\sqrt{2}p_{\text{zpf}})$, with
$p_{\text{zpf}}^{2}=m\Omega_{z}\hbar/2$. The tweezer field exerts a constant
radiation pressure force on the particle, which is captured by the second term
in Eq. (S37b). With the use of Eqs. (S19), we can appreciate that this
constant force can be rewritten as $PA/c$, where
$P=(\omega_{0}^{4}\alpha^{2}E_{0}^{2})/(12\pi\varepsilon_{0}c^{3})$ is the
optical power radiated by the induced dipole. In the following, we account for
the radiation pressure by shifting the equilibrium position of the
nanoparticle to
$q_{\text{eq}}=\sqrt{2\Gamma_{\text{qba}}\Gamma_{\text{rp}}}/\Omega_{z}$. The
quantum fluctuations associated with the input optical field are depicted by
the last term in Eq. (S37b), and act as a random force noise . We notice that
there is no dissipation term associated with these fluctuations. This is in
contrast with a relativistic treatment, which predicts the presence of damping
via radiative dissipation Novotny (2017).

Note that the Hamiltonian used in the main text was taken in the interaction
picture with respect to the free optical field. In this picture, with $t_{-}$
chosen as the initial time and neglecting the static radiation pressure, the
Hamiltonian $H^{(\mathrm{I})}$ takes the form

$H^{(\mathrm{I})}=\frac{1}{2}\hbar\Omega_{z}\left(p^{2}+q^{2}\right)-\sqrt{4\Gamma_{\mathrm{qba}}}\
qX_{\mathrm{in}}(t),$ (S38)

whose last term is the interaction Hamiltonian reported in the main text.
Equation (S38) is time dependent, however it commutes with itself at different
times. With Eq. (S38), Eqs (S37) follow directly from the Heisenberg equations
of motion for $q$ and $p$.

####  I.3.1 Heisenberg measurement-disturbance relation

Equations (S36) and (S37) describe a continuous displacement quantum
measurements Jacobs and Steck (2006); Wiseman and Milburn (2010). Similar
equations can be derived for the transverse motion, with
$\Gamma_{\text{qba}}^{(j)}=2\pi g_{j}^{2}j^{2}_{\text{zpf}}$ and $j={x,y}$.
The interacting modes along these direction are symmetric (differently from
the $z$ mode), thus do not feature any constant radiation pressure force.
Because the input optical field is in the vacuum state, the quadratures have
the following non-zero correlations

$\displaystyle\langle
X_{\mathrm{in}}(t)X_{\mathrm{in}}(t^{\prime})\rangle=\langle
Y_{\mathrm{in}}(t)Y_{\mathrm{in}}(t^{\prime})\rangle$
$\displaystyle=\frac{1}{2}\delta(t-t^{\prime}).$ (S39)

This has two profound implications for the measurement.

Firstly, the estimation of the particle position from a phase quadrature
measurement of the output field, Eq. (S36b), is affected by measurement
imprecision due to the fluctuating input phase quadrature. We calibrate this
measurement imprecision in displacement units by dividing by
$\sqrt{2\Gamma_{\text{qba}}/j_{\text{zpf}}}$. From Eqs. (S36b) and (S39), we
find the following imprecision power spectral densities (PSDs) 555Given a
hermitian operator $\mathcal{O}(t)$, we define its spectral density as
$S_{\mathcal{O}\mathcal{O}}(\Omega)=1/(2\pi)\int_{\mathbb{R}}\mathrm{d}\tau\
e^{i\Omega\tau}\langle\overline{\mathcal{O}(t)\mathcal{O}(t+\tau)}\rangle$,
where the bar indicates symmetrization.

$\begin{pmatrix}\overline{S}^{(x)}_{\mathrm{imp}}\\\
\overline{S}^{(y)}_{\mathrm{imp}}\\\
\overline{S}^{(z)}_{\mathrm{imp}}\end{pmatrix}=\frac{5}{8\pi}\frac{1}{k_{0}^{2}}\frac{\hbar\omega_{0}}{P}\begin{pmatrix}1\\\
1/2\\\ \left(2/5+A^{2}\right)^{-1}\end{pmatrix},$ (S40)

where we made again use of the power scattered by the dipole
$P=(\omega_{0}^{4}\alpha^{2}E_{0}^{2})/(12\pi\varepsilon_{0}c^{3})$. These
spectra are inversely proportional to the number of scattered photons,
$P/(\hbar\omega_{0})$, and to the (square of the) wave number $k_{0}$. This
makes sense, since the uncertainty is reduced if the number of probing photons
is increased and since the resolution of each photon increases with lower
wavelengths.

Secondly, the vacuum fluctuations in the amplitude quadrature randomly drive
the mechanical motion, as shown in Eq. (S37b). This fluctuating force, known
as photon recoil Jain  _et al._ (2016); Gonzalez-Ballestero  _et al._ (2019);
Tebbenjohanns  _et al._ (2019b), is a form of quantum backaction, a ubiquitous
aspect of quantum measurements. The photon recoil PSDs are

$\begin{pmatrix}\overline{S}^{(x)}_{FF}\\\ \overline{S}^{(y)}_{FF}\\\
\overline{S}^{(z)}_{FF}\end{pmatrix}=\frac{1}{5}\frac{\hbar^{2}k_{0}^{2}}{2\pi}\frac{P}{\hbar\omega_{0}}\begin{pmatrix}1\\\
2\\\ 2+5A^{2}\end{pmatrix}.$ (S41)

The strength of this backaction increases with the number of the scattered
photons and the momentum carried by each of them.

The imprecision and quantum backaction, respectively Eqs. (S40) and (S41), are
not independent. In fact, their product fulfills the relation

$\overline{S}^{(j)}_{\mathrm{imp}}\,\overline{S}^{(j)}_{FF}=\frac{1}{4\pi^{2}}\frac{\hbar^{2}}{4}.$
(S42)

This is an example of the Heisenberg measurement-disturbance relation
Braginsky and Khalili (1992), and stems from the quantum nature of the
measuring system, that is the electromagnetic field.

####  I.3.2 Coupling to a thermal bath

In practice, the mechanical system is also driven by fluctuating forces other
than the photon recoil. We can take these additional fluctuating forces into
account by including an effective thermal bath. In the absence of any other
forces, the mechanical motion thermalizes to the same average phonon occupancy
of this bath, $\overline{n}_{\mathrm{th}}$, at a rate $\gamma_{m}$. The
equations of motion become then

$\displaystyle\dot{q}$ $\displaystyle=\Omega_{z}p,$ (S43a)
$\displaystyle\dot{p}$
$\displaystyle=-\Omega_{z}q-\gamma_{m}p+\xi(t)+\sqrt{4\Gamma_{\text{qba}}}X_{\text{in}}(t),$
(S43b)

where we have introduced an effective thermal force with correlations
$\langle\overline{\xi(t)\xi(t^{\prime})}\rangle=2\gamma_{m}\left(\overline{n}_{\mathrm{th}}+1/2\right)\delta(t-t^{\prime})$
Giovannetti and Vitali (2001).

Equations (S36) and (S43) correspond to Eqs. (3) and (4) of the main text.

##  II Experimental details

In this section, we provide more details about experimental methods and
calibration, the spectral fittings and the reconstruction of the Wigner
functions.

###  II.1 Systematic error and calibration

We operate the balanced homodyne receiver with a local oscillator (LO) much
stronger than the backscattered beam ($P_{\mathrm{LO}}\sim
1\leavevmode\nobreak\ \mathrm{mW}$, $P_{\mathrm{bs}}\sim 2\leavevmode\nobreak\
\mathrm{\mu W}$). Excess intensity noise in the LO can manifest itself in the
measured spectra as additional imprecision noise. One can cancel this
systematic error by using identical photodiodes and balancing the optical
powers in front of them. In practice, these two conditions are never exactly
met, as photodiodes have always a slightly different responsivity and the
optical powers need to be unbalanced to measure quadratures at angles
$\theta\neq\pi/2$. The systematic error in the spectral imprecision depends on
the power unbalancing, that is, on the angle of the measured optical
quadrature Safavi-Naeini  _et al._ (2013); Chen  _et al._ (2020).

We characterize this error by blocking the backscattered field and unbalancing
the local oscillator. This is done by rotating the half-wave plate preceding
the polarizing beam-splitter in front of the photodetector, shown in Fig.
S1(a).

Figure S1: Homodyne receiver characterization. (a) Sketch of the homodyne
receiver. The LO and the scattered field are linearly polarized along
orthogonal directions. HWP, half-wave plate; PBS, polarizing beamsplitter; LO,
local oscillator. (b) Relative variance (dark blue circle) as a function of
the local oscillator power unbalance. The light blue line is a fit to a
quadratic model. The gray shaded area corresponds to the range of unbalancing
explored in the main experiment.

For a given position of the half-wave plate, which corresponds to a given
value of power unbalancing, we record the PSD of the photodetctor output
voltage. Then, we integrate the PSD from $75$ kHz to $85$ kHz to obtain the
variance $\sigma^{2}$. In Fig. S1(b) we show the results for several values of
the power unbalancing, which is expressed as the difference between the DC
voltages of the photodiodes. The measured variances are expressed relative to
the one at the balanced configuration, $\sigma^{2}_{0}$, for which we know
that the dominant source of noise is the optical shot noise Tebbenjohanns  _et
al._ (2021). We fit the data to a parabola (light blue), which models the
measured noise in the non-ideal case Safavi-Naeini  _et al._ (2013); Chen  _et
al._ (2020). The gray shaded area indicates the region of unbalancing that we
access in the experiments reported in the main text.

In the following we report the calibration procedure employed for the spectra
reported in the main text. For each angle $\theta$, we record the
corresponding PSD and the photodiodes power unbalancing shown as the voltage
value on the horizontal axis of Fig. S1(b). We repeat this measurement for a
set of angles in the range $[0,\pi]$. Then, we block the backscattered field
and we record our reference PSD corresponding to the (balanced) local
oscillator only. We compute its average spectral value around the mechanical
resonance frequency and use this value to calibrate all the other PSDs, i.e.,
to convert them from $V^{2}/\text{Hz}$ to $\text{Hz}/\text{Hz}$ (shot noise
units). Finally, we subtract from each PSD the systematic error in their
background, according to the measured power unbalancing and the
characterization in Fig. S1(b). In our regime of operation
($P_{\mathrm{LO}}\gg P_{\text{bs}}$), the excess classical noise measured by
the photodiode mainly comes from the LO rather then the backscattered beam,
thus it can be safely subtracted.
We notice that this excess noise contributes to at most $\leq 5\%$, much less
that the maximum measured squeezing of $25\%$. Thus, this systematic error
does not call into question the observation of ponderomotive squeezing.

The homodyne angle $\theta$ is inferred from the DC output voltage of the
balanced homodyne detector. This voltage can be modelled as
$V_{\text{DC}}=V_{\text{off}}-V_{\text{amp}}\cos(\theta)$. The knowledge of
the offset voltage and the peak amplitude allows us to invert this relation
and express $\theta$ in terms of the measured $V_{\text{DC}}$. To measure
these two parameters, we sweep the angle $\theta$ by modulating the local
oscillator path length with a ramp waveform and recording the oscillating
output voltage $V_{\text{DC}}$. The maximum and minimum recorded voltages can
be promptly linked to the offset and peak amplitude values.

###  II.2 Fitting the spectra

We record the PSDs at different angles in the range $[0,\pi]$, as shown in
Fig. S2(a). We simultaneously fit all these PDSs to Eq. (5) of the main text.

Figure S2: Homodyne spectra for different quadratures. (a) PSDs of the
homodyne photocurrent at various angle, from $0$ to $\pi$. The spectra are
normalized to the shot noise value. (b) Predicted spectra from the parameters
extracted from the fit to Eq. (5) of the main text.

The free parameters are the mechanical damping rate $\gamma_{m}$, the total
decoherence rate $\Gamma_{\text{tot}}$, the measurement rate
$\Gamma_{\text{meas}}$, the mechanical resonance frequency $\Omega_{m}$ and
the angle $\theta$. We constrain the former three parameters to be the same
for all the PSDs, whereas the latter two can vary across different PSDs. This
is to take into account slow drifts of the mechanical resonance frequency
across the duration of the experiment ($\sim 1\leavevmode\nobreak\
\mathrm{h}$) as well as errors in the calibration of the angle from the DC
value of the detector output voltage (termed ‘inferred angle’ from now on).

Figure S3: Results of the fits. (a) Fitted homodyne angle and (b) fitted
mechanical resonance frequency. Each blue circle corresponds to a measured
spectrum, labelled by the inferred angle $\theta$. The light blue lines in (a)
and (b) show, respectively, the expected linear and constant trend.

In Fig. S3(a) and (b) we show the fitted angles and mechanical resonance
frequencies for the different spectra, labelled by the inferred angle
$\theta$. We observe that the fitted angles are linearly proportional to the
inferred angles, as shown by the light blue line. The small offset of
$0.05\pi$ is due to a weak field reflected by the trapping lens and co-
propagating with the backscattered field, and is consistent with what was
observed in Fig. 1(c) of the main text.

The mechanical resonance frequency has a mean value of
$\Omega_{m}/(2\pi)=73.24$ kHz and a relative standard deviation of $0.05\%$.
We hypothesize that this drift arises from slow fluctuations of the power of
the trapping laser. We use the mean value of the resonance frequency and the
fitted damping rate, measurement rate and total decoherence rate to plot the
predicted spectra, shown in Fig. S2(b).

###  II.3 Wigner quasidistribution reconstruction

To reconstruct the distributions in Fig. 3 of the main text, we record a time
trace of the nanoparticle for values of the analyser angle $\theta$ spanning
the interval $[0,\pi]$. We split the time traces in chunks of
$\approx$8\text{\,}\mathrm{m}\mathrm{s}$$ Tebbenjohanns  _et al._ (2021),
apply a Hann window and perform a Fourier transform of each chunk. The real
and imaginary parts of the Fourier transform at frequency $f$ are treated
equally, in other words they are independent measurements of $X_{\theta}(f)$.
For every frequency of interest $f$, we evaluate the Fourier transform of each
chunk at $f$, concatenate real and imaginary parts and perform a histogram of
the data. The histograms of the quadrature measurements as a function of the
analyser angle $\theta$ constitute a so called sinogram, shown in Fig. S4 for
$f=$70.1\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$.

Figure S4: Sinogram of the optical mode at frequency
$f=$70.1\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$. Each column corresponds to
a histogram of the Fourier components at frequency $f$. The dashed gray (dash-
dotted orange) line corresponds to the vacuum fluctuations (variance of the
histogram).

We then perform a numerical inverse Radon transform of the sinogram based on
simultaneous algebraic reconstruction Andersen and Kak (1984); Van Der Walt
_et al._ (2014). The covariance circle of the vacuum fluctuations, which is
also used to calibrate the axes of the graphs, is estimated by blocking the
particle signal, recording the shot noise of the local oscillator and
estimating the variance of the signal. To estimate the covariance ellipse of
the particle signal, instead, we extract three vertical cuts of the sinogram,
respectively at angles $\theta_{1}=0$, $\theta_{2}=\pi/4$ and
$\theta_{3}=\pi/2$. We fit these cuts with a gaussian distribution and extract
the three respective variances $V_{1},V_{2},V_{3}$. The covariance matrix of
the gaussian Wigner distribution follows from $\langle X^{2}\rangle=V_{1}$,
$\langle Y^{2}\rangle=V_{3}$ and $\langle XY\rangle=(V_{1}+V_{3})/2-V_{2}$.

## References

  * Chang  _et al._ (2010) D. E. Chang, C. A. Regal, S. B. Papp, D. J. Wilson, J. Ye, O. Painter, H. J. Kimble, and P. Zoller, Proc. Natl. Acad. Sci. U.S.A. 107, 1005 (2010).
  * Romero-Isart  _et al._ (2011) O. Romero-Isart, A. C. Pflanzer, M. L. Juan, R. Quidant, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, Phys. Rev. A 83, 013803 (2011).
  * Rodenburg  _et al._ (2016) B. Rodenburg, L. P. Neukirch, A. N. Vamivakas, and M. Bhattacharya, Optica 3, 318 (2016).
  * Gonzalez-Ballestero  _et al._ (2019) C. Gonzalez-Ballestero, P. Maurer, D. Windey, L. Novotny, R. Reimann, and O. Romero-Isart, Phys. Rev. A 100, 013805 (2019).
  * Magrini  _et al._ (2021) L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Nature 595, 373 (2021).
  * Maurer  _et al._ (2021) P. Maurer, C. Gonzalez-Ballestero, and O. Romero-Isart, “Quantum electrodynamics with a nonmoving dielectric sphere: Quantizing lorenz-mie scattering,” (2021), arXiv:2106.07975 [quant-ph] .
  * Ricci  _et al._ (2017) F. Ricci, R. A. Rica, M. Spasenović, J. Gieseler, L. Rondin, L. Novotny, and R. Quidant, Nat Commun 8, 15141 (2017).
  * Jain  _et al._ (2016) V. Jain, J. Gieseler, C. Moritz, C. Dellago, R. Quidant, and L. Novotny, Phys. Rev. Lett. 116, 243601 (2016).
  * Tebbenjohanns  _et al._ (2019a) F. Tebbenjohanns, M. Frimmer, A. Militaru, V. Jain, and L. Novotny, Phys. Rev. Lett. 122, 223601 (2019a).
  * Tebbenjohanns  _et al._ (2020) F. Tebbenjohanns, M. Frimmer, V. Jain, D. Windey, and L. Novotny, Phys. Rev. Lett. 124, 013603 (2020).
  * Hebestreit  _et al._ (2018) E. Hebestreit, M. Frimmer, R. Reimann, and L. Novotny, Phys. Rev. Lett. 121, 063602 (2018).
  * Hoang  _et al._ (2018) T. M. Hoang, R. Pan, J. Ahn, J. Bang, H. T. Quan, and T. Li, Phys. Rev. Lett. 120, 080602 (2018).
  * Novotny (2017) L. Novotny, Phys. Rev. A 96, 032108 (2017).
  * Pflanzer  _et al._ (2012) A. C. Pflanzer, O. Romero-Isart, and J. I. Cirac, Physical Review A - Atomic, Molecular, and Optical Physics (2012), 10.1103/PhysRevA.86.013802.
  * Note (1) We define $a^{\dagger}_{j}(\omega)=\left[a_{j}(\omega)\right]^{\dagger}$.
  * Gardiner and Collett (1985) C. W. Gardiner and M. J. Collett, Phys. Rev. A 31, 3761 (1985).
  * Clerk  _et al._ (2010) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010).
  * Aspelmeyer  _et al._ (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * Note (2) The integral over the total solid angle takes the standard definition $\intop\ilimits@_{4\pi}\mathrm{d}\Omega=\intop\ilimits@_{0}^{\pi}\mathrm{d}\theta\ \intop\ilimits@_{0}^{2\pi}\mathrm{d}\varphi$.
  * Note (3) Use $[a_{\varepsilon},a_{j}^{\dagger n}]=\Lambda_{j}(\varepsilon,\mathbf{k})na_{j}^{\dagger(n-1)}$ and the Taylor expansion of $\mathop{exp}\nolimits(\beta a_{j}^{\dagger})$.
  * Tebbenjohanns  _et al._ (2019b) F. Tebbenjohanns, M. Frimmer, and L. Novotny, Phys. Rev. A 100, 043821 (2019b).
  * Note (4) Note that we use the expression interacting mode to refer to the spatial shape of the optical modes described by Eqs. (S12). Each interacting mode consists then of infinitely many frequency modes.
  * Jacobs and Steck (2006) K. Jacobs and D. A. Steck, Contemp. Phys. 47, 279 (2006).
  * Wiseman and Milburn (2010) H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_ (Cambridge University Press, 2010).
  * Note (5) Given a hermitian operator $\mathcal{O}(t)$, we define its spectral density as $S_{\mathcal{O}\mathcal{O}}(\Omega)=1/(2\pi)\intop\ilimits@_{\mathbb{R}}\mathrm{d}\tau\ e^{i\Omega\tau}\langle\overline{\mathcal{O}(t)\mathcal{O}(t+\tau)}\rangle$, where the bar indicates symmetrization.
  * Braginsky and Khalili (1992) V. B. Braginsky and F. Y. Khalili,  _Quantum Measurement_ (Cambridge University Press, Cambridge, 1992).
  * Giovannetti and Vitali (2001) V. Giovannetti and D. Vitali, Phys. Rev. A 63, 023812 (2001).
  * Safavi-Naeini  _et al._ (2013) A. H. Safavi-Naeini, S. Gröblacher, J. T. Hill, J. Chan, M. Aspelmeyer, and O. Painter, Nature 500, 185 (2013).
  * Chen  _et al._ (2020) J. Chen, M. Rossi, D. Mason, and A. Schliesser, Nat. Commun. 11, 943 (2020).
  * Tebbenjohanns  _et al._ (2021) F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Nature 595, 378 (2021).
  * Andersen and Kak (1984) A. H. Andersen and A. C. Kak, Ultrasonic Imaging (1984), 10.1177/016173468400600107.
  * Van Der Walt  _et al._ (2014) S. Van Der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager, E. Gouillart, and T. Yu, PeerJ (2014), 10.7717/peerj.453.