11institutetext: Departamento de Ciencias Integradas y Centro de Estudios
Avanzados en Física, Matemática y Computación, Universidad de Huelva, 21071
Huelva, Spain 22institutetext: Departamento de Física Atómica, Molecular y
Nuclear, Facultad de Física, Universidad de Sevilla, Apartado 1065, E-41080
Sevilla, Spain 33institutetext: Instituto Carlos I de Física Teórica y
Computacional, Universidad de Granada, Fuentenueva s/n, 18071 Granada, Spain

# At the edge of shape coexistence in the Z=40 region

E. Maya-Barbecho 11<EMAIL_ADDRESS>S. Baid 22<EMAIL_ADDRESS>J.M.
Arias 2233<EMAIL_ADDRESS>J.E. García-Ramos 1133<EMAIL_ADDRESS>
###### Abstract

In this contribution, the shape coexistence phenomenon near the proton sub-
shell closure at Z=40 is analyzed. Particular emphasis is placed on extracting
the nuclear deformation values by examining experimental B(E2) transition
rates and observing how the kinematic moment of inertia evolves in these
nuclei. Based on the analysis of these two observables, we arrive at the
conclusion that the notable effects of shape coexistence at approximately Z=40
are largely diminished in Mo isotopes and are scarcely detectable in Ru
isotopes.

##  1 Introduction

Shape coexistence is a phenomenon observed not only in atomic nuclei but also
in other systems, such as certain molecules. It is known to occur in various
nuclear mass regions, with a particular prevalence around neutron and proton
shell closures heyde11 ; Garr22 .

In nuclei near a shell closure, the tendency is for them to be spherical, and
particle-hole excitations in such regions are energetically costly. However,
the energy of these excitations can be significantly reduced due to the
enhanced residual interaction brought about by the presence of more effective
valence nucleons heyde11 . Consequently, these excitations correspond to more
deformed states. This effect becomes particularly significant when protons are
near a shell closure, while neutrons are situated in the middle of a shell, or
vice versa. In this scenario, shape coexistence arises in nuclei, where
regular states (0-particle-0-hole) coexist with intruder states
(2-particle-2-hole) within the same energy range despite having distinct
deformations. In extreme cases, intruder states can even have lower energies
than regular states Garc19 .

It is essential to note that deformation in nuclear structure is not a
directly measurable quantity but must be inferred from certain observable
data. In this contribution, we will extract information regarding the
deformation of isotopes of Sr, Zr, Mo, and Ru. We will then employ this
information to explore how the phenomenon of shape coexistence evolves as we
move further away from the Z=40 subshell closure. Of particular interest is
the examination of the onset of deformation in different isotopic chains and
its potential connection with Quantum Phase Transitions (QPTs) Cejn09 .

##  2 The “crossing” of configurations

Figure 1: Unperturbed regular (black dots and lines) and intruder (red squares
and lines) ground state energies for Sr, Zr, Mo, and Ru isotopes as a function
of the neutron number.

Identifying regular and intruder states in a given nucleus is a non-trivial
task, and except for the most evident cases, it often requires the use of a
theoretical model. One particularly suitable choice is the Interacting Boson
Model (IBM) iach87  with Configuration Mixing (IBM-CM) duval82 . It
incorporates 2p-2h excitations and is valuable for describing the key
spectroscopic properties of nuclei providing a microscopic representation of
the wave function. With this model, it becomes feasible to determine whether a
given state belongs to the regular or the intruder sector.

The IBM-CM has been employed in the study of isotopes of Sr Maya2022 , Zr
Garc19 , Mo and Ru isotopes Maya2023 . In Figure 1, we present a visualization
of the lowest energy states from both the regular and intruder sectors as
extracted from these calculations. Importantly, these energies are computed
with the coupling between the two sectors removed, providing a clear
separation of pure regular and intruder states.

One striking observation from Figure 1 is that regular and intruder states are
often quite close in energy, and in certain cases, they even cross paths.
These crossings are evident in Sr, Zr, and Mo isotopes, particularly around
the neutron number N=60. Notably, even in Ru isotopes, where crossings are
absent, the closest approach between regular and intruder states occurs around
the same neutron number, N=60. This region is well-known for its rapid
evolution of ground state structure and the presence of a QPT at N=60 Cejn09 .
Consequently, it is reasonable to connect the presence of the QPT with the
occurrence of crossings between regular and intruder configurations as
elucidated in Gavr19 ; Gavr22 .

##  3 The deformation

Figure 2: Value of the deformation parameter, $\beta$, for the ground state
band in Sr, Zr, Mo, and Ru isotopes extracted from $B(E2:2_{1}^{+}\rightarrow
0_{1}^{+})$ as a function of the neutron number.

To extract information about the quadrupole deformation, $\beta$, for a given
band we will rely on the geometrical view of the nucleus extracting the
$\beta$ parameter from a known $B(E2)$ value as,

$\beta=\frac{4\pi\sqrt{B(E2:I_{i}\rightarrow I_{f})}}{3\langle
I_{i}020|I_{f}0\rangle ZeR_{0}^{2}},$ (1)

where $R_{0}=1.2A^{1/3}$ fm and $<....|..>$ is the Clebsch-Gordan coefficient
and for the case of the yrast band we use $I_{i}=2$ and $I_{f}=0$.

In Figure 2, we present the extracted values of the deformation parameter,
$\beta$, for the ground state bands in the four isotopic chains. It is evident
that these chains exhibit distinct and contrasting systematics. In the cases
of Sr and Zr, there is a clear and rapid onset of deformation at neutron
number N=60. For Mo, the increase in deformation is more gradual. In the case
of Ru, the value of $\beta$ remains nearly constant. In the instances of Sr
and Zr, the sudden increase in deformation arises from the crossing of regular
and intruder configurations, as previously shown in Figure 1. Notably, the
intruder configuration has a significantly larger deformation, contributing to
the rapid change in $\beta$. However, in the case of Mo, although the crossing
of regular and intruder configurations exists, the strong interaction between
these configurations effectively smooths out the variation in $\beta$ Maya2023
. This smoothing effect has also been observed in Pt Garc09 . These
observations are further corroborated by the systematics of the two-neutron
separation energy Garc19 ; Maya2022 ; Maya2023 . In the language of QPTs, Sr
and Zr exhibit a first order Type II QPT, as seen in the behavior of $\beta$
Gavr19 ; Gavr22 . Similarly, Mo showcases a second order Type II QPT. In
contrast, Ru provides ample evidence of a second-order QPT of Type I Maya2023
.

##  4 The moment of inertia

Figure 3: Value of the kinematic moment of inertia for Sr, Zr, Mo, and Ru
isotopes as a function of the $\gamma$ energy, for different neutron numbers.

An alternative and valuable approach for gaining insights into the evolving
collective structure along the yrast band, with a particular focus on
deformation, is to examine the variations in the moment of inertia as one
progresses up the band structure. In this context, we will specifically
investigate the kinematic moment of inertia, denoted as J(1), which is defined
as follows Bohr75 ,

$J^{(1)}=\frac{\hbar^{2}}{2}\left(\frac{dE}{d\left(J(J+1)\right)}\right)^{-1}\approx\frac{\hbar^{2}(2J-1)}{E_{\gamma}(J\rightarrow(J-2))},$
(2)

where $E_{\gamma}(J\rightarrow(J-2))$ is the energy difference $E(J)-E(J-2)$.

In constructing the experimental values of J(1), we rely on the yrast band
energies. However, complexities are encountered at high-spin levels (8+, 10+,
and 12+). These complexities arise from the crossing of non-collective states,
such as those associated with broken-pair configurations, and the crossing of
regular and intruder configurations at specific energy levels. Such crossings
can lead to significant variations in the smooth progression of J(1) sometimes
resulting in a backbending pattern.

In Figure 3, we present the kinematic moment of inertia for the different
isotope chains. To aid in visualization, it is helpful to note that a purely
rotational system exhibits a constant moment of inertia, while a pure
vibrational structure appears as a vertical trend. For Sr isotopes a
significant change in behavior is observed from N=60 onward, where the moment
of inertia demonstrates a relatively smooth trend, indicative of a rotational
structure. However, for N=58, it appears that the lowest angular momenta
correspond to a more vibrational-like structure, while the highest angular
momenta exhibit rotational characteristics. For N=54-56, it is challenging to
draw a definitive conclusion, though the substantial values of Eγ (gamma-ray
energy) suggest a vibrational system. The results points towards the
coexistence of two configurations, with the rotational structure corresponding
to the intruder configuration and the vibrational structure to the regular
one. For Zr and Mo isotopes, similarly to Sr, the information presented in
Figure 3 supports the existence of two different configurations. From N=60 and
onwards, there is a clear presence of a rotational structure. However, for
N=54-60, the observed backbending suggests the coexistence of two structures
within the yrast band. In the Ru isotopes, a relatively consistent behavior is
observed across the entire chain. Nevertheless, the highest angular momenta
exhibit an upward trend, suggesting the presence of a vibrational structure,
while the lowest angular momenta maintain a relatively horizontal trend,
typical of a more deformed structure. In summary, the kinematic moment of
inertia serves as a suitable observable for studying the coexistence of
different configurations. For Sr, Zr, and Mo, the data suggests the presence
of two crossing configurations, while in Ru, a single configuration exists,
showing more of a rotational nature for lower angular momenta and a more
vibrational character for higher angular momenta.

##  5 Conclusions

This contribution presents a comprehensive study investigating the onset of
deformation and its relationship with shape coexistence around the Z=40
region, focusing on the isotope chains 38Sr, 40Zr, 42Mo, and 44Ru. The study
employs experimental data, including the deformation parameter, $\beta$, and
kinematic moment of inertia, derived from B(E2) values and excitation
energies. The key findings reveal the significance of shape coexistence in
explaining energy systematics for Sr and Zr, while Mo exhibits a notable
impact from intruder and regular configurations. However, this impact
diminishes for Ru isotopes, with N=60 marking a boundary between spherical and
deformed nuclei. The study also relates these results to Quantum Phase
Transitions (QPTs), with Sr and Zr experiencing first-order Type II QPTs at
N=60, Mo showing a second-order Type II QPT, and Ru undergoing a second-order
Type I QPT at N=60. This study, complementing Gavr19 ; Garc20 ; Gavr22 ;
Maya2022 ; Maya2023 , sheds light on the complex interplay between nuclear
structure, coexisting configurations, and quantum phase transitions.

##  6 Acknowledgements

This work was partially supported by grants PID2019-104002GB-C21,
PID2019-104002GB-C22, PID2022-136228NB-C21, and PID2022-136228NB-C22 funded by
MCIN/AEI/10.13039/50110001103 and “ERDF A way of making Europe”. It was also
partially supported by the ConsejerÃa de EconomÃa, Conocimiento, Empresas y
Universidad de la Junta de AndalucÃa (Spain) under Group FQM-160, and under
Projects No. P20-01247, and No. US-1380840. Resources supporting this work
were provided by the CEAFMC and Universidad de Huelva High Performance
Computer (HPC@UHU) funded by ERDF/MINECO project UNHU-15CE-2848.

## References

  * (1) K. Heyde, J.L. Wood, Rev. Mod. Phys. 83, 1467 (2011)
  * (2) P.E. Garrett, M. ZieliÅska, E. ClÃment, Progress in Particle and Nuclear Physics 124, 103931 (2022)
  * (3) J.E. García-Ramos, K. Heyde, Phys. Rev. C 100, 044315 (2019)
  * (4) P. Cejnar, J. Jolie, Progress in Particle and Nuclear Physics 62, 210 (2009)
  * (5) F. Iachello, A. Arima, _The interacting boson model_ (Cambridge University Press, Cambridge, 1987)
  * (6) P.D. Duval, B.R. Barrett, Nucl. Phys. A 376, 213 (1982)
  * (7) E. Maya-Barbecho, J.E. García-Ramos, Phys. Rev. C 105, 034341 (2022)
  * (8) E. Maya-Barbecho, S. Baid, J.M. Arias, J.E. García-Ramos, Phys. Rev. C 108, 034316 (2023)
  * (9) N. Gavrielov, A. Leviatan, F. Iachello, Phys. Rev. C 99, 064324 (2019)
  * (10) N. Gavrielov, A. Leviatan, F. Iachello, Phys. Rev. C 105, 014305 (2022)
  * (11) J. García-Ramos, K. Heyde, Nucl. Phys. A 825, 39 (2009)
  * (12) A. Bohr, B.R. Mottelson, _Nuclear Structure, Vol. II_ (W. A. Benjamin, Reading MA, 1975)
  * (13) J.E. García-Ramos, K. Heyde, Phys. Rev. C 102, 054333 (2020)


# A biased Ising model using two coupled Kerr parametric oscillators with
external force

Pablo Álvarez  Laboratory for Solid State Physics, ETH Zürich, CH-8093 Zürich,
Switzerland.  Davide Pittilini  Laboratory for Solid State Physics, ETH
Zürich, CH-8093 Zürich, Switzerland.  Filippo Miserocchi  Laboratory for Solid
State Physics, ETH Zürich, CH-8093 Zürich, Switzerland.  Sathyanarayanan
Raamamurthy  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.  Orjan Ameye  Department of Physics,
University of Konstanz, D-78457 Konstanz, Germany.  Javier del Pino  Institute
for Theoretical 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.  Quantum Center, ETH Zurich, CH-8093 Zurich,
Switzerland

###### Abstract

Networks of coupled Kerr parametric oscillators (KPOs) are a leading physical
platform for analog solving of complex optimization problems. These systems
are colloquially known as “Ising machines”. We experimentally and
theoretically study such a network under the influence of an external force.
The force breaks the collective phase-parity symmetry of the system and
competes with the intrinsic coupling in ordering the network configuration,
similar to how a magnetic field biases an interacting spin ensemble.
Specifically, we demonstrate how the force can be used to control the system,
and highlight the crucial role of the phase and symmetry of the force. Our
work thereby provides a method to create Ising machines with arbitrary bias,
extending even to exotic cases that are impossible to engineer in real spin
systems.

††preprint: APS/123-QED

The Kerr parametric oscillator (KPO) is a nonlinear resonator with a time-
dependent harmonic potential term Ryvkine and Dykman (2006); Mahboob and
Yamaguchi (2008); Wilson  _et al._ (2010); Eichler  _et al._ (2011, 2018a);
Gieseler  _et al._ (2012); Lin  _et al._ (2014); Puri  _et al._ (2017a);
Eichler  _et al._ (2018b); Nosan  _et al._ (2019); Frimmer  _et al._ (2019);
Grimm  _et al._ (2019); Wang  _et al._ (2019); Puri  _et al._ (2019); Miller
_et al._ (2019); Yamaji  _et al._ (2022). In a certain range of parameters,
this potential modulation renders the zero-amplitude solution unstable. There,
the system undergoes a spontaneous period-doubling $Z_{2}$ symmetry-breaking
phase transition and assumes a large-amplitude oscillation at a so-called
phase state. Importantly, a KPO features two such phase states with equal
amplitude but phases separated by $\pi$, shown as black dots in Fig. 1(a).

The KPO is at the focus of much research work because its two phase states are
analogous to the two polarization states of an Ising spin, “up” and “down”.
Consequently, it was proposed that networks of KPOs can be used to find the
ground state of Ising Hamiltonians, that is, the energetically preferred
configuration of a spin network Ising (1925). Such resonator-based Ising
solvers Mahboob  _et al._ (2016); Goto (2016); Puri  _et al._ (2019); Bello
_et al._ (2019); Okawachi  _et al._ (2020) are of high interest because the
corresponding calculations are NP-hard to tackle with conventional computers
Mohseni  _et al._ (2022), and yet they map to many other key optimization
problems, such as the travelling salesman problem Lucas (2014), the MAX-CUT
problem Inagaki  _et al._ (2016a); Goto  _et al._ (2019), and the number
partitioning problem Nigg  _et al._ (2017). In Fig. 1(b), we sketch a network
of two spins, where each spin takes the form of a double-well potential whose
wells corresponds the two levels (“spin up” or “spin down”).

Figure 1: KPOs as Ising spin analogues. (a) KPO phase states in a rotating
phase space spanned by the quadratures $u$ and $v$, where $x=u\cos(2\pi
f_{d}t)-v\sin(2\pi f_{d}t)$, with $f_{d}$ the driving frequency. Stationary
solutions in the absence of an external force are shown as black dots and are
labelled “up” and “down” to indicate a formal similarity to the two states of
an Ising spin, marked by arrow symbols. Solution in the presence of an
external force $F$ are shown as white dots as a function of the force phase
$\theta$ (dotted circle). (b) Representation of a 2-spin network as coupled
double wells in the absence of a bias field. The levels of each double well
($1,2$) are energetically degenerate. The spin-spin coupling is indicated by
$J$ and leads to an antisymmetric state in this example. (c) In the presence
of a bias field $B$, the degeneracy is broken and each spin has a preferred
polarization. This preference can overcome the solution favored by the
coupling $J$. (d) Coupled resonators are often described in a normal-mode
basis as illustrated by two coupled pendula. For two KPOs, the relevant
nonlinear stationary states we will refer to in this paper include states with
zero amplitude (Z), symmetric (S) and antisymmetric (A) oscillation, as well
as mixed-symmetry (M) solutions.

In the low-amplitude limit, the oscillations of a resonator network are
approximately harmonic and can be understood in terms of their normal modes. A
two-KPO system starting from the zero-amplitude (Z) solution can therefore
ring up to symmetric (S) or antisymmetric (A) phase states that map to
ferromagnetic or antiferromagnetic Ising configurations, respectively, see
Fig. 1(d) Heugel  _et al._ (2019). At large amplitudes, the nonlinearities
become significant and can lead to deviations from a normal-mode basis. The
deviations manifest for example as oscillations of mixed symmetry (M) Heugel
_et al._ (2022), which have no clear Ising counterpart. Nevertheless, the
normal-mode basis remains useful as a frame of reference to study both
strongly Heugel  _et al._ (2022) and weakly coupled KPOs Margiani  _et al._
(2023).

Figure 2: (a) Measured amplitudes $A_{i}$ and phase responses $\phi_{i}$ of
the two resonators when driven by a parametric pump tone with
$U_{p}=$3.5\text{\,}\mathrm{V}$$. The arrow symbols represent the Ising spin
analogy, and the colors indicate the relative phase configuration. (b)
Measured phase $\phi_{2}$ as a function of $f_{d}$ and $U_{p}$ for $F=0$. (c)
Measured amplitude and phase configuration of the two resonators for $F=0$.
Dashed lines are theory predictions for the outline of the characteristic
double-tongue. (d) Same as (b) with a force strength
$F=$10\text{\,}\mathrm{mV}$$ at a phase
$\theta=$55\text{\,}\mathrm{\SIUnitSymbolDegree}$$. (e) Measured amplitude and
phase configuration of the two resonators for $F=$10\text{\,}\mathrm{mV}$$ and
$\theta=$55\text{\,}\mathrm{\SIUnitSymbolDegree}$$. Dashed lines mark the
outlines calculated for (c) to highlight the shifted boundaries.

Previous experimental and theoretical studies of KPO-based Ising simulators
with bilinear coupling considered the case of unbiased Ising Hamiltonians,
where the solutions are only defined up to a global sign Goto (2016); Puri
_et al._ (2017b); Nigg  _et al._ (2017); Goto  _et al._ (2018); Dykman  _et
al._ (2018); Rota  _et al._ (2019); Calvanese Strinati  _et al._ (2019);
Heugel  _et al._ (2022); Margiani  _et al._ (2023). For example, the state
“down-up” shown in Fig. 1(b) can be identically replaced by “up-down”, as the
individual spin levels are degenerate. Many Ising problems, however, require
breaking of this degeneracy; the archetypal case being a magnetic field
applied to a spin ensemble. In Fig. 1(c), a magnetic field $B$ biases the
potential to compete with the spin-spin coupling $J$. A strong field can
overcome the coupling-induced ordering and dramatically change the
configuration that emerges as the optimal solution of the spin system. This
functionality was recently included in optical parametric systems with
dissipative coupling Takesue  _et al._ (2020), which constitute an alternative
route to KPO networks Inagaki  _et al._ (2016b); Gershenzon  _et al._ (2020).

In this paper, we demonstrate experimentally how the functionality of an
external $B$ field can be introduced in a classical KPO network. By applying
an external force $F$ to each resonator, we can displace the network’s
solutions in their phase space, as indicated in Fig. 1(a). Our experiment
features two coupled KPOs, but the concepts we present are easily extended to
larger networks. We provide a general framework to understand the role of the
external force term for coupled KPOs and resonator-based Ising solvers.
Importantly, while this forcing can emulate the effect of a simple bias field,
we show that its consequences are much richer due to the freedom of selecting
the driving amplitude and phase for each resonator individually. Our study
will therefore not only provide a practical guide for applications, but also
motivate further fundamental research on Ising networks with inhomogeneous
magnetic fields, including quantum implementations Puri  _et al._ (2017b);
Grimm  _et al._ (2019); Dykman  _et al._ (2018); Goto (2016).

Our experiment consists of two electrical resonators that feature a nonlinear
capacitance Nosan  _et al._ (2019). The system can be described by the coupled
differential equations

$\displaystyle\ddot{x}_{i}+\omega_{0}^{2}\left[1-\lambda\cos\left(4\pi
f_{d}t\right)\right]x_{i}+\beta
x_{i}^{3}+\Gamma\dot{x}_{i}-Jx_{i}=F_{i}(t)\,,$ (1)

where the displacement $x_{i}$ is the measured voltage signal of resonator $i$
Nosan  _et al._ (2019); Heugel  _et al._ (2022) as a function of time $t$. The
angular resonance frequency
$\omega_{0}/2\pi=f_{0}=$2.646\text{\,}\mathrm{MHz}$$, the damping rate
$\Gamma=\frac{\omega_{0}}{Q}=$12\text{\,}\mathrm{kHz}$$, and the nonlinearity
$\beta$ are approximately equal for both resonators, see section S1 in the
supplemental material Sup . Note that the capacitive diodes used to generate
the nonlinearity $\beta$ mainly possess a quadratic force component
(equivalent to the three-wave mixing enabled by $\chi^{(2)}$ in nonlinear
optics), which allows us to implement parametric pumping $\lambda$ by driving
the system with a voltage $U_{p}$ at the frequency $2f_{d}$ Nosan  _et al._
(2019). The effective model we obtain in a frame rotating at $2\pi f_{d}$ is
equivalent to that resulting from the well-known form in Eq. (1) Eichler and
Zilberberg (2023). Finally, $J$ quantifies the coupling between the resonators
($j\neq i$), and the $F_{i}(t)$ are external forces applied to the resonators
individually. In the following, we will study the competition between these
two effects ($J$ versus $F_{i}$) in ordering the phases of the two KPOs.

In Fig. 2(a), we demonstrate a frequency sweep by slowly increasing $f_{d}$ at
a constant pump amplitude $U_{p}$ with $F_{i}=0$. We use lock-in amplifiers
(Zurich Instruments MFLI) to measure the amplitudes
$A_{i}=(u_{i}^{2}+v_{i}^{2})^{1/2}$ and phases $\phi_{i}=\arctan(v_{i}/u_{i})$
of the two resonators, where $x=u\cos(2\pi f_{d}t)-v\sin(2\pi f_{d}t)$. We
observe that both resonators jump from zero to a finite amplitude around
$2.63\text{\,}\mathrm{MHz}$. The phase of the two resonators assumes a well-
defined value after the jump, with $\phi_{1}-\phi_{2}=\pi$. In our spin
analogy, this corresponds to an antisymmetric ordering (“up-down”). Around
$2.645\text{\,}\mathrm{MHz}$, the amplitudes jump again, accompanied by a
$\pi$-shift of $\phi_{1}$. The resulting state is analogous to a symmetric
spin state (“down-down”). Roughly at $2.66\text{\,}\mathrm{MHz}$, the
amplitudes drop to zero and the phases become random again.

We repeat frequency sweeps at different values of $U_{p}$. In Fig. 2(b), the
measured $\phi_{2}$ as a function of $f_{d}$ and $U_{d}$ is plotted. We
clearly see that in each sweep, the phase assumes a well-defined value in a
certain frequency range described by two overlapping, rounded triangles often
called “Arnold tongues” Heugel  _et al._ (2022); Eichler and Zilberberg
(2023). Importantly, the phase randomly assumes one of two values separated by
$\pi$ in each sweep. This randomness is a consequence of the symmetry
illustrated in Fig. 1(b), and the resulting spontaneous time-translation
symmetry breaking at each jump from zero to finite amplitude. It is a
fundamental feature of a KPO, and of networks thereof, in the absence of an
external force.

To map the different symmetry phases of the system, we employ the symmetric
and antisymmetric quadratures $v_{S}=\frac{1}{\sqrt{2}}(v_{1}+v_{2})$ and
$v_{A}=\frac{1}{\sqrt{2}}(v_{1}-v_{2})$. Note that $u_{S,A}$ can be defined
analogously and yield qualitatively similar results. In Fig. 2(c), we
represent the total amplitude $X=(v_{S}^{2}+v_{A}^{2})^{1/2}$ as a brightness
contrast, while the relative phase of the two resonators is shown in color
code. As in a previous work Heugel  _et al._ (2022), we observe $S$, $A$ and
$M$ phases inside the overlapping Arnold tongues. The center frequencies and
phase symmetries of the two Arnold tongues are inherited from the normal modes
of the system in the linear regime. The precise shapes of these zones are
understood to originate from an interplay of the nonlinearity $\beta$ and the
coupling $J$, and can be reproduced precisely with numerical simulations and
with an analytical solver Košata  _et al._ (2022); Sup .

In a next step, we additionally apply an external force to each resonator. For
simplicity, we select $F_{i}(t)=F\cos(2\pi f_{d}t+\theta)$ for both
resonators, with $\theta$ a global phase. In Fig. 2(d), we see a striking
change in the measured $\phi_{2}$ compared to Fig. 2(b): the phase now
exhibits a well-defined value over the entire parameter space, and the jump
between solutions follows a deterministic pattern. The spontaneous time-
translation symmetry breaking observed in Fig. 2(b) is entirely replaced by a
force-induced bias, as sketched in Fig. 1(c). The force therefore allows us to
control the “spin polarization” of the individual KPOs in each sweep.

Figure 3: Ringup experiments with different values of force phase $\theta$.
Both resonators are initialized in $u_{i}=0$ and $v_{i}=0$ (center of graphs).
First, only an external force $F=$20\text{\,}\mathrm{mV}$$ with phase $\theta$
is applied to break the symmetry, bringing the resonators to an amplitude of
$\approx$0.1\text{\,}\mathrm{mV}$$. Then the parametric pump
$U_{p}=$3.15\text{\,}\mathrm{V}$$ is added and the resonators ring up to
symmetric or mixed-symmetry solutions, as marked by white disks. Arrow symbols
indicate the corresponding spin analogy. All experiments were performed at
$f_{d}=$2.639\text{\,}\mathrm{MHz}$$. Figure 4: Stochastic sampling of KPO
network states. (a) Quadratures $v_{1}$ and $v_{2}$ measured as a function of
time in the presence of force noise. The applied noise is white up to
approximately $30\text{\,}\mathrm{MHz}$ and has a standard deviation of
$350\text{\,}\mathrm{mV}$. Separate noise generators were employed for the two
KPOs and connected via additional inductive lines. Only the $v_{i}$ are shown
here for simplicity. $U_{p}=$3.5\text{\,}\mathrm{V}$$, $F=0$, lock-in
demodulation rate $B=$500\text{\,}\mathrm{Hz}$$, and
$f_{d}=$2.586\text{\,}\mathrm{Hz}$$ (the resonance frequencies of both KPO
shifted due to the additional inductive lines). (b) Representation of the
timetrace from (a) in a so-called symmetry space spanned by
$v_{S}=\frac{1}{\sqrt{2}}(v_{1}+v_{2})$ and
$v_{A}=\frac{1}{\sqrt{2}}(v_{1}-v_{2})$. (c) A symmetric force of
$30\text{\,}\mathrm{mV}$ with $\theta=0$ (after subtracting a global phase
offset) favors the negative symmetric state. (d) With
$\theta=$180\text{\,}\mathrm{\SIUnitSymbolDegree}$$, the symmetric force
favors the positive symmetric state. (e) An antisymmetric force with
$\theta_{1}=0$ and $\theta_{2}=$180\text{\,}\mathrm{\SIUnitSymbolDegree}$$
favors the negative antisymmetric state. (f) inverting both phases favors the
positive antisymmetric state. A total of $60\text{\,}\mathrm{s}$ of data was
taken for each case. Insets show spin with grey dashed arrows for the local
magnetic field direction.

When plotting the system states in Fig. 2(e), we find several important
differences to the unbiased example in Fig. 2(c). First, the outlines of the
Arnold tongues are shifted by the force. This effect arises due to changes in
the stability conditions of the parametric oscillators, which was previously
studied for the case of a single KPO Ryvkine and Dykman (2006); Rhoads and
Shaw (2010); Papariello  _et al._ (2016); Leuch  _et al._ (2016); Eichler  _et
al._ (2018a); Nosan  _et al._ (2019). The most prominent manifestation of this
effect is a jump (of both phase and amplitude) at the right border of the $S$
phase. This jump is caused by the termination of the selected symmetric phase
state Sup . For a symmetric force, we mainly observe a shift of the $S$
region. Second, the antisymmetric state labelled $A$ is no longer visible, as
it is not favored by the symmetric force. Numerical simulations indicate that
the $A$ state still appears at higher values of $U_{p}$, where the impact of
the external force is reduced relative to the parametric pump. The symmetric
state $S$ therefore fills a greater portion of the diagram, in agreement with
the intuition that a symmetric force $F_{1}=F_{2}$ should favor this phase
configuration. Third, the $S$ state now only comprises an “up-up” component,
instead of allowing both “up-up” and “down-down” as in Fig. 2(c). The force
therefore fulfills the role of a potential bias, as anticipated.

In Fig. 2, we employ frequency sweeps with a fixed force phase $\theta$ to
study the solution space of our system. In many applications, however, it may
be more useful to start from the equilibrium state ($u_{i}=0$ and $v_{i}=0$),
and then directly access a particular solution through the correct choice of
$\theta$. We demonstrate this capability experimentally in Fig. 3, where the
system rings up to four different symmetric and mixed-symmetry states,
depending on the phase $\theta$ of the symmetric external force. In the spirit
of the two-spin analogy that we employ throughout the paper, we label the
final states with corresponding symbols. We find that small changes in
$\theta$ can cause the system to attain an entirely different state. For
instance, $\theta=$36\text{\,}\mathrm{\SIUnitSymbolDegree}$$ leads to a
symmetric state while $\theta=$54\text{\,}\mathrm{\SIUnitSymbolDegree}$$ leads
to a mixed-symmetry state. Rotating $\theta$ by
$180\text{\,}\mathrm{\SIUnitSymbolDegree}$ inverts the final state of both
resonators, leaving the relative configuration intact (e.g. “up-up” becomes
“down-down”).

The deterministic experiments we reported so far always reach the same stable
solutions for the same initial conditions and driving parameters. This leads
to a limited understanding of our system, which possesses several stable
solutions for certain positions in the $f_{d}$-$U_{p}$ diagram Sup .
Stochastic sampling allows us to explore these different solutions, and to
assess their dwell times in the long-time limit Margiani  _et al._ (2023). In
Fig. 4(a), we show how the system jumps between different solutions when
activated by white noise and in the absence of external forcing ($F_{i}=0$).
For most of the time, the system jumps between the two symmetric solutions,
which are more stable than the antisymmetric solutions due to their higher
amplitude. A high amplitude imposes a larger “momentum barrier” and makes the
state more stable against jumps Margiani  _et al._ (2023). In a so-called
symmetry space spanned by $v_{S}$ and $v_{A}$, this dynamics result in the
plot shown in Fig. 4(b). Here, symmetric and antisymmetric states appear at
the corners of a diamond at $v_{A}=0$ and $v_{S}=0$, respectively (plots for
the $u_{S}$ and $u_{A}$ look very similar). The edges of the diamond result
from points measured during transitions between the states. Below the graph,
we sketch the corresponding picture of two spins without an external field,
allowing both symmetric and antisymmetric states (and fluctuations between
them).

When applying an external force, we can change the relative weight (occupation
probability) of these four states. In Fig. 4(c) and (d), we apply symmetric
forces whose phase $\theta$ is tuned to favor either of the two symmetric
states. As a consequence, we observe the system only in the corresponding
symmetric state, and all jumps are suppressed. The spin picture with external
field below each graph shows how the spins are forced to align by the
homogeneous external field.

In contrast to a real (nanoscale) spin ensemble, our system allows us to
arbitrarily tune the amplitude and phase of the applied force (corresponding
to the strength and direction of the external field) for each KPO
individually. As a demonstration example, we show in Fig. 4(e) and (f) the
results for an antisymmetric force $F_{i}(t)=F\cos(2\pi f_{d}t+\theta_{i})$
with $\theta_{2}=-\theta_{1}$. Here, we find that the system occupies almost
exclusively the selected antisymmetric state. However, the suppression of the
symmetric states is weaker than in the opposite case in Fig. 4(c) and (d).
Again, this is due to the fact that the symmetric state is generally more
stable than the antisymmetric state for $J/\beta>0$ Margiani  _et al._ (2023).
Even larger forces would be necessary to entirely overcome this intrinsic
bias.

In our experiments and in the theory analysis, we find that the external force
can bias our system of coupled KPOs. In the simplest case, this bias is
analogous to an external magnetic field acting on an ensemble of coupled
spins. The tunable phase $\theta$ of the external force assumes the role of
the magnetic field angle, which can be different for each KPO, cf. Fig. 3 and
Fig. 4. This freedom in selecting the phases will be crucial in future
experiments that go beyond conventional spin systems, as it allows access to
unexplored, exotic networks that have no counterpart in solid state physics.
Such novel networks include, for instance, the Ising chain in the presence of
a tunable local impurity Falk (1966), the random-field Ising model Grinstein
(1976); Belanger and Young (1991); Bingham  _et al._ (2021); Yao and Jack
(2023) and corresponding avalanche models Perković  _et al._ (1995); Im  _et
al._ (2009); Field  _et al._ (1995); Lahini  _et al._ (2017), and the magnetic
Bose polaron Mistakidis  _et al._ (2022). We believe that controlled
experimental realizations of these elusive phenomena will spur new
developments in theory in many directions. At the same time, our work
demonstrates new strategies to control Ising machines, and to program such
systems to solve complex optimization tasks Mohseni  _et al._ (2022); Lucas
(2014); Inagaki  _et al._ (2016a); Goto  _et al._ (2019); Nigg  _et al._
(2017).

##  Appendix A Setup Calibration

Before starting any experiments with the two coupled resonators, we calibrate
their characteristics individually. The resonance frequencies of our circuits
can be tuned via the bias voltages applied to the diodes across
$47\text{\,}\mathrm{M\SIUnitSymbolOhm}$ resistors that act as a low-pass
filters Nosan  _et al._ (2019). We can apply a voltage to one resonator while
keeping the other resonator at zero bias, which results in a large detuning
between the resonance frequencies. Driving the biased resonator with a small
external drive enables us to extract its resonance frequency, while a sweep
with a large drive allows us to assess its nonlinear coefficient. This
procedure is then repeated with the second resonator after tuning the bias
voltage until the resonance frequencies match within the experimental
uncertainty, which is on the order of $100\text{\,}\mathrm{Hz}$. We noted,
however, that the device characteristics can change slightly when they are
resonantly coupled. For this reason, we extract the precise parameters of the
coupled system from the measured data, as described below.

##  Appendix B Numerical Simulations

The boundary of the large-amplitude region in Fig. 2(c) corresponds to two
overlapping Arnold tongues Heugel  _et al._ (2022); Eichler and Zilberberg
(2023). As an initial guess, we assume the two resonators to be identical. By
fitting a line to this boundary, we can extract the parametric threshold
$U_{\mathrm{th}}=$1.78\text{\,}\mathrm{V}$$ (the lowest tip of the tongues)
and the quality factor $Q=$227$$, as well as a resonance frequency
$f_{0}=$2.646\text{\,}\mathrm{MHz}$$ (the middle between the two tongues) and
frequency splitting $\Delta f=$9.4\text{\,}\mathrm{kHz}$$ (the separation
between the tongues). The coupling constant $J$ can then be calculated as
$J=(2\pi)^{2}f_{0}\Delta f$, while the threshold $U_{\mathrm{th}}$ and the
quality factor $Q$ allow us to transform the parametric drive amplitudes
$U_{d}$ to unit-less parametric modulation depths
$\lambda=\frac{2U_{d}}{QU_{\mathrm{th}}}$.

Figure 5: Amplitude and phase configuration (symmetry) of the two resonators
calculated with a numerical simulation tool Eichler and Zilberberg (2023). The
results were matched to the measurements shown in Fig. 2 of the main text by
fine-tuning the parameters. We obtained
$\omega_{1}/2\pi=$2.647\text{\,}\mathrm{MHz}$$,
$\omega_{2}/2\pi=$2.645\text{\,}\mathrm{MHz}$$,
$J=$0.98\text{\,}{\mathrm{MHz}}^{2}$$, $Q_{1}=234$, $Q_{2}=219$,
$\beta_{1}=$-8.4\text{\times}{10}^{17}\text{\,}{\mathrm{V}}^{-2}\text{\,}{\mathrm{Hz}}^{-2}$$,
and
$\beta_{2}=$-6.8\text{\times}{10}^{17}\text{\,}{\mathrm{V}}^{-2}\text{\,}{\mathrm{Hz}}^{-2}$$.
(a) Simulation result for $F=0$. The black line is identical with the one in
Fig. 2(c). (b) We reproduce the result in Fig. 2(e) with an effective force
$F=$1.6\text{\times}{10}^{08}\text{\,}\mathrm{V}\text{\,}{\mathrm{Hz}}^{-2}$$
and a phase $\phi=$-40\text{\,}\mathrm{\SIUnitSymbolDegree}$$. As for the
experiment, the theoretically expected Arnold tongue boundaries (dashed lines)
match in (a) but not in (b) due to the shifted boundaries.

In a second step, we use these parameters to perform numerical simulations of
frequency sweeps for different modulation depths $\lambda$. We find that the
phase transitions between antisymmetric, symmetric and mixed dynamical phases
($A$, $S$ and $M$, respectively) are highly sensitive to non-degeneracies
bewteen the resonator parameters. We therefore recursively modify the
individual resonance frequencies, quality factors, and nonlinear coefficients
slightly to optimize the agreement with the measured results. This procedure
yielded the parameters used for Fig. 5(a).

Finally, an external force term is added to the simulations. The applied force
strength and phase depends on the induction parameters between the driving
coils and the resonators coils, while the measurement gain and phase depend on
the induction parameters between the resonator coils and the pickup coils. As
these two sets of parameters could not be calibrated independently, we decided
to determine the effective strength and phase of the drive $F$ by fine-tuning
until the observed features are comparable to the measured results. This leads
us to the values used in Fig. 5(b).

##  Appendix C Number and stability of Solutions

For a numerical analysis of the KPO network, we employ the open-source
software package HarmonicBalance.jl Košata  _et al._ (2022). The package
allows us to transform the coupled equations from Eq. (1) into slow-flow
equations for the quadratures $u_{i}$ and $v_{i}$:

$\displaystyle\dot{u}_{1,2}^{\phantom{\dagger}}$ $\displaystyle=-\frac{\Gamma
u_{1,2}^{\phantom{\dagger}}}{2}-\left(\frac{3\beta}{8\omega}X_{1,2}^{2}+\frac{\omega_{0}^{2}-\omega^{2}}{2\omega}+\frac{\lambda\omega_{0}^{2}}{4\omega}\right)v_{1,2}^{\phantom{\dagger}}$
$\displaystyle+\frac{Jv_{2,1}^{\phantom{\dagger}}}{2\omega}+\frac{F_{u,1,2}}{m}\,,$
$\displaystyle\dot{v}_{1,2}^{\phantom{\dagger}}$ $\displaystyle=-\frac{\Gamma
v_{1,2}^{\phantom{\dagger}}}{2}+\left(\frac{3\beta}{8\omega}X_{1,2}^{2}+\frac{\omega_{0}^{2}-\omega^{2}}{2\omega}-\frac{\lambda\omega_{0}^{2}}{4\omega}\right)u_{1,2}^{\phantom{\dagger}}$
$\displaystyle-\frac{Ju_{2,1}^{\phantom{\dagger}}}{2\omega}+\frac{F_{v,1,2}}{m}\,,$
(2)

with $X_{i}=(u_{i}^{2}+v_{i}^{2})^{1/2}$ the total amplitude of resonator $i$
and $F_{u,i}$, and $F_{v,i}$ the quadrature components of the force acting on
it, respectively. Our main focus lies in examining the system’s steady states,
specifically the solutions of $u_{1,2}$ and $v_{1,2}$ when
$\dot{u}_{1,2}=\dot{v}_{1,2}=0$. Therefore, the task of identifying all the
steady states reduces to determining all the roots of the polynomial system
described by Eq. (C). To accomplish this, HarmonicBalance.jl utilizes Homotopy
Continuation Breiding and Timme (2018), a technique that enables the
computation of these roots.

We can evaluate all possible roots of Eq. (C) as a function of $\omega$ and
$\lambda$, and assess their stability with standard methods Eichler and
Zilberberg (2023). Comparing these solutions as a function of $f_{d}$ to an
experimental sweep provides additional insights. An example comparison in the
case with external force is plotted in Fig. 6. There, we can follow the
measured system state (thick grey line) from low to high $f_{d}$, starting in
a state near zero amplitude. Roughly at $2.631\text{\,}\mathrm{MHz}$, the
system jumps to one out of four solutions with mixed symmetry (see arrow
labelled i). These solutions are all different for non-identical resonators
and respond sensitively to any tuning of the parameters, making it difficult
to obtain a good agreement between theory and experiment without heuristic
fine-tuning (which we avoided here).

Figure 6: Comparison of measured quadratures $u$ and $v$ for resonator 1 and 2
as a function of frequency. Experimental data are shown as thick grey lines
and are extracted from the measurement shown in Fig. 2(e) for
$U_{p}=$3.5\text{\,}\mathrm{V}$$. Stable theory solutions are shown in blue
for symmetric ($S$), red for antisymmetric (A), pink for mixed-symmetry (M),
and black for the solution near zero amplitude ($Z$). Dashed lines represent
unstable solutions. Arrow symbols indicate analogy to a two-spin system for
thy symmetric solutions. The roman numbers i to iv indicate frequencies where
the experimental trace jumps between solutions.

At $2.636\text{\,}\mathrm{MHz}$, the system jumps from the mixed-symmetry
solution to a symmetric solution (ii). From the theory analysis, we understand
that a jump is necessary because the mixed-symmetry solutions all terminate in
bifurcations where they merge with unstable solutions. In principle, the
system offers two symmetric ($S$) and two antisymmetric ($A$) stable solutions
that the system could jump to. We interpret the selected symmetric state as
dictated by the external force, as already discussed in Fig. 2. Future
analysis will be dedicated to analyzing the precise rules governing the state
selection.

At $2.66\text{\,}\mathrm{MHz}$, the system jumps a third time, with both
resonators undergoing a phase change of $\pi$ (iii). This situation is
reminiscent to what was found in a single KPO in the presence of a symmetry-
breaking force Rhoads and Shaw (2010); Papariello  _et al._ (2016); Leuch  _et
al._ (2016); Eichler  _et al._ (2018a); Nosan  _et al._ (2019). There, the
external force resulted in a local amplitude minimum associated to a
bifurcation, and to a jump in a frequency sweep. In the present case, we
observe the same phenomenon for the symmetric normal mode of the system. The
third jump at iii is therefore a direct manifestation of the interplay between
the parametric pump and the external force applied to a parametric oscillator
network.

## References

  * Ryvkine and Dykman (2006) D Ryvkine and M I Dykman, “Resonant symmetry lifting in a parametrically modulated oscillator,” Physical Review E 74, 061118 (2006).
  * Mahboob and Yamaguchi (2008) I. Mahboob and H. Yamaguchi, “Bit storage and bit flip operations in an electromechanical oscillator,” Nature Nanotechnology 3, 275–279 (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).
  * Eichler  _et al._ (2011) Alexander Eichler, Julien Chaste, Joel Moser, and Adrian Bachtold, “Parametric amplification and self-oscillation in a nanotube mechanical resonator,” Nano Letters 11, 2699–2703 (2011).
  * Eichler  _et al._ (2018a) Alexander Eichler, Toni L Heugel, Anina Leuch, Christian L Degen, Ramasubramanian Chitra, and Oded Zilberberg, “A parametric symmetry breaking transducer,” Applied Physics Letters 112, 233105 (2018a).
  * Gieseler  _et al._ (2012) Jan Gieseler, Bradley Deutsch, Romain Quidant, and Lukas Novotny, “Subkelvin parametric feedback cooling of a laser-trapped nanoparticle,” Phys. Rev. Lett. 109, 103603 (2012).
  * Lin  _et al._ (2014) Z.R. 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 Communications 5, 4480 (2014).
  * Puri  _et al._ (2017a) S. Puri, S. Boutin, and A. Blais, “Engineering the quantum states of light in a kerr-nonlinear resonator by two-photon driving,” npj Quantum Information 3, 18 (2017a).
  * Eichler  _et al._ (2018b) Alexander Eichler, Toni L. Heugel, Anina Leuch, Christian L. Degen, R. Chitra, and Oded Zilberberg, “A parametric symmetry breaking transducer,” Applied Physics Letters 112, 233105 (2018b), https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.5031058/14512731/233105_1_online.pdf .
  * Nosan  _et al._ (2019) Z. Nosan, P. Märki, N. Hauff, C. Knaut, and A. Eichler, “Gate-controlled phase switching in a parametron,” Phys. Rev. E 99, 062205 (2019).
  * Frimmer  _et al._ (2019) Martin Frimmer, Toni L. Heugel, Ziga Nosan, Felix Tebbenjohanns, David Hälg, Abdulkadir Akin, Christian L. Degen, Lukas Novotny, R. Chitra, Oded Zilberberg, and Alexander Eichler, “Rapid flipping of parametric phase states,” Phys. Rev. Lett. 123, 254102 (2019).
  * Grimm  _et al._ (2019) Alexander Grimm, Nicholas E. Frattini, Shruti Puri, Shantanu O. Mundhada, Steven Touzard, Mazyar Mirrahimi, Steven M. Girvin, Shyam Shankar, and Michel H. Devoret, “Stabilization and operation of a kerr-cat qubit,” Nature 584, 205–209 (2019).
  * Wang  _et al._ (2019) Zhaoyou Wang, Marek Pechal, E. Alex Wollack, Patricio Arrangoiz-Arriola, Maodong Gao, Nathan R. Lee, and Amir H. Safavi-Naeini, “Quantum dynamics of a few-photon parametric oscillator,” Phys. Rev. X 9, 021049 (2019).
  * Puri  _et al._ (2019) Shruti Puri, Alexander Grimm, Philippe Campagne-Ibarcq, Alec Eickbusch, Kyungjoo Noh, Gabrielle Roberts, Liang Jiang, Mazyar Mirrahimi, Michel H. Devoret, and S. M. Girvin, “Stabilized cat in a driven nonlinear cavity: A fault-tolerant error syndrome detector,” Phys. Rev. X 9, 041009 (2019).
  * Miller  _et al._ (2019) James M.L. Miller, Dongsuk D. Shin, Hyun-Keun Kwon, Steven W. Shaw, and Thomas W. Kenny, “Phase control of self-excited parametric resonators,” Phys. Rev. Applied 12, 044053 (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).
  * Ising (1925) Ernst Ising, “Beitrag zur theorie des ferromagnetismus,” Zeitschrift für Physik 31, 253–258 (1925).
  * Mahboob  _et al._ (2016) Imran Mahboob, Hajime Okamoto, and Hiroshi Yamaguchi, “An electromechanical ising hamiltonian,” Science Advances 2, e1600236 (2016).
  * Goto (2016) H. Goto, “Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network,” Scientific Reports 6, 21686 (2016).
  * Bello  _et al._ (2019) Leon Bello, Marcello Calvanese Strinati, Emanuele G. Dalla Torre, and Avi Pe’er, “Persistent coherent beating in coupled parametric oscillators,” Phys. Rev. Lett. 123, 083901 (2019).
  * Okawachi  _et al._ (2020) Yoshitomo Okawachi, Mengjie Yu, Jae K. Jang, Xingchen Ji, Yun Zhao, Bok Young Kim, Michal Lipson, and Alexander L. Gaeta, “Demonstration of chip-based coupled degenerate optical parametric oscillators for realizing a nanophotonic spin-glass,” Nature Communications 11, 4119 (2020).
  * Mohseni  _et al._ (2022) Naeimeh Mohseni, Peter L McMahon, and Tim Byrnes, “Ising machines as hardware solvers of combinatorial optimization problems,” Nature Reviews Physics , 1–17 (2022).
  * Lucas (2014) Andrew Lucas, “Ising formulations of many np problems,” Frontiers in Physics 2, 5 (2014).
  * Inagaki  _et al._ (2016a) Takahiro Inagaki, Yoshitaka Haribara, Koji Igarashi, Tomohiro Sonobe, Shuhei Tamate, Toshimori Honjo, Alireza Marandi, Peter L. McMahon, Takeshi Umeki, Koji Enbutsu, Osamu Tadanaga, Hirokazu Takenouchi, Kazuyuki Aihara, Ken-ichi Kawarabayashi, Kyo Inoue, Shoko Utsunomiya, and Hiroki Takesue, “A coherent ising machine for 2000-node optimization problems,” Science 354, 603–606 (2016a).
  * Goto  _et al._ (2019) Hayato Goto, Kosuke Tatsumura, and Alexander R Dixon, “Combinatorial optimization by simulating adiabatic bifurcations in nonlinear hamiltonian systems,” Science advances 5, eaav2372 (2019).
  * Nigg  _et al._ (2017) Simon E. Nigg, Niels Lörch, and Rakesh P. Tiwari, “Robust quantum optimizer with full connectivity,” Science Advances 3, e1602273 (2017).
  * Heugel  _et al._ (2019) Toni L. Heugel, Matthias Oscity, Alexander Eichler, Oded Zilberberg, and R. Chitra, “Classical many-body time crystals,” Phys. Rev. Lett. 123, 124301 (2019).
  * Heugel  _et al._ (2022) Toni L. Heugel, Oded Zilberberg, Christian Marty, R. Chitra, and Alexander Eichler, “Ising machines with strong bilinear coupling,” Phys. Rev. Research 4, 013149 (2022).
  * Margiani  _et al._ (2023) Gabriel Margiani, Javier del Pino, Toni L. Heugel, Nicholas E. Bousse, Sebastián Guerrero, Thomas W. Kenny, Oded Zilberberg, Deividas Sabonis, and Alexander Eichler, “Deterministic and stochastic sampling of two coupled kerr parametric oscillators,” Phys. Rev. Res. 5, L012029 (2023).
  * Puri  _et al._ (2017b) Shruti Puri, Christian Kraglund Andersen, Arne L. Grimsmo, and Alexandre Blais, “Quantum annealing with all-to-all connected nonlinear oscillators,” Nature Communications 8, 15785 (2017b).
  * Goto  _et al._ (2018) H. Goto, Z. Lin, and Y. Nakamura, “Boltzmann sampling from the ising model using quantum heating of coupled nonlinear oscillators,” Scientific Reports 8, 7154 (2018).
  * Dykman  _et al._ (2018) M. I. Dykman, Christoph Bruder, Niels Lörch, and Yaxing Zhang, “Interaction-induced time-symmetry breaking in driven quantum oscillators,” Phys. Rev. B 98, 195444 (2018).
  * Rota  _et al._ (2019) Riccardo Rota, Fabrizio Minganti, Cristiano Ciuti, and Vincenzo Savona, “Quantum critical regime in a quadratically driven nonlinear photonic lattice,” Phys. Rev. Lett. 122, 110405 (2019).
  * Calvanese Strinati  _et al._ (2019) Marcello Calvanese Strinati, Leon Bello, Avi Pe’er, and Emanuele G. Dalla Torre, “Theory of coupled parametric oscillators beyond coupled ising spins,” Phys. Rev. A 100, 023835 (2019).
  * Takesue  _et al._ (2020) Hiroki Takesue, Kensuke Inaba, Takahiro Inagaki, Takuya Ikuta, Yasuhiro Yamada, Toshimori Honjo, Takushi Kazama, Koji Enbutsu, Takeshi Umeki, and Ryoichi Kasahara, “Simulating ising spins in external magnetic fields with a network of degenerate optical parametric oscillators,” Phys. Rev. Appl. 13, 054059 (2020).
  * Inagaki  _et al._ (2016b) T. Inagaki, K Inaba, R. Hamerly, K Inoue, Y Yamamoto, and H. Takesue, “Large-scale ising spin network based on degenerate optical parametric oscillators,” Nature Photonics 10, 415–420 (2016b).
  * Gershenzon  _et al._ (2020) Igor Gershenzon, Geva Arwas, Sagie Gadasi, Chene Tradonsky, Asher Friesem, Oren Raz, and Nir Davidson, “Exact mapping between a laser network loss rate and the classical xy hamiltonian by laser loss control,” Nanophotonics 9, 4117–4126 (2020).
  * (38) “See supplemental material at [url will be inserted by publisher] for additional setup descriptions, more measurement data.” .
  * Eichler and Zilberberg (2023) Alexander Eichler and Oded Zilberberg,  _Classical and Quantum Parametric Phenomena_ (Oxford University Press, 2023).
  * Košata  _et al._ (2022) Jan Košata, Javier del Pino, Toni L. Heugel, and Oded Zilberberg, “HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance,” SciPost Phys. Codebases , 6 (2022).
  * Rhoads and Shaw (2010) Jeffrey F. Rhoads and Steven W. Shaw, “The impact of nonlinearity on degenerate parametric amplifiers,” Applied Physics Letters 96, 234101 (2010), https://doi.org/10.1063/1.3446851 .
  * Papariello  _et al._ (2016) Luca Papariello, Oded Zilberberg, Alexander Eichler, and R. Chitra, “Ultrasensitive hysteretic force sensing with parametric nonlinear oscillators,” Phys. Rev. E 94, 022201 (2016).
  * Leuch  _et al._ (2016) Anina Leuch, Luca Papariello, Oded Zilberberg, Christian L. Degen, R. Chitra, and Alexander Eichler, “Parametric symmetry breaking in a nonlinear resonator,” Phys. Rev. Lett. 117, 214101 (2016).
  * Falk (1966) H. Falk, “Ising chain with a spin impurity,” Phys. Rev. 151, 304–311 (1966).
  * Grinstein (1976) G. Grinstein, “Ferromagnetic phase transitions in random fields: The breakdown of scaling laws,” Phys. Rev. Lett. 37, 944–947 (1976).
  * Belanger and Young (1991) DP Belanger and AP Young, “The random field ising model,” Journal of magnetism and magnetic materials 100, 272–291 (1991).
  * Bingham  _et al._ (2021) N. S. Bingham, S. Rooke, J. Park, A. Simon, W. Zhu, X. Zhang, J. Batley, J. D. Watts, C. Leighton, K. A. Dahmen, and P. Schiffer, “Experimental realization of the 1d random field ising model,” Phys. Rev. Lett. 127, 207203 (2021).
  * Yao and Jack (2023) Liheng Yao and Robert L Jack, “Thermal vestiges of avalanches in the driven random field ising model,” Journal of Statistical Mechanics: Theory and Experiment 2023, 023303 (2023).
  * Perković  _et al._ (1995) Olga Perković, Karin Dahmen, and James P. Sethna, “Avalanches, barkhausen noise, and plain old criticality,” Phys. Rev. Lett. 75, 4528–4531 (1995).
  * Im  _et al._ (2009) Mi-Young Im, Peter Fischer, Dong-Hyun Kim, and Sung-Chul Shin, “Direct observation of individual barkhausen avalanches in nucleation-mediated magnetization reversal processes,” Applied Physics Letters 95 (2009).
  * Field  _et al._ (1995) Stuart Field, Jeff Witt, Franco Nori, and Xinsheng Ling, “Superconducting vortex avalanches,” Phys. Rev. Lett. 74, 1206–1209 (1995).
  * Lahini  _et al._ (2017) Yoav Lahini, Omer Gottesman, Ariel Amir, and Shmuel M. Rubinstein, “Nonmonotonic aging and memory retention in disordered mechanical systems,” Phys. Rev. Lett. 118, 085501 (2017).
  * Mistakidis  _et al._ (2022) SI Mistakidis, GM Koutentakis, F Grusdt, P Schmelcher, and HR Sadeghpour, “Inducing spin-order with an impurity: phase diagram of the magnetic bose polaron,” New Journal of Physics 24, 083030 (2022).
  * Breiding and Timme (2018) Paul Breiding and Sascha Timme, “HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia,” in  _International Congress on Mathematical Software_ (Springer, 2018) pp. 458–465.