Solving combinatorial optimization problems involves finding the maximum or minimum of an objective function under certain constraints, which are common in fields such as science, engineering, and business.1^–3 These problems are categorized as non-deterministic polynomial-time (NP)-hard or NP-complete problems, meaning their computation time scales exponentially with problem size,4 making them intractable with modern computers.5 Mathematically, NP problems can be reduced to one another in polynomial time,1 which implies that solving one NP problem efficiently can potentially provide solutions to others. Therefore, proposing an efficient solution method for one class of NP problems can enable algorithmic mapping to other problems. The quest for efficient solutions to NP problems has driven the development of computing accelerators, both digital and physical.6^–10

The Ising model is a well-known NP-hard decision model, widely used to study complex dynamics in fields such as physics, computer science, biology, and social systems.11^–14 The Ising problem involves finding the configuration of binary particles (spins) that minimizes the energy of the Ising Hamiltonian, where particles interact with each other. Finding the ground state of the Ising Hamiltonian can be used to solve other combinatorial optimization problems through polynomial-time mapping.15^–17

As Moore’s law for conventional computers slows, there has been a surge in the development of unconventional computing architectures based on various physical systems for solving Ising problems, including optics10^,16^,18^–32 and magnetic devices.33^–40 Optical systems offer advantages such as parallelism, scalability, and low power consumption,41^–47 making them ideal for the development of an optical Ising machine (OIM) that can simulate Ising problems, potentially outperforming traditional von Neumann architectures.

Several OIM schemes have been proposed, including those based on spatial and temporal domains. The time-domain OIM typically relies on degenerate optical parametric oscillators (DOPOs) and long-distance optical fibers.15^,20^,48^–51 However, this scheme is implemented in long optical fibers, resulting in low integration density. In addition, DOPOs require high-power light sources for pumping, leading to high power consumption. The spatial-domain OIMs, which utilize spatial light modulators (SLMs) for encoding, offer better scalability and are more suitable for constructing large Ising networks.26^,52^–59 However, they face limitations in programming flexibility due to the restricted number of grid points, and the refresh rate of the SLMs limits computational speed. Furthermore, the architecture of discrete spatial computing systems tends to be bulky due to their reliance on discrete optical components, leading to higher operational costs. Integrated spatial-domain OIM using Mach–Zehnder interferometers (MZIs) shows some promise, but the number of MZIs needed increases quadratically with the number of spins, making them less suited for large-scale integration.60^,61

Here, we propose a wavelength-domain OIM, where different wavelengths represent different spins. This approach holds significant potential for achieving Ising simulations in an incoherent manner. Wavelength-domain OIMs are expected to leverage the wavelength resources of light, offering improved scalability and reduced experimental complexity due to the incoherence of the wavelengths. Our integrated wavelength-domain OIM has the potential for large-scale scalability and exhibits low computational complexity.

An Ising machine is a physical system designed to simulate Ising problems, which involve multi-spin systems with pairwise linear interactions. In this system, each spin is represented by a binary variable σj∈{+1,−1}, where σj indicates the state of the spin. The spin configuration refers to the arrangement of all spin states in the system. The total energy of the system is determined by the interactions between spins and the influence of an external magnetic field h. This total energy is known as the Ising Hamiltonian16H(σ)=−∑i≠jJi,jσiσj−h∑l=1Nσl.(1)

In the equation, σj and σj represent the different spin states, and Jij denotes the interaction strength between the corresponding pair of spins. This interaction strength can be expressed as an interaction matrix. In the absence of an external magnetic field, the Ising Hamiltonian can be simplified as H(σ)=−∑i≠jJi,jσiσj.(2)

According to the Ising Hamiltonian equation, the key aspects of simulating the Ising problem involve mapping the spins, Ising Hamiltonian, and interaction strengths. As schematically illustrated in Fig. 1(a), the wavelength-domain OIM is based on cascaded microrings (MRRs). The process of solving Ising problems with the proposed Ising machine consists of two main components: problem encoding and optimization. The spin variables are encoded in the resonant states of the MRRs, which are controlled by thermo-optic phase shifters. The Ising Hamiltonian is obtained by measuring the intensity of the output. The coupling coefficients among spins are achieved through intensity modulation of the input signal. Two different encoding schemes are proposed to accommodate both arbitrary coupling types and all-to-all connected coupling types.

The first encoding scheme, which corresponds to the arbitrary coupling type, is illustrated in Figs. 1(b) and 1(c). Based on Eq. (2) of the Ising Hamiltonian, each term represents a coupling between two spins, involving a coupling coefficient Jij and two binary spins. The expression of the Ising Hamiltonian can then be rewritten in a more intuitive form H(σ)=−∑n=1NJin,jnσinσjn.(3)

Equation (3) provides a clearer correspondence with Fig. 1(b), which shows the on-chip structure of Fig. 1(a). The two rows of MRRs in Fig. 1(b) are arranged in a one-to-one manner. For example, the MRRs in the nth column correspond to the same computing wavelength λn. The resonant states of the two MRRs correspond to the spins σin and σjn. The interaction strength Jin,jn is encoded by the intensity of the input light at wavelength λn. Therefore, the computing light passing through the MRR array corresponds to the summation term in the Ising Hamiltonian. The differential detection of the output corresponds to the Hamiltonian Outarb=∑n=1NJin,jnσinσjn=−H(σ).(4)

For an Ising problem with N coupling terms, 2N MRRs and N wavelengths are required for encoding. Figure 1(c) illustrates the correspondence with Fig. 1(b), where MRRs associated with the same computing term are represented in the same color. In the coupling diagram of Fig. 1(c), solid lines of corresponding colors indicate these couplings. By controlling the mapping between MRR pairs and spins, different couplings can be realized, enabling a programmable Ising model.

The second encoding scheme is applicable to Ising problems where coupling exists between any two spins, known as the all-to-all connected coupling type, as illustrated in Figs. 1(d) and 1(e). In this case, all the MRRs operate at different wavelengths, and the binary spin states are represented by the resonant and non-resonant states of the MRR. MRRs in Fig. 1(d) with different colors indicate different operating wavelengths, corresponding one-to-one with the colored spins in Fig. 1(e). The input is encoded with the intensity εj, determined by the coupling coefficient, and depicted in Fig. 1(e) by different line widths. The output, obtained through coherent detection, can then be expressed as Outall=∑jεjσj.(5)

The detection results are squared in the electrical domain, resulting in cross terms related to σ: Outall2=∑i(εiσi)2+2∑i,jεiεjσiσj.(6)

As the spin state can only be +1 or −1, the first term above remains constant, independent of the spin state. When encoding the intensity, the coupling coefficient Jij=−2εiεj, resulting in the following equation: Outall2=∑i(εiσi)2+H(σ).(7)

As a result, this processing outcome differs from the Ising Hamiltonian by a constant term, which means it can be used directly to evaluate the optimization process.

The two encoding schemes mentioned above each have distinct characteristics and are suited for different types of problems. The arbitrary coupling type offers high flexibility, as it allows for coupling between any pair of spins or the absence of coupling altogether. This makes it particularly well-suited for solving sparse Ising problems. On the other hand, the all-to-all coupling type utilizes N MRRs to implement an Ising machine with N spins. This approach reduces the number of required MRRs, making it more efficient for solving specific problems that demand full connectivity.

The optimization algorithm is used for iterative computation to find the spin configuration that minimizes the Ising Hamiltonian, which corresponds to the solution of the Ising problem. In this context, the simulated annealing (SA) algorithm is employed. A key feature of the SA algorithm is its ability to avoid local optima through stochastic search, ensuring convergence to the global optimum, i.e., the system reaches the lowest energy state. The main parameters in the SA algorithm include the initial temperature (Tini), final temperature (Tfin), and cooling rate (α). Together, these parameters determine the probability (P) of the system accepting changes during the iteration process of the optimization procedure:62P={1,H(n+1)<H(n)e−H(n+1)−H(n)T,H(n+1)≥H(n).(8)

In Eq. (8), H(n) represents the Ising Hamiltonian at the nth iteration state. The Metropolis algorithm, which accepts new states probabilistically, forms the basis of the SA algorithm. The SA algorithm, based on the Metropolis method, can be described as62

Through iterative optimization of the spin configuration, the Ising Hamiltonian gradually decreases.

Figure 2 illustrates the detailed design of the cascaded array consisting of 32 MRRs. Increasing the number of MRRs allows for a larger Ising problem to be solved on the chip. The radius and the coupling region of each MRR are carefully designed to ensure accurate optical analog computing in the C-band. All MRRs in the array share the same structural parameters, as shown in Fig. 2(c), with a zoomed-in micrograph provided in Fig. 2(b). Thermo-optic phase shifters, which utilize TiN heaters, are employed to fine-tune the MRRs. Figure 2(d) shows an overall photo of the packaged layout, which includes both wire-bonding and vertical grating coupling for electrical and optical input/output (I/O). The voltage source is a specially designed field-programmable gate array circuit equipped with digital-to-analog converters, enabling programmable voltages with 16-bit resolution. To control the environmental temperature, a thermo-electric cooler (TEC) is mounted below the chip. In addition, a wavelength-selective switch (WSS) is used for intensity encoding and filtering of the input light. The filtered input spectrum and the through-port spectrum of the MRRs after the preset voltage are shown in Fig. 2(e). The insertion loss introduced by the WSS ranges from 2 to 6.5 dB, with an insertion loss uniformity of less than 1.5 dB. Utilizing sequential channel scanning, the total time required to configure all 32 input wavelengths is less than 64 ms. During the experiment, only minor adjustments are required for the MRR voltages to achieve the resonant state. This experimental method helps minimize the impact of thermal cross-talk in the MRRs. Based on the parameter of MRR, each MRR requires an average tuning power of ∼20mW to achieve the desired phase shift corresponding to a single spin. Given a typical tuning duration of 1 ms, the energy consumption per MRR of each iteration is ∼20μJ. Therefore, for solving an Ising problem with 32 spins over 100 iterations, the energy budget is ∼64mJ.

Light from an external broadband light source is filtered and intensity-encoded using a WSS to generate multi-wavelength optical signals with a bandwidth of 0.5 nm, which serves as computing light. These optical signals are transmitted through the MRR array and detected by a multi-channel photodetector, where the differential output is calculated. The computational results are stored and processed in the electrical domain, with the SA algorithm iteratively updating the resonance states of the MRRs. The iterative computation is carried out throughout the experimental process. Once the SA process is complete, the set of resonance states of the MRRs, stored in the electrical domain, corresponds to the spin configuration that represents the solution of the Ising problem. Based on the experimental setup and control scheme used in this work, the encoded Ising problem can be solved within 2.5 ms.

The experimental validation is divided into two sections based on different coupling methods. First, we discuss the validation of the Ising machine with arbitrary coupling type, which offers high flexibility and allows for the encoding of sparse coupling matrices without excessive resource usage.

MAX-CUT is a problem that involves dividing the nodes of a given graph G into two subsets {S1,S2}. The coefficient w represents the adjacency matrix of the graph, which defines the coupling among the nodes. The solution to the MAX-CUT problem is the partitioning that maximizes the sum of the weights wi,j(i∈S1,j∈S2) of the edges connecting the two subsets. The CUT value C(σ) can be defined as C(σ)=∑i∈S1,j∈S2wij.(9)

By representing the subset membership of node j with σj=±1, the state of the edge corresponding to node j and i can be expressed as 1−σiσj2. The value is 1 when σj and σi are not equal, indicating that the edge belongs to the CUT. Conversely, the edge does not belong to the CUT and does not contribute to the CUT value calculation when σj and σi are equal. Thus, the CUT value can be expressed as C(σ)=∑i<jwij(1−σiσj)2=12[∑i<jwij−H(σ)].(10)

The coupling coefficient Ji,j is defined as −wi,j. In Eq. (10), H represents the Ising Hamiltonian, as defined in Eq. (2). By correlating the adjacency matrix with the coupling coefficient matrix of the Ising problem, the MAX-CUT problem can be solved using an Ising machine. The CUT value is maximized when the Ising Hamiltonian is minimized. The MAX-CUT problem is commonly used to address various issues, such as structural optimization in communication and computer networks, user group behavior analysis in social networks, and image processing in computer vision. The spin coupling patterns in these practical problems are flexible, necessitating the design of simulators with adaptable encoding schemes. The wavelength-domain OIM we proposed can solve the MAX-CUT problem. In our experiment, we observe the differential output from the through and drop ports of the MRR array. The spin coupling scheme we use corresponds to the first type of encoding scheme outlined in Sec. 2. It is implemented by expanding the summation equation of the Ising Hamiltonian as expressed in Eq. (2). Each input wavelength, with a specific intensity, corresponds to one of the additive terms in the Ising Hamiltonian. This encoding scheme offers flexibility and is a significant improvement over previous simulated Ising models, which have limitations in their coupling schemes. It supports arbitrary coupling types, making it particularly well-suited for sparse coupling matrices.

As shown in Fig. 3(a), a graph with six nodes was selected for the experiment, which includes eight coupling terms corresponding to eight input wavelengths. Depending on the sign of the coupling coefficient, positive and negative input ports are used to perform on-chip Ising Hamiltonian computation. As each coupling term involves two MRRs, a total of 16 computation MRRs are used. After power normalization, the differential detection values from the two output ports of the MRR array yield the Hamiltonian, which is then employed for SA optimization. In each iteration, the node configurations at intermediate states are used to calculate the CUT value theoretically, as shown in Eq. (9) and Fig. 3(c). Using the proposed wavelength-domain OINM, a solution to the MAX-CUT problem is successfully obtained, with the corresponding Ising Hamiltonian of −4 and the CUT value of 4. The node classification at this stage is shown in Fig. 3(b). Multiple node assignments result in the maximum CUT value of 4, and the solution will converge to one of these assignments. This experimental result demonstrates that the proposed on-chip OIM can solve Ising problems with arbitrary coupling types and is effective in solving the MAX-CUT problem.

In the Ising machine with the second type of coupling, there is a connection between every pair of distinct spins. Using the on-chip wavelength-domain method we proposed, N spins can be encoded with N MRRs. To demonstrate the problem-solving capability of the Ising machine, we use it to address the number-partitioning problem (NPP). The NPP has a time complexity of O(2N), and its difficulty increases with both the length of the set and the size of the elements when solved using an exhaustive approach. The NPP involves partitioning a set S into two subsets, S1 and S2, such that the sum of the elements in S1 is equal to the sum of the elements in S2. The difference between the sum of two subsets can be expressed as Diff=(∑i=1Nniσi)2.(11)

In Eq. (11), ni represents an element in the set, and σi∈{1,−1} denotes the subset to which the corresponding element belongs. Therefore, the solution to the NPP is transformed into finding the configuration of {σ} to minimize Diff. We encode these elements into input optical signals using intensity encoding, such that ni=εi. Equation (11) has the same structure as Eq. (6) and can be transformed into an Ising problem for solution. The elements in the set are encoded in the coupling coefficients. The resonant and non-resonant states of the MRRs represent different subsets of the problem. The Ising Hamiltonian, as shown in Eq. (7), can be detected from the output. During each iteration, we record the current set partition state and compute the difference in the subset sum. The NPP is then solved using the SA algorithm, which leverages the OIM we proposed. In the experiment, we focus on the NPP with a set of 32 elements. The resulting partition of the set is shown in Fig. 3(d). The iterative process of the Ising Hamiltonian in the experiment is illustrated in Fig. 3(e).

We also demonstrate the scalability of the proposed wavelength-domain Ising machine through experiments. To assess its potential for scaling to larger numbers of spins, we tested an all-to-all coupling-type Ising machine. The experimental setup for the wavelength-domain Ising machine, implemented using spatial light modulation, is shown in Fig. 4(a). The coupling coefficient is encoded through intensity modulation with an external WSS. This step is consistent with the on-chip scheme. The state of spin is encoded through a second-level WSS, with the optical signal output from two different ports depending on the spin states. The routing implemented here by the WSS is consistent with the function of the MRR array, and the two outputs match the labels in Fig. 1(a). Following differential power detection and normalization, the results match the expression given in Eq. (6). As a result, the solution can be derived using the same principle as the on-chip scheme. This architecture can realize an Ising machine with a positive coupling coefficient, suitable for solving NPP problems. The input optical spectrum, measured before encoding, is shown in Fig. 4(b). The operating spectral wavelength is in the C-band, and the single-wavelength optical signal has a 3-dB bandwidth of 0.5 nm, limited by both the light source and the WSS channel width used in the experiment. We solved an NPP problem with 256 variables, with the results presented in Fig. 4(c). In Fig. 4(c), different colors represent elements at corresponding positions within different subsets. This experimental result demonstrates the expanded implementation range of the wavelength-domain Ising machine, highlighting its ability to utilize the rich wavelength resources of light and offering excellent scalability.

Based on the architectural design of OIM, the computation time per iteration scales linearly with the number of spins. In contrast, for a GPU-based implementation, the number of operations per iteration increases quadratically with system size. This fundamental difference endows the OIM with a clear advantage in solving large-scale Ising problems, offering greater potential for high-speed computation as problem size grows. Therefore, designing OIMs with strong scalability represents a promising direction for future research. The scalability of the OIM is fundamentally limited by the working domain available for spin encoding. We classify the working domains of Ising machines into three categories: spatial, temporal, and wavelength domains. The wavelength-domain OIM we propose achieves 32 spins on-chip, limited by factors such as the number of MRRs, free spectral range (FSR), and 3-dB bandwidth. By employing a nanobeam-based scheme, which offers a narrow linewidth and is free of FSR, the wavelength-domain OIM chip could enable large-scale on-chip computation of Ising problems. State-of-the-art nanobeam devices can achieve linewidths as narrow as 0.16 nm,63 and through tailored design of parameters such as aperture size, the operational wavelength range can be further extended. For instance, using a 100-nm computing bandwidth, it would be possible to implement up to 312 spins entirely on-chip. Table 1 provides an overview of several recently proposed Ising machines and their key performance metrics, indicating the distinct advantages of different approaches.

The scheme we proposed has three distinct advantages compared with other OIMs.

We propose and experimentally demonstrate a chip-scale, reconfigurable wavelength-domain OIM. This approach leverages the wavelength resources of light to achieve large-scale scalability for the Ising machine. By controlling an array of MRRs and encoding the intensity of input optical signals at different wavelengths, the system can be optimized to reach the lowest energy state using the SA algorithm. We experimentally validate the capability of a 32-MRR on-chip array to implement the Ising machine, employing two encoding methods for different coupling types. This setup successfully solves NP problems such as the NPP and MAX-CUT. Furthermore, by modulating the spatial light intensity, we demonstrate a wavelength-domain OIM with 256 spins, further confirming the scalability of the approach. This method reduces the complexity of solving practical NP problems, thereby lowering the computational power requirements for large-scale systems. Compared with previously proposed Ising machines, our method offers significant advantages, including enhanced flexibility in coupling types and improved scalability with low complexity.

Acknowledgment. This work was supported by the National Key Research and Development Program of China (Grant No. 2022YFB2804203), the National Natural Science Foundation of China (Grant No. U21A20511), and the Knowledge Innovation Program of Wuhan - Basic Research (Grant No. 2023010201010049).

Xinyu Liu received her undergraduate training in optics engineering from the Huazhong University of Science and Technology (HUST), China, and received a PhD in 2025. She is mainly engaged in the research of optical matrix computation and its applications, including optical neural computation and Ising machine.

Wenkai Zhang received his bachelor’s degree from the HUST in 2021. Then he joined the Wuhan National Laboratory for Optoelectronics at the HUST as a PhD candidate. His research interests include photonic integrated circuits and neuromorphic photonics.

Wenguang Xu received his bachelor’s degree from the HUST in 2023. Then he joined the School of Optical and Electronic Information for Optoelectronics at the HUST as a PhD candidate. His research interests include photonic integrated circuits and neuromorphic photonics.

Hailong Zhou received his undergraduate training in optics engineering from the HUST and received a PhD in 2017 and now works as an associate professor at the Wuhan National Laboratory for Optoelectronics of HUST. He is mainly engaged in the research of optical matrix computation and its applications, including optical neural computation and optical logic computation. In recent years, he has published more than 30 SCI academic papers.

Ming Li received his PhD from the University of Shizuoka, Japan, in 2009. In 2013, he was with the Institute of Semiconductors, Chinese Academy of Sciences, as a full professor. He has authored or co-authored more than 210 high-impact journal papers in Nature Photonics, Nature Communications, Light: Science & Applications, and Physical Review Letters. His research interests include integrated microwave photonics and ultrafast optical signal processing.

Jianji Dong received his PhD in optical engineering from the HUST. He is currently a professor with the Wuhan National Laboratory for Optoelectronics. He has authored or co-authored more than 100 Journal papers, including Nature Communications, Light: Science & Applications, and Physical Review Letters. His research interests include integrated microwave photonics and photonic computing. He was honored the Fund of Excellent Youth Scholar by NSFC and the First Award of Natural Science of Hubei Province.

Xinliang Zhang: Biography is not available.