To break the limitations of conventional computing models, several new paradigms in analog computing have emerged recently [1–4]. These draw upon the fundamental physical principles or the incorporation of quantum phenomena, including the minimization of entropy, energy, quantum entanglement, and quantum annealing [5,6]. For instance, the Ising machine (IM) is investigated as a promising tool for tackling complex optimization challenges that are presently beyond the scope of traditional methods [7]. At the core of this innovative computing paradigm lies the integration of physics-based algorithms and hardware, aimed at transforming intricate optimization problems into universal spin Hamiltonians: $$H=-\sum _{ij}^N{J}_{ij}{s}_i{s}_j-\sum _i^N{h}_i{s}_i,$$(1)where $s_i=\pm 1$ represents the spin value, $J_{ij}$ denotes the interaction strength between the $i$-th and $j$-th spins, and $h_i$ represents the external magnetic field [8]. When the spin variables are not restricted to a discrete distribution but can take on continuous values in modulo $2\pi$, this continuous Ising model is referred to as an XY model [9]. The Hamiltonian of the XY model is thus written as $$H=-J\sum _{i,j}^N\cos ({\theta }_i-{\theta }_j)-h\sum _i^N\cos {\theta }_i,$$(2)where the XY spin variable ${\theta }_i$ can be identified with a point on the unit circle. Obviously, the XY model extends spin dimensionality; however, this complex energy landscape also introduces more local minima.

Among the proposed Ising machines, optical Ising machines, typically the coherent Ising machine (CIM) utilizing degenerate optical parametric oscillators [10–12] and the spatial photonic Ising machine (SPIM) based on spatial light modulators (SLMs) [13–17], have shown a significant performance advantage and will become promising hardware architecture for solving large-scale combinatorial problems. Though CIM and SPIM both demonstrate fully programmable over-10,000-spin scale Ising machines, fiber- and spatial-optics-based architectures still remain bulky, are very power demanding, and have poor stability. On-chip photonic Ising machines exhibit a diverse architecture, which can be categorized into two main groups by their utilization of either time coherence or space coherence of light. The first category encompasses integrated CIMs that are developed using a poor man’s structure design [18–21]. The other category comprises unitary matrix networks constructed using a cascaded Mach-Zehnder interferometer (MZI) mesh to enable matrix multiplication and addition operations [22,23]. The former architecture involves the allocation of numerous computational tasks onto FPGAs, and thus poses a challenge to the storage capacity of FPGAs. The latter, accompanied by intricate optical network design, lacks scalability. More importantly, both of these architectural approaches face difficulties in constructing fully connected models and are more commonly employed for sparse Ising models.

Here, to enhance integration and decrease the footprint, we introduce a novel concept, an integrated spatial photonic XY Ising sampler (ISPIS), utilizing a highly uniform multimode interferometer (MMI) implementation. This approach facilitates straightforward integration using off-the-shelf on-chip photonic components and retains the natural full connection of SPIMs. The designed MMI coupler with unequal-width tapers provides multi-optical-channel coherent computing with a low insertion loss (IL) of 0.85 dB and high uniformity (UF) of 0.42 dB. To illustrate the programmable capabilities of general spin coupling interactions, we present the process of ground state search for the XY model through continuous phase sampling of the spin state within a $2\pi$ range, utilizing our ISPIS. Based on this, we further investigate $\pm J$ Ising models and observe the paramagnetic, ferromagnetic, ferrimagnetic, and spin glass phases by simulating equilibrium systems at various temperatures. Additionally, more complex Ising models such as those with external magnetic fields, as well as weighted full rank models are also illustrated. Furthermore, instances of addressing the weighted MaxCut problems provide potential applications in solving combinatorial optimization problems, although the output converges to a stable distribution rather than a state, indicating that further statistics of the sampling results are necessary to obtain the optimal solution.

As shown in Fig. 1, the proposed ISPIS consists of four modules: the coherent input module, the spin encoding module, the coherent superposition module, and the detection and control module. The coherent input module includes a laser and a beam splitter designed to generate N coherent beams with uniform amplitude and phase. The spin encoding module constitutes the central component of the integrated Ising sampling system, comprising a phase shifter array. The ISPIS utilizes a phased array approach akin to SPIM but on the chip level, resulting in specific signal modulations and substantially enhanced performance. Its primary functions include modulating the phase within the Ising model, performing phase calibration, and loading the interaction matrix $J$. The modulated light field can be expressed as $E(\mathit{k})=E_0{\mathrm{e}}^{-\mathrm{i}({\theta }_k)}$. The coherent superposition module is a $1\times N$ MMI differing from SPIM in that the modulated beam is focused using a lens. This configuration offers an evident advantage in that it yields a singular output optical signal on the chip, thereby streamlining the process of detection. In contrast, space-light systems require the center locking of the output spot on the focal plane. The detection and control module consists of an optical power meter and the control circuit. The output light intensity measured by the detector can be expressed as $$I={|E_{\text{total}}|}^2=E_0^2\sum _{k,l}^N\cos ({\theta }_k-{\theta }_l),$$(3)where $E_{\text{total}}={\sum }_{k=1}^N{E}_0{\mathrm{e}}^{-\mathrm{i}({\theta }_k)}$ represents the coherent superposition of all N modulated light fields. Here, Eq. (3) corresponds to an XY model with the external field $h=0$ as $$H=-J\sum _{k,l}^N\cos ({\theta }_k-{\theta }_l),$$(4)where the variable at spin $k$ is represented by a real value ${\theta }_k$ within a range of $[\mathrm{0,}2\pi ]$, and the interaction energy is given by the expression $-J\cos ({\theta }_k-{\theta }_l)$. When the spin values are restricted to any pair of values with a phase difference of $\pi$, such as $\{0,\pi \}$, the formulation of Eq. (4) becomes the classical Ising model as Eq. (1). In addition, the wavefront phase of light theoretically corresponds to the spin value in the XY model, offering a more direct mapping compared to the spin values $s_i=\pm 1$ used in the classical Ising model. And the introduction of the XY model serves another purpose of achieving amplitude modulation, which will be elaborated later. Finally, an electronic processor is utilized to determine whether the current spin state is accepted or not based on the output light intensity. The above process, starting from the generation of a spin state to the acceptance of that state, is considered as a sampling. And iterative sampling is essential to collect enough statistics aligned with specific objectives.

Furthermore, the aforementioned ISPIS structure can also satisfy the mapping to weighted Ising models and the spin glass model. A simple operation is to achieve weight encoding by overlaying the amplitude-modulated phase ${\chi }_k$ on each optical pathway, as shown in Fig. 1(b). The corresponding Hamiltonian can be obtained as $$H=-J\sum _{k,l}\cos (({\theta }_k-{\chi }_k)-({\theta }_l-{\chi }_l)).$$(5)

Then we introduce a phase difference parameter ${\chi }_{k,l}={\chi }_k-{\chi }_l$ into Eq. (5); the Hamiltonian can be expressed as $$H=-J\sum _{k,l}(\cos ({\theta }_k-{\theta }_l)\cos {\chi }_{k,l}+\sin ({\theta }_k-{\theta }_l)\sin {\chi }_{k,l}).$$(6)

By observing Eq. (6), it becomes evident that $\sin ({\theta }_k-{\theta }_l)·\sin {\chi }_{k,l}$ is consistently zero when the values of ${\theta }_k$ and ${\theta }_l$ differ by an integer multiple of $\pi$. Now, by setting the interaction coefficient $J_{k,l}=\cos {\chi }_{k,l}$, we can rewrite Eq. (6) as $$H=-\sum _{k,l}^N{J}_{k,l}\cos ({\theta }_k-{\theta }_l),$$(7)which describes the spin glass model and can be used to solve various types of weighted combinatorial optimization problems.

As a proof of concept, an eight-spin ISPIS is fabricated on the silicon-on-insulator (SOI) platform, with the micro-graphs displayed in Fig. 1(c). A laser operating at a wavelength of 1.55 μm is coupled into the waveguide through an input grating coupler and evenly distributed into eight channels with a three-stage cascaded 50:50 splitter. The spin states are encoded onto the phase term of the light field through eight thermal phase shifters using 100-μm-long, 2-μm-wide TiN metal strips. To minimize thermal cross talk, the design specifies a separation of more than 100 μm between adjacent phase shifters. The modulated light is subsequently directed into a highly uniform $1\times 8$ MMI, facilitating coherent superposition among the channels. The output light is detected by optical power meters at the output interface and is then transmitted to the central processing unit (CPU). Based on the detection intensity, the CPU generates the spin state for the subsequent iteration and modulates the phase shifters through a multichannel programmable voltage source. Following processing and device-level testing, the chip is packaged to facilitate subsequent experiments. The packaging process involves the attachment of high-density electrodes to a printed circuit board (PCB) through wire-bonding, as illustrated in Fig. 1(d).

While the realization of optical discrete Fourier transform (DFT) using $N\times N$ MMI couplers has been proposed [24,25], the requirement for DFT in optical Ising machines differs somewhat. In a SPIM system, the lens takes on the DFT, but in fact only the light intensity at the focal position needs to be detected. Essentially, the Fourier device in this context primarily facilitates coherent superposition, a function that can be effectively substituted by a $1\times N$ MMI. And in ISPIS systems, the uniformity of the MMI coupler is the primary metric to be considered. In addition, a low insertion loss with phase error is also desired.

To validate the design, we fabricate and characterize the $1\times 8$ MMI splitter on an SOI platform, employing 10 chips for testing purposes. We measure the output light intensity of the eight ports of the optimized MMI at 1550 nm. After normalization, we compared them with the simulation results, as shown in Fig. 2(c). Additionally, we conduct tests on the cascaded $1\times 2$ MMI structure and the $1\times 8$ MMI coupler with equal-width tapers as a control. The test IL and UF within a wavelength range of 1530–1565 nm are presented in Figs. 2(d) and 2(e). The optimized $1\times 8$ MMI coupler with unequal-width tapers exhibits a UF of 0.42 dB and an IL of 0.85 dB at 1550 nm, which are better than the initial results of 1.10 dB and 2.05 dB. A UF improvement of 0.39 dB is also pronounced when compared to the cascaded $1\times 2$ MMI structure, highlighting the superior uniformity achieved by the former. Therefore, the optimized MMI structure provides superior performance in terms of uniformity and insertion loss across the C-band.

Due to the inter-channel random errors of the system and the intrinsic errors associated with the $1\times 8$ MMI, the Ising system exists in a stochastic initial state, where each spin exhibits a continuously adjustable phase ($\theta$) instead of a discrete phase ($\pm 1$), corresponding to the solution space of the XY model. Unlike the classical binary Ising machine, an XY Ising machine necessitates high accuracy to realize more phase states. In our ISPIS, the resolvable phase control is contingent upon the precision of the thermal phase shifters. The device tests demonstrate that a phase control precision of 6 bits is achieved. And the recorded intensities will be curve-fitted to establish the relationship between electrical power and optical phase, as illustrated in Fig. 3(a).

To search the ground spin state of the XY Ising sampler, the phase of one channel is selected as a reference, and the control voltages of the shifters in the remaining channels are adjusted sequentially to achieve a phase change to $2\pi$. A specific operation involves incrementally raising the control voltage of a designated channel with a step of 50 mV and monitoring the output light intensity. After the voltage sweeping range typically spans 1.5 to 1.7 phase cycles, the voltage of the current channel is adjusted to the value corresponding to the maximum output optical power and the aforementioned process is repeated for the next channel’s phase shifter. This method is applied sequentially to complete the phase calibration for the remaining seven channels, and ultimately fine-tuning the phase of the reference channel concludes the ground state search process. Figure 3(b) illustrates the ground state search for aligning the phase, and the IL of the system is measured to be 2.12 dB. In addition, Fig. 3(c) demonstrates the phase changes of each channel to facilitate phase encoding in subsequent experiments.

To evaluate the performance of the ISPIS, we experimentally traversed all 256 spin states of the eight-spin fully connected Ising model ($J=1$). Additionally, we construct the on-chip photonic domain utilizing the appropriate photonic devices simulated by Lumerical Interconnect software. The experimental and simulated output light intensities are as illustrated in Fig. 4(a). According to Hamiltonian quantities, all spin states can eventually be classified into five distinct merger states ($H=0,-4,-16,-32,-64$). The light intensity values detected by the optical power meter in the experiment are represented by the blue box plots, while the orange line plots depict the error-free simulation results. According to the trend of the lines in the figure, it can be noted that the output light intensity in both the experiment and simulation exhibits a logarithmic association with the theoretical Hamiltonian. While the simulation results demonstrate a precise alignment with the theoretical Hamiltonian, the experimental outcomes display a certain range of distribution [see the error bar in Fig. 4(a)]. This is caused by the weak intensity nonuniformity between the eight waveguides.

In subsequent experiments, we simulate several well-studied spin systems at various temperatures using the Monte Carlo algorithm. First, we perform 20 sets of full samples encompassing all the spin states of the different Ising models and compile a look-up table comprising all the spin states along with their corresponding detected light intensity values. Next, the Monte Carlo algorithm is employed to undertake a stochastic traversal of the spin states, utilizing a randomly selected loss function (the Hamiltonian) from the 20 sets of light intensity values corresponding to the current spin state. The length of the Markov chain will affect the result of sampling. In particular, fluctuations near the critical temperature are significant, necessitating extended equilibration times and an increased number of statistical samples to achieve a more accurate average value. Therefore, we take Markov chains with different chain lengths for sampling, specifically 100 and 1000 iterations. The experimental results for these varying chain lengths are presented in Figs. 4(b)–4(e). The parameter T indicates the relative magnitude of temperature rather than representing an absolute physical temperature value. Despite the variations in chain length, the phases exhibit a consistent trend with temperature. As the length of the Markov chain increases, the relationship between the Hamiltonian and phase transition with temperature becomes more pronounced, thereby facilitating the determination of the critical temperature of the phase transition phenomenon.

Here, we investigate the process of phase transition with different spin interactions and observe the three kinds of Ising models: fully connected $\pm J$ model, Ising model under external magnetic field, and weighted Ising model with full rank. For the $\pm J$ model, the interactions $J_{i,j}$ between spins are either 1 with a probability of $p$ or $-1$ with a probability of $1-p$. And the explicit form of the probability distribution is $p(J_{i,j})=p·\delta (J_{i,j}-J)+(1-p)·\delta (J_{i,j}+J)$. Depending on the value of probability $p$ and temperature, the spin will show four different distributions: paramagnetic phase, ferromagnetic phase, ferrimagnetic phase, and spin glass phase. As the temperature increases, the thermal motion of spins gradually surpasses the inter interactions, resulting in the substance exhibiting macroscopic paramagnetic properties. As shown in Figs. 4(b)–4(d), it is evident that at high temperatures, the magnetization of all three models is rendered zero, indicating a disordered spin state. However, as the temperature decreases, the Hamiltonian of the Ising system decreases gradually, and the system manifests ordering characteristics where the spin state is mainly determined by the inter interaction. The probability $p$ governs whether the spins primarily exhibit promoting or excluding interactions among each other. When $p>0.5$, the spin-spin interactions show strong interaction, resulting in the entire system exhibiting ferromagnetic ($p=1$) or ferrimagnetic phase ($0.5<p<1$), as shown in Figs. 4(b) and 1(c). When $p<0.5$, strong excluding interactions result in the spin glass phase, as shown in Fig. 4(d). Moreover, net overall magnetization is observed in the ferromagnetic and ferrimagnetic system. The bifurcation of the curves can be attributed to the fact that, at low temperatures, the spin clusters may remain in either 1 or $-1$ state, with equal probabilities in the absence of external field intervention. Additionally, the ferromagnetic system displays a greater magnetization susceptibility and the spins tend to line up in the same direction [see Fig. 4(b)]. Ferrimagnetic phase is somewhat like the antiferromagnetic phase in that the spin clusters align antiparallel [see Fig. 4(c)]. However, some spin clusters are larger than others, and this imbalance in cluster sizes results in a net overall magnetization. In a spin glass system, it is impossible to minimize the energy of ferromagnetic interaction and antiferromagnetic interaction at the same time. The spins are randomly oriented at low or high temperatures, resulting in an overall magnetization of zero [see Fig. 4(d)].

We also investigate the ferromagnetic phase with an external magnetic field by simulating the equilibrium systems at different temperatures. In our experimental setup, we mimic a uniform external magnetic field $h$ by fixing the spin states in specific channels. A six-spin Ising model with magnetic fields $h=4$ is presented in Fig. 4(e). The magnetization deviates from splitting into two clusters with equal probability and instead tends to align with the direction of the external field.

To demonstrate the generality of the ISPIS system, we employ a time-division multiplexing method to achieve full sampling of a weighted full-rank Ising model. The eigenvalue decomposition method is performed on the full-rank matrix, yielding multiple low-rank matrix slices [15,26]. These processed low-rank matrices are subsequently loaded into the ISPIS for full sampling. The coefficients of each slice matrix are modulated by amplitude-modulated phase $\chi$. Finally, the Hamiltonian of the original full-rank model is calculated in the electric domain by summing the output light intensities of all Ising model slices for every sample. The total output light intensity calculated in the CPU is shown as the blue solid scatter in Fig. 4(f), while the simulation results are indicated by the orange solid scatter. It can be observed that the Interconnect simulation results exhibit perfect linear correlation with the theoretical calculations, validating that our design effectively facilitates mapping and solving of the full-rank Ising models. Though the experimental results are susceptible to noise, leading to some errors, they still maintain a robust linear correlation with the theoretical results. As a result of the linear fitting of the experimental data, the coefficient of determination $R^2=0.82$.

To further highlight the capability of this work in addressing combinatorial optimization problems, searching for the optimal solution of eight-node unweighted and weighted MaxCut instances is demonstrated. The simulated annealing (SA) algorithm is conducted based on the sampled results; 100 repeated searches are implemented based on ISPIS for the unweighted and weighted instances. The evolution curves of the experimentally detected optical power are shown in Figs. 5(a) and 5(b). The evolution curves of the corresponding cut values are presented in Figs. 5(c) and 5(d). The slight difference between the detected intensity and calculated cut value mainly results from the amplitude nonuniformity of the optical channels and the small thermal drift of the whole system during operation. Figures 5(e) and 5(f) show the probability distributions of the search results. The identified solutions populate the tail of the cut value distribution, and ISPIS gives the correct maximum solutions with 69% and 91% probability, respectively.

We have proposed and implemented an ISPIS utilizing on-chip phase modulators and an optical discrete Fourier transform (DFT) device. Compared to SPIMs, our design offers the advantage of a more compact and controllable setup. The highly integrated structure avoids the coupling of spatial light, thereby enhancing the system’s stability and precision. Moreover, unlike traditional SPIM systems that are constrained by the refresh frequency of spatial light modulation (180 Hz), the ISPIS holds promising potential for achieving high-speed ground state search. Although the current design uses thermal phase modulation with a response rate of 16 μs, the main contribution of this work lies in showcasing the on-chip SPIM architecture and its capability to construct various Ising models. Recent advancements in electro-optic integrated modulators enable high modulation rates on both conventional SOI platforms and thin-film lithium niobate platforms. Replacing the thermal phase shifters with an electro-optic version can boost operation speed and reduce energy consumption.

In addition to modulation speed, the performance of the ISPIS may be influenced by factors such as signal-to-noise ratio (SNR) and phase noise. SNR is mostly determined by IL and photodetection noise. With the utilization of advanced waveguide and photodetection technologies, high SNR is expected to enable high-quality sampling.

Unlike other integrated photonic computing architectures, photonic Ising machines are particularly sensitive to phase noise, resulting in large variances in the distribution of experimental sampling results. To mitigate the impact of phase errors, we adopt various means, including the design of highly homogeneous MMI couplers, phase calibration, and sampling schemes instead of solver schemes. In future work, the incorporation of an intensity modulation module could be implemented on the chip to facilitate independent and precise intensity modulation and calibration, thereby alleviating the assignments of the phase-shifter array.

Furthermore, the ISPIS has the potential for spatial scaling to leverage the optical advantages of parallelism. Generally, to construct an N-channel integrated photonic Ising machine, the number of modulators required according to our design operates at the $O(N)$ level, whereas the architecture of the MZI network requires a complexity of $O(N^2)$ [22,27]. Consequently, our design is more favorable for achieving large-scale integration, as it requires only an expansion of the modulator array without a redesign of the network structure. Additionally, $1\times N$ MMI couplers can be cascaded to enable coherent computation across a larger number of optical channels. Recent advancements in space-division multiplexing and time-division multiplexing present promising solutions to these scalability limitations. For example, TSMC develop two sets of 1-to-512 fan-out circuits based on multi-layer silicon nitride, achieving compact, large-scale architectures with a total optical path loss of around $-35\mathrm{dB}$ [28]. The Taichi photonic chiplets incorporate 64 modulators for weight loading and 56 phase shifters for unitary matrix transformation, successfully demonstrating on-chip classification at the level of 1000 categories [29]. The 16-channel photonic solver for optimization problems integrates 16 high-speed electro-optic modulators, 88 thermo-optic phase shifters, and 16 balanced photodetectors on the SOI substrate [30]. The integrated optoelectronic IM successfully addresses the MaxCut problems involving up to 16,384 spins using just one electro-optical (EO) modulator [19].

In summary, we propose the ISPIS based on a highly uniform MMI design to facilitate the simulation of various spin phases and the solution of combinatorial optimization problems. We simulate the classical Ising model using the expression of the XY model, which primarily provides a method for phase-amplitude modulation, thereby reducing the number of on-chip modulators, as well as the complexity of control and the overall chip footprint. Furthermore, an unequal-width taper-based MMI coupler with the IL of 0.42 dB and UF of 0.85 dB enhances the uniformity of the system and ensures that the contributions from each spin Ising system are equivalent. The uniformity of the system is influenced by the phase of the optical waveguide, leading to the formulation of an XY model for optimization. The initial phase is calibrated by searching the ground state of the XY model, which facilitates (1) the assessment of uniformity in response to varying spin configurations and (2) the identification of an appropriate operating point and driving voltage. The performance of ISPIS is demonstrated by realizing fully programmable spin couplings and external magnetic fields. The linear correlation between the experimental results and theoretical calculations reached 0.82. Additionally, the (un)weighted MaxCut instances solved illustrate the capability of ISPISs to function as general solvers for optimization problems. Furthermore, thanks to recent advances in hybrid integration of electronics and photonics, it is optimistic to build a monolithic system for high-performance computing with an ISPIS and an affiliated driver board for signal control. Consequently, our design enables the miniaturization of programmable SPIMs into a single energy-efficient photonic integrated circuit, hence taking full advantage of the high bandwidth of optical systems and significantly enhancing speed compared to existing SPIM concepts.