The Kerr Parametric Oscillator (KPO) is a nonlinear resonator system that is often described as a synthetic two-level system. In the presence of noise, the system switches between two states via a fluctuating trajectory in phase space, instead of following a straight path. The presence of such fluctuating trajectories makes it hard to establish a precise count or even a useful definition, of the “lifetime” of the state. Addressing this issue, we compare several rate counting methods that allow to estimate a lifetime for the levels. In particular, we establish that a peak in the Allan variance of fluctuations can also be used to determine the levels' lifetime. Our work provides a basis for characterizing KPO networks for simulated annealing where an accurate determination of the state lifetime is of fundamental importance.

Synthetic two level systems (TLSs) generated in driven nonlinear resonators have recently caught a significant attention in the physics community.1,2 A particularly prominent example is the Kerr Parametric Oscillator (KPO, also known as parametron)3–15 whose potential energy is pumped at frequency fp close to twice its resonance frequency f0, i.e., at fp2f0. If the modulation strength λ exceeds a threshold λth, the device responds with oscillations locked to fdfp/2 within a certain detuning range. This well-known “period doubling” of the response relative to the pump gives rise to two stable “phase states” with the same amplitude but separated by a phase difference of π. The phase states can be used to encode the two polarization states (up/down) of a classical spin. This analogy leads to the idea of using networks of coupled KPOs to build noisy intermediate-scale quantum (NISQ) machines.16,17 These machines can simulate the dynamics of mathematical problems that overwhelm traditional computers, such as the ground state of an Ising Hamiltonian,18–27 or of other complex systems that can be mapped onto the same framework.28–32 

An important quantity for many applications of TLSs is their lifetime τ.33 It is the typical time spent on a level before the interaction with an environment induces a (seemingly) spontaneous “jump” from one state to the other. The rates of environmental noise-induced switching have previously been investigated for different systems, such as trapped electrons,34 cold atoms,35 micromechanical systems,36–38 and analog electronic circuits.39 

The situation is more subtle for a KPO. Here, the synthetic levels are formed by coherent bosonic states forming attractors in phase space. These attractors are not separated by an energy gap but by a phase gap.11 When switches occur on a slow timescale (relative to the resonator relaxation time) and follow narrow channels in phase space, the fluctuations are termed as “weak.” Such a setting allows for situations with negligible back action where the fluctuations during a single switch can be observed. Currently, however, there exist very few studies of the fascinating physics unfolding during individual switches.37,40–42

In this work, we study a classical micromechanical KPO and investigate its switching rates in the presence of weak fluctuations. We invoke and compare several methods previously used to characterize the rates of charge and parity state switching in cooper pair boxes and superconducting qubits.43,44 Furthermore, we propose a method to calculate the switching rate that is based on the Allan variance of the resonator displacement.45 In the final part of the paper, we compare all methods and find good agreement between several (but not all) of them.

Our KPO consists of a micro-electromechanical resonator (MEMS) in a room-temperature setup schematically shown in Fig. 1(a). The resonator is a doubly clamped beam, with the length of 200 μm, width 3 μm, and 60 μm in thickness with a lumped mass of 25.4 ng made from highly doped single crystal silicon and fabricated in a wafer-scale encapsulation process.46 Electrodes on both sides separated from the conducting beam with a gap 1 μm enable capacitive driving and sensitive detection of oscillations in the presence of a bias voltage, Vbias=10V.47 We use a Zurich Instruments HF2LI lock-in amplifier to apply the driving voltage Vin and to readout the resonator displacement xVout=ucos(ωt)vsin(ωt) with quadrature amplitudes u and v. For convenience, we drop the proportionality factor between x and Vout and identify in the following xVout.10 

FIG. 1.

Experimental setup and phase states of the KPO. (a) A Zurich Instruments HF2LI lock-in amplifier is used to apply a bias voltage Vbias to the beam and to capacitively drive and readout the voltage signal Vout=ucos(ωt)vsin(ωt) generated by the displacement x of the resonator. (b) Measured out-of-phase response v of the resonator to parametric driving as a function of detuning Δ=fdf0 with Vin=0.4V. Bright and dark dots correspond to different sweeps that showcase the amplitude-degenerate phase states of the KPO that can be interpreted as a synthetic TLS, e.g., spin-12 states. Each sweep contains 300 points measured within 685 s. (c) Switching between the phase states observed in v as a function of time with Δ = 0 Hz, Vin=0.4 V), and σV=1V. A dotted line represent the threshold between the phase states.

FIG. 1.

Experimental setup and phase states of the KPO. (a) A Zurich Instruments HF2LI lock-in amplifier is used to apply a bias voltage Vbias to the beam and to capacitively drive and readout the voltage signal Vout=ucos(ωt)vsin(ωt) generated by the displacement x of the resonator. (b) Measured out-of-phase response v of the resonator to parametric driving as a function of detuning Δ=fdf0 with Vin=0.4V. Bright and dark dots correspond to different sweeps that showcase the amplitude-degenerate phase states of the KPO that can be interpreted as a synthetic TLS, e.g., spin-12 states. Each sweep contains 300 points measured within 685 s. (c) Switching between the phase states observed in v as a function of time with Δ = 0 Hz, Vin=0.4 V), and σV=1V. A dotted line represent the threshold between the phase states.

Close modal

Our mechanical resonator can be described by the nonlinear equation of motion (in units of the measured electrical signal)

x¨+ω02[1λcos(2ωdt)]x+αx3+γẋ=ξ.
(1)

Here, dots indicate time derivatives, ω0/2π=f0=439.56kHz is the resonance frequency, α=1.47×1018V2s2 is the coefficient of the Duffing nonlinearity, γ=ω0/Q=770Hz is the resonator relaxation rate, and Q = 3580 is the quality factor of the resonator. The potential energy term (x) is pumped with the parametric modulation depth λ=2Vin/(VthQ) at the angular frequency 2ωd=4πfd where Vth=320μV is the voltage threshold for parametric oscillations for the case fd=f0 (demodulation frequency). The potential modulation arises because the electrostatic force due to Vin pulls the beam closer toward one electrode. The force is nonlinear, i.e., it grows stronger for small beam-electrode distances, which corresponds to a change in the overall spring constant that the beam experiences. As a consequence, the drive generates small frequency variations δf0Vin. The term ξ in Eq. (1) represents a fluctuating thermal bath (see the supplementary material for details.).

Figure 1(b) shows the v-quadrature response of the resonator during two sweeps of fd from positive to negative detuning Δfdf0. Close to Δ=50Hz, the response jumps from v = 0 to v=±50μV, marking a bifurcation point of the underlying nonlinear system. At the bifurcation, the resonator experiences a spontaneous Z2 symmetry breaking, also known as a period-doubling bifurcation or a discrete time-translation breaking.48,49 At this point, the resonator jumps to a positive or negative response with equal probability. The two responses belong to stable attractors (1 and 2) with opposite phases, i.e., v1=v2 (and u1=u2).27,50–52

To study switching between the phase states of our KPO, we apply white electrical noise ξ characterized by a standard deviation σV (over a bandwidth of 30 MHz) that causes the state of the resonator to fluctuate around its initial solution. If the fluctuations are large enough, they will occasionally carry the resonator across the threshold in the middle between the phase states. The resonator is then captured by the opposite attractor, corresponding to a switch of the synthetic TLS. Several such processes can be observed in Fig. 1(c). From this observation, it appears natural to attribute a lifetime to the inverse switching rate, τ=Γ1. However, calculating the switching rate is not straightforward due to the fluctuating trajectory.

For a deeper understanding of the system's transient behavior during switching events, we perform measurements with a high temporal resolution. In Figs. 2(a) and 2(b), we display a narrow time segment before, during, and after a single switch. We find many data points in the unstable zone between the two phase states. A ten-point moving average filter helps to visualize the trajectory of the system during the transition. The total switching time is roughly 10 ms, much longer than the lock-in integration time of 15 μs and the moving-average filter time of 700 μs. The measurement error of each data point is 3.7 μV, in agreement with the measured point-to-point fluctuations, but significantly smaller than the 10μV fluctuations visible on the 5 ms scale.

FIG. 2.

Phase space representation of states and switching. (a) u and v quadratures of a single phase state switch composed of 2170 data points measured with a 15 μs integration time at 14 391 samples per second. Bright dots and dark lines correspond to raw data and to a 10-point moving average, respectively, which allows to reduce the influence of detection noise. A dashed line indicates the threshold between the phase states. Δ = 0 Hz, Vin=0.4 V, and σV=0.6V. (b) Phase space representation of the data in (a). White squares indicate the attractor points measured in the absence of noise, and a dashed line indicates the threshold between the phase states. (c) Probability density of the KPO steady state calculated with a numerical evolution of a Fokker–Planck description of the system. Driven by classical force noise, the system explores its phase space stochastically. Dark blue indicates a low probability that the KPO visits a position in phase space within a finite time, and bright yellow indicates a high probability (scale not normalized).

FIG. 2.

Phase space representation of states and switching. (a) u and v quadratures of a single phase state switch composed of 2170 data points measured with a 15 μs integration time at 14 391 samples per second. Bright dots and dark lines correspond to raw data and to a 10-point moving average, respectively, which allows to reduce the influence of detection noise. A dashed line indicates the threshold between the phase states. Δ = 0 Hz, Vin=0.4 V, and σV=0.6V. (b) Phase space representation of the data in (a). White squares indicate the attractor points measured in the absence of noise, and a dashed line indicates the threshold between the phase states. (c) Probability density of the KPO steady state calculated with a numerical evolution of a Fokker–Planck description of the system. Driven by classical force noise, the system explores its phase space stochastically. Dark blue indicates a low probability that the KPO visits a position in phase space within a finite time, and bright yellow indicates a high probability (scale not normalized).

Close modal

Our observation depicted in Fig. 2 demonstrates that activated switches between the phase states are not deterministic but include prominent random elements. For instance, in the phase-space representation of the switch in Fig. 2(b), we can clearly see that the system performs a winding path close to the origin. In our device, the fluctuations generally have a slight preference for counterclockwise rotations around the phase states and clockwise ones around the origin. This can be explained by the combination of the drive and the nonlinearity, which leads to an effective detuning of the fluctuations from the lock-in amplifier clock.53 In the corresponding Fokker–Planck steady-state calculation presented in Fig. 2(c), we, therefore, find a channel with a significant probability density between the phase states.

These visualizations of the fluctuating trajectories expose a fundamental problem in estimating the lifetime τ: since transitions follow no straight lines, they can cross any point in phase space multiple times during a single switching event. An example of this can be observed in Fig. 2(b), where the averaged (dark) trajectory crosses the dotted threshold line from bottom to top, describes a clockwise winding that traverses back across the threshold, and finally crosses the line a third time before completing the switch. A simple counting algorithm will in this case register three crossing events during a single switch. In general, any counting method based on a simple threshold (such as a line) will, therefore, overestimate the switching number Nswitch during the full time T, and therefore, also Γ=Nswitch/T. This problem has been known since a long time.

The problem of overestimating the switching count can be reduced by defining multiple thresholds that have to be crossed in a particular order to constitute an event. In Fig. 3(a), we demonstrate this with the example of two circles in phase space. The count is increased by one each time a circular threshold is left and the opposing circle is entered. This method is less sensitive to small fluctuations, but it requires a subjective measure that impacts the estimated Γ, in our case, the radii of the circular thresholds. Calibrating the measured switching rate Γ as a function of the radius helps to reduce this degree of arbitrariness (see the supplementary material), but it cannot be removed entirely.

FIG. 3.

Methods used to estimate the phase state lifetime. All plots show the same 15 m dataset, an extract of which is shown in Fig. 1(c). Data were recorded at 899 samples per second with an integration time τ = 143 μs and with Δ = 0 Hz, Vin=0.4 V, and σV=1V. (a) Phase space representation of the two phase states and of switching between them (cf. Fig. 2). The white squares indicate the attractors measured in the absence of noise, and the dotted line and circles indicate different threshold methods outlined in the text. The radius of the circle in this case was set to be 70% of the distance from the center of the circle to the origin of the coordinate system. The estimated activation rates are Γline13±0.1Hz and Γcirc4.35±0.07Hz, with standard deviations calculated assuming Poisson statistics of the jumps. (b) PSD analysis of the fluctuations in terms of a telegraph noise model, cf. Eq. (2), yielding a fit result Γpsd3.60±0.01Hz with a fit value F = 5.86 × 10–5 V2. Bright and dark lines correspond to the measured data and to the fit, respectively. (c) Allan variance of the measured fluctuations (bright), with a maximum at ΓAllan1/τ=4.00±0.08Hz, where the precision is limited by the separation of points in τ. A dark line is the function expected (with arbitrary vertical scaling) for pure telegraph noise with a mean switching rate of 4 Hz [see Eq. (5)].

FIG. 3.

Methods used to estimate the phase state lifetime. All plots show the same 15 m dataset, an extract of which is shown in Fig. 1(c). Data were recorded at 899 samples per second with an integration time τ = 143 μs and with Δ = 0 Hz, Vin=0.4 V, and σV=1V. (a) Phase space representation of the two phase states and of switching between them (cf. Fig. 2). The white squares indicate the attractors measured in the absence of noise, and the dotted line and circles indicate different threshold methods outlined in the text. The radius of the circle in this case was set to be 70% of the distance from the center of the circle to the origin of the coordinate system. The estimated activation rates are Γline13±0.1Hz and Γcirc4.35±0.07Hz, with standard deviations calculated assuming Poisson statistics of the jumps. (b) PSD analysis of the fluctuations in terms of a telegraph noise model, cf. Eq. (2), yielding a fit result Γpsd3.60±0.01Hz with a fit value F = 5.86 × 10–5 V2. Bright and dark lines correspond to the measured data and to the fit, respectively. (c) Allan variance of the measured fluctuations (bright), with a maximum at ΓAllan1/τ=4.00±0.08Hz, where the precision is limited by the separation of points in τ. A dark line is the function expected (with arbitrary vertical scaling) for pure telegraph noise with a mean switching rate of 4 Hz [see Eq. (5)].

Close modal

To avoid overcounting and subjective dependencies, it is desirable to extract Γ from a method that does not require thresholds at all. Interestingly, the parity lifetime of superconducting qubits can be determined via their charge-parity power spectral density (PSD).43,54,55 Assuming that the switching is dominated by telegraph noise, the PSD of v of our KPO can be fitted to a Lorentzian function,

PSD(f)=2F2τ4+(2πfτ)2,
(2)

where the lifetime τ corresponds to the characteristic timescale between level switching events, and F is a constant related to the measurement fidelity.55 In this case, the lifetime or the switching rate is related to the width of the spectral peak in the frequency domain.56 To make the estimate quantitative in Fig. 3(b), we fit the measured displacement power spectral density with Eq. (2), yielding a third estimation for Γ=1/τ lifetime. The method can also be applied after a Fourier transform by fitting the sliding average autocorrelation with the function AC(Δt)=Ae2ΔtΓ under the assumption of stationarity and ergodicity (not shown).

Crucially, the autocorrelation is intimately related to the Allan variance (see the supplementary material for the derivation). Originally invented to characterize the fidelity of clocks, the Allan variance measures the frequency fluctuations of a resonator as a function of integration time τA. As we are interested in the time τ over which the typical fluctuations of u (or v) of our KPO are maximal, we apply the Allan variance formalism57 to the measured values,

σAllan2(τA)=12τA2(ai,22ai,1+ai,0)2i.
(3)

In this notation,

ak,l=m=0k+τAl/δtv(m)
(4)

are sums over the measured v values (or u values) and i denotes the mean over i, running from i = 1 to i=N2τA/δt, where N is the total number of data points and δt is the sampling time. Assuming that the signal is dominated by telegraph-like switching with lifetime τ and amplitude B, we obtain45 

σAllan2(τA)=B24τA/τ+e4τA/τ4e2τA/τ+34τA2/τ2.
(5)

Hence, the maximum of σAllan2(τA) should occur around the value τAτ=Γ1. In Fig. 3(c), we, indeed, find a peak at the expected value, yielding Γ4Hz. In contrast to the PSD method, the Allan variance method does not necessarily require a fitting process, as the peak can be read off directly and is easy to interpret even in the presence of noise.

We compare the results of the different methods in Fig. 4. We find excellent agreement between four out of five of the methods for values of Γ varying over more than two orders of magnitude. The only method that we wish to discard from this comparison is the simple line threshold approach, which consistently overestimates the count rate as expected from the discussion above. The method using two circles for thresholding overestimates Γ slightly for Vin<0.4V, where the separation between the attractors is small and the “clouds” start to overlap significantly, cf. the example in Fig. 3(a). Additional comparison as a function of the noise strength σV can be found in the supplementary material.

FIG. 4.

Comparison of results for Γ obtained with different rate estimation methods. Switching rate as a function of parametric drive amplitude Vin for Δ = 0 and σV=1V estimated using simple line-based thresholding (blue square), circle-based thresholding (filled circle), power spectral density of telegraph noise (triangle), autocorrelation (star), and Allan deviation (hollow circle). The radius of the circle method in this case was set to be 50% of the distance from the center of the circle to the origin of the coordinate system.

FIG. 4.

Comparison of results for Γ obtained with different rate estimation methods. Switching rate as a function of parametric drive amplitude Vin for Δ = 0 and σV=1V estimated using simple line-based thresholding (blue square), circle-based thresholding (filled circle), power spectral density of telegraph noise (triangle), autocorrelation (star), and Allan deviation (hollow circle). The radius of the circle method in this case was set to be 50% of the distance from the center of the circle to the origin of the coordinate system.

Close modal

We emphasize that there is no fundamental reason why the estimators we obtain should be identical at all. The surprisingly good agreement between most of the estimators confirms that the notion of a lifetime τ is useful to characterize the switching between phase states in a KPO, where a parametric pump generates a synthetic potential landscape.49 This approach may be useful in other systems where multi-stable potentials in dimensions higher than one are present, such as three-dimensional protein folding or other chemical reactions. For advanced applications in the future, the resonator networks could be realized through bilinear, resonant coupling of several KPOs25,53 (see the supplementary material for details). For MEMS, such as those studied here, bilinear coupling can be achieved in multiple ways, such as pairwise capacitive, inductive, optical, or mechanical coupling, or indirect all-to-all coupling through a separate radio frequency cavity.

See the supplementary material for theory derivations for the probability density, the Allan variance, and the autocorrelation of telegraph noise; experimental demonstrations of the dependence of the extracted switching rate on the circle threshold radius and on the noise strength; and a short summary of various coupling methods for parametric oscillators.

Fabrication was performed in nano@Stanford labs, which are supported by the National Science Foundation (NSF) as part of the National Nanotechnology Coordinated Infrastructure under Award No. ECCS-1542152, with the support from the Defense Advanced Research Projects Agency's Precise Robust Inertial Guidance for Munitions (PRIGM) Program, managed by Ron Polcawich and Robert Lutwak. This work was further supported by the Swiss National Science Foundation through Grant Nos. CRSII5_177198/1 and PP00P2 190078 and a Deutsche Forschungsgemeinschaft (DFG) - Project No. 449653034.

The authors have no conflicts to disclose.

Gabriel Margiani: Data curation (equal); Investigation (equal); Resources (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Nicholas E. Bousse: Resources (equal); Writing – review & editing (equal). Thomas W. Kenny: Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). Oded Zilberberg: Formal analysis (equal); Supervision (equal); Writing – original draft (equal). Deividas Sabonis: Supervision (equal); Writing – original draft (equal). Alexander Eichler: Conceptualization (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal). Sebastián Guerrero: Resources (equal); Software (equal); Writing – review & editing (equal). Toni Louis Heugel: Data curation (equal); Formal analysis (equal); Writing – review & editing (equal). Christian Marty: Resources (equal); Writing – review & editing (equal). Raphael Pachlatko: Resources (equal); Writing – review & editing (equal). Thomas Gisler: Investigation (equal); Methodology (equal); Writing – review & editing (equal). Gabrielle D. Vukasin: Resources (equal); Writing – review & editing (equal). Hyun-Keun Kwon: Resources (equal); Writing – review & editing (equal). James M Lehto Miller: Conceptualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
D.
Gottesman
,
A.
Kitaev
, and
J.
Preskill
,
Phys. Rev. A
64
,
012310
(
2001
).
2.
M. H.
Devoret
and
R. J.
Schoelkopf
,
Science
339
,
1169
(
2013
).
3.
I.
Mahboob
and
H.
Yamaguchi
,
Nat. Nanotechnol.
3
,
275
(
2008
).
4.
C. M.
Wilson
,
T.
Duty
,
M.
Sandberg
,
F.
Persson
,
V.
Shumeiko
, and
P.
Delsing
,
Phys. Rev. Lett.
105
,
233907
(
2010
).
5.
A.
Eichler
,
J.
Chaste
,
J.
Moser
, and
A.
Bachtold
,
Nano Lett.
11
,
2699
(
2011
).
6.
A.
Eichler
,
T. L.
Heugel
,
A.
Leuch
,
C. L.
Degen
,
R.
Chitra
, and
O.
Zilberberg
,
Appl. Phys. Lett.
112
,
233105
(
2018
).
7.
J.
Gieseler
,
B.
Deutsch
,
R.
Quidant
, and
L.
Novotny
,
Phys. Rev. Lett.
109
,
103603
(
2012
).
8.
Z.
Lin
,
K.
Inomata
,
K.
Koshino
,
W. D.
Oliver
,
Y.
Nakamura
,
J. S.
Tsai
, and
T.
Yamamoto
,
Nat. Commun.
5
,
4480
(
2014
).
9.
S.
Puri
,
S.
Boutin
, and
A.
Blais
,
npj Quantum Inf.
3
,
18
(
2017
).
10.
Z.
Nosan
,
P.
Märki
,
N.
Hauff
,
C.
Knaut
, and
A.
Eichler
,
Phys. Rev. E
99
,
062205
(
2019
).
11.
M.
Frimmer
,
T. L.
Heugel
,
Z.
Nosan
,
F.
Tebbenjohanns
,
D.
Hälg
,
A.
Akin
,
C. L.
Degen
,
L.
Novotny
,
R.
Chitra
,
O.
Zilberberg
, and
A.
Eichler
,
Phys. Rev. Lett.
123
,
254102
(
2019
).
12.
A.
Grimm
,
N. E.
Frattini
,
S.
Puri
,
S. O.
Mundhada
,
S.
Touzard
,
M.
Mirrahimi
,
S. M.
Girvin
,
S.
Shankar
, and
M. H.
Devoret
,
Nature
584
,
205
(
2020
).
13.
S.
Puri
,
A.
Grimm
,
P.
Campagne-Ibarcq
,
A.
Eickbusch
,
K.
Noh
,
G.
Roberts
,
L.
Jiang
,
M.
Mirrahimi
,
M. H.
Devoret
, and
S. M.
Girvin
,
Phys. Rev. X
9
,
041009
(
2019
).
14.
J. M.
Miller
,
D. D.
Shin
,
H.-K.
Kwon
,
S. W.
Shaw
, and
T. W.
Kenny
,
Phys. Rev. Appl.
12
,
044053
(
2019
).
15.
D.
Ryvkine
and
M. I.
Dykman
,
Phys. Rev. E
74
,
061118
(
2006
).
17.
T.
Albash
and
D. A.
Lidar
,
Rev. Mod. Phys.
90
,
015002
(
2018
).
19.
H.
Goto
,
Z.
Lin
, and
Y.
Nakamura
,
Sci. Rep.
8
,
7154
(
2018
).
20.
R.
Rota
,
F.
Minganti
,
C.
Ciuti
, and
V.
Savona
,
Phys. Rev. Lett.
122
,
110405
(
2019
).
21.
B.
Heim
,
T. F.
Rønnow
,
S. V.
Isakov
, and
M.
Troyer
,
Science
348
,
215
(
2015
).
22.
I.
Mahboob
,
H.
Okamoto
, and
H.
Yamaguchi
,
Sci. Adv.
2
,
e1600236
(
2016
).
23.
T.
Inagaki
,
K.
Inaba
,
R.
Hamerly
,
K.
Inoue
,
Y.
Yamamoto
, and
H.
Takesue
,
Nat. Photonics
10
,
415
(
2016
).
25.
L.
Bello
,
M.
Calvanese Strinati
,
E. G.
Dalla Torre
, and
A.
Pe'er
,
Phys. Rev. Lett.
123
,
083901
(
2019
).
26.
Y.
Okawachi
,
M.
Yu
,
J. K.
Jang
,
X.
Ji
,
Y.
Zhao
,
B. Y.
Kim
,
M.
Lipson
, and
A. L.
Gaeta
,
Nat. Commun.
11
,
4119
(
2020
).
27.
M.
Calvanese Strinati
,
L.
Bello
,
A.
Pe'er
, and
E. G.
Dalla Torre
,
Phys. Rev. A
100
,
023835
(
2019
).
29.
S. E.
Nigg
,
N.
Lörch
, and
R. P.
Tiwari
,
Sci. Adv.
3
,
e1602273
(
2017
).
30.
T.
Inagaki
,
Y.
Haribara
,
K.
Igarashi
,
T.
Sonobe
,
S.
Tamate
,
T.
Honjo
,
A.
Marandi
,
P. L.
McMahon
,
T.
Umeki
,
K.
Enbutsu
,
O.
Tadanaga
,
H.
Takenouchi
,
K.
Aihara
,
K-i
Kawarabayashi
,
K.
Inoue
,
S.
Utsunomiya
, and
H.
Takesue
,
Science
354
,
603
(
2016
).
31.
H.
Goto
,
K.
Tatsumura
, and
A. R.
Dixon
,
Sci. Adv.
5
,
eaav2372
(
2019
).
32.
T.
Honjo
,
T.
Sonobe
,
K.
Inaba
,
T.
Inagaki
,
T.
Ikuta
,
Y.
Yamada
,
T.
Kazama
,
K.
Enbutsu
,
T.
Umeki
,
R.
Kasahara
 et al,
Sci. Adv.
7
,
eabh0952
(
2021
).
33.
P.
Krantz
,
M.
Kjaergaard
,
F.
Yan
,
T. P.
Orlando
,
S.
Gustavsson
, and
W. D.
Oliver
,
Appl. Phys. Rev.
6
,
021318
(
2019
).
34.
L. J.
Lapidus
,
D.
Enzer
, and
G.
Gabrielse
,
Phys. Rev. Lett.
83
,
899
(
1999
).
35.
K.
Kim
,
M.-S.
Heo
,
K.-H.
Lee
,
H.-J.
Ha
,
K.
Jang
,
H.-R.
Noh
, and
W.
Jhe
,
Phys. Rev. A
72
,
053402
(
2005
).
36.
J. S.
Aldridge
and
A. N.
Cleland
,
Phys. Rev. Lett.
94
,
156403
(
2005
).
37.
H. B.
Chan
and
C.
Stambaugh
,
Phys. Rev. Lett.
99
,
060601
(
2007
).
38.
W. J.
Venstra
,
H. J.
Westra
, and
H. S.
Van Der Zant
,
Nat. Commun.
4
(
1
),
2624
(
2013
).
39.
D. G.
Luchinsky
,
R. S.
Maier
,
R.
Mannella
,
P. V. E.
McClintock
, and
D. L.
Stein
,
Phys. Rev. Lett.
82
,
1806
(
1999
).
40.
H. B.
Chan
,
M. I.
Dykman
, and
C.
Stambaugh
,
Phys. Rev. Lett.
100
,
130602
(
2008
).
41.
I.
Mahboob
,
M.
Mounaix
,
K.
Nishiguchi
,
A.
Fujiwara
, and
H.
Yamaguchi
,
Sci. Rep.
4
,
4448
(
2015
).
42.
M. I.
Dykman
,
C. M.
Maloney
,
V. N.
Smelyanskiy
, and
M.
Silverstein
,
Phys. Rev. E
57
,
5202
(
1998
).
43.
K.
Serniak
,
M.
Hays
,
G.
De Lange
,
S.
Diamond
,
S.
Shankar
,
L.
Burkhart
,
L.
Frunzio
,
M.
Houzet
, and
M.
Devoret
,
Phys. Rev. Lett.
121
,
157701
(
2018
).
44.
K.
Serniak
,
S.
Diamond
,
M.
Hays
,
V.
Fatemi
,
S.
Shankar
,
L.
Frunzio
,
R.
Schoelkopf
, and
M.
Devoret
,
Phys. Rev. Appl.
12
,
014052
(
2019
).
45.
C. M.
Van Vliet
and
P. H.
Handel
,
Physica A
113
,
261
(
1982
).
46.
Y.
Yang
,
E. J.
Ng
,
Y.
Chen
,
I. B.
Flader
, and
T. W.
Kenny
,
J. Microelectromech. Syst.
25
,
489
(
2016
).
47.
J. M.
Miller
,
N. E.
Bousse
,
D. B.
Heinz
,
H. J. K.
Kim
,
H.-K.
Kwon
,
G. D.
Vukasin
, and
T. W.
Kenny
,
J. Microelectromech. Syst.
28
,
965
(
2019
).
48.
L.
Landau
and
E.
Lifshitz
,
Mechanics
(
Butterworth-Heinemann
,
1976
).
49.
T. L.
Heugel
,
M.
Oscity
,
A.
Eichler
,
O.
Zilberberg
, and
R.
Chitra
,
Phys. Rev. Lett.
123
,
124301
(
2019
).
50.
M. C.
Lifshitz
and
R.
Cross
, “
Nonlinear dynamics of nanomechanical and micromechanical resonators
,” in
Reviews of Nonlinear Dynamics and Complexity
(
Wiley-VCH
,
2009
), pp.
1
52
.
51.
M.
Dykman
and
M.
Krivoglaz
,
Sov. Phys. JETP
50
,
30
(
1979
).
52.
L.
Papariello
,
O.
Zilberberg
,
A.
Eichler
, and
R.
Chitra
,
Phys. Rev. E
94
,
022201
(
2016
).
53.
T. L.
Heugel
,
O.
Zilberberg
,
C.
Marty
,
R.
Chitra
, and
A.
Eichler
, arXiv:2103.02625 (
2021
).
54.
P.
Dutta
and
P. M.
Horn
,
Rev. Mod. Phys.
53
,
497
(
1981
).
55.
D.
Ristè
,
C.
Bultink
,
M.
Tiggelman
,
R.
Schouten
,
K.
Lehnert
, and
L.
DiCarlo
,
Nat. Commun.
4
(
1
),
1913
(
2013
).
56.
W.
Demtröder
,
Laser Spectroscopy
(
Springer
,
1973
), Vol.
5
.
57.

Supplementary Material