Experimental Observation of Phase Transitions in Spatial Photonic Ising Machine
aa r X i v : . [ c ond - m a t . d i s - nn ] N ov Experimental Realization of Gauge Transformation and Observation of PhaseTransition in Spatial Optical Spin Systems
Yisheng Fang, ∗ 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,Physics Department, Zhejiang University, Hangzhou 310027, China
Statistical spin dynamics plays a key role to understand the working principle for novel opticalIsing machines. Here we propose the gauge transformations for spatial optical spin systems, where asingle spatial phase modulator simultaneously encodes spin configurations and programs interactionstrengths. Thanks to gauge transformation, we experimentally evaluate the phase diagram of high-dimensional spin-glass equilibrium system with 100 fully-connected spins. We identify that there isa spin-glass phase that the gauge-transformed spins become synchronized when disorder interactionsare rather strong. Furthermore, we exploit the gauge transformation for high accuracy and goodscalability in solving combinatorial optimization problems, for the spin number as large as N =40000. Our results show that the gauge transformation enables a promising approach for solvingNP-hard problems with large-scale and ultra-fast optical computing. I. INTRODUCTION
Recent developments in optical Ising machines havegenerated renewed interest due to the prospects of fea-sibly controlling interaction strengths between artificialspins and utilizing the high speed and parallelism of opti-cal signal processing. In the temporal domain, optimiza-tion problem solvers with remarkable performance havebeen realized with coupled parametric oscillators [1–11],injection-locked or degenerate cavity lasers [12–15], andintegrated nanophotonics circuits [16–22]. In contrast,the spatial photonic Ising machine with optical modula-tion in spatial domain has been demonstrated with reli-able large-scale Ising spin systems, even up to thousandsof spins [23–27]. These setups benefit from the high speedand parallelism of optical signal processing, thus demon-strate high efficiency in solving the ground states of Isingand spin glass Hamiltonians and therefore the combina-torial optimization problems.Generally, statistical spin dynamics plays a key roleto understand the working principle in the spin sys-tems. Although the interaction bonds are randomly dis-tributed between massive spins, the competing interac-tion strengths lead to a new kind of ordered phase, thespin glass phase [28]. Since the first theoretical studyproposed half a century ago [29], understanding the dy-namic phase transitions in spin glass models have beenhighly beneficial for investigations on broad classes of in-teracting systems ranging from statistical physics to com-puter science, including magnetism in condensed matterphysics [28–32], quantum information [33, 34], Hopfieldassociative memory in neural networks [35–37], error-correcting codes and image restoration in informationscience [38, 39], and combinatorial optimization prob- ∗ These authors contributed equally to this work. † [email protected] lems [39, 40]. In particular, a couple of combinatorial op-timization problems with non-deterministic polynomial-time (NP) hard computing complexity, including graphpartition, number partition, travelling salesman prob-lems and others, are directly mapped to the ground-state-searching problems of spin glass models [40, 41].In this Article, we propose and experimentally demon-strate gauge transformation and investigate the phasetransition in spatial optical spin systems. We note thatin order to configure the spatial optical spin system,the spins and the interaction strengths are separatelyencoded by phase and amplitude modulations [23–27],where each pixels in two modulators are required to beaccurately aligned and realize one-to-one mapping. Onthe other hand, such spin systems are strongly relatedto the gauge transformation [39]: The spin configura-tion and the interaction strengths change simultaneouslywhile the total Hamiltonian of the system remains in-variant. This property, termed as gauge invariance, en-ables relieving the disorder in interacting strengths andtransforming it to the effective spin variables. As a fun-damental property, the spin glass systems connected bygauge transformation share the same macroscopic statis-tical properties in equilibrium states, including the phasediagrams and critical behaviors, even though different indetails of microscopic interaction bonds [42–45].By implementing the process of gauge transformation,we incorporate both the spin configuration and interac-tion strengths and integrate these two modulations to-gether in only one phase modulator. There are twofoldadvantages with gauge transformations. First, by per-forming the spatial optical spin system in the equilibriumstates, we explore the phase diagram with the spins num-ber of N = 100. Interestingly, in addition to the standardparamagnetic phase and the ferromagnetic phase, weidentify the spin glass phase where the gauge-transformedspin becomes synchronized when disorder interactions arerather strong. The experimental results for the criticalpoints of phase transitions agree well with the mean-fieldtheory predictions. Second, since with gauge transfor-mation there is no additional amplitude modulation, oneof the advantages is to circumvent the difficulty of pixelalignment, which greatly improve the system stabilityand the computing fidelity. This is of essential impor-tance because measuring the field intensity on the de-tection plane further acts as the feedback to govern theevolution of the spin system [23–26]. Therefore, we ap-ply the optical spin glass system with gauge transforma-tion to solving combinatorial optimization problems. Wedemonstrate the solution of NP-hard number-partitionproblems with high accuracy where the size of the num-ber set ranges from N = 1600 to N = 40000. Our resultsshow that the proposed gauge transformation is greatlybeneficial for high accuracy and good scalability withlarge-scale and ultra-fast optical computing. II. RESULTSA. Gauge transformation in spatial optical spinsystems
To introduce the operation of gauge transformation,we consider the fully-connected spin glass systems withMattis-type interactions. Generally, the spin glassHamiltonian in absence of external field is describedby the classical Ising model H = − P ij J ij σ i σ j . Here, S = { σ j } ( j = 1 , , · · · N ) stands for the spin configura-tion where each σ j takes binary value of either +1 or − J ij , where J ji = J ij . For fully connectedspin systems with Mattis-type interactions [Fig. 1(a)],the interaction strengths are J ij = Jξ i ξ j , determined byΞ = { ξ j } related to all spins. For simplicity, here J isassumed to be fixed and ξ j s are normalized to the region − ≤ ξ j ≤ { ξ j } . Thenthrough the pixel alignment, the modulated beam im-pinges on the phase-only spatial light modulator (SLM),where the spin configuration S = { σ j } is encoded onthe beam wavefront through the phase modulation as ϕ j, SLM = σ j π . Here, ϕ j, SLM is the phase of the j th pixelon SLM, taking binary values of ± π , corresponding tothe spin variables σ j = ±
1, respectively. Subsequently, alens performs Fourier transformation of the optical field,then the field intensity at the focal plane, I ( u, v ), is de-tected by the CCD. Especially, the intensity I at thecenter where ( u, v ) = (0 ,
0) is contributed by the inter-actions between every two spins as I = P ij ξ i ξ j σ i σ j (seeSupplementary Sec. I for details). Thus the Hamiltonian L1 (b) (d) φ j,SLM π = σ j φ j,SLM π = σ j +(-1) j α j (c) feedback CCD
CCD feedback ξ j SLM
SLM1
SLM2
Gauge Transformation (a)
FIG. 1. Schematic of the gauge transformation andthe experimental optical setups for spatial optical spin sys-tems. (a,b) Without gauge transformation, the spin S = { σ j } and the interaction strengths are separately encoded ona phase spatial light modulator (SLM1) and an amplitudeone (SLM2), whose the phase and amplitude modulationare ϕ j, SLM and ξ j , respectively. (c, d) By gauge transfor-mation, the gauge-transformed effective spin configuration S ′ z = (cid:8) σ ′ jz (cid:9) is encoded through only one phase SLM, fol-lowing Eq. (4). The modulated light is detected by a CCDcamera in the back-focus plane of a lens L1. For both twocases, the SLMs are illuminated by collimated laser beams.Details of the experimental setup are presented in the Meth-ods. of such an artificial spin glass system can be defined as H = − J · I (1)Here, we propose the gauge transformation in the casethat the interaction strengths J ij are inhomogeneous andtransform the interaction strengths into the spin orien-tations. As shown in Fig. 1(c), each original spin σ j isrotated clockwise with respect to the z -axis with the an-gle α j = arccos ξ j to arrive at a new spin vector σ ′ j , then σ ′ j is projected on the z -axis to obtain the effective spin σ ′ jz = ξ j σ j and S ′ z = (cid:8) σ ′ j z (cid:9) is the gauge-transformedeffective spin configuration. As the results, the inter-actions between the z components of gauge-transformedspins become uniform in both short and long ranges, withthe strength of J . The gauge transformation operationabove is given as σ j → σ ′ j z , J ij → J. (2)The Hamiltonian remains invariant after gauge transfor-mation, H = − X ij J ij σ i σ j = − X ij Jσ ′ iz σ ′ jz . (3) FIG. 2. Phase diagram of the spatial optical spin sys-tem with 100 spins. (a) The absolute value of magnetiza-tion strength M and spin glass order parameter Q as func-tions of the effective temperature T for fixed probabilities of p = 0 , . , . , . , . , .
5. (b) The absolute value of mag-netization strength M and (c) the square-root value of spinglass order parameter Q as functions of p and T . The phasetransition points are labeled by the critical temperature T c and the critical probability p c . (d) The value of √ Q/ | M | for spin glasses in low effective temperature T = 89 . J , withrespect to the probability p . The results clearly show thetransition between FM and MSG phases. The inset is theschematic of the spin configurations of three different phases. We experimentally implement the gauge transforma-tion with the setup shown in Fig. 1(d). Instead of sepa-rately encoding the interaction strengths J ij = Jξ i ξ j andthe original spin configuration S = { σ j } on two SLMs,here we encode the gauge-transformed spin configuration S ′ z = { σ ′ jz } on a single SLM. In this case, the rotation ofspins corresponds to the modification of the phase mod- ulation on SLM as ϕ j, SLM = σ j π − j α j . (4)Equation (4) is the main result of the present work,where the interaction strengths are incorporated into thespin configuration as the second term. By performingthe gauge transformation, the requirement of amplitudemodulation as shown in Fig. 1(b) is eliminated. Thederivation of Eq. (4) is inspired by the complex encod-ing method with double-phase hologram [46–49] and thedetails are given by Supplementary Sec. II. Due to thegauge invariance, the detected Hamiltonian before andafter the gauge transformation remain invariant, as pre-sented in Eq. (3). B. Observation of phase transition
In order to show the validation of the gauge trans-formation, we first experimentally study the statisticalproperties of spatial optical spin system in the equilib-rium states. The interaction strengths J ij are determinedby Ξ = { ξ j } related to all spins, where each ξ j is ran-domly chosen following the distribution probability of p ( ξ j ) = p · δ ( ξ j + 1) + (1 − p ) · δ ( ξ j −
1) (5)i.e. each ξ j independently takes the value of either − p and 1 − p respectively.With respect to the probability p and an effective tem-perature T , the phases of this spin system are charac-terized by two statistical order parameters, the magne-tization strength M and the spin glass order parameter Q . These two statistical order parameters are definedas M = "* N P j σ j + and Q = * N P j σ j + respec-tively, where h· · · i denotes the ensemble average over thespin configurations S and [ · · · ] denotes the configura-tional average over different sets of Ξ = { ξ j } generatedfollowing the probability p . In order to show the phasediagram of spatial optical spin system, we first formulatethe statistical ensembles containing sufficient samples ofspin configuration S .The proposed gauge transformation is rather conve-nient for investigating the phase transitions of these op-tical spin models. By this way, only one set of experimentis needed to compute the full phase diagram of the spinglass system. This is because though the interactionsof the spin systems are different determined by differentΞ = { ξ j } , they share the same model after gauge trans-formation. In the experiment, it is the gauge-transformedeffective spin configurations S ′ = (cid:8) σ ′ jz (cid:9) that are encodedto formulate the statistical ensembles containing suffi-cient samples of spin configurations. With the knowledgeof any given Ξ = { ξ j } , the statistical ensemble for theoriginal spin configurations S before gauge transforma-tion can be obtained by simply performing σ j = σ ′ jz /ξ j ,where for each given probability p , 100 different sets ofΞ = { ξ j } are generated following the probability equation[Eq. (5)] for configurational averaging in determining theorder parameters M and Q of the spin glass systems.In the experiment, an arbitrarily initial effective spinconfiguration S ′ = (cid:8) σ ′ j z (cid:9) with 100 spins, where each σ ′ jz takes value of either 1 or -1 randomly, is encoded onthe phase-only SLM (Holoeye PLUTO-NIR-011, 1920 × µ m × µ m) following ϕ j, SLM = σ ′ j z π . Here an active area with the size of3200 µ m × µ m on SLM is divided in to an array of10-by-10 macropixels, where each macropixel encode aneffective spin.We generate the ensemble of equilibrium states inthe optical spin glass system with each given effectivetemperature T following the measurement and feedbackscheme. Governed by the Markov Chain Monte Carlo(MCMC) algorithm, we tentatively flip one spin on theSLM during each iteration and measure the field in-tensity I ( u, v ) on the back focal plane of the lens (fo-cal length f = 100mm) by CCD (Ophir SP620), whichthen gives the Hamiltonian as Eq. (1). Following theMetropolis-Hasting sampling procedure (see the Meth-ods for details), with the effective temperature T , weupdate the spin configuration with the knowledge ofthe optically computed feedback Hamiltonian H . Un-der each fixed effective temperature T , 1000 MonteCarlo iterations are performed to formulate a statisti-cal ensemble containing 1000 effective spin configurations { S ′ ( T )1 , S ′ ( T )2 , · · · , S ′ ( T )1000 } . These 1000 samples representthe equilibrium state of the spin glass system, where thespin configurations satisfy the Gibbs-Boltzmann distri-bution, p ( S ′ ) ∝ e − H/T . We gradually cool the effectivetemperature from T = 400 J to T = 89 . J in the simu-lated annealing manner. In this way, the whole statisticalequilibrium states in the effective temperature range areobtained. The experimentally measured order parame-ters are presented in Fig. 2, forming the phase diagramwith respect to p and T . Here five experimental runsare conducted independently and then averaged to re-duce the error. In general, following the processes de-scribed above, we can determine the full phase diagramof the spatial optical spin system with the aid of thegauge transformation.As shown in Fig. 2, these two order parameters M and Q categorize the spin glass systems into three dis-tinct phases, the paramagnetic phase (PM, M = Q = 0),the ferromagnetic phase (FM, M > Q >
0) and theMattis spin glass phase (MSG, M = 0, Q > T for fixed probabilities p . When p = 0 for example, thespin system has uniform interactions, and it experiencesthe order-disorder phase transition, at a critical tempera-ture T c = 199 J , between the PM and FM phases duringannealing. At the high effective temperature T > T c ,the original spins σ j s align randomly, thus the averagedmagnetization strength M vanishes in the paramagneticphase. At low effective temperature T < T c , the spins σ j tend to align in parallel to minimize the total energy H , resulting in spontaneous magnetizations with nonzero M in the ferromagnetic phase. This phase transition be-tween PM and FM during annealing still holds for otherspin glass systems where the degree of disorder in inter-actions is lower than a critical value p c .When p > p c however, a phase of MSG occurs at loweffective temperature, where the original spins σ j s seemto be randomly distributed, but they are in fact lockedto ξ j s, that is, all σ j s tend to align parallel to ξ j s. Thusthe gauge-transformed spins σ ′ jz s all take uniform valueof +1 or − σ j s originates from the dis-order in the interaction strengths, i.e. the disorder of ξ j s, which is different from the case in the PM phase.When the interaction strengths given by { ξ j } are fixed,there’s no randomness in the value of { σ j } and the en-semble averaged magnetization strength * N P j σ j + isnon-vanishing ( Q > M = 0).In order to determine the critical probability p c , wedistinguish the FM and MSG phases by plotting the val-ues of √ Q/ | M | with respect to the probability p , at a loweffective temperature T = 89 . J . The data are extractedfrom the left-most lines of Fig. 2(b) and (c). The re-sults are shown in Fig. 2(d), and we determine p c = 0 . √ Q/ | M | remains approximately 1 for p < p c anddiverges when p > p c , showing the characteristics of FMand MSG phases respectively. We note that the criticalpoints T c = 199 J and p c = 0 .
45 for phase transitions,determined from the phase diagram, agree well with thepredictions of the mean field theory, T c = 2( N − J and p c + (1 − p c ) = 1 / (1 + e − J/T c ) [39]. Indeed, Figs.2(b) and (c) show the fluctuations of the order param-eters in the vicinity of the critical points, which comefrom the critical-slowing-down phenomenon in statisti-cal models. These results show the validation that ourscheme can model the spatial optical spin systems withhigh accuracy, even when the optical system has someaberrations. C. Application to number-partition problems
The spatial optical spin system provides a new compu-tation platform for solving the challenging combinatorialoptimization problems. More importantly, the process ofgauge transformation circumvents the difficulty of pixelalignment and therefore greatly improves the system sta-bility and the computing fidelity. As a demonstration,here we present examples of solving a combinatorial op-timization problem, the NP-hard number-partition prob-lem [40, 41]: A set containing N real numbers Ξ = { ξ j } ( j = 1 , , · · · N ) is divided into two subsets Ξ and Ξ ,such that the difference between the summations of ele- -11 H | m ' | -11 (a) (c) (d)(b) iterations ×10 initial I ( u,v ) final I ( u,v ) -4 -2 FIG. 3. Experimental results for the number-partition prob-lem. (a) The number set Ξ = { ξ j } containing 40000 elements,where each ξ j is chosen randomly from [0 ,
1] with uniformprobability. (b) Evolutions of the Hamiltonian H and am-plitude of the gauge-transformed magnetization | m ′ | duringthe iterations. Four independent trials are conducted withuniform initial spin configurations. Inset images show the de-tected intensity I ( u, v ) on focal planes for the initial state(left) and the final state (right). The white bar correspondsto the length of 100 µ m. For clear visualization, (c) and (d)respectively presents only a part of the final configurationsof the gauge-transformed spins (cid:8) σ ′ jz (cid:9) on SLM and the finalpartition solutions { σ j } , which correspond to the red squarebox in (a). The complete solutions can be found in the Sup-plementary Sec. IV. ments in two subsets P = P ξ j ∈ Ξ ξ j and P = P ξ j ∈ Ξ ξ j is as small as possible. Without loss of generality, wesuppose all ξ j s in the set Ξ are real numbers belongingto the region [0 , S ′ = (cid:8) σ ′ j z (cid:9) where σ ′ j z = ξ j σ j . Here the absolutevalues of the gauge-transformed spins (cid:12)(cid:12) σ ′ jz (cid:12)(cid:12) = ξ j repre-sent the elements in the number set Ξ, while σ j = 1and − and Ξ respectively. The optimization is tominimize the value | P − P | , which is related to thetotal magnetization strength of the gauge-transformedspins | m ′ | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N P j σ ′ j z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) as | P − P | = N | m ′ | . Suchgoal is equivalent to minimizing the Hamiltonian H =( P j ξ j σ j ) = P ij ξ i ξ j σ i σ j = P ij σ ′ iz σ ′ j z . Therefore, solv-ing the number-partition problem can be mapped to the ground-state-searching process for the spatial optical spinsystem, where J ij = ξ i ξ j and J = − T = 0, the system definitely evolvesto the ground states, thus minimizing the Hamiltonian H , and then the combinatorial optimization problem issolved. number of spins N f i de li t y −4 −3 − − FIG. 4. The computing fidelities for number sets with sizesof N = 1600 , , , , , . We first solve the number-partition problem for a setΞ = { ξ j } with N = 40000 elements. The set Ξ is ran-domly generated that each element ξ j is a real numberrandomly chosen from [0 , σ j = 1, corresponding to the sit-uation that all elements are put into the same partitionsubset. During the iterations, σ j s are tentatively flippedand the gauge-transformed spins σ ′ jz = ξ j σ j are encodedon the active area with the size of 3200 µ m × µ mon SLM following Eq. (4). Here each σ ′ jz is encodedby a macropixel with 2-by-2 pixels on SLM with thelength of W = 16 µ m. Figure 3(b) shows the evolu-tion of the Hamiltonian H and amplitude of the gauge-transformed magnetization | m ′ | during the ground-state-searching process. The computing fidelity in solving suchproblems is defined as η = (cid:12)(cid:12)(cid:12) P − P P + P (cid:12)(cid:12)(cid:12) and it is expectedto be as small as possible. In the experiment, four inde-pendent trials are conducted with the initial state thatall spins are uniformly distributed. For all these four tri-als, η reach lower than 3 . × − within 100 iterations.The slightly non-perfect fidelities result from the fact thatthe glassy systems are easily trapped in the meta-stablestates during evolution.Here for clear visualization, corresponding to the redsquare indicated in Fig. 3(a), Figs. 3(c) and (d) presentonly a part of the final configurations of the gauge-transformed spins S ′ = (cid:8) σ ′ j z (cid:9) and the resulted partitionsolutions { σ j } , respectively. The complete solutions canbe found in the Supplementary Sec. V. These four inde-pendent trials reach four different solutions belonging tothe highly-degenerate solution sets. For the final steadystate, the detected image presents as a dark spot in thecenter, I = 0, shown as the inset of Fig. 3(b), whichindicates the minimized Hamiltonian in Eq. (3).To evaluate the scalability of the spatial opticalspin system, we investigate its performance on solvingnumber-partition problems as a function of the size ofthe number set N . We perform experiments for sizes N varying from 1600 to 40000. For each N , we per-form 10 independent experimental trials with differentsets Ξ = { ξ j } , and then the averaged computing fidelityis obtained by measuring the final states after 1000 itera-tions. From the results in Fig. 4, we can see that fidelityremains within 6 . × − , demonstrating that the ex-perimental setup works effectively for large-scale numbersets. The good scalability of spatial optical spin systemwith gauge transformation is inherited in the parallelismof the optical analog signal processing in the spatial do-main. III. CONCLUSION AND DISCUSSION
Benefiting from the gauge transformation, we charac-terize the phase transition to understand the statisticalspin dynamics in the novel optical spin system. We alsonotice a study to simulate large-scale random spin net-works with a disordered medium [26]. However, the tech-nique is distinct from the gauge transformation and can-not characterize the phase transition presented here. Fur-thermore, based on the Ising systems constructing withnonlinear processes as shown in [27], we also believe thatthe proposed gauge transformation can be extended tostudy more complex spin models with four-body interac-tions.In summary, we have proposed a spatial optical spinsystem with gauge transformation. The encoding of thespin configurations and the programming of the inter-action strengths between spins are realized by only onespatial phase modulation process, which significantly im-proves the stability and fidelity of the optical Ising ma-chine. For both studies of the statistical equilibrium stateproperties of spin glass systems and practical applicationsin solving combinatorial optimization problems, such anoptical system exhibits high accuracy and great robust-ness against noises and aberrations, intriguing for its pro-grammability and scalability in large-scale Ising machine.
IV. ACKNOWLEDGEMENT
The authors acknowledge funding through the Na-tional Natural Science Foundation of China (NSFCGrants Nos. 91850108 and 61675179), the National KeyResearch and Development Program of China (Grant No.2017YFA0205700), the Open Foundation of the StateKey Laboratory of Modern Optical Instrumentation, and the Open Research Program of Key Laboratory of 3D Mi-cro/Nano Fabrication and Characterization of ZhejiangProvince.
V. APPENDIX: METHODSA. Experimental setup and measurement of systemHamiltonian
In the experiment, we use a green laser source (wave-length λ = 532 nm) to generate a collimated beam witha planar wavefront and a uniform amplitude distribution(see the detailed setup in Fig. S1). Here the collimatedlaser with a beam waist radius of about 3.6mm is ex-panded by lenses L2(focal length is 50mm) and L3(focallength is 500mm), which generates a sufficiently wideGaussian beam to cover the phase-only SLM (HoloeyePLUTO-NIR-011), such that the amplitude distributionat the used region of SLM is rather uniform. The po-larizer P1 is used to prepare the incident beam linearlypolarized along the long display axis of the SLM. A CCD(Ophir SP620) is used to detect the optical field intensityon the back focal plane.We note that with a finite pixel-size detector, itis hard to exactly detect I at the center, due tothe restriction of the NA of lenses and the resolu-tion of CCD. So the Hamiltonian is calculated bynormalizing the detected intensity in a finite regionas H = − J RR A [ I ( u, v ) × ¯ I ( u, v )]d u d v/ RR A ¯ I ( u, v )d u d v ,where ¯ I ( u, v ) is the field intensity on the detectionplane when the SLM has uniform phase modulation of ϕ j, SLM = 0 in the squared detection region A : | u | , | v | ≤ d/ d > . fλ W . Here f is the focal length of lens L1shown in Fig. 1(d), and W corresponds to the length ofa macropixel on SLM encoding the effective spin config-urations. B. Optical Metropolis-Hasting sampling
For the purpose of studying the statistical propertiesof spin glass systems, we first formulate the statisticalensemble containing sufficient samples of spin configura-tion S . The statistical ensembles are obtained by opti-cal iterations governed by Monte Carlo algorithm withMetropolis-Hasting sampling. In each iteration, a sin-gle spin is flipped and the updated spin configuration S new is accepted with the probability of e − ∆ H/T , de-pending on both the variance of the energy function∆ H = H ( S new ) − H ( S ) and the effective temperatureas follow: ( if ∆ H > S new with the probability of e − ∆ H/T if ∆ H S new . The authors declare no conflict of interest. [1] P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly,C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Ut-sunomiya, K. Aihara, et al. , A fully programmable 100-spin coherent Ising machine with all-to-all connections,Science , 614 (2016).[2] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Ta-mate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki,K. Enbutsu, et al. , A coherent Ising machine for 2000-node optimization problems, Science , 603 (2016).[3] T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Ya-mamoto, and H. Takesue, Large-scale Ising spin networkbased on degenerate optical parametric oscillators, Na-ture Photonics , 415 (2016).[4] F. B¨ohm, G. Verschaffelt, and G. Van der Sande, A poorman’s coherent Ising machine based on opto-electronicfeedback systems for solving optimization problems, Na-ture Communications , 1 (2019).[5] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli,A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba,T. Honjo, et al. , Experimental investigation of perfor-mance differences between coherent Ising machines and aquantum annealer, Science Advances , eaau0823 (2019).[6] L. Bello, M. Calvanese Strinati, E. G. Dalla Torre, andA. Pe’er, Persistent coherent beating in coupled para-metric oscillators, Physical Review Letters , 083901(2019).[7] F. B¨ohm, T. Inagaki, K. Inaba, T. Honjo, K. Enbutsu,T. Umeki, R. Kasahara, and H. Takesue, Understandingdynamics of coherent Ising machines through simulationof large-scale 2D Ising models, Nature Communications , 1 (2018).[8] E. S. Tiunov, A. E. Ulanov, and A. Lvovsky, Annealingby simulating the coherent Ising machine, Optics Express , 10288 (2019).[9] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Ya-mamoto, Coherent Ising machine based on degenerateoptical parametric oscillators, Physical Review A ,063853 (2013).[10] Y. Yamamoto, K. Aihara, T. Leleu, K.-i. Kawarabayashi,S. Kako, M. Fejer, K. Inoue, and H. Takesue, CoherentIsing machines optical neural networks operating at thequantum limit, npj Quantum Information , 1 (2017).[11] A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Ya-mamoto, Network of time-multiplexed optical parametricoscillators as a coherent Ising machine, Nature Photonics , 937 (2014).[12] S. Utsunomiya, K. Takata, and Y. Yamamoto, Mappingof Ising models onto injection-locked laser systems, Op-tics Express , 18091 (2011).[13] K. Takata, S. Utsunomiya, and Y. Yamamoto, Transienttime of an Ising machine based on injection-locked lasernetwork, New Journal of Physics , 013052 (2012).[14] M. Babaeian, D. T. Nguyen, V. Demir, M. Akbulut, P.-A. Blanche, Y. Kaneda, S. Guha, M. A. Neifeld, andN. Peyghambarian, A single shot coherent Ising machine based on a network of injection-locked multicore fiberlasers, Nature Communications , 1 (2019).[15] C. Tradonsky, I. Gershenzon, V. Pal, R. Chriki, A. A.Friesem, O. Raz, and N. Davidson, Rapid laser solverfor the phase retrieval problem, Science Advances ,eaax4530 (2019).[16] C. Roques-Carmes, Y. Shen, C. Zanoci, M. Prabhu,F. Atieh, L. Jing, T. Dubˇcek, C. Mao, M. R. John-son, V. ˇCeperi´c, et al. , Heuristic recurrent algorithms forphotonic Ising machines, Nature Communications , 1(2020).[17] M. Prabhu, C. Roques-Carmes, Y. Shen, N. Har-ris, L. Jing, J. Carolan, R. Hamerly, T. Baehr-Jones,M. Hochberg, V. ˇCeperi´c, et al. , Accelerating recurrentIsing machines in photonic integrated circuits, Optica ,551 (2020).[18] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle,D. Englund, et al. , Deep learning with coherent nanopho-tonic circuits, Nature Photonics , 441 (2017).[19] K. Wu, J. G. De Abajo, C. Soci, P. P. Shum, andN. I. Zheludev, An optical fiber network oracle for NP-complete problems, Light: Science & Applications ,e147 (2014).[20] M. R. V´azquez, V. Bharadwaj, B. Sotillo, S.-Z. A. Lo,R. Ramponi, N. I. Zheludev, G. Lanzani, S. M. Eaton,and C. Soci, Optical NP problem solver on laser-writtenwaveguide platform, Optics Express , 702 (2018).[21] Y. Okawachi, M. Yu, J. K. Jang, X. Ji, Y. Zhao, B. Y.Kim, M. Lipson, and A. L. Gaeta, Nanophotonic spin-glass for realization of a coherent Ising machine, arXivpreprint arXiv:2003.11583 (2020).[22] M. Prabhu, C. Roques-Carmes, Y. Shen, N. Har-ris, L. Jing, J. Carolan, R. Hamerly, T. Baehr-Jones,M. Hochberg, V. ˇCeperi´c, et al. , A recurrent Ising ma-chine in a photonic integrated circuit, arXiv preprintarXiv:1909.13877 (2019).[23] D. Pierangeli, G. Marcucci, and C. Conti, Large-scalephotonic Ising machine by spatial light modulation,Physical Review Letters , 213902 (2019).[24] D. Pierangeli, G. Marcucci, D. Brunner, and C. Conti,Noise-enhanced spatial-photonic Ising machine,Nanophotonics (2020).[25] D. Pierangeli, G. Marcucci, and C. Conti, Adiabatic evo-lution on a spatial-photonic Ising machine, arXiv preprintarXiv:2005.08690 (2020).[26] D. Pierangeli, M. Rafayelyan, C. Conti, and S. Gi-gan, Scalable spin-glass optical simulator, arXiv preprintarXiv:2006.00828 (2020).[27] S. Kumar, H. Zhang, and Y.-P. Huang, Large-scale Isingemulation with four body interaction and all-to-all con-nections, Communications Physics , 1 (2020).[28] K. Binder and A. P. Young, Spin glasses: Experimentalfacts, theoretical concepts, and open questions, Reviewsof Modern Physics , 801 (1986). [29] S. F. Edwards and P. W. Anderson, Theory of spinglasses, Journal of Physics F: Metal Physics , 965(1975).[30] D. Sherrington and S. Kirkpatrick, Solvable model of aspin-glass, Physical Review Letters , 1792 (1975).[31] M. Gabay and G. Toulouse, Coexistence of spin-glass andferromagnetic orderings, Physical Review Letters , 201(1981).[32] J. D. Reger and A. Zippelius, Three-dimensional random-bond Ising model: Phase diagram and critical properties,Physical Review Letters , 3225 (1986).[33] G. De las Cuevas, W. Dur, H. J. Briegel, and M. A.Martin-Delgado, Unifying all classical spin models ina lattice gauge theory, Physical Review Letters ,230502 (2009).[34] M. C. Ba˜nuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac,M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein,S. Montangero, et al. , Simulating lattice gauge theo-ries within quantum technologies, The European PhysicalJournal D , 1 (2020).[35] J. J. Hopfield, Neural networks and physical systems withemergent collective computational abilities, Proceedingsof the National Academy of Sciences , 2554 (1982).[36] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Spin-glass models of neural networks, Physical Review A ,1007 (1985).[37] E. Agliari, A. Barra, A. De Antoni, and A. Galluzzi, Par-allel retrieval of correlated patterns: From Hopfield net-works to Boltzmann machines, Neural Networks , 52(2013).[38] N. Sourlas, Spin-glass models as error-correcting codes,Nature , 693 (1989).[39] H. Nishimori, Statistical physics of spin glasses and in-formation processing: an introduction , 111 (ClarendonPress, 2001).[40] A. Lucas, Ising formulations of many NP problems, Fron-tiers in Physics , 5 (2014).[41] S. Mertens, Phase transition in the number partitioningproblem, Physical Review Letters , 4281 (1998).[42] H. Nishimori, Exact results and critical properties ofthe Ising model with competing interactions, Journal ofPhysics C: Solid State Physics , 4071 (1980).[43] H. Nishimori, Internal energy, specific heat and correla-tion function of the bond-random Ising model, Progressof Theoretical Physics , 1169 (1981).[44] Y. Ueno and Y. Ozeki, Interfacial approach tod-dimensional ± J Ising models in the neighborhood ofthe ferromagnetic phase boundary, Journal of StatisticalPhysics , 227 (1991).[45] R. R. P. Singh, Spin-glass–ferromagnetic–paramagneticmulticritical point, Physical Review Letters , 899(1991).[46] C. Hsueh and A. Sawchuk, Computer-generated double-phase holograms, Applied Optics , 3874 (1978).[47] O. Mendoza-Yero, G. M´ınguez-Vega, and J. Lancis, En-coding complex fields by using a phase-only optical ele-ment, Optics Letters , 1740 (2014).[48] S. Ngcobo, I. Litvin, L. Burger, and A. Forbes, A digitallaser for on-demand laser modes, Nature Communica-tions , 1 (2013).[49] A. Dudley, R. Vasilyeu, V. Belyi, N. Khilo, P. Ropot,and A. Forbes, Controlling the evolution of nondiffract-ing speckle by complex amplitude modulation on a phase-only spatial light modulator, Optics Communications , 5 (2012). Supplementary Information: Experimental Realization of Gauge Transformation andObservation of Phase Transition in Spatial Optical Spin Systems
I. THE EXPERIMENTAL SETUP WITHOUT GAUGE TRANSFORMATION, AS PRESENTED INFIG. 1(B)
The spin configuration S = { σ j } is encoded on the SLM in an array of N x × N y macropixels. The size of a singlemacropixel is W × W. As presented in Fig. 1(b), a paraxial beam with amplitude modulation { ξ j } illuminates on theSLM, which has the spatial phase modulation of { σ j } . After reflected by SLM, the electric field is E ( x ) = i X j ξ j σ j rect W ( x − x j ) = [ i X j ξ j σ j δ ( x − x j )] ⊗ rect W ( x ) . (S1)Here x = ( x, y ) denotes the spatial coordinate on the SLM plane, and x j is the center position of the j th pixel onSLM, which takes the value of x j = W( m e x + n e y ) (1 ≤ m ≤ N x , ≤ n ≤ N y ).In Eq. (S1), the notation ⊗ is the convolution operation, and the rectangular functions are defined asrect W ( x ) = rect( x W ) = (cid:26) | x | , | y | ≤ W / | x | , | y | > W / x ) = (cid:26) | x | , | y | / | x | , | y | > / E ( x ), E ( x ) FT −−→ ˜ E ( k ) = π ) R + ∞−∞ E ( x ) e − i x · k d x .Specifically, P j ξ j σ j δ ( x − x j ) FT −−→ ( ) P j ξ j σ j e i k · x j rect W ( x ) FT −−→ W sinc W ( k )where sinc W ( k ) = sinc( k W2 π ) = sin W kx Wkx · sin Wky Wky , sinc( k ) = sin πk x πk x sin πk y πk y .Therefore we have ˜ E ( k ) = ( i P j ξ j σ j e i k · x j ) · sinc W ( k ).The electric field on the detection plane, i.e. the back focal plane, is E ( u, v ), where u = ( u, v ) is the spatialcoordinate on the detection plane. Here the focal length of the FT lens is f and the wavelength of the laser source is λ , and u = k x fλ π , v = k y fλ π . So E ( u, v ) = ( i X j ξ j σ j e i πfλ u · x j ) · sinc( u W f λ ) . (S2)The detected intensity image on CCD is I ( u, v ) = E ∗ ( u, v ) · E ( u, v ) = X ij ξ i ξ j σ i σ j e i πfλ ( x i − x j ) · u sinc ( u W f λ ) . (S3)0The Hamiltonian is defined as H = − J · I = − J · X ij ξ i ξ j σ i σ j , (S4)where J is a constant and I is the intensity on the center of the detection plane, ( u, v ) = (0 , FIG. S1. Experimental setup of the spatial optical spin system with gauge transformation.
II. DERIVATION OF GAUGE TRANSFORMATION IN EQ. (4)
Supposing that a collimated beam, with uniform amplitude, impinges on the phase-only SLM with the spatial phasemodulation following Eq. (4) in the text, the beam wavefront becomes E ( x ) = iM ( x )[ X j σ j e iα j rect W ( x − x j )] + i (1 − M ( x ))[ X j σ j e − iα j rect W ( x − x j )]= i [ H ( x ) M ( x ) + H ( x )(1 − M ( x ))] ⊗ rect W ( x ) (S5)where H ( x ) = P j σ j e i cos − ξ j δ ( x − x j ) H ( x ) = P j σ j e − i cos − ξ j δ ( x − x j ) H ( x ) + H ( x ) = 2 P j ξ j σ j δ ( x − x j ) H ( x ) − H ( x ) = 2 i P j q − ξ j σ j δ ( x − x j )and M ( x ) and 1 − M ( x ) are the checkerboard pattern functions as M ( x ) = + ∞ P m,n = −∞ [ δ ( x − m W , y − n W) + δ ( x − (2 m + 1)W , y − (2 n + 1)W)]1 − M ( x ) = + ∞ P m,n = −∞ [ δ ( x − m W , y − (2 n + 1)W) + δ ( x − (2 m + 1)W , y − n W)] .The values of M ( x ) ⊗ rect W ( x ) and (1 − M ( x )) ⊗ rect W ( x ) are presented in Fig. S2. Specifically, in the case that all ξ j s take binary values of either +1 or -1, the modulated optical field is simplified as E ( x ) = i X j σ ′ iz rect W ( x − x j ) . (S6)According to the Fourier optics, the lens L1(with focal length f = 100mm) performs the spatial Fourier transforma-tion(FT) on E ( x ),1 (a) (b) x y x y FIG. S2. Checkerboard patterns for functions (a) M ( x ) ⊗ rect W ( x ) and (b) (1 − M ( x )) ⊗ rect W ( x ). H ( x ) FT −−→ h ( k x , k y ) H ( x ) FT −−→ h ( k x , k y ) M ( x ) FT −−→ ( 12W ) ∞ X m,n = −∞ (1 + e ik x W e ik y W ) δ ( k x − m π W , k y − n π W )= ( 12W ) ∞ X m,n = −∞ (1 + ( − m + n ) δ ( k x − m π W , k y − n π W )1 − M ( x ) FT −−→ ( 12W ) ∞ X m,n = −∞ ( e ik x W + e ik y W ) δ ( k x − m π W , k y − n π W )= ( 12W ) ∞ X m,n = −∞ (( − m + ( − n ) δ ( k x − m π W , k y − n π W ).So the spatial spectrum is˜ E ( k ) = i { h ( k x , k y ) ⊗ [ + ∞ X m,n = −∞ (1 + ( − m + n ) δ ( k x − m π W , k y − n π W )]+ h ( k x , k y ) ⊗ [ + ∞ X m,n = −∞ (( − m + ( − n ) δ ( k x − m π W , k y − n π W )] }· sinc W ( k x , k y ) . (S7)The convolution terms state that the optical field on the focal plane is distributed periodically. This is because thepixels on the SLM have finite sizes, so ˜ E ( k ) has the discrete Fourier transform (DFT)-like phenomena, where differentdiffraction orders are centered at ( k x , k y ) = ( m, n ) π W . The zeroth diffraction order is˜ E ( k ) = i h ( k x , k y ) + h ( k x , k y )] · sinc W ( k x , k y )= ( X j ξ j σ j e i k · x j ) · sinc W ( k )= ( X j σ ′ jz e i k · x j ) · sinc W ( k ) . (S8)In fact, the higher-diffraction-order fields may overlap with the zeroth-order field, thus the total field is the interferencebetween the zeroth and higher diffraction orders,˜ E ( k ) = ˜ E ( k ) + ˜ E ′ ( k ) (S9)where ˜ E ′ ( k ) is the contribution from higher diffraction orders. The field ˜ E ′ ( k ) acts as noises that degrades theperformance of the optical setup in modeling the spin glass systems.2The electric field on the detection plane, related to ˜ E ( k ) through u = k x fλ π , v = k y fλ π , is written as E ( u, v ) = ( X j σ ′ jz e i πfλ u · x j ) · sinc( u W f λ ) . (S10)If H ( x ) and H ( x ) are band-limited such that E ( u, v ) is confined in the first Brillouin zone ( − fλ ≤ u, v ≤ fλ ),thus different diffraction orders do not overlap with each other. In the experiments, we are only interested in thedetected field intensity within the finite area ( − fλ ≤ u, v ≤ fλ ), so only the zeroth diffraction order is containedand thus E ( u, v ) = E ( u, v ).The detected intensity image on CCD is I ( u, v ) = E ∗ ( u, v ) · E ( u, v ) = X ij σ ′ iz σ ′ jz e i πfλ ( x i − x j ) · u sinc ( u W f λ ) . (S11)Then we arrive at the same result as Eq. S3. The Hamiltonian is defined as H = − J · I = − J · X ij σ ′ iz σ ′ j z (S12)where J is a constant and I is the intensity on the center of the detection plane ( u, v ) = (0 , III. EVOLUTION OF THE ARTIFICIAL SPIN GLASS SYSTEM TOWARDS GROUND STATE AT T=0
Next we show how the optical setup operates to search for the ground states of spin glass systems. We perform theoptical setup at zero effective temperature T = 0 for modeling spin glass systems in the ground state. The numberof spins is N = 100. At T = 0, the interacting system definitely evolves to its ground state during the Monte Carloiterations. The operation starts from a randomly initiated gauge-transformed spin configuration S ′ = (cid:8) σ ′ iz (cid:9) , encodedon the SLM with an array of 10-by-10 macropixels. The iteration process is presented in Fig. S3. We can see that thesystems evolve towards lower energies, and within 1000 iterations the equilibrium distribution is achieved, where allspin configurations in the statistical ensemble reside in the ground state. This is verified by the fact that the fields onSLM and CCD remain invariant after 1000 iterations. The detected images of the intensity I ( u, v ) of initial and finalstates are shown as the insets of Fig. S3 . When the ground state is arrived, the detected image becomes a brightspot, I ( u, v ) = δ ( u, v ), and the gauge-transformed spin configuration S ′ = (cid:8) σ ′ iz (cid:9) is uniform on SLM, which is exactlythe solution of minimizing the Hamiltonian H in Eq.(3). Here we use the normalized Hamiltonian by normalizing thedetected intensity in a finite region as H = − RR A [ I ( u, v ) × ¯ I ( u, v )]d u d v/ RR A ¯ I ( u, v )d u d v , as stated in the main text.Here ¯ I ( u, v ) is the field intensity on the detection plane when the SLM has uniform phase modulation of ϕ j, SLM =0in the square detection region A : | u | , | v | ≤ d/ d = fλ , which corresponds to the area for only the zerothdiffraction order.3 H | m |
500 100000 500 100000.20.40. .
810 500 100000.20.40.60.810 500 100000.20.40.60.810 500 100000.20.40.60.81 iterations initial I ( u,v ) final I ( u,v ) -4-8 ×10 ×10 FIG. S3. Experimental results of spatial optical spin glass system in the ground state. Evolution of the normalized Hamiltonian H and the gauge-transformed magnetization | m ′ | = (cid:12)(cid:12)(cid:12)(cid:12) N P i σ ′ iz (cid:12)(cid:12)(cid:12)(cid:12) during the Monte Carlo iterations for searching the ground state.Insets show the detected images on CCD, I ( u, v ), of initial and final states. The white bar corresponds to the length of 100 µ m. IV. THE SOLUTIONS FOR THE NUMBER-PARTITION PROBLEM AS PRESENTED IN FIG. 3 a be f cg dh FIG. S4. Solutions of the number partition problem for a number set containing 40000 elements. (a-d) One of the final statesof the gauge-transformed spin distribution (cid:8) σ ′ jz (cid:9) on SLM. (e-h) The resulted partition solutions { σ j } corresponding to the(a-d), respectively. V. ESTIMATION OF THE IMPACT OF THE INTEGRATION AREA SIZE ON THE DETECTIONPLANE