Studying phonon coherence with a quantum sensor

Nanomechanical oscillators hold the potential to serve as long-lived memories for computation Hann et al. (2019); Chamberland et al. (2022), transducers for communication Safavi-Naeini and Painter (2011); Vainsencher et al. (2016); Jiang et al. (2020), and high-precision quantum sensors Mason et al. (2019). Their ability to interact with superconducting qubits through the piezoelectric effect has allowed mechanical systems to be brought into the quantum regime O’Connell et al. (2010); Chu et al. (2017, 2018); Satzinger et al. (2018); Arrangoiz-Arriola et al. (2019); Sletten et al. (2019); Bienfait et al. (2019, 2020); Chu and Gröblacher (2020); Wollack et al. (2022). In order to fully realize the potential of this hybrid platform, it is crucial to understand decoherence processes affecting mechanical oscillators in the quantum limit. For mechanical states with large numbers of phonons ($\bar{n}\gg 1$), the established work has followed semiclassical spectroscopic methods to detect time-averaged energy loss Manenti et al. (2016); MacCabe et al. (2020); Wollack et al. (2021) and frequency noise Gao et al. (2008a); Tubsrinuan et al. (2022) in the resonator. More recent studies of mechanical devices in the few-phonon regime have used a superconducting qubit to perform single-phonon characterization Chu et al. (2017); Satzinger et al. (2018); Wollack et al. (2022); von Lupke et al. (2022). In this manuscript, we use a Ramsey protocol with a superconducting qubit that allows us to perform time- and phonon number-resolved measurements on a mechanical resonator Lachance-Quirion et al. (2020); von Lupke et al. (2022). This technique, which could easily be extended to other bosonic systems, affords us a highly granular picture of how larger mechanical states of an oscillator evolve in time, shedding new light on the quantum nonlinear dissipation processes.

The mechanical dissipation dynamics we observe are consistent with emission into a bath of rapidly dephasing two-level system (TLS) defects, similar to those which exist in amorphous glasses and imperfect crystalline materials Phillips (1987). TLS can couple to electromagnetic and elastic fields, creating a dissipation channel for nearby modes Martinis et al. (2005); Muller et al. (2019); McRae et al. (2020). The dynamics of TLS-induced decay and decoherence at microwave frequencies have been the subject of extensive study in both electromagnetic Gao et al. (2008b, a); Quintana et al. (2014); Klimov et al. (2018); Lisenfeld et al. (2019) and nanomechanical systems Hoehne et al. (2010); Suh et al. (2013); Manenti et al. (2016); MacCabe et al. (2020); Wollack et al. (2021); Andersson et al. (2021). Their contribution to microwave loss has been shown to depend on power and temperature, allowing the TLS-induced loss tangent to be extracted from measurements of the resonator’s frequency and quality factor Gao et al. (2008b); Wollack et al. (2021). In particular, TLS contribution to dielectric loss is known to be saturable – it is suppressed when the intra-cavity energy exceeds a critical threshold. In this regime, the TLS spontaneous excitation and emission rates equilibrate as the TLS reach a steady state in which they no longer absorb energy from the resonator mode. Prior studies of TLS saturation have involved monitoring the scattered response of a resonator for a range of drive powers Quintana et al. (2014); Gao et al. (2008b); Wollack et al. (2021). In this work, the qubit nonlinearity allows for quantum non-demolition measurement of the mechanical state via established circuit QED techniques. We leverage this approach to observe signatures of TLS-induced dissipation in the time domain, which we reproduce using a simple numerical model.

The mechanical oscillators under study are one-dimensional phononic crystal cavities made of thin-film lithium niobate (LN). While the device contains two resonators, these experiments focus on only one cavity, which supports the higher frequency mode (Fig 1a). The cavity is formed from periodic structures, acting as acoustic mirrors, which support a defect site from both ends Arrangoiz-Arriola and Safavi-Naeini (2016). The periodicity of the mirrors gives rise to a full phononic bandgap which protects the localized mode from phonon radiation channels (supplementary Fig. S1 and Section I). This design approach, and the large coupling strengths arising from the strong piezoelectric effect of LN, make these devices a compelling platform for studies in quantum acoustics, because it robustly removes the effects of clamping and other linear scattering losses. Furthermore, the small mode volume of these cavities make them an ideal candidate for studies of TLS as mechanical loss channels, since it yields stronger coupling to fewer TLS.

We assemble our device in a hybrid flip-chip architecture. The mechanical resonators are fabricated on the top chip, and the bottom chip includes a superconducting qubit, control lines, and a coplanar waveguide readout resonator (Fig. 1b). The qubit is a frequency-tunable planar transmon, with an on-chip flux line providing capability for both static and pulsed frequency control (Fig 1b, inset). The qubit fabrication is described in prior work Wollack et al. (2022); Arrangoiz-Arriola et al. (2019) and follows methods developed in Ref. Kelly (2015). As illustrated in Fig. 1c, the mechanics fabrication process begins with thin-film LN on a silicon substrate. We thermally anneal the sample for 8 hours at 500 C before thinning the LN by argon ion milling to a target final thickness of 250 nm. The nanomechanical structures are then patterned by electron-beam lithography and argon ion milling using a hydrogen silesquioxane (HSQ) mask. We remove redeposited material from the LN sidewalls in a heated bath of dilute hydrofluoric acid, followed by a piranha etch. Next, we define metallic layers with a combination of electron-beam and photolithography, including bandages to create the relevant galvanic connections. We perform an oxygen plasma descum before each deposition to minimize the presence of polymer residues at metallic interfaces. Finally, we undercut the nanomechanical structures from the substrate with a masked xenon difluoride (XeF₂) dry etch. This concludes the mechanical device fabrication, before the LN chip is aligned and bonded to an accompanying qubit chip (Fig. 1d) in a submicron die bonder Wollack et al. (2022); Satzinger et al. (2019). The qubit is capacitively coupled to the mechanical mode across a small vacuum gap, which is defined by the flip-chip separation distance (Fig. 1e).

The completed flip-chip package is cooled to a temperature of 10 mK in a dilution refrigerator. As we bias the qubit away from its maximum frequency $\omega_{\text{ge}}^{\text{max}}/2\pi$ = 2.443 GHz, we measure the qubit $T_{1}$ over a wide tuning range to find a mean value $T_{1}=4.9{\,\pm\,}2.3\,\mu\text{s}$, with $T_{2}=1.4\,\mu\text{s}$ at the flux sweet spot. When the qubit is rapidly tuned into resonance with the mechanics at $\omega_{\text{m}}/2\pi=2.339$ GHz, the resulting Rabi oscillations allow us to extract the coupling rate $g/2\pi=10.5\,\pm\,0.1$ MHz. Further details on basic device characterization can be found in a prior manuscript Wollack et al. (2022).

In the following experiments, we statically bias the qubit near the mechanical mode and use it as a probe to study dissipation and dephasing of mechanical coherent states. We extract information about the mechanical state through the qubit spectrum by means of their dispersive coupling. The qubit and mechanics are coupled through the piezoelectric effect to give an interaction Hamiltonian $\hat{H}_{\text{int}}=g(\hat{b}+\hat{b}^{\dagger})\hat{\sigma}_{x}$, with $\hat{b}$ corresponding to the mechanics and Pauli operators $\hat{\sigma}$ corresponding to the qubit. When the two modes are far detuned, the effective Hamiltonian for the system is given by Koch et al. (2007):

$\displaystyle\hat{H}_{\text{eff}}$

$\displaystyle=$

$\displaystyle\omega_{\text{m}}\hat{b}^{\dagger}\hat{b}+\frac{1}{2}\bigg{(}\omega_{\text{ge}}+2\chi\hat{b}^{\dagger}\hat{b}\bigg{)}\hat{\sigma}_{z}\,.$

In this limit, the qubit spectrum acquires a frequency shift $2\chi$ for each populated energy level in the mechanical oscillator. The magnitude of this shift depends on the qubit anharmonicity $\alpha_{\text{q}}$, mechanical coupling rate $g$, and detuning $\Delta=\omega_{\text{ge}}-\omega_{\text{m}}$ as:

$$\chi=-\frac{g^{2}}{\Delta}\frac{\alpha_{\text{q}}}{\Delta-\alpha_{\text{q}}}\,.$$

Since this value varies with the qubit frequency, it is important to choose an operating point which yields a large $\chi$ while keeping the system in the strong dispersive limit. In this regime, $\chi$ exceeds all decoherence rates of the system, but the detuning is sufficiently large to prevent significant linear coupling between the two modes. For our experiments, we statically bias the qubit to $\omega_{\text{ge}}/2\pi=2.262$ GHz, which corresponds to $|\Delta|\approx 8g$ and a measured dispersive shift $2\chi/2\pi=-1.48\pm 0.05$ MHz.

When the mechanical mode is populated, we can perform a Ramsey measurement on the qubit to resolve the array of dispersive shifts and thus obtain the full phonon number distribution of the mechanical state. The Ramsey signal $S(t)$ takes the form of a sum of oscillating terms whose frequencies differ by integer multiples of $2\chi$, indexed by the Fock number $n$, with an exponentially decaying envelope:

$\displaystyle S(t)=\sum_{n=0}^{n_{\text{max}}}A_{n}e^{-\kappa t}\cos\bigg{(}(\omega_{0}+2\chi n)t+\varphi_{n}\bigg{)}\,.$

(1)

We fit the Ramsey signal to Eq. 1 with fit parameters $A_{n}$, $\kappa$, and $\chi$. We extract the phonon number distribution from the normalized amplitudes, $P(n)=A_{n}/\sum_{n}A_{n}$. In contrast to previous work Wollack et al. (2022) where the exponential decay rate in Eq. 1 was taken to depend linearly on $n$, we choose here to model decay as independent of the phonon occupation number (supplementary Section II).

We use the Ramsey protocol to study decay of mechanical coherent states. We first apply a displacement pulse $\hat{D}_{\alpha}$ at the mechanical frequency to the qubit’s XY line (Fig. 2a), which drives the mechanics into a coherent state. After a variable delay $\tau$, we perform a Ramsey measurement on the qubit to extract the mechanical $P(n)$, as described above (Fig. 2b,c). From this, we calculate the average number of phonons, $\bar{n}=\sum_{n=0}n\cdot P(n)$. Repeating this sequence for a range of $\tau$ values constitutes a phonon-number resolved ringdown measurement of a mechanical coherent state (Fig. 2d). The measurement provides a high degree of resolution, revealing not only the decaying mean phonon number in the resonator, but also how the distribution evolves in time (Fig. 2e and Fig. S2).

We observe that the mechanical dissipation follows a double-exponential trajectory, with an initial fast decay followed by significantly slower relaxation. Similar multi-exponential behavior appears in single-phonon measurements in a prior study of the same device Wollack et al. (2022). By varying the duration and amplitude of the initial displacement pulse, we study how the energy decay and saturation dynamics depend on the size of the initial mechanical state (Fig. 2d). We fit each decay curve to a double-exponential form:

$\displaystyle\bar{n}(\tau)=a_{1}e^{-\kappa_{1}\tau}+a_{2}e^{-\kappa_{2}\tau}.$

(2)

We examine the behavior of the fast decay $\kappa_{1}/2\pi\simeq 30-70$ kHz, slow decay $\kappa_{2}/2\pi\simeq 1-3$ kHz, and their corresponding weights $a_{1}$, $a_{2}$ as we vary the initial mean phonon number $\bar{n}_{0}$.

Nonexponential decay has been observed in superconducting qubits due to fluctuations in quasiparticle population Gustavsson et al. (2016) as well as in graphene mechanical resonators operated in a strongly nonlinear regime Singh et al. (2016). Neither mechansism is plausible in our experiments. We argue that the multi-exponential energy decay in our system is caused by resonant interaction of the mechanics with a handful of rapidly dephasing TLS. The initial fast decay occurs as the mechanical mode emits energy into the TLS. We note that we do not observe any coherent oscillations between the two coupled systems. As the TLS become populated, they lose their ability to absorb more phonons, thus becoming saturated. Consequently, the rapid decay ends and slower decay emerges, caused by either energy relaxation of the TLS or other dissipation processes affecting the resonator. We perform numerical modeling to reproduce the salient features of this process, described later in this paper.

A modification of our measurement protocol allows us to extract the mechanical dephasing time $T_{\text{2m}}$ using coherent states of motion with varying amplitude. As shown in Fig. 3a, we now displace the mechanical state twice: once to initialize the measurement with a coherent state, and again after a variable delay $\tau$. The two displacement pulses have identical amplitude and duration, but the second pulse has a programmed phase $\phi$ that is swept from 0 to $2\pi$ for each value of $\tau$. The protocol concludes with a Ramsey measurement to characterize the final mechanical state.

We illustrate the principle behind this displacement interferometry protocol in Fig. 3b. The first displacement initializes a coherent state (Fig. 3b,i). After the second displacement, the resulting mechanical states lie on a circle in the complex plane. In the absence of dephasing ($\tau\ll T_{\text{2m}}$), maximal displacement $2\alpha$ occurs at $\phi=0$, and the final state is returned to the origin when the pulses are out of phase by $\phi=\pi$ (Fig. 3b,ii). In the presence of dephasing ($\tau\gtrsim T_{\text{2m}}$), this resulting path is shifted in the complex plane and no longer intersects the origin (Fig. 3b,iii). In these measurements, the final state size $\bar{n}$ depends sinusoidally on $\phi$, and the effect of dephasing can be seen in the reduced oscillation amplitude for larger delays (Fig. 3c). These amplitudes decay exponentially in time at a rate $\gamma_{\text{2m}}=1/T_{\text{2m}}$ (supplementary Section III). We fit the interference data to the following form to extract the amplitude $C$ and offset $\bar{n}_{\text{off}}$:

$\displaystyle\bar{n}(\phi)=C\hskip 1.42271pt\text{cos}(\phi+\phi_{0})+\bar{n}_{\text{off}}.$

(3)

The phase offset $\phi_{0}$ is dependent on the delay $\tau$, and could arise due to frequency uncertainty.

A few representative interferometry traces are shown in Fig. 3d for an initial state $\bar{n}_{0}=2.29$. Each interferogram is plotted relative to its offset $\bar{n}_{\text{off}}$ for visual clarity. We note that the offsets also decay in time at a rate $\gamma_{\text{1m}}=1/T_{\text{1m}}$, but the timescale of the interferometry data is insufficient to accurately determine this value. In Fig. 3e we plot the extracted offsets $\bar{n}_{\text{off}}$ (inset), amplitudes $C$ and exponential fit for this data set, from which we extract $T_{\text{2m}}=2.2\,\pm\,0.2\,\mu$s.

The effect of dephasing is also visible in the phonon number distribution of the mechanical oscillator. In principle, a coherent state undergoing amplitude decay should remain a coherent state and therefore maintain a Poisson distribution. However, in the presence of dephasing, the state evolves into something other than a coherent state whose $P(n)$ may deviate significantly from the Poissonian form. This effect is evident in Fig. 3f, where we examine $P(n)$ extracted from data for two representative states: one with $\tau=0$ (top) and another with $\tau=5.0\,\mu$s (bottom). Each distribution is compared to a coherent state $|\beta\rangle$ whose $P(n)$ is given by a Poisson distribution with the same average phonon number as the data, $|\beta|^{2}=\bar{n}$. The top panel shows good agreement between the data (red bars) and the coherent state distribution (shaded blue), while the bottom panel shows a pronounced divergence between the two. We observe a similar effect in a simulation of the dephasing process (Fig. 4f).

We repeat the dephasing measurement for three distinct initial states and find, perhaps surprisingly, that the dephasing rate is reduced for larger phonon states. In a naive model, we expect the mechanical frequency jitter to be either primarily due to fluctuating off-resonant TLS, which would not be saturated by weakly driving the mechanics, or from excitation and relaxation of resonant TLS, which occur even when they are saturated. Saturation implies that the rates of phonon emission and absorption into the TLS bath are roughly equal, cancelling their effect on energy decay and leading to longer $T_{\text{1m}}$ as we observe. However, dephasing should still occur if emission and absorption of phonons are occuring incoherently. An increase in $T_{\text{2m}}$ may therefore point to coherence of the TLS bath with the driving field.

To better understand these decay and decoherence signatures, we develop and study numerically Johansson et al. (2011) a highly simplified model of a mechanical mode interacting with a small number of TLS Remus et al. (2009). Our model includes a collection of $N$ identical two-level oscillators, with annihilation operators $\hat{a}_{k}$, coupled at a rate $g_{\text{TLS}}$ to a harmonic oscillator (Fig. 4a). The TLS have intrinsic decay and dephasing rates $\gamma_{1}$ and $\gamma_{2}$, respectively, and are detuned from the mechanical frequency by $\Delta_{\text{TLS}}$. We assume there is no TLS-TLS coupling. In the frame of the TLS, the system Hamiltonian is given by:

$\displaystyle\hat{H}=\Delta_{\text{TLS}}\hat{b}^{\dagger}\hat{b}+\sum_{k=1}^{N}g_{\text{TLS}}(\hat{a}_{k}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a}_{k}).$

(4)

We initialize the mechanics and TLS in thermal states, each populated to a level $n_{\text{th}}=0.05$. Next, we prepare the mechanical mode in a coherent state with an instantaneous displacement $\hat{D}_{\alpha}$ with variable amplitude. We then allow the system to freely evolve under the action of the Hamiltonian in Eq. 4 for a total duration of 50 $\mu$s, with collapse operators acting on the TLS. The operators $\sqrt{\gamma_{1}(n_{\text{th}}+1)}\hat{a}_{k}$ and $\sqrt{\gamma_{1}n_{\text{th}}}\hat{a}_{k}^{\dagger}$ account for spontaneous emission and absorption in the TLS, while $\sqrt{\gamma_{2}/2}\hat{a}_{k}^{\dagger}\hat{a}_{k}$ accounts for their dephasing. Dissipation and decoherence of the phonon mode are induced by interactions with the TLS; we do not include collapse operators for the mechanical mode itself.

The numerical solution to this Lindblad master equation returns the density matrix of the total system for every point in time. We examine the evolution of the mechanical mode’s $P(n)$ and state size $\bar{n}$. Similar to experiment, the model shows approximately double-exponential energy decay in Fig. 4b, with weights and rates that depend on the initial state size. We fit each simulated decay profile to Eq. 2; to constrain the fit, we fix $\kappa_{1}$ to the mean of the experimentally observed values while allowing $a_{1}$, $a_{2}$ and $\kappa_{2}$ to vary.

Although the measurements are insufficient to precisely determine the system parameters $N$, $g_{\text{TLS}}$, $\Delta_{\text{TLS}}$, $\gamma_{1}$ and $\gamma_{2}$, they provide guidance in our choice of values, and the resulting model helps to qualitatively understand the experimental data. It is possible to estimate the microscopic Hamiltonian parameters of the TLS from properties of the host material and the mechanical eigenmode (supplementary Section IV and Fig. S3). From the expected distribution of these parameters, the model produces reasonable agreement with experiments for a range of values. In Fig. 4b,c we plot the results for $N=5$, $g_{\text{TLS}}/2\pi=33$ kHz and $\Delta_{\text{TLS}}/2\pi=100$ kHz. The lack of evidence of coherent mechanics-TLS energy exchange suggests that the TLS are rapidly dephasing, $\gamma_{2}/g_{\text{TLS}}\simeq 20-30$. For this simulation, we use $\gamma_{2}/2\pi=660$ kHz and $\gamma_{1}/2\pi=4.0$ kHz. In Fig. 4b we plot the simulated energy decay trajectories (black) and corresponding fits (red) for a few initial states, and Fig. 4c shows good agreement between the fit metrics extracted from these simulations (lines) and the experimental data (points). See Fig. S3 for the corresponding TLS behavior.

To model the dephasing experiment, we apply a second displacement operation with phase $\phi\in[0,2\pi]$ at a range of delays in the time-evolved mechanical state trajectory. Following the same procedure as before, we extract the sinusoidal amplitudes from the phonon interference fringes (Fig. 4d) and determine the TLS-induced $\gamma_{\text{2m}}$ from the decaying amplitudes (Fig. 4d, inset). The results shown in Fig. 4d-f are computed with a stronger coupling rate ($g_{\text{TLS}}/2\pi=0.33$ MHz, $N=5$, $\gamma_{1}/2\pi=4.0$ kHz, $\gamma_{2}=20g_{\text{TLS}}$ and $\Delta_{\text{TLS}}=3g_{\text{TLS}}$). In Fig. 4e, we compare the extracted $\gamma_{\text{2m}}$ for these parameters to the experimental data, as well as to the weakly-coupled simulation results. For the weak-coupling parameters which match the energy decay experiments well, we find that the extracted $\gamma_{\text{2m}}$ (blue line) is significantly lower than the experimentally observed values (points). With the stronger coupling rate, $g_{\text{TLS}}/2\pi=0.33$ MHz, the extracted $\gamma_{\text{2m}}$ values (green line) fall within the range of our experimental observations. We note that for both values of $g_{\text{TLS}}$, the model predicts a reduced dephasing rate for larger phonon states. Finally, in Fig. 4f we show the simulated $P(n)$ for a short delay (top) and longer delay (bottom), where the dephasing effects can be seen as a deviation from the coherent state distribution.

The discrepancy in optimal simulation parameters suggests that our model provides an incomplete picture of the decoherence processes. For instance, it is possible that the mechanical mode suffers from additional dephasing channels, such as thermal excitation of low frequency TLS. It is also possible that more complicated bath dynamics, such as a distribution of non-identical defects with TLS-TLS interactions, are necessary to more accurately describe our data.

In conclusion, we have applied a dispersive phonon number measurement to study coherent states in a phononic crystal resonator. We have examined how the energy decay and dephasing of these states depend on the initial phonon state size, and have reproduced some of these signatures using a simple model incorporating an ensemble of saturable TLS. Future studies would benefit from a more complex computational model, as well as additional measurements to probe the saturation dynamics, such as spectral hole burning Andersson et al. (2021). Our results have direct relevance for bosonic error-correction schemes involving coherent states of phonons Chamberland et al. (2022). This measurement technique builds a foundation for further explorations of fundamental physics in quantum acoustic systems, including observing quantum jumps in mechanical states Sun et al. (2014).

The authors would like to thank M. I. Dykman and M. L. Roukes for useful discussions. The authors also thank R. G. Gruenke, K. K. S. Multani, N. R. Lee, T. Makihara, and O. Hitchcock. We acknowledge the support of the David and Lucille Packard Fellowship. This work was funded by the U.S. government through the Office of Naval Research (ONR) under grant No. N00014-20-1-2422, the U.S. Department of Energy through Grant No. DE-SC0019174, and the National Science Foundation CAREER award No. ECCS-1941826. A.Y.C. was supported by the Army Research Office through the Quantum Computing Graduate Research Fellowship as well as the Stanford Graduate Fellowship. E.A.W. was supported by the Department of Defense through the National Defense & Engineering Graduate Fellowship. Device fabrication was performed at the Stanford Nano Shared Facilities (SNSF), supported by the National Science Foundation under grant No. ECCS-1542152, and the Stanford Nanofabrication Facility (SNF). The authors wish to thank NTT Research for their financial and technical support.

A.Y.C. and E.A.W. designed and fabricated the device and performed the experiments. A.Y.C. analyzed the data with E.A.W. assisting. A.Y.C. performed the simulations and wrote the manuscript with A.H.S-N. assisting. A.H.S.-N. supervised all efforts.

E. Alex Wollack is currently a research scientist at Amazon, and A. H. Safavi-Naeini is an Amazon Scholar. The other authors declare no competing financial interests. Correspondence and requests for materials should be addressed to A. H. Safavi-Naeini ([email protected]).

The datasets generated and analysed for the current study are available from the corresponding author on reasonable request.

To optimize the mechanical resonator design, we perform finite-element simulations using COMSOL Multiphysics, similar to those detailed in Ref. Arrangoiz-Arriola et al. (2019). Our model includes a small sidewall angle $\theta_{\text{sw}}=10^{\circ}$ and corners rounded with radius $r=80$ nm to account for expected fabrication imperfections. We use $X$-cut, MgO-doped lithium niobate, with the crystal extraordinary axis oriented perpendicular to the cavity propagation axis. The frequency placement and bandwidth of the phononic bandgap depend on the geometry of the phononic crystal mirror cells. Our choice of design yields a simulated bandgap extending from approximately 1.9 to 2.5 GHz (Fig. S1a).

We also perform simulations of the mechanical cavity, which includes the defect site suspended by 4 mirror cells on each side. This is chosen for computational simplicity; the experimental device includes 8.5 mirror cells on each side. The resonant frequency of the localized mode is determined by the defect width, $L_{x}$. Target dimensions for both the defect and mirror cells are reported in Table S1. From the full structure simulations, we can extract the electromechanical admittance $Y(\omega)$ by simultaneously solving the electrostatic and elastic constitutive relations, coupled by the piezoelectric effect Arrangoiz-Arriola and Safavi-Naeini (2016). By applying an oscillating voltage boundary condition at the electrode surface and computing the induced current, we obtain the admittance as the ratio of the two. The simulated results are shown in Fig. S1d with a fit to an equivalent circuit model. From this fit, we extract the Butterworth-van Dyke circuit parameters $C_{0}=213.5$ aF, $C_{\text{m}}=51.4$ aF, and $L_{\text{m}}=90.9$ $\mu$H, corresponding to the equivalent circuit shown in the inset.

Though a decay rate proportional to phonon number is theoretically justified for linear decay of a bosonic mode, it is less suitable for our case, as coupling to TLS leads to more complex behavior. We attempted several approaches, including the original form from Ref. Wollack et al. (2022) as well as having independent $\kappa_{n}$ for each Fock state. We found that for the types of phonon states being analyzed in this study, which have larger average phonon numbers than our previous work, a constant decay of $\kappa$ assumed for all Fock states yields the best fit quality and stability. We roughly justify this by noting that since TLS are saturated by driving, larger Fock states will tend to have a relatively smaller decay than would be expected from a nonsaturable linear dissipation channel, and so a constant decay across all states is a useful heuristic for fitting the Ramsey signal. Importantly, we only assume this state-independent decay to fit the 700 nanosecond duration of the Ramsey signal, $t$. The time-dependent phonon occupations are extracted from data over the longer period of measurement, $\tau$.

A harmonic oscillator with resonant frequency $\omega_{0}$ and annihilation operator $\hat{b}$ can undergo frequency fluctuations due to a stochastic process, $\eta(t)$, leading to a modulated resonant frequency $\omega(t)=\omega_{0}+\eta(t)$. In the frame rotating at $\omega_{0}$, the quantum master equation which describes this process is given by:

$\displaystyle\frac{d\rho}{dt}=\mathcal{D}\big{[}\frac{1}{\sqrt{T_{\text{2m}}}}\hat{b}^{\dagger}\hat{b}\big{]}\rho$

(S1)

where $\mathcal{D}$ represents the Lindblad dissipator. Classically, these fluctuations of the harmonic oscillator’s frequency lead to an equation of motion for a coherent state:

$\displaystyle\frac{d\alpha}{dt}=-i\eta(t)\alpha$

(S2)

which, when integrated, leads to a coherent state amplitude $\alpha(t)=\alpha(0)e^{i\theta(t)}$ where $\theta(t)=\int_{0}^{t}\eta(t)dt$. This is reflected in the density matrix as

$\displaystyle\rho(0)=|\alpha\rangle\langle\alpha|\rightarrow\rho(t)=\sum_{\theta}p(\theta)|\alpha e^{i\theta}\rangle\langle\alpha e^{i\theta}|.$

(S3)

The second displacement operation in the protocol of Fig. 3, with amplitude and phase $\alpha e^{i\phi}$, leads to:

$\displaystyle D\rho(t)D^{\dagger}=\sum_{\theta}p(\theta)|\alpha e^{i\theta}+\alpha e^{i\phi}\rangle\langle\alpha e^{i\theta}+\alpha e^{i\phi}|.$

(S4)

The resulting phonon number in the oscillator is then given by

	$\displaystyle\langle n\rangle$	$\displaystyle=\sum_{\theta}p(\theta)|\alpha(e^{i\theta}+e^{i\phi})|^{2}$
		$\displaystyle=|\alpha|^{2}\sum_{\theta}p(\theta)(2+e^{i\phi-i\theta}+e^{i\theta-i\phi})$
		$\displaystyle=|\alpha|^{2}\big{(}2+e^{i\phi}\langle e^{-i\theta}\rangle+e^{-i\phi}\langle e^{i\theta}\rangle\big{)}$		(S5)

The expectation values can be expanded as:

$\displaystyle\langle e^{\pm i\theta}\rangle=\langle 1\pm i\theta-\frac{\theta^{2}}{2}+...\rangle$

(S6)

For a stochastic process $\eta(t)$, the expected value for terms linear in $\theta$ are zero: $\langle\theta\rangle=\langle-\theta\rangle=0$. For short times $\Delta t$, the quadratic terms are given by the correlation function $S_{\theta\theta}$:

$\displaystyle\langle\theta^{2}\rangle=S_{\theta\theta}(\omega=0)\Delta t=\frac{2\Delta t}{T_{\text{2m}}}$

(S7)

More generally, $\langle\theta^{n}\rangle=0$ for odd $n$, and $\langle\theta^{2n}\rangle\propto(\Delta t)^{n}$. This allows us to write Eq. S6 in terms of the time increment $\Delta t$ and decoherence time $T_{\text{2m}}$:

$\displaystyle\langle e^{\pm i\theta}\rangle=1-\frac{\Delta t}{T_{\text{2m}}}+\mathcal{O}(\Delta t)^{2}=e^{-t/T_{\text{2m}}}$

(S8)

Making this substitution in Eq. S5 brings us to the expected phonon number relation:

$\displaystyle\langle n\rangle=2|\alpha|^{2}\big{(}1+e^{-t/T_{\text{2m}}}\text{cos}\hskip 1.42262pt\phi)$

(S9)

The interaction Hamiltonian describing a single TLS coupled to a strain field is given in Refs. Gao (2008); Behunin et al. (2016) as:

$\displaystyle\hat{H}_{\text{int}}=\bigg{(}\frac{\Delta_{\text{as}}}{\varepsilon}\hat{\sigma}_{z}+\frac{\Delta_{0}}{\varepsilon}\hat{\sigma}_{x}\bigg{)}\gamma\cdot\xi$

(S10)

Here, $\hat{\sigma}_{x}$ and $\hat{\sigma}_{z}$ are Pauli operators for the TLS, $\gamma$ is the elastic dipole moment, and $\xi$ is the strain field. $\Delta_{0}$ and $\Delta_{\text{as}}$ represent, respectively, the tunneling energy and asymmetry energy from the bare TLS Hamiltonian, and $\varepsilon=\sqrt{\Delta_{0}^{2}+\Delta_{\text{as}}^{2}}$ gives the associated eigenenergies, $\pm\varepsilon/2$. We can use this relation to estimate the expected coupling rate between the TLS and the strain field, $g_{\text{TLS}}=\gamma\cdot\xi/\hbar$. The elastic dipole moment can be extracted from material constants and the measured loss tangent for our material platform:Behunin et al. (2016); Muller et al. (2019)

$\displaystyle\gamma=\sqrt{\delta^{0}_{\text{TLS}}\frac{\rho v^{2}}{\pi P_{0}}}.$

(S11)

Here $v$ is the speed of sound in the host material, $\rho$ the mass density, and $P$ the spectral and spatial density of TLS states. For simplicity, we assume the longitudinal and transverse components of $\gamma$ to be equivalent, and we use the RMS average of all strain components $\xi_{ij}$ in our calculation.

From a finite-element simulation of the mechanical eigenmode, it is useful to calculate the effective mass

$\displaystyle m_{\text{eff}}\equiv\frac{\int dV\rho|u(\vec{r})|^{2}}{\text{max}|u(\vec{r})|^{2}}$

(S12)

from which we can extract the displacement amplitude $x_{\text{zpf}}=\sqrt{\hbar/2m_{\text{eff}}\omega_{\text{m}}}$, which can be understood as the maximum zero-point fluctuation displacement. For the mode in question, we find $m_{\text{eff}}=440$ fg and $x_{\text{zpf}}=2.9$ fm. We can also calculate the average zero-point strain amplitude over the volume of the resonator

$\displaystyle\bar{\xi}_{\text{zpf}}=\frac{x_{\text{zpf}}}{\text{max}|u(\vec{r})|}\bar{\xi}$

(S13)

where $\bar{\xi}$ represents the volume-averaged RMS value of all strain components,

$\displaystyle\bar{\xi}=\sqrt{\frac{1}{6V}\int dV(\xi_{xx}^{2}+\xi_{xy}^{2}+\xi_{xz}^{2}+\xi_{yy}^{2}+\xi_{yz}^{2}+\xi_{zz}^{2})}.$

(S14)

For the phononic crystal resonator mode, we find $\bar{\xi}_{\text{zpf}}=1.6\times 10^{-9}$.

The TLS-induced loss tangent at zero temperature $\delta^{0}_{\text{TLS}}$ can be determined experimentally. Previous measurements of the temperature-dependent frequency shift of resonators with a similar geometry in this material platform provide an estimated value of the $F\delta^{0}_{\text{TLS}}$ product Wollack et al. (2021). For our system, where we approximate that the TLS are distributed through the entire volume of the cavity, the filling factor is simply $F=1$.

For the TLS density of states $P_{0}$, we assume a value of order $10^{44}-10^{46}$ 1/J$\cdot$m³ based on literature values extracted for silica whose measured loss tangents are comparable to ours Behunin et al. (2016); Phillips (1987). We calculate the expected distributions for $\gamma$ and $g_{\text{TLS}}$ using Eq. S11 with randomly sampled $P_{0}$ and $\delta^{0}_{\text{TLS}}$. For each iteration, we randomly sample $P_{0}=10^{\lambda_{1}}$ with $\lambda_{1}$ uniformly distributed between $[44,46]$ and $\delta^{0}_{\text{TLS}}=10^{\lambda_{2}}$ with $\lambda_{2}$ uniformly distributed between $[-4.5,-4]$. We calculate these values using 10,000 iterations to compute the expected distribution for elastic dipole $\gamma$ and coupling rate $g_{\text{TLS}}$, shown in Fig. S3d and S3e.

We can also use $P_{0}$ to estimate $N$, the number of TLS interacting with the mechanical mode. The relevant number to include in our numerical model is the number of TLS within the spatial volume of the resonator $V$ and within a pertinent bandwidth $\delta\omega$ of the mechanical frequency. This bandwidth can be evaluated as $\delta\omega=\text{max}(g_{\text{TLS}},\gamma_{2})$ where $\gamma_{2}$ is the decoherence rate of the TLS. We must also consider the temperature-dependent distribution of the TLS, $P\neq P_{0}$. For the operating temperature of our experiment ($T=10$ mK), this is found to be $P\simeq 10P_{0}$ Duquesne and Bellessa (1985). This allows us to extract

$\displaystyle N=10P_{0}\times V\times\hbar\delta\omega$

(S15)

for a known $P_{0}$ and $\gamma_{2}$. The calculations in existing TLS literature predict a much smaller $\gamma_{2}$ arising from TLS-TLS interactions Gao (2008). Our simulation findings are not consistent with this prediction; the model reproduces our experimentally observed effects only in the limit of much faster $\gamma_{2}$. For these values $\gamma_{2}/2\pi\gtrsim 500$ kHz, Eq. S15 predicts $N\simeq 1-5$ TLS. Note that the faster than expected $\gamma_{2}$ is also consistent with the fact that we do not observe coherent oscillations between TLS and the mechanical resonance.

We use Monte Carlo error propagation to determine the uncertainties for all values extrapolated from $\bar{n}$ data. First, we extract $P(n)$ from the Ramsey signal $S(t)$ by nonlinear least squares regression, which returns an estimate for the $P(n)$ uncertainties from the fit parameters’ covariance matrix. To determine the uncertainty on each calculated phonon number $\bar{n}$, we randomly resample $P(n)$ using the statistical uncertainty in each Fock level population to generate normally-distributed, zero-centered random noise. We repeat this process for 2,000 iterations and calculate $\bar{n}$ at each iteration to build a distribution of values. The reported uncertainty represents the standard deviation of this distribution. The same procedure is applied to determine the uncertainty in extracted values at each level of extrapolation; for fitted values, the resampled data are re-fit at each iteration. Fig. S4 shows the distributions generated by this procedure for fitted $T_{\text{2m}}$ in a displacement interferometry experiment as well as $\kappa_{1}$ and $\kappa_{2}$ from a representative nonlinear ringdown measurement.