Wavelength-division multiplexing optical Ising simulator enabling fully programmable spin couplings and external magnetic fields

Ising model is an archetypal model widely used for investigations of complex dynamics in physics, computer science, biology, and even social systems. Owing to Moore’s law for conventional computers, there has been tremendous interest and a boost in the development of unconventional computing architectures for simulating Ising Hamiltonians, for example, based on optical parametric oscillators [1, 2, 3, 4, 5], lasers [6, 7, 8, 9, 10], polariton [11, 12, 13], trapped ions [14], atomic and photonic condensates [15, 16], electronic memorisers [17], superconducting qubits [18, 19, 20], and nanophotonics circuits [21, 22, 23, 24, 25, 22]. Despite different approaches and technologies, it is worth noting that the error probability and time-to-solution metrics of these Ising machines have similar scaling trend as a function of the number of spins [26]. Practically, the difficulty to implement the spin coupling interactions with the proposed hardware has become the main factor limiting scalability and performance for unconventional Ising simulators [27, 28]. Also complete characterization of possible stable phases is necessary for estimating the Ising description and plays a key role to understand the working principle in spin systems [29, 30, 31, 32, 33].

In this regard, by encoding spins on the phase terms of a monochromatic field with spatial light modulators (SLMs), spatial photonic Ising machines (SPIMs) benefit with reliable large-scale Ising spin systems, even up to thousands of spins by exploring the spatial degrees of freedom [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. Like other optical analog computations [45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55], the calculation of spin system energy is just by instantaneously measuring the light intensity, therefore with ultrafast speed and high power efficiency [38]. Moreover, to eliminate the impact of pixel alignment, the gauge transformation is proposed to simultaneously encode spin configurations and interaction strengths with a single spatial phase modulator [40]. However, the original proposed SPIM [38] is only applicable to Mattis-type coupling interactions [56]. Even with scattering medium, tunable SPIM was demonstrated based on multiple light scattering, while the Ising spin system is still limited to specific fully connected couplings [43]. Therefore, how to realize completely programmable spin couplings is a primary target to develop SPIMs for general Ising models.

Here we report a wavelength-division multiplexing SPIM to enable fully programmable spin couplings and external magnetic fields as well. Beyond Mattis-type interaction, we propose a gauge transformation for general Ising models with arbitrary spin interactions. More importantly, with wavelength-division multiplexing, we optically compute the general spin Hamiltonian with the advantages of large scale and ultrafast speed. We experimentally demonstrate such a wavelength-division multiplexing SPIM with a single SLM, where the gauge transformation is implemented to eliminate the impact of pixel alignment. To show the programmable capability of general spin coupling interactions, we explore three kinds of Ising models: fully connected $\pm J$ model, Sherrington-Kirkpatrick (SK) model, and only locally connected ${{J}_{1}}\texttt{-}{{J}_{2}}$ model. For $\pm J$ and SK models, we investigate the process of phase transition with different spin interactions and experimentally observe the spin-glass (SG), ferromagnetic (FM), and paramagnetic (PM) phases by simulating the equilibrium systems at different temperatures. We verify that the critical temperatures evaluated by the wavelength-division multiplexing SPIM are consistent with the predictions of the mean-field theory. We also demonstrate the phase transition from PM to stripe-antiferromagnetic phase in the ${{J}_{1}}\texttt{-}{{J}_{2}}$ model with competing interactions and experimentally observe the presence of SG phase by increasing the stochasticity of the next-nearest-neighbor interactions. These results show the great programmable flexibility of spin couplings and external magnetic fields by the wavelength-division multiplexing SPIM, which provides important potential applications in solving combinatorial optimization problems.

The wavelength-division multiplexing SPIM, designed for full programmability of spin interactions and magnetic fields, is shown in Fig. 1(a). The illumination components consist of a collimated super-continuum laser, a diffraction grating, and a cylindrical lens. The diffraction grating and cylindrical lens diffract light with different wavelengths onto a SLM along the $x$-axis, while the $y$-axis pixels are illuminated coherently by the same wavelength. By adjusting the diffraction angle, the operation wavelengths are selected within the quasi-stationary region of super-continuum light [57]. The optical field modulated by the SLM is then transformed by a Fourier lens, resulting in an incoherent intensity summation of different wavelengths and coherent interference for each wavelength at the back focus plane.

Next we show that the wavelength-division multiplexing SPIM can optically compute the Hamiltonian of any Ising model by detecting the intensity at the center position of the back focus plane. The Ising model is given by $H=-\sum\limits_{i\neq j}{{J_{ij}}{\sigma_{i}}{\sigma_{j}}}-\sum\limits_{i}^{N}{{h_{i}}}{\sigma_{i}}$, where ${J_{ij}}$ and ${{h}_{i}}$ are the interaction strength and the external magnetic field strength for $N$ spins, respectively, and ${{\sigma}_{i}}$ can take binary values $+1$ or $-1$. By adding an auxiliary spin with ${{\sigma}_{N+1}}=1$, we use a Cholesky-like decomposition to transform the Hamiltonian into multiple Mattis models ${{H}_{k}}=-\sum_{i,j=k}^{N+1}{J\kappa_{i}^{k}\kappa_{j}^{k}{{\sigma}_{i}}{{\sigma}_{j}}}$, so that $H=\sum_{k=1}^{N}{{{H}_{k}}}+{{H}_{0}}$. The interaction strength for the $k$-th Mattis model is given by $J_{ij}^{k}=J\kappa_{i}^{k}\kappa_{j}^{k}$, where $\left|\kappa_{i}^{k}\right|\leq 1$, $J$ and ${{H}_{0}}$ are constants with the unit of energy. The transformation is summarized in Sec. I of the Supplemental Material (SM), and we note that despite the transformation, the degrees of freedom of $\kappa_{i}^{k}$ remain equal to $N(N+1)/2$, which is the summation of those of $\{J_{ij}\}$ and $\{h_{i}\}$ for arbitrary interactions and magnetic fields.

Inspired by the encoding scheme for Mattis model [40, 41], here we propose a general one for arbitrary Ising models by a gauge transformation. This transformation allows to encode the spin configurations and program the interaction strengths on a single phase-only SLM, as shown in Fig. 1(b). The transformation rotates each original spin about the $z$-axis by an angle ${\alpha}_{j}^{k}=\arccos\kappa_{j}^{k}$ to arrive at a new spin vector, which is then projected onto the $z$-axis to obtain the effective spin $\sigma_{j}^{\prime k}=\kappa_{j}^{k}\sigma_{j}$. By the gauge transformation given as

the transformed Hamiltonian remains invariant, $H=-J\sum_{k=1}^{N}{\sum_{i,j=k}^{N+1}{\sigma_{i}^{{}^{\prime}k}\sigma_{j}^{{}^{\prime}k}}}+{H_{0}}$, where the interactions of the transformed spins are uniform with a strength of $J$ in both short and long ranges.

We encode the gauge-transformed spin configurations $\sigma_{j}^{{}^{\prime}k}$ on a single phase-only SLM with wavelength-division multiplexing. For the $k$-th Mattis model, we assume a uniform spectrum intensity light illuminating on SLM and apply a phase modulation on the $y$-directional pixels for each spin as

We note that for each Mattis model, the phase modulation should account for the calibration of different wavelengths when encoding $\varphi_{j}^{k}$ in the $x$ direction on the SLM. The normalized intensity $\tilde{I}$ at the center position of the back focus plane is the summation of the incoherent field intensities for light with different wavelengths, $\tilde{I}=\sum_{k=1}^{N}{\sum_{i,j=k}^{N+1}{\sigma_{i}^{{}^{\prime}k}\sigma_{j}^{{}^{\prime}k}}}$. The Ising Hamiltonian for the general spin interactions is thus optically computed as $H=-J\tilde{I}+{{H}_{0}}$. The details of the gauge transformation and the encoding for general illumination cases are described in SM Sec. II.

To evaluate the performance of the wavelength-division multiplexing SPIM, we conduct an experiment (detailed experiment setup can be found in SM Sec.III) and simulate three well-studied spin systems: the $\pm J$ model [Fig. 2(a)], the SK model [Fig. 2(d)] and the SK model under external magnetic field [Fig. 2(g)]. The interactions between spins in the $\pm J$ model are either $1$ with a probability of $p$ or $-1$ with a probability of $1-p$. The probability distribution is given as $P\left({{J}_{ij}}\right)=p\delta\left({{J}_{ij}}-J\right)+(1-p)\delta\left({{J}_{ij}}+J\right)$, where $\delta(x)$ is a Kronecker delta function equal to 1 for $x=0$ and 0 for other cases. In the SK model, the probability distribution of ${{J}_{ij}}$ is Gaussian, with the distribution given as $P\left({{J}_{ij}}\right)=\mathcal{N}({{J}_{ij}};\frac{{{J}_{0}}}{N},\frac{{{{\Delta J}}^{2}}}{N})$, where ${J}_{0}$ and ${\Delta J}$ are two constant parameters, such that the energy can be extensive and proportional to $N$ [32].

We use the wavelength-division multiplexing SPIM to examine the phase transition with 80 spins. For each model, a quenched realization is generated for the interactions ${{J}_{ij}}$ randomly assigned based on their respective probability distributions. By varying the temperature and utilizing the Markov chain Monte Carlo algorithm, we generate 100 replicas from the final configuration, where each replica is obtained by randomly initializing spin configurations and through 800 iterations of the optical Metropolis Hasting sampling procedure [40, 58, 39]. The spin overlap is then calculated as ${{q}_{\alpha\beta}}=\frac{1}{N}\sum_{i=1}^{N}{\sigma_{i}^{\alpha}}\sigma_{i}^{\beta}$, which measures the similarity between replicas $\alpha$ and $\beta$. The phase transition is characterized by the probability density function $P({q})$ of the overlap.

Figures 2(b) and (c) present the experimental results for the $\pm J$ model with the parameters $p=0.7$ and $0.55$, respectively. At high temperatures, both of the two parameters $p=0.7$ and $p=0.55$ result in randomly arranged spins and little correlation between replicas, indicating the PM phase, where $P({{q}})$ has a peak at around zero. At low temperatures, however, $P({{q}})$ displays the distinct density distributions for the two values of $p$. For $p=0.7$, since the interactions are mostly positive, the spins attract each other and the replicas remain in only two ground states of the FM phase. Consequently, $P({{q}})$ at low temperature has two peaks around $1$ and $-1$ as shown in Fig. 2(b). On the other hand, for $p=0.55$, the interactions are composed of both positive and negative values, causing frustration during the energy minimization process at low temperatures. This frustration results in a multi-valley energy landscape and $P({{q}})$ takes on a wide range of values at low temperature, as depicted in Fig. 2(c). This feature of widely distributed ${{q}}$ is the hallmark of the SG phase. These results demonstrate that the SG phase transition indeed emerges in the systems simulated by the wavelength-division multiplexing SPIM.

To further study the phase transition with the wavelength-division multiplexing SPIM, we examine the SK model without magnetic fields and compare the estimated critical temperature with the mean-field theory. Fig. 2(e) and (f) show the results for ${{J}_{0}}=40$ and $8$, respectively, with $\Delta J=\sqrt{80}$. At high temperatures, both figures show the PM phase with $P({{q}})$ dominated around zero, due to the randomly arranged spin configurations. However, as shown in Figs.  2(e) and (f), $P({{q}})$ at low temperatures exhibit the FM and SG phases for the different values of ${{J}_{0}}$ and occur around $T_{c}={{J}_{0}}$ and $T_{c}=\Delta J$, respectively. The experiment results are consistent with the mean-field theory [32] with the critical temperatures for the PM-FM transition and for the PM-SG transition.

Additionally, we also investigate the stability of the complex multi-valley energy landscape of the SG phase in the presence of a uniform external magnetic field in the SK model [Fig. 2(g)]. As shown in Fig. 2(h), for the SK model with ${{J}_{0}}=8$ and $\Delta J=\sqrt{80}$, a weak magnetic field with $h=0.2$ aligns spins with the field direction even at high temperature, causing most ${{q}}$ values to cluster around positive values. Upon decreasing the temperature, the probability distribution of ${{q}}$ resembles that of the SG phase covering a wide range, indicating that weak magnetic fields are not strong enough to completely alter the multi-valley energy landscape. However, when subjected to a stronger magnetic field with $h=2$ [Fig. 2(i)], at low temperatures ${{q}}$ are more clustered, which suggests that the magnetic field flattens out some valleys in the energy landscape.

The wavelength-division multiplexing SPIM can be utilized to study Ising systems beyond the full-coupling models. To demonstrate its full programmability of spin couplings, we examine a locally connected ${{J}_{1}}\texttt{-}{{J}_{2}}$ model [Fig.3(a)]. The Hamiltonian only includes nearest-neighbor ferromagnetic and next-nearest-neighbor antiferromagnetic interactions on a square lattice with cyclic boundary conditions: $H=-{{J}_{1}}\sum\limits_{\langle ij\rangle}{{{\sigma}_{i}}}{{\sigma}_{j}}-{{J}_{2}}\sum\limits_{\langle\langle ij\rangle\rangle}{{{\sigma}_{i}}}{{\sigma}_{j}}$, where $\langle\rangle$ and $\langle\langle\rangle\rangle$ denote the nearest and next-nearest neighbors, respectively. The nearest-neighbor ferromagnetic interactions align adjacent spins with $J_{1}>0$ (represented by solid lines in Fig.3(a)), while the next-nearest-neighbor antiferromagnetic interactions drive two adjacent rows and columns to have opposite orientations with $J_{2}<0$ (represented by double parallel lines). When the antiferromagnetic interaction $J_{2}$ is strong enough to overcome the ferromagnetic interaction $J_{1}$, a striped phase is produced, characterized by a two-component order parameter $(m_{x},m_{y})$, where $m_{x}=\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}(-1)^{x_{i}}$, and ${{m}_{y}}=\frac{1}{N}\sum_{i=1}^{N}{{{\sigma}_{i}}}{{(-1)}^{{{y}_{i}}}}$, and $(x_{i},y_{i})$ are the coordinates of spin $\sigma_{i}$. The ground states of the ${{J}_{1}}\texttt{-}{{J}_{2}}$ model are ${{\mathbb{Z}}_{4}}$ ordered when $|J_{2}|/J_{1}>1/2$, and can be represented by the order parameters of $(m_{x},m_{y})=(\pm 1,0)$ and $(m_{x},m_{y})=(0,\pm 1)$, corresponding to two longitudinal and transverse striped states with opposite spin directions, respectively [59].

Using the wavelength-division multiplexing SPIM, we simulate the transition between the striped states and the PM phases. Fig. 3(b) shows the simulation results for the parameters ${{J}_{1}}=0.2$ and ${{J}_{2}}=-1$ on an $8\times 8$ lattice. Above the critical temperature of $T_{c}=20J_{1}$, the order parameters $(m_{x},m_{y})$ are both close to zero due to the thermal fluctuations destroying the long-range correlation between spins as a sample of spin configurations shown in Fig. 3(c). However, below $T_{c}$, the order parameters split into four clusters located in different quadrants, corresponding to the four striped states [Fig. 3(b)]. Figs. 3(d-g) show four samples of the spin configurations, which exhibit a long stripe spatial distribution. These results for the ${{J}_{1}}\texttt{-}{{J}_{2}}$ model demonstrate the programmability of spin couplings.

The ${{J}_{1}}\texttt{-}{{J}_{2}}$ model plays a crucial role in comprehending the low-temperature phase of short-range spin glass. With the help of the wavelength-division multiplexing SPIM, we consider the next-nearest-neighbor interactions are Gaussian as $P\left({{J}_{2}}\right)=\mathcal{N}({{J}_{2}};{{J}_{0}},{\Delta J}^{2})$ with a mean value of ${{J}_{0}}=-1$ and a standard deviation $\Delta J$, while maintaining a fixed value of nearest-neighbor interactions ${{J}_{1}}=1$. For small variances of $\Delta J=0.2$ [Fig. 4(a)], the density distribution of ${{q}}$ exhibits a single peak at high temperatures, signifying the PM phase. However, as the temperature decreases, $P({{q}})$ transforms into three clusters of peaks, indicative of the ${{\mathbb{Z}}_{4}}$ ordered striped ground states. The peak around ${{q}}=0$ is twice as high as those around ${{q}}=\pm 1$. In contrast, by increasing the standard deviation to $\Delta J=2$, the increased disorder in the next-nearest-neighbor interactions results in a SG phase at low temperatures. The SG phase is characterized by many pairs of ground states and results in multiple sharp peaks of $P({{q}})$ as Fig.4(b).

We also experimentally measure the susceptibility $\chi=\frac{1}{N}\sum_{i=1}^{N}{\frac{{1-m_{i}^{2}}}{T}}$[60] to investigate the transition temperature in the ${{J}_{1}}\texttt{-}{{J}_{2}}$ model, where ${{m}_{i}}$ represents the time ensemble average of ${{\sigma}_{i}}$. At high temperatures, thermodynamic fluctuations cause the spin orientation to constantly change, resulting in ${{m}_{i}}$ being close to zero, as shown in Fig. 4(c). The susceptibility in both two cases decreases with increasing temperature and scales with $\frac{1}{T}$. The results also indicate that the critical temperatures are different, with ${{T}_{c}}=14{{J}_{1}}$ and ${{T}_{c}}=10{{J}_{1}}$ for $\Delta J=0.2$ and $\Delta J=2$, respectively. In the case of $\Delta J=0.2$, the interactions of the next-nearest-neighbors are relatively small and close to $-1$. As the temperature decreases below ${{T}_{c}}$, the system reaches the striped ground states with ${{m}_{i}}$ approaching $\pm 1$, causing the average susceptibility $\chi$ to converge to zero. In contrast, for $\Delta J=2$, at low temperatures, the local spin orientation still varies slowly due to the presence of numerous ground states in the SG phase, resulting in $\chi$ being close to a non-zero constant. These experimental results clearly demonstrate the spin-glass phase transition in the ${{J}_{1}}\texttt{-}{{J}_{2}}$ model with short-range interactions.

In summary, we propose the wavelength-division multiplexing SPIM to realize fully programmable spin couplings and external magnetic fields. With the full-coupling $\pm J$ models, SK models and the only locally connected ${{J}_{1}}\texttt{-}{{J}_{2}}$ models, we demonstrate the great programmable flexibility of spin couplings and external magnetic fields with the wavelength-division multiplexing SPIM. By simulating the equilibrium systems at different temperatures, we experimentally observe the phase transitions among the spin-glass, the ferromagnetic, the paramagnetic and the stripe-antiferromagnetic phases with the wavelength-division multiplexing SPIM. The exhibited rich phase diagrams are beneficial to searching for the true ground state of the Ising model, and thus the wavelength-division multiplexing approach provides important potential applications in solving combinatorial optimization problems.

The authors acknowledge funding through the National Key Research and Development Program of China (2022YFA1405200) and the National Natural Science Foundation of China (12174340).

Supplementary Material: Wavelength-division multiplexing optical Ising simulator enabling fully programmable spin couplings and external magnetic fields

Li Luo, Zhiyi Mi, Junyi Huang, and Zhichao Ruan

Interdisciplinary Center of Quantum Information,
State Key Laboratory of Modern Optical Instrumentation,
and Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou 310027, China