Microwave Photonic Ising Machine

Ising machines based on analog systems have the potential of acceleration in solving ubiquitous combinatorial optimization problems. Although some artificial spins to support large-scale Ising machine is reported, e.g. superconducting qubits in quantum annealers and short optical pulses in coherent Ising machines, the spin coherence is fragile due to the ultra-low equivalent temperature or optical phase sensitivity. In this paper, we propose to use short microwave pulses generated from an optoelectronic parametric oscillator as the spins to implement the Ising machine with large scale and also high coherence under room temperature. The proposed machine supports 10,000 spins, and the high coherence leads to accurate computation. Moreover, the Ising machine is highly compatible with high-speed electronic devices for programmability, paving a low-cost, accurate, and easy-to-implement way toward to solve real-world optimization problems.

and j-th spins, σi and σj denote the z projection of the spins with eigenvalues of either 1 or −1.
Finding the optimal answer is equivalent to searching the ground state of an Ising machine (11), and computing acceleration is obtained by utilizing the system intrinsic convergence property since convergence to the ground state is often at a high speed.
To solve the real-world optimization problems, large-scale, high-coherence Ising machines are highly desired (1)(2)(3)(4)(5). Many Ising machines, e.g. based on trapped ions (12,13), superconducting circuits (14), molecules (15,16), and optical systems (17)(18)(19)(20)(21) have been reported. The key features of these Ising spins are shown in Fig. 1. Although an Ising machine based on trapped ions has been successfully demonstrated with advantages such as high fidelity, long trapping and strong spin-spin couplings, it must operate under an ultra-low temperature of millikelvin and is constrained to small scale networks of about 300 spins due to the difficulties of constructing the requisite optical and electronic control systems (13). A large-scale Ising machine operates under room temperature can be achieved by using molecular spin, but it takes a long time of hundreds of seconds to optimize due to the limited speeds of molecular motions and chemical reactions (16). By the use of quantum superconducting or degenerate optical parametric oscillator (DOPO) spins, large-scale Ising machines have been implemented with computation acceleration. For example, the state-of-the-art quantum annealers based on superconducting circuits can achieve 5000 qubits/spins (22) and computation speedup can be realized thanks to the parallel search process (23,24). However, maintaining the spins coherence requires a temperature near absolute zero and extremely clean electromagnetic environment, which is costly and technically difficult. When spin values are represented by the phase of short optical pulses, Ising machines based on DOPOs can operate at room temperature, and more than 10,000 spins have been reported (20). Moreover, DOPO-based Ising machines can also search the answer parallelly (25) and is superior to quantum annealers (26) and digital computer (27) in dense graphic optimization. Nevertheless, as a long fiber cavity is used to store the large number of spins, it weakens the spin coherence since the optical phase is very sensitive to the cavity delay variation. According to 2 o ft   =  , where   is the phase change of the spin, o f is the optical frequency (around 200 THz) and t  is the time jitter, a few femtoseconds of jitter would reverse the sign of the spin, causing the system Hamilton evolves away from the expectation. In short, the above-mentioned approaches have difficulties in implementing a long-term-operating, large-scale and high-coherence Ising machine under room temperature.
Here, we propose a novel Ising machine, whose artificial Ising spins are represented by the phases of the short microwave pulses generated from an optoelectronic parametric oscillator (OEPO) (28). Similar to the DOPO, the OEPO can also operate at degenerate state with relative phase of 0 or π, but with a much longer wavelength. Since the microwave wavelength (2 cm) is about 4 orders of that of the lightwave (1.5 μm), the microwave phase is much less sensitive to the cavity variation, resulting in high spin coherence. Moreover, by using low-loss, long optical fiber and time-division multiplexing technique, large scale spins can be obtained, proving enough space to code the realworld problems. Figure 2 shows the experimental setup and operation principle of the proposed microwave photonic Ising machine. One may figure out that the fabric of OEPO in Fig. 2A is similar to an optoelectronic oscillator (OEO) where a hybrid microwave photonics feedback loop has also been used. In essence, the conventional OEO is a continuous-wave microwave generator using a lowloss optical fiber as energy storage media to produce high-coherence microwave signals with an uncertain initial phase (29)(30)(31)(32). In our scheme, the microwave phase is locked to a local oscillation (LO1) through electrical second-order degenerate parametric process. If the central frequency and the bandwidth of the microwave filter are properly designed, the generated microwave signal has a frequency half of LO1. In this case, the parametric process is a phase conjugate operator that reverse the phase of the input microwave signal (28). If the frequency of the LO1 is also carefully designed and synchronized to the cavity resonant frequency, the phase conjugate operator would lock oscillating frequency in a relative phase of either 0 or π (see supplementary materials for details). The binary-phase oscillation is then used to simulate the artificial Ising spin, e.g. 0-phase/π-phase oscillation represents up/down spin. By using long fiber combined with time-division multiplexing, large-scale, discrete pulse oscillations are obtained in a single cavity. In such an oscillation network, the spin-spin interaction can be implemented by delay lines (20) or measurement-feedback circuit (11). The delay-line scheme, which is implemented by a wavelength-division multiplexing (WDM) system and optical delay lines (ODLs), is adopted in our demonstration. As shown in Fig. 2B, specific spin-spin interaction, defined by matrix J, would change the global loss of the OEPO network and result in particular phase configuration (see supplementary materials for details). Phase configuration with minimum loss has the maximum possibility to operate (19). By gradually increasing cavity gain, the OEPO network can operate around the minimum-loss state and outputs the corresponding microwave signal whose phase configuration is the optimal or suboptimal solution of the given Ising problem.

Generation of 10,000 spins with high coherence
Large scale and high coherence of the proposed system are verified in a non-interaction (J is a unit matrix) artificial spin network and results are shown in Fig. 3. The oscillation frequency of microwave photonic Ising spins is 10 GHz, and the wavelength is about 2 cm in the optical fiber with a refractive index of 1.47. The spin repetition period is 10 nanoseconds, so that 10,000 spins are supported by using 20-km fiber. Since the ratio between the cavity length and the oscillation wavelength is dramatically decreased from 10 9 (20) to 10 6 , the spin values in the OEPO are much less sensitive to the inevitable temperature fluctuation and ambient vibration, obtaining a much higher coherence compared to the counterpart in DOPO.
Another external microwave source (LO2), whose frequency is the same as the oscillating signal and synchronized to the LO1, is used to demodulate the short microwave pulses to the baseband to extract the microwave phase. By tuning the phase of LO2, the demonstrated baseband pulses can reach the maximum peak-peak value, and the peak value in a single test (each test starts from noise to stable oscillation) is measured and shown in the histogram in Fig. 3A. The absolute value varies around a small range but may have opposite sign, which indicates the spins oscillate with stable amplitude and either 0 or π phase. To further evaluate the spins coherence, 60,000 overlapped waveforms are recorded within 30 minutes by using a high-speed oscilloscope under persistence mode and results are shown in Fig. 3B. The "open eye" in Fig. 3B indicates that spins are highly stable both in amplitude and in phase, while the "eye" would be closed in an unstable system since its amplitude and/or phase varies over time. The high coherence guarantees the system evolves toward the expectation, thus accurate computation can be expected. This is a unique advantage over other Ising machines, such as based on DOPO, quantum superconducting circuits, or trapped ions, because maintaining high spin coherence is challenging in these networks. Since oscillation starts from shot noise and thermal noise, each spin is independent and global phase configuration would be chaotic if no interaction is implemented. In 100 tests, the ratio between the positive and negative spins of each test falls into the range from 0.97 to 1.02 (mean 0.9998), suggesting equal probabilities of 0-phase and π-phase oscillations. Furthermore, we extract the spin sign and calculate the correlation, as shown in Fig. 3C. In the autocorrelation, only one single peak is observed at the zerobit delay, suggesting the spins are independent in each test. And no peak is observed along the whole span in the cross-correlation between two tests, which indicates that each test is independent.

Simulation of 1D/2D Ising model
Two Ising models, i.e., the closed one-dimensional (1D) chain and the 2D square lattice (33), are first used to verify the feasibility of the proposed Ising machine. To implement a 1D Ising simulator, a second channel with an additional 1-bit delay is used to interact the neighboring spins, as shown in Fig. 4A. The formed 1D Ising model is a closed loop and energy coupling is unidirectional (Ji,i+1 ≠0, Ji+1,i=0), from (i−1)-th spin to i-th spin for i≤10,000, and 10,000-th spin to the first spin. Figure   4B shows the evolutions of the spin phases and Ising energy as a function of roundtrip number in positive coupling 1D chain. The blue/yellow dot represents up/down spin. In our experiment, the cavity gain is gradually increased by enlarging the LO1 power. At the beginning, the cavity net gain is far below the oscillation threshold, so that system noise dominates the spin phases and leads to a chaotic configuration. As a result, the domains, namely different regions where spins have the same sign, are relatively short. As the LO1 power increases, cavity gain exceeds its loss, and the signal starts to oscillate along with the adjusting of the spin phases. From the noise to stable oscillation, some domains shrink and disappear, while others grow correspondingly. As the last domain wall (topological defect that separate different domains) falls down, all domains merges into one and all the spins have the same sign. This is because the positive energy couplings between spins forces the neighboring spins to share the same phase, so that the collaborative loss can be minimal. We calculate the frequencies of domains with different lengths as a function of the roundtrip number, and result is presented in Fig. 4C. In the chaotic stage, the short domains (domain length is below 10) have a relatively high frequency while the long domains are rare. As the roundtrip number increases, the spins are getting more and more orderly, and the frequencies of the short domain decrease. As the short domains collapse and merge into longer domains, the frequencies of the long domain increase in a certain period (roundtrip number from 2,000 to 4,000). As evolution continues, long domain continuously grows, swallowing up other short domains. The evolution stops at around 9,000 roundtrips with the spins oscillate stably in 0-phase state and the system reaches the ground state. Noted that the domain structure of the DOPO-based Ising machine with the same spin number would "freeze out" and will not reach the ground state unless one waits for an exponentially long time (33).
As shown in Fig. 5A, another channel with additional 100-bit delay is used to implement the vertical coupling. Combined with the above 1-bit horizontal coupling, we can simulate a 100×100 2D Ising square lattice. The (i−1)-th and (i−100)-th spins are coupled into the i-th spin with boundary conditions that the 10,000-th spin is coupled into the first spin, where i≤10,000, and the (9900+k)-th spin is coupled into the k-th spin, where 1 ≤ k ≤ 100. Thanks to the additional constraints brought by the interaction of the vertical dimension, the machine takes much less time to evolve from the chaotic state to the ordered state compared to the 1D Ising simulator, as shown in Fig. 5C. A much steeper curve can also be observed in the Ising energy evolution. Only tens of roundtrips are taken to reach a low energy close to the ground state, suggesting a clear phase transition between the chaos and the order state. Some snapshots at specific roundtrips are shown in Fig. 5B to reveal the evolution process. Here, the coupling coefficients in two dimensions are both positive and unidirectional. As a result, only one domain remains and all spins fall into the same phase configuration at the end. The frequencies of domains with different sizes as a function of roundtrip number are shown in Fig. 5D. As the roundtrip number increases, the large domain gradually swallows up other domains and the machine takes less than 1,200 roundtrips to the ground state.

Solving a Max-cut problem
Max-cut problem is an NP-hard problem that can be mapped onto Ising formulation (11). Here, we program an unweighted, unidirectional Möbius ladder graph with 20 vertices into our machine, as shown in insect of Fig. 6A. In this demonstration, the 20-km fiber is replaced by a short one, and the 100-bit delay line is also replaced by a 10-bit one while other parameters are maintained. In a graph with 20 vertices, the number of possible combinations is 2 20 . Although the possible Ising energy has only 27 different values which vary from −26 to 30, finding the lowest states is still very difficult since there are only twenty out of 2 20 combinations meet the requirement, as shown in Fig.  6A. By tuning the according ODLs, the desired interaction matrix is implemented. Figure 6B shows the demodulated baseband pulses train. The Ising energy is −26 when the oscillation is stable, suggesting the lowest energy under the given interaction matrix is found. We performed 100 tests and find that lowest Ising energy is obtained in each test. The minimum Ising energy corresponds to the max cut, which is 28. Figure 6C and D present each spin evolutions and the Ising energy/graph cut size as a function of round-trip number. The spins have random phase at the beginning, corresponding to the relatively high Ising energy and small cut. As the round-trip number increases, spins evolve towards certain, either 0 or π phase. Meanwhile, the Ising energy drops down while the cut size increase. At around 60 roundtrips, the machine reaches the lowest-energy phase configuration and the given Max-cut problem is solved.

Discussion
In order to solve combinatorial optimization problems in real-world applications, the Ising machine should be qualified with performances such as large scale, high coherence, and also flexible programmability. Large-scale spin network provides enough space for coding real-world complex problems. High coherence ensures accurate computation of the Ising machine. Programmability makes the machine be common to use for various optimization problems. Thanks to the use of lowloss, long optical fiber in the optoelectronic cavity of OEPO, a large-scale Ising network can be easily implemented. Since the spins in the proposed machine are represented by short microwave pulses, which are less sensitive to the environment disturbance and maintain high coherence, the machine evolves fast toward the ground state. Moreover, flexible programmability can also be easily realized by utilizing measurement-feedback schemes, since the microwave spin can be easily measured and controlled with high-speed electrical devices (21). In a word, the microwave photonic Ising machine promises an optimizer with a large scale, high coherence, and programmability, paving a way to solve real-world problems.

Experimental details
A multi-channel laser (IDPhotonics, CoBriteDX4) is used in our experiment. In the noninteraction oscillation, only one channel is used. In the simulation of 1D Ising model, two channels are used. In 2D Ising model simulation and solving the Max-cut problem, three channels are used. These channels have a 100-GHz frequency spacing with the wavelength from 1549.32 nm to 1550.92 nm. The continuous lightwave from the laser are coupled together through a WDM and then shaped into optical pulse trains by using electrical pulse trains and a low-biased intensity modulator. The electrical pulse trains are generated from a programmable pulse pattern generator (PPG, Keysight N4960A) with a repetition rate of 100 MHz and a duty cycle of 20%. The generated optical pulse trains are then used as the carrier of the oscillating microwave signal in the optoelectronic cavity. A multiplexer and a demultiplexer are used in the cavity to enable independent delay control for each channel. Three tunable ODLs are inserted between the multiplexer and the demultiplexer. The tunable range of the ODLs is 500 ps, which is much larger than the period of the oscillating microwave signal (100 ps). As a result, the sign of coupling coefficients in interaction matrix J can be easily adjusted, while their magnitudes are also alterable by tuning the laser power. To minimize the dispersion, a 20-km DSF is used to store the microwave spins. The PD (DSC40S) in our experiment has 3-dB bandwidth of 18 GHz and responsivity of 0.75 A/W at 1550 nm. The signal after PD is filtered by an 8-12 GHz BPF and amplified by an electrical amplifier, then is divided into two parts: one is used for measurement and another one is launched into an electrical mixer and interacts with an external oscillation (LO1) with 20-GHz frequency and 10-dBm power. A second 8-12 GHz BPF is used to suppress the sum-frequency signal and leaked LO1. The filtered IF signal is then amplified and feedback to a quadrature-biased intensity modulator (IM). The LO1 frequency LO1 , the carrier frequency of the microwave pulses s , and the free spectrum range (FSR) of the optoelectronic cavity FSR satisfy LO1 = 2 s = 2 FSR , where M is an integer.
The microwave pulses (Ising spins) oscillate from noise in the OEPO with random phase, then iterate in the cavity with increasing amplitude and their phases converge to either 0 or π state. Since the signal bandwidth is about 1 GHz, which is much less than the channel frequency spacing, the lightwave interference between different channels has no contribution to the microwave pulses. Another external microwave source (LO2), which is synchronized to LO1, with a frequency of      The cavity net gain is calculated to be 2 dB. The results of ten spins are shown in Fig. S1. After less than 30 roundtrips, the oscillation amplitude reaches the maximum and then stays stable, as can be seen in Fig.1SA. Different from the monotonic increasing in the amplitude, the spin phase jumps back and forth at the beginning because of the phase conjugate operator, as shown in Fig. S1B. As the roundtrip number increase, the wobbles become smaller and finally the spin phases converge to either 0 or π. The stability of the oscillation is also evaluated by the calculating the jitter and results are shown in Fig. S1C. One can conclude that the spin is highly stable and jitters drops down to around 10 -6 after 100 roundtrips. The cavity net gains are around 2 dB at the beginning and decrease to one when the spin amplitudes are stable.

Mapping the Ising problem onto the OEPO network
In an N-spin, non-interaction oscillation network, assuming each of spin has a normalized loss of one, we obtain the global loss N. If the i-th spin is injected by the j-th spin with a coupling coefficient of Ji,j, the loss of the i-th spin would be changed. If the sign of Ji,jσiσj is positive/negative, the loss would decrease/increase from 1 to 1−Ji,jσiσj. As a result, all the couplings from other spins have an integrated loss contribution of 1, , and the total loss of the i-th spin is then given by that we can conclude that the potential phase configurations in an N-spin oscillation network under the given interaction matrix J corresponds to the global loss and maps onto the Ising energy landscape of a given Ising model. Based on the minimum-loss principle, the oscillation network most likely operates at this state. To solve an Ising problem, one can program the corresponding matrix J into the oscillation network, and then gradually increase the cavity gain to search the minimum-loss state. As the cavity gain exceeds its loss, oscillation starts with an increase in amplitude and adjustment of phase. When the oscillation network runs stably, one can measure the corresponding phase configuration, from which the given problems can be solved.

Other experimental results
The temporal waveforms of the oscillating microwave pulses and the demodulated baseband pulses are shown in Fig. S2A. The microwave pulse lasts several periods. In the non-interaction oscillation, 100 tests are performed. The ratios between the numbers of the positive and negative pulses are calculated and the results fall into the range from 0.97 to 1.02, as shown in Fig. S2B. The power spectrum of the non-interaction oscillation is also measured by an electrical spectrum analyzer (ESA, Keysight N9030A) and the results are shown in Fig. S3C. One can find that the power spectrum is actually a microwave frequency comb whose frequency spacing equals to the FSR of the optoelectronic cavity. Clear comb line under high-resolution bandwidth (RBW) suggests the microwave signal is highly stable, which corresponds to high spin coherence.
The temporal waveform of the demodulated baseband pulses in the 1D positive coupling Ising simulation is shown in Fig. S3A. Since the energy couplings are positive, all the spins share the same value. The period of the signal is 10 ns, which equals to the interval of two neighboring spins. This results in a much sparser comb line in the frequency domain, as shown in Fig. S3B and C.