# Testing spontaneous wavefunction collapse with quantum electromechanics:
Supplemental Material

Germain Tobar1,2, Stefan Forstner3, Arkady Federov1, and Warwick P. Bowen1
1Australian Centre for Engineered Quantum Systems, School of Mathematics and
Physics, The University of Queensland, Australia  2Department of Applied
Mathematics and Theoretical Physics, Centre for Mathematical Sciences,
Wilberforce Road, Cambridge, CB3 OWA, UK  3ICFO-Institut de Ciencies
Fotoniques, The Barcelona Institute of Science and Technology, Castelldefels
(Barcelona) 08860, Spain

##  I Heating rate from spontaneous collapse

###  I.1 CSL-Linblad term for an internal mechanical mode

Here we calculate the decoherence rate of a bulk acoustic wave resonator by
deriviing the CSL master equation for the internal mechanical mode of a
resonator. We provide the most accurate calculation of the CSL cross-section
of an internal mechanical mode of a resonator, and derive the generalisation
of the CSL cross-section for higher order internal mechanical modes.

In order to determine the CSL master equation for an internal mechanical mode,
we consider the CSL master equation for a body of N masses $m_{n}$ with
position operators $\vec{r}_{n}$ [1]

$\mathcal{L}_{\mathrm{csl}}\hat{\rho}=\frac{\lambda_{\mathrm{csl}}}{\pi^{3/2}r_{\mathrm{csl}}^{3}\mathrm{amu}^{2}}\int\mathrm{d}^{3}s\left(m(\vec{s})\hat{\rho}m(\vec{s})-\frac{1}{2}\left\\{\hat{\rho},m^{2}(\vec{s})\right\\}\right).$
(1)

Here, $\mathcal{L}_{\mathrm{CSL}}$ is the CSL-Lindblad operator, $\hat{\rho}$
is the system’s density matrix, $\mathrm{amu}$ is the atomic mass unit and
$\lambda_{c}$, $r_{c}$ are the CSL collapse rate and noise correlation length
respectively.

The Gaussian-smeared mass density operator for CSL is defined as [1, 2]

$m(\vec{s})=\sum_{n}m_{n}e^{-\frac{(\vec{s}-\vec{r}n)^{2}}{2r_{c}^{2}}}=\int_{V}\mathrm{~{}d}^{3}\vec{r}\rho(\vec{r})e^{-\frac{(\vec{s}-\vec{r})^{2}}{2r_{c}^{2}}},$
(2)

where $\vec{r}_{n}$ is the coordinate of the $n$-th mass $m_{n}$ and
$\rho(\vec{r})$ is the mass density. By taking a Fourier transform, the mass
density can be expressed as

$m(\vec{s})=\frac{r_{c}^{3}}{(2\pi)^{3/2}}\int\mathrm{d}^{3}\vec{k}e^{-\frac{1}{2}r_{c}^{2}k^{2}}\int_{V}\mathrm{~{}d}^{3}\vec{r}_{n}\rho\left(\vec{r}_{n}\right)e^{i\vec{k}\cdot\left(\vec{s}-\vec{r}_{n}\right)}.$
(3)

The position operator for an infinitesimal volume element of the continuous
structure is

$\vec{r}=\vec{r}_{\mathrm{CoM}}+\vec{r}^{(0)}+\Delta\vec{r},$ (4)

where $\Delta\vec{r}$ is the displacement of the volume element relative to
its equilibrium position, $\vec{r}^{(0)}$ is the equilibrium position of the
volume element relative to the cenre of mass and $\vec{r}_{\mathrm{CoM}}$ is
the centre-of-mass (CoM) operator. For an internal mechanical mode of a
resonator there is no oscillation of the system’s centre-of-mass requiring
meaning that $\vec{r}_{\mathrm{CoM}}=0$ [2]. Therefore, we must examine more
closely the effects of $\Delta\vec{r}$. In this work, we will model
$\Delta\vec{r}$ as an ideal breathing mode. A breathing mode is an expansion
of a structure along one or multiple axes with contraction along the other
axes [2]. The amount of contraction is characterised by its Poisson ratio [2].

####  I.1.1 1D breathing mode of a cylinder

We start by assuming that expansion and contraction of the breathing motion of
a cylindrical cavity occurs along a single direction - the vertical direction
of phonon polarisation. This assumption holds for materials with low Poisson
ratios such as silicon or quartz and simple resonator geometries [2]. We
consider a cylindrical cavity of height $\lambda/2$ and radius $R$ with
periodic expansion and contraction solely in the vertical $x$-direction with
no motion in the transverse directions. We define the variable $x_{s}$ to
quantify the displacement at each breathing surface. The position operator of
a volume element in a $1-\mathrm{D}$ breathing mode in cylindrical coordinates
is

$\vec{r}_{n}=\underbrace{x_{n}^{(0)}\hat{e}_{x}+y_{n}^{(0)}\hat{e}_{y}+z_{n}^{(0)}\hat{e}_{z}}_{\vec{r}_{n}^{(0)}}+\underbrace{x_{s}\sin\left(\frac{2\pi
x_{n}^{(0)}}{\lambda}\right)\hat{e}_{x}}_{\Delta\vec{r}_{n}},$ (5)

where the amplitude of motion of the volume element follows the sinusoidal
mode profile. We assume that the material is homogeneous, allowing us to
approximate the density of the material as a constant
$\rho\left(\vec{r}_{n}\right)\approx\rho_{0}$ and integrate over
$\vec{r}^{(0)}$ instead of $\vec{r}$. Therefore, the expression for the mass
density operator becomes

$\begin{array}[]{r}m(\vec{s})=\frac{\rho_{0}r_{c}^{3}}{(2\pi)^{3/2}}\int\mathrm{d}^{3}ke^{-\frac{1}{2}r_{c}^{2}k^{2}}e^{i\vec{k}\cdot\vec{s}}\cdot\int_{-\lambda/4}^{\lambda/4}e^{-ik_{x}\left(x^{(0)}+x_{s}\sin\left(\frac{2\pi
x_{n}^{(0)}}{\lambda}\right)\right)}\mathrm{d}x^{(0)}\\\
\cdot\int_{-R}^{R}\int_{-\sqrt{R^{2}-y^{2}}}^{\sqrt{R^{2}-y^{2}}}\exp\left(-ik_{y}y_{n}^{(0)}\right)\exp\left(-ik_{z}z_{n}^{(0)}\right)\mathrm{d}z_{n}^{(0)}\mathrm{d}y_{n}^{(0)}.\end{array}$
(6)

Here we have broken up the 3-D integral into a 1-D integral in the
longitudinal coordinate $x_{n}^{(0)}$, and a 2-D integral in the transverse
coordinates $y_{n}^{(0)}$ and $z_{n}^{(0)}$. The 2-D integral in transverse
coordinates can be more easily solved by converting to polar coordinates:

$\begin{split}\tilde{\rho}_{\perp}\left(k_{\perp}\right)=\tilde{\rho}_{\perp}\left(k_{y},k_{z}\right)&=\int_{-R}^{R}\int_{-\sqrt{R^{2}-y^{2}}}^{\sqrt{R^{2}-y^{2}}}\exp(-ik_{y}y_{n}^{(0)})\exp(-ik_{z}z_{n}^{(0)})\mathrm{d}z^{(0)}_{n}\mathrm{d}y^{(0)}_{n}\\\
&=\int_{0}^{2\pi}\int_{0}^{R}r\exp(-ik_{y}r\sin(\theta))\exp(-ik_{z}r\cos(\theta))\mathrm{dr}\mathrm{d\theta}\\\
&=2\pi\int_{0}^{R}r\mathrm{J}_{0}\left(r\sqrt{k_{y}^{2}+k_{z}^{2}}\right)\mathrm{dr}\\\
&=2\pi\int_{0}^{R}r\mathrm{J}_{0}\left(rk_{\perp}\right)\mathrm{dr}\\\
&=\frac{2\pi R}{k_{\perp}}\mathrm{J}_{1}(Rk_{\perp}).\end{split}$ (7)

Here, $\mathrm{J}_{\alpha}$ are Bessel functions of the first kind. After
making the substitution $u=\frac{2\pi}{\lambda}x_{n}^{(0)}$, the longitudinal
integral simplifies to

$\tilde{\rho}_{x}\left(k_{x}\right)=\left(\frac{\lambda}{2\pi}\right)\cdot\int_{-\pi/2}^{\pi/2}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}+x_{s}\sin(u)\right)}\mathrm{d}u.$ (8)

we approximate $\sin(x)$ as a quadratic for $0\leq x\leq\pi$:

$\sin(x)\approx-\frac{4}{\pi^{2}}x^{2}+\frac{4}{\pi}x,$ (9)

while for $-\pi\leq x\leq 0$:

$\sin(x)\approx\frac{4}{\pi^{2}}x^{2}+\frac{4}{\pi}x.$ (10)

After making these approximations, the longitudinal integral becomes

$\tilde{\rho}_{x}\left(k_{x}\right)=\left(\frac{\lambda}{2\pi}\right)\cdot\left(\int_{0}^{\pi/2}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}+x_{s}(-\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u)\right)}\mathrm{d}u+\int_{-\pi/2}^{0}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}+x_{s}(\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u)\right)}\mathrm{d}u\right).$
(11)

Now expressing the complex exponential
$e^{-ik_{x}\left(x_{s}(-\frac{4}{\pi}u^{2}+\frac{4}{\pi}u)\right)}$ in terms
of $\mathrm{sin}$ and $\mathrm{cos}$, taking the Taylor series of the
trigonometric funtions and eliminating all terms up to first order in
$k_{x}x_{s}$ we obtain:

$\exp\left(-ik_{x}\left(x_{s}(-\frac{4}{\pi}u^{2}+\frac{4}{\pi}u)\right)\right)\approx
1-ik_{x}x_{s}\left(-\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u\right).$ (12)

After a similar calculation for
$e^{-ik_{x}\left(x_{s}(\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u)\right)}$ we
obtain:

$\exp\left(-ik_{x}\left(x_{s}(\frac{4}{\pi}u^{2}+\frac{4}{\pi}u)\right)\right)\approx
1-ik_{x}x_{s}\left(\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u\right).$ (13)

We now substitute these expressions into the mass density operator:

$\begin{split}m(\vec{s})=\frac{\rho_{0}r_{c}^{3}}{(2\pi)^{3/2}}&\int\mathrm{d}^{3}ke^{-\frac{1}{2}r_{c}^{2}k^{2}}e^{i\vec{k}\cdot\vec{s}}\left(\frac{\lambda}{2\pi}\right)\\\
&\cdot(\int_{0}^{\pi/2}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}\right)}\left(1-ik_{x}x_{s}\left(-\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u\right)\right)\mathrm{d}u\\\
&+\int_{-\pi/2}^{0}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}\right)}\left(1-ik_{x}x_{s}\left(\frac{4}{\pi^{2}}u^{2}+\frac{4}{\pi}u\right)\right)\mathrm{d}u)\tilde{\rho}_{\perp}\left(k_{y},k_{z}\right).\end{split}$
(14)

Once substituted into the CSL master equation, only terms proportional to
$x_{s}$ are non-zero, allowing the mass density operator to be expressed as

$\begin{split}m(\vec{s})=\frac{\rho_{0}r_{c}^{3}}{(2\pi)^{3/2}}&\int\mathrm{d}^{3}ke^{-\frac{1}{2}r_{c}^{2}k^{2}}e^{i\vec{k}\cdot\vec{s}}\left(\frac{2\lambda}{\pi^{2}}\right)\\\
&\cdot(\int_{0}^{\pi/2}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}\right)}\left(-ik_{x}x_{s}\left(-\frac{1}{\pi}u^{2}+u\right)\right)\mathrm{d}u\\\
&+\int_{-\pi/2}^{0}e^{-ik_{x}\left(\frac{\lambda
u}{2\pi}\right)}\left(-ik_{x}x_{s}\left(\frac{1}{\pi}u^{2}+u\right)\right)\mathrm{d}u)\tilde{\rho}_{\perp}\left(k_{y},k_{z}\right).\end{split}$
(15)

We can now express the mass density operator as

$m(\vec{s})=\frac{r_{\mathrm{csl}}^{3}\rho_{0}}{(2\pi)^{3/2}}\int\mathrm{d}^{3}\vec{k}e^{-\frac{1}{2}r_{\mathrm{csl}}^{2}k^{2}}e^{i\vec{k}\cdot\vec{s}}\tilde{\rho}_{x}\left(k_{x}\right)\tilde{\rho}_{\perp}\left(k_{\perp}\right),$
(16)

substituting this expression into the CSL master equation, it becomes:

$\mathcal{L}_{\mathrm{csl}}\hat{\rho}\approx\frac{\lambda_{\mathrm{csl}}r_{\mathrm{csl}}^{3}\rho_{0}^{2}}{\pi^{3/2}\mathrm{amu}^{2}}\int\mathrm{d}^{3}ke^{-r_{\mathrm{csl}}^{2}k^{2}}\left|\tilde{\rho}_{\perp}\right|^{2}\left(\tilde{\rho}_{x}\hat{\rho}\tilde{\rho}_{x}-\frac{1}{2}\left\\{\hat{\rho},\left|\tilde{\rho}_{x}\right|^{2}\right\\}\right),$
(17)

here we have applied $\int
e^{i\left(\vec{k}-\vec{k}^{\prime}\right)\cdot\vec{s}}\mathrm{d}^{3}\vec{s}=(2\pi)^{3}\delta\left(\vec{k}-\vec{k}^{\prime}\right)$
to reduce the number of variables. Simplifying with the identity

$\tilde{\rho}_{x}\hat{\rho}\tilde{\rho}_{x}-\frac{1}{2}\left\\{\hat{\rho},\left|\tilde{\rho}_{x}\right|^{2}\right\\}=-\frac{1}{2}\left[\tilde{\rho}_{x},\left[\tilde{\rho}_{x},\hat{\rho}\right]\right],$
(18)

and expressing the commutator in terms of the position operator for the
displacement

$-\frac{1}{2}\left[\tilde{\rho}_{x},\left[\tilde{\rho}_{x},\hat{\rho}\right]\right]=|\tilde{\rho}^{\prime}_{x}\left(k_{x}\right)|^{2}\cdot(-1)\cdot\left[x_{s},\left[x_{s},\hat{\rho}\right]\right]$
(19)

where

$|\tilde{\rho}^{\prime}_{x}\left(k_{x}\right)|^{2}=\left(\frac{\lambda}{2\pi}\right)^{2}k_{x}^{2}\frac{\left(-8+\left(8+a^{2}\pi^{2}\right)\cos\left(\frac{a\pi}{2}\right)\right)^{2}}{4a^{6}\pi^{2}},$
(20)

and $a=k_{x}\lambda/2\pi$, we can finally express the CSL Master equation as

$\begin{split}\mathcal{L}_{\mathrm{csl}}\hat{\rho}&=-\frac{\lambda_{c}r_{c}^{3}\rho_{0}^{2}}{\pi^{3/2}\mathrm{amu}^{2}}\left(\frac{16}{\pi\lambda}\right)\int\exp(-r_{c}^{2}\left(\frac{2\pi
a}{\lambda}\right)^{2})\frac{\left(-8+\left(8+a^{2}\pi^{2}\right)\cos\left(\frac{a\pi}{2}\right)\right)^{2}}{4a^{4}\pi^{2}}\mathrm{d}a\\\
&\cdot\int\int
e^{-r_{\mathrm{csl}}^{2}(k_{x}^{2}+k_{y}^{2})}\frac{4\pi^{2}R^{2}}{(k_{y}^{2}+k_{z}^{2})}\mathrm{J}_{1}^{2}\left(R\sqrt{k_{y}^{2}+k_{z}^{2}}\right)\mathrm{d}k_{y}\mathrm{d}k_{z}(-1)\left[x_{s},\left[x_{s},\hat{\rho}\right]\right].\end{split}$
(21)

The 2-D integral in transverse coordinates can be more easily solved by
converting to polar coordinates [1]:

$\begin{split}&\iint
e^{-r_{\mathrm{csl}}^{2}\left(k_{\mathrm{x}}^{2}+k_{y}^{2}\right)}\frac{4\pi^{2}R^{2}}{\left(k_{y}^{2}+k_{z}^{2}\right)}\mathrm{J}_{1}^{2}\left(R\sqrt{k_{y}^{2}+k_{z}^{2}}\right)\mathrm{d}k_{y}\mathrm{d}k_{z}\\\
&=\int_{0}^{2\pi}\int_{0}^{\infty}e^{-r_{\mathrm{csl}}^{2}\left(k_{r}^{2}\right)}k_{r}\frac{4\pi^{2}R^{2}}{\left(k_{r}^{2}\right)}\mathrm{J}_{1}^{2}\left(Rk_{r}\right)\mathrm{d}k_{r}\mathrm{d}\theta\\\
&=8\pi^{3}R^{2}\int_{0}^{\infty}\frac{e^{-r_{\mathrm{csl}}^{2}\left(k_{r}^{2}\right)}}{k_{r}}\mathrm{J}_{1}^{2}\left(\mathrm{Rk}_{r}\right)\mathrm{d}k_{r}\\\
&=4\pi^{3}R^{2}\left[\frac{\exp(\frac{R^{2}}{2r_{\mathrm{CSL}}^{2}})-\mathrm{I}_{0}\left(\frac{R^{2}}{2r_{\mathrm{CSL}}^{2}}\right)-\mathrm{I}_{1}\left(\frac{R^{2}}{2r_{\mathrm{CSL}}^{2}}\right)}{2\exp(\frac{R^{2}}{2r_{\mathrm{CSL}}^{2}})}\right].\end{split}$
(22)

Here $\mathrm{I}_{n}$ are the modified Bessel functions of the first kind.
Therefore, the CLS-Lindblad operator is

$\begin{split}\mathcal{L}_{\mathrm{csl}}\hat{\rho}=-\lambda_{c}\eta_{\mathrm{new}}\left[x_{s},\left[x_{s},\hat{\rho}\right]\right].\end{split}$
(23)

where

$\begin{split}\eta_{\mathrm{new}}=\frac{r_{c}^{3}\rho_{0}^{2}R^{2}}{\mathrm{amu}^{2}}\left(\frac{64\pi^{1/2}}{\lambda}\right)&\cdot\left[\frac{\exp\left(\frac{R^{2}}{2r_{\mathrm{c}}^{2}}\right)-\mathrm{I}_{0}\left(\frac{R^{2}}{2r_{\mathrm{c}}^{2}}\right)-\mathrm{I}_{1}\left(\frac{R^{2}}{2r_{\mathrm{c}}^{2}}\right)}{2\exp\left(\frac{R^{2}}{2r_{\mathrm{c}}^{2}}\right)}\right]\\\
&\cdot\int_{-\infty}^{\infty}\exp(-r_{c}^{2}\left(\frac{2\pi
a}{\lambda}\right)^{2})\frac{\left(-8+\left(8+a^{2}\pi^{2}\right)\cos\left(\frac{a\pi}{2}\right)\right)^{2}}{4a^{4}\pi^{2}}\mathrm{d}a,\end{split}$
(24)

We now take the mode expansion for $x_{s}$ by substituting
$x_{s}=x_{0}\left(b+b^{{\dagger}}\right)$ where
$x_{0}=\sqrt{\frac{\hbar}{2m_{\mathrm{eff}}\omega}}$, such that
$\tilde{D}=\lambda_{c}x_{0}^{2}\eta_{\mathrm{new}}$ is the CSL-induced
decoherence rate. We can neglect terms of the form
$\hat{b}^{{\dagger}}\hat{b}^{{\dagger}}$ and $\hat{b}\hat{b}$ as long as the
mechanical oscillator’s resonance frequency $\omega$ is much larger than the
decoherence rate $\tilde{D}$ and the oscillator’s damping rate $\gamma_{r}$:
$\tilde{D},\gamma_{r}\ll\omega$. In this regime after taking the mode
expansion, the CSL-Lindblad term has the form
$-\lambda_{c}\eta_{\mathrm{new}}\left[\delta_{x},\left[\delta_{x},\hat{\rho}\right]\right]\approx
2\tilde{D}(\mathcal{D}[\hat{b}]+\mathcal{D}[\hat{b}^{{\dagger}}])\hat{\rho}$.
If the mechanical oscillator is in its ground state only the
$\mathcal{D}[\hat{b}^{{\dagger}}]$ term applies. The coefficient of the
$\mathcal{D}[\hat{b}^{{\dagger}}]$ superoperator quantifies the rate at which
phonons are generated in the resonator.

We now define the dimensionless the dimensionless CSL cross-section of a
cylindrical bulk acoustic wave resonator to be:

$\begin{split}D&=\frac{2\tilde{D}}{\lambda_{c}},\\\
&=2x_{0}^{2}\eta_{\mathrm{new}}.\end{split}$ (25)

For quartz with a density of $\rho_{0}\approx 2648\;\mathrm{kgm^{-3}}$ [3],
and a cylindrical resonator of radius $R\approx 35\mathrm{\mu m}$ with a
phonon wavelength of $\lambda\approx 9.5\times 10^{-7}\;\mathrm{\ m}$, which
corresponds to a mode frequency of $\omega/2\pi\approx 6.65\;\mathrm{GHz}$ and
assuming the expected collapse noise correlation length to be $r_{c}\approx
10^{-7}\;\mathrm{m}$, we obtain

$D\approx 7.5\times 10^{5}$ (26)

at $r_{c}=10^{-7}\;\mathrm{m}$.

Eq. (24) shows that the CSL cross-section is maximised for phonon wavelengths
$\lambda\approx 6r_{c}$. Therefore, a phonon wavelength of $\lambda\approx
9.5\times 10^{-7}\;\mathrm{\ m}$ is chosen to enhance the sensitivity to
collapse-induced heating, while ensuring that the mode frequency is within the
tunable frequency range of the qubit.

####  I.1.2 Generalising to multiple anti-nodes

Here we generalise the calculation of the CSL cross-section to higher-order
breathing modes.

We consider the resonator to contain $l$ breathing modes side-by-side in the
longitudinal direction, where each breathing mode has an equilibrium height of
$\lambda/2$. Therefore, there are $l$ breathing modes between $x=0$, and
$x=h$. The direct result is that the effect of generalising to $l$
longitudinal breathing modes is that there are $l$ of the masses which make up
the mass density operator. In this case the mass density operator becomes:

$\begin{split}m_{\mathrm{new}}(\vec{s})&=\sum_{i=0}^{i=l}\frac{r_{\mathrm{csl}}^{3}}{(2\pi)^{3/2}}\int\mathrm{d}^{3}\vec{k}\int_{(i-\frac{1}{2})\frac{\lambda}{2}}^{(i+\frac{1}{2})\frac{\lambda}{2}}\exp(-ik_{x}(x_{n}^{(0)}-i\frac{\lambda}{2}+x_{s}\sin(\frac{2\pi(x_{n}^{(0)}-i\frac{\lambda}{2})}{\lambda})))\mathrm{d}x_{n}^{(0)}\\\
&\cdot\int_{-R}^{R}\int_{-\sqrt{R^{2}-y^{2}}}^{\sqrt{R^{2}-y^{2}}}\exp(-ik_{y}y_{n}^{(0)})\exp(-ik_{z}z_{n}^{(0)})\mathrm{d}z^{(0)}_{n}\mathrm{d}y^{(0)}_{n}e^{-\frac{1}{2}r_{\mathrm{csl}}^{2}k^{2}}\\\
&=l\frac{r_{\mathrm{csl}}^{3}}{(2\pi)^{3/2}}\int\mathrm{d}^{3}\vec{k}e^{-r_{\mathrm{csl}}^{2}k^{2}}\int_{V}\mathrm{d}^{3}\vec{r}_{n}\rho\left(\vec{r}_{n}\right)e^{+i\vec{k}\left(\vec{s}-\vec{r}_{n}\right)}.\end{split}$
(27)

Each integral in the sum simply has the position operator’s equilibrium
position offset $\lambda/2$ from its nearest neighbour. This allows the
simplification of the expression in the first line to the expression in the
final line by recognising that it is the sum of $l$ identical integrals.
Therefore, increasing to multiple anti-nodes causes the mass density operator
to pick up a factor of the number of anti-nodes $l$.

The increase in the mass density operator by a factor of $l$ causes the total
CSL-Lindblad term to pick up a factor of $l^{2}$:

$\begin{split}\mathcal{L}^{\mathrm{total}}_{\mathrm{csl}}\hat{\rho}&=\frac{\lambda_{\mathrm{csl}}}{\pi^{3/2}r_{\mathrm{csl}}^{3}\mathrm{amu}^{2}}\int\mathrm{d}^{3}s\left(m_{\mathrm{new}}(\vec{s})\hat{\rho}m_{\mathrm{new}}(\vec{s})-\frac{1}{2}\left\\{\hat{\rho},m_{\mathrm{new}}^{2}(\vec{s})\right\\}\right),\\\
&=l^{2}\frac{\lambda_{\mathrm{csl}}}{\pi^{3/2}r_{\mathrm{csl}}^{3}\mathrm{amu}^{2}}\int\mathrm{d}^{3}s\left(m_{\mathrm{new}}(\vec{s})\hat{\rho}m_{\mathrm{new}}(\vec{s})-\frac{1}{2}\left\\{\hat{\rho},m_{\mathrm{new}}^{2}(\vec{s})\right\\}\right).\end{split}$
(28)

Using this modified result in the derivation of the CSL master equation, it is
easy to see that this results in the coefficient of the CSL-Lindblad term
displayed in Eq. (24) being scaled by a factor of $l^{2}$ due to the presence
of additional anti-nodes in the resonator:

$\begin{split}\eta^{\prime}=l^{2}\eta_{\mathrm{new}}.\end{split}$ (29)

However, the zero-point motion of the resonator is also increased due to the
presence of additional nodes. For the case of a single anti-node resonator in
a box of height $\lambda/2$, the effective mass is
$m_{\mathrm{eff}}=\frac{\lambda}{4}\rho d^{2}$. However, increasing the height
of the box to $l\frac{\lambda}{2}$, the effective mass becomes
$m_{\mathrm{eff}}^{\prime}=l\frac{\lambda}{4}\rho d^{2}=lm_{\mathrm{eff}}$.
Therefore the decrease in the zero-point motion for this larger effective mass
reduces the factor of $l^{2}$ in $\eta_{\mathrm{new}}$ to simply be a factor
of $l$:

$(x_{0}^{\prime})^{2}=\frac{\hbar}{m_{\mathrm{eff}^{\prime}}\omega}=\frac{\hbar}{lm_{\mathrm{eff}}\omega}=\frac{1}{l}x_{0}^{2}.$
(30)

Therefore, the CSL cross-section scales linearly with the number of anti-nodes
($l$):

$D=2(x_{0}^{\prime})^{2}\eta^{\prime}=2lx_{0}^{2}\eta_{\mathrm{new}}.$ (31)

This result provides the first derivation of the CSL cross-section of a non-
fundamental internal mechanical mode of a resonator. The CSL cross-section
demonstrates a mode number dependent scaling of the CSL-induced collapse rate
at fixed acoustic frequency. The mode number dependent scaling is a direct
result of the CSL-Lindblad term’s dependence on the system’s mass. Increasing
the number of anti-nodes by a factor of $l$ while simultaneously increasing
the height of the resonator to keep the frequency constant increases the
effective mass of the cavity that is susceptible to mechanical motion on the
order of the oscillator’s zero point motion by a factor of $l$. This explains
the increase in the CSL-induced heating rate as an impact of the increased
amount of mass susceptible to mechanical motion from coupling to the collapse
noise field.

The estimate for the CSL cross-section of $l=503$ breathing modes confined to
a cylindrical box can be found by simply scaling the previously calculated
result for $D$ by the longitudinal mode number $l$. The result is

$D_{l=503}\approx 3.8\times 10^{8}.$ (32)

For the lower bounds to the collapse rate postulated by Adler [4] and Bassi et
al. [5], our calculated CSL cross-section gives collapse-induced heating rates
of $\dot{n}_{c}=\lambda_{c}D\approx 3.8\times 10^{\pm 2}\;s^{-1}$ and
$\dot{n}_{c}=\lambda_{c}D\approx 3.8\times 10^{-2\pm 2}\;s^{-1}$,
respectively.

##  II Measurement Scheme

As shown in Figure 1, after the first swap operation and dispersive
measurement, the qubit is then tuned on-resonance with the second mechanical
resonator on a timescale fast compared to the swap operation period. This is
followed by another swap operation and dispersive measurement. This process is
repeated for all mechanical modes in the array. Once all mechanical modes have
had one interaction with the qubit, the qubit is then decoupled from all
resonators for $T\sim\gamma_{r}^{-1}\;\mathrm{s}$ in order to allow the
resonators to cool and collapse-induced heating to re-accumulate. This process
is repeated to build up statistics. A given CSL collapse rate can be ruled out
if it predicts a higher average phonon occupation than is observed for the
relevant correlation length.

Figure 1:  Experimental protocol. Swap operations are sequentially performed
between the qubit and an array of mechanical resonances, with mechanical
resonator $j$ having an electomechanical coupling rate $g_{j}$. The qubit
state is measured after each swap via interaction with a microwave cavity with
coupling rate $g$. The first few measurements, both at the beginning of the
protocol and after a detection event (red lines) are discarded, with the
associated swap operations instead used to initialize the qubit to near its
ground state. After all mechanical resonators have been measured, the qubit is
decoupled for a time of $\sim\gamma_{r}^{-1}\;\mathrm{s}$ in order to allow
the resonators to cool and collapse-induced heating to accumulate. The
protocol is then repeated until sufficient measurement statistics have been
obtained.

##  III Interaction Hamiltonian and coupling strength

The interaction energy between the qubit and the BAW resonator is given by

$H=-\int\sigma(\vec{x})s(\vec{x})d\mathrm{V}.$ (33)

Here, $\sigma(\vec{x})$ is the stress generated by the piezoelectric
transducer, $E(\vec{x})$ is the qubit’s electric field. Now, the stress
generated by the AIN is,

$\sigma(\vec{x})=c_{33}d_{33}(\vec{x})E(\vec{x}).$ (34)

Here $d_{33}(\vec{x})$ is the dominant piezoelectric tensor component
(component pointing in the longitudinal direction), $c_{33}$ is the dominant
stiffness tensor component of the resonator. We approximate that only the
dominant components of the stiffness and piezoelectric tensors contribute to
the interaction energy. Now, we quantise the qubit mode as
$E(\vec{x})(\hat{a}+\hat{a}^{{\dagger}})$, and the phonon mode as
$s(\vec{x})(\hat{b}+\hat{b}^{{\dagger}})$, we obtain:

$H=-c_{33}d_{33}\int
E(\vec{x})(\hat{a}+\hat{a}^{{\dagger}})s(\vec{x})(\hat{b}+\hat{b}^{{\dagger}})dV.$
(35)

We have approximated $d_{33}(\vec{x})$ as a constant inside the transducer
disk and zero everywhere else.

Assuming longitudinally polarised thickness modes confined below the PZT
transducer, the strain profile of the phonon mode is

$s_{l,m}(\vec{x})=\beta_{l,m}\sin\left(\frac{l\pi
z}{h}\right)s_{\perp}\left(r,\theta\right),$ (36)

where $s_{\perp}\left(r,\theta\right)$ is the phonon’s transverse mode profile
and $\beta_{l,m}$ is a normalisation constant so that the total energy of the
mode equal $\hbar\omega_{l,m}$. Here the mode frequency is given by

$\omega_{l,m}=\sqrt{\left(\frac{l\pi}{h}\right)^{2}v_{l}^{2}+\left(\frac{2j_{0,m}}{d}\right)^{2}v_{t}^{2}},$
(37)

and the normalisation factor is

$\beta_{l,m}=\sqrt{\frac{\hbar\omega_{l,m}}{\pi hc_{33}\int
J_{0}^{2}\left(\frac{2j_{0,m}r}{d}\right)rdr}}.$ (38)

We now assume that the electric field is the constant $E(\vec{x})=E_{0}$
throughout the transducer, and then apply the rotating wave approximation at
the fundamental transverse mode $\omega_{l,0}$ and equate the interaction
energy to the Jaynes-Cummings Hamiltonian $H_{\mathrm{int}}~{}=~{}-\hbar
g(\hat{a}\hat{b}^{{\dagger}}+\hat{a}^{{\dagger}}\hat{b})$, where $g$ is the
coupling strength between the qubit and the BAW resonator. This gives us the
following expression for the coupling strength:

$\begin{split}\hbar
g&=c_{33}d_{33}E_{0}\beta_{l,m}\int\sin\left(\frac{l\pi}{h}z\right)s_{\perp}\left(r,\theta\right)dV\\\
&=c_{33}d_{33}E_{0}\beta_{l,m}\int_{0}^{h}\sin\left(\frac{l\pi}{h}z\right)\mathrm{d}z\\\
&\cdot\int_{0}^{2\pi}\int_{0}^{d_{e}/2}rs_{\perp}\left(r,\theta\right)\mathrm{d}r\mathrm{d}\theta.\end{split}$
(39)

Here, $d_{e}$ is the diameter of the qubit electrode which can be chosen to
selectively couple to the fundamental transverse mode. The expression for the
coupling strength in Eq. (39) is exactly true for a cylindrical flat-cap
cavity geometry and it is an excellent approximation for a plano-convex cap
geometry as long as the radius of curvature is large compared to the cavity
height $R\gg h$, where $h$ height of the cylindrical cavity.

Henceforth we will assume that the plano-convex radius of curvature is large
such that we can approximate the transverse modes as Bessel functions. Using
the Bessel function solutions for the cylindrical flat-cap cavity geometry [6]
$s_{\perp}\left(r,\theta\right)=J_{0}\left(\frac{2j_{0,m}r}{d}\right)$, we
obtain the following expression for the coupling strength:

$\hbar
g_{l,m}=2\left(\frac{2h}{l}\right)c_{33}d_{33}E_{0}\beta_{l,m}\int_{0}^{d_{e}/2}J_{0}\left(\frac{2j_{0,m}r}{d}\right)r\mathrm{dr}.$
(40)

To take a specific example, targeting a correlation length on the order of
$r_{c}\approx 10^{-7}\;\mathrm{m}$, we consider a BAW resonator internal mode
with a frequency of $\omega/2\pi=6.65\;\mathrm{GHz}$ as have been previously
coupled to a qubit in Ref. [6]. In Ref. [6], the $l=503$ longitudinal and
$m=0$ transverse mode of a 420 $\mathrm{\mu m}$ thick sapphire BAW resonator
was used. Here, we propose to use quartz instead of sapphire due to the ultra-
high Q factors observed in quartz BAW resonators [7, 3]. We use parameters
based on Ref. [6] by proposing to couple to the $l=503$ longitudinal and $m=0$
transverse mode, however a thinner 240 $\mu m$ thick BAW resonator is required
to obtain the same resonance frequency as in Ref. [6] due to the lower density
of quartz, which corresponds to a phonon wavelength of $\lambda\approx
9.5\times 10^{-7}~{}\mathrm{m}$. As discussed earlier, in a BAW resonator the
CSL cross-section is maximised for phonon wavelengths $\lambda\approx 6r_{c}$,
such that any resonator with $\lambda\sim r_{c}$ has a large CSL cross-section
at $r_{c}$, however phonon wavelengths that are slightly larger than $r_{c}$
are more ideal, allowing our choice of phonon wavelength to target the
correlation length $r_{c}\approx 10^{-7}\;\mathrm{m}$, while providing a
resonance frequency within the tunable frequency range of the qubit.

In contrast to recent experiments which used AIN as the transducer material
[6, 8], we propose to use PZT. Reasonably conservative estimates of the
piezoelectricity of single crystal PZT (${d_{33}\approx 80\;\mathrm{pm}/V}$)
[9] offer a significant improvement on the piezoelectricity of AIN
(${d_{33}\approx 5\;\mathrm{pm}/V}$), allowing us to estimate an improved
coupling strength in our proposed experiment.

In Fig. 2 we show coupling strengh as a function of mode frequency
$\omega_{l,m}$ for a longitudinal mode number $l=503$ and the lowest order
transverse mode numbers of $m=0$ to $m=3$ for a cylindrical geometry bulk
acoustic wave resonator. We use $d=70\;\mathrm{\mu m}$ for the cavity
diameter, $d_{e}=55\mathrm{\mu m}$ for the electrode diameter,
$E_{0}=2.8\times 10^{-2}\;\mathrm{V/m}$, $c_{33}=106\;\mathrm{GPa}$ and a
cavity vertical height of $h=240\;\mu\mathrm{m}$. With these parameters, we
find the fundamental transverse mode at $l=503$ has frequency
$\omega_{l,0}/2\pi=6.65\;\mathrm{GHz}$ and its on-resonance coupling strength
with the qubit to be $g\approx 2\pi\times 3.05\;\mathrm{MHz}$ for a quartz BAW
resonator. This coupling strength is much larger than typical qubit and
mechanical resonator decay rates [6, 7], allowing the coupled system to
operate in the regime of strong coupling and therefore efficient swap
operations to be performed.

Figure 2:  Qubit-BAW resonator coupling strengths for the $l=503$ longitudinal
mode for a quartz resonator of height $h=240\;\mu\mathrm{m}$. The modes go to
zero at the boundary of the cavity of diameter $d=70\;\mathrm{\mu m}$. The
qubit electrode diameter $d_{e}=55\;\mathrm{\mu m}$ is slightly smaller to
create selective coupling to the fundamental transverse mode.

We can now model the idealised coupled dynamics In the rotating frame the
standard Jaynes-Cummings Hamiltonian with qubit coupling strengths $g_{l,m}$
and frequencies of the discrete phonon modes $\omega_{l,m}$. Therefore, the
Hamiltonian describing the coupled system is

$\begin{split}H=&\delta_{q}\sigma_{+}\sigma_{-}+\sum_{l,m}\delta_{l,m}b_{l,m}^{\dagger}b_{l,m}\\\
&+g_{l,m}\left(\sigma_{+}b_{l,m}+\sigma_{-}b_{l,m}^{\dagger}\right).\end{split}$
(41)

Here $\delta_{q}$ and $\delta_{l,m}$ is the detuning from the frequency
$\omega_{l,0}$ of the qubit and $l,m$ resonator mode respectively, $g_{l,m}$
is the coupling strength to the $l,m$ resonator mode, $\sigma_{+}$ and
$\sigma_{-}$ are the raising and lowering operators for the qubit and
$b_{l,m}^{\dagger},$ $b_{l,m}$ are the creation and annihilation operators
respectively for the $l,m$ resonator mode.

It is sufficient to only consider the $m=0$ to $m=3$ transverse modes and
adjacent longitudinal modes ($l=502,m=0$) and ($l=504,m=0$) in eq. (41) since
the coupling strength to other modes is negligible.

In the experimentally relevant regime where the mechanical resonator has a
high Q-factor (${\omega_{l,m}\gg\gamma_{r},\gamma_{q}}$), the dynamics of the
system are described with the phenomenological Linblad master equation [10]

$\begin{split}\dot{\varrho}&=\frac{-i}{\hbar}[H,\varrho]+\gamma_{q}\mathcal{D}\left[\sigma_{-}\right]\varrho\\\
&+(\Gamma_{QP}+\Gamma_{P})\mathcal{D}\left[\sigma_{+}\right]\varrho+\frac{\gamma_{\varphi}}{2}\mathcal{D}\left[\sigma_{z}\right]\varrho\\\
&+\sum_{l,m}\left[\gamma_{r}(1+\bar{n}_{\mathrm{th}})\mathcal{D}[b_{l,m}]\varrho+(\gamma_{r}\bar{n}_{\mathrm{th}}+\dot{n}_{c})\mathcal{D}[b^{{\dagger}}_{l,m}]\varrho\right],\end{split}$
(42)

which accounts for all known sources of decoherence.
$H=\sum_{l,m}g_{l,m}\left(\sigma_{+}b_{l,m}+\sigma_{-}b_{l,m}^{\dagger}\right)$
is the standard Jaynes-Cummings Hamiltonian (see App. III), $\varrho$ is the
density matrix,
$\bar{n}_{\mathrm{th}}=\left(e^{\hbar\omega_{l,m}/k_{B}T}-1\right)^{-1}$ is
the mechanical mean thermal occupancy (where $k_{B}$ is the Boltzmann constant
and $T$ is the temperature of the resonator), $\mathcal{D}$ is the dissipating
superoperator
$\mathcal{D}[L]\hat{\rho}=L\hat{\rho}L^{\dagger}-\frac{1}{2}\left(L^{\dagger}L\hat{\rho}+\hat{\rho}L^{\dagger}A\right)$.
The subscripts l,m refer to the $l,m$ mechanical modes, which are for
simplicity, assumed to have equal dissipation rates rates $\gamma_{r}$.
$\Gamma_{QP}/2\pi=30\;\mathrm{Hz}$ and $\Gamma_{P}/2\pi=0.5\;\mathrm{Hz}$ are
the heating rates in the qubit due to quasiparticle poisoning and the Purcell
effect respectively, while $\gamma_{\varphi}/2\pi=0.3\;\mathrm{MHz}$ is the
qubit dephasing rate [11] and $\gamma_{q}/2\pi=27\;\mathrm{kHz}$ is the qubit
dissipation rate [6].

In order to solve this master equation, we assume a qubit decay rate of
$\gamma_{q}/2\pi\approx 27\;\mathrm{kHz}$, a mechanical decay rate of
$\gamma_{r}/2\pi\approx 300\;\mathrm{Hz}$ and a dephasing rate of
$\gamma_{\varphi}/2\pi=0.3\;\mathrm{MHz}$, consistent with existing
experiments [11]. We assume the joint system is initially in the state
$\ket{g}\ket{n_{l=503,m=0}=1}$, where $\ket{n_{l,m}=1}$ denotes the $l,m$ mode
in the $n=1$ Fock state and all other modes in the ground state. Using our
previously calculated coupling strength of $g\approx 2\pi\times
3.05\;\mathrm{MHz}$, as shown in Fig. 3, we find the system has an efficiency
of $\eta_{\mathrm{swap}}\approx 0.8$ after a single swap operation, with
evidence of minor off-resonance coupling to adjacent modes.

Figure 3:  Occupancy probabilities for the qubit (red), BAW resonator $l=503$,
$m=0$ mode (blue), and off-resonant BAW modes (green). The x-axis is
normalised to the coupling strength $g$.

##  IV Gaussian measurement distribution and state discrimination error

During dispersive readout, quantum noise on the microwave field and dispersive
measurement inefficiency leads to a Gaussian statistical distribution for the
measured phase for each qubit state, which is known as state discrimination
error. The centre of the Gaussian for the ground and excited state is $x_{0}$
and $x_{1}$ respectively with each having the same standard deviation
$\sigma$. The discrimination point
$\left(x_{d}=x_{0}+\frac{(x_{1}-x_{0})}{a}\right)$ is the threshold phase that
separates measurements of the qubit state as $\ket{e}$ from those of
$\ket{g}$. Typically $a=2$ is chosen such that the probability of a false
measurement of the qubit in the $\ket{e}$ state (false positives) is equal to
the probability of a false measurement of the qubit in the $\ket{g}$ state
(false negatives) [12]. In this case the state discrimination error is
$\epsilon_{\text{sde
}}=\frac{1}{2}\operatorname{erfc}\left[\frac{\left|x_{0}-x_{1}\right|}{2\sqrt{2\sigma^{2}}}\right]$
and we can identify the signal to noise ratio as
$\mathrm{SNR}~{}=~{}\frac{\left(x_{0}-x_{1}\right)^{2}}{2\sigma^{2}}$ such
that $\epsilon_{\text{sde }}$ can be used as a measure of the $\mathrm{SNR}$.
As mentioned in the main text, we require false positives to be suppressed as
much as possible. This can be done by shifting the discrimination point as
shown in Fig. 4, such that $a<2$.

Figure 4:  Gaussian measurement distributions in the IQ plane for the qubit in
the ground state (blue) and excited state (red). The discrimination line (Fig.
4) distinguishes points on the IQ plane that are assigned to be a measurement
of the qubit in the $\ket{g}$ state from points that are assigned to be a
measurement of the qubit in the $\ket{e}$ state. By moving the discrimination
line upwards, we can suppress the probability of falsely identifying the qubit
as being in the excited state, at the expensive of also suppressing the
probability that we correctly identify the qubit in the excited state.

Shifting the discrimination point from the traditional choice of $a=2$ leads
to the following modified expressions for the probability of detecting a false
positive and true positive:

$\begin{split}P(\mathrm{False}\;\ket{e})&=\frac{1}{2}\operatorname{erfc}\left[\frac{\sqrt{\mathrm{SNR}}}{a}\right]\\\
P(\mathrm{True}\;\ket{e})&=\frac{1}{2}\operatorname{erfc}\left[\frac{\sqrt{\mathrm{SNR}}(1-a)}{a}\right].\end{split}$
(43)

In order to distinguish collapse-induced heating above state discrimination
error, we need to suppress $P(\mathrm{False}\;\ket{e})$ while maintaining a
reasonably high $P(\mathrm{True}\;\ket{e})$. Fig. 5 displays compares
$P(\mathrm{False}\;\ket{e})$ and $P(\mathrm{True}\;\ket{e})$ for a large range
of the applicable $(a,\mathrm{SNR})$ parameter space. We find that there is a
regime in the parameter space of $(a,\mathrm{SNR})$ for which it is possible
to exponentially suppress the error in false positives, while maintaining a
reasonable $P(\mathrm{True}\;\ket{e})$. For $a\geq 0.8$ and $\mathrm{SNR}\geq
10$ we find that $P(\mathrm{False}\;\ket{e})$ can be suppressed below
$P_{\mathrm{CSL}}\approx 10^{-6}$, while maintaining
$P(\mathrm{True}\;\ket{e})\gtrapprox 0.16$.

Figure 5:  (a) Contour plot of the probability of a false positive for the
parameter space of $\mathrm{SNR}$ and $a$. (b) Contour plot of the probability
of successfully distinguishing a true positive for the parameter space of
$\mathrm{SNR}$ and $a$. The red line displayed in the parameter space
corresponds to the values for $a$ and $\mathrm{SNR}$ for which the false
positive probability is suppressed to $10^{-6}$. Crucially there is a region
of the $\mathrm{SNR}$ and $a$ parameter space where the probability of a false
positive is suppressed below $P_{\mathrm{CSL}}\approx 10^{-6}$ while the
probability of a true positive still remains sufficiently large
$P(\mathrm{True}\;\ket{1})\sim 1$.

Ref. [13] reports the demonstration of rapid ($\sim 50\;\mathrm{ns}$) high SNR
single shot dispersive readout of a transmon qubit. Using similar experimental
parameters to those we assumed in the main text, a $\tau=56\;\mathrm{ns}$
measurement of the qubit state resulted in a state discrimination error of
only $\epsilon_{\mathrm{SDE}}\approx 0.002$. This corresponds to a relatively
high signal to noise ratio of $\mathrm{SNR}\approx 17$. Therefore, rapid
single shot dispersive readout that sufficiently suppresses
$P(\mathrm{False}\;\ket{e})$ below the collapse-induced occupation number
($P_{\mathrm{CSL}}\approx 10^{-6}$) is feasible with state-of-the-art quantum
technology.

##  V Heating due to quasiparticles and the Purcell effect

Recent measurements found a residual excited state population due to
quasiparticle excitations in a transmon qubit of $P(\ket{e})\approx 3\times
10^{-3}$ [14]. From their qubit decay rate of $\gamma_{q}/2\pi\approx
100\;\mathrm{kHz}$, we estimate the quasiparticle induced excitation rate to
be $\Gamma_{\mathrm{QP}}/2\pi\approx 300\;\mathrm{Hz}$. As discussed in the
main text, we expect current technology to mitigate quasiparticles should be
able to reduce this to $\Gamma_{\mathrm{QP}}/2\pi\approx 300\;\mathrm{Hz}$.

The Purcell effect is the enhancement of the qubit’s spontaneous excitation
rate due to the coupling of the qubit to the microwave resonator which is in
turn coupled to a transmission line [15]. The coupling to the transmission
line produces spontaneous excitation of the qubit at the rate [15]

$\Gamma_{P}\approx\frac{\kappa
g^{2}}{16\Delta^{2}}\left(\frac{\bar{n}}{n_{\mathrm{crit}}}\right)^{2}.$ (44)

For the experimental parameters we have proposed in previous sections, we
estimate the Purcell rate to be $\Gamma_{P}/2\pi\approx 26\;\mathrm{Hz}$. A
Purcell filter can be employed to reduce this rate by a further factor of
$\sim 50$ to $\Gamma_{P}/2\pi\approx 0.5\;\mathrm{Hz}$ [16]. A swap operation
will switch the spurious excited state population of the qubit into the
resonator and switch the CSL-induced occupation number of the resonator into
the qubit excited state population.

Here we solve the master equation displayed in Eq. (42) with the addition of
the term
$\left(\Gamma_{QP}+\Gamma_{P}\right)\mathcal{D}\left[\sigma_{+}\right]\hat{\rho}$
to account for quasiparticle and Purcell induced heating of the qubit.

Consider the qubit to be in equilibrium with quasiparticle poisoning induced
excitation rate $\Gamma_{\mathrm{QP}}$, decay rate
$\gamma_{\mathrm{\mathrm{QP}}}$, Purcell induced heating rate $\Gamma_{P}$,
decay rate $\gamma_{P}$ and thermally induced excitation rate
$\dot{n}_{\mathrm{th},q}$ and decay rate $\gamma_{q}$. In this case, the
equilibrium qubit occupancy $P(\ket{e})$, relates the sum of the excitation
and decay rates according to

$\Gamma_{\mathrm{QP}}+\Gamma_{P}+\dot{n}_{\mathrm{th},q}=P(\ket{e})\left(\gamma_{\mathrm{QP}}+\gamma_{\mathrm{P}}+\gamma_{q}\
\right).$ (45)

As discuss in the main text, $\dot{n}_{\mathrm{th},q}$ is suppressed at milli-
Kelvin temperatures. Furthermore the dominant source of decay is $\gamma_{q}$,
especially with the addition of a quasiparticle trap and Purcell filter to
suppress $\gamma_{\mathrm{QP}}$ and $\gamma_{\mathrm{P}}$. Therefore,
$\dot{n}_{\mathrm{th},q}$, $\gamma_{\mathrm{QP}}$ and $\gamma_{\mathrm{P}}$
can be neglected, and using $(\Gamma_{\mathrm{QP}}+\Gamma_{P})/2\pi\approx
30\;\mathrm{Hz}$ and again assuming a qubit decay rate of
$\gamma_{q}/2\pi=27\;\mathrm{kHz}$, the spurious excited state population of
the qubit is

$P(\ket{e})\approx\frac{\Gamma_{\mathrm{QP}}+\Gamma_{P}}{\gamma_{q}}\approx
1\times 10^{-3}.$ (46)

Therefore, the initial state of the joint BAW resonator-qubit coupled system
is

$\begin{split}\hat{\rho}(0)&\approx 1\times
10^{-3}\ket{e}\ket{\mathrm{vac}}\bra{e}\bra{\mathrm{vac}}\\\ &+(1-1\times
10^{-3})\ket{g}\ket{\mathrm{vac}}\bra{g}\bra{\mathrm{vac}}.\end{split}$ (47)

Here, we have assumed that there are is no occupancy of the $n=1$ Fock state
in the resonator due to CSL, since we are interested in the spurious noise
rather than the CSL signal. Solving the master equation for this initial state
with a qubit-BAW resonator coupling strength of $g/2\pi\approx
3\;\mathrm{MHz}$, we find that the qubit excited state population is only
suppressed to $P(\ket{e})\sim 10^{-4}$ due to dephasing noise. However if the
qubit is tuned to another mode on a timescale short compared to the period of
a Rabi oscillation and another swap operation is performed, the spurious
heating can be further suppressed. In other words, the qubit can be cooled by
swapping spurious heating due to quasiparticle poisoning and the Purcell
effect into the mechanical resonator. The saturation limit to how
significantly the qubit can be cooled using this protocol is the ratio of the
spurious excitation rate and the qubit-resonator coupling strength
$\frac{\Gamma_{\mathrm{QP}}+\Gamma_{P}}{g_{l,m}}$. The initial population of
the qubit is $\frac{\Gamma_{\mathrm{QP}}+\Gamma_{P}}{\gamma_{q}}$, therefore
our protocol is effective in cooling the qubit in the strong-coupling regime
$g_{l,m}>\gamma_{q}$. We predict this relation through a detailed balance
equation similar to Eq. (45) in which coupling to the mechanical resonator is
modelled as providing a decay channel for the qubit. After each swap operation
the qubit is tuned on-resonance with another mechanical mode and another swap
operation is performed in order to cool the qubit further, until the
saturation limit is obtained.

We now numerically confirm the predicted saturation limit of how much the
qubit can be cooled using our proposed procedure. If no heating additional to
the noise level corresponding to the collapse rate $\lambda_{\mathrm{min}}$ is
observed in the proposed experiment, then any value for the collapse rate of
$\lambda_{c}\geq\lambda_{\mathrm{min}}$ is strictly excluded. Therefore, we
will firstly solve the master equation for the initial state

$\begin{split}\hat{\rho}(0)&\approx
10^{-5}\ket{e}\ket{\mathrm{vac}}\bra{e}\bra{\mathrm{vac}}\\\
&+(1-10^{-5})\ket{g}\ket{\mathrm{vac}}\bra{g}\bra{\mathrm{vac}},\end{split}$
(48)

which assumes the dominant source of heating in the qubit is quasiparticle
poisoning. We will then compare to the solution of the master equation with an
initial state that assumes no quasiparticle or Purcell-induced heating
$\Gamma_{\mathrm{QP}}=\Gamma_{P}=0$, such that the only source of heating in
the coupled system is CSL-induced heating:

$\begin{split}\hat{\rho}(0)&\approx 1.2\times
10^{-6}\ket{g}\ket{l,m}\bra{g}\bra{l,m}\\\ &+(1-1.2\times
10^{-6})\ket{g}\ket{\mathrm{vac}}\bra{g}\bra{\mathrm{vac}}.\end{split}$ (49)

In this way, we compare the expected excited state population due to CSL alone
to the excited state population due to quasiparticle poisoning and the Purcell
effect. Assuming Bassi et. al’s lowest estimate for the collapse rate lower
bound of $\lambda_{c}=10^{-12}\;\mathrm{s}^{-1}$, the probability a higher-
order mechanical mode $l\sim 503$ will be in the $n=1$ Fock state due to
collapse-induced heating is $P=1\times 10^{-6}$. This probability scales
linearly with the collapse rate $\lambda_{c}$. Therefore, we identify the
minimum testable collapse rate as the collapse rate for which the excited
state probability due to CSL intersects the excited state probability due to
noise sources after a single swap operation. The quasiparticle-induced heating
is suppressed to its minimum possible value after a single qubit-mechanical
resonator swap operation ($\tau\approx 85\;\mathrm{ns}$). A measurement should
therefore be made $60-100\;\mathrm{ns}$ after switching on the coupling, and
the absence of any qubit excited state population above $P(\ket{e})\approx
6\times 10^{-6}$ would conclusively falsify CSL for collapse rates
$\lambda_{c}\geq 6\times 10^{-12}\;\mathrm{s}^{-1}$ at the expected
correlation length ($r_{c}\approx 100\;\mathrm{nm}$).

Figure 6:  Black lines: Qubit occupancy due to CSL as a function of time for
$\lambda_{c}=10^{-10}\;\mathrm{s}^{-1}$ and
$\lambda_{c}=10^{-12}\;\mathrm{s}^{-1}$, corresponding to full exclusion of
Adler and Bassi’s respective proposals. Dashed blue and light blue lines:
occupancy due to quasiparticle poisoning and the Purcell effect at the rates
of $(\Gamma_{\mathrm{QP}}+\Gamma_{\mathrm{P}})/2\pi\approx 30\;\mathrm{Hz}$
and $(\Gamma_{\mathrm{QP}}+\Gamma_{\mathrm{P}})/2\pi\approx 3\;\mathrm{Hz}$,
respectively.

In Fig 6 the dominant noise source is compared to the CSL heating
corresponding to a full exclusion of Adler and Bassi et al.’s lower bounds.
The qubit population due to quasi-particle poisoning (blue lines) and the
qubit population due to CSL heating (grey lines) are displayed in Fig 6, which
shows how these populations vary over a single swap operation. While the
heating rates that are in existing quasiparticle trap technology allow for a
full exclusion of Adler’s proposed lower bounds, quasiparticle heating reduced
by just one order of magnitude allows a full exclusion of Bassi’s proposal
after just one swap operation is complete ($gt/{2\pi}\approx 0.25$). The
crucial feature for suppressing quasiparticle is the swap operation, after a
single swap operation the quasiparticle heating in the qubit is minimised,
while the CSL heating in the qubit is maximised.

##  VI Measurement Induced heating

Non-QND effects arise due to higher order terms in $g/\Delta$ in the
dispersive approximation. Most significantly, the drive Hamiltonian acts as
creation and annihilation operators on eigenstates of the dispersive
Hamiltonian, while the eigenstates of the system in the lab frame are instead
eigenstates of the full interaction Hamiltonian
$H=\hbar\omega_{r}a^{\dagger}a-\hbar\omega_{q}\frac{\sigma_{z}}{2}+\hbar
g\left(a^{\dagger}\sigma_{-}+a\sigma_{+}\right)$. Mismatch between these
dressed eigenstates and the bare eigenstates of the dispersive Hamiltonian
causes the drive Hamiltonian ($H_{d}$) to induce qubit transitions in the
dressed basis:

$H_{d}=\varepsilon^{*}a+\varepsilon
a^{\dagger}\approx\left(\epsilon^{*}\bar{a}+\epsilon\bar{a}^{\dagger}\right)-\frac{g}{\Delta}\left(\epsilon^{*}\bar{\sigma}_{-}+\epsilon\bar{\sigma}_{+}\right),$
(50)

where $a$ and $\bar{a}$ are the bare and dressed creation operators
respectively. This mismatch between the bare and dressed creation and
annihilation operators in the readout Hamiltonian generates qubit excited
state population in the dressed basis given by
$P_{H}\approx\left|\frac{\epsilon g}{\Delta^{2}}\right|^{2}$ [17]. For a
measurement drive of $\epsilon/2\pi=25\;\mathrm{MHz}$ (allowing the readout
resonator to be populated with 2.5 photons within the first 10 $\mathrm{ns}$
being driven), a coupling strength and detuning of $g/2\pi=100\;\mathrm{MHz}$
and $\Delta/2\pi=2\;\mathrm{GHz}$ respectively we obtain $P_{H}\approx 4\times
10^{-7}.$

##  VII Thermal Noise

A large thermal noise background has precluded previous experiments from
resolving collapse-induced heating within the expected CSL parameter range
[18, 19]. This is because these experiments have operated with low frequencies
in the high temperature limit ($k_{B}T\gg\hbar\omega$), requiring them to
resolve CSL-induced heating above this thermal noise background. This
necessitates a precise knowledge of the thermal noise level, and can lead to
inconsistent interpretations [18]. Here this problem is solved by accessing
the regime where $k_{B}T\ll\hbar\omega$ [2]. Taking
$\omega=\omega_{q}=\omega_{r}$ and $k_{B}T\ll\hbar\omega$, the thermal heating
rates are
$\dot{n}_{\mathrm{th,r}}\approx\frac{\gamma_{r}}{2\pi}\exp\left(-\frac{\hbar\omega}{k_{B}T}\right)$
and
$\dot{n}_{\mathrm{th,q}}\approx\frac{\gamma_{q}}{2\pi}\exp\left(-\frac{\hbar\omega}{k_{B}T}\right)$
in the mechanical resonator and the qubit respectively. An upper bound on the
minimum testable collapse rate limited by thermal noise is therefore
$\lambda_{c,\mathrm{min}}=(\dot{n}_{\mathrm{th},r}+\dot{n}_{\mathrm{th},q})/D$.
Assuming the system is thermalised to the base-plate of a dilution fridge at
$T=10\;\mathrm{mK}$, while taking $\omega/2\pi=6.65\;\mathrm{Hz}$, we obtain a
minimum testable collapse rate of $\lambda_{c,\mathrm{min}}\approx
10^{-19}\;\mathrm{s}^{-1}$, many orders of magnitude below both Adler and
Bassi’s predicted lower bounds.

##  VIII Coloured CSL

One of the weaknesses of the white noise CSL model is that it is incompatible
with the collapse noise field having a physical origin [20]. If the collapse
noise field has a cosmological origin, it should have a non-white spectrum
with a high frequency cut-off [21]. Therefore, it is important to consider the
implications of a non-white spectrum for the bounds on the CSL model’s
phenomenological parameters $\lambda_{c}$ and $r_{c}$. In order to examine how
a non-white spectrum influences the decoherenece rate we examine the
formulation of CSL as a random potential [20]:

$\hat{V}_{\mathrm{cCsL}}(t)=-\frac{\hbar\sqrt{\lambda}}{m_{0}}\int\mathrm{d}\mathbf{z}\hat{M}(\mathbf{z})w(\mathbf{z},t).$
(51)

Here, $w(\mathbf{z},t)$ is a classical Gaussian noise characterised by

$\mathbb{E}[w(\mathbf{z},t_{1})]=0,\quad\mathbb{E}[w(\mathbf{z},t_{1})w(\mathbf{x},t_{2})]=\delta^{(3)}(\mathbf{z}-\mathbf{x})f(t_{1}-t_{2}),$
(52)

where $\mathbb{E}[\cdot]$ denotes the stochastic average. For the correlation
function $f(t)=\delta(t)$ the simple white-noise CSL model is obtained. In
order to model the effects of a coloured CSL noise spectrum, we consider an
exponentially decaying noise correlation function:

$f(t_{1}-t_{2})=\frac{\Omega_{c}}{2}e^{-\Omega_{c}|t_{1}-t_{2}|}.$ (53)

Here, $\Omega_{c}$ is the high frequency cut-off, and for
$\Omega_{c}\rightarrow\infty$ we obtain the delta function correlation
function corresponding to the simple white-noise CSL model. For an oscillator
of resonance frequency $\omega$, the coloured CSL noise spectrum modifies the
CSL-induced noise spectral density for the oscillator to be [20]

$S_{\mathrm{CSL}}=S_{\mathrm{c}}\cdot\frac{\Omega_{\mathrm{c}^{2}}}{\Omega_{\mathrm{c}}^{2}+\Omega^{2}},$
(54)

while in the simple white noise CSL model, the CSL-induced heating rate is
$\dot{n}_{c}=\lambda_{c}D$ (where $D$ is the CSL cross-section of the
oscillator). In the coloured CSL model the CSL-induced heating rate is

$\dot{n}_{c}=\lambda_{c}D\frac{\Omega_{\mathrm{c}^{2}}}{\Omega_{\mathrm{c}}^{2}+\Omega^{2}}.$
(55)

If the collapse noise field has a cosmological origin, it is expected to be at
$\sim 10^{10}-10^{11}\;\mathrm{Hz}$. In this case, we note that for
$\Omega_{\mathrm{c}}=10^{10}\;\mathrm{Hz}$ and
$\Omega_{\mathrm{c}}=10^{11}\;\mathrm{Hz}$ there is only a minor reduction
from the white noise CSL-induced heating rate of $30\%$ and $0.4\%$
respectively (assuming a mechanical mode frequency of
$\omega/2\pi=6.62\;\mathrm{GHz}$). Therefore, our proposed experiment has the
capability to test coloured CSL models.

##  IX Sample Experimental design

Here we propose a sample experimental design of our proposal similar to the
design used in Ref. [8]. The hybrid quantum system we are proposing to use
consists of a quartz phononic Fabry-Perot resonator that can support
longitudinally polarised acoustic modes (Fig. 7). We propose to perform phonon
counting using a superconducting qubit as a means of reading out phonons in a
bulk acoustic wave (BAW) resonator. The resonator is coupled to a frequency-
tunable Aluminium transmon (similar to that used in Ref. [6]) with a thin disk
of piezoelectric lead zirconate titanate (PZT). The dominant source of the
piezoelectricity is the PZT which generates stress $\sigma(\vec{x})$ due to
the transmon’s electric field $E(\vec{x})$ producing the phonon strain field
$s(\vec{x})$. The proposed flip-chip device consisting of a qubit on a non-
piezoelectric sapphire chip and the mechanical resonator on a separate quartz
chip is similar to the design used in Ref. [8]. Quartz is used for the chip
containing the mechanical resonator due to recent demonstrations of decay
rates as low as $\gamma_{r}/2\pi\approx 300\;\mathrm{Hz}$ in $\mathrm{GHz}$
frequency modes of curved surface quartz BAW resonators [7]. Sapphire is used
for the chip containing the qubit since it is a non-piezoelectric material,
preventing the generation of mechanical modes in the lower chip from producing
unwanted interference effects.

Figure 7:  Side view schematics of the proposed HBAR device. The proposed
flip-chip device is similar to the device used in Ref. [8]. The qubit is on a
non-piezoelectric sapphire chip while the acoustic resonator is on a separate
quartz chip. The phonon mode is confined to a cylindrical cavity of height $h$
with a curved surface in order to reduce the mechanical decay rate. The phonon
mode is coupled to a frequency-tunable transmon qubit with a thin disk of PZT
(green) of height $h_{\mathrm{PZT}}$. The electrode diameter $d_{e}$ is
smaller than the curved surface resonator diameter $d_{c}$ in order to produce
selective coupling to the fundamental transverse mode.

## References

##

  * Nimmrichter  _et al._ [2014] S. Nimmrichter, K. Hornberger, and K. Hammerer, Phys. Rev. Lett. 113, 020405 (2014).
  * Forstner  _et al._ [2020] S. Forstner, M. Zych, S. Basiri-Esfahani, K. E. Khosla, and W. P. Bowen, Optica 7, 1427 (2020).
  * Kharel  _et al._ [2018] P. Kharel, Y. Chu, M. Power, W. H. Renninger, R. J. Schoelkopf, and P. T. Rakich, APL photonics 3, 066101 (2018).
  * Adler [2007] S. L. Adler, Journal of Physics A: Mathematical and Theoretical 40, 2935 (2007).
  * Bassi  _et al._ [2010] A. Bassi, D.-A. Deckert, and L. Ferialdi, EPL (Europhysics Letters) 92, 50006 (2010).
  * Chu  _et al._ [2017] Y. Chu, P. Kharel, W. H. Renninger, L. D. Burkhart, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, Science 358, 199 (2017).
  * Renninger  _et al._ [2018] W. Renninger, P. Kharel, R. Behunin, and P. Rakich, Nature Physics 14, 601 (2018).
  * Chu  _et al._ [2018] Y. Chu, P. Kharel, T. Yoon, L. Frunzio, P. Rakich, and R. Schoelkopf, Nature 563, 666 (2018).
  * Fraga  _et al._ [2014] M. A. Fraga, M. A. Fraga, H. Furlan, H. Furlan, R. S. Pessoa, R. S. Pessoa, M. Massi, and M. Massi, Microsystem Technologies 20, 9 (2014).
  * Scala  _et al._ [2007] M. Scala, B. Militello, A. Messina, J. Piilo, and S. Maniscalco, Phys. Rev. A 75, 013811 (2007).
  * Boissonneault  _et al._ [2009] M. Boissonneault, J. Gambetta, and A. Blais, Phys. Rev. A 79 (2009).
  * Jeffrey  _et al._ [2014] E. Jeffrey, D. Sank, J. Y. Mutus, T. C. White, J. Kelly, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Megrant, P. J. J. O’Malley, C. Neill, P. Roushan, A. Vainsencher, J. Wenner, A. N. Cleland, and J. M. Martinis, Phys. Rev. Lett. 112, 190504 (2014).
  * Walter  _et al._ [2017] T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, and A. Wallraff, Phys. Rev. Applied 7 (2017).
  * Kulikov  _et al._ [2020] A. Kulikov, R. Navarathna, and A. Fedorov, Phys. Rev. Lett. 124, 240501 (2020).
  * Sete  _et al._ [2014] E. A. Sete, J. M. Gambetta, and A. N. Korotkov, Phys. Rev. B 89, 104516 (2014).
  * Sete  _et al._ [2015] E. A. Sete, J. M. Martinis, and A. N. Korotkov, Phys. Rev. A 92, 012325 (2015).
  * Khezri  _et al._ [2016] M. Khezri, E. Mlinar, J. Dressel, and A. N. Korotkov, Phys. Rev. A 94 (2016).
  * Vinante  _et al._ [2017] A. Vinante, R. Mezzena, P. Falferi, M. Carlesso, and A. Bassi, Phys. Rev. Lett. 119, 110401 (2017).
  * Vinante  _et al._ [2020] A. Vinante, M. Carlesso, A. Bassi, A. Chiasera, S. Varas, P. Falferi, B. Margesin, R. Mezzena, and H. Ulbricht, Phys. Rev. Lett. 125, 100404 (2020).
  * Carlesso  _et al._ [2018] M. Carlesso, L. Ferialdi, and A. Bassi, The European Physical Journal D 72, 1 (2018).
  * Adler and Bassi [2007] S. L. Adler and A. Bassi, Journal of Physics A: Mathematical and Theoretical 40, 15083 (2007).