**Operation principle and experimental setup**

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 low-loss optical fiber as energy storage media to produce high-coherence microwave signals with an uncertain initial phase (*29-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 109 (*20*) to 106, 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 zero-bit 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 (*J*i,i+1≠0, *J*i+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 220. 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 220 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.