Adiabatic evolution on a spatial-photonic Ising machine

Combinatorial optimization problems are crucial for widespread applications but remain diﬃcult to solve on a large scale with conventional hardware. Novel optical platforms, known as coherent or photonic Ising machines, are attracting considerable attention as accelerators on optimization tasks formulable as Ising models. Annealing is a well-known technique based on adiabatic evolution for ﬁnding optimal solutions in classical and quantum systems made by atoms, electrons, or photons. Although various Ising machines employ annealing in some form, adiabatic computing on optical settings has been only partially investigated. Here, we realize the adiabatic evolution of frustrated Ising models with 100 spins programmed by spatial light modulation. We use holographic and optical control to change the spin couplings adiabatically, and exploit experimental noise to explore the energy landscape. Annealing enhances the convergence to the Ising ground state and allows to ﬁnd the problem solution with probability close to unity. Our results demonstrate a photonic scheme for combinatorial optimization in analogy with adiabatic quantum algorithms and enforced by optical vector-matrix multiplications and scalable photonic technology.


INTRODUCTION
Ising machines are physical devices aimed to accelerate the minimization of Ising Hamiltonians. Scalable implementations of Ising machines are of paramount importance because many of the most challenging combinatorial optimization problems in science, engineering, and social life can be cast in terms of an Ising model [1,2]. Finding the ground state of the Ising spin system gives the solution to the optimization, but requires resources growing exponentially with the problem size. For this reason, intense research focuses on unconventional architectures that use computational units such as light pulses [3], superconducting [4] and magnetic junctions [5], electromechanical modes [6], lasers and nonlinear waves [7][8][9][10][11][12], or polariton and photon condensates [13,14].
Adiabatic computing is a valuable technique to solve combinatorial optimizations by slowly evolving an easyto-prepare initial configuration towards the ground state of a target Hamiltonian, which encodes the combinatorial problem [15,16]. Examples are adiabatic quantum computing using nuclear magnetic resonance [17] and superconducting gates [18], as well as quantum annealing with superconducting circuits [19], and simulated annealing on CMOS networks [20]. Annealing is a form of adiabatic computing at non-zero temperatures in which classical, quantum, or nonlinear perturbations enable the exploration of the complex energy landscape [21][22][23][24]. Quantum annealers designed to solve classical Ising problems succeed in optimization tasks ranging from protein folding [25] to prime factorization [26]. If there are advantages from entanglement and quantum tunneling is still debated [27][28][29]. Therefore, recently large interest is also centred on the the realization of non-electronic annealing devices that can exploits classical nonlinear and photonic properties.
Optical Ising machines use multiple frequency or spa-tial channels to process data at high speed and in parallel. Coherent Ising machines (CIM) employ optical parametric oscillators [30][31][32], fiber lasers [33], or optoelectronic oscillators [34] to solve optimization problems with remarkable performance for hundreds of spins [35]. These machines exploit a gain-dissipative principle [36][37][38][39][40][41]: the search for the minimum energy configuration is conducted in an upward direction by gradually raising the gain. Other platforms based on integrated nanophotonic circuits [42][43][44][45][46][47] operate as optical recurrent neural networks that converge to Ising ground states. Recurrent feedback also has a key role in the recently-demonstrated large-scale photonic Ising machine by spatial light modulation [48][49][50]. In all these optical settings, the possibility of performing a evolution of the machine's parameters to improve the computation remains largely unexplored.
In this Article, we demonstrate adiabatic evolution on a spatial photonic Ising machine (SPIM). We exploit the features of the SPIM, which encodes the spins and their couplings via spatial light modulators (SLMs). Starting from the energy minimum of a simple Hamiltonian, and slowly changing the system, we find the low-energy ground state of a target model. For the adiabatic transformation of the Ising Hamiltonian, we vary the spin couplings at fixed experimental noise level [ Fig. 1(a)]. The annealing protocol occurs optically by amplitude modulation. We also test a different holographic annealing scheme, which uses the image formed during light propagation to control the instantaneous Hamiltonian. We consider Mattis spin glasses with 100-spins initially prepared in a uniform ferromagnetic state. For a sufficiently slow evolution of the Hamiltonian, the success probability approaches unity. Good performance are maintained also when increasing the number of spins. Our findings demonstrate a novel approach that allow applying adiabatic computing principles on photonic devices.  The energetic landscape with several minima varies as the Hamiltonian evolves from the initial H0 to HP corresponding to the target combinatorial problem. When the dynamics is slow enough (adiabatic condition), the system remains in the lowenergy state towards the final ground state, which gives the optimal solution. (b) Realization of the adiabatic evolution on a spatial photonic Ising machine. In the experimental setup, the spins σi are encoded by SLM2 into binary optical phases (the inset shows a representative configuration), whereas SLM1 modulates the beam amplitudes ξi to control spin interaction (top panels). Adiabatic evolution is implemented by the controlled variation of the couplings via the amplitude-modulated light, as illustrated by the random intensity matrices in the top panels (M1-M2, mirrors; L1, lens; CCD, camera).

Spin dynamics on a spatial-photonic Ising machine
In a SPIM, a coherent wavefront encodes binary spin variables by spatial light modulation [48,49]. The device takes advantage of optical vector-matrix multiplications and of the large pixel density of SLMs, properties that enable implementing large-scale photonic computing and machine learning [51][52][53][54][55][56][57]. Figure 1(b) shows the SPIM with an optical path with two SLMs: SLM1 fixes the couplings of the Hamiltonian by amplitude modulation, and SLM2 controls the spin variables (see Methods). Ising spins σ i = ±1 are imprinted on a continuous beam by 0-π phase-delay values (SLM2). The spin interaction occurs by interference on the detection plane. Spatial modulation of the input intensity (SLM1) fixes the interaction strength. Minimizing the difference between the image detected on the camera and a chosen target image I T is equivalent to minimizing an Ising Hamiltonian with couplings determined by the amplitude values set by SLM1 and by I T [48].
Classical annealing is a well-known strategy to minimize problems having many local minima. A timedependent Hamiltonian H(t) is made evolving adiabatically from a simple H(0) = H 0 to the target problem H(T ) = H P [ Fig. 1(a)], with T the annealing time. If the evolution is slow enough, the system remains trapped in the ground state of H(t), and a final measurement provides the spin configuration minimizing H P . In our implementation, we first perform a state preparation phase in which the system reaches the minimum of a Hamiltonian H 0 with homogeneous couplings; then we have the adiabatic evolution, during which the Hamiltonian reaches the target model H P .
State preparation. The optical machine works in a measurement and feedback loop, with all the parameters kept constant once initialized. Recurrent feedback from the detected intensity allows the phase distribution on the SLM2 to converge towards the ground state of an Ising [48]. The quenched variable ξ i is the optical amplitude impinging on the i-th spin;Ĩ T is the Fourier transform of a pre-determined target image. The difference between I T and the image detected on the CCD is the cost function. At each measurement and feedback cycle, we update the spin configuration {σ i } to minimize the cost function. Differently from CIMs [31,34], once initialized, the spin state is updated without electronically computing the energy as well as the field on each spin.
Adiabatic evolution. To realize annealing, we vary the Hamiltonian with time. During the adiabatic phase, H(t) changes towards the target H P . If the evolution is slow enough to prevent high energy excitations, the final state gives the optimal solution. We implement the

time-dependent Ising Hamiltonian
with the conditions H(0) = H 0 and H(T ) = H P . The spin couplings J ij can be changed by the target image I T ("holographic annealing"), and also by varying the modulated amplitudes ξ i ("optical annealing" ). This feature enables different and versatile adiabatic evolution protocols. We remark that, due to the experimental noise, the spin system is coupled to an effective thermal bath [49].

Holographic annealing
A random spin state is prepared in a low-temperature homogeneous ferromagnetic configuration and evolves toward a Mattis spin glass [58,59]. This target H P belongs to a class of frustrated Ising models whose zerotemperature ground states is characterized by ferromagnetic and anti-ferromagnetic domain blocks. During the evolution, time is discretized in M steps. In each step, the Ising machine performs a fixed number of iterations n i using the instantaneous Hamiltonian, until the condition H = H P is reached (see also Methods). The parameter M determines the speed of the adiabatic trajectory, i.e., the total annealing time is T = n i M , following a typical protocol [17].
In the holographic annealing scheme, we vary the target image I T to perform the temporal evolution in Eq. (1), as shown in Fig. 2(a). The graphs for the initial and final problem are inset in Fig. 2(b), where we report the time evolution of the magnetization m = σ i for M = 20. The results correspond to 32 realizations with N = 64 spins. The final ground state always has zero mean magnetization, as expected for a spin glass with vanishing mean interaction. Fluctuations around the average trajectory (thick purple line in Fig. 2(b)) unveil the influence of the experimental noise, which helps the exploration of the various energy configurations during the annealing.
The considered Mattis models enable to directly test the minimization trajectory, as these specific Ising Hamiltonians admit exact zero-temperature ground states [59]. Therefore, we can verify the obtained solutions by inspecting the spin configurations. Figure 2(c) shows the dynamics of the spins. Figure 2(d) reports the final ground state configuration and its comparison with the exact solution. Remarkably, the two spin states have the same energy and differ only for two spin flips. This discrepancy is due to the effective thermal fluctuations. The result indicates that the global minimum is successfully found.

Optical annealing
We implement Eq. (1) by varying the Ising machine parameters optically. Adiabatic evolution is performed by evolving the spatially-modulated intensity inside the photonic machine while keeping constant the target image I T . In each of the M steps, the machine operates with a different amplitude mask ξ i (see also Methods). The instantaneous Hamiltonian H(t) is specified by J ij (t) = ξ i (t)ξ j (t) (within some multiplicative constant).
The initial Hamiltonian H(t) = H 0 corresponds to a homogeneous amplitude distribution [ Fig. 3(a)]. We choose a random set of amplitudes ξ, and gradually decrease their values towards zero. The corresponding system is an instance of a target Mattis model H P with randomly coupled clusters of spin, as represented by the graph inset in Fig. 3(b). Figure 3(b) shows the time-evolving magnetization when approaching the Mattis ground state during the optical annealing. Each trajectory is affected by noise-driven fluctuations, and, at the end of the annealing, the expected mean magnetization m = ±0.5 is reached. This value arises from fact that half of the spins have J ij = 0 in H P . Figure 3(c) shows the evolution of a configuration with N = 64 spins for a representative case. The final ground state averaged over thermal fluctuations coincides with minimum energy state of the programmed Hamiltonian H P . Figure 3(d) shows the remarkable agreement of the obtained spin configuration with the corresponding theoretical zero-temperature solution. This demonstrates that the spins maintain their state at the lowest energy during the optical change of To test the performance of the holographic and optical annealing, we vary the number of discretization steps M , i.e., the annealing time. The number of machine iterations in each step is kept constant. Figure 4(a) shows the success probability p s (see Methods) versus M for the holographic annealing; results refer to the problems in Fig. 2 (N = 64). At a small M , when the evolution occurs rapidly, the excitation of high-energy states reduces the effectiveness of the minimization. The probability of converging to the optimal solution increases with the annealing time. The adiabatic condition is identified by p s reaching a plateau with values exceeding 90% when averaging over thermal fluctuations (magenta bars in Fig. 4(a), see Methods). Figure 4(b) shows the results for the optical annealing scheme. The adiabatic condition is reached on much shorter annealing time with respect to the holographic case in Fig. 4(a). The ground state is found also for rapid processes (M 10). The quality of each solution is given by its Hamming distance, i.e., the number of spins that should be changed to obtain the known zero-temperature ground state [19]. The mean Hamming distance we measure in adiabatic condition is h = 8 ± 1 (h = 7.5 ± 0.8) for the optical (holographic) annealing. In both methods, we found that adiabatic evolution provides a substantial enhancement of the success probability with respect to a "no-annealing" strategy, in which the Hamiltonian H(t) = H P is kept constant since the initial instant [dashed lines in Fig. 4(a,b)]. Moreover, Fig. 4(a) and Fig. 4(b) show that -in all the considered cases -the problem ground state can be found more efficiently if the spin fluctuations at equilibrium are observed and averaged out. This circumstance indicates that a large part of the SPIM error can be ascribed to effective thermal noise.
We also investigate the scaling properties of the photonic setting by varying the number of spins. In Fig. 4(c) and Fig. 4(d), we report the scaling of the success probability for holographic and optical annealing, respectively. For a fast evolution protocol [M = 10, Fig. 4(b)], a degrading effect when increasing the size is observed in the holographic case. However, as the adiabatic condition is reached, we found a remarkable property: good performance are maintained as the problem size grows. For N = 100, values of p s close to unity correspond to a measured mean Hamming fraction (h/N ) of 0.11. The results suggest that, on these specific problems, adiabatic computing can be extended to larger scales with comparable performance.

CONCLUSION
Annealing is one of the most general and consolidated heuristic approaches to solve combinatorial optimization problems and complex physical models. Its experimental implementation on unconventional physical systems operating at room temperature can impact future computing architectures. We have realized adiabatic computing schemes on a spatial-photonic Ising machine showing that ground states are found with enhanced success probability. Computing devices based on spatial light modulation are scalable to larger sizes and can potentially host systems consisting of millions of spins. They can also benefit from temporal fluctuations to speed-up the search for the optimal minimum [49]. Recent SLM technologies can reduce the operation time of our scheme to milliseconds [60]. Other developments may include the use of nonlinear media [50], or metasurfaces, to implement a larger class of combinatorial optimization problems, beyond Mattis instances of the Ising model, and for realizing compact devices without free space propagation. Exploiting optical matrix multiplications, which can be performed efficiently for large sizes [57], our spatial-photonic Ising machine represents a route to tackle hard optimization problems at an unprecedented scale, and opens the route to the experimental demonstration of various minimization strategies.

Experimental setup and feedback method
Light from a continuous-wave laser source with wavelength λ = 532nm (max. power 1.5W) is expanded, polarization controlled and spatially filtered. The beam is thus spatially modulated in amplitude by a first spatial light modulator (SLM1) and then it is independently phase modulated by the second modulator (SLM2). The optical path shown in Fig. 1(b) is realized by a single nematic liquid crystal reflective modulator (Holoeye LC-R 720, 1280 × 768 pixels, pixel pitch 20 × 20µm). A section of the modulator is employed in amplitude mode to generate controlled intensity distributions ξ i , which are imaged by a 4-f system [not shown in Fig. 1(b)] and mirrors M1-M2 on the second section, which perform binary phase modulation. By a combination of incident and analyzed polarizations phase-modulation occurs with less than 10% residual intensity variations. To maintain the setting optically stable, an active area of approximatively 200 × 200 SLM pixels is divided into N optical spins by grouping several pixels. Modulated light is separated using an holographic grating and focused by a lens L1 (f= 500mm) on a CCD camera. The intensity is detected on a region of interest composed of n×m = 18×18 spatial modes, where the signal in each mode is obtained averaging over 10 × 10 camera pixels. The measured intensity pattern determines the feedback signal. At each iteration a spin is randomly flipped, the recorded pattern is compared with a reference image I T on the same number of modes [see Fig. 2(a)], and the spin configuration on the SLM2 is updated to minimize the difference between the two images. Due to intensity fluctuations, as well as to the grouping procedure at the readout, exists a finite probability to update the spin configuration in any case. We experimentally evaluate this probability to p ≈ 0.1. These are the classical fluctuations through which various energy configurations are explored. The value of noise can be also tuned by applying a post-processing step to the recorded intensity [49]. However, in all the presented results it is kept fixed at the same level.

Adiabatic evolution schemes
Holographic annealing. The machine starts from a random configuration of N spins, and it is prepared on a uniform ferromagnetic state. Specifically, the initial Hamiltonian is H 0 = − ijJ σ i σ j , which is implemented using a plane wave of constant amplitude ξ i = E 0 and a target image that is composed by a single spot [ Fig.  2(a)], so thatĨ T = c, being c an arbitrary constant, and J = cE 2 0 . By varying in time the target image I T , being J ij (t) = ξ i ξ jĨT (t), we implement a linear evolution protocol of the form H(t) = (1 − t/T )H 0 + (t/T )H P , with H P an Ising problem where each spin-spin interaction can have only a positive or a negative value (frustrated Mattis model). Therefore, the time-dependent Hamiltonian is where G ij is a block matrix with element values ±1. For example, a four-block matrix G = [G 11 , G 12 ; G 21 , G 22 ], where G 11 = G 22 is a all-ones matrix and G 12 = G 21 = −G 11 , corresponds to a target image I T with two horizontal intensity spots [ Fig. 2(a)]. In this case, pairs of negatively-coupled spins correspond to points of the optical field resulting in destructive interference on the central mode of the detection plane. To discretize the evolution, we divide the total time interval [0, T ] in M equal intervals. For a fixed system size, we choose to keep constant the number of iterations n i performed by the machine in each step. The annealing time thus in-creases with M . In Fig. 2 and Fig. 3, we normalize the evolving time t (machine iterations) with respect to n i . Optical annealing. Starting from a random spin state, a uniform ferromagnetic configuration is prepared by implementing H 0 as in the holographic protocol. However, in the optical evolution method, the target image is maintained fixed to the initial single-spot profile and the whole dynamics is performed by varying the amplitude distribution ξ i via intensity modulation on the SLM1. This corresponds to individually changing the spin-spin coupling values. The evolution protocol we implement is thus where c is an arbitrary constant. Since the interaction matrix is given in any case as a product of two variables J ij ∝ ξ i ξ j , the problem Hamiltonian H P still maintains the form of a Mattis model. In this case, each instance corresponds to a random set of positively interacting spins [ Fig. 3(b)]. The time dependence of the optical amplitudes ξ i (t) can also be selected arbitrarily. We choose to implement couplings that change linearly in time, that is, to variate the optical intensity with a constant rate. This linear schedule corresponds to ξ i (t) = E 0 − (t/T )E 0 , and it is shown in Fig. 3(a) for a random set of evolving ξ i . The total run time T is divided in M equal intervals. Experimentally, the dynamics is implemented by applying an amplitude mask varied in each step. Digital modulation over 256 grey values (8-bit) corresponds to a maximum intensity modulation depth of the order of 10 3 .

Ground states analysis
To quantify the ground state probability at finite temperature, we compute the correlation coefficients between the measured spin configuration and the known optimal solutions of the programmed instance, which for a Mattis graph is identical to the interaction configuration ξ i , or to its sign reversal [59]. The correlation reads as C = i σ i ξ i /|ξ i |, being C = ±1 for the ideal spin system in the lowest energy state. It is successful an evolution whose final equilibrium state gives |C| > 0.75; the success probability p s is the fraction of runs converging to the correct ground state over a set of experiments with different random initial conditions. Instantaneous values of p s are obtained from a single measurement at a fixed time (machine iteration) during the equilibrium stage; values averaged over thermal fluctuations (magenta bars in Fig. 4(a,b)) result from averaging on a time interval the spin configuration measured at equilibrium and identifying the spin with the local magnetization sign. The Hamming distance h is the number of spins that need to be flipped to reach the minimum energy configuration.
The Hamming fraction is the Hamming distance normalized to the total number of spin, h/N .