# Hamiltonian reconstruction via ringdown dynamics

Vincent Dumont<EMAIL_ADDRESS>Laboratory for Solid State Physics, ETH
Zürich, CH-8093 Zürich, Switzerland  Markus Bestler  Department of Physics,
University of Konstanz, D-78457 Konstanz, Germany  Letizia Catalini
Laboratory for Solid State Physics, ETH Zürich, CH-8093 Zürich, Switzerland
Gabriel Margiani  Laboratory for Solid State Physics, ETH Zürich, CH-8093
Zürich, Switzerland  Oded Zilberberg  Department of Physics, University of
Konstanz, D-78457 Konstanz, Germany  Alexander Eichler  Laboratory for Solid
State Physics, ETH Zürich, CH-8093 Zürich, Switzerland

###### Abstract

Many experimental techniques aim at determining the Hamiltonian of a given
system. The Hamiltonian describes the system’s evolution in the absence of
dissipation, and is often central to control or interpret an experiment. Here,
we theoretically propose and experimentally demonstrate a method for
Hamiltonian reconstruction from measurements over a large area of phase space,
overcoming the main limitation of previous techniques. A crucial ingredient
for our method is the presence of dissipation, which enables sampling of the
Hamiltonian through ringdown-type measurements. We apply the method to a
driven-dissipative system – a parametric oscillator – observed in a rotating
frame, and reconstruct the (quasi-)Hamiltonian of the system. Furthermore, we
demonstrate that our method provides direct experimental access to the so-
called symplectic norm of the stationary states of the system, which is tied
to the particle- or hole-like nature of excitations of these states. In this
way, we establish a method to unveil qualitative differences between the
fluctuations around stabilized minima and maxima of the nonlinear out-of-
equilibrium stationary states. Our method constitutes a versatile approach to
characterize a wide class of driven-dissipative systems.

††preprint: APS/123-QED

##  I Introduction

The evolution of any physical system is governed by an interplay between
conservative and nonconservative forces. The former are generated by a
Hamiltonian ‘landscape’, i.e., the sum of all energy terms within a closed
system. Understanding and controlling an experiment usually require knowledge
of the system’s Hamiltonian. In a realistic setting, however, it can be
difficult to compute the Hamiltonian from first principles, as it requires
full insight into all microscopic constituents. Alternatively, extracting the
effective Hamiltonian from a measurement provides direct access to the system
dynamics even when its theoretical model is incomplete. Experimentalists
therefore invest much effort in calibrating the parameters of their system in
order to extract and fit their effective model [1]. A prime example is the
calibration of qubits, whose gate operations rely on precise Hamiltonian
estimations [2, 3, 4, 5].

A vast majority of systems are inherently ‘open’, i.e., the dissipative
coupling to an environment cannot be ignored. An open system experiences
fluctuations that cause it to sample the available landscape over time. This
allows estimating a system’s Hamiltonian by measuring the probability of
finding it in a certain state while it is subject to fluctuations [6, 7, 8, 9,
10, 11, 12, 13]. This method was used in experimental studies of the escape
dynamics [13] and Kramers turnover [11] of a particle trapped in an optical
potential, investigations of the force fields acting on these particles [14,
15], or the stability of coupled nonlinear systems [16, 17]. Probabilistic
methods have the drawback that certain states are rarely explored. As such,
the methods are often inadequate in situations where regions of interest (in
state space) are separated by large energy differences. In addition, these
methods are insensitive to temporal correlations in the system, which pertain
to the notion of causality of excitations, as quantified by e.g. out-of-time-
ordered correlators [18, 19].

In the presence of dissipative coupling to an environment, it is also possible
to probe a system’s Hamiltonian from its deterministic relaxation. This
becomes feasible when the measured variables of interest are large enough to
neglect the impact of fluctuations. By initializing a system in a well-defined
state using an external drive and then turning off the drive, the system
decays into a stationary state due to dissipative coupling to a large and,
most often, Markovian reservoir. By measuring such a ‘ringdown’ into a
stationary state (attractor), both the nonlinearity of a system close to a
stable solution [20, 21, 22, 23] and the curvature of engineered potential
landscapes [24] have been measured in the absence of a drive.

In principle, such ringdown-type experiments can be extended to driven-
dissipative systems. This is particularly interesting for driven nonlinear
systems with multiple stable oscillation states. In a rotating frame, such
oscillations can appear as stationary states within their respective basins of
attraction [25, 26]. There, the system dynamics can be described by a
rotating-frame quasi-Hamiltonian that, for instance, allows understanding out-
of-equilibrium phase transitions [27, 28, 29, 30, 31] in lattices of cold
atoms [24], optical oscillators [32], or cavity magnonic systems [33].
Interestingly, such systems can stabilize out-of-equilibrium phases where
seemingly anti-causal excitations can manifest [34, 29]. However, to our
knowledge, no full reconstruction of such a rotating-frame Hamiltonian has
been reported to date.

In this paper, we demonstrate a precise and deterministic method of
reconstructing the full rotating frame quasi-Hamiltonian of a driven-
dissipative nonlinear system using systematic ringdown measurements. Notably,
the presence of dissipation enables sampling of a large section of a
Hamiltionian from a limited number of ringdown measurements. The method
provides high resolution even far from stationary solutions, where stochastic
approaches typically fail. Furthermore, we obtain a direct measurement of the
‘symplectic norm’ of each solution [29, 30], providing a qualitative
understanding of the different phases the system can enter, including a
dissipation-stabilized maximum with hole-like quasiparticle excitations.
Importantly, our method can be extended and applied to both undriven and
driven-dissipative oscillating systems far from equilibrium, making it a
valuable tool in many contemporary fields of physics.

We first present in Sec. II our electro-mechanical resonator system. In Sec.
III, we introduce the general method to reconstruct the Hamiltonian from
ringdown measurements. In Sec. IV, we reconstruct the effective Hamiltonian of
our resonator in the absence of a drive – i.e., demonstrating the method for a
standard damped harmonic oscillator. Then, in Sec. V, we subject the resonator
to a large parametric drive to allow multiple stable oscillation solutions and
compare the reconstructed rotating-frame Hamiltonian to theoretical
predictions. Finally, in Sec. VI, we discuss the extraction and implications
of the symplectic norm as a tool for exploring and understanding out-of-
equilibrium stationary solutions.

##  II Electromechanical Device

Figure 1: Micro-electromechanical resonator and $H_{\mathrm{rot}}$
reconstruction of a simple damped harmonic oscillator. (a) A cantilever,
biased at voltage $U_{\mathrm{bias}}$, with oscillation displacement $x$ and
natural frequency $\omega_{0}$, is capacitively coupled to two conductors.
Typically, an input voltage $U_{\mathrm{in}}$ with frequency
$\omega_{\mathrm{r}}$ ($2\omega_{\mathrm{r}}$) can drive the mechanical
resonator resonantly (parametrically), when
$\omega_{\mathrm{r}}\approx\omega_{0}$. The mechanical motion induces a
voltage $U_{\mathrm{out}}$, which we read out. (b) Several ringdown
measurements for a few initialization points (dots) in rotating (i.e., at
demodulated frequency $\omega_{\mathrm{r}}$) phase space $(u,v)$ for blue
detuned demodulation, $(\omega_{\mathrm{r}}-\omega_{0})/2\pi\approx 91$ Hz.
(c) Several ringdown measurements in phase space for red detuned demodulation,
$(\omega_{\mathrm{r}}-\omega_{0})/2\pi\approx-209$ Hz. (d) Reconstructed
Hamiltonian for positive detuning, i.e., from measurements of panel (b), using
Eq. (5). (e) Reconstructed Hamiltonian for negative detuning, i.e. from the
measurements of panel (c).

Our device is a micro-electromechanical resonator, as illustrated in Fig.
1(a). It consists of a mechanical cantilever [35] biased with voltage
$U_{\mathrm{bias}}\approx 32$ V, used to tune the mechanical resonance
frequency $\omega_{0}$ [36], and which induces a Duffing non-linearity $\beta$
(due to the non-linearity of the electrostatic force between the mechanical
element and the electrodes [35]). To generate a weak near-resonant forcing
term $F$, we apply a voltage
$U_{\mathrm{in}}=U_{\mathrm{r}}\cos(\omega_{\mathrm{r}}t+\theta)$ to one of
the electrodes, where $U_{\mathrm{r}}$ is the amplitude, $\theta$ is a phase
offset, and $\omega_{\mathrm{r}}\approx\omega_{0}$. A large off-resonant
voltage $U_{\mathrm{in}}=U_{\mathrm{p}}\cos(2\omega_{\mathrm{r}}t+\psi)$
periodically modulates $\omega_{0}$ with a modulation depth $\lambda\propto
U_{\mathrm{p}}$, making it possible to parametrically drive the resonator
[37]. We read out the mechanical displacement $x$ via the output voltage
$U_{\mathrm{out}}\propto x$ with a lock-in amplifier (Zurich Instruments
MFLI), and forgo the proportionality coefficient for convenience (i.e., we
define $U_{\mathrm{out}}\equiv x$) [38]. The device displacement is described
by the equation of motion (EOM)

$\displaystyle\frac{d^{2}x}{dt^{2}}+\omega_{0}^{2}[1-\lambda\cos(2\omega_{\mathrm{r}}t+\psi)]x+\Gamma\frac{dx}{dt}+\beta
x^{3}=F,$ (1)

with time $t$. The device has a mechanical resonance frequency
$\omega_{0}/2\pi\approx 1.12$ MHz, energy decay rate $\Gamma/2\pi\approx 112$
Hz (quality factor $Q=\omega_{0}/\Gamma\approx 10^{4}$), and Duffing non-
linearity $\beta\approx-9.9\times 10^{16}$ V-2s-2. The forcing term (in units
of Vs-2) is $F=AU_{\mathrm{in}}$ with a conversion factor $A\approx 16\times
10^{6}$ s-2 and with the oscillator’s mass absorbed in the definition of $F$.
The parametric drive voltage amplitude $U_{\mathrm{p}}$ can be converted to a
parametric modulation depth $\lambda=2U_{\mathrm{p}}/QU_{\mathrm{th}}$ by
measuring the parametric drive threshold amplitude $U_{\mathrm{th}}\approx
1.98$ V. Above the parametric threshold, the system can be driven into an out-
of-equilibrium stationary oscillation via a spontaneous time-translation
symmetry breaking [39, 40, 8, 41, 42, 43, 44, 4, 45, 46, 37, 47]. The full
device characterization is shown in Appendix A.

##  III Hamiltonian Reconstruction

We now introduce our method for reconstructing the Hamiltonian from measured
coordinates along a trajectory.The evolution of a classical lossless system in
the laboratory frame with coordinate $x$ and its conjugate $p$ is given by
Hamilton’s EOMs

$\displaystyle\begin{aligned} \frac{dx}{dt}&=+\frac{\partial H(x,p)}{\partial
p},\\\ \frac{dp}{dt}&=-\frac{\partial H(x,p)}{\partial x},\end{aligned}$ (2)

where $H(x,p)$ is the Hamiltonian. Here, for any initial conditions and
without drive, trajectories form closed loops along equipotential lines in
coordinate space $(x,p)$ since there is no energy loss. This means that only a
single loop of the Hamiltonian can be sampled when following any given
trajectory. To sample all potential elevations, one has to initialize the
system at infinitely many different initial conditions.

In many cases, we are interested in the slowly-varying in-phase and out-of-
phase quadratures $u(t)$ and $v(t)$ observed in a frame rotating at
$\omega_{\mathrm{r}}\approx\omega_{0}$, which are defined via
$x(t)=u(t)\cos(\omega_{\mathrm{r}}t)-v(t)\sin(\omega_{\mathrm{r}}t)$. The
demodulated quadrature we measure with our lock-in amplifier at frequency
$\omega_{\mathrm{r}}$ are precisely $u$ and $v$. By applying the averaging
method [48] on an equation of motion such as Eq. (1), we obtain the so-called
slow-flow equations (see Appendix B). We open the system by adding dissipative
terms $\propto\Gamma$, i.e., we consider rotating Lagrangian dynamics, such
that the equations of motion can be formulated in a similar structure as Eqs.
(2):

$\displaystyle\begin{aligned} \frac{du}{dt}&=+\frac{\partial
H_{\mathrm{rot}}(u,v)}{\partial v}-\frac{\Gamma}{2}u,\\\
\frac{dv}{dt}&=-\frac{\partial H_{\mathrm{rot}}(u,v)}{\partial
u}-\frac{\Gamma}{2}v.\end{aligned}$ (3)

In this frame, the energy-conserving evolution of the system is governed by a
rotating-frame quasi-Hamiltonian $H_{\mathrm{rot}}$. For our resonator,
subjected to both parametric and external drives, $H_{\mathrm{rot}}$ reads
111Note that $H_{\mathrm{rot}}$ may not have dimensions of energy in this
notation because $u$ and $v$ have the same units. The quasi-Hamiltonian can be
expressed in units of energy by changing $H_{\mathrm{rot}}\rightarrow
H_{\mathrm{rot}}/m\omega_{\mathrm{r}}$ in Eqs. (3) and (III).

$\displaystyle H_{\mathrm{rot}}=$
$\displaystyle\frac{\omega_{\mathrm{r}}^{2}-\omega_{0}^{2}}{4\omega_{\mathrm{r}}}(u^{2}+v^{2})$
$\displaystyle-\frac{3\beta}{32\omega_{\mathrm{r}}}(2u^{2}v^{2}+u^{4}+v^{4})$
$\displaystyle+\frac{\lambda\omega_{0}^{2}}{8\omega_{\mathrm{r}}}\left(2uv\sin\psi+(u^{2}-v^{2})\cos\psi\right)$
$\displaystyle+\frac{AU_{\mathrm{r}}}{2\omega_{\mathrm{r}}}\left(u\cos\theta+v\sin\theta\right).$
(4)

This Hamiltonian $H_{\mathrm{rot}}$ is a function in the rotating phase space
spanned by $u$ and $v$ [20], and is the quantity we reconstruct in this paper.

The addition of dissipation in Eq. (3) is crucial. In the presence of
dissipation, a trajectory is no longer confined to a single closed loop, but
samples different energy elevations, ending up in one of the stable stationary
states of the system, see Fig. 1(b). This allows us to probe
$H_{\mathrm{rot}}$ by experimentally measuring $u(t)$ and $v(t)$ using
rotating-frame ringdown measurements. By initializing the system in an initial
state $(u_{\mathrm{i}},v_{\mathrm{i}})$ and letting it evolve to a final state
$(u_{\mathrm{f}},v_{\mathrm{f}})$, we can extract the change in the
Hamiltonian $\Delta H_{\mathrm{rot}}$ at any point $(u_{j},v_{j})$ along this
ringdown’s trajectory. We isolate $H_{\mathrm{rot}}$ in Eqs. (3) and integrate
over the slow coordinates from $(u_{\mathrm{i}},v_{\mathrm{i}})$ to
$(u_{j},v_{j})$ to obtain

$\displaystyle\Delta
H_{\mathrm{rot}}=\int_{v_{\mathrm{i}}}^{v_{j}}\left[\frac{du}{dt}+\frac{\Gamma}{2}u\right]dv-\int_{u_{\mathrm{i}}}^{u_{j}}\left[\frac{dv}{dt}+\frac{\Gamma}{2}v\right]du.$
(5)

In practice, the values $(u_{j},v_{j})$ are measured in discrete steps,
allowing to compute $\Delta H_{\mathrm{rot}}$ only at these points. Measuring
multiple ringdowns with different initial $(u_{\mathrm{i}},v_{\mathrm{i}})$,
$\Delta H_{\mathrm{rot}}$ can be deterministically sampled and reconstructed
over a large area of phase space, with a resolution limited by the measurement
uncertainty or fluctuations (e.g. thermal or quantum noise) in $u$ and $v$.
Note that Eq. (5) does not provide the relative change of $\Delta
H_{\mathrm{rot}}$ between different ringdown measurements.

To compare different traces, we make use of the fact that
$H_{\mathrm{rot}}(u_{\mathrm{f}},v_{\mathrm{f}})$ should be single-valued at
stationary points. This means that all traces sharing the same final
coordinates $(u_{\mathrm{f}},v_{\mathrm{f}})$ should have the same final value
$H_{\mathrm{rot}}(u_{\mathrm{f}},v_{\mathrm{f}})$. We thus find the relative
Hamiltonian offset between different ringdown traces sharing the same
$(u_{\mathrm{f}},v_{\mathrm{f}})$ by comparing
$H_{\mathrm{rot}}(u_{\mathrm{f}},v_{\mathrm{f}})$. Finding the offset between
traces is more complicated when they have different values of
$(u_{\mathrm{f}},v_{\mathrm{f}})$, i.e., they do not share a common end point.
Here, the offset can be calculated by making the Hamiltonian continuous, i.e.
by finding the offset that minimizes the difference between nearby starting
points $(u_{\mathrm{i}},v_{\mathrm{i}})$ that decay into different final
$(u_{\mathrm{f}},v_{\mathrm{f}})$, see Appendix C.

Note that a Hamiltonian reconstruction analogous to Eq. (5) can also be
performed in the nonrotating frame, i.e., for a dissipative Eq. (2). We
concentrate here on the case of a rotating-frame Hamiltonian to be in line
with the theory. In the following, we will test the Hamiltonian reconstruction
before discussing what information we can extract from it.

##  IV Harmonic Oscillator Case

As a first demonstration, we reconstruct the rotating-frame Hamiltonian of a
damped harmonic oscillator. In Figs. 1(b) and (c), we show multiple measured
trajectories for the case $\lambda=0$, and with amplitudes that are small
enough to neglect the effect of the Duffing non-linearity, i.e. only
considering the first line in Eq. (III). For each of those ringdown
trajectories, we first displace the resonator in phase space using a near-
resonant drive with fixed $U_{\mathrm{r}}$ and $\omega_{\mathrm{r}}$, and with
an individually selected $\theta$. In a frame rotating at
$\omega_{\mathrm{r}}$, the corresponding initial resonator coordinates
$(u_{\mathrm{i}},v_{\mathrm{i}})$, shown as dots, are stationary under the
resonant drive. Then, we switch off the drive ($U_{\mathrm{r}}=F=0$), and
track the ringdown trajectory while the system decays to the state
$(u_{\mathrm{f}}=0,v_{\mathrm{f}}=0)$. During this decay, the angle of the
state in phase space evolves in time since the rotating frame’s frequency is
detuned from the resonance frequency [21]. For
$\omega_{\mathrm{r}}-\omega_{0}>0$ in Fig. 1(b), the trajectories spiral
towards $(0,0)$ with a clockwise orientation, while for
$\omega_{\mathrm{r}}-\omega_{0}<0$ in Fig. 1(c), they spiral in
counterclockwise direction.

We reconstruct the rotating-frame Hamiltonian $H_{\mathrm{rot}}$ from the
ringdown trajectories using Eq. (5). For the reconstructions in Fig. 1(d) and
(e), the resolution was improved by using more ringdown traces than shown in
Figs. 1(b) and (c). The resulting $H_{\mathrm{rot}}$ is a paraboloid whose
sign of curvature depends on $\omega_{\mathrm{r}}-\omega_{0}$. This is
predicted in Eq. (III) and reflects the fact that $H_{\mathrm{rot}}$ describes
energy relative to a rotating frame. Indeed, for
$\omega_{\mathrm{r}}-\omega_{0}=0$, the rotating potential would be entirely
flat. Note that the origin $(0,0)$, which appears as a maximum of
$H_{\mathrm{rot}}$ in Fig. 1(e), remains the only stable solution of the
system irrespective of the detuning. As a (local) maximum can appear as a
stationary state, this example demonstrates that the rotating-frame
Hamiltonian cannot be interpreted as a simple energy function as is the case
in the non-rotating frame [48]. After successfully reconstructing the rotating
frame potential of a harmonic oscillator, we now apply our method to a more
complex system.

##  V Parametric Oscillator Case

The reconstruction of $H_{\mathrm{rot}}$ can also be applied to driven
nonlinear systems with multiple stable solutions. Here, we demonstrate this
principle on the example of the electromechanical resonator described in Sec.
II when subjected to a parametric drive, cf. Eq. (III) and Fig. 1(a).

We first analyze a single rotating-frame ringdown into an out-of-equilibrium
stationary state in phase space. In Fig. 2, we study an example trajectory of
our system in the presence of parametric driving beyond the instability
threshold ($\lambda>\lambda_{\mathrm{th}}$) [48], where two time-translation
symmetry-broken high-amplitude solutions are stabilized due to the interplay
between the drive and the nonlinearity, see Appendix A. As shown schematically
in Fig. 2(a), the system is initialized in $(u_{\mathrm{i}},v_{\mathrm{i}})$
by a near-resonant force at $\omega_{\mathrm{r}}$, followed by a parametric
drive tone whose amplitude $\lambda\propto U_{\mathrm{p}}$ and frequency
$2\omega_{\mathrm{r}}$ define the stationary solutions of the system in the
readout frame rotating at frequency $\omega_{\mathrm{r}}$ [50, 51, 48]. In
Fig. 2(b), we show a typical ringdown trajectory of our resonator in phase
space, noting that $(u_{\mathrm{f}}\neq 0,v_{\mathrm{f}}\neq 0)$.

Figure 2: Ringdown of a parametric oscillator. (a) For a ringdown measurement,
the resonator is first displaced with a near-resonant drive
$U_{\mathrm{in}}(\omega_{\mathrm{r}})$ with
$\omega_{\mathrm{r}}\approx\omega_{0}$. This resonant drive is turned off and
a parametric drive $U_{\mathrm{in}}(2\omega_{\mathrm{r}})$ is immediately
turned on, while the slowly-varying (i.e., demodulated at
$\omega_{\mathrm{r}}$) quadratures $u$ and $v$ are read out. (b) Ringdown in
phase space for fixed parametric drive $U_{\mathrm{p}}=10$ V ($\lambda\approx
10^{-3}$) and frequency $(\omega_{\mathrm{r}}-\omega_{0})/2\pi\approx-253$ Hz.
The dark blue curve is the predicted ringdown with the linearized model in
Appendix B [Eq. (19)]. (c) In panels (i) the in-phase $\tilde{u}$ and out-of-
phase $\tilde{v}$ displaced quadratures, (ii) amplitude
$\tilde{R}=\sqrt{\tilde{u}^{2}+\tilde{v}^{2}}$, and (iii) instantaneous
rotating-frame frequency $\tilde{\omega}=\partial_{t}\tilde{\phi}/2\pi$ with
phase $\tilde{\phi}=\arctan(\tilde{v}/\tilde{u})$ are plotted as a function of
time for the same ringdown as in (b). These measurements are referenced to the
attractor in which they ringdown into, i.e., $\tilde{u}\equiv
u-u_{\mathrm{f}}$ with
$u_{\mathrm{f}}\approx$-538\text{\,}\mathrm{\SIUnitSymbolMicro V}$$ the final
value of $u$ (the position of the attractor), and similarly for $\tilde{v}$
with $v_{\mathrm{f}}\approx$-572\text{\,}\mathrm{\SIUnitSymbolMicro V}$$.

To facilitate the Hamiltonian reconstruction, we define displaced quadratures
$\tilde{u}\equiv u-u_{\mathrm{f}}$ and $\tilde{v}\equiv v-v_{\mathrm{f}}$,
with the final coordinates of the ringdown $(u_{\mathrm{f}},v_{\mathrm{f}})$
(i.e, the attractor position). In this shifted reference frame, the relaxation
into the stable solution resembles a ringdown process in an equilibrium
system, with $\tilde{u}$ and $\tilde{v}$ performing damped oscillations
towards a final state $(\tilde{u}_{\mathrm{f}}=0,\tilde{v}_{\mathrm{f}}=0)$,
see Fig. 2(c)(i). The damping rate associated with these trajectories is
quantified by the same $\Gamma$ as we introduced in Eq. (1), while the
rotation around $(u_{\mathrm{f}},v_{\mathrm{f}})$ depends on
$H_{\mathrm{rot}}$. The dynamics of these shifted quadratures are well
captured by linearizing the equations of motion [Eqs. (3)] around the
attractors, see derivation in Appendix B.

Figure 3: Ringdowns and Hamiltonian reconstruction in parametric phase space.
(a) Multiple ringdowns in phase space under different parametric drive
strenght and frequency. Each ringdown’s starting point is indicated by a dot,
and color-coded for the attractor in which it rings down: purple the zero-
amplitude, orange the top-right, and green the bottom-left attractor. The
attractors extracted from the ringdown measurements are indicated by black
dots. The parametric drive is (i) below threshold, (ii) above threshold, and
(iii) above threshold and far detuned, cf. Appendix A for the locations of the
measurements relative to the device’s phase diagram. (b) Rotating-frame
Hamiltonian reconstruction from the ringdown measurements in (a), cf. Eq. (5).
(c) Theoretical plots of the rotating-frame Hamiltonian [Eq. (III)] using
potential angle (i) $\psi\approx-2.1$ rad, (ii) $\psi\approx-1.76$ rad, and
(iii) $\psi\approx-1.86$ rad as the only fit parameter, while fixing the
independently measured mechanical parameters, see Appendix A. Note that the
attractors are not necessarily at the minima of the quasi-Hamiltonian due to
dissipation.

In contrast to undriven relaxation processes in the lab frame $(x,p)$, the
amplitude $\tilde{R}=\sqrt{\tilde{u}^{2}+\tilde{v}^{2}}$ does not decay
monotonically. Instead, we observe oscillations imposed on top of the
exponential decay in Fig. 2(c)(ii). These oscillations stem from the fact that
the shape of $H_{\mathrm{rot}}$ around $(u_{\mathrm{f}},v_{\mathrm{f}})$ is
not rotationally symmetric, such that the system samples different
$\partial_{u}H_{\mathrm{rot}}$ and $\partial_{v}H_{\mathrm{rot}}$ as it moves
around the Hamiltonian landscape. We emphasize that despite the fact that we
observe oscillations growing in amplitude at certain times, this does not
violate any law of conservation, as we are dealing with a driven system. The
non-rotationally symmetric Hamiltonian also manifests in Fig. 2(c)(iii), where
an instantaneous rotating-frame frequency
$\tilde{\omega}=\partial_{t}\tilde{\phi}/2\pi$ is defined by the phase
$\tilde{\phi}=\arctan(\tilde{v}/\tilde{u})$ relative to the attractor.

There are clear oscillations in the frequency $\tilde{\omega}$ due to the
resonator sampling different $\partial_{u}H_{\mathrm{rot}}$ and
$\partial_{v}H_{\mathrm{rot}}$ as it moves in phase space. This behavior is
also well captured by the averaged and linearized dynamics near the attractors
(cf. Eqs. (3) and Appendix B). Note that for smaller parametric drives or for
dynamics further away from the attractors, one needs to consider a model that
goes beyond the linearization we employ (i.e. expand the equations to higher-
order terms) to fully capture the dynamics.

In a next step, we measure multiple trajectories with various initialization
conditions, but fixed parametric drive strength and frequency, allowing us to
reconstruct the rotating frame Hamiltonian $H_{\mathrm{rot}}$ of our driven
nonlinear system. In Fig. 3(a), we show such sets of trajectories for three
different parametric drives: (i) below threshold with a single squeezed state,
(ii) above threshold with two stable phase states, and (iii) at large driving
and large detuning with a combination of phase states and a zero-amplitude
state (see Appendix A for a phase diagram of the device). Each ringdown is
color-coded by the state it eventually approaches, allowing us to identify the
corresponding attractor pools. From ringdown measurements, we can thus
directly identify the number of stable solutions in the given phase-space
area, as well as the separatrices of the system – the regions where the color
of nearby ringdown measurements changes. In addition, we obtain a
visualization of the stream flow in phase space, as governed by Eq. (3). We
stress that the main assumption of this reconstruction is that the damping is
linear, cf. Eqs. (3), though these equations could be adapted to include other
types of damping. The reconstruction does not assume linearized dynamics.

The reconstructed Hamiltonians $H_{\mathrm{rot}}$ are shown in Fig. 3(b)
alongside theoretical Hamiltonians in Fig. 3(c), as calculated from Eq. (III)
with independently measured parameters $\omega_{0}$, $\beta$, and $\lambda$,
see Appendix A. We find excellent qualitative agreement between measurement
and theory. Crucially, the quality and resolution of the Hamiltonian
reconstruction is consistently high over the entire sampled phase space, which
would be hard to achieve with statistical methods [6, 7, 8, 9, 10, 11, 12,
13]. In the reconstructed $H_{\mathrm{rot}}$, we can clearly see the
appearance of one, two, and three stable states in the three cases (i)-(iii),
respectively, indicating different phases of the driven system [27, 28, 52].
Our reconstruction method confirms the theoretical prediction that the
stationary solutions can qualitatively differ in the rotating frame; while the
phase states at finite amplitude in the cases (ii) and (iii) are marked by a
minimum in $H_{\mathrm{rot}}$, the stable state appearing at $(u=0,v=0)$ in
(iii) corresponds to a maximum, signalling a fundamentally different type of
solution, i.e., a dissipation-stabilized state [29, 30]. In the following, we
analyze this difference using the symplectic norm of the individual solutions,
and we show that this quantity yields valuable insights into the behavior of
driven-dissipative systems.

##  VI Symplectic Norm

In sections IV and V, we successfully reconstructed rotating-frame
Hamiltonians $H_{\mathrm{rot}}$ using ringdown measurements. One of the most
prominent distinctions of a rotating-frame Hamiltonian, compared with a
Hamiltonian in a laboratory frame such as in Eq. (2), is that both maxima and
minima of $H_{\mathrm{rot}}$ can constitute stable oscillation states of the
system, unlike what we expect from the minimal action principle in equilibrium
systems. This counter-intuitive feature can be clearly observed in row (iii)
of Fig. 3, where two minima at finite amplitudes are separated by a stable
local maximum at $(u=v=0)$. Importantly, the minima in Fig. 3 are stable due
to the nonlinearity $\beta$, while the maximum is stabilized by dissipation
[30]. This fundamental difference, however, is difficult to quantify in
standard measurements, such as frequency sweeps and stability diagrams [41,
53].

A method tailored to classify and distinguish minima and maxima in
$H_{\mathrm{rot}}$ is the symplectic norm $ds^{2}$ [29, 30]. This quantity
indicates if the excitations (i.e., out-of-equilibrium phonons) of a system
around an attractor are more hole-like or particle-like: when an excitation
with a negative (positive) symplectic norm is created on top of a stationary
system, it reduces (increases) the energy of the system relative to the
rotating frame, that is, relative to an excitation at the driving frequency
$\omega_{\mathrm{r}}$. This difference manifests in the Hamiltonian:
stationary solutions with a negative (positive) symplectic norm appear as
maxima (minima) in $H_{\mathrm{rot}}$ and are formally associated with a hole-
like (particle-like) excitation, see Appendix B. We can therefore use the
reconstructed Hamiltonians in Fig. 3 to directly determine the symplectic norm
of the stable oscillations states, and to classify the corresponding
excitations.

A careful study of Fig. 3 reveals additional information. Namely, the
trajectories leading to Hamiltonian extrema assume clockwise or counter-
clockwise rotations. For our system, we show in Appendix B that the symplectic
norm is directly linked to the sense of rotation of the trajectories close to
an attractor. This allows us to extract the symplectic norm $ds^{2}$ of the
different attractors directly from the measured ringdown trajectories, even
without a full reconstruction. To capture this observation mathematically, we
define a correlator characterizing the direction of rotation of a ringdown,

$\displaystyle
G^{c}(t^{\prime}-t)=\Theta(t^{\prime}-t)\langle\tilde{v}(t)\tilde{u}(t^{\prime})-\tilde{u}(t)\tilde{v}(t^{\prime})\rangle_{t},$
(6)

with the Heavyside step-function $\Theta(t-t^{\prime})$, and
$\langle...\rangle_{t}$ indicating an average over all times $t$. This
quantity $G^{c}(t-t^{\prime})$ is a classical analogue to the quantum Green’s
function used in Refs. [54, 30, 34]. Following these works, we calculate the
corresponding spectral response $\mathcal{A}$ of a stable oscillation state

$\displaystyle\mathcal{A}(\omega)$
$\displaystyle=-2\mathrm{Im}[G^{c}(\omega)]$ (7)
$\displaystyle=\frac{|\zeta|^{2}ds^{2}}{2}\left[\frac{1}{(\omega-\omega_{\mathrm{lin}})^{2}+\frac{\Gamma^{2}}{4}}-\frac{1}{(\omega+\omega_{\mathrm{lin}})^{2}+\frac{\Gamma^{2}}{4}}\right],$

where $G^{c}(\omega)$ is the Fourier transform of $G^{c}(t-t^{\prime})$,
$\zeta$ is a constant related to the ringdown starting conditions, and
$\omega_{\mathrm{lin}}$ is the oscillation frequency of the quadratures. The
linear response function in the second line is only valid near the attractor
of interest.

Looking at $\mathcal{A}$, we see that the spectral response has peaks at
$\pm\omega_{\mathrm{lin}}$, with width $\Gamma$, and an overall sign which is
determined by the symplectic norm $ds^{2}$. Therefore, if the resonator rings
down with a counter-clockwise (clockwise) rotation in phase space, the final
steady state of the system has a negative (positive) symplectic norm and a
negative (positive) peak in $\mathcal{A}$.

Figure 4: Spectral response of a parametric oscillator. With $\tilde{u}$ and
$\tilde{v}$ obtained during a ringdown measurement, we can extract the
spectral response $\mathcal{A}$ using Eqs. (6) and (7) (dots) [cf. Figs. (2)
and (3)]. These spectral responses are in good agreement with the fitted
theoretical curves of Eq. (7), with $\omega_{0}$ left as a free parameter to
account for small frequency drifts. All measurements are done with the same
parametric drive strenght $U_{\mathrm{p}}=10$ V ($\lambda\approx 10^{-4}$),
and different demodulation frequencies $\omega_{\mathrm{r}}$. The detuning is
$(\omega_{\mathrm{r}}-\omega_{0})/2\pi=170$ Hz (orange) and
$(\omega_{\mathrm{r}}-\omega_{0})/2\pi=-222$ Hz (green) for ringdowns into
high amplitude states, and $(\omega_{\mathrm{r}}-\omega_{0})/2\pi=-865$ Hz
(purple) for the zero amplitude state.

With this knowledge, we can extract the sign of $ds^{2}$ from the spectral
response of ringdown measurements, as shown in Fig. 4. We record $\tilde{u}$
and $\tilde{v}$ during one trajectory, calculate the correlator
$G^{c}(t-t^{\prime})$ [Eq. (6)], take its Fourier transform, and compute the
spectral response $\mathcal{A}$ according to Eq. (7). In Fig. 4, we present
the measured result for each of the three stable states of a parametric
oscillator at different detunings. We compare these measured spectral
responses with those expected from theory, calculated using independently
measured parameters and the linearized slow-flow equations [Eq. (7)], and
obtain very good agreement. As expected, for a stable solution corresponding
to a maximum in $H_{\mathrm{rot}}$ (purple curve), the spectral response is a
dip, while the two other solutions, which are minima in $H_{\mathrm{rot}}$
(orange and green curves), correspond to peaks in the spectral response. These
measurements show the link between the orientation of rotation, the symplectic
norm, the spectral function, and the rotating frame Hamiltonian
$H_{\mathrm{rot}}$.

We emphasize that the notion of maxima and minima of the rotating quasienergy
potential depend on the chosen rotating frame frequency, and that stabilized
maxima are a manifestation of out-of-equilibrium stationary states. As such,
the fact that the excitations on top of the maxima are hole-like implies that
their dynamics are slower than the clock in the rotating frame as manifested
by the chirality of their ringdown. In other words, we can think of the
response in this case as non-causal relative to the clock. This is analogous
to how antiparticles, that exhibit a negative mass dispersion, behave in
relativistic quantum mechanics [55].

##  VII Outlook

We report a precise method to reconstruct the Hamiltonian of a system via
ringdown measurements. The method allows for a full characterization of the
energy potential, including multiple stable solutions, saddle points, and
attractor pools. The method is particularly suited for studies in the growing
field of driven-dissipative nonlinear systems, where a Hamiltonian
characterization from first principles is often very difficult. Furthermore,
it bestows the ability to characterize the symplectic norm of different stable
oscillation solutions in the rotating frame with a connection to relativistic
quantum mechanics and causality. We expect that this approach will allow the
classification of a broad variety of systems, including nanomechanics,
superconducting circuits, light-matter systems, and nonlinear optics.

## Acknowledgements

We acknowledge Christian Marti and Sebastián Guerrero for help building the
measurement setup. We thank Nicholas E. Bousse and T. W. Kenny for providing
the MEMS device, and Matteo Fadel and Tobias Donner for useful discussions. V.
D. acknowledges support from the ETH Zurich Postdoctoral Fellowship Grant No.
23-1 FEL-023. O. Z. acknowledges funding from the Deutsche
Forschungsgemeinschaft (DFG) via project number 449653034 and through SFB143.
A. E. and O. Z. acknowledge financial support from the Swiss National Science
Foundation (SNSF) through the Sinergia Grant No. CRSII5_206008/1.

##  Appendix A Experimental Characterization

###  A.1 Parametric Sweeps

In this section, we describe how we measure the parametric response of our
resonator in order to extract its resonance frequency $\omega_{0}$, energy
decay rate $\Gamma$, and Duffing nonlinearity $\beta$.

We send a parametric tone
$U_{\mathrm{in}}=U_{\mathrm{p}}\cos(2\omega_{\mathrm{r}}t)$ to our device and
read out the induced displacement at frequency $\omega_{\mathrm{r}}$. In Fig.
5(a), we show the system’s response while sweeping the frequency of the
parametric drive from high to low frequencies for fixed parametric drive
voltage $U_{\mathrm{p}}=4.54$ V. By repeating this measurement for different
parametric drive voltages, we obtain the ‘phase diagram’ of our parametric
oscillator, shown in Fig. 5(b). Now, performing the same measurement but
instead sweeping the frequency from low to high frequencies (i.e., “against”
the Duffing nonlinearity), we obtain a different diagram commonly referred to
as ‘Arnold Tongue’, shown in Fig. 5(c).

Comparing with theoretical predictions [48], Figs. 5(b) and (c) allow us to
read out the different phases and number of solutions of our resonator
depending on the parametric drive strength and frequency. Indeed, the outline
of the Arnold tongue, shown as a red dashed line, indicates the (frequency-
dependent) parametric threshold, above which the resonator has exactly two
stable states (the parametric phase states), both of which have finite-
amplitude but opposite phases. The orange region outside the Arnold tongue in
Fig. 5(b) features a zero-amplitude state in addition to the two parametric
phase states. Which state is selected depends on the initial condition of the
resonator. Finally, in the white region in Fig. 5(b) only the zero-amplitude
state is stable.

To find the equation predicting the parametric response of our device, we use
the ‘slow-flow’ equations [cf. Eqs. (3) and (B.1)], and find the steady state
amplitude response $R=\sqrt{u^{2}+v^{2}}$ by setting $\dot{u}=\dot{v}=0$,
yielding [48]

$\displaystyle\left[-\frac{\lambda^{2}}{4}+\left(\frac{\Gamma\omega}{\omega_{0}^{2}}\right)^{2}+\left(1-\frac{\omega^{2}}{\omega_{0}^{2}}+\frac{3}{4}\frac{\beta
R^{2}}{\omega_{0}^{2}}\right)^{2}\right]R^{2}=0\,.$ (8)

Equation (8) has the trivial solution $R=0$. For the case $R\neq 0$, we can
divide Eq. (8) by $R^{2}$ to obtain the non-trivial solutions

$\displaystyle
R(\omega)=\left(\frac{4\omega_{0}^{2}}{3\beta}\left[\left(\frac{\omega^{2}}{\omega_{0}^{2}}-1\right)\pm\left(\frac{\lambda^{2}}{4}-\frac{\Gamma^{2}\omega^{2}}{\omega_{0}^{4}}\right)^{1/2}\right]\right)^{1/2}\,,$
(9)

when the radicands of both square roots are positive. Thus, the parametric
amplitude response allows (amongst other things) to extract the Duffing
nonlinearity.

Figure 5: Parametric response. (a) Response of our resonator to a parametric
drive $U_{\mathrm{in}}=U_{\mathrm{p}}\cos(2\omega_{\mathrm{r}}t)$ with voltage
amplitude $U_{\mathrm{in}}=4.54$ V. The frequency is swept from high to low
values. (b) Response to a parametric drive as in (a) for different values of
$U_{\mathrm{p}}$. The dashed grey line indicates the measurement of panel (a),
the red dashed line corresponds to the fitted ‘Arnold tongue’ outline [see
panel (c)], while the dots indicate the parameters used in Fig. 3. (c)
Response to a parametric drive as in (a) for different values of
$U_{\mathrm{p}}$ when sweeping from low to high frequencies.

The outline of the Arnold tongue corresponds to the limit $R\to 0$ in Eq. (9),
leading to

$\displaystyle\left(\omega^{2}-\omega_{0}^{2}\right)^{2}-\frac{\lambda^{2}\omega_{0}^{4}}{4}+\Gamma^{2}\omega^{2}=0\,.$
(10)

Solving for the parametric drive yields

$\displaystyle\lambda=2\sqrt{\frac{\Gamma^{2}\omega^{2}}{\omega_{0}^{4}}+\left(1-\frac{\omega^{2}}{\omega_{0}^{2}}\right)^{2}},$
(11)

where we kept the (physical) positive solution.

To extract our device parameters, we start by fitting Eq. (11) to the outline
of the Arnold tongue, see red dashed line in Fig. 5(c). For the fit, we
replace $\lambda=U_{\mathrm{p}}/C$, where $C$ is a constant converting the
unitless parametric drive strength $\lambda$ to the applied parametric voltage
$U_{\mathrm{p}}$. Doing so, we obtain the resonance frequency
$\omega_{0}/2\pi=1.1198294(3)$ MHz and the energy decay rate
$\Gamma/2\pi=112(2)$ Hz of our resonator, as well as the conversion constant
$C=9.88(5)\times 10^{3}$ V.

In a second step, we use the measured parametric response in Fig. 5(a) and fit
it to Eq. (9) in order to extract the Duffing nonlinearity
$\beta=-9.894(4)\times 10^{16}$ (V$\cdot$s)-2, while fixing $\omega_{0}$,
$\Gamma$, and $C$ to the previously fitted values.

###  A.2 Resonant Sweep

We now extract the factor $A$, converting from an input voltage
$U_{\mathrm{in}}=U_{\mathrm{r}}\cos(\omega_{\mathrm{r}}t)$ to an applied force
$F$ (in units Vs-2), from the resonant response of our resonator.

To calibrate this, we consider the EOM for our resonator under near-resonant
(but small, i.e., neglecting the Duffing term) drive [cf. Eq. (1)],

$\displaystyle\ddot{x}(t)+\omega_{0}^{2}x(t)+\Gamma\dot{x}(t)=AU_{\mathrm{in}}(t),$
(12)

where we replaced the forcing term by $F=AU_{\mathrm{in}}$. Taking the Fourier
transform of Eq. (12) yields

$\displaystyle
x(\omega)=AU_{\mathrm{in}}(\omega)\left[\omega_{0}^{2}-\omega^{2}-i\omega\Gamma\right]^{-1},$
(13)

which provides a direct link between the voltage $U_{\mathrm{in}}$ applied to
our device and the voltage $U_{\mathrm{out}}\equiv x$ we read out.

We can thus extract $A$ by sweeping the frequency $\omega_{\mathrm{r}}$ of a
driving tone across the mechanical resonance of our resonator and measuring
its response. In Fig. 6, we fit the measured response with Eq. (13) and obtain
$A=16.22(1)\times 10^{6}$ s-2.

Figure 6: Linear resonant response. An input voltage
$U_{\mathrm{in}}=U_{\mathrm{r}}\cos(\omega_{\mathrm{r}}t)$ with
$U_{\mathrm{r}}=10$ mV is swept across the resonance frequency and the
amplitude response $R=\sqrt{u^{2}+v^{2}}$ at frequency $\omega_{\mathrm{r}}$
is measured (dots). The data is fitted (line) according to Eq. (13) with fixed
$\Gamma/2\pi=112$ Hz, while $\omega_{0}$ is left as a free parameter to
account for small frequency drifts. We extract a value $A=16.22(1)\times
10^{6}$ s-2.

##  Appendix B Linearized Model

In this Appendix, we describe the method we use to obtain the EOMs for the
rotating frame coordinates $u$ and $v$, and the link between the symplectic
norm and the ringdown trajectories. We first obtain the steady state solutions
under parametric drive, then we linearize the EOMs in their vicinity. We solve
the linearized EOMs to get analytical expressions for the resonators rotation
frequency and direction around the attractors. The latter reveals a direct
link between the ringdown rotation and the symplectic norm of the
corresponding final steady state. Finally, we show how to extract the
symplectic norm from a ringdown measurement using the spectral response
function.

###  B.1 Averaging Method

As described in the main text, starting from Eq. (1), we move to a rotating
frame with out-of-phase quadratures $u(t)$ and $v(t)$, defined via
$x(t)=u(t)\cos(\omega_{\mathrm{r}}t)-v(t)\sin(\omega_{\mathrm{r}}t)$. Assuming
that $u(t)$ and $v(t)$ vary slowly in time compared to the mechanical
oscillations with frequency $\omega_{\mathrm{r}}$, we average the EOMs over
the time $2\pi/\omega_{\mathrm{r}}$ to obtain approximated EOMs for $u(t)$ and
$v(t)$ [48],

$\displaystyle\dot{u}$ $\displaystyle=-\frac{\Gamma}{2}u-\delta
v-\frac{3\beta}{8\omega_{\mathrm{r}}}\left(u^{2}v+v^{3}\right)-\frac{\lambda\omega_{0}^{2}}{4\omega_{\mathrm{r}}}v\,,$
$\displaystyle\dot{v}$ $\displaystyle=-\frac{\Gamma}{2}v+\delta
u+\frac{3\beta}{8\omega_{\mathrm{r}}}\left(u^{3}+uv^{2}\right)-\frac{\lambda\omega_{0}^{2}}{4\omega_{\mathrm{r}}}u\,,$
(14)

where we introduce the detuning
$\delta\equiv(\omega_{0}^{2}-\omega_{\mathrm{r}}^{2})/2\omega_{\mathrm{r}}$.
For simplicity, we also set $\psi=0$ when moving from Eq. (1) to Eq. (B.1).
This phase $\psi$ rotates the potential in phase space around the origin with
an angle of $\psi/2$, so the solutions for $u$ or $v$ calculated for $\psi=0$
in this section can be rotated by an angle of $\psi/2$ to match the
experimental results.

###  B.2 Linearization of the EOMs

To predict the motion of the resonator near an attractor, we linearize the
slow-flow equations [Eqs. (B.1)] near the attractors [48, 56]. To do so, we
first determine the coordinates of the attractors (i.e., the stationary states
of the system) in the rotating frame by setting $\dot{u}=\dot{v}=0$ in Eq.
(B.1). Solving the resulting coupled polynomial equations leads to five
possible solutions [57, 58]:

$\displaystyle(u_{\mathrm{f}},v_{\mathrm{f}})_{1}$ $\displaystyle=(0,0)\,,$
$\displaystyle(u_{\mathrm{f}},v_{\mathrm{f}})_{2,3}$ $\displaystyle=\left(\pm
u_{+},\pm\frac{u_{+}(3\beta\lambda\omega_{0}^{2}u_{+}^{2}-2\lambda^{2}\omega_{0}^{4}+4\Gamma^{2}\omega_{\mathrm{r}}^{2}+4\lambda\omega_{0}^{4}-4\lambda\omega_{0}^{2}\omega_{\mathrm{r}}^{2}}{2\Gamma\omega_{\mathrm{r}}\left(2\omega_{\mathrm{r}}^{2}-2\omega_{0}^{2}+\lambda\omega_{0}^{2}\right)}\right)\,,$
$\displaystyle(u_{\mathrm{f}},v_{\mathrm{f}})_{4,5}$ $\displaystyle=\left(\pm
u_{-},\pm\frac{u_{-}(3\beta\lambda\omega_{0}^{2}u_{-}^{2}-2\lambda^{2}\omega_{0}^{4}+4\Gamma^{2}\omega_{\mathrm{r}}^{2}+4\lambda\omega_{0}^{4}-4\lambda\omega_{0}^{2}\omega_{\mathrm{r}}^{2}}{2\Gamma\omega_{\mathrm{r}}\left(2\omega_{\mathrm{r}}^{2}-2\omega_{0}^{2}+\lambda\omega_{0}^{2}\right)}\right)\,,$
(15)

where

$\displaystyle u_{\pm}$
$\displaystyle=\left[\frac{\beta\lambda\omega_{0}^{2}\left(2\lambda\omega_{0}^{2}\omega_{\mathrm{r}}^{2}+\lambda^{2}\omega_{0}^{4}-2\lambda\omega_{0}^{4}-4\Gamma^{2}\omega_{\mathrm{r}}^{2}\right)\pm\sqrt{\beta^{2}\lambda^{2}\omega_{0}^{4}\left(2\omega_{\mathrm{r}}^{2}-2\omega_{0}^{2}+\lambda\omega_{0}^{2}\right)^{2}\left(\lambda^{2}\omega_{0}^{4}-4\Gamma^{2}\omega_{\mathrm{r}}^{2}\right)}}{3\beta^{2}\lambda^{2}\omega_{0}^{4}}\right]^{1/2}\,.$
(16)

For any given set of parameters, a maximum of three of these solutions are
stable, i.e., they act as attractors, cf. Fig. 3. To see if a solution is
stable, we now introduce the displaced quadratures
$(\tilde{u},\tilde{v})_{k}\equiv(u-u_{\mathrm{f}},v-v_{\mathrm{f}})_{k}$,
relative to a given solution position $(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$,
and linearize Eq. (B.1) around one attractor by neglecting higher-order terms
in $\tilde{u}$ and $\tilde{v}$. This procedure yields

$\displaystyle\begin{pmatrix}\dot{\tilde{u}}\\\ \dot{\tilde{v}}\end{pmatrix}$
$\displaystyle=\underbrace{\left.\begin{pmatrix}\frac{\partial\dot{u}}{\partial
u}&\frac{\partial\dot{u}}{\partial v}\\\\[4.26773pt]
\frac{\partial\dot{v}}{\partial u}&\frac{\partial\dot{v}}{\partial
v}\end{pmatrix}\right|_{(u,v)=(u_{\mathrm{f}},v_{\mathrm{f}})_{k}}}_{\mathbf{J}_{\mathrm{f}}}\begin{pmatrix}\tilde{u}\\\
\tilde{v}\end{pmatrix}\,.$ (17)

Here, $\mathbf{J}_{\mathrm{f}}$ is the Jacobian of the slow-flow equations
evaluated around the solution position $(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$,
which takes the form

$\displaystyle\mathbf{J}_{\mathrm{f}}$
$\displaystyle=\begin{pmatrix}-\frac{\Gamma}{2}-\frac{3\beta
u_{\mathrm{f}}v_{\mathrm{f}}}{4\omega_{\mathrm{f}}}&-\delta-\frac{3\beta\left(u_{\mathrm{f}}^{2}+3v_{\mathrm{f}}^{2}\right)}{8\omega_{\mathrm{r}}}-\frac{\lambda\omega_{0}^{2}}{4\omega_{\mathrm{r}}}\\\
\delta+\frac{3\beta\left(3u_{\mathrm{f}}^{2}+v_{\mathrm{f}}^{2}\right)}{8\omega_{\mathrm{r}}}-\frac{\lambda\omega_{0}^{2}}{4\omega_{\mathrm{r}}}&-\frac{\Gamma}{2}+\frac{3\beta
u_{\mathrm{f}}v_{\mathrm{f}}}{4\omega_{\mathrm{r}}}\end{pmatrix}.$ (18)

The Jacobian $\mathbf{J}_{\mathrm{f}}$ describes the linearized forces acting
in the rotating frame near a given solution. Solving the linear first-order
differential equation given by Eq. (17) leads to

$\displaystyle\begin{pmatrix}\tilde{u}\\\
\tilde{v}\end{pmatrix}=\zeta\mathbf{w}_{+}e^{\mu_{+}t}+\zeta^{*}\mathbf{w}_{-}e^{\mu_{-}t}\,.$
(19)

Here,

$\displaystyle\zeta=$
$\displaystyle\frac{1}{2}(v_{\mathrm{i}}-v_{\mathrm{f}})-i\bigg{(}\frac{3\beta
u_{\mathrm{f}}v_{\mathrm{f}}}{8\omega_{\mathrm{r}}\omega_{\mathrm{lin}}}(v_{\mathrm{i}}-v_{\mathrm{f}})$
(20) $\displaystyle+\frac{(9\beta u_{\mathrm{f}}^{2}+3\beta
v_{\mathrm{f}}^{2}+4\omega_{0}^{2}-4\omega_{\mathrm{r}}^{2}-2\lambda\omega_{0}^{2})}{16\omega_{\mathrm{r}}\omega_{\mathrm{lin}}}(u_{\mathrm{i}}-u_{\mathrm{f}})\bigg{)}$

is a constant which depends on the initial position
$(u_{\mathrm{i}},v_{\mathrm{i}})$ of the resonator in the rotating frame.
Furthermore, we use in Eq. (19) the Jacobian’s eigenvectors

$\displaystyle\mathbf{w}_{\pm}=\begin{pmatrix}\frac{-6\beta
u_{\mathrm{f}}v_{\mathrm{f}}\pm
8i\omega_{\mathrm{r}}\omega_{\mathrm{lin}}}{9\beta u_{\mathrm{f}}^{2}+3\beta
v_{\mathrm{f}}^{2}+4\omega_{0}^{2}-4\omega_{\mathrm{r}}^{2}-2\lambda\omega_{0}^{2}}\\\
1\end{pmatrix},$ (21)

and corresponding eigenvalues

$\displaystyle\mu_{\pm}=-\Gamma/2\pm i\omega_{\mathrm{lin}},$ (22)

with the complex frequency

$\displaystyle\omega_{\mathrm{lin}}=\frac{1}{8\omega_{\mathrm{r}}}\Bigl{[}$
$\displaystyle 27\beta^{2}u_{\mathrm{f}}^{4}+6\beta
u_{\mathrm{f}}^{2}\left(9\beta
v_{\mathrm{f}}^{2}+8\omega_{0}^{2}-8\omega_{\mathrm{r}}^{2}+2\lambda\omega_{0}^{2}\right)$
$\displaystyle+\left(3\beta
v_{\mathrm{f}}^{2}+4\omega_{0}^{2}-4\omega_{\mathrm{r}}^{2}-2\lambda\omega_{0}^{2}\right)$
$\displaystyle\times\left(9\beta
v_{\mathrm{f}}^{2}+4\omega_{0}^{2}-4\omega_{\mathrm{r}}^{2}+2\lambda\omega_{0}^{2}\right)\Bigl{]}^{1/2}\,.$
(23)

If the real part of one eigenvalue in Eq. (22) is positive, Eq. (19) diverges
for $t\rightarrow\infty$ and $(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$ is an
unstable state. If both real parts are negative, $(\tilde{u},\tilde{v})$
converge to $(0,0)$ and $(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$ is an attractor
[48, 56]. The theory results depicted in Fig. 2(b), (c)(i), and (c)(ii) are
directly calculated using Eq. (19).

The theory prediction for the frequency of the rotation around an attractor in
Fig. 2(c)(iii) is analytically calculated using

$\displaystyle\frac{1}{2\pi}\frac{d\tilde{\phi}}{dt}$
$\displaystyle=\frac{\partial}{\partial
t}\arctan\left(\frac{\tilde{v}}{\tilde{u}}\right)$
$\displaystyle=\frac{1}{2\pi}\frac{\tilde{u}(d\tilde{v}/dt)-\tilde{v}(d\tilde{u}/dt)}{\tilde{u}^{2}+\tilde{v}^{2}}$
$\displaystyle=-\frac{|\zeta|^{2}\omega_{\mathrm{lin}}}{\pi}\frac{ds^{2}}{\tilde{u}^{2}+\tilde{v}^{2}}e^{-\Gamma
t}\,,$ (24)

with the ‘symplectic norm’ defined as

$\displaystyle ds^{2}\equiv i(w_{+,1}w_{-,2}-w_{-,1}w_{+,2})\,,$ (25)

and where $w_{\pm,l}$ is the $l$-th entry of $\mathbf{w}_{\pm}$. The sign of
$ds^{2}$ therefore decides the rotational sense of the ringdown. We now show
that this quantity is identical to the symplectic norm derived for quantum
driven-dissipative systems, which allows to classify stable states in driven-
dissipative systems [29, 30, 59].

###  B.3 Symplectic Norm

To see that the symplectic norm [cf. Eq. (25)] is the same as the one defined
in quantum driven-dissipative systems, we transform our definition of the
symplectic norm to the form commonly used in quantum optics [29, 30, 59].
First, we express the rotating frame coordinates $\tilde{u}(t)$ and
$\tilde{v}(t)$ in terms of the complex coordinate $\alpha(t)$ and its complex
conjugate $\alpha^{*}(t)$ via the transformation matrix $\mathbf{S}$:

$\displaystyle\begin{pmatrix}\tilde{u}\\\ \tilde{v}\end{pmatrix}$
$\displaystyle=\underbrace{\sqrt{\frac{\hbar}{2m\omega_{0}}}\begin{pmatrix}1&1\\\
i&-i\end{pmatrix}}_{\equiv\mathbf{S}^{-1}}\begin{pmatrix}\alpha\\\
\alpha^{*}\end{pmatrix}\,.$ (26)

This basis change can be interpreted as moving to the mean-field limit of
bosonic creation and annihilation operators of a quantum harmonic oscillator
representing the quadratures $\tilde{u}(t)$ and $\tilde{v}(t)$ (i.e.,
$\alpha=\langle a\rangle$ with a bosonic annihilation operator $a$ and
expectation value $\langle...\rangle$).

Note that the $\sqrt{\hbar/2m\omega_{0}}$ prefactor is the same for both
quadratures due to our definition of $\tilde{u}$ and $\tilde{v}$ sharing the
same units. This prefactor leads to

$\displaystyle\mathbf{S}^{-1}=\frac{\hbar}{m\omega_{0}}\mathbf{S}^{\dagger},$
(27)

implying that $\mathbf{S}$ is not unitary and therefore not a norm-preserving
transformation. This will later be accounted for by scaling the vectors in the
new basis, and has no consequences for the sign of the symplectic norm, which
is the relevant quantity for us.

We use the transformation introduced in Eqs. (17) and (26) to find the
equations of motion for $\alpha$ and $\alpha^{*}$:

$\displaystyle i\,\mathbf{S}\begin{pmatrix}\dot{\tilde{u}}\\\
\dot{\tilde{v}}\end{pmatrix}$
$\displaystyle=i\,\mathbf{S}\mathbf{J}_{\mathrm{f}}\begin{pmatrix}\tilde{u}\\\
\tilde{v}\end{pmatrix}$ $\displaystyle\Rightarrow
i\begin{pmatrix}\dot{\alpha}\\\ \dot{\alpha}^{*}\end{pmatrix}$
$\displaystyle=\underbrace{i\mathbf{S}\mathbf{J}_{\mathrm{f}}\mathbf{S}^{-1}}_{\equiv\mathbf{D}}\begin{pmatrix}\alpha\\\
\alpha^{*}\end{pmatrix}\,.$ (28)

We multiplied both sides of Eq. (26) with a complex factor $i$ in order to
bring Eq. (B.3) into a form where the so-called dynamic matrix $\mathbf{D}$
[60] of the system can be directly read out. This is crucial as the symplectic
norm defined in Ref. [30] is formulated in terms of the eigenvectors of the
dynamic matrix.

The eigenvectors $\mathbf{v}_{\pm}$ of the dynamic matrix $\mathbf{D}$ can be
calculated using the eigenvectors $\mathbf{w}_{\pm}$ of the Jacobian
$\mathbf{J}_{\mathrm{f}}$:

$\displaystyle\mathbf{J}_{\mathrm{f}}\mathbf{w}_{\pm}$
$\displaystyle=\mu_{\pm}\mathbf{w}_{\pm}$ $\displaystyle\Rightarrow
i\,\mathbf{S}\mathbf{J}_{\mathrm{f}}\mathbf{S}^{-1}\mathbf{S}\mathbf{w}_{\pm}$
$\displaystyle=i\,\mu_{\pm}\mathbf{S}\mathbf{w}_{\pm}$
$\displaystyle\Rightarrow i\,\mathbf{D}\mathbf{S}\mathbf{w}_{\pm}$
$\displaystyle=i\,\mu_{\pm}\mathbf{S}\mathbf{w}_{\pm}$
$\displaystyle\Rightarrow\mathbf{D}\left(\sqrt{\frac{m\omega_{0}}{\hbar}}\mathbf{v}_{\pm}\right)$
$\displaystyle=i\mu_{\pm}\left(\sqrt{\frac{m\omega_{0}}{\hbar}}\mathbf{v}_{\pm}\right)\,.$
(29)

In the last line, we used the definition of $\mathbf{D}$ from Eq. (B.3) and
introduced an additional factor of $\sqrt{m\omega_{0}/\hbar}$ to ensure that
$|\mathbf{w}_{\pm}|=|\mathbf{v}_{\pm}|$, accounting for the fact that
$\mathbf{S}$ is not a unitary transformation, cf. Eq. (27).

Since $\mathbf{w}_{+}=\mathbf{w}_{-}^{*}$, this implies that
$\mathbf{v}_{+}=\mathbf{v}_{-}^{*}$, and we have
$\mathbf{w}_{+}^{\dagger}=\mathbf{w}_{-}^{T}$, as well as
$\mathbf{v}_{+}^{\dagger}=\mathbf{v}_{-}^{T}$ (with $T$ denoting transpose),
allowing us to rewrite the classical symplectic norm defined in Eq. 25 as

$\displaystyle ds^{2}$
$\displaystyle=\left(w_{-,1},w_{-,2}\right)\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\begin{pmatrix}w_{+,1}\\\ w_{+,2}\end{pmatrix}$
$\displaystyle=\left(w_{-,1},w_{-,2}\right)\mathbf{S}^{\dagger}\left(\mathbf{S}^{\dagger}\right)^{-1}\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\mathbf{S}^{-1}\mathbf{S}\begin{pmatrix}w_{+,1}\\\
w_{+,2}\end{pmatrix}$
$\displaystyle=\left(\mathbf{S}\mathbf{w}_{+}\right)^{\dagger}\left(\mathbf{S}^{\dagger}\right)^{-1}\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\mathbf{S}^{-1}\mathbf{S}\mathbf{w}_{+}$
$\displaystyle=\mathbf{v}_{+}^{\dagger}\mathbf{I}_{-}\mathbf{v}_{+}\,,$ (30)

with $\mathbf{I}_{-}=\mathrm{diag}(1,-1)$, and where we used
$\mathbf{v}_{+}=\mathbf{S}\mathbf{w}_{+}$. Since the symplectic norm is real
valued, one could also calculate it using $\mathbf{v}_{-}$, yielding the same
result:

$\displaystyle
ds^{2}=\left(ds^{2}\right)^{\dagger}=\left(\mathbf{v}_{+}^{\dagger}\mathbf{I}_{-}\mathbf{v}_{+}\right)^{\dagger}=\mathbf{v}_{-}^{\dagger}\mathbf{I}_{-}\mathbf{v}_{-}\,.$
(31)

Equations (B.3) and (31) show that the (classical) symplectic norm defined in
Eq. (25) coincides exactly with the symplectic norm introduced in Ref. [30].
As the sign of Eq. (B.2) is determined by $ds^{2}$, the rotational sense of a
ringdown is directly linked with the sign of the symplectic norm of the
attractor. Since the sign of the symplectic norm of an attractor is positive
(negative) for a minimum (maximum) of the underlying potential [29, 30], the
rotational sense of the ringdown can be used to differ between maxima and
minima in the rotating frame potential.

###  B.4 Spectral Response

Having established a link between the rotational sense of ringdowns and the
symplectic norm, we derive the ‘spectral response’ of fluctuations around an
attractor. This spectral response also allows us to extract the symplectic
norm from ringdown measurements.

In analogy to driven-dissipative quantum systems [29, 30], we first introduce
the correlator function

$\displaystyle G^{c}(t^{\prime}-t)$
$\displaystyle=\Theta(t^{\prime}-t)\langle\tilde{v}(t)\tilde{u}(t^{\prime})-\tilde{u}(t)\tilde{v}(t^{\prime})\rangle_{t}$
(32)
$\displaystyle=\Theta(t^{\prime}-t)\frac{2ds^{2}|\zeta|^{2}}{\Gamma}e^{-\frac{\Gamma}{2}(t^{\prime}-t)}\sin(\omega_{\mathrm{lin}}(t^{\prime}-t)),$

where ‘cc.’ denotes the complex conjugate, and $\Theta(t-t^{\prime})$ is the
Heavyside step-function. The function $G^{c}$ correlates the displaced
quadratures $\tilde{u}(t)$ and $\tilde{v}(t^{\prime})$ at different times to
capture the rotational sense of a ringdown near an attractor. This correlator
is reminiscent of Eq. (25), and is analogous to a retarded Green’s function in
driven-dissipative quantum systems, correlating at different times the
creation and annihilation of bosons on top of a stable solution. Note that the
second equality in Eq. (32) holds in the limit of the linearization procedure
presented earlier.

The sign of the symplectic norm can now be read out directly by inspecting the
peaks and dips of the spectral function [29, 30, 34]. This spectral response
is related to the imaginary part of the Fourier transform of the Green’s
function $G^{c}(\omega)$ via

$\displaystyle\mathcal{A}(\omega)$
$\displaystyle=-2\mathrm{Im}[G^{c}(\omega)].$ (33)

Explicitly computing Eq. (33) using Eq. (32) yields the spectral response of
our system:

$\displaystyle\mathcal{A}(\omega)=\frac{|\zeta|^{2}ds^{2}}{2}\left[\frac{1}{(\omega-\omega_{\mathrm{lin}})^{2}+\frac{\Gamma^{2}}{4}}-\frac{1}{(\omega+\omega_{\mathrm{lin}})^{2}+\frac{\Gamma^{2}}{4}}\right]\,.$
(34)

A peak (dip) of $\mathcal{A}$ at positive frequencies therefore corresponds to
a positive (negative) symplectic norm and hence to a minimum (maximum) of the
rotating frame potential $H_{\mathrm{rot}}$ at
$(u_{\mathrm{f}},v_{\mathrm{f}})$.

##  Appendix C Details of the Hamiltonian Reconstruction

In this section, we provide further details on the Hamiltonian reconstruction
process. We start the reconstruction by using a single ringdown measurement
with the quadratures $u_{j}(t_{j})$ and $v_{j}(t_{j})$ measured at discrete
times $t_{j}$. Assuming that our device is linearly damped, and using the
extracted damping $\Gamma$, we calculate the change in the Hamiltonian $\Delta
H_{\mathrm{rot}}(u_{j},v_{j})$ at each coordinates along the ringdown path by
numerically integrating Eq. (5).

We then find the attractor $k$, e.g. $k={1,2,3}$ for the case of three
attractors, in which the resonator rings down into by looking at the final
value $(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$. We fix the Hamiltonian
$H_{\mathrm{rot}}(u_{\mathrm{f}},v_{\mathrm{f}})=C_{k}$ at the final
coordinate $(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$ to an arbitrary value
$C_{k}$. To do so, we add a constant offset to the Hamiltonian change $\Delta
H_{\mathrm{rot}}$ calculated for the ringdown. This will insure that the
Hamiltonian is singled-valued at each attractor, i.e.,
$H_{\mathrm{rot}}(u_{\mathrm{f}},v_{\mathrm{f}})$ is the same for each
ringdown ending up in the same attractor. At this stage, we simply set
$C_{k}=0$.

As a second step, we repeat the integration procedure for all ringdown
measurements and check if the final value of each ringdown
$(u_{\mathrm{f}},v_{\mathrm{f}})_{k}$ is the same as the first ringdown, or if
it ends up at another attractor. If it is another attractor, we label the new
attractor as $k+1$, and add an arbitrary value $C_{k+1}$ to the trace, which
we will later adjust. We repeat this procedure for all measured ringdowns,
allowing to find the Hamiltonian change for all the ringdown paths, as shown
in Fig. 7 for an example with three attractors.

Figure 7: Hamiltonian reconstruction before adjusting the offset $C_{k}$ of
each attractor. After using Eq. (5) to recover the relative Hamiltonian
potential along each ringdown trajectory and setting the end point of each
ringdown to $C_{k}=0$, we obtain a discontinuous Hamiltonian. The clear
discontinuity happens at the separatrix of the system since the offset between
the different attractors $C_{k}-C_{k+1}$ has yet to be found.

We generally find discontinuities in the Hamiltonian reconstruction at the
separatrices between different attractors, see Fig. 7. Such discontinuities
are unphysical and stem from our choice to set all offsets to $C_{k}=0$. As a
last step of the reconstruction, we thus need to find the relative height
$C_{k}-C_{k+1}$ between the attractors such that the discontinuities
disappear. In practice, we fix one of the constants (e.g., $C_{1}=0$), and
vary the constant offsets $C_{k}$ of the other attractors ($C_{2}$ and $C_{3}$
in Fig. 7) until the Hamiltonian is smooth (i.e., such that the second
differentials of $H_{\mathrm{rot}}$ with respect to $u$ and $v$ are
minimized). We furthermore set the minimum of $H_{\mathrm{rot}}=0$, yielding
Fig. 3(a)(iii). For our demonstration examples, this optimization was done
manually.

## References

  * Cuairan  _et al._ [2022] M. T. Cuairan, J. Gieseler, N. Meyer, and R. Quidant, “Precision calibration of the duffing oscillator with phase control,” Phys. Rev. Lett. 128, 213601 (2022).
  * Bairey  _et al._ [2019] E. Bairey, I. Arad, and N. H. Lindner, “Learning a local hamiltonian from local measurements,” Phys. Rev. Lett. 122, 020504 (2019).
  * Qi and Ranard [2019] X.-L. Qi and D. Ranard, “Determining a local hamiltonian from a single eigenstate,” Quantum 3, 159 (2019).
  * Lin  _et al._ [2014] Z. Lin, K. Inomata, K. Koshino, W. D. Oliver, Y. Nakamura, J. S. Tsai, and T. Yamamoto, “Josephson parametric phase-locked oscillator and its application to dispersive readout of superconducting qubits,” Nature Commu. 5, 4480 (2014).
  * Siva  _et al._ [2023] K. Siva, G. Koolstra, J. Steinmetz, W. P. Livingston, D. Das, L. Chen, J. Kreikebaum, N. Stevenson, C. Jünger, D. Santiago, I. Siddiqi, and A. Jordan, “Time-dependent hamiltonian reconstruction using continuous weak measurements,” PRX Quantum 4, 040324 (2023).
  * Florin  _et al._ [1998] E. Florin, A. Pralle, E. Stelzer, and J. Hörber, “Photonic force microscope calibration by thermal noise analysis,” Appl. Phys. A 66, S75 (1998).
  * McCann  _et al._ [1999] L. I. McCann, M. Dykman, and B. Golding, “Thermally activated transitions in a bistable three-dimensional optical trap,” Nature 402, 785 (1999).
  * Chan and Stambaugh [2007] H. B. Chan and C. Stambaugh, “Activation barrier scaling and crossover for noise-induced switching in micromechanical parametric oscillators,” Phys. Rev. Lett. 99, 060601 (2007).
  * Jeney  _et al._ [2004] S. Jeney, E. H. Stelzer, H. Grubmüller, and E.-L. Florin, “Mechanical properties of single motor molecules studied by three-dimensional thermal force probing in optical tweezers,” ChemPhysChem 5, 1150 (2004).
  * Tanaka  _et al._ [2013] Y. Tanaka, S. Kaneda, and K. Sasaki, “Nanostructured potential of optical trapping using a plasmonic nanoblock pair,” Nano Lett. 13, 2146 (2013).
  * Rondin  _et al._ [2017] L. Rondin, J. Gieseler, F. Ricci, R. Quidant, C. Dellago, and L. Novotny, “Direct measurement of kramers turnover with a levitated nanoparticle,” Nat. Nanotechnol. 12, 1130 (2017).
  * Zijlstra  _et al._ [2020] N. Zijlstra, D. Nettels, R. Satija, D. E. Makarov, and B. Schuler, “Transition path dynamics of a dielectric particle in a bistable optical trap,” Phys. Rev. Lett. 125, 146001 (2020).
  * Militaru  _et al._ [2021] A. Militaru, M. Innerbichler, M. Frimmer, F. Tebbenjohanns, L. Novotny, and C. Dellago, “Escape dynamics of active particles in multistable potentials,” Nat. Comm. 12, 2446 (2021).
  * Wu  _et al._ [2009] P. Wu, R. Huang, C. Tischer, A. Jonas, and E.-L. Florin, “Direct measurement of the nonconservative force field generated by optical tweezers,” Phys. Rev. Lett. 103, 108101 (2009).
  * Pérez García  _et al._ [2018] L. Pérez García, J. Donlucas Pérez, G. Volpe, A. V. Arzola, and G. Volpe, “High-performance reconstruction of microscopic force fields from brownian trajectories,” Nat. Comm. 9, 5166 (2018).
  * Heugel  _et al._ [2023] T. L. Heugel, R. Chitra, A. Eichler, and O. Zilberberg, “Proliferation of unstable states and their impact on stochastic out-of-equilibrium dynamics,” arXiv:2307.13718 (2023).
  * Álvarez  _et al._ [2023] P. Álvarez, D. Pittilini, F. Miserocchi, S. Raamamurthy, G. Margiani, O. Ameye, J. del Pino, O. Zilberberg, and A. Eichler, “A biased ising model using two coupled kerr parametric oscillators with external force,” arXiv:2307.13676 (2023).
  * Swingle [2018] B. Swingle, “Unscrambling the physics of out-of-time-order correlators,” Nat. Phys, 14, 988 (2018).
  * Blocher  _et al._ [2022] P. D. Blocher, S. Asaad, V. Mourik, M. A. I. Johnson, A. Morello, and K. Mølmer, “Measuring out-of-time-ordered correlation functions without reversing time evolution,” Phys. Rev. A 106, 042429 (2022).
  * Bachtold  _et al._ [2022] A. Bachtold, J. Moser, and M. I. Dykman, “Mesoscopic physics of nanomechanical systems,” Rev. Mod. Phys. 94, 045005 (2022).
  * Antoni  _et al._ [2012] T. Antoni, K. Makles, R. Braive, T. Briant, P.-F. Cohadon, I. Sagnes, I. Robert-Philip, and A. Heidmann, “Nonlinear mechanics with suspended nanomembranes,” EPL 100, 68005 (2012).
  * Polunin  _et al._ [2016] P. M. Polunin, Y. Yang, M. I. Dykman, T. W. Kenny, and S. W. Shaw, “Characterization of mems resonator nonlinearities using the ringdown response,” J. Microelectromechanical Syst. 25, 297 (2016).
  * Catalini  _et al._ [2021] L. Catalini, M. Rossi, E. C. Langman, and A. Schliesser, “Modeling and observation of nonlinear damping in dissipation-diluted nanomechanical resonators,” Phys. Rev. Lett. 126, 174101 (2021).
  * Morales  _et al._ [2018] A. Morales, P. Zupancic, J. Léonard, T. Esslinger, and T. Donner, “Coupling two order parameters in a quantum gas,” Nat, Mater. 17, 686 (2018).
  * Kozinsky  _et al._ [2007] I. Kozinsky, H. C. Postma, O. Kogan, A. Husain, and M. L. Roukes, “Basins of attraction of a nonlinear nanomechanical resonator,” Phys. Rev. Lett. 99, 207201 (2007).
  * Yan  _et al._ [2022] Y. Yan, X. Dong, L. Huang, K. Moskovtsev, and H. Chan, “Energy transfer into period-tripled states in coupled electromechanical modes at internal resonance,” Phys. Rev. X 12, 031003 (2022).
  * Woo and Landauer [1971] J. Woo and R. Landauer, “Fluctuations in a parametrically excited subharmonic oscillator,” IEEE J. Quantum Electron 7, 435 (1971).
  * Dykman  _et al._ [1998] M. I. Dykman, C. M. Maloney, V. N. Smelyanskiy, and M. Silverstein, “Fluctuational phase-flip transitions in parametrically driven oscillators,” Phys. Rev. E 57, 5202 (1998).
  * Soriente  _et al._ [2020] M. Soriente, R. Chitra, and O. Zilberberg, “Distinguishing phases using the dynamical response of driven-dissipative light-matter systems,” Phys. Rev. A 101, 023823 (2020).
  * Soriente  _et al._ [2021] M. Soriente, T. L. Heugel, K. Omiya, R. Chitra, and O. Zilberberg, “Distinctive class of dissipation-induced phase transitions and their universal characteristics,” Phys. Rev. Res. 3, 023100 (2021).
  * Ferri  _et al._ [2021] F. Ferri, R. Rosa-Medina, F. Finger, N. Dogra, M. Soriente, O. Zilberberg, T. Donner, and T. Esslinger, “Emerging dissipative phases in a superradiant quantum gas with tunable decay,” Phys. Rev. X 11, 041046 (2021).
  * Roy  _et al._ [2021] A. Roy, S. Jahani, C. Langrock, M. Fejer, and A. Marandi, “Spectral phase transitions in optical parametric oscillators,” Nat. Comm. 12, 835 (2021).
  * Zhang  _et al._ [2021] G.-Q. Zhang, Z. Chen, W. Xiong, C.-H. Lam, and J. Q. You, “Parity-symmetry-breaking quantum phase transition via parametric drive in a cavity magnonic system,” Phys. Rev. B 104, 064423 (2021).
  * Scarlatella  _et al._ [2019] O. Scarlatella, A. A. Clerk, and M. Schiro, “Spectral functions and negative density of states of a driven-dissipative nonlinear quantum resonator,” New J. Phys. 21, 043040 (2019).
  * Agarwal  _et al._ [2008] M. Agarwal, S. A. Chandorkar, H. Mehta, R. N. Candler, B. Kim, M. A. Hopcroft, R. Melamud, C. M. Jha, G. Bahl, G. Yama, T. W. Kenny, and B. Murmann, “A study of electrostatic force nonlinearities in resonant microstructures,” Appl. Phys. Lett. 92, 104106 (2008).
  * Schmid  _et al._ [2016] S. Schmid, L. G. Villanueva, and M. L. Roukes,  _Fundamentals of nanomechanical resonators_ , Vol. 49 (Springer, 2016).
  * Miller  _et al._ [2019] J. M. Miller, D. D. Shin, H.-K. Kwon, S. W. Shaw, and T. W. Kenny, “Phase control of self-excited parametric resonators,” Phys. Rev. Applied 12, 044053 (2019).
  * Nosan  _et al._ [2019] Ž. Nosan, P. Märki, N. Hauff, C. Knaut, and A. Eichler, “Gate-controlled phase switching in a parametron,” Phys. Rev. E 99, 062205 (2019).
  * Marthaler and Dykman [2006] M. Marthaler and M. I. Dykman, “Switching via quantum activation: A parametrically modulated oscillator,” Phys. Rev. A 73, 042108 (2006).
  * Ryvkine and Dykman [2006] D. Ryvkine and M. I. Dykman, “Resonant symmetry lifting in a parametrically modulated oscillator,” Phys. Rev. E 74, 061118 (2006).
  * Mahboob and Yamaguchi [2008] I. Mahboob and H. Yamaguchi, “Bit storage and bit flip operations in an electromechanical oscillator,” Nat. Nanotechn. 3, 275 (2008).
  * Wilson  _et al._ [2010] C. M. Wilson, T. Duty, M. Sandberg, F. Persson, V. Shumeiko, and P. Delsing, “Photon generation in an electromagnetic cavity with a time-dependent boundary,” Phys. Rev. Lett. 105, 233907 (2010).
  * Leuch  _et al._ [2016] A. Leuch, L. Papariello, O. Zilberberg, C. L. Degen, R. Chitra, and A. Eichler, “Parametric symmetry breaking in a nonlinear resonator,” Phys. Rev. Lett. 117, 214101 (2016).
  * Gieseler  _et al._ [2012] J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny, “Subkelvin parametric feedback cooling of a laser-trapped nanoparticle,” Phys. Rev. Lett. 109, 103603 (2012).
  * Grimm  _et al._ [2020] A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar, and M. H. Devoret, “Stabilization and operation of a kerr-cat qubit,” Nature 584, 205 (2020).
  * Wang  _et al._ [2019] Z. Wang, M. Pechal, E. A. Wollack, P. Arrangoiz-Arriola, M. Gao, N. R. Lee, and A. H. Safavi-Naeini, “Quantum dynamics of a few-photon parametric oscillator,” Phys. Rev. X 9, 021049 (2019).
  * Yamaji  _et al._ [2022] T. Yamaji, S. Kagami, A. Yamaguchi, T. Satoh, K. Koshino, H. Goto, Z. R. Lin, Y. Nakamura, and T. Yamamoto, “Spectroscopic observation of the crossover from a classical duffing oscillator to a kerr parametric oscillator,” Phys. Rev. A 105, 023519 (2022).
  * Eichler and Zilberberg [2023] A. Eichler and O. Zilberberg,  _Classical and Quantum Parametric Phenomena_ (Oxford University Press, 2023).
  * Note [1] Note that $H_{\mathrm{rot}}$ may not have dimensions of energy in this notation because $u$ and $v$ have the same units. The quasi-Hamiltonian can be expressed in units of energy by changing $H_{\mathrm{rot}}\rightarrow H_{\mathrm{rot}}/m\omega_{\mathrm{r}}$ in Eqs. (3) and (III).
  * Lifshitz [2009] M. C. Lifshitz, R. Cross, “Nonlinear dynamics of nanomechanical and micromechanical resonators,” in  _Reviews of Nonlinear Dynamics and Complexity_ (Wiley, 2009) pp. 1–52.
  * Dykman [2012] M. Dykman,  _Fluctuating Nonlinear Oscillators_ (Oxford University Press, 2012).
  * Heo  _et al._ [2010] M.-S. Heo, Y. Kim, K. Kim, G. Moon, J. Lee, H.-R. Noh, M. I. Dykman, and W. Jhe, “Ideal mean-field transition in a modulated cold atom system,” Phys. Rev. E 82, 031134 (2010).
  * Karabalin  _et al._ [2010] R. Karabalin, S. Masmanidis, and M. Roukes, “Efficient parametric amplification in high and very high frequency piezoelectric nanoelectromechanical systems,” Appl. Phys. Lett. 97, 183101 (2010).
  * Soriente  _et al._ [2018] M. Soriente, T. Donner, R. Chitra, and O. Zilberberg, “Dissipation-induced anomalous multicritical phenomena,” Phys. Rev. Lett. 120, 183603 (2018).
  * Dirac [1932] P. A. M. Dirac, “Relativistic quantum mechanics,” Proc. R. Soc. Lond. A 136, 453 (1932).
  * Ricardo [2021] H. J. Ricardo,  _A Modern Introduction to Differential Equations (Third Edition)_, edited by H. J. Ricardo (Academic Press, 2021) pp. 361–420.
  * Košata  _et al._ [2022] J. Košata, J. del Pino, T. L. Heugel, and O. Zilberberg, “HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance,” SciPost Phys. Codebases , 6 (2022).
  * Borovik  _et al._ [2024] V. Borovik, P. Breiding, J. del Pino, M. Michałek, and O. Zilberberg, “Khovanskii bases for semimixed systems of polynomial equations – approximating stationary nonlinear newtonian dynamics,” J. Math. Pures Appl. 182, 195 (2024).
  * Fan and Jia [2023] J. Fan and S. Jia, “Collective dynamics of the unbalanced three-level dicke model,” Phys. Rev. A 107, 033711 (2023).
  * Xiao [2009] M. Xiao, “Theory of transformation for the diagonalization of quadratic hamiltonians,” arXiv:0908.0787 (2009).