## Abstract

Solving intractable mathematical problems in simulators composed of atoms, ions, photons, or electrons has recently emerged as a subject of intense interest. We extend this concept to phonons that are localized in spectrally pure resonances in an electromechanical system that enables their interactions to be exquisitely fashioned via electrical means. We harness this platform to emulate the Ising Hamiltonian whose spin 1/2 particles are replicated by the phase bistable vibrations from the parametric resonances of multiple modes. The coupling between the mechanical spins is created by generating two-mode squeezed states, which impart correlations between modes that can imitate a random, ferromagnetic state or an antiferromagnetic state on demand. These results suggest that an electromechanical simulator could be built for the Ising Hamiltonian in a nontrivial configuration, namely, for a large number of spins with multiple degrees of coupling.

- electromechanical resonator
- Phonon
- Non-linear
- parametric resonance
- non-degenerate parametric amplification
- simulator
- Ising Hamiltonian

## INTRODUCTION

Physical simulators have emerged as a novel means for solving problems in quantum physics that are beyond the capacity of classical computers (*1*, *2*). This notion can also be applied to mathematical problems, for example, the Ising Hamiltonian used in describing a spin glass, which is characterized as a nondeterministic polynomial-time (NP)–hard problem (*3*). However, despite this formidable challenge, the ground-state spin configuration of an Ising system could yield solutions to optimization problems that can be mapped onto its spin lattice, and thus, its efficient extraction is highly desired (*4*). One approach to rapidly determine the ground state of the Ising Hamiltonian is to use quantum annealing where quantum tunneling is harnessed to search the energy landscape corresponding to a given problem programmed into the underlying spin lattice (*5*–*7*). Recently, an alternative and apparently classical approach to this problem has emerged where a time-multiplexed optical parametric resonator network is programmed with the Ising Hamiltonian. The ground state is then determined by slowly activating the network, which preferentially resonates in its global potential minima (*8*, *9*). Here, a variation of this concept is developed with phonons in a frequency-multiplexed electromechanical parametric resonator (*10*, *11*), and we show that this platform could readily be extended to multiple parametric resonances with arbitrary degrees of coupling, namely, all the ingredients necessary to solve the Ising Hamiltonian in a nontrivial configuration (*12*–*14*).

The Ising model conceived using the tools of statistical physics describes ferromagnetism, and in the absence of a magnetic field, its Hamiltonian is given by with *N* particles each having two spin states σ_{i} = ± 1 where the coupling between the *i*th and *k*th particles is parameterized by *J*_{ik}. In the first steps to building a phonon-based simulator, the fundamentals of the Ising Hamiltonian need to be replicated in the electromechanical domain, namely, a mechanical spin that encodes σ, multiple mechanical spins that play the role of an *N* particle bath, and finally coupling between the mechanical spins *J* whose magnitude and polarity can be controlled and extended beyond nearest-neighbor particles.

## RESULTS

### Mapping the electromechanical system onto the Ising Hamiltonian

The prototype electromechanical system in which these concepts are investigated is shown in Fig. 1 (see Materials and Methods), and it consists of two strongly coupled mechanical beams that sustain symmetric (*S*) and asymmetric (*A*) vibration modes at ω_{S}/2π ≈ 298.4 kHz and ω_{A}/2π ≈ 310.3 kHz with bandwidths of Δω_{S}/2π ≈ 130 Hz and Δω_{A}/2π ≈ 138 Hz, respectively (*15*). The mechanical elements are integrated with piezoelectric transducers, which enable the underlying harmonic potential of both modes to be parametrically modulated, leading to a system Hamiltonian given by(1)where the summation expresses the kinetic and potential energies from both modes in terms of their position *q*_{n} and canonically conjugate momentum *p*_{n} (*15*–*17*). The potential energy term contains three contributions, with the second contribution arising from the periodic modulation of the mechanical spring constant with amplitude Γ_{n} at twice the natural mode frequency, which yields degenerate parametric amplification and parametric resonance (*18*, *19*), and the third contribution arising from the Duffing nonlinearity β_{n}, a well-known anharmonicity that emerges at large displacements (*20*). Γ_{n} can be experimentally activated by modulating the spring constant of either mode with the application of voltage *V*_{n}(2ω_{n}) to induce stress from the piezoelectric transducers, and at sufficiently large amplitude, it results in a parametric resonance as detailed in the Supplementary Materials (*16*, *19*). Projecting the parametric resonances into phase space under pulsed driving, via the demodulation circuit in Fig. 1 (which records the in-phase *Q*_{n} and quadrature *P*_{n} components of position), reveals that they can vibrate with only two phases separated by π radians, as shown in Fig. 2 (A and B) (*16*, *21*). These phase bistable vibrations provide the ideal means to encode a classical spin in the mechanical domain where the positive in-phase component is defined as spin up, that is, σ_{n} = 1, and vice versa as depicted by the arrows in Fig. 2 (A and B) (*8*, *12*).

The last term in Eq. 1 describes nondegenerate parametric down-conversion from the pump with amplitude Λ at the sum frequency of both modes , which results in the symmetric and asymmetric modes becoming correlated, yielding a two-mode squeezed state (*15*). Λ can be experimentally activated by piezoelectrically pumping the spring constant of both modes with voltage , which simultaneously amplifies their thermomechanical fluctuations, as shown in Fig. 2 (C and D) (*22*). The concurrent generation of phonons in this process leads to the vibrations of both modes becoming correlated, which can be identified by reconstructing their cross-quadratures in phase space, namely, the in-phase component of the asymmetric mode versus the quadrature component of the symmetric mode and vice versa, as shown in Fig. 2 (E and F). The resultant squashed distributions in phase space imply that the motion of both modes is perfectly intertwined and is statistically confirmed by their unity correlation coefficient as detailed in the Supplementary Materials (*15*).

These observations suggest that correlations generated between the two harmonic modes from nondegenerate parametric down-conversion, as schematically depicted in Fig. 3A with Γ_{n} = 0, and the double-well potential underpinning the phase bistable vibrations of a parametric resonance (*17*, *23*) provide all the ingredients necessary to realizing a phonon-based simulator for the Ising Hamiltonian. The key to implementing this vision is the ability to generate correlations between the double-well potentials underlying the parametric resonances from both modes through parametric down-conversion as visualized in Fig. 3A with Γ_{n} ≠ 0.

First, to theoretically confirm the viability of this approach, the above Hamiltonian is transformed in the rotating-frame approximation with the introduction of a new canonical position coordinate *Q*_{n} and a conjugate momentum *P*_{n}, as detailed in Materials and Methods, which yields(2)where the first term describes the quasi-energy separating the two stable oscillation phases of the parametric resonances at *P*_{n} = 0 and with their two phases encoding σ_{n} = ±1 (see the Supplementary Materials) and the second term quantifies the coupling between the modes. This electromechanical system can therefore be formally mapped onto the Ising Hamiltonian composed of two particles with spins corresponding to the phase bistable parametric resonances of both modes, which can be coupled via nondegenerate parametric down-conversion. This analysis also reveals that the polarity of the coupling can be tuned by the pump phase ϕ (*24*).

### Experimental implementation and analysis

To verify this concept, the protocol depicted in Fig. 3A is developed where a two-mode thermally squeezed state is initially created with as visualized by the correlated fluctuations of the balls (signifying mechanical motion) in the harmonic potentials of the symmetric and asymmetric modes. Next, both modes are simultaneously and slowly (that is, << Δω_{S}/2π and Δω_{A}/2π) activated via *V*_{n}(2ω_{n}), which results in their harmonic potentials evolving to the double-well potentials of their parametric resonance, as shown in Fig. 3A with Γ_{n} ≠ 0, where the balls have now stabilized in one of the two potential minima corresponding to either a spin up or a spin down. In the slow transition from the harmonic to the double-well potentials, an intermediate regime exists where thermal fluctuations (in addition to dissipation) can drive transitions between the two oscillation phases, namely, *E*_{n} < *k*_{B}*T*, where *k*_{B} is the Boltzmann constant and *T* is the temperature, which can stimulate the search for the global potential minima in an electromechanical Ising simulator. The spin information is then deleted by deactivating both parametric resonances, and the protocol is repeated for 2000 s. The spin information encoded in the two modes is identified by the polarity of the in-phase component of their parametric resonance, and it yields a train of switching outputs from both modes, with a period defined by this sequence, namely, 0.1 Hz, as shown in Fig. 3A. Implementation of this protocol over this duration provides a statistical ensemble from which the nature of the correlation between the two mechanical spins can be reliably and quantitatively evaluated.

Experimentally, in the case when (peak-to-peak voltage), the outputs from the symmetric and asymmetric modes reveal spin polarities, which are independent of each other, implying the absence of coupling between them or, in other words, *J*_{SA} = 0 as shown in Fig. 3B. The corresponding correlation coefficient extracted from this measurement is almost 0, as shown in Fig. 3C, which quantitatively confirms this observation. On the other hand, implementing this sequence with yields the output shown in the lower panel of Fig. 3B, which reveals that the mechanical spins always exhibit parallel alignment, namely, ferromagnetic coupling with *J*_{SA} >> *k*_{B}*T*. The corresponding correlation coefficient confirms this observation, yielding a value of almost 1, as shown in Fig. 3C. Next, implementing this protocol as a function of pump amplitude and extracting the resultant correlation coefficients yield the response shown in Fig. 3C. This indicates that the ferromagnetic state can be controllably created by the pump, with the corresponding correlation coefficient smoothly transitioning from 0 to 1 with a profile that is consistent with correlations between two hypothetical spins interacting via the Ising Hamiltonian (as detailed in Materials and Methods and shown by the solid line in Fig. 3C). Consequently, , namely, Λ, can convincingly play the role of *J*_{SA} in the electromechanical Ising simulator.

To control the sign of *J*_{SA} in the electromechanical Ising simulator, ϕ is adjusted as suggested by Eq. 2 and experimentally depicted in Fig. 1 (see the Supplementary Materials). To this end, the protocol outlined in Fig. 3A is reimplemented but now as a function of ϕ with , which ensures that the mechanical spins are perfectly coupled, as shown in Fig. 3C. The results of this measurement in terms of the extracted correlation coefficient are shown in Fig. 4A, which reveal that it can be continuously tuned from 1 to −1 crossing 0, where the latter two time series are shown explicitly in Fig. 4B. In other words, starting from a ferromagnetic state, the mechanical spins transition to an antiferromagnetic state via an uncoupled state; that is, the *J*_{SA} → −*J*_{SA} operation can be executed on demand via ϕ.

## DISCUSSION

Electro-optomechanical systems (*10*, *25*) have emerged as a versatile platform where ultra-precise sensors can be developed (*26*, *27*), dynamically engineered nonlinearities can be harnessed (*28*, *29*), and quantum mechanics within a macroscopic context can even be studied (*30*–*33*). The notion of solving mathematical problems with phonons localized in spectrally pure resonances, which is advanced here, offers a new chapter to the mechanical resonator narrative.

The results detailed in Figs. 3 and 4 confirm that the requisite features of the Ising Hamiltonian can be reproduced by phonons sustained by an electromechanical system. However, this implementation with two mechanical modes is trivial, and therefore, it is instructive to examine the prospects of a more useful Ising machine based on these ideas, as visualized in Fig. 5. Here, an array of mechanical elements with different frequencies are weakly mechanically coupled to their neighbors (*34*). The elements encode spin information via the bistable phase of their piezoelectrically activated parametric resonance in their fundamental mode, via the right clamping point, from where this information can be programmed and read out (*16*). The key difference here is that spin information stored in a given mode is concentrated within each mechanical element in contrast to the above demonstration. On the left clamping point, a coupling gate electrode is defined, which interconnects the mechanical spins via piezoelectrically activated parametric down-conversion at the sum frequency of two elements, for instance, . Naturally, the reduced mechanical coupling between the elements will require a more intense pump to compensate, which is readily available to this architecture (*15*). Uniquely, in this scheme, each mechanical spin can then be easily coupled to all its neighbors by using a frequency division–multiplexed (FDM) pump composed of multiple sum frequencies, where the availability of arbitrary coupling between the mechanical elements via the FDM pump is depicted in Fig. 5. The compact and highly flexible form of this coupling is in stark contrast to the optical Ising machine, which requires delay lines that increase both in number for more spins and in length for higher-order couplings (*9*). The FDM piezoelectric pump, in principle, can enable a large number of mechanical spins to sustain multiple degrees of coupling, thus permitting the electromechanical Ising simulator to be programmed to explore problems that challenge conventional computers. However, in practice, the ultimate number of spins, with maximal coupling (that is, *N*^{2}), will be limited by the global coupling gate’s ability to sustain the sum of the *N*^{2} coupling voltages before its Schottky barrier breaks down and neutralizes the piezoelectric transduction.

Although this platform could probe NP-hard Ising problems (*8*, *9*), it does not offer speed-up because it uses classical annealing where the thermomechanical fluctuations of the mechanical elements drive the search in the underlying potential energy landscape for the ground state. However, if all the mechanical elements can be operated in their ground state, then quantum effects could be potentially harnessed to explore the possibility of increasing the speed of search (*30*–*32*). Alternatively, the universal nature of this concept allows for it to be exploited by any kind of resonator, even from superconductors that sustain both parametric resonances (*35*) and nondegenerate parametric down-conversion (*36*), thus enabling the quantum dynamics of this concept to be explored.

## MATERIALS AND METHODS

### Experimental

The electromechanical system was synthesized from a GaAs/AlGaAs heterostructure sustaining a 100-nm-thick, uniformly doped GaAs layer located 300 nm below the surface (shown in blue in Fig. 1) via conventional micromachining processes (*15*). The mechanical elements have a length, width, and thickness of 80 μm, 20 μm, and 800 nm, respectively, which are strongly intercoupled via two 40-μm-wide overhangs with the same thickness and a length of 16 μm. Piezoelectric transducers were incorporated into the clamping points of both mechanical elements, which are composed of the doped layer and gold electrodes sandwiching an insulating GaAs/AlGaAs superlattice.

The parametric resonances and nondegenerate parametric down-conversion were piezoelectrically activated with multiple AC signal generators (NF Wave Factory 1974) connected to the gold electrodes on the left element, as depicted in Fig. 1, which modulated the stress and hence the spring constant of the symmetric and asymmetric modes. The thermal motion and the parametric resonances of both modes were probed with a HeNe laser reflected from the right mechanical element with a 3-μW input and detected in a Doppler interferometer (Neoark MLD-221). The output from the interferometer was demodulated either in a spectrum analyzer (Agilent 89410A) or in two phase-sensitive detectors (Stanford Research Systems SR844) that were mixed with two local oscillators (NF Wave Factory 1974) locked onto the resonances of both modes.

### Theory

The Hamiltonian for a resonantly excited degenerate parametric oscillator with natural frequency ω_{0} can be expressed as (*37*)(3)This Hamiltonian can be translated into the rotating frame following the standard procedure (*38*) with the introduction of canonically conjugate variables *P* and *Q* defined as(4)The effective Hamiltonian in these variables is then given by(5)where all the off-resonant and rapidly oscillating terms have been omitted (*37*). This Hamiltonian has two minima, extracted via , at(6)separated by a saddle point at *P* = *Q* = 0 (see the Supplementary Materials). At steady state, the bistable oscillating states of the parametric resonator are energetically degenerate with , and they can mimic a classical two-level or a spin 1/2 system via *Q*, as experimentally confirmed in Fig. 2 (A and B).

The Hamiltonian for two degenerate parametric oscillators that are coupled via nondegenerate parametric down-conversion is given in Eq. 1, and it can also be translated into the rotating frame via Eq. 4, which yields(7)If the nondegenerate parametric coupling strength Λ is weak (that is, in the amplification regime as detailed in fig. S1C), then the steady state in Eq. 6 can be used to approximate the lowest-order solution yielding Eq. 2.

### Correlations between spins

The theoretical correlations between two spins, whose interaction is mediated by the Ising Hamiltonian, were evaluated to verify the experimental variation of the correlation coefficient as a function of pump amplitude and phase to confirm the viability of the electromechanical simulator. For a pair of spins, the probability of finding a particular spin configuration is given by(8)where the partition function is defined as(9)with α = 1/*k*_{B}*T*. The correlation between nearest-neighbor spins can then be expressed as(10)where this metric is analogous to the correlation coefficient. This equation can reproduce the experimental response using a fitting parameter defined as , which yields the line in Fig. 3C. This analysis reveals that as the pump intensity is increased, a ferromagnetic state can be created, the activation of which is consistent with the Ising Hamiltonian where the correlation generated via parametric down-conversion competes with the random thermal fluctuations and is balanced at , where *J*_{SA} = *k*_{B}*T*. A similar analysis can also be performed as a function of pump phase with *J*_{SA}cos(ϕ), yielding 〈σ_{S}σ_{A}〉 = tanh(α*J*_{SA}cos(ϕ)), which reproduces the experimental result shown by the line in Fig. 4A.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/6/e1600236/DC1

I. Degenerate and nondegenerate parametric amplification

II. The double-well potential

III. Pump phase

fig. S1. Experimentally measured degenerate and nondegenerate parametric amplification of both modes in the electromechanical system.

fig. S2. The double-well potential underpinning a parametric resonance.

fig. S3. The pump phase dependence of the two-mode squeezing.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank Y. Ishikawa and K. Onomitsu for growing the heterostructure and H. Takesue and T. Inagaki for discussions.

**Funding:**This work was partly supported by MEXT (Ministry of Education, Culture, Sports, Science and Technology) KAKENHI grant no. 15H05869.

**Author contributions:**I.M. and H.Y. conceived the idea. I.M. performed the measurements, analyzed the data, and wrote the paper. H.O. fabricated the device. H.Y. developed the analytical model and planned the project.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors