# Limit cycles as stationary states of an extended Harmonic Balance ansatz

Javier del Pino  Institute for Theoretical Physics, ETH Zürich, 8093 Zürich,
Switzerland  Jan Košata  Institute for Theoretical Physics, ETH Zürich, 8093
Zürich, Switzerland  Oded Zilberberg  Department of Physics, University of
Konstanz, 78464 Konstanz, Germany

###### Abstract

A limit cycle is a self-sustained periodic motion appearing in autonomous
ordinary differential equations. As the period of the limit cycle is a-priori
unknown, it is challenging to find them as stationary states of a rotating
ansatz. Correspondingly, their study commonly relies on brute-force time-
evolution or on circumstantial evidence such as instabilities of fixed points.
Alas, such approaches are unable to account for the coexistence of multiple
solutions, as they rely on specific initial conditions. Here, we develop a
multifrequency rotating ansatz with which we find limit cycles as stationary
states. We demonstrate our approach and its performance in the simplest case
of the Van der Pol oscillator. Moving beyond the simplest example, we show
that our method can capture the coexistence of all fixed-point attractors and
limit cycles in a modified nonlinear Van der Pol oscillator. Our results
facilitate the systematic mapping of out-of-equilibrium phase diagrams, with
implications across all fields of natural science.

Limit cycles (LCs) are periodic solutions of nonlinear ordinary differential
equations (ODEs) that persist indefinitely in time. Commonly, the oscillation
period is a-priori unknown, with oscillations reflecting the autonomous
balance between the various forces acting on the system. Limit cycles are
ubiquitous in all scientific fields, including physics, biology, and
engineering [1], with broad applications, such as developing optical frequency
combs and sensors [2, 3, 4], understanding biological oscillators [5, 6, 7],
designing control systems [8], as well as realizing neuromorphic computing
[9]. Limit cycles have recently gained significant attention in the realm of
out-of-equilibrium physics with a variety of realizations in light-matter
ensembles [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], parametric oscillators
[21, 22, 23, 24], electromechanical systems [25, 26, 27, 28], superconducting
circuits [29], and optomechanics [30, 31, 32, 33, 34, 35].

The concept of LCs is defined differently in different fields, leading to
ambiguity in the literature [36]. Most definitions require the LC to be an
isolated, self-sustained periodic motion, where the system is stationary,
where its Fourier transform amplitudes remain time-invariant and are
asymptotically stable. While some definitions allow for more complex solutions
such as quasiperiodic or chaotic motions, the key feature of a LC is the
absence of a preferred time origin and phase locking to external forces. This
“resilience" against external perturbations or noise is crucial for the role
LCs play in synchronization [37], and in instabilities in fluid dynamics [38,
39].

The study of LCs poses enduring challenges, due to their initially unknown
oscillation period and the nonlinear nature of the problem. Various methods
have been proposed for their study, including brute-force ODE integration [40]
and bifurcation theory of unstable fixed points [41, 42]. While useful, these
methods crucially depend on the choice of initial conditions or are applicable
only for specific types of LCs. As such, they often fail to capture the
coexistence of multiple LCs in the system. Moreover, to date, the upper bound
for the number of LCs in a polynomial vector field remains unknown [43].

The Harmonic Balance Method (HBM) [44, 45, 46, 47, 48] stands as a powerful
technique for solving stationary (static and periodic) solutions of ODEs. The
method relies on an ansatz assuming that the solution takes the form of a
multi-frequency Fourier series. Then, a nonlinear optimization algorithm is
used to determine the coefficients of the series that satisfy the boundary
conditions and equations of motion. The HBM connects to the Krylov-Bogoliubov
averaging method for solving nonlinear ODEs with harmonic dependence [49, 50],
effectively describing various effects in nanomechanics [51, 52, 53], optics
[54], optomechanics [55] and nuclear magnetic resonance [56].

Figure 1: Gain versus loss and Hopf bifurcations. Stream plots illustrating
the behavior of a linear resonator with (a) positive loss, and (b) positive
gain around the origin of parameter space (marked with black dot), cf. Eq. (1)
with $\alpha_{3,5}=0$ and $(\gamma-\mu)>0$ or $(\gamma-\mu)<0$, respectively.
(c) A plot of the Jacobian eigenvalues around the origin $(0,0)$ as a function
of linear damping for a Van der Pol oscillator [cf. Eq. (1)]. Their imaginary
part (orange, dashed-dotted lines) marks the oscillation frequency around the
origin (excitation). Their real part (blue lines) denotes the inverse lifetime
of the excitations (when negative) or their rate of expansion (when positive).
From right to left, we observe transitions from overdamped to underdamped
stable excitation (see the exceptional point at $\gamma-\mu=2\omega_{0}$, star
marker), to unbound (unstable) spiral (Hopf bifurcation at 0, cross marker),
and to an unbound non spiraling expansion (see the exceptional point at
$\gamma-\mu=-2\omega_{0}$, star marker). Vertical dashed lines mark cases (a)
and (b).

In this work, we introduce a variant of the HBM that can efficiently find LCs
as stationary states. Specifically, we employ the HBM ansatz, and through
gauge-fixing allow for a variational search of the LC frequencies using a
powerful Homotopy continuation method [57, 58, 59]. A substantial advantage of
our approach is the ability to find the complete set of both LCs and fixed
point stationary states under a generic ansatz, even when they coexist. We
first illustrate the merit of our approach on the prototypical Van der Pol
(VdP) oscillator, which is the simplest case study for LCs. We then generalize
the model to include nonlinearities, where LCs and fixed point attractors can
coexist. Crucially, our approach can be applied to efficiently explore the
phase diagrams of nonlinear, coupled, driven-dissipative systems [22, 60, 61,
62], including those commonly encountered in nanomechanics, cold atoms, and
nonlinear photonics.

We consider a modified VdP oscillator that includes additional polynomial
potential terms in the equation. The oscillator displacement $x(t)$ evolves
according to Newton’s equation $\ddot{x}+F(x,\dot{x})=0$ with

$F(x,\dot{x})=\omega_{0}^{2}x+\gamma\dot{x}-\mu(1-x^{2})\dot{x}+\alpha_{3}x^{3}+\alpha_{5}x^{5}\,,$
(1)

where the natural angular frequency is $\omega_{0},\gamma$ controls the linear
damping, $\mu>0$ is the standard VdP parameter that controls the linear gain
channel alongside nonlinear damping, and $\alpha_{3,5}$ are cubic and quintic
nonlinearities, respectively. The limit of $\gamma,\alpha_{3,5}=0$ in Eq. (1)
recovers the standard VdP equation, where the gain channel, $-\mu\dot{x}$,
drives the oscillator to higher amplitudes, and the nonlinear damping, $\mu
x^{2}\dot{x}$, saturates the motion into a limit cycle (cf. the sign change of
the overall term based on the value of $|x|$). In an opposite limit where only
$\mu=0$, we have a nonlinear oscillator with quartic and sextic potential
terms. This system can have multiple energy minima, which will be studied in
the final part of the manuscript.

Our model has two distinct regimes related to the sign of $(\gamma-\mu)$. A
positive sign signals a damped system that will relax into one of the
potential minima depending on the initial boundary conditions, see Fig. 1(a).
In the opposite case (negative sign), gain overcomes damping and the system is
incoherently driven to higher amplitudes, see Fig. 1(b). Generally, for each
point in configuration space, ${\mathbf{x}}_{0}=(x_{0},\dot{x}_{0})$, we can
evaluate the Jacobian matrix $\mathcal{J}({\mathbf{x}}_{0})$ around that point
to see where excitations flow to in the linear vicinity of ${\mathbf{x}}_{0}$.
In particular, by diagonalizing the Jacobian at the origin,
${\mathbf{x}}_{0}=(0,0)$, we can discern between the two regimes above: the
real part of the Jacobian eigenvalues
$\epsilon_{0}=\frac{1}{2}\left(-\left(\gamma-\mu\right)\pm\sqrt{(\gamma-\mu)^{2}-4\omega_{0}^{2}}\right)$
corresponds to whether the origin is a sink or a source (cf. discussion on the
sign of $\gamma-\mu$ above); the imaginary part parametrizes the oscillation
frequency for small excitations near the origin, see Fig. 1(c). These small
excitations are stable and overdamped when $\gamma-\mu>2\omega_{0}$, and
stable and underdamped when $0<\gamma-\mu<2\omega_{0}$. The under- to over-
damped transition is marked by an exceptional point [63, 64]. The source-to-
sink transition occurs via a Hopf bifurcation, where the excitations lose
stability while maintaining their oscillation frequency. At higher gain, no
notion of unstable excitations remains as the oscillation frequency is
removed.

Figure 2: The van der Pol limit cycle and the harmonic balance method. (a)
Stream plots depicting the dynamics in the VdP oscillator [cf. Eq. (1) with
$\alpha_{3,5},\gamma=0$] for $\mu=0.1\omega_{0}$, $\mu=0.5\omega_{0}$, and
$\mu=\omega_{0}$, from left-to-right, respectively. The solid (red) curve
marks the LC obtained from direct time propagation, while the dashed (black)
curve shows the result from the extended HBM with $M=5$ harmonics [cf. Eqs.
(4)]. (b) Fourier transforms of the LC time evolutions in (a) unveil the
emergence of a frequency comb. Red dots show normalized squared amplitudes at
odd frequencies $(2k+1)\omega$ $k\in\mathbb{N}$, captured by the e-HBM. We use
windowed Fourier transforms with a Hamming function to mitigate spectral
leakage caused by finite integration times. (c) Overlap between e-HBM and
exact time evolution result as a function of $\mu$ for $M=5$ [cf. Eq. (5)].
Dots show a linear fit. (inset) Same as (b) as a function of $\mu$ marking the
change in the comb frequencies. e-HBM frequencies are shown as dashed lines,
overlaid on the positive part of the Fourier spectrum.

The local (sink-to-source) instability around the origin does not provide
sufficient information regarding whether a LC forms once the Hopf bifurcation
occurs. Indeed, without any nonlinearity in the model, the excitations will
expand under the gain to infinity. However, in the standard VdP oscillator
($\alpha_{3,5},\gamma=0$), the nonlinear damping prevents the system from
spiraling to $x,\dot{x}\rightarrow\infty$, by stabilizing a LC that encircles
the origin at a distance, see Figs. 2(a). Interestingly, the geometry of the
LC is determined by the value of $\mu$. For smaller values, the motion
exhibits a circular orbit in configuration space, which gradually morphs into
a more complex pattern as $\mu$ increases. The periodic character of the LC
makes it suitable for a Fourier decomposition, see Figs. 2(b). We observe that
the oscillation frequency is not simply that of the bare resonator
($\omega_{0}$). Specifically, the deviation from a single harmonic behavior (a
circular motion in configuration space) is captured by the Fourier amplitudes
$|\tilde{x}(\omega)|^{2}$, which reveal the formation of a frequency comb with
harmonics $m\omega$. Note that the LC, in our case, has a half-wave symmetry
in the time domain, such that only odd harmonics contribute.

Due to the intricate motion of LCs, with unknown frequencies, it is
challenging to describe them as stationary states of frequency expansion,
i.e., a rotating ansatz. Instead, LCs are commonly obtained by brute-force
numerical integration, or by seeding them out of a Hopf bifurcation, cf. Fig.
1. However, both methods have their limitations. The former may not capture
the LC if the initial condition is not within its basin of attraction. The
latter approach is not universal in capturing LCs: (i) a parameter deformation
of the system that leads to a sink-to-source transition –– a necessary
condition for a Hopf bifurcation –– is required; in the standard VdP with
$\mu>0$, this transition is omitted. In addition, (ii) LCs can arise from
various types of bifurcations, and not just Hopf, e.g., orbits connecting
saddle equilibria, and even from infinity [50, 65]. Consequently, in the
following, we devise a variant of HBM [47] that captures LCs as stationary
states of a varying rotating frame, i.e., at a time-evolving frequency to be
discovered.

In the standard HBM ansatz, the motion

$x(t)=\sum_{l=1}^{M}u_{l}(t)\cos(\omega_{l}t)+v_{l}(t)\sin(\omega_{l}t)\,,$
(2)

is split into $M$ fast harmonics with frequencies
$\boldsymbol{\omega}=(\omega_{1},\ldots,\omega_{M})$ that have slowly evolving
amplitudes ${\mathbf{u}}=(u_{1}(t),\ldots,u_{M}(t))$ and
${\mathbf{v}}=(v_{1}(t),\ldots,v_{M}(t))$. The ansatz (2) approximates the
time evolution with its leading frequency components. Finding the amplitudes
${\mathbf{u}}$ requires inserting it into Eq. (1), which for $\alpha_{3,5}=0$
yields

$\displaystyle\sum_{l}\cos\left(\omega_{l}t\right)\left[\ddot{u}_{l}+2\omega_{l}\dot{v}_{l}-\mu\dot{u}_{l}+\left(\omega_{0}^{2}-\omega_{l}^{2}\right)u_{l}-\mu\omega_{l}v_{l}\right]+\sum_{l}\sin\left(\omega_{l}t\right)\left[\ddot{v}_{l}-2\omega_{l}\dot{u}_{l}-\mu\dot{v}_{l}+\mu\omega_{l}u_{l}+\left(\omega_{0}^{2}-\omega_{l}^{2}\right)v_{l}\right]=$
(3)
$\displaystyle\sum_{m,n,l}\mu\left[u_{m}\cos\left(\omega_{m}t\right)+v_{m}\sin\left(\omega_{m}t\right)\right]\left[u_{n}\cos\left(\omega_{n}t\right)+v_{n}\sin\left(\omega_{n}t\right)\right]\cdot\hskip
5.69054pt\left[\sin\left(\omega_{l}t\right)\left(\omega_{l}u_{l}-\dot{v}_{l}\right)-\cos\left(\omega_{l}t\right)\left(\dot{u}_{l}+\omega_{l}v_{l}\right)\right]\,,$

We search for stationary motion implying fixed harmonic amplitudes, i.e.,
$\ddot{\mathbf{u}}=\ddot{\mathbf{v}}=\dot{\mathbf{u}}=\dot{\mathbf{v}}=0$.
Furthermore, the premise of the HBM is to “balance” the harmonics at both
sides of Eq. (3), i.e., the prefactors of each harmonic satisfy the equation
independently. Note that nonlinear terms in Eq. (1) lead to frequency mixing
among different harmonics through ’multi-photon’ resonance conditions,
$\omega_{l}=\omega_{l}\pm\omega_{m}\pm\omega_{n}$. These conditions lead to
third-degree nonlinear coupling among the prefactors of various harmonics. As
a result, we obtain from Eq. (1) a system of $2M$ coupled polynomial equations
$\mathcal{F}({\mathbf{u}},{\mathbf{v}})$. When the stationary state
frequencies $\omega_{l}$ are known, as in stationary motion fixed by an
external driving, the HBM simplifies the problem to finding the roots
$\mathcal{F}({\mathbf{u}},{\mathbf{v}})=0$ [66].

Figure 3: Coexistence of limit cycles and fixed points. (a) Phase diagram of
the model (1) in the $\alpha_{3},\mu$ parameter space, showing the count of
stable stationary states, including LCs and fixed points. Dashed lines mark
additional boundaries where the number of unstable stationary states changes.
Panels (b)-(g) illustrate distinct behaviors in each region through stream
plots [cf. Fig. 1]]. Stable (unstable) fixed points derived from the HBM are
depicted by dots (crosses). The LC found by the e-HBM is shown by a red dashed
line. Regions (b) and (d) exist below the Hopf bifurcation, i.e.,
$\mu<\gamma$, and only sustain fixed points. When $\mu>\gamma$, regions (c)
and (e) host a LC around the origin, that coexists with stable fixed points in
(e). The transition between (b) and (d) [respectively (c) and (e)] is
associated with the appearance of two additional fixed points. The LC
destabilizes in (f) and vanishes in (g) by merging with the unstable fixed
points, and leaves behind solely fixed points. In all panels, $\omega_{0}=1$,
$\mu=0.1$, $\gamma=0.03$, $\alpha_{5}-0.01$ and $M=2$.

Limit cycles evolve with unknown self-oscillation frequencies. We, therefore,
modify the standard HBM approach and introduce a slow time variation to the
frequencies, $\omega_{l}\mapsto\omega_{l}(t)$, on top of the slowly evolving
amplitudes $u_{l}(t),v_{l}(t)$ in the ansatz (2). In our extended HBM (e-HBM)
approach, the variational prefactors evolve similarly to Eqs. (3), with
additional terms arising from time derivatives of $\omega_{l}(t)$. Now, fixing
the harmonic amplitudes in the generalized stationary state, where
$\ddot{{\mathbf{u}}}=\dot{{\mathbf{u}}}=\ddot{{\mathbf{v}}}=\dot{{\mathbf{v}}}=\ddot{\boldsymbol{\omega}}=\dot{\boldsymbol{\omega}}=0$,
results in a system of $2M$ coupled polynomials
$\tilde{\mathcal{F}}({\mathbf{u}},{\mathbf{v}},\boldsymbol{\omega})=0$ in the
$3M$ variables ${\mathbf{u}},{\mathbf{v}},\boldsymbol{\omega}$. The harmonic
balance equations are thus insufficient to determine unequivocally the LCs.
Fortunately, we can deduce the missing $M$ conditions from the spontaneous,
self-oscillating nature of the LC: With no external drive to lock onto, the
LC’s phases are free from constraints. This freedom leaves the complex
harmonic amplitude vector ${\mathbf{a}}=({\mathbf{u}}+i{\mathbf{v}})/2$ with a
free $U(1)$ gauge symmetry under global transformations $a_{i}\mapsto
e^{i\varphi_{i}}a_{i}$, where $\varphi_{i}\in(0,2\pi)$ [66]. As the stationary
solution of the LC does not rely on the gauge, we can fix the $M$ gauges by
redefining the time origin in Eq. (2). One possible gauge choice yields, e.g.,
the following system of $3M$ equations

$\\{\mathcal{F}({\mathbf{u}},{\mathbf{v}},\boldsymbol{\omega})=0,{\mathbf{u}}=0\\}\,.$
(4)

Note that our treatment simplifies when introducing relationships within
$\boldsymbol{\omega}$, e.g. commensurate frequencies. For example, the LC in
Eq. (2)(b) is described by a frequency comb with a single unknown frequency
$\omega$, i.e., $\omega_{l}=l\omega$ with $l\in\mathbb{N}$. This requires only
fixing the gauge of a single harmonic to attain a solution, e.g., by setting
$u_{1}=0$.

Our e-HBM embeds LCs as fixed points in a high-dimensional phase space,
rotating with variational frequencies $\omega_{l}$. The complexity of solving
the problem at this stage is delineated to the proliferation of roots in the
system, which complicates the use of straightforward root-finding (e.g.,
Newton-Raphson) methods. Inspired by a similar challenge encountered in the
standard HBM for driven motion [47], we solve Eqs. (4) using Homotopy
Continuation [57, 58, 48, 67]. Homotopy Continuation combines root finding
with continuous deformation of an exactly solvable algebraic system into the
target coupled polynomials (4) to find all the roots in a single run [57, 58].
Applying our methodology to the VdP oscillator ($\alpha_{3,5}=\gamma=0$), we
find a single real solution that corresponds to the amplitudes
${\mathbf{u}},{\mathbf{v}}$ and frequencies $\boldsymbol{\omega}$ that well
estimate the LC [Figs. 2(a)]. This same agreement manifests more clearly in
the frequency domain [see Fig. 2(b)], and persists as a function of $\mu$ for
various LC frequencies, see inset of 2(c). To further benchmark the quality of
the estimated solution, we compare the estimated e-HBM stationary amplitudes,
and those obtained from the Fourier transform of the exact time-evolution. To
this end, we evaluate the overlap

$\displaystyle\mathcal{O}={\mathbf{a}}_{\text{HBM}}^{*}\cdot{\mathbf{a}}_{\text{FT}}\,,$
(5)

between the normalized harmonic amplitudes
${\mathbf{a}}_{\text{HBM}}=({\mathbf{u}}+i{\mathbf{v}})/\sqrt{{\mathbf{u}}^{2}+{\mathbf{v}}^{2}}$,
and ${\mathbf{a}}_{\text{FT}}$, respectively. The overlap decreases with
increasing $\mu/\omega_{0}$, approximately as $\sim(\mu/\omega_{0})^{2}$ [cf.
Fig. 2(c)] reaching few-percent deviations for $\mu=2\omega_{0}$. The
relationship between the overlap and the number of harmonics, $M$, is not
strictly monotonic due to the specific orders retained in the HBM truncation
[66].

The combination of HBM and Homotopy Continuation captures all solutions
simultaneously, unlike direct time evolution, which only explores a single
solution from a specific initial condition (a basin of attraction). While the
system with $\alpha_{3,5}=0$ sustains a unique limit cycle, as guaranteed by
Lienard’s theorem [68], when $\alpha_{3,5}\neq 0$ a coexistence of different
stationary states manifests, see Fig. 3(a). Such coexistence arises from the
conservative potential that hosts multiple energy minima, whereas around the
origin the LC physics manifests. To describe this coexistence using e-HBM, we
apply a two-step procedure: (i) we first solve for stationary fixed-points
using a zero-frequency ansatz $x(t)=a_{0}(t)$, and then (ii) solve for the LCs
using the e-HBM ansatz (2). Here, convergence is achieved by incorporating
harmonics at frequencies $0$ and $(2l+1)\omega$ with $l\in\mathbb{N}$. Note
also that the solution for $a_{0}$ recovers also the trivial fixed point in
the VdP oscillator (the origin). The two steps are needed because the e-HBM
ansatz leads to infinitely-degenerate solutions arising around fixed points,
as it corresponds to zero amplitude solutions for arbitrary
$\boldsymbol{\omega}\neq 0$. Importantly, by combining the results from the
two steps, we fully estimate the expected outcome. Crucially, we can now use
this combined e-HBM to capture the stationary solutions of the system
(including LCs) for different parameters, e.g., as a function of the cubic
potential and nonlinear gain, and obtain the phase diagram of the system, see
Fig. 3(a) highlighting labeled distinct regions. Remarkably, a stable LC can
coexist with fixed point attractors at region (e). For the purpose of this
work, we focus on the parameter regime where LCs coexist with fixed points,
and do not study the potential mixing between the LC and attractors in the
system. This coexistence exemplifies the loose connection between Hopf
bifurcations and LCs, as no bifurcation appears, see e.g. transition between
regions (e) to (f).

Our method facilitates the investigation of LCs in a breadth of fields and is
readily available as part of HarmonicBalance.jl [48]. Instead of relying on
numerical evolution and experimentation, it finds the LCs as fixed points of
an extended rotating ansatz. As such, our treatment can be readily extended to
non-autonomous systems, particularly those subjected to periodic forcing,
where dynamics can undergo spontaneous shifts from periodic to quasi-periodic
behavior [21]. By employing this approach, we expect to uncover the
coexistence of solutions with stable phases, and also identify the presence of
multiple limit cycles within the same system [69, 70]. Crucially, our approach
is well-suited for the analysis of many-body phases in mean-field problems
[71, 72, 73], and of the conditions for synchronization in classical and
quantum oscillator arrays [74, 75, 76, 77, 78, 79, 80, 81, 82, 83]. Notably,
in all these fields, the ability to describe the LCs as a fixed point provides
the framework for studying stochastic dynamics and activation between LCs.
Last, while the e-HBM serves as a promising approach, its effectiveness is
ultimately constrained by the solution of $2M$ coupled polynomial equations
[84]. Progress towards efficiently solving the equations (4) will allow us to
solve complex systems with more degrees of freedom [85].

Note added: While preparing our manuscript, we encountered a manuscript that
proposes a related approach for LCs in aerodynamics [39], where a related HBM
approach is used followed by a single solution continuation.

We thank A. Eichler, A. Nunnenkamp, and C. W. Wächtler for fruitful
discussions. J.d.P. was supported by the ETH Fellowship program (grant no.
20-2 FEL-66). J.K. was supported by the Swiss National Science Foundation
through grant CRSII5 $177198/1$. O.Z. acknowledges funding through SNSF grant
CRSII5_206008/1, and the Deutsche Forschungsgemeinschaft (DFG) through project
number 449653034, and via SFB1432.

## References

  * Strogatz [2000] S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering_ (Westview Press, 2000).
  * Kim  _et al._ [2020] S. Kim, Y. G. Rubo, T. C. Liew, S. Brodbeck, C. Schneider, S. Höfling, and H. Deng, Emergence of microfrequency comb via limit cycles in dissipatively coupled condensates, Physical Review B 101, 085302 (2020).
  * Hamerly and Mabuchi [2015] R. Hamerly and H. Mabuchi, Optical devices based on limit cycles and amplification in semiconductor optical cavities, Phys. Rev. Appl. 4, 024016 (2015).
  * Piergentili  _et al._ [2021] P. Piergentili, W. Li, R. Natali, D. Vitali, and G. Di Giuseppe, Absolute determination of the single-photon optomechanical coupling rate via a hopf bifurcation, Physical Review Applied 15, 034012 (2021).
  * May [1972] R. M. May, Limit cycles in predator-prey communities, Science 177, 900 (1972).
  * Jewett and Kronauer [1998] M. E. Jewett and R. E. Kronauer, Refinement of limit cycle oscillator model of the effects of light on the human circadian pacemaker, Journal of Theoretical Biology 192, 455 (1998).
  * Turchin [2003] P. Turchin,  _Complex Population Dynamics: A Theoretical/Empirical Synthesis (MPB-35)_ (Princeton University Press, 2003).
  * Isidori [2013] A. Isidori,  _Nonlinear control systems II_, Communications and Control Engineering (Springer, 2013).
  * Chalkiadakis and Hizanidis [2022] D. Chalkiadakis and J. Hizanidis, Dynamical properties of neuromorphic josephson junctions, Physical Review E 106, 044206 (2022).
  * Seibold  _et al._ [2019] K. Seibold, R. Rota, and V. Savona, A dissipative time crystal in an asymmetric non-linear photonic dimer, Phys. Rev. A 101, 1 (2019), arXiv:1910.03499 .
  * Carlon Zambon  _et al._ [2020] N. Carlon Zambon, S. R. Rodriguez, A. Lemaître, A. Harouri, L. Le Gratiet, I. Sagnes, P. St-Jean, S. Ravets, A. Amo, and J. Bloch, Parametric instability in coupled nonlinear microcavities, Phys. Rev. A 102, 23526 (2020), arXiv:1911.02816 .
  * Marconi  _et al._ [2020] M. Marconi, F. Raineri, A. Levenson, A. M. Yacomotti, J. Javaloyes, S. H. Pan, A. E. Amili, and Y. Fainman, Mesoscopic Limit Cycles in Coupled Nanolasers, Phys. Rev. Lett. 124, 213602 (2020), arXiv:1911.10830 .
  * Keeling  _et al._ [2010] J. Keeling, M. J. Bhaseen, and B. D. Simons, Collective Dynamics of Bose-Einstein Condensates in Optical Cavities, Phys. Rev. Lett. 105, 043001 (2010), arXiv:1002.3108 .
  * Bhaseen  _et al._ [2011] M. J. Bhaseen, J. Mayoh, B. D. Simons, and J. Keeling, Dynamics of Nonequilibrium Dicke Models, Phys. Rev. A - At. Mol. Opt. Phys. 85, 1 (2011), arXiv:1110.1348 .
  * Lee and Sadeghpour [2013a] T. E. Lee and H. R. Sadeghpour, Quantum Synchronization of Quantum van der Pol Oscillators with Trapped Ions, Phys. Rev. Lett. 111, 234101 (2013a), arXiv:1306.6359 .
  * Piazza and Ritsch [2015] F. Piazza and H. Ritsch, Self-organised Limit-Cycles, Chaos and Phase-Slippage with a Superfluid inside an Optical Resonator, Phys. Rev. Lett. 115, 1 (2015), arXiv:1507.08644 .
  * Chiacchio and Nunnenkamp [2019] E. I. R. Chiacchio and A. Nunnenkamp, Dissipation-Induced Instabilities of a Spinor Bose-Einstein Condensate Inside an Optical Cavity, Phys. Rev. Lett. 122, 193605 (2019), arXiv:1901.06996 .
  * Keßler  _et al._ [2019] H. Keßler, J. G. Cosme, M. Hemmerling, L. Mathey, and A. Hemmerich, Emergent limit cycles and time crystal dynamics in an atom-cavity system, Phys. Rev. A 99, 1 (2019), arXiv:1905.02751 .
  * Peters and Rodriguez [2022] K. J. Peters and S. R. Rodriguez, Limit cycles and chaos induced by a nonlinearity with memory, arXiv:2108.02680 (2022), arXiv:2108.02680 .
  * Kongkhambut  _et al._ [2022] P. Kongkhambut, J. Skulte, L. Mathey, J. G. Cosme, A. Hemmerich, and H. Keßler, Observation of a continuous time crystal, Science (80-. ). 377, 670 (2022), arXiv:2202.06980 .
  * Bello  _et al._ [2019] L. Bello, M. C. Strinati, E. G. D. Torre, and A. Pe’er, Persistent coherent beating in coupled parametric oscillators, Phys. Rev. Lett. 123, 83901 (2019), arXiv:1901.06202 .
  * Strinati  _et al._ [2020] M. C. Strinati, I. Aharonovich, S. Ben-Ami, E. G. D. Torre, L. Bello, and A. Pe’er, Coherent dynamics in frustrated coupled parametric oscillators, New J. Phys. 22, 85005 (2020), arXiv:2005.03562 .
  * Inagaki  _et al._ [2021] T. Inagaki, K. Inaba, T. Leleu, T. Honjo, T. Ikuta, K. Enbutsu, T. Umeki, R. Kasahara, K. Aihara, and H. Takesue, Collective and synchronous dynamics of photonic spiking neurons, Nat. Commun. 12, 2325 (2021).
  * Eichler and Zilberberg [2023] A. Eichler and O. Zilberberg,  _Classical and quantum parametric phenomena._ (Oxford University Press, 2023).
  * Chen  _et al._ [2016] C. Chen, D. H. Zanette, J. R. Guest, D. A. Czaplewski, and D. López, Self-Sustained Micromechanical Oscillator with Linear Feedback, Phys. Rev. Lett. 117, 1 (2016), arXiv:1612.03245 .
  * Houri  _et al._ [2019] S. Houri, D. Hatanaka, M. Asano, R. Ohta, and H. Yamaguchi, Limit cycles and bifurcations in a nonlinear MEMS resonator with a 1:3 internal resonance, Appl. Phys. Lett. 114, 103103 (2019).
  * Houri  _et al._ [2021] S. Houri, M. Asano, H. Okamoto, and H. Yamaguchi, Self-Sustained Libration Regime in Nonlinear Microelectromechanical Devices, Phys. Rev. Appl. 16, 064015 (2021).
  * Brenig [2022] W. Brenig, Spinless fermions in a Z2 gauge theory on a triangular ladder, Phys. Rev. B 105, 245105 (2022).
  * Khan and Türeci [2018] S. Khan and H. E. Türeci, Frequency combs in a lumped-element josephson-junction circuit, Physical review letters 120, 153601 (2018).
  * Metzger  _et al._ [2008] C. Metzger, M. Ludwig, C. Neuenhahn, A. Ortlieb, I. Favero, K. Karrai, and F. Marquardt, Self-Induced oscillations in an optomechanical system driven by bolometric backaction, Phys. Rev. Lett. 101, 26 (2008).
  * Lörch  _et al._ [2013] N. Lörch, J. Qian, A. Clerk, F. Marquardt, and K. Hammerer, Laser Theory for Optomechanics: Limit Cycles in the Quantum Regime, Phys. Rev. X 4, 1 (2013), arXiv:1310.1298 .
  * Bakemeier  _et al._ [2015] L. Bakemeier, A. Alvermann, and H. Fehske, Route to Chaos in Optomechanics, Phys. Rev. Lett. 114, 013601 (2015), arXiv:1407.5529 .
  * Navarro-Urrios  _et al._ [2017] D. Navarro-Urrios, N. E. Capuj, M. F. Colombano, P. D. Garc\́mathrm{i}a, M. Sledzinska, F. Alzina, A. Griol, A. Mart\́mathrm{i}nez, and C. M. Sotomayor-Torres, Nonlinear dynamics and chaos in an optomechanical beam, Nature Communications 8, 14965 (2017).
  * Roque  _et al._ [2019] T. F. Roque, F. Marquardt, and O. M. Yevtushenko, Nonlinear dynamics of weakly dissipative optomechanical systems, Nat. Rev. Phys. 2, 229 (2019), arXiv:1909.03126 .
  * Li  _et al._ [2022] J. Li, Z.-H. Zhou, S. Wan, Y.-L. Zhang, Z. Shen, M. Li, C.-L. Zou, G.-C. Guo, and C.-H. Dong, All-optical synchronization of remote optomechanical systems, Phys. Rev. Lett. 129, 063605 (2022).
  * Jenkins [2013] A. Jenkins, Self-oscillation, Physics Reports 525, 167 (2013).
  * Pikovsky  _et al._ [2001] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization - A Universal Concept in Nonlinear Sciences_ , Cambridge Nonlinear Science Series, Vol. 12 (Cambridge University Press, 2001) pp. 1–411.
  * Thomas  _et al._ [2004] J. P. Thomas, E. H. Dowell, and K. C. Hall, Modeling Viscous Transonic Limit Cycle Oscillation Behavior Using a Harmonic Balance Approach, J. Aircr. 41, 1266 (2004).
  * Petrov [2012] E. P. Petrov, Analysis of Flutter-Induced Limit Cycle Oscillations in Gas-Turbine Structures With Friction, Gap, and Other Nonlinear Contact Interfaces, J. Turbomach. 134, 10.1115/1.4006292 (2012).
  * Parker and Chua [2012] T. S. Parker and L. Chua,  _Practical numerical algorithms for chaotic systems_ (Springer Science & Business Media, 2012).
  * Han [2006] M. Han, Bifurcation theory of limit cycles of planar systems, in  _Handbook of Differential Equations: Ordinary Differential Equations_ , Vol. 3 (Elsevier, 2006) pp. 341–433.
  * Guckenheimer and Holmes [2013a] J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ , Vol. 42 (Springer Science & Business Media, 2013).
  * Hilbert [1900] D. Hilbert, Mathematische probleme, Nachrichten von der Koniglichen Gesellschaft der Wissenschaften zu Gottingen (1900).
  * Genesio and Tesi [1992] R. Genesio and A. Tesi, Harmonic balance methods for the analysis of chaotic dynamics in nonlinear systems, Automatica 28, 531 (1992).
  * Guskov  _et al._ [2008] M. Guskov, J.-J. Sinou, and F. Thouverez, Multi-dimensional harmonic balance applied to rotor dynamics, Mechanics Research Communications 35, 537 (2008).
  * Luo and Huang [2012] A. C. Luo and J. Huang, Approximate solutions of periodic motions in nonlinear systems via a generalized harmonic balance, J. Vib. Control 18, 1661 (2012).
  * Krack and Gross [2019] M. Krack and J. Gross,  _Harmonic Balance for Nonlinear Vibration Problems_ (Springer, 2019) p. 159.
  * 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), arXiv:2202.00571 .
  * Bogoliubov and Mitropolsky [1961] N. N. Bogoliubov and Y. A. Mitropolsky,  _Asymptotic methods in the theory of non-linear oscillations_ , Vol. 1 (Gordon and Breach New York, 1961).
  * Rand [2005] R. H. Rand, Lecture Notes on Nonlinear Vibrations, http://audiophile.tam.cornell.edu/randdocs/nlvibe52.pdf (2005).
  * 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, 22201 (2016), arXiv:1608.08896 .
  * Heugel  _et al._ [2022] T. L. Heugel, O. Zilberberg, C. Marty, R. Chitra, and A. Eichler, Ising machines with strong bilinear coupling, Phys. Rev. Res. 4, 013149 (2022).
  * Bachtold  _et al._ [2022] A. Bachtold, J. Moser, and M. I. Dykman, Mesoscopic physics of nanomechanical systems, Rev. Mod. Phys. 94, 045005 (2022).
  * Lugiato [1984] L. A. Lugiato, Ii theory of optical bistability (Elsevier, 1984) pp. 69–216.
  * Burgwal  _et al._ [2020] R. Burgwal, J. del Pino, and E. Verhagen, Comparing nonlinear optomechanical coupling in membrane-in-the-middle and single-cavity systems, New Journal of Physics 22, 113006 (2020).
  * Buishvili and Menabde [1979] L. Buishvili and M. Menabde, Application of the averaging method in problems involving high-resolution nuclear magnetic-resonance in solid bodies, ZHURNAL EKSPERIMENTALNOI I TEORETICHESKOI FIZIKI 77, 2435 (1979).
  * Sommese and Wampler [2005] A. J. Sommese and C. W. Wampler,  _Numer. Solut. Syst. Polynomials Aris. Eng. Sci._ (World Scientific, 2005) pp. 1–402.
  * Breiding and Timme [2018] P. Breiding and S. Timme, HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia, in  _Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics)_, Vol. 10931 LNCS (Springer, 2018) pp. 458–465, arXiv:1711.10911 .
  * Timme [2021] S. Timme,  _Numerical nonlinear algebra_ , Ph.D. thesis, Technische Universität Berlin (Germany) (2021).
  * Alaeian  _et al._ [2021] H. Alaeian, G. Giedke, I. Carusotto, R. Löw, and T. Pfau, Limit cycle phase and Goldstone mode in driven dissipative systems, Phys. Rev. A 103, 013712 (2021), arXiv:2007.10680 .
  * Chiacchio  _et al._ [2023] E. I. R. Chiacchio, A. Nunnenkamp, and M. Brunelli, The non-reciprocal dicke model (2023), arXiv:2302.06386 [quant-ph] .
  * Dreon  _et al._ [2022] D. Dreon, A. Baumgärtner, X. Li, S. Hertlein, T. Esslinger, and T. Donner, Self-oscillating pump in a topological dissipative atom–cavity system, Nature 608, 494 (2022).
  * Fernández [2018] F. M. Fernández, Exceptional point in a simple textbook example, European Journal of Physics 39, 045005 (2018).
  * Khedri  _et al._ [2022] A. Khedri, D. Horn, and O. Zilberberg, Fate of exceptional points in the presence of nonlinearities (2022), arXiv:2208.11205 [cond-mat.str-el] .
  * Guckenheimer and Holmes [2013b] J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ (Springer, 2013).
  * [66] See Supplemental Material at [URL will be inserted by publisher] for further details on the Extended Harmonic Balance approach and its convergence, in the case of a Van der Pol oscillator.
  * Borovik  _et al._ [2023] V. Borovik, P. Breiding, J. del Pino, M. Michałek, and O. Zilberberg, Khovanskii bases for semimixed systems of polynomial equations – a case of approximating stationary nonlinear newtonian dynamics (2023), arXiv:2306.07897 [math.AG] .
  * Liénard [1928] A. Liénard, Etude des oscillations entretenues, Revue Generale de l’Elactricite 23, 901 (1928).
  * Kuznetsov  _et al._ [2013] N. V. Kuznetsov, O. A. Kuznetsova, and G. A. Leonov, Visualization of four normal size limit cycles in two-dimensional polynomial quadratic system, Differential Equations and Dynamical Systems 21, 29 (2013).
  * Yan  _et al._ [2022] Y. Yan, X. Dong, L. Huang, K. Moskovtsev, and H. B. Chan, Energy transfer into period-tripled states in coupled electromechanical modes at internal resonance, Phys. Rev. X 12, 031003 (2022).
  * Aldana  _et al._ [2013] S. Aldana, C. Bruder, and A. Nunnenkamp, Equivalence between an optomechanical system and a Kerr medium, Phys. Rev. A - At. Mol. Opt. Phys. 88, 10.1103/PhysRevA.88.043826 (2013), arXiv:1306.0415 .
  * Chitra and Zilberberg [2015] R. Chitra and O. Zilberberg, Dynamical many-body phases of the parametrically driven, dissipative Dicke model, Phys. Rev. A - At. Mol. Opt. Phys. 92, 23815 (2015), arXiv:1501.07098 .
  * 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, 41046 (2021), arXiv:2104.12782 .
  * Wächtler and Platero [2023] C. W. Wächtler and G. Platero, Topological synchronization of quantum van der pol oscillators, Phys. Rev. Res. 5, 023021 (2023).
  * Moreno  _et al._ [2023] J. N. Moreno, C. W. Wächtler, and A. Eisfeld, Synchronized states in dissipatively coupled harmonic oscillator networks (2023), arXiv:2301.13614 [nlin.CD] .
  * Lee and Sadeghpour [2013b] T. E. Lee and H. R. Sadeghpour, Quantum synchronization of quantum van der pol oscillators with trapped ions, Phys. Rev. Lett. 111, 234101 (2013b).
  * Walter  _et al._ [2014] S. Walter, A. Nunnenkamp, and C. Bruder, Quantum synchronization of a driven self-sustained oscillator, Phys. Rev. Lett. 112, 094102 (2014).
  * Walter  _et al._ [2015] S. Walter, A. Nunnenkamp, and C. Bruder, Quantum synchronization of two van der pol oscillators, Annalen der Physik 527, 131 (2015).
  * Xu  _et al._ [2014] M. Xu, D. A. Tieri, E. C. Fine, J. K. Thompson, and M. J. Holland, Synchronization of two ensembles of atoms, Phys. Rev. Lett. 113, 154101 (2014).
  * Roulet and Bruder [2018a] A. Roulet and C. Bruder, Quantum synchronization and entanglement generation, Phys. Rev. Lett. 121, 063601 (2018a).
  * Sonar  _et al._ [2018] S. Sonar, M. Hajdušek, M. Mukherjee, R. Fazio, V. Vedral, S. Vinjanampathy, and L.-C. Kwek, Squeezing enhances quantum synchronization, Phys. Rev. Lett. 120, 163601 (2018).
  * Roulet and Bruder [2018b] A. Roulet and C. Bruder, Synchronizing the smallest possible system, Phys. Rev. Lett. 121, 053601 (2018b).
  * Arosh  _et al._ [2020] L. B. Arosh, M. C. Cross, and R. Lifshitz, Quantum limit-cycles and the Rayleigh and van der Pol oscillators, Phys. Rev. Res. 3, 013130 (2020), arXiv:2011.02706 .
  * Cox  _et al._ [2013] D. Cox, J. Little, and D. O’Shea,  _Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra_ (Springer, 2013).
  * Breiding  _et al._ [2022] P. Breiding, M. Michałek, L. Monin, and S. Telen, The algebraic degree of coupled oscillators (2022), arXiv:2208.08179 [math.AG] .
  * Landau [1944] L. D. Landau, On the problem of turbulence, in  _Dokl. Akad. Nauk USSR_ , Vol. 44 (1944) p. 311.
  * Stuart [1958] J. T. Stuart, On the non-linear mechanics of hydrodynamic stability, Journal of Fluid Mechanics 4, 1 (1958).

Supplemental Material: Limit cycles as stationary states of an extended
Harmonic Balance ansatz

##  I Extended Harmonic Balance method

The Harmonic Balance Method (HBM) is commonly applied to handle harmonically
driven nonlinear systems. Such systems are described by general nonlinear
system of $N$ second-order ordinary differential equations (ODEs) of the form

$\dot{\textbf{x}}(t)=G(\textbf{x}(t),t),$ (S1)

where $\textbf{x}(t)=(x_{1}(t),x_{2}(t),\cdots,x_{N}(t))^{T}$, and $x_{i}(t)$
are real variables with $i=1,2,...,N$ and time $t$ as the independent
variable. The function $G(\textbf{x}(t),t)$ is assumed to depend harmonically
on $\textbf{x}(t)$ and its time derivatives111Note that in a linear system,
$\mathbf{G}(\mathbf{x},t)=\ddot{\mathbf{x}}+\mathbf{M}(t)\mathbf{x}+\mathbf{b}(t)$,
where $\mathbf{M}(t)$ contains spring constants and linear couplings, and
$\mathbf{b}(t)$ represents external forces. Diagonalizing $\mathbf{M}(t)$
yields the normal modes of the system.. The system (S1) eventually reaches a
stationary state, typically exhibiting an oscillatory response that remains
constant over time. The HBM starts with the premise of representing such a
solution with an expansion of the form [see main text Eq. (2)]

$\displaystyle
x_{i}(t)=\sum_{l=1}^{M_{i}}u_{l,i}(t)\cos(\omega_{l,i}t)+v_{l,i}(t)\sin(\omega_{l,i}t)\,,$
(S2)

where we introduced $M_{i}$ harmonics oscillating at frequencies
$\omega_{l,i}$ with slowly-varying envelopes $u_{l,i}$ and $v_{l,i}$, in order
to describe the evolution of each variable $x_{i}(t)$. An approach to finding
such amplitudes is to perform numerical experiments, where ODE solvers
propagate the evolution for a long time and Fourier analysis of the solutions
is performed. Once a given set of solutions has been found for a parameter
set, additional solutions can be obtained by continuation techniques [47].
This methodology, however, can be challenging to apply, since nonlinear
systems can have multiple stationary states, resulting in different responses
depending solely on the initial conditions. Obtaining a complete stationary-
state diagram requires time-consuming sampling of initial conditions, which
may not guarantee the discovery of all branches of stationary solutions.

To simplify the process, an alternative approach involves using Eq. (S2) as an
ansatz of the ODE system. Replacement of the expansion (S2) into (S1) leads to
harmonic oscillating terms at both sides of the equation. The stationary state
problem can be transformed into an algebraic one, within the truncated ansatz
(S2), by ‘balancing’ harmonics at both sides of Eq. (S1). This procedure
results in equations equivalent to those found by Krylov-Bogoliubov averaging
[50] (see main text and Ref. [48] for further details):

$\mathcal{F}(\textbf{u},\textbf{v})=0,$ (S3)

where $\textbf{u}=(u_{1,1},\cdots,u_{M_{N},N})$,
$\textbf{v}=(v_{1,1}\cdots,v_{M_{N},N})$.

In typical scenarios, an externally driven system settles into stationary
motion states that oscillate at frequencies that are equal or commensurate
with the driving frequencies. Therefore, the frequencies $\omega_{l,i}$ can be
assumed a priori (e.g., they are prescribed by the driving frequencies), while
the variables $\textbf{u},\textbf{v}$ remain unknown. In limit cycle (LC)
problems, however, this is no longer the case, since the complex interplay of
driving, dissipation, and nonlinearity determines the oscillation frequencies.
In order to solve such frequencies, we resort to the concept of gauge freedom
in LCs. Gauge freedom fundamentally stems from the spontaneous emergence of
LCs, which causes a striking change in the time-translation symmetry of the
problem: a time-dependent solution can arise even in a time-independent
system. We illustrate how this concept allows to extract limit cycles as the
non-redundant stationary solutions from an HBM ansatz, upgrading the HBM into
the extended HBM (e-HBM).

In the main text, the concept of gauge freedom enables the determination of
the system of polynomial equations for the stationary state by identifying a
specific, non-redundant solution. Here, we delve into the emergence of gauge
freedom and discuss in the case of $N=1$ to keep the discussion simple. We
deduce the notion of gauge freedom from the LC symmetry under time
translations with the ODE’s period $T_{s}$, namely the operation
$\mathcal{T}(T_{s})x(t)=x(t+T_{s})$, or equivalently,
$\mathcal{T}(T_{s})a_{m}=a_{m}e^{im\omega_{m}T_{s}}$. A periodic nonlinear
system is invariant under the discrete operation $\mathcal{T}(T_{s})$, so it
is a time-independent (autonomous) system with a continuous, arbitrary
$T_{s}$. However, their long-term solutions can spontaneously break these
symmetries. For instance sub-harmonic oscillations with commensurate periods
$T_{s}q/p$ ($p,q\in\mathbb{N}$) break $\mathcal{T}(T_{s})$ symmetry. This
allows all solutions returning to themselves after $q$ iterations of the
original translation $\mathcal{T}(T_{s})$, which are related via discrete
rotations $a_{m}\rightarrow e^{i2\pi j/q}a_{m}$ for
$j\in\\{0,1,\cdots,q-1\\}$. However, a continuous time-translation-symmetry
broken LC has an incommensurate period $T_{s}r$ with $r$ irrational;
mathematically, LC orbits only return to themselves after infinite actions of
$\mathcal{T}(T_{s})$ with fixed $T_{s}$. This property frees oscillation
phases in LCs from constraints and yields a continuous set of solutions
related by full $U(1)$ symmetry $a_{i}\mapsto e^{i\varphi_{i}}a_{i}$, where
$\varphi_{i}\in(0,2\pi)$.

##  II Convergence of the e-HBM for the Van der Pol oscillator

Here, we provide details on the application of the e-HBM to the Van der Pol
(VdP) oscillator in the main text. For the sake of simpler algebraic
manipulation, we adopt a complex notation where
$\textbf{a}=\textbf{u}+i\textbf{v}$, thus writing the HBM ansatz as

$x(t)=\sum_{m=1}^{M}a_{m}(t)e^{-im\omega
t}+\mathrm{H.c.=}\sum_{k=0}^{\frac{M-1}{2}}a_{k}(t)e^{-i(2k+1)\omega
t}+\mathrm{H.c.}\,.$ (S4)

The first step of the HBM involves inserting Eq. (S4) into the VdP equation,
namely

$\ddot{x}-\mu(1-x^{2})\dot{x}=0\,,$ (S5)

as we illustrated in Eq. (3) of the main text. Next, we match oscillatory
terms with a specific frequency $\omega_{k}$ and isolate the corresponding
frequency terms. Such frequency matching procedure filters each oscillating
contribution in an equivalent way as averaging the equation over the period
$2\pi/\omega_{k}$, such that $\omega_{k}=k\omega$.

The nonlinearity in the VdP equation leads to couplings between different
harmonic amplitudes. We focus on the impact of this nonlinear damping in the
VdP oscillator. The corresponding term in the Harmonic Balance equation at a
frequency $\omega_{k}$ follows from the average integral:

$\langle
x^{2}\dot{x}\rangle^{c}=-i\frac{\omega_{k}}{2\pi}\int_{-2\pi/\omega_{k}}^{2\pi/\omega_{k}}\mathrm{d}t\sum_{m,n,p}\omega_{m}\left[a_{m}e^{-i\omega_{m}t}-a_{m}^{*}e^{i\omega_{m}t}\right]\left[a_{n}e^{-i\omega_{n}t}+a_{n}^{*}e^{i\omega_{n}t}\right]\left[a_{p}e^{-i\omega_{p}t}+a_{p}e^{i\omega_{p}t}\right]e^{-i\omega_{k}t}\,.$
(S6)

with $a_{m}(t)$ assumed constant for the purpose of the time integration [see
[48]].

The equations for the harmonic amplitudes $u_{k}$ and $v_{k}$ can be derived
by separating the real and imaginary terms from the aforementioned expression.
The key point is that, for each $k$, all terms with $m,n\leq M$ are present,
and explicitly, the surviving terms are

$\displaystyle\langle
x^{2}\dot{x}\rangle_{k}^{c}=i\sum_{m,n,p}\omega_{m}\left[a_{m}a_{n}\left[a_{p}\delta_{m+n+p+k}+a_{p}^{*}\delta_{m+n-p+k}\right]-a_{m}^{*}a_{n}\left[a_{p}\delta_{-m+n+p+k}+a_{p}^{*}\delta_{-m+n-p+k}\right]\right]+$
$\displaystyle
i\sum_{m,n,p}\omega_{m}\left[a_{m}a_{n}^{*}\left[a_{p}\delta_{m-n+p+k}+a_{p}^{*}\delta_{m-n-p+k}\right]-a_{m}^{*}a_{n}^{*}\left[a_{p}\delta_{-m-n+p+k}+a_{p}^{*}\delta_{-m-n-p+k}\right]\right]\,,$
(S7)

where $\delta_{m-n}=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{i(m-n)\theta}\,d\theta$
stands for the Kronecker symbol. The preceding step assumes the
commensurability of frequencies, which has been previously assumed in the Van
der Pol ansatz Eq. (S4). The ansatz insertion, expansion, and frequency
balance are facilitated in the main text by the symbolic layer of the
HarmonicBalance.jl package [48].

Equation (II) reveals explicit nonlinear coupling between different harmonics,
weighted by ‘selection rules’ arising from the Kronecker delta factors. Under
the assumption of convergence of the expansion Eq. (S4), we expect some of the
harmonics to play a more predominant role than others: in practice, we observe
the harmonic amplitudes to monotonically decay, cf. Fig. 2(b) in main text.
Therefore, to assess convergence, it is instrumental to analyze how higher-
frequency harmonics impact low-frequency ones. For this analysis, it is
important to also notice that the LC in the VdP oscillator exhibits a series
of step functions in the time domain that oscillate with half-wave symmetry
with period $T$, ($x(t)=-x(t\pm T/2)$), i.e., only odd harmonics contribute to
the ansatz. The sums $\sum_{m,n,p}$ must be thus restricted to odd harmonics
only, so we replace $m\rightarrow 2m+1,n\rightarrow 2n+1,p\rightarrow
2p+1,k\rightarrow 2k+1$ in the harmonic frequencies and set the origin in the
indices to zero: $m,n,p,k\geq 0$.

Let as analyze the terms that arise from a given truncation $M$. Keeping only
the first harmonic in the ansatz, i.e., setting $M=1$ and collecting terms
with $m=n=p=k=0$, the only non-vanishing contributions in section II are

$\displaystyle-
a_{m}^{*}a_{n}^{*}a_{p}\delta_{2(k-m-n+p)}=-a_{0}^{*}|a_{0}|^{2},$
$\displaystyle-
a_{m}^{*}a_{n}a_{p}^{*}\delta_{2(k-m+n-p)}=-a_{0}^{*}|a_{0}|^{2},$
$\displaystyle
a_{m}a_{n}^{*}a_{p}^{*}\delta_{2(k+m-n-p)}=a_{0}^{*}|a_{0}|^{2}\,.$ (S8)

At this approximation level, equivalent to the lowest-order averaging method
[see [50] for further details], the nonlinear dissipation in the VdP
oscillator appears as an amplitude-dependent dispersive shift. It simplifies
to the Stuart-Landau equation [86, 87], commonly used to investigate nonlinear
oscillator systems near a Hopf bifurcation.

We can further analyze the corrections to this lower harmonic $a_{0}$ from the
interaction between $a_{0}$ with $a_{1}$. Keeping $m,n,p\leq M=1$ leads to

$\displaystyle\langle x^{2}\dot{x}\rangle_{k=0}^{c}=$ $\displaystyle-
a_{0}^{*}|a_{0}|^{2}+a_{0}^{2}a_{1}^{*}-2a_{1}a_{0}^{*}a_{1}\,.$ (S9)
Similarly, keeping $m,n,p\leq M=2$, introduces additional corrections,
$\displaystyle\langle x^{2}\dot{x}\rangle_{k=0}^{c}=$ $\displaystyle
a_{0}^{2}a_{1}^{*}-a_{0}\left(a_{0}^{*}{}^{2}-2a_{1}a_{2}^{*}\right)-2a_{1}a_{0}^{*}a_{1}^{*}-a_{2}\left(a_{1}^{*}{}^{2}+2a_{0}^{*}a_{2}^{*}\right)\,,$
(S10) where $a_{0}$ and $a_{2}$ are directly coupled, but also $a_{0}$, where
$a_{0}$, $a_{1}$ and $a_{2}$ have a threefold interaction. At increasing
orders, nonlinear contributions between all harmonics proliferate, see
$\displaystyle\langle x^{2}\dot{x}\rangle_{k=0}^{c}=$ $\displaystyle
a_{0}^{2}a^{*}_{1}+a_{0}\left(2(a_{1}a^{*}_{2}+a_{2}a^{*}_{3}+a_{3}a^{*}_{4})-a^{*}_{0}{}^{2}\right)-a_{2}a^{*}_{1}{}^{2}-a_{4}a^{*}_{2}{}^{2}-2a_{2}a^{*}_{0}a^{*}_{2}-2a_{3}a^{*}_{1}a^{*}_{2}+a_{1}^{2}a^{*}_{3}-$
$\displaystyle
2a^{*}_{3}(a_{3}a^{*}_{0}+a_{4}a^{*}_{1})+a_{1}(2a_{2}a^{*}_{4}-2a^{*}_{0}a^{*}_{1})\,.$
(S11)

The Fourier analysis [cf. Fig.2(b) in main text] shows $a_{k}\mapsto 0$ as $k$
increases. Given the prefactors in that the terms above are comparable, the
convergence of the ansatz entails that the terms $a_{m}a_{n}a_{p}$ are
increasingly weaker as they contain harmonics for larger $m,n,p$, and the
result is expected to converge. However, further analysis is needed to
understand if such convergence is monotonic. In Fig. S1 we show how, in fact,
this is not the case: Increasing the number of harmonics in the ansatz $M$
does not necessarily increase the fidelity of the solution.

Figure S1: Fidelity of e-HBM ansatz, further details: (a)-(c) Overlap between
e-HBM and exact time evolution result as a function of the number of harmonics
$M$ for varying $\mu$, shown at each panel [cf. Fig. 2 in the main text]. The
overlap is computed by matching frequency peaks in $|\tilde{x}(\omega)|^{2}$,
from the quasi-exact evolution $\textbf{a}_{\text{FT}}$ (cf. main text) to the
nearest frequency estimates obtained by the HBM $\textbf{a}_{\text{HBM}}$,
with percentile deviations below $2\%$, shown in (d).

## References