# Proliferation of unstable states
and their impact on stochastic out-of-equilibrium dynamics

Toni L. Heugel  Institute for Theoretical Physics, ETH Zürich, CH-8093 Zürich,
Switzerland.  R. Chitra  Institute for Theoretical Physics, ETH Zürich,
CH-8093 Zürich, Switzerland.  Alexander Eichler  Laboratory for Solid State
Physics, ETH Zürich, CH-8093 Zürich, Switzerland.  Quantum Center, ETH Zurich,
CH-8093 Zurich, Switzerland  Oded Zilberberg  Department of Physics,
University of Konstanz, D-78457 Konstanz, Germany

###### Abstract

Networks of nonlinear parametric resonators are promising candidates as Ising
machines for annealing and optimization. These many-body out-of-equilibrium
systems host complex phase diagrams of coexisting stationary states. The
plethora of states manifest via a series of bifurcations, including
bifurcations that proliferate purely unstable solutions, which we term “ghost
bifurcations”. Here, we demonstrate that the latter take a fundamental role in
the stocahstic dynamics of the system in the presence of noise. Specifically,
they determine the switching paths and the switching rates between stable
solutions. We demonstrate experimentally the impact of ghost bifurcations on
the noise-activated switching dynamics in a network of two coupled parametric
resonators.

##  I Introduction

Statistical physics has provided valuable insights into the phenomenon of
noise-induced switching between local energy minima, as exemplified by the
well-known Kramers double-well problem Kramers (1940); Hänggi  _et al._
(1990); Rondin  _et al._ (2017a). Such stochastic dynamics is highly relevant
for a wide variety of phenomena spanning protein folding Best and Hummer
(2006); Chung  _et al._ (2015), chemical reactions García-Müller  _et al._
(2008), as well as stability in mechanical Badzey and Mohanty (2005); Rondin
_et al._ (2017b) and electrical systems Fulton and Dunkleberger (1974);
Silvestrini  _et al._ (1988). While the noisy dynamics of equilibrium systems
have been extensively studied, that of systems driven far out of equilibrium
remains largely unexplored Dykman  _et al._ (1998); Dykman (2007); Tadokoro
_et al._ (2020).

An important class of out-of-equilibrium systems are driven systems. Such
systems are characterized by stationary oscillation states that manifest when
the conserving and nonconserving forces in the system are in balance. For
nonlinear systems, there may be several such stationary states that act as
attractors, just like potential wells do in equilibrium systems. In a rotating
frame, the resulting dynamics can resemble that of an equilibrium potential
landscape Hänggi  _et al._ (1990). Extending the analogy, noise can induce
stochastic switching between the attractors, and the switching rate can be
treated with the abstract notion of a potential activation barrier Dykman  _et
al._ (1998); Luchinsky  _et al._ (1999a); Lapidus  _et al._ (1999); Kim  _et
al._ (2005); Aldridge and Cleland (2005); Chan and Stambaugh (2007); Chan  _et
al._ (2008); Venstra  _et al._ (2013); Mahboob  _et al._ (2014); Margiani  _et
al._ (2022). However, it is important to emphasize that this activation
barrier is not related to a gap in free energy, but instead to a “phase gap”
that separates the attractors Frimmer  _et al._ (2019). Understanding how
often the system switches between attractors, and which path it selects during
the switch, requires different methods and can be cumbersome Hänggi  _et al._
(1990); Dykman  _et al._ (1998).

A paradigmatic example of a bistable out-of-equilibrium system is the Kerr
parametric oscillator (KPO) Ryvkine and Dykman (2006); Mahboob and Yamaguchi
(2008); Wilson  _et al._ (2010); Eichler  _et al._ (2011); Leuch  _et al._
(2016); Gieseler  _et al._ (2012); Lin  _et al._ (2014); Puri  _et al._
(2017a); Eichler  _et al._ (2018); Nosan  _et al._ (2019); Frimmer  _et al._
(2019); Grimm  _et al._ (2019); Puri  _et al._ (2019); Miller  _et al._
(2019); Dykman (2012); Eichler and Zilberberg (2023). Recently, networks of
coupled, driven KPOs have been proposed as a simulation platform to solve
complex problems optimally. Mahboob  _et al._ (2016); Inagaki  _et al._
(2016); Goto (2016); Puri  _et al._ (2017b); Nigg  _et al._ (2017); Dykman
_et al._ (2018); Okawachi  _et al._ (2020). Such networks typically possess a
large number of stationary states, analogous to an multi-well potential with
rich phase transitions Calvanese Strinati  _et al._ (2021); Bello  _et al._
(2019); Heugel  _et al._ (2022a, b). Studying activated switching between
states in the presence of fluctuations is crucial for understanding their
stability and lifetimes. It will influence how these networks are operated,
and it can also provide a characterization method that is unaffected by the
danger of local trapping Margiani  _et al._ (2023).

In this work, we address the physics of stochastic activation in a tractable
system of two strongly coupled, classical KPOs. The system possesses various
stable and unstable stationary states Heugel  _et al._ (2019, 2022a), and we
observe stochastic switching between two such states in the presence of
fluctuations. Surprisingly, the switching rate $\Gamma$ deviates significantly
from the exponential model expected for a single KPO Dykman  _et al._ (1998).
Seeking to explain this deviation, we calculate the dominant transition paths
between the states with the Onsager-Machlup function Lehmann  _et al._ (2003);
Wio (2013). Our analysis shows that new kinds of unstable states, termed ghost
states emerge in the two KPO system. These states do not manifest in the
stationary deterministic dynamics of the system. Interestingly, however, they
offer new transition paths and contribute significantly to the transition
paths and the corresponding transition rates. We thus identify a striking
example of out-of-equilibrium statistical physics in a nonlinear multistable
system. Our work paves the way for the exploration of larger systems,
especially in view of KPO networks as solvers for complex optimization tasks.

##  II System

In the following, we analyze noise-induced switching dynamics using an
experimental setup composed of two electrical KPOs with capacitive coupling
Nosan  _et al._ (2019); Heugel  _et al._ (2022a). Each KPO consists of a coil
with inductance $L$ and a diode that provides a nonlinear capacitance $C$, cf.
Appendix A. The resonance frequency of each KPO can be tuned by applying a DC
voltage across the diode. We drive and measure the resonators inductively
through auxiliary coils. Our electrical circuits are well described by the
following coupled equations of motion

$\ddot{x}_{i}+\omega_{i}^{2}\left[1-\lambda\cos\left(2\omega_{d}t\right)\right]x_{i}+\alpha_{i}x_{i}^{3}\\\
+\gamma_{i}\dot{x}_{i}-\sum_{j\neq i}J_{ij}x_{j}=\xi_{i}(t)\,,$ (1)

where dots indicate time derivatives, $x_{i}=u_{i}\cos(\omega
t)-v_{i}\sin(\omega t)$ is the measured voltage with quadrature amplitudes
$u_{i}$ and $v_{i}$, $\omega_{i}=2\pi f_{i}$ is the angular eigenfrequency,
and $J_{ij}$ ($i\neq j$) denotes the inter-circuit linear coupling strength.
Each resonator has an effective Duffing (Kerr) nonlinearity with coefficient
$\alpha_{i}$ and a damping rate $\gamma_{i}=\omega_{i}/Q_{i}$, with $Q_{i}$
the quality factor. Our resonators are constructed and tuned to be (nearly)
identical in their bare characteristics, $\omega_{i}\approx\omega_{0}=2\pi
f_{0}$. The same parametric pumping with angular modulation frequency
$2\omega_{d}=4\pi f_{d}\approx 2\omega_{0}$ and modulation depth
$\lambda\propto U_{d}$ is applied to all resonators.

Crucially, beyond a frequency-dependent driving threshold $U_{\mathrm{th}}$,
the KPO has exactly two stable solutions that we refer to as ‘phase states’,
which have identical amplitudes but differ in phase by $\pi$ McLachlan (1951);
Lifshitz (2009); Eichler and Zilberberg (2023). These are the two attractors
of a single KPO in a frame rotating at $\omega_{d}$. In order to induce
switching events between the attractors, we add an artificial noise $\xi_{i}$
generated by a fluctuating voltage with white power spectral density $S_{n}$.
The noise $\xi_{i}$ simulates a thermal force noise with
$\langle\xi_{i}(t_{1})\xi_{j}(t_{2})\rangle=\varsigma^{2}\delta_{ij}\delta(t_{1}-t_{2})$
and power spectral density calibrated to be
$\varsigma^{2}=$4.93\text{\times}{10}^{-20}\text{\,}\mathrm{Hz}^{4}$S_{n}$,
see Appendix B.

In our experiments, we use a lock-in amplifier to measure the quadratures
$(u_{i},v_{i})$ which vary on timescales much longer than $1/\omega_{0}$. The
evolution of these quadratures is well captured in the rotating-frame picture
obtained by applying the averaging method Guckenheimer and Holmes (1990);
Papariello  _et al._ (2016); Heugel  _et al._ (2019). to our model. Equation
(1) then leads to the following slow-flow equations:

$\displaystyle\dot{u}_{i}$ $\displaystyle=-\frac{\gamma
u_{i}}{2}-\left(\frac{3\alpha}{8\omega_{d}}X_{i}^{2}+\frac{\omega_{0}^{2}-\omega_{d}^{2}}{2\omega_{d}}+\frac{\lambda\omega_{0}^{2}}{4\omega_{d}}\right)v_{i}+\frac{Jv_{j}}{2\omega_{d}}+\Xi_{u_{i}}\,,$
$\displaystyle\dot{v}_{i}$ $\displaystyle=-\frac{\gamma
v_{i}}{2}+\left(\frac{3\alpha}{8\omega_{d}}X_{i}^{2}+\frac{\omega_{0}^{2}-\omega_{d}^{2}}{2\omega_{d}}-\frac{\lambda\omega_{0}^{2}}{4\omega_{d}}\right)u_{i}-\frac{Ju_{j}}{2\omega_{d}}+\Xi_{v_{i}}\,,$
(2)

where $X_{i}^{2}=u_{i}^{2}+v_{i}^{2}$, and we have additive uncorrelated noise
terms $\Xi_{u_{i}}$, $\Xi_{v_{i}}$, whose power spectral densities are given
by $\sigma^{2}=~{}\varsigma^{2}/2\omega_{d}^{2}$ Khas’minskii (1966); Roberts
and Spanos (1986). As shown in previous works, the averaging method fully
captures the physics of our networks in the regime where $\lambda$,
$\gamma/\omega_{0}$, $J/\omega_{0}^{2}$, $(\alpha/\omega_{0}^{2})x_{i}^{2}$,
and $\alpha\varsigma^{2}/\omega_{0}^{5}$ are all of order $\epsilon$ with
$0<\epsilon\ll 1$ Nayfeh and Mook (2008); Khas’minskii (1966); Papariello  _et
al._ (2016); Eichler  _et al._ (2018). Here and in the following, we assume
identical dissipation rates $\gamma_{i}=\gamma$, nonlinearities
$\alpha_{i}=\alpha$, and coupling $J_{ij}=J$.

Figure 1: Noise-induced switching in a single KPO. (a) Time trace of the
rotating quadratures $u_{1}$ and $v_{1}$, showing noisy fluctuations around
and switching between the two phase states for $U_{d}=$3\text{\,}\mathrm{V}$$
and
$S_{n}=$2.75\text{\times}{10}^{-8}\text{\,}{\mathrm{V}}^{2}\text{\,}{\mathrm{Hz}}^{-1}$$
at $f_{d}=$2.681\text{\,}\mathrm{MHz}$$. Arrows indicate the position of the
phase states. (b) Heatmap of normalized measurement counts $c$ versus $u_{1}$
and $v_{1}$, showing two dominant attractors around the phase states (colored
dots). The data stem from a time trace measurement as in (a) taken over
$8\text{\,}\mathrm{s}$. Counts close to the origin arise due to inter-state
switching with rate $\Gamma$. The orange line shows the most probable
switching path through the origin (grey square) predicted from the Onsager-
Machlup formalism, cf. Eq. (3) and Appendix A. (c) Switching rate $\Gamma$ as
a function of $f_{d}$ with the same parameters as in (a) obtained from the
experiment (black dots) Margiani  _et al._ (2022). The orange line corresponds
to the approximate analytical result from Ref. Dykman  _et al._ (1998), cf.
Appendix C, and the orange dots to numerical minimization of Eq. (3).

##  III Transition Rates

A switching experiment with KPO 1 (while KPO 2 is detuned) is shown in Fig.
1(a). There, we observe that the system resides in each phase state for a
certain dwell time before switching to the opposite state. The average dwell
time $\tau$ can be expressed as a rate of activated switching
$\Gamma=\tau^{-1}$. In the rotating phase space, the same measurement data can
be represented as a density of count rates, see Fig. 1(b). As discussed in
Ref. Dykman  _et al._ (1998), the logarithm of the switching rate is inversely
proportional to the distance between the two phase states. In Fig. 1(c), this
leads to an exponential decrease in $\Gamma$ with decreasing $f_{d}$ (the
direction in $f_{f}$ in which $\Gamma$ increases depends on the sign of the
nonlinearity).

Figure 2: Noise-induced switching in the two-KPO system. (a) Quasi-static
frequency sweep with a parametric drive $U_{d}=$3.7\text{\,}\mathrm{V}$$
applied to both resonators. The full data (including the $u_{1,2}$ and the
stability diagram) are shown in Fig. 4. Arrows indicate the sweep direction.
Vertical dashed lines indicate the bifurcation frequencies, c.f. (e). (b)
Switching between the two symmetric phase states at
$f_{d}=$2.37\text{\,}\mathrm{MHz}$$ in the presence of applied voltage noise
(see Appendix B for details). (c) Switching rate $\Gamma$ as a function of
$f_{d}$ for a noise PSD
$S_{n}=$1.1\text{\times}{10}^{-7}\text{\,}{\mathrm{V}}^{2}\text{\,}{\mathrm{Hz}}^{-1}$$.
Dark (light) blue corresponds to experimental (simulated) data with error bars
estimated from textbook Poisson statistics. A thick red line reflects the
exponential trend. (d) Switching between the two symmetric states visualized
in the rotating frame in terms of $v_{S}$ and $v_{A}$. I:
$2.37\text{\,}\mathrm{MHz}$, II: $2.3675\text{\,}\mathrm{MHz}$, III:
$2.3625\text{\,}\mathrm{MHz}$, IV: $2.36\text{\,}\mathrm{MHz}$. Greyscale heat
maps represent the normalized counts $c$ from multiple switches. Circles and
squares indicate stable and unstable solutions, respectively. The ghost states
are marked in yellow. The solid lines are the theoretically predicted
switching paths between the stable symmetric states using the Onsager-Machlup
formalism, cf. Eq. (3). The color indicates whether a path passes the unstable
ghost state (red) or the unstable 0-amplitude state (blue-gray). (e)
Schematics of the steady states and bifurcation points as a function of
$f_{d}$. Solid and dashed lines are for stable and unstable solutions,
respectively. Squares mark ghost bifurcations that proliferate additional
unstable states (ghost states) from an already unstable attractor. Ghost
states are marked in yellow. (f) Activation $W$ as a function of $f_{d}$. Dark
(light) blue corresponds to experimental (simulated) data with error bars
estimated from textbook Poisson statistics. The solid gray and red lines are
obtained by minimizing $S_{\rm OM}$ for paths via the 0-amplitude state and
the antisymmetric state, respectively. The system parameters are: $Q=265$,
$f_{0}=$2.3670\text{\,}\mathrm{MHz}$$,
$\alpha=$-6.5\text{\times}{10}^{17}\text{\,}{\mathrm{V}}^{-2}\text{\,}{\mathrm{s}}^{-2}$$,
$U_{th}=$1.73\text{\,}\mathrm{V}$$, and
$J=$-1.28\text{\,}\mathrm{M}\mathrm{Hz}^{2}$$.

We now consider the case that the two KPO are tuned to have the same
frequency. The experimental phase diagram of the system as a function of
$\lambda$ and $\omega_{p}$ was characterized in Ref. Heugel  _et al._ (2022a).
In Fig. 2(a), we show the results of a measured frequency sweep passing
through various regions containing up to four different types of stable two-
oscillator states: a state with both resonators having amplitude zero; (S)
symmetric states, i.e., the two oscillators have the same phase; (A)
antisymmetric states where the two oscillators have $\pi$-shifted phases; and
(M) mixed-symmetry states that are neither symmetric nor antisymmetric. The
symmetric/antisymmetric solutions can be interpreted as the parametrically
driven symmetric/antisymmetric normal modes of the two resonators. See Fig. 4
for more information regarding the stability diagram.

Typically, when initialized in one such state, the resonator explores the
vicinity of the attractor under the influence of the noise terms $\Xi_{u_{i}}$
and $\Xi_{v_{i}}$ Heugel  _et al._ (2022b). Occasionally, the noise activates
the system over the quasi-potential barrier to reach a different attractor, as
seen for instance in Fig. 2(b) Margiani  _et al._ (2022). The system can
switch back and forth in time and is characterized by a switching rate
$\Gamma$. In the following, we analyze $\Gamma$ along the sweep in Fig. 2(a)
for $f_{d}>$2.36\text{\,}\mathrm{MHz}$$. In this frequency range, only the
symmetric solution is stable Heugel  _et al._ (2022a). We therefore expect
$\Gamma(f_{d})$ to manifest an exponential behaviour analogous to that seen in
a single KPO in Fig. 1. Motivated by this idea, we set out to confirm the
hypothesis that each normal mode has the same activation rate scaling as a
single KPO.

In Fig. 2(c), we show the measured transition rate $\Gamma$ for switches
between the two symmetric states as a function of $f_{d}$ while the parametric
driving strength $U_{d}$ is fixed. With decreasing $f_{d}$, the transition
rate $\Gamma$ decreases exponentially as expected from previous single KPO
studies Dykman  _et al._ (1998); Ryvkine and Dykman (2006); Stambaugh and Chan
(2006); Chan  _et al._ (2008); Margiani  _et al._ (2022). This holds for a
distinct range marked as I and II until
$f_{d}\approx$2.365\text{\,}\mathrm{MHz}$$ is reached. Surprisingly, below
this frequency, we observe substantial deviations from the simple exponential
model even though the symmetric state remains the only stable configuration.
First, a kink manifests at $f_{d}=$2.365\text{\,}\mathrm{MHz}$$, implying a
relative enhancement of the switching in region III. At even lower
frequencies, the slope of $\Gamma$ changes sign, signaling a significant
increase of the switching with decreasing frequency in region IV. This
behaviour is fundamentally at odds with the standard expectation of decreasing
rates.

To obtain deeper insights into the curious switching behaviour of the two-KPO
system, we look at several measured transition events and systematically
collect the 4-quadrature state vector $\mathbf{Y}=(u_{1},v_{1},u_{2},v_{2})$,
cf. Appendix D. To visualize the transitions in a 2-dimensional space, we use
symmetric and antisymmetric coordinates, $v_{\rm S}=(v_{1}+v_{2})/\sqrt{2}$,
$v_{\rm A}=(v_{1}-v_{2})/\sqrt{2}$ (and analogous for $u$). By plotting
several time traces containing multiple transition events in the phase space
spanned by $v_{S}$ and $v_{A}$, we obtain the corresponding probability
distribution.

Comparing the distributions at four representative frequencies, we find
striking differences, see Fig. 2(d). Regions I and II are characterized by two
attractors with high probabilities, corresponding to the two phase states of
the symmetric mode. By contrast, for $f_{d}<$2.365\text{\,}\mathrm{MHz}$$ in
the regions III and IV, we find the emergence of a substantial probability
centered around the attractors of the antisymmetric modes. This indicates that
the appearance of the kink at $f_{d}=$2.365\text{\,}\mathrm{MHz}$$ follows a
drastic change of the transition path chosen by the system, and of the
underlying quasi-potential landscape. Such a deviation in region IV is to be
expected as the activation dynamics now involves four different stable
attractors. In other words, the naive model of activation between two states
is insufficient in region IV. This is confirmed by a numerical simulation of
the noisy time evolution of the EOMs, given in Eq. (II), which is in accord
with the experimental results, see bright blue data in Fig. 2(c). Crucially,
however, we note that the antisymmetric state is not stable in region III, as
indicated by a red square (instead of a circle). This observation raises the
question: why should the unstable solutions of the syste influence the
transition paths? In the following, we precisely address this question through
an in-depth theoretical study of the transition dynamics.

##  IV Transition Paths

The weak noise-induced switching between stable oscillation states is
analogous to noise-activated jumping over a barrier $W$ studied in an
equilibrium system Stambaugh and Chan (2006); Hänggi  _et al._ (1990). The
principal difference is that in our driven system, the barrier $W$ between two
stable attractors resides in a quasipotential structure in a rotating frame.
We use the Onsager-Machlup formalism to identify the optimal transition paths
in phase space between two stable attractors, whose corresponding action then
provides an estimation for the barrier $W$. We first define the Onsager-
Machlup function Lehmann  _et al._ (2003); Wio (2013)

$S_{\rm
OM}[\mathbf{Y}]=\int_{t_{i}}^{t_{f}}\frac{1}{4}\left(\dot{\mathbf{Y}}-\mathbf{f}(\mathbf{Y})\right)^{2}dt\,,$
(3)

where $t_{i}$ ($t_{f}$) is the initial (final) time of the trajectory of a
system composed of $N$ resonators,
$\mathbf{Y}=(u_{1},v_{1},...,u_{N},v_{N})^{T}$, and $\mathbf{f}(\mathbf{Y})$
is the right hand side of Eq. (II) without noise terms written as a column
vector. The switching probability density between two stable states connected
by a path $\mathbf{Y}(t)$ is given by $e^{-2S_{\rm
OM}[\mathbf{Y}]/\sigma^{2}}$. The total switching probability $P_{if}$ from a
stable attractor at $\mathbf{Y}_{i}$ to one at $\mathbf{Y}_{f}$ is obtained by
integration over all allowed trajectories connecting them. From this
probability one can derive the switching rate $\Gamma$ Lehmann  _et al._
(2003), which in the weak-noise limit scales as
$\Gamma\propto\exp(-2W/\sigma^{2})$ with barrier $W$ Stambaugh and Chan
(2006).

At low noise, the switching rate $\Gamma$ is dominated by the path
$\mathbf{Y}_{\rm min}$ that minimizes $S_{\rm OM}$ Lehmann  _et al._ (2003).
Hence, we can neglect the integration over all possible paths and the
switching rate is approximately given by

$\Gamma\approx\Gamma_{\rm min}\equiv\Gamma_{0}e^{-2S_{\rm OM}[\mathbf{Y}_{\rm
min}]/\sigma^{2}}\,,$ (4)

where $\Gamma_{0}$ is an overall prefactor and we identify the effective
activation barrier $W_{\rm eff}=S_{\rm OM}[\mathbf{Y}_{\rm min}]$ Lehmann  _et
al._ (2003). If $S_{\rm OM}$ has multiple local minima, one needs to find the
contribution from all relevant minimizing paths and weigh their relative
contributions.

The landscape of $S_{OM}$ is characterized by a bifurcation diagram of the
system equations, which categorize both stable and unstable solutions as
effective minima and saddles/maxima, respectively. In Fig. 2(d), we
schematically show this bifurcation diagram along the frequency sweep in Fig.
2(a). Various pitchfork bifurcations lead to the emergence of stable and
unstable states. Interestingly, some of the bifurcations involve purely
unstable states and are not visible in a frequency sweep. We therefore dub
them “ghost bifurcations” as they leave the deterministic steady state physics
unaffected. Such ghost bifurcations separate regions I and II, as well as II
and III. At each ghost bifurcation, additional unstable states emerge. To
emphasize their origin, we refer to these unstable states as “ghost states”.
From comparing Fig. 2(c) and (e), it becomes clear that the emergence of the
ghost states is accompanied by new switching paths, and that these changes
impact $\Gamma$. Specifically, rather than acting as additional obstacles in
the switching paths, the ghost states appear to favor an increase in $\Gamma$.

To elucidate the surprising role of the ghost states, we study the transition
paths $\mathbf{Y}_{\rm min}$ and the corresponding barriers $W$ for the
different representative frequencies along the sweep shown in Fig. 2(a). To
this end, we need to minimize Eq. (3) which is a complex task. A simple
variational scheme with equal timesteps gives inconsistent results and more
advanced methods as the sgMAM-method Grafke  _et al._ (2017) are necessary to
obtain the correct physical paths (see details in Appendix C).

In Fig. 2(d), we show the calculated switching paths corresponding to the
experimental parameters, which are representative for the regions (I-IV). In
region I where the symmetric states are the only stable attractors, we find
only one switching path, which passes through the intermediate 0-amplitude
state in agreement with our experimentally observed distributions. As
expected, this is equivalent to the single KPO case in Fig. 1, confirming that
we should expect a monotonic decrease transition rate $\Gamma(f_{d})$ with
decreasing frequency as the stable attractors move apart Dykman  _et al._
(1998); Chan  _et al._ (2008).

In region II, we find two additional switching paths that avoid the
0-amplitude state and instead pass through emergent unstable antisymmetric
states. This is in line with the experimental data that exhibits a broader
distribution around the 0-amplitude state, extending to the two unstable
antisymmetric states. Similar paths also arise in regions III and IV, where
the unstable states provide transient ledges where the system can hover during
switching events, cf. Appendix D. The additional switching paths are visible
in the experiment and seem to be the dominant paths in regions III and IV.
These alternative switching paths bring forth a complexity to the system
dynamics, a feature intrinsically related to the existence of multiple normal
modes in KPO networks. In our system, we observe that the minimal path,
$\mathbf{Y}_{\rm min}$, connects the stable states always via an unstable one
in see Fig. 2(d) Maier and Stein (1992); Tang  _et al._ (2017), although
exceptions have been observed Luchinsky  _et al._ (1999b); Feng  _et al._
(2014). This underscores the importance of investigating stable as well as
unstable states of the system in order to understand the stochastic dynamics
of Ising networks.

Based on our theoretical analysis of optimal switching paths, we can now
obtain $W$ and compare it with that extracted from the experimental data, cf.
Appendix C. From this, we identify the dominant transition path, see Fig.
2(d). Within regions I-III, we find good qualitative and approximate
quantitative agreement (save for a small overall shift). Interestingly, in
region II, both the paths contribute equally to the activation: one via the
antisymmetric ghost states and the other via the origin. However, this effect
is not sufficiently strong to be observed in Fig. 2(c). In region III, the
former overtakes the latter, which manifests as the observed kink in Fig.
2(c). This indicates that the ghost states support the antisymmetric
switching, and markedly participate and modify the expected stochastic
dynamics of the system. In region IV, both symmetric and antisymmetric phase
states are stable. Here, the analytical calculation deviates from the
experimental and numerical results. This deviation is likely due to the fact
that the Onsager-Machlup method only considers switches that connect the two
symmetrical states, while the counting algorithm that we used for the
experimental and numerical data includes all possible switches between states.
The latter includes repeated switches between a symmetric and an antisymmetric
state as individual events.

##  V Discussion

At a fundamental level, our results demonstrate unambiguously the existence of
an unusual type of bifurcation arising from the coupling between the
individual resonators. Although these ghost bifurcations remain undetected in
a deterministic system characterization Heugel  _et al._ (2022a), they impact
the inter-state switching path, switching rate $\Gamma$, and the relative
dwell times in the symmetric and antisymmetric states. This inevitably affects
the stochastic switching processes in KPO networks, and their characterization
via stochastic sampling Margiani  _et al._ (2023).

Generalizing our results to $N$ identical resonators with identical all-to-all
coupling, we theoretically predict that the ghost bifurcations, as well as the
altered stability of the symmetric/antisymmetric regions, persist. Such a
network has a single symmetric lobe and a $N-1$-fold degenerate antisymmetric
lobe Heugel  _et al._ (2019). The boundaries delineating the overlapping
region once more involve ghost bifurcations. However, one of the ghost
bifurcations is now $N-1$ fold degenerate. As demonstrated above for two
resonators, the antisymmetric states only become stable when they undergo
another bifurcation, described by the general equation

$\lambda_{A}=\frac{\sqrt{4\gamma^{2}\omega^{2}+\left(J(N+2)-2\left(\omega^{2}-\omega_{0}^{2}\right)\right)^{2}}}{\omega_{0}^{2}}\,,$
(5)

resulting in an extended symmetric region. We thus expect a KPO network with
identical all-to-all coupling to manifest the four different regimes of
activation seen in the two-KPO case. When the coupling coefficient $J_{ij}$ is
different for each resonator pair, the degeneracy is lifted. However, the
instability lobes still overlap and generally form ghost bifurcations which
will also impact the switching behaviour of the system. In analogy to our
observations in Fig. 2, many competing switching paths open up. Finding the
dominant switching paths is a complex and demanding task. Our work provides
the motivation for such future investigations.

All of our observations have important consequences for logic networks built
from KPOs and nonlinear resonators in general, because they impact the
solution that a network will choose after a finite transition time. With
proper modelling and calibration of the ghost bifurcations, a network can be
operated at a position in $f_{d}-\lambda$ space where the many-body character
of the network is preserved and the annealing speed is optimized. It becomes
clear from our work that large networks bear very complex switching dynamics
and a careful analysis of the bifurcation topology is very important. Future
work might find ways to use this complexity in an advantageous manner to
perform faster calculations.

##  VI Summary and Outlook

In summary, we experimentally and theoretically investigated the noise-induced
dynamics of a system of two coupled nonlinear Kerr parametric oscillators
(KPOs). Our study implements the smallest form of a KPO network, and tests the
switching behavior in so-called Ising machines. We found that ghost
bifurcations play an important role, with consequences for the switching
dynamics of the system as it progresses towards its most stable configuration.
A better understanding of such effects can be very helpful for the calibration
to and for stochastic logic protocols, such as simulated annealing. As coupled
networks of parametric resonators are one of the main candidates for future
parallel computation architectures, our study provides crucial input for a
growing subcommunity working towards classical and quantum analog computation
Gottesman  _et al._ (2001); Devoret and Schoelkopf (2013); Mahboob  _et al._
(2016); Inagaki  _et al._ (2016); Goto (2016); Puri  _et al._ (2019); Nigg
_et al._ (2017); Dykman  _et al._ (2018). Furthermore, it provides additional
incentives for the fundamental exploration of complex driven-dissipative
nonlinear networks in a multitude of fields Dykman (2012).

###### Acknowledgements.

This work received financial support from the Swiss National Science
Foundation through grants (CRSII5_206008/1) and (PP00P2_163818), and the
Deutsche Forschungsgemeinschaft (DFG) through project number 449653034 and
SFB1432. We thank Peter Märki, Žiga Nosan and Christian Marty for technical
help.

##  Appendix A Single KPO

For $J=0$, each resonator can be driven into parametric resonance when
$U_{d}\geq U_{th}$ Landau and Lifshitz (1976); Lifshitz (2009). We
characterize each resonator using frequency sweeps as described in Leuch  _et
al._ (2016); Heugel  _et al._ (2019); Nosan  _et al._ (2019) and obtain the
values $Q_{1}=295$, $f_{0}=$2.6784\text{\,}\mathrm{MHz}$$,
$\alpha_{1}=$-9\text{\times}{10}^{17}\text{\,}{\mathrm{V}}^{-2}\text{\,}{\mathrm{s}}^{-2}$$,
and $U_{th}=$1.21\text{\,}\mathrm{V}$$. Using Eqs. (II), we can describe the
steady-state of the single KPO by applying the condition
($\dot{u}_{1}=\dot{v}_{1}=0$) Papariello  _et al._ (2016); Leuch  _et al._
(2016); Eichler and Zilberberg (2023); Košata  _et al._ (2022). This yields a
quintic characteristic polynomial with up to three different stable solutions
(attractors) in phase space, cf. Fig. 3(b). As a function of $f_{d}$, the
number of stable solutions changes at specific bifurcation points. In the
single KPO, we only observe pitchfork bifurcations, which involve at least one
stable solution, cf. Fig. 3(c). In Fig. 3(d) we show the characteristic
parametric instability lobe. Region (i) accommodates only one stable state
with amplitude 0. Inside the region marked as (ii), the linear resonator
becomes unstable, bifurcates, and settles into one of the two steady states
that are stabilized by $\alpha$ Lifshitz (2009). These phase states have the
same amplitude but are $\pi$-shifted in phase, cf. Fig. 3(b). In region (iii),
the phase states coexist with the amplitude 0 solution.

Single KPO –  In Fig. 1 we inspect the noise-induced switching of a single
KPO, whose properties are well known Dykman  _et al._ (1998); Ryvkine and
Dykman (2006); Chan  _et al._ (2008); Dykman  _et al._ (2018). As expected, we
find that $\Gamma$ decreases monotonically with increasing separation between
the phase states, which we control here through $f_{d}$ Dykman  _et al._
(1998). Similar results have been previously measured in other KPO
implementations Chan  _et al._ (2008).

The monotonic decrease of $\Gamma$ in Fig. 1(c) is derived using the Onsager-
Machlup approach Dykman  _et al._ (1998). Specifically, at low noise, the
switching rate $\Gamma$ is dominated by the path $\mathbf{Y}_{\rm min}$ that
minimizes $S_{\rm OM}$ Lehmann  _et al._ (2003). Repeating this estimation as
a function of $f_{d}$ and calculating $\Gamma_{\rm min}$ yields good agreement
with the experimentally observed $\Gamma$, cf. Fig. 1(c). Note that the
prefactor $\Gamma_{0}$ is not obtained by this method but reused from Ref.
Dykman  _et al._ (1998), leading to a slight overall shift towards larger
$\Gamma$. The analytical formula derived in Ref. Dykman  _et al._ (1998)
produces a similarly good agreement, cf. Eq. (6) in Appendix C.

Figure 3: (a) Schematics of a parametric RLC resonator with resistance $R$,
inductance $L$, nonlinear capacitance $C$, and tuning voltage $U_{1}$. (b)
Schematic representation of the parametric phase states (wine red, orange),
and the 0-amplitude state (grey) in phase space. (c) Schematics of the steady
states and bifurcation points as a function of $f_{d}$. Solid (dashed) lines
correspond to stable (unstable) solutions. (d) Stability phase diagram: (i)
White: only the 0-amplitude solution is stable; (ii) orange: only the phase
states are stable; (iii) light orange: 0-amplitude and phase states are both
stable.

##  Appendix B Fluctuating versus Coherent Signal Amplitude

To obtain an optimal agreement between the measured switching rates and the
theoretical predictions, we consistently found that the noise power spectral
density in the model had to be a factor $\approx 4.2$ smaller than the value
applied in the experiment. This discrepancy is likely due to an additional
attenuation of a factor 2 in the path of the fluctuating voltage, for instance
a voltage division at a $50\text{\,}\mathrm{\SIUnitSymbolOhm}$ matched input
port. The fluctuating signal with power spectral density $S_{n}$ was provided
by two dedicated voltage sources with the same output intensity and added to
the coherent signal via the ADD channel of the Zurich Instruments HF2LI lock-
in amplifier. The resulting noise process $\xi_{i}$ acting on our system has a
power spectral density $\varsigma^{2}=C_{in}S_{n}$, where the coefficient for
the signal in-coupling efficiency is
$C_{in}=$4.93\text{\times}{10}^{-20}\text{\,}\mathrm{Hz}^{4}$$ for the single-
KPO experiment. For the two-KPO experiment, we find best agreement for a
slightly lower value for $C_{in}$, which is probably due to differences in the
coil geometry between the devices or between the experimental runs.

Figure 4: (a) stability diagram of the system of two coupled KPOs in a space
spanned by $f_{d}$ and $U_{p}$. The red and blue lines mark the instability
lobes of the parametrically driven antisymmetric and symmetric normal modes,
respectively. The frequency sweep in Fig. 2(a) and in (b) of this figure was
taken at the $U_{p}$ value indicated by a dashed horizontal line, with dots
and squares corresponding to the bifurcations shown in Fig. 2(e). (b) Full
data of the frequency sweep shown in Fig. 2(a). Arrows indicate the sweep
direction. Figure 5: Switching rate $\Gamma$ between the symmetric states of
the two-KPO system as a function of inverse noise strength, $1/\sigma^{2}$,
obtained from simulation (black) for $U_{d}=$3.7\text{\,}\mathrm{V}$$ and
$f_{d}=$2.3725\text{\,}\mathrm{MHz}$$. The optimal fit (gray) corresponds to
Eq. (4) for the fitting parameters
$\Gamma_{0}=$2\text{\times}{10}^{4}\text{\,}\mathrm{Hz}$$ and
$W=$0.024\text{\,}\mathrm{Hz}\mathrm{V}^{2}$$.

##  Appendix C Details on the calculations for noise-induced switching

Determination of the switching rate –  The experimental determination of the
switching rate $\Gamma$ in Fig. 2(c) was performed with a lock-in amplifier
(Zurich Instruments HF2LI). We used a sampling rate of
$450\text{\,}\mathrm{Hz}$ and a total measurement time of
$300\text{\,}\mathrm{s}$ for $f_{d}\leq$2.3696\text{\,}\mathrm{Hz}$$ and
$60\text{\,}\mathrm{s}$ for $f_{d}>$2.3696\text{\,}\mathrm{Hz}$$. Counting of
the switching events was done with a numerical algorithm that compared the
amplitudes and phases of successive measurement points for an entire time
trace measurement. Concretely, the program increased the switching counter by
1 if the phase difference of two successive points was above a ‘phase
threshold’ ($130\text{\,}\mathrm{\SIUnitSymbolDegree}$) while at least one of
the points was above an ‘amplitude threshold’ ($0.5\text{\,}\mathrm{mV}$), or
if exactly one out of two successive points was above the amplitude threshold.
The same algorithm was used to evaluate the switching rate in numerical
simulations that emulated the measurements (including the effective sampling
rate). Similar results were obtained by finding the maximal turning point of
Allan deviations of the phase Margiani  _et al._ (2022).

Figure 6: Examples of switching events in the two-KPO system. (a) For
$U_{d}=$3.7\text{\,}\mathrm{V}$$ and $f_{d}=$2.37\text{\,}\mathrm{MHz}$$
(marked as I in Fig. 2(c)), all four measured coordinates
($u_{1},v_{1},u_{2},v_{2}$) switch simultaneously on our sampling timescale.
The system switches from one symmetric (S) configuration to the opposite one
via the 0-amplitude (0) state. (b) For $U_{d}=$3.7\text{\,}\mathrm{V}$$ and
$f_{d}=$2.36\text{\,}\mathrm{MHz}$$ (marked as IV in Fig. 2(c)), the
coordinates ($u_{1},v_{1}$) switch first, followed by ($u_{2},v_{2}$) after a
delay of roughly $0.5\text{\,}\mathrm{ms}$. In the time span between the
jumps, the system occupies the antisymmetric (A) state.

Determination of the activation barrier $W$ –  We simulate the coupled KPO
system at different noise strengths, see Fig 5, and we verify that our
analysis is in the low-noise limit by showing that Eq. (4) is obeyed. The
optimal fit (gray) yields
$\Gamma_{0}=$2\text{\times}{10}^{4}\text{\,}\mathrm{Hz}$$ and
$W=$0.024\text{\,}\mathrm{Hz}\mathrm{V}^{2}$$. In this limit, $\Gamma_{0}$ and
$W$ are independent of the noise and purely depend on the properties of the
non-stochastic system, and on the switching paths in phase space. This
procedure allows us to extract $\Gamma_{0}$ and $W$ for the experimental and
the numerical data in Fig. 2(d), as well as calculate $\Gamma$ at different
noise strengths. For the numerical switching rate in Fig. 2(c), we used this
scaling law to convert the numerical data, simulated at $1.3$ times stronger
noise, to the experimentally applied noise.

Analytic expression for a single KPO –  For a single parametric Kerr
oscillator, the switching rate was calculated in Dykman  _et al._ (1998) and
is given by

$\Gamma=\frac{\left(\gamma\sqrt{\frac{\lambda^{2}\omega_{0}^{4}}{\gamma^{2}\omega^{2}}-4}-4\omega+4\omega_{0}\right)\sqrt{\left|1-\frac{\lambda^{2}\omega_{0}^{4}}{4\gamma^{2}\omega^{2}}\right|}\exp\left(-\frac{\gamma^{2}\omega^{3}\left(\gamma\sqrt{\frac{\lambda^{2}\omega_{0}^{4}}{\gamma^{2}\omega^{2}}-4}-4\omega+4\omega_{0}\right)^{2}\sqrt{\left|1-\frac{\lambda^{2}\omega_{0}^{4}}{4\gamma^{2}\omega^{2}}\right|}}{3\alpha\lambda^{2}\sigma^{2}\omega_{0}^{4}}\right)}{2\sqrt{2}\pi}\,.$
(6)

Details on the path optimization –  The switching rate $\Gamma$ can be
described by $S_{\rm OM}[\mathbf{Y}_{\rm min}]$. Obtaining the minimal path
$\mathbf{Y}_{\rm min}$ is a complex task for coupled parametric resonators. As
a simple variation of a discretized path with equal timesteps fails to obtain
correct results, we use the sgMAM method Grafke  _et al._ (2017). It is an
improved path optimization scheme based on scaled time, leading to consistent
converged results. We start with a guessed initial path that connects two
stable states via an unstable attractor. Then, we perform numerical
minimization of $S_{\rm OM}$ by varying the path in phase space between the
chosen end points. For the single KPO, we choose the two phase states as
initial and final state, see Fig 1(b). In the coupled system with
$\mathbf{Y}=(u_{1},v_{1},u_{2},v_{2})^{T}$, we choose one of the symmetric
states as the initial point and the other one as the final point, and try
different unstable states as intermediate points. We thus obtain the
corresponding locally-minimizing switching paths $\mathbf{Y}_{\rm min}$, see
Fig. 2(d).

##  Appendix D Switching via the antisymmetric state

In Fig. 6, we show examples of timetraces during noise-induced switching
between symmetric states. Figure 6(a) corresponds to
$f_{d}=$2.37\text{\,}\mathrm{MHz}$$ in region I of Fig. 2(c), where switches
occur via the unstable 0-amplitude state because both resonators switch
synchronously. In Fig. 6(b), we show an example for
$f_{d}=$2.36\text{\,}\mathrm{MHz}$$ in region IV, where the two resonators
switch with a finite delay. In the short time interval between the two
switches, the system dwells in the antisymmetric state.

## References

  * Kramers (1940) H. Kramers, Physica 7, 284 (1940).
  * Hänggi  _et al._ (1990) P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
  * Rondin  _et al._ (2017a) L. Rondin, J. Gieseler, F. Ricci, R. Quidant, C. Dellago, and L. Novotny, Nature nanotechnology 12, 1130 (2017a).
  * Best and Hummer (2006) R. B. Best and G. Hummer, Phys. Rev. Lett. 96, 228104 (2006).
  * Chung  _et al._ (2015) H. S. Chung, S. Piana-Agostinetti, D. E. Shaw, and W. A. Eaton, Science 349, 1504 (2015), https://www.science.org/doi/pdf/10.1126/science.aab1369 .
  * García-Müller  _et al._ (2008) P. L. García-Müller, F. Borondo, R. Hernandez, and R. M. Benito, Phys. Rev. Lett. 101, 178302 (2008).
  * Badzey and Mohanty (2005) R. L. Badzey and P. Mohanty, Nature 437, 995 (2005).
  * Rondin  _et al._ (2017b) L. Rondin, J. Gieseler, F. Ricci, R. Quidant, C. Dellago, and L. Novotny, Nature Nanotechnology 12, 1130 (2017b).
  * Fulton and Dunkleberger (1974) T. A. Fulton and L. N. Dunkleberger, Phys. Rev. B 9, 4760 (1974).
  * Silvestrini  _et al._ (1988) P. Silvestrini, S. Pagano, R. Cristiano, O. Liengme, and K. E. Gray, Phys. Rev. Lett. 60, 844 (1988).
  * Dykman  _et al._ (1998) M. I. Dykman, C. M. Maloney, V. N. Smelyanskiy, and M. Silverstein, Phys. Rev. E 57, 5202 (1998).
  * Dykman (2007) M. I. Dykman, Phys. Rev. E 75, 011101 (2007).
  * Tadokoro  _et al._ (2020) Y. Tadokoro, H. Tanaka, and M. Dykman, Scientific Reports 10, 10413 (2020).
  * Luchinsky  _et al._ (1999a) D. G. Luchinsky, R. S. Maier, R. Mannella, P. V. E. McClintock, and D. L. Stein, Phys. Rev. Lett. 82, 1806 (1999a).
  * Lapidus  _et al._ (1999) L. J. Lapidus, D. Enzer, and G. Gabrielse, Phys. Rev. Lett. 83, 899 (1999).
  * Kim  _et al._ (2005) 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).
  * Aldridge and Cleland (2005) J. S. Aldridge and A. N. Cleland, Phys. Rev. Lett. 94, 156403 (2005).
  * Chan and Stambaugh (2007) H. B. Chan and C. Stambaugh, Phys. Rev. Lett. 99, 060601 (2007).
  * Chan  _et al._ (2008) H. B. Chan, M. I. Dykman, and C. Stambaugh, Phys. Rev. Lett. 100, 130602 (2008).
  * Venstra  _et al._ (2013) W. J. Venstra, H. J. Westra, and H. S. Van Der Zant, Nature communications 4, 1 (2013).
  * Mahboob  _et al._ (2014) I. Mahboob, M. Mounaix, K. Nishiguchi, A. Fujiwara, and H. Yamaguchi, Scientific Reports 4, 4448 (2014).
  * Margiani  _et al._ (2022) G. Margiani, S. Guerrero, T. L. Heugel, C. Marty, R. Pachlatko, T. Gisler, G. D. Vukasin, H.-K. Kwon, J. M. Miller, N. E. Bousse, _et al._ , Applied Physics Letters 121, 164101 (2022).
  * Frimmer  _et al._ (2019) M. Frimmer, T. L. Heugel, i. c. v. 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).
  * Ryvkine and Dykman (2006) D. Ryvkine and M. I. Dykman, Physical Review E 74, 061118 (2006).
  * Mahboob and Yamaguchi (2008) I. Mahboob and H. Yamaguchi, Nature Nanotechnology 3, 275 (2008).
  * Wilson  _et al._ (2010) C. M. Wilson, T. Duty, M. Sandberg, F. Persson, V. Shumeiko, and P. Delsing, Phys. Rev. Lett. 105, 233907 (2010).
  * Eichler  _et al._ (2011) A. Eichler, J. Chaste, J. Moser, and A. Bachtold, Nano Letters 11, 2699 (2011), pMID: 21615135.
  * Leuch  _et al._ (2016) A. Leuch, L. Papariello, O. Zilberberg, C. L. Degen, R. Chitra, and A. Eichler, Phys. Rev. Lett. 117, 214101 (2016).
  * Gieseler  _et al._ (2012) J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny, Phys. Rev. Lett. 109, 103603 (2012).
  * Lin  _et al._ (2014) Z. Lin, K. Inomata, K. Koshino, W. D. Oliver, Y. Nakamura, J. S. Tsai, and T. Yamamoto, Nature Communications 5, 4480 (2014).
  * Puri  _et al._ (2017a) S. Puri, S. Boutin, and A. Blais, npj Quantum Information 3, 18 (2017a).
  * Eichler  _et al._ (2018) A. Eichler, T. L. Heugel, A. Leuch, C. L. Degen, R. Chitra, and O. Zilberberg, Applied Physics Letters 112, 233105 (2018), https://doi.org/10.1063/1.5031058 .
  * Nosan  _et al._ (2019) Z. Nosan, P. Märki, N. Hauff, C. Knaut, and A. Eichler, Phys. Rev. E 99, 062205 (2019).
  * Grimm  _et al._ (2019) 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 (2019).
  * Puri  _et al._ (2019) 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).
  * Miller  _et al._ (2019) J. M. Miller, D. D. Shin, H.-K. Kwon, S. W. Shaw, and T. W. Kenny, Phys. Rev. Applied 12, 044053 (2019).
  * Dykman (2012) M. Dykman,  _Fluctuating Nonlinear Oscillators_ (Oxford University Press, 2012).
  * Eichler and Zilberberg (2023) A. Eichler and O. Zilberberg,  _Classical and Quantum Parametric Phenomena_ (Oxford University Press, 2023).
  * Mahboob  _et al._ (2016) I. Mahboob, H. Okamoto, and H. Yamaguchi, Science Advances 2, e1600236 (2016).
  * Inagaki  _et al._ (2016) T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, Nature Photonics 10, 415 (2016).
  * Goto (2016) H. Goto, Scientific Reports 6, 21686 (2016).
  * Puri  _et al._ (2017b) S. Puri, C. K. Andersen, A. L. Grimsmo, and A. Blais, Nature Communications 8, 15785 (2017b).
  * Nigg  _et al._ (2017) S. E. Nigg, N. Lörch, and R. P. Tiwari, Science Advances 3 (2017), 10.1126/sciadv.1602273.
  * Dykman  _et al._ (2018) M. I. Dykman, C. Bruder, N. Lörch, and Y. Zhang, Phys. Rev. B 98, 195444 (2018).
  * Okawachi  _et al._ (2020) Y. Okawachi, M. Yu, J. K. Jang, X. Ji, Y. Zhao, B. Y. Kim, M. Lipson, and A. L. Gaeta, Nature Communications 11, 4119 (2020).
  * Calvanese Strinati  _et al._ (2021) M. Calvanese Strinati, L. Bello, E. G. Dalla Torre, and A. Pe’er, Phys. Rev. Lett. 126, 143901 (2021).
  * Bello  _et al._ (2019) L. Bello, M. Calvanese Strinati, E. G. Dalla Torre, and A. Pe’er, Phys. Rev. Lett. 123, 083901 (2019).
  * Heugel  _et al._ (2022a) T. L. Heugel, O. Zilberberg, C. Marty, R. Chitra, and A. Eichler, Physical Review Research 4, 013149 (2022a).
  * Heugel  _et al._ (2022b) T. L. Heugel, A. Eichler, R. Chitra, and O. Zilberberg, arXiv preprint arXiv:2203.05577 (2022b).
  * Margiani  _et al._ (2023) G. Margiani, J. del Pino, T. L. Heugel, N. E. Bousse, S. Guerrero, T. W. Kenny, O. Zilberberg, D. Sabonis, and A. Eichler, Phys. Rev. Res. 5, L012029 (2023).
  * Heugel  _et al._ (2019) T. L. Heugel, M. Oscity, A. Eichler, O. Zilberberg, and R. Chitra, Phys. Rev. Lett. 123, 124301 (2019).
  * Lehmann  _et al._ (2003) J. Lehmann, P. Reimann, and P. Hänggi, physica status solidi (b) 237, 53 (2003), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.200301774 .
  * Wio (2013) H. S. Wio,  _Path Integrals for Stochastic Processes_ (WORLD SCIENTIFIC, 2013) https://www.worldscientific.com/doi/pdf/10.1142/8695 .
  * McLachlan (1951) N. McLachlan,  _Theory and application of Mathieu functions_ (Clarendon, 1951).
  * Lifshitz (2009) M. C. Lifshitz, R. Cross, “Nonlinear dynamics of nanomechanical and micromechanical resonators,” in  _Reviews of Nonlinear Dynamics and Complexity_ (Wiley-VCH, 2009) pp. 1–52.
  * Guckenheimer and Holmes (1990) J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ , Applied mathematical sciences (Springer-Verlag, 1990).
  * Papariello  _et al._ (2016) L. Papariello, O. Zilberberg, A. Eichler, and R. Chitra, Phys. Rev. E 94, 022201 (2016).
  * Khas’minskii (1966) R. Z. Khas’minskii, Theory of Probability & Its Applications 11, 390 (1966), https://doi.org/10.1137/1111038 .
  * Roberts and Spanos (1986) J. Roberts and P. Spanos, International Journal of Non-Linear Mechanics 21, 111 (1986).
  * Nayfeh and Mook (2008) A. H. Nayfeh and D. T. Mook,  _Nonlinear Oscillations_ , Physics textbook (Wiley, 2008).
  * Stambaugh and Chan (2006) C. Stambaugh and H. B. Chan, Phys. Rev. B 73, 172302 (2006).
  * Grafke  _et al._ (2017) T. Grafke, T. Schäfer, and E. Vanden-Eijnden, “Long term effects of small random perturbations on dynamical systems: Theoretical and computational tools,” in  _Recent Progress and Modern Challenges in Applied Mathematics, Modeling and Computational Science_, edited by R. Melnik, R. Makarov, and J. Belair (Springer New York, New York, NY, 2017) pp. 17–55.
  * Maier and Stein (1992) R. S. Maier and D. L. Stein, Phys. Rev. Lett. 69, 3691 (1992).
  * Tang  _et al._ (2017) Y. Tang, R. Yuan, G. Wang, X. Zhu, and P. Ao, Scientific Reports 7, 15762 (2017).
  * Luchinsky  _et al._ (1999b) D. G. Luchinsky, R. S. Maier, R. Mannella, P. V. E. McClintock, and D. L. Stein, Phys. Rev. Lett. 82, 1806 (1999b).
  * Feng  _et al._ (2014) H. Feng, K. Zhang, and J. Wang, Chem. Sci. 5, 3761 (2014).
  * Gottesman  _et al._ (2001) D. Gottesman, A. Kitaev, and J. Preskill, Phys. Rev. A 64, 012310 (2001).
  * Devoret and Schoelkopf (2013) M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013).
  * Landau and Lifshitz (1976) L. Landau and E. Lifshitz,  _Mechanics_ , Butterworth-Heinemann (1976).
  * Košata  _et al._ (2022) J. Košata, J. del Pino, T. L. Heugel, and O. Zilberberg, SciPost Phys. Codebases , 6 (2022).