Extracting the lifetime of a synthetic two-level system

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][4][5][6][7][8][9][10][11][12][13][14][15] whose potential energy is pumped at frequency f p close to twice its resonance frequency f 0 , i.e. at f p ≈ 2f 0 . If the modulation strength λ exceeds a threshold λ th , the device responds with oscillations locked to f d ≡ f p /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][19][20][21][22][23][24][25][26][27], or of other complex systems that can be mapped onto the same framework [28][29][30][31][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][37][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 * Corresponding author: eichlera@ethz.ch 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 "weak". Such a setting allows for situations with negligible backaction 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][41][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 doublyclamped 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, V bias = 10 V [47]. We use a Zurich Instruments HF2LI lock-in amplifier to apply the driving voltage V in and to measure the resonator displacement x ∝ V out = u cos(ωt) − v sin(ωt) with quadrature amplitudes u and v. For convenience, we drop the proportionality factor between x and V out and identify in the 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 V bias to the beam, and to capacitively drive and read-out the voltage signal Vout = u cos(ωt)− v sin(ω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 ∆ = f d − f0 with Vin = 0.4 V. 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-1 2 states. Each sweep contains 300 points measured within 685 seconds. (c) Switching between the phase states observed in v as a function of time with ∆ = 0 Hz, Vin = 0.4 V, and σV = 1 V. A dotted line represent the threshold between the phase states.
Our mechanical resonator can be described by the nonlinear equation of motion (in units of the measured electrical signal) Here, dots indicate time derivatives, ω 0 /2π = f 0 = 439.56 kHz is the resonance frequency, α = 1.47 × 10 18 V −2 s −2 the coefficient of the Duffing nonlinearity, γ = ω 0 /Q = 770 Hz the resonator relaxation rate, and Q = 3580 the quality factor of the resonator. The potential energy term (∝ x) is pumped with the parametric modulation depth λ = 2V in /(V th Q) at the angular frequency 2ω d = 4πf d , and where V th = 320 mV is the voltage threshold for parametric oscillations for the case f d = f 0 (demodulation frequency). The potential modulation arises because the electrostatic force due to V in pulls the beam closer towards 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 δf 0 ∝ V in . The force term ξ in Eq. 1 represents a fluctuating thermal bath [see Supplemental Material (SM) for details].   shows the v-quadrature response of the resonator during two sweeps of f d from positive to negative detuning ∆ ≡ f d − f 0 . Close to ∆ = 50 Hz, 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 Z 2 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., v 1 = −v 2 (and u 1 = −u 2 ) [27,[50][51][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 behaviour during switching events, we perform measurements with a high temporal resolution. In Fig. 2(a)-(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 10-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.
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 N switch during the full time T , and therefore also Γ = N switch /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 SM), but it cannot be removed entirely.
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, where the lifetime τ corresponds to the characteristic time scale 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) = Ae −2∆tΓ under the assumption of stationarity and ergodicity (not shown). Crucially, the autocorrelation is intimately related to the Allan variance (see SM 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 formalism [57] to the measured values, In this notation, are sums over the measured v values (or u values) and ... i denotes the mean over i, running from i = 1 to i = N − 2τ 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 obtain [45]: Hence, the maximum of σ 2 Allan (τ A ) should occur around the value τ A ≈ τ = Γ −1 . In Fig. 3(c), we indeed find a peak at the expected value, yielding Γ ≈ 4 Hz. 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 V in < 0.4 V, 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 SM. 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 threedimensional protein folding or other chemical reactions. For advanced applications in the future, the resonator networks could be realized through bilinear, resonant coupling of several KPOs [25,53] (see SM 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.
Data availability statement: The data that support the findings of this study are available from the corresponding author upon reasonable request.
The authors have no conflicts to disclose.

I. SUPPLEMENTARY MATERIAL
The supplementary material contains 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.

ACKNOWLEDGMENTS
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 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 grants (CRSII5 177198/1) and (PP00P2 190078), and by the Deutsche Forschungsgemeinschaft (DFG) -project number 449653034.

PROBABILITY DENSITY
In this section, we briefly detail our analysis steps yielding the probability density shown in Fig. 2(c) in the main text. First, we find the deterministic equations of motion for the "slow" quadrature amplitudes u and v. To this end, we use the so-called van der Pol transformation and lowest-order Krylov-Bogoliubov averaging method [1][2][3], to replace the full time-dependent equation of motion [cf. Eq. (1) in the main text] by time-independent averaged equations of motion. The white noise term ξ is transformed in a similar way using the method described in [4,5] which yieldsu where ω d = 2πf d is the angular demodulation frequency and X 2 = u 2 + v 2 . Ξ u and Ξ v are new independent white noise processes of strength σ. Equations (S1) provide a good description of the system when λ, γ/ω 0 , (α/ω 2 0 )x 2 , and ασ 2 /ω 2 0 are small. We use the stochastic differential equations (S1) to derive the corresponding Fokker-Planck equation [6,7] that describes the time evolution of the probability density p(u, v, t): We then solve this partial differential equation (PDE) numerically for the steady state (ṗ = 0) which is reached after long times and plot the outcome for the experimental parameters in Fig. 2(c).

ALLAN VARIANCE FOR A SYSTEM DOMINATED BY TELEGRAPH NOISE
In this section, we provide a short derivation of the Allan variance of a system subject to telegraph noise [cf. Eq. (5) in the main text]. The Allan variance can be calculated via the autocorrelation function. For a random telegraph signal x which switches between B and −B the autocorrelation function is given by [8]: where 1/τ is the mean rate of transitions. The Allan variance is obtained by the following expectation value for arXiv:2112.03357v2 [cond-mat.mes-hall] 26 Oct 2022 2 the cumulative amplitudes y(t) = t 0 x(t )dt We can now look for the maximum of the obtained expression by setting the first derivative equal to 0, i.e. dσ 2 A /dτ A ≡ 0. Expanding the numerator in τ A around τ up to forth order, we find that the maximum is located at τ A ≈ 0.946τ ≈ τ .

CALCULATION OF POWER SPECTRAL DENSITY AND AUTOCORRELATION
The data was recorded as two quadratures u and v that result from the signal demodulation at half the parametric drive frequency. The quadratures were then digitized at a sampling rate of 899.46 Sa/s. The total length of the time trace was 899.94 s with a total number of 809472 recorded points. The two quadratures were then transformed into the spectral domain by using the complex (u+iv) Welch power spectral density function in Python. Afterwards, the resulting double-sided power spectral density was converted to a single sided power spectral density. As a last step, Eq. 2 in the main text was fitted to the data.
The autocorrelation AC (T ) for correlation times T = nδt where n is an integer and δt is the sampling time, is calculated from the measured quadrature data points u k and v k in the following way. Let a k = u k +iv k be the complex data points, then is the standard deviation, where N the total number of data points. Here (C) denotes the complex conjugate and Re(C) the real part of a complex value C. The autocorrelation is then: We fit Eq. (S3) to the resulting AC (T ) to extract the switching rate Γ = 1/τ .

SWITCHING RATE AS A FUNCTION OF σV
Below, we show the switching rate that was estimated with all methods as a function of noise strength in a range where a two-level system exists. The parametric drive strength for this set was chosen to be V in = 0.436 V measured on resonance at f = 0.438 MHz. Switching rate as a function of noise amplitude σV for ∆ = 0 and Vin = 0.436 V estimated using simple average-based threshold (blue square), circle-based threshold (circles), power spectral density of telegraph noise (triangle), autocorrelation (star) and Allan deviation (hollow circle). Note that for low values of noise (roughly σV < 0.6 V) the larger discrepancies between the four rate estimation methods (circle-based, PSD, autocorrelation and Allan deviation) appear due to a small switching rate (small number of switching events).
We next comment on small deviations between switching rates estimated using different methods in the Fig. S2. In the parametron, we have two attractors in a rotating two-dimensional phase-space. One can think of the situation to be similar to a double-well potential well along one phase-space axis, alongside a single confining well along the perpendicular axis. As such, the system exhibits "fast" oscillations in the basins of phase-space as well as "slower switching" event where a phase-slip occurs and the oscillator moves from one phase state to the other. This is the reason behind the fact that the noise exhibits more involved dynamics than in a 1D double-well potential.
Our study here focuses on comparing methods by which to extract the different fluctuation sources, and specifically estimate the slow activation dynamics. To this end, we use the phenomenological observation that the activation rates are slower than those around the attractor, and evaluate the PSD, auto-correlation, and Allan-variance based on the assumption of having solely telegraph noise. We search for clear signatures of the slow dynamics, where we can differentiate them from the fast ones.

COUPLING OF MULTIPLE MICROMECHANICAL RESONATORS
Resonator networks for advanced applications could be realized through bilinear, resonant coupling of several KPOs [9,10]. For micromechanical resonators such as those studied here, bilinear coupling can be achieved in multiple ways, the most common being: pairwise (I) capacitive (II) inductive, (III) optical or, (IV) mechanical coupling, or (V) indirect all-to-all coupling through a separate radio-frequency cavity. Commonly, in the design of networks, the 4 most complex task will be the design of multiple connections with tunable interaction strength. The following classes of coupling architectures can be envisaged: Proximity coupling. Exchange of photons/phonons between resonators leads to nearest-neighbor coupling. For mechanical resonators, vibrations transmitted evanescently through the substrate or electrostatic interactions through the vacuum are sources of proximity coupling. The coupling can be set via geometric design but generally not tuned in the finished device.
Connector lines. Electrical striplines or optical waveguides that are proximity coupled to two separate resonators can act as a relay to achieve coupling between distant devices. An added advantage is that variable coupling can be obtained through changes in the transmission coefficient of the line. For instance, electrical connector lines can contain an inductive element or a tunable resonator whose transmission can be tuned.
The Lechner-Hauke-Zoller (LHZ) model. All-to-all coupling between more than two or three devices poses a severe challenge due to the necessity of a complex network of connector lines. In the minor embedding scheme [11,12], this problem is partially mitigated by creating chains of strongly coupled KPOs that are weakly coupled to other chains, where each chain effectively corresponds to a single logic unit. The LHZ model, on the other hand, presents a more radical possibility to build an all-to-all network using individually adjustable couplings, but without the complexity of a physical connector network [13,14]. In the LHZ model, the pairwise orientation of network states is encoded in additional parametron devices, and there is no need for physical all-to-all connections.
Nondegenerate resonators with parametric coupling. Many implementations of KPOs rely on degenerate parametric devices, i.e., they all have nominally identical resonance frequencies ω i = ω j = ω 0 and driving frequencies at 2ω 0 . In general, nondegenerate resonators have negligible coupling. However, three-wave mixing [15,16] offers a method for sizable and tunable coupling between nondegenerate parametrons; devices with different resonance (ω i = ω j for all pairs) and driving frequencies (2ω i , 2ω j , etc.) are coupled to each other through a common medium, for instance a low-Q optical cavity [17] that is proximity-coupled to all resonators. Modulating the cavity field at the frequency differences (ω i − ω j ) generates parametric frequency up-and down-conversion between device pairs [18,19]. The parametrons are thus in effect resonantly coupled.