Antiferromagnetic spatial photonic Ising machine through optoelectronic correlation computing

Recently, spatial photonic Ising machines (SPIM) have been demonstrated to compute the minima of Hamiltonians for large-scale spin systems. Here we propose to implement an antiferromagnetic model through optoelectronic correlation computing with SPIM. Also we exploit the gauge transformation which enables encoding the spins and the interaction strengths in a single phase-only spatial light modulator. With a simple setup, we experimentally show the ground state search of an antiferromagnetic model with $40000$ spins in number-partitioning problem. Thus such an optoelectronic computing exhibits great programmability and scalability for the practical applications of studying statistical systems and combinatorial optimization problems.

Recently, like optical analog computations exploring the spatial degrees of freedom [44][45][46][47][48][49][50][51][52], the spatial photonic Ising machine (SPIM) has been proposed with reliable large-scale Ising spin systems, even up to thousands of spins [53]. With spatial light modulations, these spatial Ising spin setups benefit from the high speed and parallelism of optical signal processing [53][54][55][56][57]. Although the modeling of ferromagnetic and spin glass systems has been demonstrated [58], how to implement antiferromagnetic models in SPIM has not been proposed yet. In particular, the antiferromagnetic Ising models are important and extensive in research fields like oxide materials [59] and giant magnetoresistance [60,61]. Also, the combinatorial optimization problems with antiferromagnetic Ising models have many real-world applications such as multiprocessor scheduling, minimization of circuit size and delay, cryptography, and logistics analysis [62,63].
In this Work, we propose to implement the antiferromagnetic model through optoelectronic correlation computing. We show that an antiferromagnetic Hamilto- * These authors contributed equally to this work. † zhichao@zju.edu.cn nian can be evaluated through the correlation between a distribution function and the measured optical intensity with SPIM. We experimentally demonstrate the groundstate-search process with a number-partitioning problem, where 40000 spins are connected with random antiferromagnetic interaction strengths. Our results show that the proposed antiferromagnetic model in SPIM can evolve toward the ground state, exhibiting an efficient approach by scalable degrees of freedom in spatial light modulation.

A.
Optoelectronic correlation computing for Mattis-type Ising model We first introduce the gauge transformation for SPIM, which encodes the spins and the interaction strengths in a single phase-only spatial light modulator (SLM) [58]. We consider a Mattis-type spin glass system with the Ising model Hamiltonian H = − jh J jh σ j σ h , where the spin configuration is S= {σ j } (j = 1, 2, · · · N ) and σ j takes binary value of either +1 or −1, representing the spin-up or spin-down state respectively. The interaction strengths can be expressed as J jh = Jξ j ξ h , where J is an interaction strength as a function of the distance between two spins with the unit of energy, and the amplitude modulation ξ j is limited as −1 ≤ ξ j ≤ 1. By the gauge transformation shown in Fig. 1(a), when rotating each original spin σ j with the angle α j = arccos ξ j , the new spin vector σ ′ j is projected on the z-axis to obtain the effective spin σ ′ The gauge invariance property promises that the experimental implementation only needs a single phase- only SLM, with uniform illumination by a collimated uniform laser beam [ Fig. 1(b)]. This setup circumvents the difficulty of pixel alignment in the previously proposed SPIM [53][54][55][56], and therefore greatly improves the system stability and the computing fidelity. Since the spins are loaded through two-dimensional spatial modulation, the j-th spin is distributed in a square lattice at [58], after the gauge transformation, each spin is encoded by a macropixel with phase modulation ϕ m,n such that σ ′ j z = exp(iϕ m,n ) and Then with Lens L1 of the focal length f performing Fourier transformation, we detect the band-limited intensity distribution I(u) confined within the first diffraction order zone A on the focal plane, and where λ is the wavelength, f is the focal length of lens L1, W is the length of each macropixel, x j = W j is the center position of the j-th pixel, u = (u, v) is the spatial coordinate in the focal plane, and sinc(u) = sin πu πu sin πv πv . Suppose that we preset a distribution function g c (u) and evaluate the correlation function F as Here that is, G(k) is the Fourier transformation of g c (u) sinc 2 ( W u f λ ). Indeed, Eq. (2) shows that by presetting an appropriate g c , a Mattis-type Ising Hamiltonian can be evaluated as We note that the distribution function g c (u) is distinct from the target intensity I T (u) proposed in Ref. [53].
Here g c (u) can be an arbitrary real function to guarantee an even function of the interaction strength [c.f. Eq. (3)], which has either positive or negative values, while the target intensity I T always has non-negative values. In particular, for the antiferromagnetic model, all the interaction strengths J jh < 0. In the case that ξ j s are positive, G(j − h) should be negative to ensure antiferromagnetic interactions between all the spins. It leads that g c must be negative for some values of u, otherwise Eq. (3) shows G > 0 for j = h. Moreover, when g c has both positive and negative values, the Hamiltonian can be evaluated through the correlation function as Eq. (2), while it cannot be implemented by the target-intensity approach.

B. Number-partitioning problem with the antiferromagnetic Hamiltonian
The antiferromagnetic model in SPIM provides a new computation platform for studying the challenging combinatorial optimization problems. As a demonstration, here we present the ground-state-search process of a combinatorial optimization problem, the NP-hard numberpartitioning problem [13,64]: One would like to divide a set Ξ = {ξ j }, containing N real numbers (j = 1, 2, · · · N ), into two subsets Ξ 1 and Ξ 2 , such that the difference between the summations of elements in two subsets 1 = ξj ∈Ξ1 ξ j and 2 = ξj ∈Ξ2 ξ j is as small as possible. Without loss of generality, we suppose all ξ j s in the set Ξ are real numbers belonging to the range (0, 1]. Specifically, when the number set Ξ has parity symmetry, the spins of such models can be analytically optimized. For instance, for the set with an even total number of ξ j s, when ξ j = ξ N +1−j , the spin should be σ j = −σ N +1−j to ensure the equivalence of two subsets. In general, by labeling the elements belonging to two different subsets Ξ 1 and Ξ 2 with σ j = 1 and −1 respectively, the optimization is equivalent to minimizing the antiferromagnetic Hamiltonian H = ( To implement such an antiferromagnetic Hamiltonian in SPIM, we explore the gauge transformation to search a spin configuration S ′ = σ ′ j z where σ ′ j z = ξ j σ j while keeping the interaction strength between any two spins G(k) = −1. Due to the gauge invariance, the Hamiltonian is H = jh σ ′ j z σ ′ h z and the optimized value of During the experimental iterations, the spin configuration is updated gradually [53], and the system definitely evolves to the ground states, indicating the process of solving the optimization problem.

C. Experimental ground state search
We experimentally demonstrate the ground-statesearch process with the number-partitioning problem. As shown in Fig. 1(b), a collimated Gaussian beam (wavelength λ = 532 nm) is expanded by two confocal lenses L2 (50mm focal length) and L3 (500mm focal length). After expansion, the waist radius of the collimated beam is about 36mm. Then a polarizer P is used to prepare the incident beam linearly polarized along the long display axis of the SLM (Holoeye PLUTO-NIR-011). The SLM is calibrated through the two-shot method based on generalized spatial differentiator [65]. Lens L1 with the focal length f = 100mm performs Fourier transformation, where a CCD beam profiler (Ophir SP620) is used to detect the optical field intensity on the back focal plane.
As an example, we demonstrate the searching process with the number-partitioning problem of a set Ξ = {ξ j } with N = 40000 elements. The set Ξ is randomly generated that each element ξ j is a real number randomly chosen from (0, 1], as presented in Fig. 3(a) in the form of a 200-by-200 array. Here each spin is encoded by a macropixel with 2-by-2 pixels on SLM with the length of W = 16µm. Since the beam size is much larger than that of the array, we assume that light illuminates each macropixel with a uniform amplitude.
In order to evaluate the correlation function F , we first numerically calculate the distribution function g c (u) through Eq. (3). Given the interaction strength between any two spins G(k) = −1, Fig. 2(a) shows the calculated distribution function g c (u). We note that g c (u) have both positive and negative values, which means such a case cannot be implemented by the target-intensity approach as Ref. [53]. We also need to calibrate the intensity measurement such that the distribution function g c (u) has the same origin as the intensity distribution I(u). In order to reduce the impact of optical aberrations, we measure the intensity distribution on CCD plane [ Fig. 2(b)] by setting the SLM with uniform phase modulation. Therefore, the origin of I(u) is marked at the maximal intensity through the numerical fitting.
Next we start with the initial state that all spins are uniformly distributed as σ j = 1 and update the spin configuration for ground state search. Here we utilize the Markov chain Monte Carlo algorithm, where σ j s are tentatively updated during each iteration, and the gaugetransformed spins σ ′ j z = ξ j σ j are encoded on SLM following Eq. (1). Then we measure the intensity I(u) on the CCD and evaluate the system Hamiltonian through the correlation function F as Eq. (2). The updated spin configuration is accepted only when the Hamiltonian decreases. Figure 3(b) shows the evolution of the Hamiltonian H and the amplitude of the gauge-transformed magnetization |m ′ | during the ground-state-search process. In the experiment, four independent trials are conducted with the initial state that all spins are uniformly distributed. For all the four cases, the Hamiltonian H and the magnetization |m ′ | decrease rapidly at the beginning of the iteration, because the initial spin configuration strongly deviates from the ground state. As the number of iterations increases, the Hamiltonian tends to be stable, while the |m ′ | starts to fluctuate. We attribute it to the too weak intensity distribution I(u), which is strongly affected by the noise during the measurement. As a result, the spin configuration may be incorrectly updated with the distribution resulting in a larger |m ′ |. The accuracy can be improved by using a wide dynamic-range detector for intensity measurement, or by adjusting the input light intensity in real time within a suitable range.
Here for clear visualization, corresponding to the red square in Fig. 3(a) Overall, for all these four trials, |m ′ | reaches lower than 1.7 × 10 −3 within 100 iterations. Thus during the ground state search, the magnetization |m ′ | decreases by nearly three orders of magnitude, which indicates the validity of the gauge transformation method for antiferromagnetic model.
To evaluate the scalability of the ground state search for number-partitioning problems, we define the computing fidelity as , which is expected to be as small as possible like |m ′ |. We investigate its performance as a function of the size N of the number set and perform experiments for sizes N varying from 1600 to 40000. For each N , we perform 10 independent experimental trials with different sets Ξ = {ξ j }, and then the averaged computing fidelity is obtained by measuring the final states after 1000 iterations. From the results in Fig. 4, we observe that fidelity remains within 6.9 × 10 −3 , demonstrating that the experimental setup works effectively for large-scale number sets. The good scalability of SPIM with gauge transformation on number-partitioning problems is inherited in the parallelism of the optical analog signal processing in the spatial domain.

III. CONCLUSION AND DISCUSSION
We propose to implement antiferromagnetic model in SPIM. By gauge transformation, an antiferromagnetic Hamiltonian can be evaluated through the correlation between the distribution function and the measured optical intensity with SPIM. To improve the processing speed of the system, the ultrafast SLM and CCD at gigahertz rates with the most recent technologies [66,67] is helpful and practical. Also, the computing accuracy can be improved with a more sensitive CCD camera. We note that our proposed method can be applied to the ground-statesearch process, e.g., adiabatic evolution and simulated annealing algorithms [68,69].
In summary, we optically demonstrate the groundstate-search process of an antiferromagnetic Mattis model with thousands of spins, as well as the numberpartitioning problem. With the improved accuracy resulting from gauge transformation, we successfully reduce the total magnetization strength of the gaugetransformed spins |m ′ | by nearly three orders of magnitude. Thus for practical applications in modeling statistical systems and studying combinatorial optimization problems, such an optoelectronic computing exhibits great programmability and scalability in large-scale systems. The authors declare no conflict of interest.