# Macroscopic quantum test with bulk acoustic wave resonators

Björn Schrinski  Center for Hybrid Quantum Networks (Hy-Q), Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17, DK-2100 Copenhagen,
Denmark  Yu Yang  Department of Physics, ETH Zürich, 8093 Zürich, Switzerland
Uwe von Lüpke  Department of Physics, ETH Zürich, 8093 Zürich, Switzerland
Marius Bild  Department of Physics, ETH Zürich, 8093 Zürich, Switzerland
Yiwen Chu  Department of Physics, ETH Zürich, 8093 Zürich, Switzerland  Klaus
Hornberger  University of Duisburg-Essen, Faculty of Physics, Lotharstraße 1,
47048 Duisburg, Germany  Stefan Nimmrichter  Naturwissenschaftlich-Technische
Fakultät, Universität Siegen, Walter-Flex-Straße 3, 57068 Siegen, Germany
Matteo Fadel<EMAIL_ADDRESS>Department of Physics, ETH Zürich, 8093
Zürich, Switzerland

###### Abstract

Recently, solid-state mechanical resonators have become a platform for
demonstrating non-classical behavior of systems involving a truly macroscopic
number of particles. Here, we perform the most macroscopic quantum test in a
mechanical resonator to date, which probes the validity of quantum mechanics
by ruling out a classical description at the microgram mass scale. This is
done by a direct measurement of the Wigner function of a high-overtone bulk
acoustic wave resonator mode, monitoring the gradual decay of negativities
over tens of microseconds. While the obtained macroscopicity of $\mu=11.3$ is
on par with state-of-the-art atom interferometers, future improvements of mode
geometry and coherence times could test the quantum superposition principle at
unprecedented scales and also place more stringent bounds on spontaneous
collapse models.

Understanding the quantum-classical transition is one of the main challenges
of modern physics. Is the Schrödinger equation valid all the way from the
microscopic to the macroscopic world, with quantum effects increasingly hard
to observe due to environmental decoherence [1]? Or do the laws of quantum
mechanics break down at some point, so as to reinstate macrorealism in our
everyday life [2]? This question is not philosophical, as it can be tackled by
demonstrating genuine quantum effects in ever more macroscopic systems.

Single atoms have been delocalized over macroscopic length and time scales,
reaching metres and seconds in state-of-the-art experiments [3, 4, 5, 6]. On
the other hand, molecule interferometry has pushed the mass boundaries by
confirming the wave nature for compounds of more than $10^{4}$ atomic mass
units [7, 8]. Future Earth- and space-based experiments with freely falling or
levitated nanoparticles aim at even higher mass and time scales [9, 10], but
in such schemes it is a single rigid-body degree of freedom that is placed
into a superposition.

This is to be contrasted with growing efforts in electro- and opto-mechanics
toward demonstrating quantum states in the elastic deformation of condensed
matter systems of substantially higher mass. Recent experimental milestones
include the successful preparation of nonclassical states of surface and bulk
acoustic modes, demonstrated through the observation of Wigner function
negativities [11, 12, 13, 14, 15]. Such nonclassical states of motion involve
a truly macroscopic number of atoms, on the order of $10^{16}$, vibrating in
unison. At the same time, each atom is delocalized only over a tiny fraction
of the size of an atomic nucleus, which makes unclear how they compare to
previous tests of the quantum-classical transition.

Using an objective measure for the degree of macroscopicity of nonclassical
states [16], we report here the most macroscopic quantum test performed with a
solid-state mechanical resonator. The macroscopicity value obtained in our
experiment is surprisingly large, motivating us to optimize the design
specifications for future electromechanical experiments to become competitive
even with matter-wave superposition tests.

Figure 1: Illustration of the HBAR device used in this work (not to scale, see
[14] for details). A superconducting qubit couples to a mechanical acoustic
mode localized in the bulk of a sapphire crystal through a piezoelectric
material. This allows preparing and measuring the quantum state of the phonon
mode.

To justify the macroscopicity measure used, note that any genuine quantum test
probes the validity of quantum mechanics against the hypothesis of
macrorealism: that a modification of quantum theory destroys coherence above a
certain size or mass scale, rendering superposition phenomena unobservable in
practice. Hence, the degree to which a quantum experiment falsifies a broad,
generic class of minimally invasive, macrorealistic modifications of quantum
mechanics gives rise to an objective benchmark [16]: Using the measurement
data to rule out a range of values of the coherence time parameter $\tau_{e}$
quantifying the macrorealist modification, an experiment is justifiably the
more macroscopic the greater the $\tau_{e}$-values falsified through Bayesian
inference [17]. The macroscopicity measure $\mu$ is then determined by the
logarithm of the greatest ruled out $\tau_{e}$-value 111Note that the
macrorealist modification hypothesis can be probed also in experiments that to
do not involve quantum superpositons. Only observations of genuine quantum
signatures can therefore be assigned a benchmark value $\mu$ [17]. This
requires the interrogated state not to have an underlying classical
description, excluding all oscillator states with positive Wigner function..
To date, the highest macroscopicity $\mu=14$ has been achieved in molecule
interferometry [8].

In our experiment, we observe Wigner function negativities for up to
$40\,\mu$s in an bulk acoustic vibration mode of a sapphire crystal with an
effective mass of $1\,\mu$g, which is prepared in a single-phonon state and in
a superposition of Fock states. This quantum signature is detected by
measuring the displaced parity [19] of the quantum state using a
superconducting qubit, see Fig. 1. Our analysis shows that it can be
associated with a large macroscopicity $\mu=11.3$, despite the sub-nuclear
delocalization of the atoms involved. Based on this, we explain how future
experiments with bulk acoustic resonators may realize the most macroscopic
quantum test.

Figure 2: Measured Wigner functions of a bulk acoustic mode initially prepared
in the Fock state $\ket{1}$ (top row) and in the superposition state
$(\ket{0}+\ket{1})/\sqrt{2}$ (bottom row). Both states are initially non-
classical, as indicated by the phase space regions with negative quasi-
probability (blue). As time evolves (left to right), both states approach the
ground state $\ket{0}$ due to relaxation. Comparing the measured data to the
theoretical model Eq. (Macroscopic quantum test with bulk acoustic wave
resonators) allows us to extract the maximal compatible diffusion rate
$\Gamma$ through Bayesian inference, yielding the macroscopicity $\mu$ for the
considered experiment.

_Macroscopicity of bulk resonator modes._ – To evaluate the macroscopicity of
a demonstrated quantum effect, one should consider the class of macrorealistic
models possibly falsified by the experiment. Keeping those models minimally
invasive in that they preserve linearity, the semigroup nature of the time
evolution, invariance under the chosen non-relativistic reference frame, and
other basic consistency principles [16], their observable impact can be
formulated as a Lindblad generator $\mathcal{L}\rho$ modifying the quantum
evolution of the system state $\rho$. For a large compound of $N\gg 1$ atomic
constituents with masses $m_{1},\ldots,m_{N}$, the generator reduces to [16]

$\mathcal{L}\rho=\frac{1}{\tau_{e}}\int\mathrm{d}^{3}q\,g(q,\sigma_{q})\mathcal{D}[\mathsf{L}(\bm{q})]\rho\;,$
(1)

with
$\mathcal{D}[\mathsf{L}]\rho=\mathsf{L}\rho\mathsf{L}^{\dagger}-\\{\mathsf{L}^{\dagger}\mathsf{L},\rho\\}/2$.
The strength of the modification is specified by the coherence time parameter
$\tau_{e}$ and the standard deviation $\sigma_{q}$ of the isotropic Gaussian
momentum distribution
$g(q,\sigma_{q})=e^{-q^{2}/2\sigma_{q}^{2}}/(2\pi\sigma_{q}^{2})^{3/2}$. The
Lindblad operators are mass-weighted sums of single-particle momentum kicks by
$\bm{q}$,
$\mathsf{L}(\bm{q})=\sum_{n=1}^{N}(m_{n}/m_{e})e^{i\bm{q}\cdot\bm{\mathsf{r}}_{n}/\hbar}$,
with $\bm{\mathsf{r}}_{n}$ the particle position operators and $m_{e}$ a
reference mass chosen to be that of an electron.

We note that the theory of continuous spontaneous collapse (CSL) can be
described by (1), but it is usually parameterized in terms of a rate per
atomic mass unit $\lambda_{\mathrm{CSL}}=(1\,\mathrm{u}/m_{e})^{2}/\tau_{e}$,
and a localisation length $r_{\mathrm{CSL}}=\hbar/\sqrt{2}\sigma_{q}$ [20].
The strongest bounds on these CSL parameters are currently achieved by so-
called non-interferometric tests [21], where one measures the noise
experienced by a classical system by exploiting the fact that macrorealistic
modifications are accompanied by a heating effect. However, a genuine quantum
signal, such as interference fringes or a negative Wigner function, is
required for a CSL test to also probe the validity of quantum mechanics in the
system.

In the following, we consider a collective acoustic excitation in the bulk
condensed matter, taking the form

$\bm{u}(\bm{r})=\exp\left(-\frac{y^{2}+z^{2}}{w_{0}^{2}}\right)\cos\left(\pi\frac{\ell
x}{L}\right)\mathbf{e}_{x}\;$ (2)

of waist $w_{0}$, mode index $\ell$, and length $L$. (The mode vanishes for
$x\notin\left[0,L\right]$, and is laterally well confined within the
material.) The effective volume of the mode is defined as
$V_{\mathrm{eff}}=\int\mathrm{d}^{3}r\,u^{2}(\bm{r})=\pi w_{0}^{2}L/4$ and,
given the mass density $\varrho$, the effective mass is defined as
$m_{\mathrm{eff}}=\int\mathrm{d}^{3}r\,\varrho(\bm{r})u^{2}(\bm{r})$. This way
one consistently obtains the total bulk volume and mass in the limiting case
of a center of mass motion, corresponding to a uniform displacement with
$u(\bm{r})=1$, so that $m_{\mathrm{eff}}$ can be understood as the mass
fraction taking part in the displacement [22].

The mode is quantized by means of the ladder operator
$\mathsf{a}=(\mathsf{X}+i\mathsf{P})/\sqrt{2}$, involving a conjugate pair of
dimensionless position and momentum quadrature operators,
$[\mathsf{X},\mathsf{P}]=i$. A single phononic excitation thus displaces the
atoms by amplitudes of the order of the zero-point fluctuation
$x_{0}/\sqrt{2}$, with $x_{0}=\sqrt{\hbar/m_{\mathrm{eff}}\omega}$, around
their equilibrium positions $\bm{r}_{n}$ according to the mode displacement
field, i.e.
$\bm{\mathsf{r}}_{n}=\bm{r}_{n}+\bm{u}(\bm{r}_{n})x_{0}\mathsf{X}$. The
Lindblad operators thus become

$\mathsf{L}(\bm{q})=\frac{1}{m_{e}}\int\mathrm{d}^{3}r\,\varrho(\bm{r})e^{i\bm{q}\cdot\bm{r}/\hbar}e^{i\bm{q}\cdot\bm{u}(\bm{r})x_{0}\mathsf{X}/\hbar}\;.$
(3)

In the relevant regime of modification length scales $\hbar/\sigma_{q}$ larger
than the interatomic distances, it is safe to assume a homogeneous mass
density $\varrho(\bm{r})=\mathrm{bar}{\varrho}$ and to expand the operator in
(3) to first order in $\mathsf{X}$ [17]. This results in momentum diffusion,
$\mathcal{L}\rho=2\Gamma\mathcal{D}[\mathsf{X}]\rho$, with a diffusion rate
amplified by the effective oscillator mass,

$\displaystyle\Gamma$
$\displaystyle=\frac{\mathrm{bar}{\varrho}^{2}x_{0}^{2}}{2m_{e}^{2}\tau_{e}\hbar^{2}}\int\mathrm{d}^{3}q\,g(q,\sigma_{q})\,q_{x}^{2}\left|\int\mathrm{d}^{3}r\,u_{x}(\bm{r})e^{i\bm{q}\cdot\bm{r}/\hbar}\right|^{2}$
$\displaystyle=\frac{m_{\mathrm{eff}}^{2}}{m_{e}^{2}\tau_{e}}\frac{(4x_{0}/L)^{2}}{1+\sigma_{w}^{2}}\int\mathrm{d}\zeta\,\frac{e^{-\zeta^{2}/2\sigma_{L}^{2}}}{\sqrt{2\pi}\sigma_{L}}\frac{1-(-)^{\ell}\cos\zeta}{\left(1-\pi^{2}\ell^{2}/\zeta^{2}\right)^{2}}\;.$
(4)

The rate depends on how the modification length scale $\hbar/\sigma_{q}$
compares to the geometric length scales of the oscillator mode, as determined
by $\sigma_{w}=w_{0}\sigma_{q}/\hbar$, $\sigma_{L}=L\sigma_{q}/\hbar$,
$x_{0}$, and $\ell$. The integral can be given analytically, leading to a
lengthy expression for $\Gamma$. In the limit $\sigma_{L},\sigma_{w}\ll 1$ and
for $\ell$ even, we find $\Gamma\propto
m_{\mathrm{eff}}\sigma_{q}^{2}\sigma_{L}^{4}/\ell^{4}$, whereas for
$\sigma_{L}\gg\pi\ell$ and $\sigma_{w}\gg 1$, we get $\Gamma\propto
m_{\mathrm{eff}}/\sigma_{w}^{2}$. The maximum $\Gamma$ with respect to
$\sigma_{q}$ is located in between at
$\hbar/\sigma_{q}\approx\sqrt{3}L/\pi\ell$, provided that $\pi\ell\gg 1$ and
$\pi\ell w_{0}/\sqrt{3}L\gg 1$ [22],

$\max_{\sigma_{q}}\Gamma\approx\sqrt{\frac{3\pi}{2e^{3}}}\frac{6\hbar\mathrm{bar}{\varrho}}{m_{e}^{2}\omega\tau_{e}}\frac{L}{\ell}\;.$
(5)

Given a fixed sound velocity $v$, and noting $\omega=2\pi v\ell/L$, the
maximum rate value is proportional to the square of the mode wavelength
$L/\ell$.

In the co-rotating frame of the oscillator, and including energy decay at rate
$\gamma_{\downarrow}$, the modified time evolution of the oscillator state
$\rho$ is given by the master equation

$\displaystyle\dot{\rho}$
$\displaystyle=\Gamma\mathcal{D}\left[\mathsf{a}e^{-i\omega
t}+\mathsf{a}^{\dagger}e^{i\omega
t}\right]\rho+\gamma_{\downarrow}\mathcal{D}[\mathsf{a}]\rho$
$\displaystyle\approx\left(\Gamma+\gamma_{\downarrow}\right)\mathcal{D}[\mathsf{a}]\rho+\Gamma\mathcal{D}[\mathsf{a}^{\dagger}]\rho\;.$
(6)

In the second line we average over the rapidly oscillating cross-terms, as we
are interested in decoherence on time scales much longer than the oscillation
period.

The coarse-grained isotropic master equation (6) admits a compact solution in
the Wigner function representation

$\displaystyle W(X,P;t)=\frac{e^{\gamma_{\downarrow}t}}{\pi
S(t)}\int\mathrm{d}X^{\prime}\mathrm{d}P^{\prime}$ $\displaystyle
W\left(X^{\prime}e^{\gamma_{\downarrow}t/2},P^{\prime}e^{\gamma_{\downarrow}t/2};0\right)$
$\displaystyle\times e^{-[(X-X^{\prime})^{2}+(P-P^{\prime})^{2}]/S(t)}\;,$ (7)

with $S(t)=(1+2\Gamma/\gamma_{\downarrow})(1-e^{-\gamma_{\downarrow}t})$. This
expression can be computed analytically for several states of interest, such
as Fock states and their superposition [22]. Bayesian parameter estimation
based upon Eq. (Macroscopic quantum test with bulk acoustic wave resonators)
allows us to identify from the measured Wigner functions the threshold value
$\Gamma_{5\%}$ that corresponds to the most conservative $5\,\%$ quantile of
the coherence time $\tau_{e}$ at any $\sigma_{q}$ [22, 17]. Greater $\tau_{e}$
are also compatible with the data, as the observed decoherence might as well
be due to unspecified conventional noise sources, but smaller times can be
ruled out with confidence. Macroscopicity is then defined as follows: We
convert $\Gamma_{5\%}$ into the corresponding macrorealistic coherence time
parameter $\tau_{e}=\tau_{e}(\sigma_{q},\Gamma_{5\%})$ by virtue of (4),
maximize over $\sigma_{q}$ to obtain the greatest excluded $\tau_{e}$, and
assign a macroscopicity value $\mu=\log_{10}(\tau_{e}/1\,\mathrm{s})$.

By modelling the state evolution with (6), we thus assume a zero-temperature
environment and attribute any heating or decoherence beyond energy relaxation
at the rate $\gamma_{\downarrow}$ to the impact of a macrorealistic
modification described by infinite-temperature diffusion at the rate $\Gamma$.
This amounts to a _conservative_ assessment of the macroscopicity, since, by
overestimating the hypothetical contribution of macrorealistic diffusion to
the observed decoherence, we overestimate $\Gamma$ and thus underestimate
$\tau_{e}$. We characterise the energy relaxation rate $\gamma_{\downarrow}$
independently, by a standard measurement of the phonon decay time $T_{1}$ and
the steady-state population [22].

_Experiment._ – Our measurements are carried out on a high-overtone bulk-
acoustic wave resonator (HBAR) device cooled down at millikelvin temperature,
see Fig. 1, which makes use of a superconducting qubit to prepare, control,
and read-out the quantum excitation of a mechanical mode localized in the bulk
of a sapphire substrate [14].

The acoustic vibration mode is described to a very good approximation by Eq.
(2), as also verified by numerical simulations. It oscillates over the bulk
length $L=$435\text{\,}\mu\mathrm{m}$$ at angular frequency $\omega=$2\pi\cdot
5.961\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ with wavelength
$\lambda=$1.8\text{\,}\mu\mathrm{m}$$, corresponding to the mode index
$\ell=2L/\lambda=486$. Assuming a homogenous mass density
$\mathrm{bar}{\varrho}=3.98\,$g/cm3, and given the transverse waist
$w_{0}=$27\text{\,}\mu\mathrm{m}$$, we obtain the effective oscillator mass
$m_{\mathrm{eff}}=\pi w_{0}^{2}L\mathrm{bar}{\varrho}/4=1.0\,\mu$g and a zero-
point fluctuation of $x_{0}/\sqrt{2}=$1.2\cdot 10^{-18}\text{\,}\mathrm{m}$$.
For these parameters, we find from Eq. (4) a maximal diffusion rate of
$\Gamma_{\mathrm{max}}=3.5\cdot 10^{13}/\tau_{e}$ at a critical length scale
$\hbar/\sigma_{q}\simeq 0.5\,\mu$m, which is of the same order as the mode
wavelength.

The coherence properties of the phonon mode are characterized by a small one-
phonon steady-state population of $1.6\pm 0.2\,\%$, by the relaxation time
$T_{1}=85.8{\pm 1.5}\,\mu$s, and by the Ramsey dephasing time $T_{2}=147.3{\pm
2.6}\,\mu$s [22]. These imply a much longer pure dephasing time of
$T_{\phi}=(1/T_{2}-1/2T_{1})^{-1}=$1.0{\pm
0.2}\text{\,}\mathrm{m}\mathrm{s}$$. Relaxation, which originates from the
coupling to the environment through e.g. surface scattering or diffraction
loss, thus dominates the degradation of quantum features over time, and we can
safely assume $\gamma_{\downarrow}\approx 1/T_{1}$ within error tolerance.

Using the toolbox of circuit quantum acoustodynamics (cQAD), we can use the
qubit to prepare and measure the HBAR in nonclassical states of motion [11].
In particular, an arbitrary superposition with the first excited state,
$a\ket{0}+b\ket{1}$, can be prepared by swapping the qubit state
$a\ket{\downarrow}+b\ket{\uparrow}$ to the mechanical mode through the
resonant Jaynes–Cummings interaction. For any $b\neq 0$ we obtain a
nonclassical state of motion with negativities in the Wigner function, that
will degrade over time due to relaxation. To monitor this effect, we let the
system evolve for a time $t$, and then obtain the Wigner function of the
mechanical state by measuring its displaced parity [14, 22]. Fig. (2) shows
how the measured Wigner function degrades over time for an initial Fock state
$\ket{1}$ (top row) and an initial superposition state
$(\ket{0}+\ket{1})/\sqrt{2}$ (bottom row). Note the disappearance of negative
(blue) regions in the quasi-probability distribution, turning it into a
classical phase-space distribution.

In the first part of the experiment, we prepare the resonator in the single-
phonon Fock state $\ket{1}$ and measure its Wigner quasi-probability
distribution as a function of time (Fig. 2, top row). We observe that the
region of negative values around the phase space origin shrinks gradually,
persisting for tens of microseconds, before the Wigner function turns purely
positive and becomes indistinguishable from a classical distribution. The
longer and more pronounced the initial negativity remains, the stronger the
experiment falsifies macrorealist modifications of quantum mechanics. We
estimate the impact of such modifications by first fitting the initial
measurement at $t=0$ with the Wigner function of an incoherent mixture of
$\ket{0}$ and $\ket{1}$, thus accounting for imperfections in the state
preparation and measurement. The measured Wigner functions at $t=$ 10 $\mu$s,
20 $\mu$s, and 40 $\mu$s are then compared to the theoretical expectation
(Macroscopic quantum test with bulk acoustic wave resonators) including energy
relaxation with the measured rate ${\gamma_{\downarrow}\approx}1/T_{1}$ as
well as modification-induced diffusion with variable rate $\Gamma$. Bayesian
parameter estimation and subsequent maximization with respect to $\sigma_{q}$
shows that any $\Gamma>\Gamma_{5\%}=$1.6\cdot 10^{2}\text{\,}\mathrm{s}^{-1}$$
can be excluded with 95$\%$ confidence (at
$\hbar/\sigma_{q}\simeq$0.5\text{\,}\mu\mathrm{m}$$) [22]. This corresponds to
a macroscopicity of $\mu=11.3$, by far the highest value reported on
mechanical resonators, see Tab. I. We note that the obtained diffusion rate
$\Gamma_{5\%}=$1.6\cdot 10^{2}\text{\,}\mathrm{s}^{-1}$$ is consistent with
what is found by using our system to perform a so-called non-interferometric
test of macrorealistic models based on the observed residual excitation of the
mode [22]. Thus our analytic results can be straightforwardly adopted in such
schemes, as recently proposed in Ref. [23].

In the second part of the experiment, we prepare the resonator in the
superposition state $(\ket{0}+\ket{1})/\sqrt{2}$. As shown in the bottom row
of Fig. 2, the measured Wigner function again takes negative values,
confirming the nonclassicality of the initial state. The asymmetry around the
origin, whose phase can be controlled by the qubit state, rotates due to a
small detuning between the oscillator and the displacement pulse frequencies.
This rotation is accounted for in our data processing. Following the previous
analysis, we find that the superposition state data is compatible with a
maximal diffusion rate of $\Gamma_{5\%}=$6.4\cdot
10^{2}\text{\,}\mathrm{H}\mathrm{z}$$, which is a factor of four larger than
the rate found from the Fock state $\ket{1}$ and from the non-interferometric
test [22]. This discrepancy can be explained by the presence of an additional
pure dephasing channel with a rate of
$1/T_{\phi}\sim$10^{3}\text{\,}\mathrm{H}\mathrm{z}$$. While circularly
symmetric Wigner functions such as single Fock states are unaffected by this
decoherence channel, which corresponds to a random rotation around the phase-
space origin, a superposition of Fock states, which is necessarily asymmetric,
gets washed out at rate $1/T_{\phi}$, removing any negativity. However, to
fairly assess the macroscopicity of the Fock state superposition, we do _not_
include this decoherence channel into the master equation (6), so that all
decoherence is attributed to the macrorealistic modification. This results in
a greater diffusion rate, and consequently in a lower macroscopicity of
$\mu=10.7$. If we included the $1/T_{\phi}$ dephasing effect in the Bayesian
inference the resulting rate would be consistent with both the Fock state
$\ket{1}$ analysis and the non-interferometric test.

_Beating the macroscopicity record._ –Based on our experimental results, and
using the analytical result (4), we can now establish how to optimize future
quantum experiments with bulk acoustic resonators to become competitive with
the most macroscopic quantum tests to date. In contrast to quantum tests based
on the interference of center-of-mass degrees of freedom, the diffusion rate
scales here only weakly with the mass and length scales of the material: in
the realistic case of comparable lateral and longitudinal mode extension
($w_{0}\approx L$), we find that macrorealistic diffusion is strongest at a
length scale $\hbar/\sigma_{q}$ comparable to the mode wavelength $L/\ell$. To
attain greater macroscopicity in a given material, one should strive to
increase the mode wavelength since the maximum diffusion rate (5) grows like
$(L/\ell)^{2}$. Apart from these geometric considerations, the macroscopicity
also benefits from longer relaxation and coherence times $T_{1},T_{2}$, and a
better characterization of the measurement noise. The latter can be
incorporated into the Bayesian parameter estimation to extract better bounds,
instead of adding a broad Gaussian noise channel as done in the analysis
above.

As an example, keeping the same device and mode geometry as in the present
experiment, but using a mode with
$\omega=2\pi\cdot$2\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ and
$T_{1}=$10\text{\,}\mathrm{m}\mathrm{s}$$, would result in $\mu=14.4$,
surpassing all reported quantum tests, see Tab. I. We obtain this estimate by
setting $\ell=160$ and rescaling the observed coherence time in proportion to
$T_{1}$; the macroscopicity grows logarithmically with increasing coherence
time. Beyond that, out-of-plane drum modes with frequencies in the low MHz
range and huge quality factors have been reported, corresponding to
$T_{1}>$100\text{\,}\mathrm{m}\mathrm{s}$$ [24]. If genuine quantum effects
could be demonstrated in these devices, they might reach substantially greater
macroscopicities $\mu>17$ 222This rough estimate is based on [17], by
approximating the mode in [24] as the fundamental Bessel mode
$\bm{u}(\bm{r})=J_{0}(4.81\sqrt{y^{2}+z^{2}}/d)\mathbf{e}_{x}$, for a disc
with diameter $d=100\,\mu$m, $14\,$nm thickness, and $m_{\mathrm{eff}}\simeq
68\,$pg. Precise $\mu$-values would depend on a combination of effective mass,
mode profile, and observed duration of Wigner negativities..

| experiment | year | $\mu$
---|---|---|---
mechanical resonators | Bulk acoustic waves [this work] | 2022 | 11.3
Phononic crystal resonator [13] | 2022 | $\sim 9.0^{*}$
Surface acoustic waves [12] | 2018 | $\sim 8.6^{*}$
matter-wave interference | Molecule interferometry [8] | 2019 | 14.0
Atom interferometry [6] | 2019 | 11.8
BEC interferometry [5] | 2017 | 12.4
Table 1: Macroscopicities of solid-state mechanical resonator experiments
reporting Wigner function negativities ($*$: estimated), and most macroscopic
matter-wave interference experiments assessed in [17, 26]. The values with
asterix are estimates of what could have been observed based on the stated
$T_{1}$ time [22].

_Conclusion._ –We demonstrated the most macroscopic quantum test with bulk
acoustic wave resonators, based on monitoring the time evolution of
nonclassical Wigner functions. The reported macroscopicity is comparable to
the one obtained by an experiment holding cesium atoms in a spatial
superposition at $4\text{\,}\mu\mathrm{m}$ separation over 20 seconds [6]. In
the future, we envisage improved resonator designs that may surpass the most
macroscopic matter wave interferometers, testing the validity of quantum
mechanics at unprecedented scales.

_Acknowledgements._ – We would like to thank Georg Enzian and Eric Planz for
helpful discussions. BS was supported by Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation), Grant No. 449674892. UvL is funded by the SNF
project number 200021_204073. MF was supported by The Branco Weiss Fellowship
– Society in Science, administered by the ETH Zürich.

## References

  * Schlosshauer [2019] M. Schlosshauer, Quantum decoherence, Phys. Rep. 831, 1 (2019).
  * Leggett [2002] A. J. Leggett, Testing the limits of quantum mechanics: motivation, state of play, prospects, J. Phys.: Condens. Matter 14, R415 (2002).
  * Kovachy  _et al._ [2015] T. Kovachy, P. Asenbaum, C. Overstreet, C. A. Donnelly, S. M. Dickerson, A. Sugarbaker, J. M. Hogan, and M. A. Kasevich, Quantum superposition at the half-metre scale, Nature 528, 530 (2015).
  * Abend  _et al._ [2016] S. Abend, M. Gebbe, M. Gersemann, H. Ahlers, H. Müntinga, E. Giese, N. Gaaloul, C. Schubert, C. Lämmerzahl, W. Ertmer, W. P. Schleich, and E. M. Rasel, Atom-chip fountain gravimeter, Phys. Rev. Lett. 117, 203003 (2016).
  * Asenbaum  _et al._ [2017] P. Asenbaum, C. Overstreet, T. Kovachy, D. D. Brown, J. M. Hogan, and M. A. Kasevich, Phase shift in an atom interferometer due to spacetime curvature across its wave function, Phys. Rev. Lett. 118, 183602 (2017).
  * Xu  _et al._ [2019] V. Xu, M. Jaffe, C. D. Panda, S. L. Kristensen, L. W. Clark, and H. Müller, Probing gravity by holding atoms for 20 seconds, Science 366, 745 (2019).
  * Eibenberger  _et al._ [2013] S. Eibenberger, S. Gerlich, M. Arndt, M. Mayor, and J. Tüxen, Matter–wave interference of particles selected from a molecular library with masses exceeding 10000 amu, Phys. Chem. Chem. Phys. 15, 14696 (2013).
  * Fein  _et al._ [2019] Y. Y. Fein, P. Geyer, P. Zwick, F. Kiałka, S. Pedalino, M. Mayor, S. Gerlich, and M. Arndt, Quantum superposition of molecules beyond 25 kda, Nat. Phys. 15, 1242 (2019).
  * Stickler  _et al._ [2018] B. A. Stickler, B. Papendell, S. Kuhn, B. Schrinski, J. Millen, M. Arndt, and K. Hornberger, Probing macroscopic quantum superpositions with nanorotors, New J. Phys. 20, 122001 (2018).
  * Kaltenbaek  _et al._ [2022] R. Kaltenbaek _et al._ , Research campaign: Macroscopic quantum resonators (MAQRO), Quantum Sci. Technol. 10.1088/2058-9565/aca3cd (2022).
  * Chu  _et al._ [2018] Y. Chu, P. Kharel, T. Yoon, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, Creation and control of multi-phonon Fock states in a bulk acoustic-wave resonator, Nature 563, 666 (2018).
  * Satzinger  _et al._ [2018] K. J. Satzinger, Y. P. Zhong, H. S. Chang, G. A. Peairs, A. Bienfait, M.-H. Chou, A. Y. Cleland, C. R. Conner, É. Dumur, J. Grebel, I. Gutierrez, B. H. November, R. G. Povey, S. J. Whiteley, D. D. Awschalom, D. I. Schuster, and A. N. Cleland, Quantum control of surface acoustic-wave phonons, Nature 563, 661 (2018).
  * Wollack  _et al._ [2022] E. A. Wollack, A. Y. Cleland, R. G. Gruenke, Z. Wang, P. Arrangoiz-Arriola, and A. H. Safavi-Naeini, Quantum state preparation and tomography of entangled mechanical resonators, Nature 604, 463 (2022).
  * von Lüpke  _et al._ [2022] U. von Lüpke, Y. Yang, M. Bild, L. Michaud, M. Fadel, and Y. Chu, Parity measurement in the strong dispersive regime of circuit quantum acoustodynamics, Nat. Phys. 18, 794 (2022).
  * Bild  _et al._ [2022] M. Bild, M. Fadel, Y. Yang, U. von Lüpke, P. Martin, A. Bruno, and Y. Chu, Schrödinger cat states of a 16-microgram mechanical oscillator, arXiv preprint arXiv:2211.00449 (2022).
  * Nimmrichter and Hornberger [2013] S. Nimmrichter and K. Hornberger, Macroscopicity of mechanical quantum superposition states, Phys. Rev. Lett. 110, 160403 (2013).
  * Schrinski  _et al._ [2019] B. Schrinski, S. Nimmrichter, B. A. Stickler, and K. Hornberger, Macroscopicity of quantum mechanical superposition tests via hypothesis falsification, Phys. Rev. A 100, 032111 (2019).
  * Note [1] Note that the macrorealist modification hypothesis can be probed also in experiments that to do not involve quantum superpositons. Only observations of genuine quantum signatures can therefore be assigned a benchmark value $\mu$ [17]. This requires the interrogated state not to have an underlying classical description, excluding all oscillator states with positive Wigner function.
  * Royer [1977] A. Royer, Wigner function as the expectation value of a parity operator, Phys. Rev. A 15, 449 (1977).
  * Bassi  _et al._ [2013] A. Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht, Models of wave-function collapse, underlying theories, and experimental tests, Rev. Mod. Phys. 85, 471 (2013).
  * Carlesso  _et al._ [2022] M. Carlesso, S. Donadi, L. Ferialdi, M. Paternostro, H. Ulbricht, and A. Bassi, Present status and future challenges of non-interferometric tests of collapse models, Nat. Phys. 18, 243 (2022).
  * [22] See Supplementary Information for details on the experimental setup and the theoretical calculations, including Refs. [27, 28, 29, 30, 31, 32, 33].
  * Tobar  _et al._ [2022] G. Tobar, S. Forstner, A. Federov, and W. P. Bowen, Testing spontaneous wavefunction collapse with quantum electromechanics, arXiv preprint arXiv:2206.14531 (2022).
  * Seis  _et al._ [2022] Y. Seis, T. Capelle, E. Langman, S. Saarinen, E. Planz, and A. Schliesser, Ground state cooling of an ultracoherent electromechanical system, Nat. Comm. 13, 1507 (2022).
  * Note [2] This rough estimate is based on [17], by approximating the mode in [24] as the fundamental Bessel mode $\bm{u}(\bm{r})=J_{0}(4.81\sqrt{y^{2}+z^{2}}/d)\mathbf{e}_{x}$, for a disc with diameter $d=100\,\mu$m, $14\,$nm thickness, and $m_{\mathrm{eff}}\simeq 68\,$pg. Precise $\mu$-values would depend on a combination of effective mass, mode profile, and observed duration of Wigner negativities.
  * Schrinski  _et al._ [2020] B. Schrinski, S. Nimmrichter, and K. Hornberger, Quantum-classical hypothesis tests in macroscopic matter-wave interferometry, Phys. Rev. Res. 2, 033034 (2020).
  * Jeffreys [1946] H. Jeffreys, An invariant form for the prior probability in estimation problems, Proc. R. Soc. Lond. A 186, 453 (1946).
  * Arrangoiz-Arriola [2019] P. Arrangoiz-Arriola,  _Quantum acoustics with lithium niobate nanostructures_ , Ph.D. thesis, Stanford University (2019).
  * Satzinger [2018] K. J. Satzinger,  _Quantum control of surface acoustic wave phonons_ , Ph.D. thesis, UC Santa Barbara (2018).
  * [30] K. J. Satzinger, private communication.
  * Chu  _et al._ [2017] Y. Chu, P. Kharel, W. H. Renninger, L. D. Burkhart, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, Quantum acoustics with superconducting qubits, Science 358, 199 (2017).
  * Armano  _et al._ [2016] M. Armano, H. Audley, G. Auger, J. Baird, M. Bassan, P. Binetruy, M. Born, D. Bortoluzzi, N. Brandt, M. Caleno, _et al._ , Sub-femto-g free fall for space-based gravitational wave observatories: Lisa pathfinder results, Phys. Rev. Lett. 116, 231101 (2016).
  * Carlesso  _et al._ [2016] M. Carlesso, A. Bassi, P. Falferi, and A. Vinante, Experimental bounds on collapse models from gravitational wave detectors, Phys. Rev. D 94, 124036 (2016).

## Supplementary information

##  Supplementary Section A Geometric factor

The considered macrorealistic modifications lead to momentum diffusion in the
oscillator state. The momentum diffusion rate $D=\hbar^{2}U/\tau_{e}$ is
inversely proportional to the coherence time parameter $\tau_{e}$ and depends
on the oscillator mode function $\bm{u}(\bm{r})$ through the geometric factor
[17]

$U=\frac{1}{2\hbar^{2}}\int
d^{3}q\,g(\bm{q},\sigma_{q})|\tilde{\bm{u}}(\bm{q})\cdot\bm{q}|^{2}\;,\qquad\text{with}\qquad\tilde{\bm{u}}(\bm{q})=\dfrac{1}{m_{e}}\int
d^{3}r\,\varrho(\bm{r})\bm{u}(\bm{r})e^{-\bm{r}\cdot\bm{q}/\hbar}\;.$ (8)

Here, $m_{e}$ denotes the reference electron mass, and
$g(\bm{q},\sigma_{q})=\exp(-|\bm{q}|^{2}/2\sigma_{q}^{2})/(2\pi\sigma_{q}^{2})^{3/2}$
a Gaussian momentum distribution of standard deviation $\sigma_{q}$. The
resonator mode employed in this work can be modelled to great accuracy with
the function

$\displaystyle\bm{u}(\bm{r})=\mathcal{N}e^{-(y^{2}+z^{2})/w_{0}^{2}}\cos\left(\pi\dfrac{\ell
x}{L}\right)\mathbf{e}_{x}\;,$ (9)

where $w_{0}$ is the beam waist ($1/e$ amplitude radius), $L$ is the height,
and $\ell$ the mode index. Based on this mode function, we can define an
effective mode volume and mass,
$V_{\mathrm{eff}}=\int\mathrm{d}^{3}r\,|\bm{u}(\bm{r})|^{2}$ and
$m_{\mathrm{eff}}=\int\mathrm{d}^{3}r\,\varrho(\bm{r})|\bm{u}(\bm{r})|^{2}$.
They both depend on the square of the amplitude prefactor $\mathcal{N}$, but
the relevant physical displacement amplitudes of the individual atoms,
$x_{0}\bm{u}(\bm{r}_{n})=\sqrt{\hbar/m_{\mathrm{eff}}\omega}\bm{u}(\bm{r}_{n})$,
as well as the dimensionless macrorealistic diffusion rate,
$\Gamma\tau_{e}=D(x_{0}/\hbar)^{2}=Ux_{0}^{2}$, do not depend on
$\mathcal{N}$. Hence we may set $\mathcal{N}=1$, without loss of generality,
as done in the main text.

Given that we consider a displacement field with only a longitudinal component
in $x$-direction, $\bm{u}(\bm{r})=u_{x}(\bm{r})\mathbf{e}_{x}$, and that the
mass density distribution is assumed to be uniform (i.e.
$\varrho(\bm{r})=\mathrm{bar}{\varrho}$), Eq. (8) can be simplified further:

$\displaystyle U=\frac{1}{2\hbar^{2}}\int
d^{3}\bm{q}\,g(\bm{q},\sigma_{q})|\tilde{u}_{x}(\bm{q})q_{x}|^{2}$
$\displaystyle=\frac{\mathrm{bar}{\varrho}^{2}}{2\hbar^{2}m_{e}^{2}}\int
d^{3}\bm{q}\,g(\bm{q},\sigma_{q})q_{x}^{2}\left|\int
dy\,dz\,e^{-i(q_{y}y+q_{z}z)/\hbar-(y^{2}+z^{2})/w_{0}^{2}}\int_{0}^{L}dx\,\cos\left(\pi\frac{\ell
x}{L}\right)e^{-iq_{x}x/\hbar}\right|^{2}$
$\displaystyle=\frac{\mathrm{bar}{\varrho}^{2}\pi^{2}w_{0}^{4}L^{2}}{2\hbar^{2}m_{e}^{2}}\int
d^{3}\bm{q}\,g(\bm{q},\sigma_{q})q_{x}^{2}e^{-(q_{y}^{2}+q_{z}^{2})w_{0}^{2}/2\hbar^{2}}\left|\int_{0}^{1}\mathrm{d}\phi\cos(\pi\ell\phi)e^{-i\phi
q_{x}L/\hbar}\right|^{2}.$ (10)

This allows us to first integrate out the momenta $\bm{q}$, change to
cylindrical coordinates and express the end result as

$\displaystyle
U=\frac{\mathrm{bar}{\varrho}^{2}}{2m_{e}^{2}}\frac{\pi^{2}w_{0}^{4}}{1+\sigma_{w}^{2}}\int\mathrm{d}\zeta\frac{\zeta^{2}}{\sqrt{2\pi}\sigma_{L}}e^{-\zeta^{2}/2\sigma_{L}^{2}}\left|\int_{0}^{1}\mathrm{d}\phi\cos(\pi\ell\phi)e^{-i\phi\zeta}\right|^{2}$
$\displaystyle=\frac{\mathrm{bar}{\varrho}^{2}\pi^{2}w_{0}^{4}}{m_{e}^{2}(1+\sigma_{w}^{2})}\int\mathrm{d}\zeta\,\frac{e^{-\zeta^{2}/2\sigma_{L}^{2}}}{\sqrt{2\pi}\sigma_{L}}\frac{1-(-)^{\ell}\cos\zeta}{\left(1-\frac{\pi^{2}\ell^{2}}{\zeta^{2}}\right)^{2}}\;.$
(11)

with $\sigma_{L}=L\sigma_{q}/\hbar$ and $\sigma_{w}=w_{0}\sigma_{q}/\hbar$.
Once we identify $m_{\mathrm{eff}}=\mathrm{bar}\varrho\pi w_{0}^{2}L/4$, we
arrive at the expression for the diffusion rate $\Gamma=Ux_{0}^{2}/\tau_{e}$
given in the main text. A non-tight upper bound follows by pulling the
absolute value into the $\phi$-integral in (11). The $\phi$-integral then
evaluates to $2/\pi$, so that
$U<2\mathrm{bar}{\varrho}^{2}w_{0}^{4}\sigma_{L}^{2}/m_{e}^{2}(1+\sigma_{w}^{2})$,
and
$\Gamma<(32/\pi^{2})(m_{\mathrm{eff}}^{2}/m_{e}^{2}\tau_{e})(\sigma_{q}x_{0}/\hbar)^{2}/(1+\sigma_{w}^{2})$.
The scaling of this bound with the mode mass and geometry is deceivingly
similar to previous results for center-of-mass oscillations. However, we will
see in the following that the bound is rather loose and the actual diffusion
rate for bulk acoustic modes exhibits a completely different scaling.

In the limiting case $\sigma_{L},\sigma_{w}\ll 1$, the Gaussian function is
sharply peaked around $\zeta=0$, and we can expand the remainder of the
integrand to lowest non-vanishing order; for even $\ell$, this yields
$U\approx
15\mathrm{bar}{\varrho}^{2}w_{0}^{4}\sigma_{L}^{6}/2\pi^{2}m_{e}^{2}\ell^{4}$,
and $\Gamma\approx
120(m_{\mathrm{eff}}^{2}/m_{e}^{2}\tau_{e})(\sigma_{q}x_{0}/\hbar)^{2}(\sigma_{L}/\pi\ell)^{4}$.
For decreasing $\sigma_{q}$, this rate drops as $\sigma_{q}^{6}$. Notice that
odd $\ell$ would result in a $\sigma_{q}^{4}$-scaling, $U\approx
3\mathrm{bar}{\varrho}^{2}w_{0}^{4}\sigma_{L}^{4}/\pi^{2}m_{e}^{2}\ell^{4}$.

Greater care must be taken in the opposite limit, $\sigma_{L},\sigma_{w}\gg
1$. Although the non-Gaussian part of the integrand in (11) is clearly peaked
around $\zeta=\pm\pi\ell$, it is asymptotically of order $\zeta^{0}$ and thus
non-integrable without the regularizing presence of the Gaussian. Hence a
dominant contribution to the $\zeta$-integral in (11) does not come from the
vicinity of the peaks, but from $|\zeta|\gg\pi\ell$. For even greater values
$\sigma_{L}\gg\pi\ell$, we can approximate to lowest order

$\displaystyle U$ $\displaystyle\approx
U^{(0)}:=\frac{\mathrm{bar}{\varrho}^{2}\pi^{2}w_{0}^{4}}{m_{e}^{2}(1+\sigma_{w}^{2})}\int\mathrm{d}\zeta\,\frac{e^{-\zeta^{2}/2\sigma_{L}^{2}}}{\sqrt{2\pi}\sigma_{L}}\left[1-(-)^{\ell}\cos\zeta\right]=\frac{\mathrm{bar}{\varrho}^{2}\pi^{2}w_{0}^{4}}{m_{e}^{2}(1+\sigma_{w}^{2})}\left[1-(-)^{\ell}e^{-\sigma_{L}^{2}/2}\right]\approx\frac{\mathrm{bar}{\varrho}^{2}\pi^{2}w_{0}^{4}}{m_{e}^{2}\sigma_{w}^{2}}.$
(12)

The corresponding diffusion rate reads as $\Gamma\approx
16(m_{\mathrm{eff}}^{2}/m_{e}^{2}\tau_{e})(x_{0}/\sigma_{L}w_{0})^{2}$, which
does not depend on $\ell$, is proportional to $m_{\mathrm{eff}}$, and is again
monotonic in $\sigma_{q}$. In order to access the most relevant intermediate
$\sigma_{q}$-regime of maximum diffusion, consider the correction
$U=U^{(0)}+U^{(1)}$, which can be expressed as

$\displaystyle
U^{(1)}=\frac{\mathrm{bar}{\varrho}^{2}\pi^{4}\ell^{2}w_{0}^{4}}{m_{e}^{2}(1+\sigma_{w}^{2})}\int\mathrm{d}\zeta\,\frac{e^{-\zeta^{2}/2\sigma_{L}^{2}}}{\sqrt{2\pi}\sigma_{L}}\left[\frac{2\zeta^{2}-\pi^{2}\ell^{2}}{2(\zeta+\pi\ell)^{2}}\mathrm{sinc}^{2}\left(\frac{\zeta-\pi\ell}{2}\right)\right]\;.$
(13)

The non-Gaussian part of this integrand (square-bracketed) is not only sharply
peaked at $\zeta=\pm\pi\ell$, but also asymptotically converges to zero like
$\zeta^{-2}$. The integral of this expression evaluates to $\pi/2$. Assuming
$\pi\ell\gg 1$, we can now replace $\zeta^{2}$ in the Gaussian by
$\pi^{2}\ell^{2}$, so that

$\displaystyle
U^{(1)}\approx\frac{\mathrm{bar}{\varrho}^{2}\pi^{5}\ell^{2}w_{0}^{4}}{2m_{e}^{2}(1+\sigma_{w}^{2})}\frac{e^{-\pi^{2}\ell^{2}/2\sigma_{L}^{2}}}{\sqrt{2\pi}\sigma_{L}}\approx\frac{\mathrm{bar}{\varrho}^{2}\pi^{5}\ell^{2}w_{0}^{4}e^{-\pi^{2}\ell^{2}/2\sigma_{L}^{2}}}{2\sqrt{2\pi}m_{e}^{2}\sigma_{w}^{2}\sigma_{L}}\;,\quad
U\approx\frac{\mathrm{bar}{\varrho}^{2}\pi^{2}w_{0}^{2}L^{2}}{m_{e}^{2}\sigma_{L}^{2}}\left(1+\frac{\pi^{3}\ell^{2}e^{-\pi^{2}\ell^{2}/2\sigma_{L}^{2}}}{2\sqrt{2\pi}\sigma_{L}}\right)=\frac{16m_{\mathrm{eff}}^{2}}{m_{e}^{2}\sigma_{L}^{2}w_{0}^{2}}\left(1+\frac{\pi^{3}\ell^{2}e^{-\pi^{2}\ell^{2}/2\sigma_{L}^{2}}}{2\sqrt{2\pi}\sigma_{L}}\right)$
(14)

Presuming $\sigma_{w}\gg 1$, the local maximum of $U$ coincides approximately
with that of $U^{(1)}$, at $\sigma_{L}\approx\pi\ell/\sqrt{3}$. It marks the
$\sigma_{q}$-value at which the experiment is most sensitive to macrorealistic
diffusion. The maximum rate reads as

$\max_{\sigma_{q}}\Gamma\approx\frac{48m_{\mathrm{eff}}^{2}x_{0}^{2}}{m_{e}^{2}\tau_{e}\pi^{2}\ell^{2}w_{0}^{2}}\left(1+\sqrt{\frac{3\pi}{8e^{3}}}\pi\ell\right)\stackrel{{\scriptstyle\pi\ell\gg
1}}{{\approx}}\sqrt{\frac{3\pi}{2e^{3}}}\frac{24}{\pi\ell\tau_{e}}\left(\frac{m_{\mathrm{eff}}x_{0}}{m_{e}w_{0}}\right)^{2}=\sqrt{\frac{3\pi}{2e^{3}}}\frac{6\hbar\mathrm{bar}{\varrho}}{m_{e}^{2}\omega\tau_{e}}\frac{L}{\ell}\;.$
(15)

For numerical evaluation, an exact expression of the diffusion rate can be
given in terms of

$U=\frac{\mathrm{bar}{\varrho}^{2}}{2m_{e}^{2}}\frac{\pi^{2}w_{0}^{4}\sigma_{L}^{2}}{1+\sigma_{w}^{2}}f_{\ell}(\sigma_{L})\;,\qquad\text{with}\qquad
f_{\ell}(\xi)=(1+\xi\partial_{\xi})\int_{[0,1]^{2}}\\!\\!\mathrm{d}\phi_{1}\mathrm{d}\phi_{2}\,\cos\ell\pi\phi_{1}\cos\ell\pi\phi_{2}\,e^{-\xi^{2}(\phi_{1}-\phi_{2})^{2}/2}\;.$
(16)

The function $f_{\ell}$ can be evaluated to

$\displaystyle
f_{\ell}(\sigma_{L})=\frac{1}{\sqrt{2}\sigma_{L}^{5}}\left\\{-\ell\pi
h(\sigma_{L},0,\ell)+\sqrt{2}\sigma_{L}\left[(-)^{\ell}e^{-\sigma_{L}^{2}/2}-1\right]\left(\ell^{2}\pi^{2}-2\sigma_{L}^{2}\right)+(-)^{\ell}\ell\pi^{2}\mathrm{Re}\left[h(\sigma_{L},\sigma_{L},\ell)\right]\right\\}\;,$
(17)

where

$\displaystyle h(a,b,\ell)=-i\frac{\sqrt{\pi}}{2}[6a^{2}-2\ell\pi(\ell\pi-
ib^{2})]\mathrm{erf}\left[\frac{i(\ell\pi/a-ib)}{\sqrt{2}}\right]\exp\left[-\frac{(\ell\pi/a-ib)^{2}}{2}\right]\;.$
(18)

In Fig. S3, we plot the resulting diffusion rate $\Gamma=D/\hbar
m_{\mathrm{eff}}\omega=\hbar U/m_{\mathrm{eff}}\omega\tau_{e}$ as a function
of the modification parameter $\sigma_{q}$. We compare the mode addressed in
our experiment to one with a much lower frequency and mode index. Based on the
above reasoning, the maxima are expected to be at $\hbar/\sigma_{q}\approx
0.5\,\mu$m ($\ell=486$) and $17\,\mu$m ($\ell=8$), which roughly agrees with
the exact positions on the logarithmic scale. The dotted line represents
$U^{(1)}$ from (14), a good approximation in the relevant regime around
maximum diffusion. The dashed lines show the asymptotic approximations for
$\sigma_{L}\gg\pi\ell$ (left) and for $\sigma_{L}\ll 1$ (right).

Figure S3: Dimensionless diffusion rate $\Gamma\tau_{e}$ as function of the
modification parameter $\sigma_{q}$ for the device studied in the main text
(blue curve). We compare it to the correction $U^{(1)}$ from (14) (dotted
line), as well as to the asymptote $U^{(0)}$ for $\sigma_{L}\gg\pi\ell$ (left
dashed) and to the asymptote for $\sigma_{L}\ll 1$ (right dashed). Much
stronger diffusion would arise (red curve) if instead of the mode at
$\ell=486$, we could excite, say, the mode with index $\ell=8$ at the
frequency $\omega=2\pi\cdot$98\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$.

##  Supplementary Section B Initial states and time evolution

In our setup, macrorealistic collapse models are put to the test by monitoring
the decay of Fock and superposition states over a finite amount of time. To
assess the influence of those models, it is convenient to solve the time
evolution in the Wigner phase space representation of the oscillator state.
The one-dimensional Wigner function is defined as

$\displaystyle W(x,p)=\frac{1}{2\pi\hbar}\int\mathrm{d}s\,e^{ips/\hbar}\langle
x-s/2|\rho|x+s/2\rangle\;,$ (19)

where $x$ and $p$ are the position and momentum coordinates of the harmonic
oscillator. For better readability we now change to dimensionless units
$X=\sqrt{m_{\mathrm{eff}}\omega/\hbar}x$ and $P=\sqrt{1/\hbar
m_{\mathrm{eff}}\omega}p$. Then the Wigner functions (19) of the ground state
$|0\rangle$, first exited Fock state $|1\rangle$, and the superposition state
$(|0\rangle+|1\rangle)/\sqrt{2}$, respectively, read as

$\displaystyle W_{0}(X,P)=\frac{1}{\pi}e^{-X^{2}-P^{2}},\quad
W_{1}(X,P)=\frac{1}{\pi}\left[2X^{2}+2P^{2}-1\right]e^{-X^{2}-P^{2}},\quad\mathrm{and}\quad
W_{10}(X,P)=\frac{1}{2\pi}\left[(1+\sqrt{2}X)^{2}-1+2P^{2}\right]e^{-X^{2}-P^{2}}.$
(20)

The time evolution of the state described by the master equation (6) in the
main text describes the impact of the macrorealistic diffusion at the rate
$\Gamma$, averaged over the free oscillation, and it includes also an
environmental decay channel at the rate $\gamma_{\downarrow}$. The two
separate Lindblad generators for decay and amplification,
$\partial_{t}^{\mathrm{dec}}\rho=\Gamma\mathcal{D}[\mathsf{a}]\rho$ and
$\partial_{t}^{\mathrm{amp}}\rho=\Gamma\mathcal{D}[\mathsf{a}^{\dagger}]\rho$
with $\mathsf{a}=(X+iP)/\sqrt{2}$, translate into phase space as

$\displaystyle\partial_{t}^{\mathrm{dec}}W(X,P)$
$\displaystyle=\frac{\Gamma}{2\pi}\int\mathrm{d}s\,e^{iPs}\langle
X-s/2|\mathsf{a}\rho\mathsf{a}^{\dagger}-\frac{1}{2}\\{\mathsf{a}^{\dagger}\mathsf{a},\rho\\}|X+s/2\rangle=\frac{\Gamma}{4}\left[\partial_{X}^{2}+\partial_{P}^{2}+2\partial_{X}X+2\partial_{P}P\right]W(X,P)\;,$
(21) $\displaystyle\partial_{t}^{\mathrm{amp}}W(X,P)$
$\displaystyle=\frac{\Gamma}{4}\left[\partial_{X}^{2}+\partial_{P}^{2}-2\partial_{X}X-2\partial_{P}P\right]W(X,P)\;.$
(22)

Both are classical Fokker-Planck equations, and the first one describes
thermalization at the (dimensionless) energy $1/2$, i.e. the ground state
energy of the quantum oscillator. The master equation then yields

$\displaystyle\partial_{t}W(X,P)=\frac{\Gamma}{2}\left[\partial_{X}^{2}+\partial_{P}^{2}\right]+\frac{{\gamma_{\downarrow}}}{4}\left[\partial_{X}^{2}+\partial_{P}^{2}+2\partial_{X}X+2\partial_{P}P\right]W(X,P).$
(23)

The solution to this Fokker-Planck equation is obtained with help of a Fourier
transform into the characteristic function representation,

$\displaystyle\chi(x,p)=\int\mathrm{d}x^{\prime}\mathrm{d}p^{\prime}\,W(x^{\prime},p^{\prime})e^{i(xp^{\prime}-px^{\prime})/\hbar}\;,$
(24)

which evolves according to

$\displaystyle\partial_{t}\chi(X,P)=\frac{1}{2}\left[-\left(\frac{{\gamma_{\downarrow}}}{2}+\Gamma\right)(X^{2}+P^{2})-{\gamma_{\downarrow}}(X\partial_{X}+P\partial_{P})\right]\chi(X,P)\;.$
(25)

We can now separate the simple relaxation dynamics (an exponential contraction
toward the origin) from the diffusion. Let
$X=\tilde{X}e^{{\gamma_{\downarrow}}t/2}$,
$P=\tilde{P}e^{{\gamma_{\downarrow}}t/2}$, and
$\Phi(\tilde{X},\tilde{P};t)=\chi(\tilde{X}e^{{\gamma_{\downarrow}}t/2},\tilde{P}e^{{\gamma_{\downarrow}}t/2};t)$,
which then obeys

$\displaystyle\partial_{t}\Phi(\tilde{X},\tilde{P};t)$
$\displaystyle=-\frac{{\gamma_{\downarrow}}+2\Gamma}{4}e^{{\gamma_{\downarrow}}t}(\tilde{X}^{2}+\tilde{P}^{2})\Phi(\tilde{X},\tilde{P};t)\;,$
(26) as solved by $\displaystyle\Phi(\tilde{X},\tilde{P};t)$
$\displaystyle=\exp\left[-\left(1+\frac{2\Gamma}{{\gamma_{\downarrow}}}\right)(e^{{\gamma_{\downarrow}}t}-1)\frac{\tilde{X}^{2}+\tilde{P}^{2}}{4}\right]\Phi(\tilde{X},\tilde{P};0)\;,$
(27) so that $\displaystyle\quad\chi(X,P;t)$
$\displaystyle=\exp\left[-\left(1+\frac{2\Gamma}{{\gamma_{\downarrow}}}\right)(1-e^{-{\gamma_{\downarrow}}t})\frac{X^{2}+P^{2}}{4}\right]\chi(Xe^{-{\gamma_{\downarrow}}t/2},Pe^{-{\gamma_{\downarrow}}t/2};0)\;.$
(28)

From this expression it is evident that the time evolution of an initial
$\chi(X,P;0)$ is simply given by a multiplication with a Gaussian exponential
combined with an isotropic rescaling of the arguments. Noting that
$\chi(0,0;t)=1$ at any time, one can check that the asymptotic steady state is
a Gaussian,
$\chi(X,P;\infty)=\exp[-(1+2\Gamma/{\gamma_{\downarrow}})(X^{2}+P^{2})/4]$,
which for $\Gamma=0$ is the ground state of the harmonic oscillator.

Back in the Wigner representation, the solution turns into a rescaling
combined with a Gaussian convolution. Introducing the dimensionless quantitiy
$\tilde{T}=1/2+\Gamma/{\gamma_{\downarrow}}$, we get

$\displaystyle W(X,P;t)$
$\displaystyle=\frac{1}{4\pi^{2}}\int\mathrm{d}\tilde{X}\mathrm{d}\tilde{P}\,e^{i\tilde{P}X-iP\tilde{X}}\exp\left[-\tilde{T}(1-e^{-{\gamma_{\downarrow}}t})\frac{\tilde{X}^{2}+\tilde{P}^{2}}{2}\right]\chi(\tilde{X}e^{-{\gamma_{\downarrow}}t/2},\tilde{P}e^{-{\gamma_{\downarrow}}t/2};0)$
$\displaystyle=\frac{1}{4\pi^{2}}\int\mathrm{d}\tilde{X}\mathrm{d}\tilde{P}\mathrm{d}X_{0}\mathrm{d}P_{0}\,W(X_{0},P_{0};0)\exp\left[-\tilde{T}(1-e^{-{\gamma_{\downarrow}}t})\frac{\tilde{X}^{2}+\tilde{P}^{2}}{2}+i\tilde{X}(P_{0}e^{-{\gamma_{\downarrow}}t/2}-P)-i\tilde{P}(X_{0}e^{-{\gamma_{\downarrow}}t/2}-X)\right]$
$\displaystyle=\frac{1}{2\pi\tilde{T}(1-e^{-{\gamma_{\downarrow}}t})}\int\mathrm{d}X_{0}\mathrm{d}P_{0}\,W(X_{0},P_{0};0)\exp\left[-\frac{(X-X_{0}e^{-{\gamma_{\downarrow}}t/2})^{2}+(P-P_{0}e^{-{\gamma_{\downarrow}}t/2})^{2}}{2\tilde{T}(1-e^{-{\gamma_{\downarrow}}t})}\right]$
$\displaystyle=\frac{e^{{\gamma_{\downarrow}}t}}{2\pi\tilde{T}(1-e^{-{\gamma_{\downarrow}}t})}\int\mathrm{d}X_{0}\mathrm{d}P_{0}\,W(X_{0}e^{{\gamma_{\downarrow}}t/2},P_{0}e^{{\gamma_{\downarrow}}t/2};0)\exp\left[-\frac{(X-X_{0})^{2}+(P-P_{0})^{2}}{2\tilde{T}(1-e^{-{\gamma_{\downarrow}}t})}\right]\;.$
(29)

The two last lines are alternative ways to write the convolution. Given the
simple Gaussian form of the initial Wigner functions (20), the integrals can
be calculated analytically and result in

$\displaystyle W_{0}(X,P;t)=$ $\displaystyle\frac{1}{\pi
R(t)}\exp\left[{\gamma_{\downarrow}}t-\frac{e^{{\gamma_{\downarrow}}t}(X^{2}+P^{2})}{R(t)}\right]\;,$
(30) $\displaystyle W_{1}(X,P;t)=$
$\displaystyle\frac{4\tilde{T}^{2}+2e^{{\gamma_{\downarrow}}t}(X^{2}+P^{2}+2(e^{{\gamma_{\downarrow}}t}-2)\tilde{T}^{2})-1}{\pi
R^{3}(t)}\exp\left[\frac{{\gamma_{\downarrow}}t(1-2\tilde{T})-e^{{\gamma_{\downarrow}}t}(X^{2}+P^{2}-2\tilde{T}{\gamma_{\downarrow}}t)}{R(t)}\right]\;,$
(31) $\displaystyle W_{10}(X,P;t)=$
$\displaystyle\frac{\sqrt{8}e^{3{\gamma_{\downarrow}}t/2}X\tilde{T}+4e^{2{\gamma_{\downarrow}}t}\tilde{T}^{2}-\sqrt{2}e^{{\gamma_{\downarrow}}t/2}X(2\tilde{T}-1)+2\tilde{T}(2\tilde{T}-1)+e^{{\gamma_{\downarrow}}t}(P^{2}+X^{2}+2\tilde{T}-8\tilde{T}^{2})}{\pi
R^{3}(t)}$
$\displaystyle\times\exp\left[\frac{{\gamma_{\downarrow}}t(1-2\tilde{T})-e^{{\gamma_{\downarrow}}t}(X^{2}+P^{2}-2\tilde{T}{\gamma_{\downarrow}}t)}{R(t)}\right]\;,$
(32)

with $R(t)=1+2(e^{{\gamma_{\downarrow}}t}-1)\tilde{T}$. These expressions can
be directly compared to the experimentally reconstructed Wigner function,
which allows us to perform Bayesian parameter estimation of the macrorealistic
$\Gamma$-values that are compatible with the measured data.

In Fig. S4, we use Eqs. (31), (32) to illustrate the evolution of states
$\ket{1}$ and $(\ket{0}+\ket{1})/\sqrt{2}$ according to either pure
relaxation, pure diffusion, or a combination of both. Notice the significant
differences in the dynamics between the different cases.

Figure S4: Comparison between relaxation, diffusion, and a combination of
both. Rows 1 and 4: State $\ket{1}$ and $(\ket{0}+\ket{1})/\sqrt{2}$ evolved
according to a relaxation process with
$1/\gamma_{\downarrow}=$40\text{\,}\mu\mathrm{s}$$. Rows 2 and 5: same initial
states, but evolved according to a purely dissipative process with
$\Gamma=$10^{4}\text{\,}\mathrm{H}\mathrm{z}$$. Rows 3 and 6: evolution of the
initial states according to a combination of both relaxation and diffusion
with same rates. Note that in the actual experiment we have
${\gamma_{\downarrow}}\gg\Gamma$, thus diffusion contributes very little
compared to relaxation.

##  Supplementary Section C Bayesian parameter estimation

By measuring the Wigner function at several instances of time and comparing it
to the theoretical predictions in Sect. B, one can infer the most likely
values for the macrorealistic diffusion rate $\Gamma$ and, subsequently, for
the values of the time parameter $\tau_{e}=\hbar U/\Gamma
m_{\mathrm{eff}}\omega$ by virtue of Bayes’ rule. The posterior distribution
of $\tau_{e}$-values given the measured data $d$, a fixed value of the
parameter $\sigma_{q}$, and an independently inferred decay time
${1/\gamma_{\downarrow}\approx T_{1}}$, updates as [17]

$\displaystyle p(\tau_{e}|d;\sigma_{q},T_{1})\propto
p(d|\tau_{e};\sigma_{q},T_{1})p(\tau_{e}|\sigma_{q},T_{1}),$ (33)

Here, $d=[d_{ij}(t)]$ subsumes all recorded pixel values
$d_{ij}(t)\in\mathbb{R}$ of the two-dimensional Wigner function
reconstructions at the times $t=10,20,$40\text{\,}\mu\mathrm{s}$$, as plotted
in the main text. The likelihood $p(d|\tau_{e};\sigma_{q},T_{1})$ for $d$
follows from the theoretical model for the Wigner function evaluated at the
recorded pixel coordinates and times, $W(X_{i},P_{j};t)$, and the assumption
of additional Gaussian noise in each data point. That is, we model the
likelihood for each pixel value by a Gaussian distribution of standard
deviation $s$, so that

$\displaystyle p(d|\tau_{e};\sigma_{q},T_{1})=\prod_{i,j,t}\frac{1}{\sqrt{2\pi
s^{2}}}\exp\left\\{-\frac{\left[d_{ij}(t)-W(X_{i},P_{j};t)\right]^{2}}{2s^{2}}\right\\}\;.$
(34)

We extract the overall noise level $s$ by taking a sample of all pixel values
for the Fock-$\ket{1}$ measurements at $t=0,10,20,$40\text{\,}\mu\mathrm{s}$$
and subtracting the respective model Wigner functions $W_{1}(X,P;t)$ with
decay time $T_{1}$ and without macrorealistic diffusion ($\Gamma=0$). A
histogram of the resulting pixel deviations
$\Delta_{ij}(t)=d_{ij}(t)-W_{1}(X_{i},P_{j};t)$ is shown in Fig. S5(a),
together with a Gaussian fit that yields the noise level $s=0.034$.

The prior $p(\tau_{e}|\sigma_{q},T_{1})$ is chosen as Jeffreys’ prior [27],
which is maximally objective and fair in a comparison with other completely
different experiments [17]. It is defined as

$\displaystyle
p(\tau_{e}|\sigma_{q},T_{1})\propto\sqrt{\mathcal{I}(\tau_{e}|\sigma_{q},T_{1})}=\sqrt{\left\langle\left(\frac{\partial}{\partial\tau_{e}}\log
p(d|\tau_{e};\sigma_{q},T_{1})\right)^{2}\right\rangle_{d}},$ (35)

where $\mathcal{I}(\tau_{e}|\sigma_{q},I)$ is the Fisher information of the
likelihood with respect to $\tau_{e}$ and $\langle\cdot\rangle_{d}$ denotes
the expectation value over all possible measurement results,
$d_{ij}(t)\in\mathbb{R}$.

We perform the Bayesian update (33) with the data obtained for the three time
snapshots at $10,20,$40\text{\,}\mu\mathrm{s}$$. Figures S5(b) and (c) show
the resulting posterior distributions over $\Gamma$-values for the
Fock-$\ket{1}$ and the superposition state, respectively. The upper five-
percent quantiles for $\Gamma$, marked by vertical lines, correspond to lower
five-percent quantiles of $\tau_{e}=\hbar U/\Gamma m_{\mathrm{eff}}\omega$ at
a given value for $\sigma_{q}$. Smaller time parameters are thus ruled out
with 95% confidence. We find the threshold values of $\Gamma=$1.6\cdot
10^{2}\text{\,}\mathrm{s}^{-1}$$ and $6.4\cdot 10^{2}\text{\,}\mathrm{s}^{-1}$
in (b) and (c), respectively. More conservatively, the respective bounds would
be $\Gamma=$3.1\cdot 10^{2}\text{\,}\mathrm{s}^{-1}$$ and $8.3\cdot
10^{2}\text{\,}\mathrm{s}^{-1}$ at the confidence level $1-10^{-3}$ (roughly
corresponding to $3\sigma$ for normally distributed estimates), or
$\Gamma=$5.5\cdot 10^{2}\text{\,}\mathrm{s}^{-1}$$ and $1.1\cdot
10^{3}\text{\,}\mathrm{s}^{-1}$ at the confidence level $1-10^{-7}$ (roughly
$5\sigma$). We attribute the greater threshold value in the superposition
measurement to the omitted influence of environmental dephasing, which occurs
at a comparable rate of $1/T_{\phi}=1.0\cdot$10^{3}\text{\,}\mathrm{s}^{-1}$$,
as discussed in the main text. The remaining discrepancy between this rate and
the inferred $\Gamma$ can be explained by the different impact of diffusion
and pure dephasing on the oscillator state.

In order to obtain the macroscopicity values, we take the maximum of the ruled
out $\tau_{e}$-range over $\sigma_{q}$, which coincides with the local maximum
of the dimensionless diffusion rate at
$\hbar/\sigma_{q}=$0.5\text{\,}\mu\mathrm{m}$$, as depicted in Fig. S3.
Employing the same 95% confidence level as in previous assessments, this leads
to the macroscopicity values $\mu=11.3$ and $10.7$ for the two measured
states.

Figure S5: Bayesian inference of the diffusion rate. (a): Histogram of pixel
deviations $\Delta_{ij}(t)$ between recorded data and model Wigner function,
sampled from four time snapshots of the Fock state measurement. The histogram
is closely matched by a Gaussian fit (solid line) with standard deviation
$s=0.034$. (b): Posterior distribution for the diffusion rate $\Gamma$ for the
single Fock state data up at $t=10,20,40\,\mu$s. For comparison, the dashed
curve depicts Jeffreys’ prior, which has been updated by the posterior toward
significantly lower $\Gamma$-values. The use of Jeffreys’ prior leads to a
more conservative estimate of excluded $\Gamma$-values than, say, the often
used flat prior. The vertical line and the yellow part mark the upper five
percent quantile of the posterior. (c) Posterior distribution for the
superposition state data.

##  Supplementary Section D Macroscopicity estimates for other mechanical
resonators

We report here the parameters used to estimate the macroscopicities reachable
in the two other reported experiments with mechanical resonators, where Wigner
function negativities have been measured in the initial state [12, 13]. These
experiments involve resonators with acoustic modes of complicated geometry
[28, 29], so that it is not easy to derive simple analytical expressions as
for the Laguerre-Gaussian modes we considered. For this reason, we proceed
with making a conservative analysis, based on approximating such devices to
cuboids hosting a sinusoidal displacement field,
$\bm{u}(\bm{r})=\cos\left(\pi/h\right)\mathbf{e}_{x}$ with $h$ the
longitudinal extension of the mode. This yields an analytic expression for the
geometric factor (8) determining macrorealistic diffusion, which would
otherwise require numerical integration of the displacement field resulting
from a finite-element simulation of the actual device. Our approximation is
conservative in the sense that it likely overestimates the effective
oscillator mass, but we do not expect it to result in macroscopicities that
deviate much from the actual one (given that $\mu$ is in logarithmic scale).

The experiments in Refs. [12, 13] show measurements of Wigner function
negativities with the sake of demonstrating nonclassical state preparation,
but they do not monitor the disappearance of such quantum features with time.
We can therefore merely estimate the potential of these setups to perform a
test of macrorealist models from the limitations imposed by the reported
$T_{1}$-times. To this end, we simply assume that our data were obtained with
their resonators at the same _relative_ times $t/T_{1}$, which amounts to
rescaling our inferred $\Gamma$-values by the ratio between our $T_{1}$ time
and theirs.

Phononic crystal resonator Ref. [13]: We consider a phonon mode with frequency
$\omega=$2\pi\cdot 2\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, wavelength
$\lambda=$1\text{\,}\mu\mathrm{m}$$, and decay time
$T_{1}=$1\text{\,}\mu\mathrm{s}$$. The resonator itself is approximated by a
lithium niobate (density $4.65\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$)
cuboid of size $1\times 1\times 0.25\,$\text{\,}\mu\mathrm{m}$$ and total mass
$m=$1.16\cdot 10^{-15}\text{\,}\mathrm{k}\mathrm{g}$$, with an effective mass
of $m_{\text{eff}}=5.8\cdot$10^{-16}\text{\,}\mathrm{k}\mathrm{g}$$. We assume
a hypothetical experiment that excludes diffusion rates
$\Gamma>$13.7\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$, which results in
$\mu=9.0$ after maximizing with respect to $\sigma_{q}$.

Surface acoustic waves Ref. [12]: We consider a phonon mode with frequency
$\omega=$2\pi\cdot 4\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, wavelength
$\lambda=$1\text{\,}\mu\mathrm{m}$$, and decay time
$T_{1}=$150\text{\,}\mathrm{n}\mathrm{s}$$. The resonator is approximated by a
lithium niobate (density $4.65\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$)
cuboid of size [30] $75\times 50\times 1\,$\text{\,}\mu\mathrm{m}$$, with
total mass $m=$1.744\cdot 10^{-11}\text{\,}\mathrm{k}\mathrm{g}$$ and
effective mass
$m_{\text{eff}}=8.7\cdot$10^{-12}\text{\,}\mathrm{k}\mathrm{g}$$. The
hypothetical experiment would exclude
$\Gamma>$91.5\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$, resulting in
$\mu=8.6$.

##  Supplementary Section E Experimental sequences

Figure S6: Characterization of phonon coherence. Left: population decay
measurement used to estimate the relaxation time $T_{1}$. Right: phonon Ramsey
measurement used to estimate the dephasing time $T_{2}$.

###  E.1 Phonon mode characterization

To characterize the phonon mode under investigation, we perform two standard
measurements [31]. The first is a population decay (“ring-down”) measurement,
which gives the relaxation rate $1/T_{1}$. This is done by preparing the
phonon mode in Fock state $\ket{1}$, and then monitoring its population decay
over time, see Fig. (S6). Fitting the measured data with a decaying
exponential $e^{-t/T_{1}}$ plus a constant background gives us the time
parameter $T_{1}=$85.8{\pm 1.5}\text{\,}\mu\mathrm{s}$$, where the error is
one standard deviation given by the fitting function. In the assumed model
(6), the measured relaxation rate corresponds to
$1/T_{1}=2\Gamma+\gamma_{\downarrow}$.

The second measurement consists of a Ramsey sequence, which monitors dephasing
effects. This is done by preparing the phonon mode in state
$(\ket{0}+\ket{1})/2$, and then monitoring its phase evolution over time, see
Fig. (S6). Fitting the measured data to an exponentially decaying oscillation
gives us the time constant $T_{2}=$147.3{\pm 2.6}\text{\,}\mu\mathrm{s}$$,
where the error is one standard deviation given by the fitting function. As
this contrast reduction is partially due to the relaxation effect previously
investigated, it is convenient to identify the pure dephasing time as
$T_{\phi}=(1/T_{2}-1/(2T_{1}))^{-1}=$1.0{\pm
0.2}\text{\,}\mathrm{m}\mathrm{s}$$. This pure dephasing could be described by
a Lindblad generator of the form
$2\mathcal{D}[\mathsf{a}^{\dagger}\mathsf{a}]/T_{\phi}$.

We can also estimate the ground state population and the effective temperature
of the phonon mode, which can be used for a classical noise test of
macrorealistic collapse models. To this end, we tune the qubit frequency to be
$0.8\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ below the one of the phonon mode,
operating our system in the strongly dispersive regime [14]. In this regime,
the qubit frequency depends on the phonon mode occupation number, and
therefore on its effective temperature. In Fig. S7 we present a spectroscopic
measurement of the qubit frequency with the phonon mode in the stationary
state, showing a prominent peak at the bare qubit frequency. From the data it
is possible to note a slight asymmetry between the left and right tails of the
peak, which is highlighted by comparing to a Voigt profile fit (blue line).
This is expected at finite temperature, as any population in higher Fock
states will result in additional peaks on the left-hand side of the main one.
As the position of the peak for $\ket{1}$ is known, we fit again the data with
a sum of two Voigt profiles of fixed frequencies (red dashed line). From the
relative amplitude between the two peaks we estimate the thermal population of
the 1-phonon excited state to be $1.6\pm 0.2\,\%$. In our model (6) in the
main text, this small population corresponds to
$\Gamma/(2\Gamma+\gamma_{\downarrow})$ to a good approximation. Hence
$\Gamma\ll\gamma_{\downarrow}$, and we can safely identify the measured
relaxation rate with the model rate, $\gamma_{\downarrow}\approx 1/T_{1}$,
within the uncertainty of the fit.

Figure S7: Phonon mode thermometry. Spectroscopy of the qubit resonance
frequency, while coupled to the phonon mode in the strong dispersive regime.
In this regime, the qubit frequency depends on the phonon occupation number.
Observing an asymmetry in the measured curve thus indicates a non-zero phonon
number, that we associate to an effective temperature (see text for details).

###  E.2 State preparation and read-out

To acquire the data presented in Fig. 2 of the main text, we adopt the
sequence illustrated in Fig. S8. First, the phonon mode is prepared in state
$\ket{1}$ or $(\ket{0}+\ket{1})/\sqrt{2}$ by swapping the corresponding qubit
states $\ket{\uparrow}$ or $(\ket{\downarrow}+\ket{\uparrow})/\sqrt{2}$ with
the phonon ground state through the resonant Jaynes–Cummings interaction.
Then, the system evolves freely for a variable time interval $t$, after which
a Wigner function measurement is performed. The latter consists of a displaced
parity measurement, where a displacement of the phonon mode $D(\beta)$ is
followed by a parity measurement implemented through a strong dispersive
interaction between the qubit and the phonon [14].

The above sequence results in a single measurement of the parity at point
$-\beta$ in phase space. The value of the Wigner function is then obtained by
averaging $\sim 4\cdot 10^{3}$ of such parity measurements. Finally, to record
a two-dimensional Wigner function like the ones in Fig. 2, we repeat the
described procedure for many different values of $\beta$.

Figure S8: Wigner function measurement of the time-evolved state. The
experimental sequence used to record the data in Fig. 2 consists of the
following steps: 1) phonon state preparation, 2) wait time $t$, and 3)
displaced parity measurement. Further details about the parity measurement
sequence are given in Ref. [14].

##  Supplementary Section F Non-interferometric thermalization test

A straightforward and widely used approach to study collapse models are so-
called non-interferometric tests, which focus on monitoring the noise or the
energy changes in the studied system [21]. This puts bounds on modifications
to quantum mechanics without the need of preparing a quantum state. In our
experiment, we can do a similar analysis by monitoring the resonator
occupation number. The Wigner function evolution under relaxation at the rate
${\gamma_{\downarrow}}$ and macrorealistic diffusion predicts a steady-state
energy $E_{\mathrm{therm}}=\hbar\omega(1+2\Gamma/{\gamma_{\downarrow}})/2$.
Hence, by measuring the steady state phonon population, we impose bounds on
the hypothetical diffusion rate $\Gamma$.

To improve the bounds resulting from this analysis, one could also take into
account additional known decoherence processes and a finite temperature
characterized from independent measurements. Instead, we follow the most
conservative estimation: we attribute all diffusion effects to the
macrorealistic modification, thus overestimating its effect and
underestimating the resulting bounds. Then the equilibrium temperature
$T_{\mathrm{therm}}=\hbar\omega\Gamma/{\gamma_{\downarrow}}k_{B}$ is a result
of the competing relaxation rate ${\gamma_{\downarrow}}$ and macrorealistic
diffusion rate $\Gamma$. The decay rate ${\gamma_{\downarrow}}$ can be
reliably measured far from equilibrium by preparing a thermal state at
$T(t=0)\gg T(t=\infty)=T_{\mathrm{therm}}$ and verifying a constant decay
rate. This leaves only one free parameter $\Gamma$ to estimate, which then
gives non-interferometric exclusion bounds on collapse models based on the
observed amount of classical heating.

With the maximum dimensionless diffusion rate $\Gamma_{\mathrm{max}}=3.5\cdot
10^{13}/\tau_{e}$ and a thermal population of $1.6\%$, this analysis leads to
$\Gamma=0.016{\gamma_{\downarrow}}$. The measured $\gamma_{\downarrow}\approx
1/T_{1}$ then results in an exclusion of $\tau_{e}<1.9\cdot 10^{11}\,$s at the
length scale $\hbar/\sigma_{q}=5\cdot 10^{-7}\,$m. This is weaker than other
non-interferometric collapse tests performed with _ad hoc_ experiments [21],
or with large-scale endeavors like the LISA pathfinder experiment [32, 33],
but it could be further improved by carefully optimizing the HBAR device
geometry and the experimental sequences. Note however that no genuine quantum
feature is actually verified in a collapse test based on classical noise
bounds, meaning that it does not qualify for a benchmark based on the
macroscopicity of a superposition.

Finally, as a side remark, let us mention that a new approach for performing
non-interferometric tests of collapse models using bulk resonators has been
recently proposed in Ref. [23]. Our analytical results can readily be applied
to this proposal, leading to an exact expression of the diffusion rate. Given
the mode function $u(\bm{r})=\cos(\pi\ell x/L)\mathbf{e}_{x}$ hosted by a
homogeneous cylinder of radius $R$ and length $L$, we obtain the diffusion
rate

$\displaystyle\Gamma=$
$\displaystyle\lambda_{\mathrm{C}}\frac{x_{0}^{2}\mathrm{bar}{\varrho}^{2}\pi^{2}r_{\mathrm{C}}^{5}R^{2}}{2m_{e}^{2}L^{3}}\left[\frac{\exp(R^{2}/2r_{\mathrm{C}}^{2})-I_{0}(R^{2}/2r_{\mathrm{C}})-I_{1}(R^{2}/2r_{\mathrm{C}})}{\exp(R^{2}/2r_{\mathrm{C}}^{2})}\right]$
$\displaystyle\times\left[-\ell\pi
h(L/\sqrt{2}r_{\mathrm{C}},0,n)+\sqrt{2}\left((-1)^{n}e^{-L^{2}/4r_{\mathrm{C}}^{2}}-1\right)\sigma_{L}\left(\ell^{2}\pi^{2}-L^{2}/r_{\mathrm{C}}^{2}\right)+(-1)^{n}\ell\pi^{2}\mathrm{Re}\left\\{h(L/\sqrt{2}r_{\mathrm{C}},L/\sqrt{2}r_{\mathrm{C}},\ell)\right\\}\right]\;.$
(36)

For convenience, we have adopted the parameter notation of the CSL collapse
model here, $\lambda_{\mathrm{C}}=m_{0}^{2}/m_{e}^{2}\tau_{e}$ and
$r_{\mathrm{C}}=\hbar/\sqrt{2}\sigma_{q}$. Our formula must be multiplied by
two for direct comparison with the approximate momentum diffusion rate term
given by the authors,

$\displaystyle 2\Gamma\approx$
$\displaystyle\lambda_{\mathrm{C}}\frac{x_{0}^{2}r_{\mathrm{C}}^{3}\mathrm{bar}{\varrho}^{2}R^{2}}{2m_{0}^{2}}\left(\frac{64\sqrt{\pi}}{L/\ell}\right)\left[\frac{\exp(R^{2}/2r_{\mathrm{C}}^{2})-I_{0}(R^{2}/2r_{\mathrm{C}})-I_{1}(R^{2}/2r_{\mathrm{C}})}{2\exp(R^{2}/2r_{\mathrm{C}}^{2})}\right]$
$\displaystyle\times\int_{-\infty}^{\infty}\mathrm{d}a\,\exp\left(-r_{\mathrm{C}}^{2}\left(\frac{2\pi
a}{L/\ell}\right)^{2}\right)\frac{(-8+(8+a^{2}\pi^{2})\cos(a\pi/2))^{2}}{4a^{2}\pi^{2}}\frac{\sin^{2}(\ell\pi
a/2)}{\cos^{2}(\pi a/2)}.$ (37)

A comparison between both expressions is shown in Fig. S9 for the parameters
chosen in Ref. [23],
$\mathrm{bar}{\varrho}=$3210\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$$,
$R=$35\text{\,}\mu\mathrm{m}$$,
$\omega=2\pi\times$6.33\text{\,}\mathrm{H}\mathrm{z}$$, and for the two
settings at $\ell=1,L=$1.5\text{\,}\mu\mathrm{m}$$ and at
$\ell=40,L=$60\text{\,}\mu\mathrm{m}$$, showing good agreement.

Figure S9: Comparison between diffusion rates. Our results (red), compared to
the results of Ref. [23] (black). Solid curve is for a mode index of $\ell=1$,
while the dashed curve is for $\ell=40$. We chose to plot $2\Gamma=2\eta
x_{0}^{2}$ for a direct comparison.