a A four-spin coupled optical Ising machine, consisting of four on-chip OEC nonlinear units and a 4 × 4 on-chip symmetric coupling matrix. Once the coupling matrix is deployed, the nonlinear coupled oscillators evolve from an initial state to the ground state of the Ising problem. b On-chip symmetric coupling matrix, comprising an optical matrix connected end-to-end with 3 dB couplers. The optical matrix, highlighted in orange, is constructed using a single-drive MZI in a rectangular unitary mesh configuration. c On-chip OEC nonlinear unit, where the differential photocurrent drives the forward-biased MRM. d The nonlinearity of the output optical power as a function of the photocurrent. e Diagram of the electric current flow in the differential photodetectors.

Figure 1b illustrates the details of the on-chip symmetric coupling matrix, which can be considered as multiple coupled Sagnac interferometers. In this system, input light is split into forward and backward paths within an optical mesh, where the light converges after completing a round-trip. Due to the reciprocity of optical linear systems, the forward and backward transmission matrices of the optical mesh are mutually transposed²⁸. We can derive the output optical fields from Ports 1 and Ports 2 as:

where P and P^T are the forward and backward complex amplitude transmission matrices of the MZI mesh, respectively, and X is the complex amplitude of the input optical vector. If the incident optical signals are mutually incoherent (e.g., operating at different wavelengths), the output light intensity will be the linear superposition of the input light intensities, and the intensity transmission matrix can be written as ref. ²⁵:

where |A|² denotes a matrix formed by the modulus squared of each element of matrix A. It is evident that both M₁ and M₂ are positive symmetric matrices. To construct a real-valued symmetric matrix, we can express the coupling matrix as:

Thanks to the inherent symmetry of the proposed optical matrix, its configurational accuracy and robustness to perturbations are superior to conventional one-way optical meshes. The optical matrix (highlighted in orange) in our design is based on an (N + 1) × (N + 1) rectangular unitary MZI mesh with a single phase shifter in each MZI (right side of Fig. 1b; see Supplementary Information S1 for derivations)^25,29. The first N out of the N + 1 ports of the MZI mesh are designated as the input and output ports of the coherent optical matrix. The total number of phase shifters is N(N + 1)/2, which satisfies the minimum degree of freedom required to achieve an N × N real-valued symmetric matrix. The reduced number of phase shifters significantly lowers the energy consumption needed for configuring and maintaining the optical matrix.

Figure 1c illustrates the detailed structure of the on-chip optical-electrical coupled (OEC) nonlinear unit. Optical signals from the on-chip coupling matrix enter the germanium/silicon photodetectors, which are arranged in a differential configuration (I₁-I₂). The identity matrix in Eq. (6) is realized by splitting 50% of the OEC nonlinear output into another differential signal (I₃). The total differential photocurrent (I₁-I₂-I₃) directly drives the MRM, producing the corresponding nonlinear response. The intensity of the supply light determines the strength of the feedback (α). Figure 1d shows the working principle of the nonlinear transfer function. When the injected current I_n > 0, the MRM is forward-biased, and the transmission spectrum experiences a significant blue shift with increasing I_n. Conversely, when I_n < 0, the MRM is reverse-biased, and the spectrum shifts slightly to the red due to the much lower carrier-depleted modulation efficiency³⁰. The nonlinear response is akin to a linearized sigmoid function. Although this function is not commonly used in optical Ising machines, we theoretically prove that it can effectively approximate the Ising problem-solving process (see Methods and Supplementary Information S2 for the capabilities of the general OEC nonlinear transfer function). In this case, the normalized output optical intensity of the OEC unit indicates the state of the corresponding spin, with a value of 1 representing ‘spin-up’ and 0 representing ‘spin-down’. Figure 1e depicts the electric current flow within the OEC nonlinear unit, where a bias current I_b provides a horizontal shift of the nonlinear transfer function.

Numerical evaluation of the proposed optical Ising machine

To highlight the advantages of the proposed optical symmetric matrix architecture, we perform a comparison of the matrix fidelity for the two structures shown in Fig. 2a. The optical MZI mesh on the right side of Fig. 1b is selected as the reference optical matrix highlighted in orange for both cases. The positive and real-valued transmission matrices are calculated using direct detection and differential detection, respectively²⁵. For 1000 randomly generated positive and real-valued symmetric matrices, we optimize the phase shifters to maximize the correlation between the current physical matrix and the target matrix (see “Methods”). Figure 2a presents a histogram of the optimized correlations. The mean correlation of the proposed optical mesh exceeds 0.9999 for both positive and real-valued symmetric matrices, while the conventional one-way optical mesh exhibits significantly lower fidelity. The imperfect fidelity of the proposed optical mesh is partly due to the non-universality of the mesh and the local convergence of the optimization algorithm. In this simulation, all MZIs and phase shifters are ideal with high tuning accuracy. However, considering the finite programming accuracy and system noise, the experimental results show slightly worse fidelity than the simulation. We further validate the performance of the OECO Ising machine with 100 spins and a linearized sigmoid function as an approximate nonlinearity. First, we set β = 0, α = 1.1, and τ = 0.1 to verify its bifurcation capability. As shown in Fig. 2b, the 100 spins evolve from their random initial states to either ‘spin-up’ (1) or ‘spin-down’ (0) states, satisfying the first basic condition of the Ising machine. For an initial intensity I_init > 0.5, the spin evolves to ‘spin-up’; otherwise, it evolves to ‘spin-down’. Next, we examine the performance of the OECO Ising machine in solving a 100-spin Ising problem for a random graph g05_100.5, included in the Biq Mac Library, whose ground energy is known to be -397³¹ (see Supplementary Information S3 for tests on other graphs). Figure 2c shows the evolution of Ising energy and spin intensity when α = 0.23 and β = 0.21. After 400 roundtrips, the Ising energy reaches the ground state and stabilizes. It is noted that the spin intensity will not always evolve to exactly 0 or 1 in a physical Ising machine. Therefore, a threshold of 0.5 is used to determine the final state of the spin. Finally, we investigate the influence of α and β on the success rate of the Ising machine. The sweeping step for both parameters is 0.01, and the success rate is counted for 100 random initial states. Figure 2d shows the sweeping results, where the maximum success rate exceeds 0.2. Compared to other nonlinearities, the success rate of our model is nearly half, which is attributed to higher amplitude inhomogeneity. This can be suppressed through engineering of the nonlinear response²⁷ (see Supplementary Information S4 for the influence of noise).

Fig. 2: Numerical investigation of the proposed optical Ising machine.

a Comparison between the proposed paradigm and the conventional approach for constructing the symmetric matrix through numerical simulations. The optical matrix, highlighted in orange, is the same in both setups and is shown on the right side of Fig. 1b. The final correlation of 1,000 random symmetric matrices is tested and analyzed in the histogram. b Evolution of the proposed OECO system above the bifurcation point (α = 1.1), where the system evolves into either the ‘spin-up’ (1) or ‘spin-down’ (0) state. The parameter I_init represents the initial intensity of the spin. c Evolution of Ising energy and spin amplitude for α = 0.23 and β = 0.21. The test graph is g05_100.5 from the Biq Mac Library. d Success rate of g05_100.5 under a sweep of α and β.

Experiment investigation of the proposed optical Ising machine

The OECO Ising machine chip is fabricated in the standard silicon photonics foundry, with a silicon layer of 220 nm in thickness. A microscope image of the chip is shown in Fig. 3a. The germanium/silicon photodetector and MRM are adopted from the process design kits, and they exhibit typical 3 dB bandwidths exceeding 40 GHz and 50 GHz (in carrier-depleted mode), respectively. Figure 3b shows the packaged chip with optical fiber array and electrical wire-bonding. Due to the forward-biased operating mode of the MRM, along with low-frequency packaging and power supply limitations, the actual operating bandwidth of the OEC unit is approximately 100 MHz, resulting in a corresponding latency of 1.6 ns (see Supplementary Information S5). The optical feedback length is approximately 8 mm, corresponding to a time constant of τ ≈ 0.11 ns, which is much shorter than the latency determined by the system’s bandwidth. However, it has been shown that both latencies can be considered together as a single effective latency without compromising system functionality (see Supplementary Information S6). Figure 3c displays the measured response of the photodetector under varying bias voltages, where the linear responsivity is 0.8 A/W. Figure 3d provides a close-up photograph of the OEC nonlinear unit, where the bias voltage V is set to 3 V to ensure the photodetector operates in the linear region. Figure 3e shows the transmission spectrum of the differentially driven MRM, excited by incremental differential optical power. The transmission spectrum for zero input is highlighted separately, and the corresponding nonlinear transfer function at the zero-input resonant wavelength is measured, as shown in Fig. 3f, which aligns well with our theoretical predictions. The optical power required to excite the spin state from 0 to 1 is less than 1 mW, indicating that the on-chip OEC unit is more energy-efficient compared to other nonlinearities used in previous Ising machines.

Fig. 3: Chip fabrication and basic characterization.

a Microscopic photograph of the optical Ising machine. b Photograph of the packaged optical Ising machine chip. c Photocurrent of the on-chip photodetector under different bias voltages. d Zoom-in photograph of the OEC unit. e Drift of the transmission spectrum of the MRM with increasing differential optical power. f Measured nonlinear transfer function of the OEC unit. P_diff denotes the differential input power, and P_out denotes the output power. g Evolution of matrix correlation given a random target positive-valued symmetric matrix. h Evolution of matrix correlation given a random target real-valued symmetric matrix.

The fidelity of the on-chip symmetric matrix is experimentally evaluated in the next step. We randomly select one positive and one real-valued symmetric matrix as targets and configure the on-chip optical mesh using a gradient descent algorithm (see “Methods”). The iteration curves for correlation are shown in Figs. 3g and h. The final optimized correlations are 0.9992 and 0.9979 for the two cases, respectively, while the mean correlations for positive and real-valued symmetric matrices are 0.998 and 0.986, as detailed in Supplementary Information S7. These results demonstrate that the proposed symmetric optical matrix can effectively represent the necessary symmetric matrix in the optical Ising machine.

The nonlinear coupled oscillators start with all spins in the ‘spin-down’ state. Therefore, it cannot automatically evolve to the ground state without additional excitations to provide an appropriate initial state. To emulate an arbitrary initial state, four independent channels of excitation signals would be necessary. In the subsequent experiment, only one excitation optical pulse is injected into the positive channel of the first OEC nonlinear unit, as shown in Fig. 4a. The excitation can also be integrated on-chip to facilitate testing of arbitrary initial states. In our setup, the supply light intensity is set to 4 dBm, while the maximum intensity of the excitation light is -6 dBm. The excitation pulse intensity determines the initial spin state. After the excitation pulse is withdrawn, the system begins its evolution towards the ground state from the corresponding initial state.

Fig. 4: Experimental verification of bifurcation in the proposed optical Ising machine.

a Experimental setup of the optical Ising machine, where the initial states of all spins are ‘spin-down’. A positive excitation pulse is sent to one of the OEC units, triggering the physical evolution of the optical Ising machine. b Influence of the nonlinear function (voltage applied to the thermal phase shifter of MRM), feedback strength (intensity of supply light), and initial state (intensity of excitation pulse) on the oscillation of the nonlinear cavity. The olive zones represent the working voltages where the cavity can oscillate under the current initial state and feedback strength. BP denotes the bifurcation point. c Bifurcation of the nonlinear oscillator under random excitation signal intensity. d Time response of the oscillator with a 40 kHz excitation signal and a 0.1% duty cycle. The upper and lower waveforms correspond to strong and weak excitation, respectively.

To start, we train the on-chip MZI mesh to act as an identity matrix to examine the bifurcation of a single OECO, with only the supply light of the first OEC nonlinear unit turned on. The control parameter for bifurcation is the self-feedback strength α, which corresponds to the power of the supply light. Above the bifurcation point, the spin will evolve to either ‘spin-up’ or ‘spin-down’, determined by the initial state. Below the bifurcation point, the spin remains in a fixed point regardless of the initial state. In addition, the bifurcation point is also influenced by the voltage applied to the thermal phase shifter of MRM, which controls the horizontal shift of the nonlinear function (see Eq. (8) in Methods).

We modify the working point of the MRM by sweeping the voltage, as shown in Fig. 4b. The sweeping curve shows an abrupt drop followed by a slow rise, attributed to the thermal effect within the ring. This occurs at high supply light power to generate a sufficient gain coefficient (α). The working region with lower voltage is above the bifurcation point, while the region with higher voltage falls below it. When an excitation pulse is applied, the MRM, which operates in the shadow region (0.58 V to 0.69 V), evolves to ‘spin-up’. As the intensity of the excitation pulse decreases, the working region (0.6 V to 0.69 V) eventually evolves to ‘spin-down’, consistent with the simulation results in Fig. 2b. Increasing the intensity of the supply light widens the working region above the bifurcation point.

Next, we select a suitable operating point for the MRM (0.63 V) to balance the probabilities of the binary states during bifurcation. We then input a sequence of excitation signals with random intensities to observe the bifurcation process. It is important to note that once the spin evolves into the ‘spin-up’ state, subsequent positive excitation pulses maintain the oscillation in the cavity, resulting in no change to the spin output intensity, regardless of the excitation pulse intensity. To clearly demonstrate the bifurcation, we reset the MRM to the ‘low-pass’ state before each new excitation pulse. This is achieved by applying a large voltage (0.8 V) to the thermal phase shifter and then returning it to 0.63 V. In practical Ising problem-solving, the reset operation is not restricted to the method described above and can be adapted based on the need for repeated testing. Figure 4c shows the captured waveform of the OEC output signal, where the time span between the adjacent excitation and reset pulses is about 50 ms. This long interval is chosen to allow for the observation of the binary states as the spin intensity stabilizes. The spikes in the waveform are attributed to the excitation and reset pulses. The deterministic evolution towards either ‘spin-up’ or ‘spin-down’ meets the basic requirements of the Ising machine.

To measure the characteristic time of spin evolution, we use excitation and reset pulses with a frequency of 40 kHz, a duty cycle of 0.1% (corresponding to a 25 ns pulse width), and rise-fall times of 48.8 ns. Figure 4d shows the OEC output signal waveforms for strong (upper panel) and weak (lower panel) excitation pulses, respectively. Upon withdrawal of the excitation pulse, the spins in the upper and lower panels evolve to the ‘spin-up’ and ‘spin-down’ states, respectively. These two cases correspond to the high and low initial intensity (I_init) scenarios from the simulation in Fig. 2b. The evolution time is approximately 150 ns, indicating that the spins take only a few round trips to reach the steady state. The waveforms of the excitation and reset pulses are provided in Supplementary Information S8.

Finally, we solve a benchmark Ising problem with four coupled spins using the optical Ising machine chip. As shown in Fig. 5a, Spins 1 and 3 are positively correlated, as are Spins 2 and 4, while Spins 2 and 3 are negatively correlated. Using the coupling matrix in Fig. 5a and selecting α = 0 (i.e., setting the diagonal elements of the on-chip matrix to zero), we train the on-chip symmetric matrix, as depicted in Fig. 5b. The excitation signal here is injected into the second channel of spins to demonstrate the versatility of our chip. The excitation intensity is kept consistent with that of previous experiments, ensuring that the final state of the second spin is 1. Under this excitation, the system evolves to the ground state of the Ising problem, which is known to be [-1 1 -1 1]. As in the previous setup, excitation and reset signals are alternately applied to the system with a time span of 10 ms. Figure 5c shows the evolution of the four OEC nonlinear units. The four spins successfully evolve to the ground state through physical annealing. In Fig. 5d, the evolutions of the four spins over time are shown together. Spins that end up in the same state have different final intensity values, partially due to the amplitude inhomogeneity of the Ising machine and slight differences in the power of the supply light. In addition, some of the observed waveform oscillations may be attributed to microwave oscillations occurring within the OECO feedback loop. The four reset and one excitation signals are supplied by low-speed electrical sources due to our limited channel of high-speed signals. Therefore, the evolution time of several microseconds is partly limited by the establishing time of the excitation signal (see Supplementary Information S8). To further demonstrate the versatility of the fabricated Ising machine chip, we solve two additional non-trivial Ising problems on-chip (see Supplementary Information S9). The results confirm that the proposed optical Ising machine is well-suited for solving a variety of Ising problems (The success rate is tested in Supplementary Information S10).

Fig. 5: Experiment demonstration of four-spin Ising problem solving.

a The Ising problem with its corresponding coupling matrix. b The actual coupling matrix deployed on the optical matrix. c Evolution of the four spin states under a specific excitation. d Zoomed-in view of the evolution process of the four spins.

Discussion