Deterministic and stochastic sampling of two coupled Kerr parametric oscillators

The vision of building computational hardware for problem optimization has spurred large efforts in the physics community. In particular, networks of Kerr parametric oscillators (KPOs) are envisioned as simulators for finding the ground states of Ising Hamiltonians. It was shown, however, that KPO networks can feature large numbers of unexpected solutions that are difficult to sample with the existing deterministic (i.e., adiabatic) protocols. In this work, we experimentally realize a system of two classical coupled KPOs, and we find good agreement with the predicted mapping to Ising states. We then introduce a protocol based on stochastic sampling of the system, and we show how the resulting probability distribution can be used to identify the ground state of the corresponding Ising Hamiltonian. This method is akin to a Monte Carlo sampling of multiple out-of-equilibrium stationary states and is less prone to become trapped in local minima than deterministic protocols.

The Kerr parametric oscillator (KPO) is a nonlinear system whose potential energy is modulated at a frequency f p close to twice its resonance frequency, f p ≈ 2f 0 [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. When the modulation depth λ exceeds a threshold λ th , the system features two stationary oscillation solutions. The solutions have a frequency f p /2, an amplitude X determined by λ relative to the Kerr nonlinearity, and phases that differ by π. These 'phase states' can be mapped to the two states σ ∈ {−1, 1} of an Ising spin. Building on that analogy, it was proposed that networks of KPOs can be utilized to simulate the ground state of coupled spin ensembles, as captured by the Ising model Hamiltonian [17]: where K i,j is the coupling coefficient between two spins with states σ i,j . Interestingly, finding this ground state is equivalent to many computational problems that are nearly intractable with conventional computers [18], such as the number partitioning problem [19], the MAX-CUT problem [20,21], and the famous traveling salesman problem [22].
Various physical implementations have been proposed or realized as "Ising solvers" [14,[23][24][25][26][27][28][29]. A well known example is the Coherent Ising Machine, a network of degenerate optical parametric oscillators (DOPOs) that are coupled through a programmable electronic feedback element [20,25,28,[30][31][32] (note that DOPOs differ from KPOs in that their amplitude is not determined by their Kerr nonlinearity but rather by two-photon loss). The feedback breaks the energy conservation of the network and imparts dissipative coupling (mutual damping) between the oscillators. As a consequence, different network configurations corresponding to different Ising solutions become stable at different driving thresholds. The opti-mal solution is assumed to possess the lowest threshold and therefore to appear as the solution when driving the network.
A different line of investigation focuses on energyconserving, bilinear coupling between KPOs [19,27,[33][34][35][36][37]. In this works, the coupled oscillator solutions can be approximated as decoupled normal modes with split resonance frequencies. In contrast to the case of dissipative coupling, the lowest threshold that is encountered when ramping up the driving strength depends here on the detuning ∆ between the external drive and the resonance frequency. This control parameter opens up the possibility for specific protocols to find different solutions [19,27,33]. It was experimentally demonstrated, however, that the combination of nonlinearities and strong bilinear coupling can give rise to a rich solution space beyond that of a simple Ising model [37]. Careful validation and testing of small systems is therefore important before larger networks can be understood and operated correctly.
In this paper, we experimentally test the validity of the Ising analogy for a system of two classical coupled KPOs. In a first step, we apply an adiabatic ramping protocol to find one particular solution for each selected combination of ∆ and λ in a deterministic way. In a second step, we use strong force noise to explore the solution space of the system: this method is based on transitions between different stationary KPO solutions [38][39][40][41][42], and it allows for the visualization of a probability distribution for all accessible states. Such "stochastic sampling" presents a way of quantifying the occupation probability of each solution, and thus selecting the optimal state. Surprisingly, we find that the oscillation state corresponding to the expected Ising ground state has not the highest but the lowest occupation probability over the entire parameter space. We reconcile this result with theory and predict 0 20  how the method can be used for larger networks.
Our system comprises two microelectromechanical resonators (MEMS) made from highly-doped single-crystal silicon. Both resonators are fabricated on the same chip in a wafer-scale encapsulation process [43], and they have the shape of double-ended tuning forks with branches 200 µm long and 6 µm thick; see Fig. 1(a) and Appendix A. They have resonance frequencies of roughly f 0 ≈ 1.124 MHz and quality factors of Q ≈ 13500. Bias voltages U b can be used to fine-tune the resonator frequencies by a few kHz and induce negative Kerrnonlinear coefficients of β ≈ −63.2 × 10 17 V −2 s −2 due to the nonlinear electrostatic forces between the biased tuning fork and the electrodes next to it [44]. Those electrodes capacitively transduce the motion into elec-trical signals that are measured with a Zurich Instruments HF2LI lock-in amplifier. The capacitive driving and measurement allows us to write effective equations of motion [10] as where ω 0 /2π = f 0 , γ = ω0 Q , x i is the measured voltage signal of resonator i, J quantifies the coupling to resonator j, and U ξ,i indicate uncorrelated white noise sources. The electrical tuning effect allows us to parametrically modulate (drive) the resonator potentials at frequency f p [42]. The required oscillating driving voltage is U d = λU 0 th Q/2, where U 0 th ≈ 2.4 V is the measured parametric threshold voltage on resonance. As a function of detuning ∆ = f p /2 − f 0 , the driving threshold for parametric oscillation, U th , is described by a socalled 'Arnold tongue', see Fig. 1(b). Outside the Arnold tongue, a resonator is stable at zero amplitude, while inside the tongue the zero-amplitude solution becomes unstable and the resonator oscillates at f p /2 with a finite effective amplitude X in one out of two possible phase states.
The resonators are mechanically coupled via their common substrate [45]. We calibrate the coupling strength from the normal-mode splitting and find ∆ f = (−2.6 ± 0.3) Hz, corresponding to a coupling coefficient J = 4π 2 ∆ f f 0 = −(113 ± 13) × 10 6 Hz 2 between the two resonators; see Appendix A for details. Even though the coupling is weak, |∆ f | f 0 /Q, we can use a normalmode basis of antisymmetric and symmetric oscillations to describe our system in the following.
When both resonators are operated as KPOs with a parametric drive voltage U d > U th , each of them selects one of its two phase states. The resonators can respond either in the same (symmetric) or in opposite (antisymmetric) phases. Which of those two solutions is preferred depends on the signs of J and ∆. In Figs. 1(c) and 1(d), we show the experimental results of sweeping the driving frequency from negative to positive ∆ at a fixed U d . The system first rings up into the antisymmetric state at ∆ = −35 Hz before it flips to a symmetric configuration close to ∆ = 25 Hz.
The ordering observed in Figs. 1(c) and 1(d) reflects the fact that for J < 0, the antisymmetric normal mode has a lower eigenfrequency than the symmetric mode. We can therefore expect to find separate normal-mode Arnold tongues for symmetric and antisymmetric oscillations with a splitting in frequency; see Fig. 2(a) [37,46]. This splitting has important consequences for the driven system: when the drive is suddenly activated at a specific detuning ∆ above threshold, the KPOs should preferentially select the normal-mode oscillation state with the lowest threshold. We experimentally test this prediction in Fig. 2(b) by measuring the chosen state once for each pixel individually, and we find very good agreement with the schematic in Fig. 2 The tongue for the antisymmetric and symmetric normal mode are shown in orange and purple, respectively. ∆ = 0 is shown as a dashed grey line and arrows indicate the solution with the lowest threshold for positive and negative ∆. (b) Measured response to a parametric drive initialized with fixed strength U d and detuning ∆. For each drive-frequency pair, the parametric drive was activated and the system was given approx. 400 ms to settle into a solution before the result was recorded with an integration time of 16 ms. Afterwards, the drive was switched off for 400 ms before moving on to the next coordinate. Symmetric (purple) and antisymmetric (orange) solutions were identified by comparing the measured phases of the two KPOs, while the shape of the Arnold tongue was extracted from their amplitudes. (c) Measured response to a parametric drive initialized with detuning ∆ and a strength that is slowly ramped upwards from U d = 0 with no initial oscillation (white). The color coding is the same as in (b).
normal-mode splitting. Figure 2(a) offers a simple interpretation of several recent theory proposals for simulating the Ising ground state with KPO networks [19,27,33]. The proposals predict that for quasiadiabatic ramping of U d with ∆ < 0, the stable solution emerging beyond the lowest instability threshold corresponds to the correct Ising ground-state solution. As we see in Fig. 2(a), such a protocol can effectively be reduced to finding the normal mode with the lowest eigenfrequency in the case of two coupled KPOs. We confirm the protocol in Fig. 2(c). Our system rings up to symmetric and antisymmetric states for ∆ > 0 and ∆ < 0, respectively. As the Ising ground state we seek is the antiferromagnetic one, the latter presents the correct solution. Note that the protocol would still work for J > 0, as both the Ising ground state and the normalmode ordering would then be reversed. For strong coupling and larger networks, additional nonlinear effects can make the solution space more complex to map [37]. This will be the scope of future experiments.
The statistical spread of states in Fig. 2(b) can be understood from the competition between deterministic ordering due to the coupling term J on the one hand, and stochastic (thermal) force noise on the other hand. When the drive is switched on suddenly, noise in a wide frequency range participates in this competition. Note, however, that the noise intensity is always finite, resulting in a finite statistical bias towards the solution preferred by the coupled system. Slow ramping of the drive additionally low-pass filters the noise and favors a deterministic ordering. This is why the outcome in Fig. 2(c) is neatly divided into two halves.
Instead of suppressing the noise, it can be interesting to enhance it in order to enable activated jumps between all states. In this way, we gain a "stochastic sampling" map of the solution space at particular values of U d and ∆, rather than just a single solution [47]. Activated escape involves a random walk from the initial state (first quasipotential well) over a quasienergy barrier, and a deterministic decay into the opposite state (second quasipotential well) [48,49]. In analogy to a thermodynamic system, we naively expect that the optimal solution occupies the deepest quasipotential well and therefore has the longest average dwell time τ between switches [35,42]. We test the stochastic sampling protocol in Fig. 3. White voltage noise with a standard deviation U ξ,i applied to each drive electrode enhances the force noise enough to cause activated jumps between the symmetric and antisymmetric solutions as a function of time; see Fig. 3(a) [42]. Plotting such a data set as a function of the variables u s = (u 1 + u 2 )/ √ 2 and u a = (u 1 − u 2 )/ √ 2, where we define x i (t) = u i cos(πf p t)+v i sin(πf p t), allows us to identify symmetric and antisymmetric solutions, respectively, cf. Fig. 3(b).
We repeat the stochastic sampling for different values of ∆ and U d , and we summarize the results in Fig. 4. The relative dwell times measured for the symmetric and antisymmetric states are shown in cake diagrams. Surprisingly, the driven out-of-equilibrium system favors the symmetric state in the entire range of parameters. This is in contrast to the behavior expected from the corresponding (equilibrium) Ising Hamiltonian H Ising , where the antisymmetric ground state would always have the highest occupation.
To investigate this enigma, we compare the switching paths that can carry the system from the symmet- (a) Extract from a measured sample path of both KPOs under the influence of applied force noise. Switches between the quasistable states are visible as jumps from −25 mV to +25 mV and vice versa. Only the u quadratures and their difference are shown for simplicity. A high difference indicates that the system is in an antisymmetric state (orange dots), while differences close to zero correspond to a symmetric state (purple dots). (b) Representation of the entire sample path in a "super phase space" spanned by us = (u1 + u2)/ √ 2 and ua = (u1 − u2)/ √ 2. Orange and purple dots mark the antisymmetric (us ≈ 0) and symmetric (ua ≈ 0) phase state configurations, respectively. The graph also shows how switches mostly involve a change in symmetry from symmetric to antisymmetric or vice versa, while symmetry-preserving transitions through the origin of the plot are rare. We use U d = 3 V and ∆ = 0. The total measurement time was 6 min, measured with an integration time of 138 µs at 3597 samples per second.
ric to the antisymmetric state or vice versa. For small switching probabilities, the main contribution of the random walk is concentrated in a narrow channel in phase space [40], such that we can approximate the total switching rate from the optimal switching path Y min for each transition (see Appendix C).
The theory results obtained for our device parameters agree with our experimental findings. As shown in Fig. 4(c), the results consistently predict longer dwell times for the symmetric configuration than for the antisymmetric one over the entire range of ∆, meaning that the time required to escape the quasipotential well of the symmetric state is longer on average. The ordering of the dwell time coincides with that of the state amplitudes shown in Fig. 4(b). Interestingly, when repeating the calculations for a positive nonlinearity β, the result changes and the antisymmetric state is found to be more stable in general; see Fig. 4(d). Note that switching the sign of β changes the ordering of the amplitudes of the normal modes in Fig. 4(b) as well. We conclude that for the system we study, the most stable state is always the one with the largest amplitude. We can understand the experimental and theoretical results in a straightforward way. To switch from one normal mode to another, one of the participating resonators must switch to the opposite phase state, which can be achieved without energy transfer [11,50] but implies a momentum reversal. With all other parameters being approximately equal, a larger normal-mode oscillation amplitude corresponds to a larger momentum to be overcome by the stochastic process. For this reason, the system typically remains trapped for a longer time in the solution with the highest amplitude, as seen in Fig. 4.
Our experimental confirmation of Ising simulation protocols [19,27,33] is a first step towards understanding large systems. For N > 2, the number of normal modes does not match the expected size of the Ising solution space, and it will be important to understand how the system evolves far beyond the parametric threshold as a function of detuning [37] and in the presence of persistent beating between solutions [51]. Mapping all solutions of a complex system can be difficult with deterministic drives due to hysteresis. Stochastic sampling, as demonstrated here, may be a way to overcome this limitation, providing a direct tomography of the solution space for a given parameter set. More generally, stochastic sampling or related experiments such as simulated annealing can be useful for understanding the analogy between an out-ofequilibrium nonlinear resonator network and thermodynamic systems. As we show in this work, the connection between the two paradigms can be counterintuitive.
Finally, we include a brief outlook on performing Ising simulation protocols with quantum-coherent systems. Quantum systems are predicted to be more efficient in finding the Ising ground state than the corresponding classical system [27], which makes them a valuable resource for solving optimization problems [18]. There, the competition between the timescales set by decoherence on the one hand and energy exchange between the KPOs on the other hand implies that strong coupling is a crucial requirement for quantum adiabatic evolution. However, the combination of strong coupling and nonlinearity was shown to impact the Ising solution space in surprising ways, calling for careful calibration [37,47]. For this reason, the development of quantum systems will benefit from methods such as stochastic sampling that allow visualizing the complete solution space.  Detailed schematics showing the two tuning-fork resonators and the electrodes for capacitive driving are shown in Fig. 5.
The resonators' properties were extracted from fits to their driven response. Applying a small driving force to only one of the resonators allows us to extract the precise resonance frequency as well as the quality factor. The same measurement can give a clear indication of coupling between the resonators when looking at the response of the nondriven resonator, cf. Fig. 6. Due to the frequencydependent driving amplitude and phase that the second resonator experiences, its amplitude response is narrower than usual and its phase changes by 180°when crossing the resonance.
The parametric threshold was found by measuring a full Arnold tongue, whose tip (lowest point) corresponds to the effective threshold on resonance U 0 th while the outer shape is determined by cf. Fig. 1(b) [52]. Equation (A1) is a reformulation of the formula for f ± shown in the caption to Fig. 2(a). The parametric modulation depth is calibrated as λ = 2U d QU 0 th . The amplitudes in a single line above threshold can then be used to determine the Duffing nonlinearity factor by fitting the parametric frequency response amplitude: where ω = πf p , b = ω 2 ω 2 0 and c = −λ 2 4 + γ 2 ω 2 ω 4 0 , cf. Fig. 1(c). Finally, we extract the coupling strength from the data in Fig. 1(d). With the frequency difference ∆ f between the end of the antisymmetric response and the end of the symmetric response, the coupling coefficient can be calculated as J = 2π∆ f ω 0 . The sign of the coupling then follows from the ordering of symmetric and antisymmetric responses. A comparable result can be extracted from the resonant amplitudes in Fig. 6. For X d and X n being the amplitudes at resonance of the driven and coupled, nondriven resonator respectively, we find J = − Xn

Appendix B: Stochastic Sampling
To test the stochastic sampling protocol as shown in Fig. 4, we measured a set of time traces for different frequencies ∆ and parametric driving strengths U d . Each measured time trace is 360 s long with 3597 samples per second. We heuristically adjusted the standard deviation σ ξ of the white noise U ξ,i applied to the resonators together with the drive as σ ξ ≈ 4 × 10 4 × (1.57X − 0.01) to reach a reasonable switching rate over the parameter range of the measurement.
In a first analysis step, each data point from a single KPO was assigned to one of the phase states labeled 0 or π. In order to do so, we defined circles in phase space whose center point were the stable attractors and whose radii were 0.6 times the state amplitude. A data point was assigned to a phase state when it was within the respective circle. If the data point was outside both circles, it was assigned to the same state as the previous data point. See Ref. [42] for details.
In a second step, we then compared the relative states of both KPOs, counting how often the two phase states are equal or opposite. The results of those polls are shown in the cake diagrams in Fig. 4. They give a qualitative estimate of which state the system prefers to be in.

Appendix C: Onsager-Machlup Function
The switching process between stable oscillation states induced by weak-noise can be described analogously to noise-activated jumping over a barrier W in equilibrium systems [53,54]. However, the barrier W between two quasistable solutions of a driven system resides in a quasipotential structure in a rotating frame. The Onsager-Machlup formalism can be used to obtain an estimate for the barrier W [55,56]. The Onsager-Machlup action is defined by: where t i (t f ) is the initial (final) time of the trajectory of the N resonator system, Y(t) = (u 1 , v 1 , ..., u N , v N ), which obeys the equation of motionẎ = f (Y). For the switching rate Γ in the weak-noise limit one can derive the scaling Γ ∝ exp(−2W/σ 2 ) with noise variance σ 2 and barrier W = S OM [Y min ], where Y min is the optimal transition path minimizing S OM [53,55]. We can thus conclude that the stable state with higher barrier will be more likely in an noisy environment. Phase and amplitude of both resonators while only one is driven by an external force. Left: resonator 1 is driven with U f = 20 mV and resonator 2 follows. Right: resonator 2 is driven while 1 follows. The non-driven resonator generally shows a narrower peak in its amplitude response, resulting from the combination of its own response function with the frequencydependent force amplitude it experiences from the externally driven resonator. Similarly, the phase of the non-driven resonator changes by 180°when crossing the resonance, which is the sum of its own harmonic response phase and the phase of the driven partner.