Controlled asymmetric Ising model implemented with parametric micromechanical oscillators

Asymmetric Ising model, in which coupled spins affect each other differently, plays an important role in diverse fields, from physics to biology to artificial intelligence. We show that coupled parametric oscillators provide a well-controlled and fully characterizable physical system to implement the model. Such oscillators are bistable. The coupling changes the rate of interstate switching of an oscillator depending on the state of other oscillators. Our experiment on two coupled micromechanical resonators reveals unusual features of asymmetric Ising systems, including the onset of a probability current that circulates in the stationary state. We relate the asymmetry to the exponentially strong effect of a periodic force on the switching rates of an individual parametric oscillator, which we measure. Our findings open the possibilities of constructing and exploring asymmetric Ising systems with controlled parameters and connectivity.


I. INTRODUCTION
Parametric oscillator is one of the best-known examples of a bistable system.It has two vibrational states with equal amplitudes and opposite phases [1].These states emerge when the oscillator eigenfrequency is periodically modulated.They have a period equal to twice the modulation period and can be associated with classical bits or Ising spin states, providing a basis for classical logic operations [2,3].Superpositions of the oppositephase coherent states of an oscillator can also encode a qubit [4,5].Coupled parametric oscillators can serve as Ising machines for classical and quantum annealing [6][7][8][9][10][11][12][13][14].Besides computation, various other applications of parametric oscillators have been studied, from force and mass sensing [15,16] to rare events in classical and quantum systems far from thermal equilibrium [17][18][19][20][21][22] and phase transitions into a time-symmetry-broken (timecrystal) state [23][24][25].
An important aspect of Ising systems pointed out by Hopfield [26] is the possibility to use coupled spins to model neural networks which memorize multiple patterns.This possibility has been attracting increasing interest over the years, particularly in view of the progress in machine learning [27,28].In the Hopfield model the spin coupling energy has the conventional form of J ij σ i σ j , where σ i , σ j take on values ±1 and J ij = J ji , and the network dynamics can be analyzed using the methods of statistical physics.The model is symmetric in the sense that the effect of spin i on spin j is the same as the effect of spin j on spin i.
However, most neuron networks are presumably asymmetric: neuron i can affect neuron j stronger than neuron j affects neuron i.If neurons are associated with spins, one can think formally that J ij ̸ = J ji and then the coupling may not be described by the coupling energy.The corresponding model is called an asymmetric Ising model.It has attracted much attention as one of the leading models of neural networks [29][30][31][32][33][34] and gene regulatory networks [35] and has been used to describe experiments on neurons, cf.[36,37] and references therein.
In spite of the importance of the asymmetric Ising model, there have been no studies that relate the spin coupling parameters to the parameters of the underlying system.Understanding the dynamics of this system enables one to examine to what extent the mapping on coupled spins is adequate, in the first place.Determining the relationship between the parameters of the system and the effective spins is essential for implementing and exploring asymmetric Ising models.
In the present paper we demonstrate that coupled parametric oscillators in the presence of noise provide a system that can be described by an asymmetric Ising model.The description is based on the Glauber picture [38] in which the rate of switching between the states of a spin depends on the states of the spins to which it is coupled.In the case of oscillators, the relevant quantity is the rate of switching between the period-two vibrational states of an oscillator that depends on which vibrational states are occupied by other oscillators.We describe the mapping of the oscillators on spins and independently measure the parameters of the system that enter the model.In particular, we measure an important characteristic of driven oscillators in the presence of fluctuations, the logarithmic susceptibility [39], which describes the exponentially strong effect of a periodic force on the switching rates of an individual uncoupled parametric oscillator.
The parametric oscillators we study are micro-electromechanical resonators modulated close to twice their eigenfrequencies.Such resonators enable exquisite con-trol of their eigenfrequencies and the coupling.Since the decay rates of our resonators are small, the modulation needed to excite parametric vibrations is comparatively weak, so that the vibrations are nearly sinusoidal.
With micromechnical resonators, we demonstrate that the asymmetric Ising model does not have detailed balance.An immediate consequence is the emergence of a probability current that circulates in the system in the stationary state.We measure this current for a system of two coupled non-identical parametric oscillators.The measurements are in excellent agreement with the theory.
We consider the case where the coupling of the oscillators is weak, so that each oscillator still has two stable vibrational states, and their amplitudes and phases are only weakly changed by the coupling.However, the coupling can significantly change the rates of noise-induced switching between the states.To gain an intuitive understanding, consider a Brownian particle in a symmetric double-well potential.Because of thermal fluctuations, the particle switches between the wells with the rate W ∝ exp(−∆U/k B T ), where ∆U is the barrier height and T is temperature [40].If the potential is tilted, the barrier heights are incremented by ±δU in the opposite wells, breaking the symmetry of the interwell switching rates.The rates acquire extra factors exp(±δU/k B T ).Even for a small tilt, the ratio δU/k B T can be large, for low temperatures.In that case the stationary populations of the wells become significantly different.
Consider now a set of weakly interacting particles in double-well potentials.A particle exerts force on other particles that depends on which well it occupies.This force tilts the potentials of the other particles and breaks the symmetry of the interwell switching rates, reminiscent of the effect of the spin-spin coupling in the Glauber model.The change of the switching rates of the spins is fully determined by the coupling energy, which in turn depends only on the relative spin orientations.For example, for two coupled spins, the change ∝ exp(−J 12 σ 1 σ 2 /k B T ) is the same for both of them.In other words, the two spins affect each other symmetrically.As we show, the picture extends to coupled parametric oscillators, even though there are no static doublewell potentials.However, a major difference is that the oscillators can affect each other asymmetrically.
If the oscillators are identical and the coupling is weak, the changes of the switching rates are equal within each pair of coupled oscillators.The system is mapped onto the symmetric Ising model.On the other hand, if the oscillators have different parameters, we show that the coupling-induced changes of the switching rates are different.The picture of a change in potential barriers no longer applies.Instead, the system is mapped onto the asymmetric Ising model.In our system, switching between the period-two vibrational states is activated by noise with controlled intensity, which allows us to fully characterize the switching rates. .c. Vibration amplitude of resonators 1 (red) and 2 (blue) subjected to identical parametric modulation as functions of ωp/2.There is no coupling between the resonators and the eigenfrequencies are tuned to be almost equal.

II. RESULTS
We present experimental results for a system of two micromechanical torsional resonators.They are shown in Fig. 1a.Each resonator consists of a movable polysilicon top plate (200 µm × 200 µm × 3.5 µm) supported by two torsional rods, with two fixed electrodes underneath.The resonators are located side by side.Their vibrations can be excited and detected independently.For resonator i (i = 1, 2), dc voltages V dc L,i , V dc R,i and V top i are applied to the left electrode, the right electrode and the top plate respectively.Application of an ac voltage on the left electrode generates a periodic electrostatic torque that excites vibrations of the top plate.The vibrations are detected by measuring the current flowing out of the top plate induced by the capacitance change between the plate and the two underlying electrodes.
In this study, only the fundamental modes of torsional vibrations are used.The eigenfrequencies of the resonators are almost identical, with ω 1 /2π ≈ 15860.562Hz for resonator 1 and ω 2 /2π ≈ 15860.598Hz for resonator 2; they can be fine tuned by adjusting dc potential difference ∆V R,i = V top i − V dc R,i between the plate and the corresponding right electrode (Supplementary Note 1).The damping constants are Γ 1 /2π ≈ 0.064 Hz and Γ 2 /2π ≈ 0.063 Hz for resonators 1 and 2 respectively.
The spring constants of the both resonators are modulated electrostatically together at frequency near 2ω 1 ≈ 2ω 2 , leading to parametric excitation of the vibrations.We also inject broadband Gaussian voltage noise for each resonator that leads to occasional switching between the period-two vibrational states.
As shown in Fig. 1a, the adjacent edges of the plates form interdigitated comb-shaped electrodes to allow the plates to couple electrostatically when there is a potential difference between them.When V cpl = 0 V, we verify that there is no coupling between the plates.We keep V cpl small as we focus on the regime of the weak coupling that only weakly perturbs the dynamics in the absence of noise.
All measurements are performed at room temperature at pressure below 10 µtorr.The eigenfrequencies and the coupling between the resonators can be tuned independently (Supplementary Note 1), which is crucial for revealing the features of the asymmetric Ising model.
The equations of motion of coupled parametric oscillators have the form For our pair of torsional resonators, i = 1, 2. The coordinate q i is the rotation angle of the ith resonator, M i is its moment of inertia, γ i is the Duffing nonlinearity parameter, F p and ω p are the amplitude and frequency of the parametric modulation, respectively, and In the experiments on the effect of the coupling, D i determines the effective temperature of the noise.We set We use resonant modulation, |ω p − 2ω i | ≪ ω i , which allows us to parametrically excite vibrations even with small F p .In the absence of resonator coupling and noise the two stable vibrational states of ith resonator are where A i and φ i are the vibration amplitude and phase, and σ i = ±1.The values of A i , φ i depend on the resonator parameters; for small damping |φ i | ≪ 1.For brevity, and where it may not cause confusion, we use ↑ and ↓ for σ i = 1 and In what follows we associate the vibrational states (2) with spin states.This association is justified provided the change of these states because of coupling the oscillators to each other, i.e., the change of the amplitudes A i and phases φ i , is small.The weakness of the coupling is thus a major condition of the mapping of the system of oscillators on the system of coupled spins.
Classical and quantum noise causes transitions between the states σ i = ±1 of an isolated oscillator.By symmetry, the rates W i (σ i ) of transitions σ i → −σ i of the ith oscillator are the same for the both states.For weak noise, the transitions are rare, W i (σ i ) ≪ Γ i , and the dependence of the switching rate on the noise intensity is given by the activation law [17,19].For classical noise where R i (σ i ) = R i (−σ i ) is the effective activation energy and C i ∼ Γ i .Activated switching in single parametric oscillators has been measured in a number of systems [18,20,23].
In our experiment, the switching rate of each resonator is extracted from the Poisson distribution of the residence times (Appendix B).Due to slight difference in the damping constants, the switching rates for the two resonators are measured to be different.We verify that, when the coupling is zero, in each resonator the two coexisting states with opposite phases are equally occupied and the populations of all 4 states σ 1,2 = ±1 are equal (Supplementary Note 3).
As seen from Eqs. ( 1) and ( 2), if the noise and the coupling of the oscillators are weak, to describe the effect of the coupling to the jth oscillator on the dynamics of the ith oscillator, one can replace the coordinate of the jth oscillator q j (t) in the equation of motion of the ith oscillator (1) by q j (σ j ; t).In this approximation, the ith oscillator is driven by a force at frequency ω p /2 exerted by the oscillators to which it is coupled.The force changes when the jth oscillator switches between its vibrational states.
The effect of weak coupling can be understood if one considers the dynamics of an isolated parametric oscillator driven by a weak extra force F d cos[(ω p /2)t + ϕ d ] that mimics the force from other oscillators [24].Such force breaks the symmetry of the vibrational states σ i = ±1.A major consequence of the symmetry lifting for weak force is the change of the switching rates W (σ i ).To leading order in F d this change has been predicted [41] to be described by an increment of the activation energy that is linear in Here Ri is the value of R(σ i ) in the absence of the drive.The parameters χ i and δ i are the magnitude and phase of the logarithmic susceptibility, i.e., the susceptibility of the logarithm of the switching rate.They strongly depend on the parameters of the oscillator and the parametric modulation, but are independent of F d and ϕ d [41].
As seen from Eqs. (3) and (4), for small noise intensity even a weak drive can significantly change the switching rates.It therefore can significantly change the stationary populations of the states w st (σ i ): from the balance equation for these populations ẇ(σ A strong population change that periodically depends on phase ϕ d was indeed seen in experiments [42].However, the general effect of the linear dependence of log[W (σ i )] on the drive amplitude of a periodic force in bistable systems has not been demonstrated other than in simulations [43].This effect may be responsible for the deviation of the escape rate from the expected quadratic dependence on the drive amplitude in Josephson junctions [44].
We measure the logarithmic susceptibility of each resonator in our two-resonator system.By setting V cpl = 0 V we ensure there is no coupling between the two resonators.For each resonator, we apply a resonant drive F d cos[(ω p /2)t + ϕ d ] on top of the parametric modulation at ω p .The drive phase ϕ d is chosen to be 3.3 • so that the results can be compared to the case of coupled os-cillators when coupling is later re-introduced.Figure 2a shows the random switches of the phase of resonator 1 as a function of time at a constant F d of 1.04 × 10 −17 Nm.The ratio of populations w st (σ 1 = +1)/w st (σ 1 = −1) ≡ w st (↑)/w st (↓) is obtained by measuring the residence time in the two states σ 1 = ±1.Figure 2b shows that this ratio deviates from 1 as F d is increased.
Next, the switching rates are measured by fitting to the Poisson distribution of the residence times (Appendix B). Figure 2c shows the effect of 1/D (which mimics the inverse noise temperature) on the logarithm of the ratio of switching rates with the symmetry breaking drive turned on and off.The upper and lower branches represent decrease and increase of the activation energy respectively, corresponding to opposite signs of σ 1 in Eq. ( 4).We obtain the increment |∆R 1 | from the average of the magnitude of the slopes of the two linear fits through the origin.The linear dependence of log[W 1 (σ 1 )/ W1 ] on 1/D in Fig. 2c confirms that the effect of a weak symmetrybreaking drive is primarily a change ∆R 1 (σ 1 ) of the activation energy of interstate switching.If D is small compared to |∆R 1 |, the change of the switching rate can be substantial.As shown in Fig. 2d, |∆R 1 | is indeed linear in F d for a weak drive.The factor χ 1 cos(ϕ d + δ 1 ) for resonator 1 is given by the slope of the linear fit (solid red line).Measurements are then repeated for resonator 2 (Supplementary Note 4) to yield χ 2 cos(ϕ d + δ 2 ).
In Fig. 2d the measurements are compared with the results of simulations of the switching rate.There is excellent agreement between measurement and the general expressions (3) and (4).However, for stronger drive the dependence of log[W i (σ i )] on F d becomes nonlinear (Supplementary Note 6).

Switching rates in the system of coupled oscillators
The above results suggest that, if we now consider coupled oscillators, the rate of switching σ i → −σ i of the ith oscillator depends on the states {σ j } of the oscillators coupled to it.From Eqs. ( 2) -( 4), for weak coupling it has the form where Wi = C i exp(− Ri /D i ) is the switching rate in the absence of coupling.The change of the activation energy Equation ( 5) has the form of the expression for the switching rates of coupled Ising spins.In the standard Ising model K ij is given by the ratio of the coupling energy J ij to k B T [38].Therefore K ij = K ji .In our case, if all oscillators are identical, we also have K ij = K ji , as seen from Eq. ( 6).Therefore the system of coupled identical parametric oscillators maps onto the standard Ising model of coupled spins.
If the oscillators are different, 6), the difference originates from both the vibration amplitudes and logarithmic susceptabilities.For K ij ̸ = K ji , the system is mapped onto the asymmetric Ising model.As seen from the known expressions for the vibration amplitudes and phases as well as the logarithmic susceptibilities (cf.[41]), the difference between K ij and K ji can be already large if, for example, the oscillator eigenfrequencies are slightly different: |ω i −ω j | ≪ ω i , but the ratio |ω i − ω j |/Γ i is not small and, most importantly, the noise intensity is small.
The stationary probability distribution w st ({σ n }) is generally not known for the asymmetric Ising model.An important feature of the model is the lack of detailed balance (Appendix D).It leads to the onset of a probability current in the stationary state.An elementary transition is a flip of a single spin, with the rate that depends on other spins.The current associated with a flip of the ith spin for a given configuration of other spins {σ j̸ =i } is, For symmetric coupling, K ij = K ji , the current ( 7) is zero (Supplementary Note 7).

Measurement of asymmetric coupling constant and probability current
We demonstrate the asymmetry in the coupling coefficients and the existence of a probability current using our system of two coupled parametric oscillators (i = 1, 2).Weak coupling between the two resonators is introduced by applying V cpl = 0.3 V. We adjust ∆V R,1 and ∆V R,2 to tune the resonant frequencies to be close but non-identical, with ω 1 − ω 2 = 0.4 Hz.The two resonators are subjected to parametric modulation of the same amplitude and the same frequency ω p .As shown in Fig. 3a, when ω p is swept up, resonator 2 undergoes a subcritical bifurcation first, followed by resonator 1.The electrostatic coupling between the two plates favors the configuration where the phases of the resonators are opposite to each other.In the absence of injected noise, resonator 1 adopts a vibration phase opposite to resonator 2 as ω p is increased.Correlations in the phase were previously observed in two nanomechanical parametric resonators [16] undergoing supercritical bifurcations.Unlike Ref. [16] where the amplitude increases from zero in a continuous fashion, in our measurement the amplitudes jump sharply from zero in a subcritical bifurcation.
Next, we fix ω p at 2ω 2 and increase the noise intensity while maintaining the same effective temperatures in the two resonators, D 1 = D 2 = D.The noise induces switching of each of the two resonators at random times.We measure the time intervals during which the 4 states are occupied, and obtain the stationary probability distributions w st (σ 1 , σ 2 ).For brevity we indicate the states σ = 1 and σ = −1 by ↑ and ↓, respectively, as we also did in Fig. 2. Therefore the 4 states are ↑↑, ↑↓, ↓↑ and ↓↓.The areas of the circles in Fig. 3b are proportional to the measured stationary probability distribution.We find that w st (↑↓) and w st (↓↑) exceed w st (↑↑) and w st (↓↓), consistent with notion that the electrostatic coupling favors opposite vibration phases in the two resonators, so that K 12 , K 21 < 0. The measured w st are in good agreement with Eq. (11).
The change of the state populations comes from the change of the switching rates.From Eq. ( 5) applied to two resonators, the rate of switching from the state σ i of resonator i is changed by exp(−K ij ) if σ j = σ i , i.e., the phases of the two resonators are almost equal, and by exp(K ij ) if the phases are opposite.For two coupled resonators, there are a total of 8 transitions, as illustrated in Fig. 3b.In the experiment, each of the 8 switching rates is individually measured, by fitting to the Poisson distribution of the residence times.Measurements are performed both before and after the coupling is turned on to give Wi and W i (σ i , σ j ) respectively [for two resonators, we use the notation W i (σ i , σ j ) rather than W i (σ i , {σ j })].The ratio W i (σ i , σ j )/ Wi represents the modification of the switching rate of resonator i due to coupling.
Figure 3c plots the logarithm of the ratio W i (σ i , σ j )/ Wi for the two resonators as a function of 1/D, where red and blue results correspond to switching of resonator 1 and 2 respectively.For the upper branches where the phases are identical, the switching rates are increased due to coupling, and vice versa for the lower branches.The lines are linear fits through the origin from which the change of the activation barriers ∆R i can be obtained by taking the negative values of the slopes.
We observe that, in agreement with Eq. ( 5), within the measurement uncertainty.Therefore in Fig. 3c we show the logarithm of the ratio of the average values of W i (σ i , σ j ) and W i (−σ i , −σ j ) to Wi .
Figure 3d shows that |∆R 1 | and |∆R 2 | are proportional to the square of potential difference V cpl between the two vibrating plates, with different proportionality constants for the two resonators.The measured values of |∆R i | are compared in Fig. 3d with Eq. ( 6) evaluated with the numerically simulated values of the logarithmic susceptibility and the vibration amplitudes and phases A j , ϕ j independently calculated for each resonator in the absence of coupling.There is good agreement between the entirely independent measurements with the coupling (circles) and the simulations with no coupling [the lines based on Eq. ( 6)].In turn, the simulations with no coupling are in excellent agreement with the measurement of the logarithmic susceptibility with no coupling, as seen from Fig. 2. The linear dependence of log[W i (σ i , σ j )/ Wi ] on 1/D in Fig. 3c confirms the proposed mechanism of the strong effect of even weak coupling, provided the noise is also weak.
As discussed earlier, a difference between |∆R 1 | and |∆R 2 |, and hence K 12 and K 21 , implies that detailed balance is broken, giving rise to a net probability current.For two resonators, the stationary probability distribution can be calculated (Appendix C), and then Eq. (7) gives In the experiment, the probability currents are obtained by taking the product of the measured stationary probability distribution and the measured switching rate out of the specific state.They are represented by block arrows in Fig. 3b.The lengths of the arrows are chosen to be proportional to the product of the measured stationary probability distribution and the measured switching rate.
Our measurement demonstrates the lack of detailed balance, as evident from the difference in length of each pair of arrows.The magnitude of the net probability current for the four branches are identical within the measurement uncertainty (Supplementary Note 5).As denoted by the purple arrow in Fig. 3b, the net probability current flows in the clockwise direction for ω 2 − ω 1 = -0.4Hz.
In our system of two coupled resonators with near identical damping, the sign and magnitude of the probability current are largely determined by the frequency mismatch ∆ω = ω 2 − ω 1 if the coupling and the noise intensity are fixed.When ∆ω is changed to 0.4 Hz by adjusting V R,2 , we find that the sign of the net probability current is reversed.Figure 4a plots the net probability current averaged over the four branches as a function of ∆ω.The line represents the probability current predicted by Eq. ( 9) with K 12 and K 21 given by the simulated value of the logarithmic susceptibilty of a single resonator using Eq. ( 6).
The difference between K 12 and K 21 , and hence the probability current, can be tuned to zero by choosing ∆ω.In our system, choosing ∆ω equal to zero makes the probability current vanish within measurement uncertainty.Detail balance is restored.The two resonators therefore map to the symmetric Ising model.The stationary distribution w st found in the experiment in this case coincides with the standard expression w st ({σ i }) ∝ exp( K ij σ i σ j /2) (Supplementary Note 8).We further show in the SM that while w st (↑↓) = w st (↓↑) exceed w st (↑↑) = w st (↓↓) due to the coupling, the switching rates given by Eq. ( 7) lead to vanishing of the net probability current.Probability current for two non-identical coupled parametric resonators Dependence of probability current on the frequency mismatch ∆ω between the two resonators at V cpl = 0.3 V. Purple circles are measurement.Calculations for two resonators based on the simulated logarithmic susceptibility of individual units are plotted in black.The straight line is a linear fit through the origin.Inset: For the considered weak coupling the frequency anticrossing as a function of ω2 − ω1 is undetectable.The color represents the sum of the amplitudes of forced vibrations of the two modes in nm; x-axis is the bias VR,1 (V) that controls ω1, whereas y-axis is the frequency of resonant drive (Hz) applied to both resonators.Red squares and blue circles mark the values of ω1 and ω2 used in the main figure.

III. DISCUSSION
Our results demonstrate that a system of slightly different parametric oscillators provides a long-sought inorganic implementation of an asymmetric Ising model.The parameters of the model are determined by the oscillator parameters, including the eigenfrequencies and the coupling, as well as the amplitude and frequency of the parametric modulation.These parameters can be controlled in a broad range.For oscillators based on microand nanomechanical resonators, this opens a way of creating asymmetric Ising networks with variable coupling strength and variable connectivity, which is the problem of interest for diverse disciplines, from biology to artificial intelligence.Besides these applications, such networks provide a conceptually simple setting for studying features of many-body dynamics away from thermal equilibrium.One of the major generic features is the lack of detailed balance, which leads to the onset of a probability current in the stationary state.
We measured the stationary probability current in an asymmetric Ising system.The magnitude of the current depends exponentially strongly on the interrelation between the coupling of the oscillators and the intensity of the noise in the range where both are small.Our analysis and measurement are done in the regime where the coupling-induced change of the oscillator frequencies is much smaller than the frequencies themselves, and the noise-induced spread of the vibration amplitudes is much smaller than the amplitudes themselves.Yet the ratio of the properly scaled coupling and noise intensity can be arbitrary.We note that, for a parametrically excited oscillator, noise necessarily comes along with relaxation, so that it is present even in the quantum regime.
The experiment shows that the effect of weak coupling of parametric oscillators can be quantitatively described in terms of an entirely independent effect of an additional drive at half the modulation frequency applied to an individual oscillator.It is demonstrated that, in a broad range of the drive amplitudes, the drive leads to a change of the logarithm of the rate of switching between the vibrational states of the oscillator, which is linear in the drive amplitude.The corresponding logarithmic susceptibility was measured and found to be in excellent agreement with simulations.
The stationary state of an asymmetric Ising model is not known, generally.This is not a consequence of disorder.A simple "ordered" system that maps onto an asymmetric Ising system is a chain of parametric oscillators where the coupling to the nearest neighbors for the oscillators on even and odd sites is different.The coupling parameters take on two values, K e = K 2n 2n±1 and K o = K 2n+1 2n+1±1 .For small |K e,o | one can analyze the spin dynamics similar to how it was done by Glauber [38] for a symmetric chain (Supplementary Note 9).In particular, we find that there are two spin-diffusion waves for a periodic chain; in an asymmetric model the wave frequencies, rather than being purely imaginary, can be complex, generally.The probability current in the stationary state is ∝ K o − K e .This model immediately extends to a square lattice, which can address the question of the possibility of an Onsager-type transition for an asymmetric Ising model with nearest-neighbor coupling.We note that, for an asymmetric Ising model, the eigenvalues of the balance equation can be complex in the general case, in contrast to a symmetric Ising model.For each resonator i (i = 1, 2) shown in Fig. 1b, the top plate is subjected to electrostatic torques exerted by the left and right electrodes underneath.If the potential difference between the two top plates

IV. ACKNOWLEDGEMENT
there is also an electrostatic attraction between the two resonators.Each top plate is connected to the input of an amplifier that is a virtual ground for ac voltages.On the left electrodes, the ac component V p,i cos(ω p t) controls the modulation of the spring constant via electrostatic springs softening.When a symmetry breaking torque is needed to measure the logarithmic susceptibility, a second ac component V d,i cos[(ω p /2)t+ϕ i ] is added.The noise voltage V n,i (t) generates the noise torque to induce transitions between the two states.
Voltages on the right electrodes only contain dc components.They are adjusted for fine tuning of the resonant frequencies of the two resonators in order to maintain the desirable difference of the oscillator frequencies ∆ω.Coupling between the two resonators is controlled by the potential difference between the top plates via V cpl (Supplementary Note 1).
Vibrations in each resonator are detected by measuring the change of capacitance between the top plate and the two underlying electrodes.The dc voltages described above lead to build up of charges on the top plates.As each of the top plates rotates, the capacitances with the two underlying electrodes change.Charge flowing out of the two top plates are detected independently by two separate amplifiers.The outputs of each amplifier is fed into a lock-in amplifier referenced at ω p /2.

APPENDIX B: MEASUREMENT OF SWITCHING RATES.
To measure the switching rate of an individual resonator, its oscillation phase φ is recorded as a function of time using a lockin amplifier.Figure 2a shows part of a record for resonator 1.If the resonator initially resides in the state σ = −1 with φ ≈ π, we identify that it has switched to the σ = +1 state with φ ≈ 0 when the phase goes over the threshold ε, where π/4 < ε < π/2.In switching from the initial state σ = +1 with φ ≈ 0 the phase with overwhelming probability jumps to −π ≡ π(mod2π).In this case the threshold is −π + ε.As the resonator switches back and forth between the two states, we record two sequences of residence times for the two states separately.The residence times in each state are plotted as a histogram.A typical histogram is shown in Fig. 5.The exponential decrease in the histogram confirms that the transitions are random and uncorrelated in time.An exponential fit to the histograms yields the switching rate.Fitting to a separate histogram gives the rate of switching from another state.We check that the measured switching rate does not depend on the choice of ε.For example, in Fig. 2a, the dark and light lines indicates two difference choices of threshold ε.They yield measured switching rates that are equal within the error bar of the fitting.
For uncoupled oscillators in the absence of the symmetry breaking drive, the measured switching rates out of the two states of each resonator are identical to within experimental uncertainty.Their value gives Wi for resonator i.Moreover, the stationary probability distributions w st (↑↓), w st (↓↑), w st (↑↑) and w st (↓↓) are measured to be equal to within measurement uncertainty (Supplementary Note 3).
To measure the logarithmic susceptibility of a single resonator, the switching rates are measured after the symmetry breaking drive is turned on.The fractional change of the switching rates for the two states are opposite in sign, as illustrated for resonator 1 in Fig. 2c.Logarithmic susceptibility can be calculated using the method of optimal fluctuation [39] or found from simulations [43].The results have been established to be in excellent agreement.Therefore here we directly used simulations to find the magnitude χ and the phase δ of the logarithmic susceptibility.To do this we incorporated the drive F d cos(ω p t/2 + ϕ d ) into the equation of motion (1) of resonator 1 and set the coupling parameters V ij equal to zero.We then switched to the rotating frame and used the standard rotating wave approximation to reduce the problem to a set of equations for the quadratures of q 1 (t).Forced vibrations at frequency ω p /2 in the lab frame correspond to stable stationary solutions of the equations for the quadratures in the absence of noise.Noise causes switching between these states.The residence times are identified and used to calculate the switching rate in a manner similar to the measurement procedure described above.This allowed us to avoid simulating multiple (≳ 10 7 − 10 9 in our case) oscillations of the parametric oscillator in the lab frame.
To measure the Ising model parameters K 12 and K 21 , the switching rates are measured before and after the coupling is turned on.The fractional change of the switching rates for resonators 1 and 2 are plotted in red and blue respectively in Fig. 3c.The dynamics of the chain of coupled parametric oscillators is mapped on the dynamics of Ising spins by associating the stable vibrational states of the oscillators with spin-1/2 states.Fluctuations lead to random switching of the spins.The evolution of the distribution w(σ 1 , σ 2 , ...) ≡ w({σ i }) over the spin states is described by the balance equation, which can be written in the form with the switching rates given by Eq. ( 5).We note that, even if the rates Wi are different for different spins (different parametric oscillators), but the model is symmetric, K ij = K ji , Eq. ( 10) has the stationary solution which is just the thermal distribution of the conventional symmetric Ising model.
For a system of N spins (oscillators) Eq. ( 10) is a system of 2 N equations.For the case of 2 oscillators it can be solved (see Supplementary Note 7 for more details).The stationary probability distribution is This expression was used to obtain Eq. ( 8) for the probability current and also in Fig. 3.For a symmetric system, K 12 = K 21 = K, we have w st (1, 1)/w st (1, −1) = exp(−2K) independent of the values of W1,2 , whereas for an asymmetric system the populations depend on the interrelation between W1 and W2 .
APPENDIX D: BREAKING OF THE DETAILED BALANCE.
The lack of detailed balance, and thus the onset of the probability current in the asymmetric Ising model can be shown without knowing the stationary distribution.One has to compare the ratio of flipping an ith spin back and forth directly or with a kth spin flipped back and forth on the way.For a system with detailed balance the result should be the same.We now compare these ratios.To shorten the notations, we keep in the expressions for the rates only the spins σ i and σ k and explicitly indicate which of them is flipped; other spins are not flipped.
The detailed balance condition reads For asymmetric Ising model the equality does not hold: the right-hand side has an extra factor exp [4(K ik − K ki )σ i σ k ].We note that the result is independent of the switching rates Wi , Wk in the absence of coupling.
Supplemental Material: "Controlled asymmetric Ising model implemented with parametric micromechanical oscillators"

SUPPLEMENTARY NOTE 1: EXCITATION AND DETECTION SCHEMES
The top plate of each resonator i (i = 1, 2) is subjected to electrostatic torques exerted by the left and right electrodes.If the potential difference between the top plates is non-zero, there is also an electrostatic attraction between the two top plates.Each top plate is connected to the input of an amplifier that is a virtual ground for ac voltages.
For resonator i, the total electrostatic torque from the two electrodes is given by the sum , where the subscript α denotes the left (L) or the right (R) electrode, C L,i (C R,i ) is the capacitance between the top plate and the left (right) electrode.V L,i and V R,i are the voltages on the left and right electrodes, respectively.Expanding the torque about the initial angle θ 0,i gives: where θ i is the rotation angle of resonator i measured from θ i,0 .C ′ i , C ′′ i , C ′′′ i and C IV i denotes the derivatives of C i with respect to θ i , evaluated at θ 0,i .Higher order terms are neglected.
V R,i contains only a dc component that is used to adjust ω i via the electrostatic spring softening effect.V L,i consists of small ac voltages on top of a dc component V dc L,i .It is used to control the parametric modulation, the symmetry-breaking torque and the noise torque.Considering only the dc component and the term responsible for the parametric excitation where To induce transitions between the two coexisting vibration states, a noise voltage V n,i (t) is added to V L,i .When the logarithmic susceptibility is measured, an additional ac component V d cos( ωp 2 t) is added.Throughout the experiments we use the modulation frequency ω p = 2ω 2 when measuring noise induced switching in two coupled resonators.
As seen from Eqs. ( 13) and ( 14), the time-independent component of τ i modifies the system parameters.The term independent of θ i shifts the equilibrium position of θ i,0 .The term 1 2 α C ′′ α,i θ i ∆V 2 α,i , which is linear in θ i , leads to electrostatic spring softening due to the static potential differences, producing a shift in the resonant frequency. Here α C IV i θ 3 i ∆V 2 α,i modify the Duffing nonlinearity.The time-dependent components of τ i are responsible for the parametric modulation term and noise term in Eq. ( 13) of the main text, with For the symmetry breaking drive, Throughout the experiment, the dc voltages V top 2 and V dc L,2 for resonator 2 are fixed at 0 V and -1.00 V respectively, so that ∆V L,2 is maintained constant at 1.00 V.For resonator 1, V top,1 is changed to control the coupling with resonator 2 as explained in more details below.V dc L,1 is then adjusted to maintain ∆V L,1 constant at -1.28V.F p is set to be identical for the two resonators, by choosing The voltages V R,1 and V R,2 applied to the right electrodes are used to control ω 1 and ω 2 respectively, as mentioned earlier.Unlike V L,i , V R,i contains no ac components.Initially, ∆V R,1 and ∆V R,2 are chosen to be -0.5 V and 0.5 V respectively to bring ω 1 and ω 2 close to each other.Subsequently, small changes to ∆V R,1 and ∆V R,2 allow the fine tuning of ∆ω to the desired value via the electrostatic spring softening effect discussed above.Furthermore, V R,1 and V R,2 are adjusted regularly to compensate for long term drifts in the resonant frequencies.The adjustment is performed by first setting V cpl to 0 V. Subsequently, the vibration amplitude in response to the parametric modulation is measured at a specific ω p .The changes in amplitude from the value recorded at the beginning of the experiment are used to infer the shifts in the eigenfrequencies.Small changes to V R,1 and V R,2 are sufficient to bring the vibration amplitudes, and hence the eigenfrequencies, back to the original values for the experiment to continue.
The two plates have identical width of 200 µm and thickness of 2 µm.Coupling between them is generated when is applied to the interdigitated comb shaped electrodes with separation ranging from 3 µm to 5 µm at different locations.Even though the equilibrium positions of θ 1,2 , as counted from the horizontal axis in Fig. 1b in the main text, are non-zero, the static rotations are small.When there are no vibrations, the sidewalls on the two plates remain largely parallel and aligned with each other.The electrostatic potential energy between the two top plates is given by 1 2 C 12 (θ 1 , θ 2 )V 2 cpl , where C 12 is the capacitance between the two plates.If we disregard the small misalignment of the plates, C 12 contains a term λ(θ 1 − θ 2 ) 2 , where λ is a proportionality constant.This term determines the coupling energy in Eq. ( 13) of the main text.Specifically, V 12 = V 21 = λV 2 cpl .As shown in the inset in Fig. 4 of the main text, level anticrossing does not occur for the small V cpl used in this paper.The term ∝ λ(θ 2 1 + θ 2 2 ) in the coupling energy leads to changes of ω 1 and ω 2 , by δω 1 and δω 2 respectively, with δω 1 ∼ δω 2 .These changes are measured at the beginning of the experiment from the shifts in the peaks of the linear resonant response from their value at V cpl = 0 V.As described earlier, in the procedure to compensate for the long term drifts in ω 1,2 , V cpl is set to 0 V.After the compensation is completed, V cpl needs to be changed back to the target value.The previously recorded values of δω 1,2 are used to adjust the applied voltages on the electrodes.In particular, ω p is changed by 2δω 2 so that the parametric drive frequency coincides with ω 2 .V R,1 is also adjusted to account for the small difference between δω 1 and δω 2 so that ∆ω ≡ ω 2 − ω 1 is maintained at the desired value.

SUPPLEMENTARY NOTE 2: GENERATION OF NOISE
Voltage noise V n,i (t) is applied on the left electrode of resonator i to induce transitions between the two coexisting states of opposite phases.The noise voltage originates from the Johnson noise of a 50 Ω resistor at room temperature.After amplification, the noise voltage is bandpass filtered with center frequency f 0 = 3000 Hz and bandwidth f bd = 40Hz.For resonator i, the filtered noise voltage is then mixed with a carrier voltage at frequency f c,i to generate two sidebands centered at f c,i ± f 0 .For instance, in Fig. 2 of the main text, f c1 and f c2 are 12841.8Hz and 12851.8Hz respectively.The resonant frequencies ω 1,2 /2π lie within the corresponding upper sideband.The frequency difference f c,2 − f c,1 is chosen to be much larger than the frequencies Γ i /2π, 3γ i A 2 i /2πω p (i = 1, 2) that characterize the motion of the modes in the rotating frame (here A i is the vibration amplitude), so that the two resonators are effectively subjected to independent noise voltages, because the relevant spectral components originate from different frequencies of the pass band (note that the above characteristic frequencies are much smaller than 40 Hz).Finally, the noise voltages are multiplied by a factor c i proportional to √ Γ i /∆V L,i to give V n,i (t) so that the effective temperature is identical in the two resonators.In Eq. ( 1) of the main text, The noise correlation time is ∼ 2π/f bd .Therefore on the time scale of slow motion in the rotating frame the noise is effectively δ-correlated.

SUPPLEMENTARY NOTE 4: LOGARITHMIC SUSCEPTIBILITY OF RESONATOR 2
The theory curves for coupled resonators are generated using the simulated logarithmic susceptibilities of individual resonators.Figure 2 of the main text shows how the results of the measurements and simulations of the logarithmic susceptibility of resonator 1. Calculation of the logarithmic susceptibility of resonator 2 follows a similar procedure.First, the changes of the activation barrier are determined as a function of 1/D in a fashion similar to Fig. 2c of the main text.Figure S2a shows the ratio of switching rates with and without the symmetry breaking drive for the two states of resonator 2 as a function of 1/D at ω p = 2ω 2 .The slope of the linear fit yields the change of the activation barrier due to the symmetry breaking drive.Figure S2b shows the dependence of the change in activation barrier on F d .The slope of the linear fit gives the logarithmic susceptibility of resonator 2.

SUPPLEMENTARY NOTE 5: MEASUREMENT OF PROBABILITY CURRENT
For two coupled resonators, there are 8 transitions as shown in Fig. 3b of the main text.In the experiment, we measure the switching rate for each transition as well as the stationary populations of the four states.The probability  current from one state A to another state B is obtained using the following procedure.First, we evaluate the product of the measured stationary population of state A and the measured transition rate out of state A to state B.Then, we calculate the product of the measured stationary population of state B and the measured transition rate out of state B to state A. In general, these two products are not equal to each other.The net probability current from state A to state B is the difference between these two products.We find that the net probability currents for the four pairs of states are identical to within measurement uncertainty, as shown in Fig. S3 for three different values of ∆ω.The values of the probability current plotted in Fig. 4 in the main text represents the mean value of the measured probability current for the four pairs of states in Fig. S3 for ∆ω = −0.4Hz.

SUPPLEMENTARY NOTE 6: CHANGES OF ACTIVATION BARRIER BEYOND THE LINEAR REGIME
The effect of the coupling on the switching rate was seen in Ref. [14], but the logarithmic-susceptibility regime was not identified there.We emphasize that the mapping on the asymmetric Ising model in our study applies in the regime where the change of the activation barrier has the form ∆R i (σ i , {σ j̸ =i }) = D i σ i j̸ =i K ij σ j .In particular, for two oscillators ∆R i (σ i , σ j ) = −∆R i (−σ i , σ j ) = −∆R i (σ i , −σ j ).The value of ∆R i (σ i , σ j ) in this regime is determined by the logarithmic susceptibility, and in our experiment we have measured it independently.
The logarithmic-susceptibility regime refers to weak coupling.If the coupling is stronger, even where it does not lead to the onset of new vibrational states, not only are the switching rates modified, but the very stable states of individual oscillators, i.e., their amplitudes and phases, significantly change depending on the states of other oscillators.Therefore one cannot think of an oscillator as having just two states, which underlies the mapping on a spin system.When the amplitude F d of the symmetry-breaking drive is small for an individual, uncoupled resonator, the change in the vibration amplitudes of the two states is also small.The changes of the activation barriers ∆R i (σ i = +1) and ∆R i (σ i = −1) are of opposite signs but near identical magnitude.|∆R i | is proportional to the amplitude of the symmetry breaking drive, with a proportionality constant given by the logarithmic susceptibility.As F d increases, the difference in the vibration amplitudes and phases of the two states cannot be ignored.Figure S4a extends the range of F d of Fig. 2d in the main text to show that, for resonator 1, |∆R 1 (σ 1 = +1)| and |∆R 1 (σ 1 = −1)| are no longer equal for large F d .Furthermore, the dependence of ∆R 1 (σ 1 ) on F d does not follow a linear relationship.Figure S4b shows a similar plot for resonator 2. As described in the main text, when weak coupling is turned on between two resonators, the effects on resonator 1 from the coupling to resonator 2 can be understood in terms of the logarithmic susceptibility of resonator 1 subjected to a symmetry breaking drive, and vice versa for resonator 2. Beyond the weak coupling regime, the increase in switching rates for initial states with identical phases is no longer equal in magnitude to the decrease of switching rates for initial states with opposite phases.In fact, the phase difference between the two vibration states in each oscillator deviates considerably from π. Figure S5 shows the measured and simulated results for the change of the activation barrier when the range of V cpl is extended beyond that in Fig. 3d in the main text.While Fig. 3d in the main text plots |∆R 1,2 |, Figure S5a and Figure S5b show ∆R 1 (σ 1 , σ 2 ) and ∆R 2 (σ 1 , σ 2 ) separately and without the absolute value.Since W i (σ i , σ j ) = W i (−σ i , −σ j ) as described by Eq. ( 8) in the main text, it follows that ∆R i (σ i , σ j ) = ∆R i (−σ i , −σ j ).In Fig. S5, the up triangles represent the average of the measured values of ∆R i (↑↓) and ∆R i (↓↑).The down triangles represent the average of ∆R i (↑↑) and ∆R i (↓↓), both of which are negative.The results for Fig. 3d in the main text are obtained by averaging the difference in the magnitude of ∆R 1,2 for initial states of the same and opposite phases becomes more apparent.In Fig. S5a and Fig. S5b, the sum of the top and bottom branches are shown as the thick lines, the deviation of which from zero increases with V 2 cpl .

SUPPLEMENTARY NOTE 7: POPULATION AND TRANSITION RATES OF COUPLED IDENTICAL RESONATORS
In this section, we tune the difference between K 12 and K 21 to near zero by choosing ∆ω/2π = 0 Hz.As a result, the probability current vanishes to within experimental uncertainty.Detailed balance is restored and our system maps onto the symmetric Ising model.
Figure S6a shows the vibration amplitudes of the two resonators with ∆ω/2π = 0 Hz under parametric modulation in the absence of injected noise.The slight difference between the response curves is due to the difference in the decay rates of the modes.Application of noise leads to switching in both resonators, as illustrated in Fig. S4b.As a result of the coupling (V cpl = 0.3 V) that favors opposite phases in the two resonators, the transition rates W i (σ i , {σ j = σ i }) and W i (σ i , {σ j = −σ i }) are increased and decreased compared to the uncoupled values Wi respectively.Moreover, the fractional changes for resonators 1 and 2 are identical to within experimental uncertainty, as shown by the adjacent red and blue circles in Fig. S6d.Figures S6f and S6g  We examine how detailed balance is manifested in the probability current, taking the transitions between the ↑↑ state and ↑↓ state as an example.The transition rates W 1 (↑, ↑) and W 1 (↑, ↓) are changed from W1 by factors exp (−K) and exp(K) respectively.The stationary probability distribution for K ij = K ji is w st = Z −1 exp 1 2 i,j(i̸ =j) Kσ i σ j where Z is the normalization factor.In the case of two resonators we see that w st (↑↑) and w st (↑↓) are changed from their value 1  4 in the absence of coupling by Z −1 exp(K) and Z −1 exp(−K) respectively, where Z = 4 cosh K (Fig. S6c).
The products W 1 (↑, ↑)w st (↑↑) and W 1 (↑, ↓)w st (↑↓) remain equal in magnitude, as represented by the same length of the two block arrows in Fig. S6e.The net probability current between the states ↑↑ and ↑↓ is therefore zero.

SUPPLEMENTARY NOTE 8: BALANCE EQUATION FOR THE STATE POPULATIONS
As discussed in the main text, on the time scale long compared to the relaxation times of the modes, the dynamics of the chain of coupled parametric oscillators is mapped on the dynamics of Ising spins, which switch at random between their states.Two states of an individual spin σ i = ±1 correspond to two values of the phase of parametrically excited vibrations.The switching rate W i (σ i , {σ j̸ =i }) of the transition σ i → −σ i of spin i depends on the state of spins j ̸ = i with which the spin is coupled.The spins are switching at random, and the evolution of the distribution w(σ 1 , σ 2 , ...) ≡ w({σ i }) over the spin states is described by the balance equation, which can be written in the form The switching rates have the form and for the spins that model differing oscillators K ij ̸ = K ji , generally, cf.Eqs. ( 5) and ( 6) of the main text.The evolution of the distribution w({σ i }) is determined by the eigenvalues of the 2 N × 2 N matrix W i (σ i , {σ j̸ =i }), where N is the number of spins.For a symmetric Ising model all eigenvalues of this matrix are real.To see this, we change from the distribution w to w, The off-diagonal elements of the matrix in the right-hand side, Wi , are independent of the spin state and are thus the same for σ i = 1 and σ i = −1.Therefore all eigenvalues of Eq. ( 17), and thus also of Eq. ( 15), are real.
For an asymmetric Ising model the eigenvalues of the matrix of the transition rates can be complex.This is another significant difference between symmetric and asymmetric Ising models.
The stationary distribution w st ({σ i }) is given by the solution of Eq. ( 15) for ẇ = 0. Formally, the condition ẇ = 0 makes Eq. (15) a system of 2 N linear equations.Only 2 N − 1 of them are linearly independent, as seen from the condition {σi} ẇ({σ i }) ≡ 0 that follows from Eq. (15).Therefore they have to be complemented by the normalization condition {σi} w st ({σ i }) = 1.Only 2 N −1 populations w st ({σ i }) are linearly independent.Indeed, as seen from Eq. ( 16), the rates W i (σ i , {σ j̸ =i }) do not change if we simultaneously change the signs of all spins.Therefore, for any spin configuration {σ i } in a finite system, w st ({σ i }) = w st ({−σ i }) (at this time we do not consider spontaneous symmetry breaking).
The distribution w st can be easily found if the spin system is symmetric, K ij = K ji .From Eq. ( 15), in this case where Z is the normalization factor, which is essentially given by the standard partition function of the Ising system with the appropriately defined coupling constants scaled by the temperature.The relation (18) holds even if the chain is disordered, i.e., the values of the prefactor in the switching rate Wi depend on site i.It immediately follows from Eqs. ( 16) and ( 18) that the probability current defined by Eq. ( 7) in the main text is equal to zero for a symmetric system.

Figure 1 .
Figure 1.Two coupled parametric torsional resonators.a. Scanning electron micrograph of two torsional resonators located side-by-side.The scale bar measures 100 µm.b.Schematic of the actuation scheme and measurement circuitry.Voltages applied to the left electrodes generate the parametric modulation at ωp, the drive at ωp/2 and the noise.Voltages V dc R,i applied to the fixed electrodes allow fine tuning of the resonant frequencies.The dc voltage differences between each top plate and the underlying electrodes leads to an ac current flowing out of the top plate as it rotates.Capacitive coupling between the two plates is controlled by the voltage difference V top 1 − V top 2

cFigure 2 .
Figure 2. Measurement of the logarithmic susceptibility of a single resonator Coupling between the two resonators is turned off.We present results for resonator 1 and indicate the states σ1 = 1 and σ1 = −1 by ↑ and ↓, respectively.a.In the presence of noise, resonator 1 randomly switches between two coexisting vibration states with opposite phase.The two light grey lines are thresholds for identifying phase switches.The dark grey lines represent another choice of threshold.A drive at half the modulation frequency with amplitude F d = 1.04 × 10 −17 Nm breaks the symmetry and renders the residence times, and thus the stationary populations of the states ↑ and ↓ different.b.The ratio wst(↑)/wst(↓) increases as F d increases.Circles are measured results for the chosen drive phase ϕ d = 3.3 o .The solid line represents theory calculated using the simulated logarithmic susceptibility.Inset: same data shown in semilog scale.c.Logarithm of the ratio of the switching rates from states ↑ and ↓ with the resonant drive turned on, W1(↑) and W1(↓) (up and down triangles, respectively), to the rate with no drive W1 = C1 exp(− R1/D), plotted as a function of 1/D.The switching rates are modified by different amounts for the two states according to Eq. (4) The increments of the effective activation energies ∆R1(↑) and ∆R1(↓) are obtained from the slopes of the linear fits through the origin.d.Increment |∆R1| as a function of F d for resonator 1.The slope of the linear fit through the origin yields χ1 cos(ϕ d + δ1) defined in Eq. (4).Measurements are shown in red.Numerical simulations are shown in pink.

Figure 3 .
Figure 3. Asymmetric Ising model implemented with two coupled parametric oscillators a. Vibration amplitudes of resonators 1 (red) and 2 (blue), with ∆ω/2π = -0.4Hz and V cpl = 0.3 V, under identical parametric modulation with no noise added.The arrow marks ωp/2 for measuring noise-induced switching for the rest of the figure.b.Switchings between the four states of the two-resonator system in the presence of noise.The areas of the circles are proportional to the measured stationary populations wst(σ1, σ2) (the first arrow from the left refers to σ1 and the second arrow refers to σ2).The lengths of the straight arrows between the circles are proportional to the products of the measured switching rates Wi(σ1, σ2) and the corresponding populations wst(σ1, σ2).The purple arrow represents the net probability current.c.Logarithm of the measured changes of the switching rates due to coupling as a function of 1/D.The values of ∆Ri(σi, σj) ≡ −D log[Wi(σi, σj)/ Wi] are determined by the slopes of the linear fits.The difference between |∆R1(σ1, σ2)| and |∆R2(σ2, σ1)| for the same pairs (σ1, σ2) is identified from the different magnitudes of the slopes.This difference determines the asymmetry of the Ising model.d.Dependence of |∆R1| (red) and |∆R2| (blue) on V 2 cpl that is proportional to the coupling constant.The values of |∆Ri| are the average values of |∆Ri(σi, σj)| for σi = σj and σi = −σj.The pink and light blue lines are obtained from theory based on the independently simulated logarithmic susceptibilities of individual uncoupled resonators.
Figure 4.Probability current for two non-identical coupled parametric resonators Dependence of probability current on the frequency mismatch ∆ω between the two resonators at V cpl = 0.3 V. Purple circles are measurement.Calculations for two resonators based on the simulated logarithmic susceptibility of individual units are plotted in black.The straight line is a linear fit through the origin.Inset: For the considered weak coupling the frequency anticrossing as a function of ω2 − ω1 is undetectable.The color represents the sum of the amplitudes of forced vibrations of the two modes in nm; x-axis is the bias VR,1 (V) that controls ω1, whereas y-axis is the frequency of resonant drive (Hz) applied to both resonators.Red squares and blue circles mark the values of ω1 and ω2 used in the main figure.
This work is supported by the Research Grants Council of Hong Kong SAR (Project No. 16304620) and partially supported by Project No. HKUST C6008-20E.MID was supported in part by the National Science Foundation through Grant No. DMR-2003815.

Figure 5 .
Figure 5. Histogram of the residence times.The residence times are recorded for resonator 1 switching out of the σ1 = +1 state at F d = 0, D = 3.01 × 10 −6 N 2 kg −2 Hz −1 and ωp/2 = ω1.The slope of the linear fit gives the rate of switching out of this state.
APPENDIX C: THE BALANCE EQUATION.
Fig. S1.a. Stationary populations of the 4 states of two uncoupled resonators.b.Measured rates of switching out of the states.For the red(blue) circles, only resonator 1 (2) switches.
Fig. S2.a. Dependence of the logarithm of the ratio of switching rates with and without the symmetry breaking drive on 1/D at ωp = 2ω2, F d = 2.12 × 10 −17 Nm and φ d = 3.3 o for the two states in resonator 2. Up and down triangles represent initial states of σ2 = +1 and -1 respectively.The straight lines are linear fits through the origin.The negative values of their slopes give the change in activation barrier.b.The change in activation barrier |∆R2| is plotted vs. F d .The line is a linear fit through the origin.The slope of |∆R2| gives the logarithmic susceptibility for resonator 2. Measurements are shown in blue.Numerical simulations are shown in light blue.
Fig.S3.The net probability current for each of the four pairs of transitions in the opposite directions for two coupled resonators for ∆ω = -0.4Hz (solid squares), 0 Hz (hollow circles) and 0.4 Hz (solid circles).
Fig. S4.a. Increments |∆R1(σ1 = +1)| (up triangles) and |∆R1(σ1 = −1)| (down triangles) of the activation energy as a function of F d for resonator 1.The range of F d is about 5 times larger than in Fig. 2d in the main text.Measurements and numerical simulations are shown in red and pink respectively.The lines are fits to the simulated results by parabolae.They have the same linear term as the line in Fig. 2d.Different quadratic terms are used for σ1 = ±1.b.Similar plot for resonator 2. Measurements and numerical simulations are shown in dark blue and light blue respectively.

Fig
Fig. S5.a. Dependence of the change ∆R1 of activation barriers for switching of resonator 1 on V 2 cpl .Initial states with the phase of resonator 1 identical (opposite) to resonator 2 are represented by up triangles (down triangles).Measurements are shown in red.The thin pink lines are obtained from theory based on the effect of a symmetry breaking drive on individual uncoupled resonators using numerical simulations from Figures S4.The thick pink line is a sum of the thin lines.b.Similar plot for resonator 2. Measurements and numerical simulations are shown in blue and light blue respectively.
Fig.S6.a. Vibration amplitudes of resonator 1 (red) and 2 (blue), with ∆ω/2π = 0 Hz and V cpl = 0.3 V, under identical parametric modulation with no noise added.b.Typical records of the phase of the two resonators as a function of time when noise is added.c.Stationary probability distributions of the four states.d.Change of the transition rates from the four initial states due to coupling.For red (blue) circles, only resonator 1 (2) switches.e. Switchings between the 4 states.The areas of circles are proportional to the measured stationary probability distributions wst(σ1, σ2).The lengths of the arrows are proportional to the product of the measured switching rates and the corresponding initial stationary probability distribution.f.Logarithm of the measured changes of the transition rates due to coupling as a function of 1/D for resonator 1.The magnitude of the slopes of the linear fits yields |∆R1|.g.Similar plot for resonator 2, with |∆R2| given by the magnitude of the slope of the linear fit.