Spatial Photonic Ising Machine (SPIM)

• SPIM is an optical analog computing platform where Ising spins are encoded in light’s spatial phases and evaluated through free-space propagation and Fourier measurements.
• It employs gauge transformations and multiplexing techniques to expand from rank-one interactions to fully programmable, high-rank Hamiltonian models.
• SPIM integrates iterative annealing protocols and measurement-feedback loops to efficiently solve large-scale combinatorial and spin-glass optimization problems.

Spatial Photonic Ising Machines (SPIMs) are optical analog computing platforms in which Ising spins are encoded in the spatial degrees of freedom of light, the energy of a candidate spin configuration is evaluated by free-space propagation and Fourier-plane intensity measurement, and a measurement-feedback loop is used to search for low-energy states of an Ising Hamiltonian. In the canonical formulation, the target problem is mapped to

$H=-\sum_{i,j}J_{ij}\sigma_i\sigma_j,$

with binary spins $\sigma_i\in\{-1,1\}$, and the optical system provides the quadratic interaction terms through interference. Across the literature, SPIMs have evolved from large-scale Mattis-type optical Ising solvers into architectures with gauge-transformed single-SLM implementations, multiplexed and decomposed higher-rank Hamiltonians, direct arbitrary-coupling encodings, and hybrid optical-digital optimization schemes for combinatorial optimization, spin-glass simulation, and learning tasks (Pierangeli et al., 2019, Pierangeli et al., 2020, Zheng et al., 22 May 2026).

1. Foundational formulation and physical principle

The original SPIM paradigm encodes each spin as a binary optical phase on a spatial light modulator (SLM), typically with $\phi_j\in\{0,\pi\}$ and $\sigma_j=\exp(i\phi_j)=\pm1$, while optical propagation and camera-based detection transform pairwise phase correlations into a measurable intensity. In the primitive formulation, the focal-plane center intensity gives a Mattis-type quadratic form,

$$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$

so the interaction matrix is effectively rank one. In target-image formulations, minimizing the discrepancy between measured intensity and a prescribed target image is equivalent to minimizing an Ising Hamiltonian with couplings determined by the amplitude profile and the Fourier content of the target image (Pierangeli et al., 2019, Shimomura et al., 26 Feb 2025).

A defining feature of SPIMs is that the Hamiltonian evaluation is optical, whereas the search itself is usually embedded in an iterative measurement-feedback loop. In one widely used protocol, a single spin is randomly selected and flipped, the resulting optical pattern is detected, and the update is accepted or rejected according to the measured cost or a Metropolis rule. This distinguishes SPIMs from architectures that explicitly compute local effective fields electronically: the optical system directly produces the global quadratic observable from which the optimization signal is inferred (Pierangeli et al., 2020).

From the outset, scale was a central motivation. Early experiments reported configurations with thousands of spins, including an initial demonstration at $N=4\times10^4$ spins for the fully connected case and a size range from $N=74\times74$ to $N=274\times274$, establishing SPIM as a large-scale free-space photonic Ising platform rather than a small-spin proof of concept (Pierangeli et al., 2019).

2. Optical architectures, phase encoding, and gauge transformation

The hardware core of a SPIM is simple in concept but varied in implementation: a coherent or incoherent source illuminates an SLM, a lens or 4-f system performs the Fourier transform, and a CCD, CMOS, or qCMOS detector reads out the optical observable. Early architectures often separated amplitude and phase functions across two SLMs, with one modulator encoding the amplitudes $\xi_i$ and a second encoding the binary spins $\sigma_i$ as phase delays. This arrangement supported optical annealing and adiabatic deformation of the Hamiltonian, but it also introduced alignment and calibration burdens (Pierangeli et al., 2020).

A major architectural inflection was the introduction of gauge transformation for SPIM. In the gauge-transformed formulation, both spin configuration and interaction structure are absorbed into a single phase-only SLM by mapping amplitude-like coefficients into transformed spin variables. For a Mattis-type fully connected spin glass, the Hamiltonian is invariant under the transformation to effective spins $\sigma_i\in\{-1,1\}$0, allowing one to encode the transformed problem with a single spatial phase pattern. This eliminated amplitude modulation in several implementations, simplified the optics, improved stability and fidelity, and removed misalignment errors associated with dual-SLM designs (Fang et al., 2020, Huang et al., 2021, Zheng et al., 22 May 2026).

Gauge transformation subsequently became a unifying device across multiple SPIM branches. It was used in the gauge-transformation time-division multiplexing SPIM (GT-SPIM), where a single phase-only SLM encodes both $\sigma_i\in\{-1,1\}$1 and $\sigma_i\in\{-1,1\}$2 and sequentially measures different rank components on the same hardware, and in wavelength-division multiplexed SPIMs, where it removes the pixel-alignment problem that would otherwise arise when different wavelength channels must encode different coupling terms (Zheng et al., 22 May 2026, Luo et al., 2023).

More recent work has emphasized optical fidelity as a systems issue rather than a purely algorithmic one. A full-aperture calibration scheme based on wavefront retrieval and correction with $\sigma_i\in\{-1,1\}$3 accuracy, combined with interaction normalization to compensate amplitude curvature, restored faithful phase encoding across the entire SLM area. In that study, spin-flip diagnostic errors were reduced from as much as 45% to less than 10%, and the energy correlation with the intended Hamiltonian improved from $\sigma_i\in\{-1,1\}$4 to $\sigma_i\in\{-1,1\}$5, enabling sampled lattices up to 9216 spins on a fixed $\sigma_i\in\{-1,1\}$6 active region (Karanikolopoulos et al., 14 Feb 2026). This suggests that scalability in SPIM is conditioned not only by encoding formalism, but also by wavefront control and amplitude uniformity.

3. Interaction classes and the expansion from rank-one to general Hamiltonians

The primitive SPIM naturally realizes only rank-one interaction matrices of Mattis form, a limitation made explicit in several later analyses:

$\sigma_i\in\{-1,1\}$7

This restricted early SPIMs to a narrow class of Ising systems, including number-partitioning-like objectives and Mattis spin glasses, even though the optical Hamiltonian evaluation itself was already highly scalable (Yamashita et al., 2023, Wang et al., 2023).

One response was the multicomponent low-rank model, in which the coupling matrix is decomposed as

$\sigma_i\in\{-1,1\}$8

so the same hardware performs $\sigma_i\in\{-1,1\}$9 optical evaluations and sums them with weights $\phi_j\in\{0,\pi\}$0. This preserves constant-time optical evaluation per component while extending the architecture to rank-$\phi_j\in\{0,\pi\}$1 problems such as knapsack and low-rank Boltzmann machines (Yamashita et al., 2023). A closely related route is interaction-matrix eigendecomposition,

$\phi_j\in\{0,\pi\}$2

implemented through time-division multiplexing. In the general spatial photonic Ising machine based on eigendecomposition, the truncated Hamiltonian Representative Value uses only the largest $\phi_j\in\{0,\pi\}$3 eigenmodes, and for a 20-vertex Max-Cut problem the paper reports that about 65% of the eigenvalues are sufficient to preserve the best optimal-solution probability, whereas performance drops significantly below about 60% (Wang et al., 2023).

A second line of development pursued full programmability rather than low-rank approximation. Wavelength-division multiplexing decomposed a general Ising model into a sum of Mattis-like terms distributed across spectral channels, enabling programmable spin couplings and external magnetic fields on a single SLM and supporting $\phi_j\in\{0,\pi\}$4, Sherrington-Kirkpatrick, and locally connected $\phi_j\in\{0,\pi\}$5 models (Luo et al., 2023). Direct arbitrary-coupling encoding was achieved by assigning SLM pixels to spin products $\phi_j\in\{0,\pi\}$6 and using an ancillary-spin subtraction protocol to cancel unwanted four-spin terms, thereby realizing sparse Hamiltonians with arbitrary couplings and connectivity; the required optical resource scales as $\phi_j\in\{0,\pi\}$7 for sparse graphs rather than $\phi_j\in\{0,\pi\}$8 (Sakellariou et al., 2024). Focal-plane-division SPIM then used eigendecomposition plus simultaneous readout of multiple Mattis energies on different camera regions to produce a fully programmable SPIM for general symmetric Ising matrices (Veraldi et al., 2024).

Other formulations broadened the representable class in different ways. The amplitude-only modulated rank-free SPIM (AR-SPIM) reformulated an arbitrary Ising Hamiltonian with external fields as a sum of Hadamard products and encoded it directly into an incoherent light field with an aligned amplitude SLM and digital micromirror device, reporting 797-spin capacity, nearly 9-bit precision over the range $\phi_j\in\{0,\pi\}$9 to $\sigma_j=\exp(i\phi_j)=\pm1$0, and global encoding metrics of $\sigma_j=\exp(i\phi_j)=\pm1$1 and Pearson correlation coefficient $\sigma_j=\exp(i\phi_j)=\pm1$2 (Zheng et al., 25 Dec 2025). By contrast, the spQUBO formulation argued that even the original non-multiplexed SPIM natively represents a broader class than rank-one if one adopts spatially convolutional QUBOs,

$\sigma_j=\exp(i\phi_j)=\pm1$3

and proved that any spQUBO can be reduced to a two-dimensional periodic spQUBO efficiently implementable on SPIM without multiplexing (Yamashita et al., 30 Jun 2025). Taken together, these works undermine the common misconception that SPIM is intrinsically limited to a rank-one optical Hamiltonian class.

4. Search dynamics, annealing protocols, and algorithmic co-design

Although SPIM evaluates energies optically, most realizations remain hybrid systems in which an external optimization algorithm updates the spin configuration. Simulated annealing (SA) and Metropolis-Hastings updates are particularly common. In one standard protocol, a proposed state $\sigma_j=\exp(i\phi_j)=\pm1$4 is obtained by flipping a single randomly chosen spin and accepted with the Metropolis probability

$\sigma_j=\exp(i\phi_j)=\pm1$5

with $\sigma_j=\exp(i\phi_j)=\pm1$6 and a decreasing temperature schedule (Zheng et al., 22 May 2026). A related adiabatic strategy slowly deforms the Hamiltonian itself,

$\sigma_j=\exp(i\phi_j)=\pm1$7

either by changing the target image (“holographic annealing”) or by changing the optical amplitudes (“optical annealing”), so that the machine follows a low-energy branch from an easy initial Hamiltonian to a frustrated target problem (Pierangeli et al., 2020).

Noise has been treated both as a liability and as a computational resource. In frustrated Mattis spin glasses, controlled readout noise in the feedback loop increases the success probability at a nonzero optimum, with $\sigma_j=\exp(i\phi_j)=\pm1$8 for dense Mattis instances and $\sigma_j=\exp(i\phi_j)=\pm1$9 for sparse ones, whereas in the non-frustrated mean-field model noise is monotonically detrimental (Pierangeli et al., 2020). A related large-scale Max-Cut study added digital white noise to the detected optical signal and interpreted the improvement as noise-assisted escape from local minima (Ye et al., 2023). These results do not imply that noise is universally beneficial; rather, they support a problem-dependent effective-temperature interpretation.

Algorithmic co-design has become increasingly explicit. In a one-SLM number-partitioning SPIM, a direct comparison on a $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$0 lattice found that Metropolis-Hastings reduced the cost to about 1/100 of its initial value within about 800 iterations, whereas a genetic algorithm reduced it only to about half its initial value after about 800 generations and required about twice as long to create a generation; the authors therefore used Metropolis-Hastings for the large-scale number-partitioning experiments (Ramesh et al., 2021). Parallel SPIM (pSPIM) addressed the serial feedback bottleneck by evaluating several Hamiltonians simultaneously through spatial multiplexing; for a 100-spin Max-Cut problem, the annealing time required to reach optimization probability 0.5 decreased from 445 at $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$1 to 164 at $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$2, and the ground state was found 8 times in 50 trials for $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$3 but not found for $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$4 (Shimomura et al., 26 Feb 2025).

The most explicit hybridization appears in the optical genetic-simulated annealing scheme for GT-SPIM. There, GA performs a global coarse-grained search in the early stage and SA performs fine-grained local refinement later. On full-rank Max-Cut simulations, the hybrid reached a Max-Cut of 104.22 in 232 steps at 100 spins, compared with 103.84 in 544 steps for SA and 103.82 in 682 steps for GA; at 500 spins, it reached 1303.39 in 2932 steps, compared with 1202.65 in 1662 steps for SA and 1075.61 in 3507 steps for GA (Zheng et al., 22 May 2026). Space-division multiplexed SPIM further introduced Dynamic Coefficient Search, in which optical coefficients are changed step-by-step during annealing; in the reported knapsack simulations, this outperformed fixed-coefficient SA by steering the search from high-value regions toward feasible regions (Sakabe et al., 2023). A plausible implication is that SPIM performance is increasingly determined by the co-design of optics, encoding, and metaheuristics rather than by optical Hamiltonian evaluation alone.

5. Statistical-physics functionality and phase behavior

SPIMs are not only heuristic optimizers; they have also been used as experimental platforms for equilibrium statistical mechanics. In a 100-spin fully connected SPIM with binary couplings $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$5, the phase diagram was mapped using the disorder-averaged order parameters

$$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$6

The experiment identified paramagnetic, ferromagnetic, and Mattis spin-glass phases, with a reported critical temperature $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$7 and critical disorder probability $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$8, both in good agreement with the corresponding mean-field predictions (Fang et al., 2020). The same work linked the critical region to slower optimization dynamics, because the system requires more iterations to re-equilibrate near the phase transition.

Programmable SPIMs expanded this thermodynamic scope. The wavelength-division multiplexing SPIM implemented $$\mathcal{H}=I(0,0)=\sum_{j,h}\xi_j\xi_h\sigma_j\sigma_h,\qquad J_{jh}=\xi_j\xi_h,$$9 models, Sherrington-Kirkpatrick models, and locally connected $N=4\times10^4$0 models with external magnetic fields, and used the overlap distribution

$N=4\times10^4$1

to resolve paramagnetic, ferromagnetic, spin-glass, and stripe-antiferromagnetic phases. In the SK case, the reported critical temperatures followed the mean-field relations $N=4\times10^4$2 for the PM-FM transition and $N=4\times10^4$3 for the PM-SG transition; in the $N=4\times10^4$4 case, the paper reported a critical temperature around $N=4\times10^4$5 for one parameter set (Luo et al., 2023).

The AR-SPIM study extended phase-transition observations to a 797-spin SK model with arbitrary-rank couplings and external fields. With $N=4\times10^4$6, $N=4\times10^4$7, and $N=4\times10^4$8, the spin-overlap distribution clustered near 0 above $N=4\times10^4$9, became bimodal below $N=74\times74$0, and turned unimodal under a uniform external field $N=74\times74$1 (Zheng et al., 25 Dec 2025). This supports the interpretation that increasingly programmable SPIMs can function as experimental many-body simulators as well as optimization accelerators.

6. Optimization applications, quantitative demonstrations, and present constraints

Combinatorial optimization remains the dominant application domain. Antiferromagnetic SPIM implemented through optoelectronic correlation computing solved a 40,000-spin number-partitioning problem on a single phase-only SLM, reduced the transformed magnetization below $N=74\times74$2 within 100 iterations, and maintained computing fidelity within $N=74\times74$3 from $N=74\times74$4 to $N=74\times74$5 (Huang et al., 2021). A later one-SLM SPIM for two-way number partitioning reached 16,384 spins with about 9 minutes runtime for the maximum reported case, average fidelity $N=74\times74$6, and direct benchmarks against Gurobi and a D-Wave 5000+ Advantage system; the reported D-Wave embedding ceiling for this dense problem was 121 spins, and the paper stated that the achieved fidelities represented a marked improvement of 2 orders of magnitude over current state-of-the-art SPIMs (Ramesh et al., 2021).

Max-Cut has served as the main dense-graph benchmark. A quadrature photonic spatial-Euler Ising machine experimentally solved weighted and unweighted, non-fully connected and fully connected 20,736-node Max-Cut instances on a single phase-only SLM, reporting average cut values about 34% above a Sahni-Gonzales greedy heuristic and about 49–50% improvement at density $N=74\times74$7; for the fully connected case, the paper reported a 122× speedup over the SG CPU comparison (Ye et al., 2023). Direct arbitrary-coupling SPIM solved graph partitioning on $N=74\times74$8, edge density $N=74\times74$9, with reported costs $N=274\times274$0, $N=274\times274$1, and $N=274\times274$2 (Sakellariou et al., 2024). Fully programmable focal-plane-division SPIM reported ground-state success probabilities of 95% for a 16-spin Möbius ladder, 50% for a 16-spin Max-Cut graph, 55% for a 16-spin biased Max-Cut graph, and 10% for a 32-spin fully connected biased Max-Cut graph, while hardware simulations to 100 spins retained favorable approximate-solution statistics (Veraldi et al., 2024). The GA-SA GT-SPIM study extended high-rank Max-Cut experiments to long-iteration runs on 529-spin and 1024-spin instances, obtaining best/mean values of 3,489,410/3,489,347 and 13,109,829/13,109,678, respectively, under a fixed 10,000-iteration budget (Zheng et al., 22 May 2026).

SPIM has also moved into statistical learning. The multicomponent low-rank SPIM was used as a low-rank Boltzmann machine, with a 13-item knapsack optimum appearing 304 times out of 150,000 samples and MNIST classification at $N=274\times274$3 remaining almost flat and comparable to full-rank performance for $N=274\times274$4 (Yamashita et al., 2023). A hybrid optical-digital implementation of Equilibrium Propagation used SPIM for recurrent interaction energy evaluation, achieved $N=274\times274$5 test accuracy on the Wine dataset in experiment versus $N=274\times274$6 in simulation, and numerically reached $N=274\times274$7 on MNIST with $N=274\times274$8 and rank $N=274\times274$9 (Abeele et al., 11 Jun 2026). AR-SPIM additionally reported successful solution of biased Max-Cut problems with arbitrary ranks and weights at less than 0.3% error rate and noted applicability to other QCBO problems including maximum independent set (Zheng et al., 25 Dec 2025).

The present constraints are repeatedly identified across the literature. Conventional SPIM usage can be slow because iterative feedback requires many SLM updates even when each Hamiltonian evaluation is optical (Shimomura et al., 26 Feb 2025). SLM refresh rates remain explicit bottlenecks: one GT-SPIM experiment operated at about 50 Hz, limited by the SLM refresh rate, while another single-SLM number-partitioning implementation reported an effective rate below 7 fps with an SLM response time of about 150 ms (Zheng et al., 22 May 2026, Ramesh et al., 2021). Optical aberrations and non-uniform illumination have historically forced operation to restricted SLM regions (Karanikolopoulos et al., 14 Feb 2026). Rank and precision remain structural bottlenecks for arbitrary Ising problems; analyses of low-rank and circulant constraints emphasize that generic full-rank matrices are expensive, that precision demands can dominate scalability, and that low-rank approximation is problem dependent (Wang et al., 2024). This suggests that future SPIM progress will likely continue along three coupled directions already visible in the literature: improved optical fidelity and calibration, more expressive but implementation-aware Hamiltonian encodings, and tighter integration of optical hardware with annealing, sampling, and learning algorithms (Karanikolopoulos et al., 14 Feb 2026, Luo et al., 2023, Zheng et al., 22 May 2026).