Large-scale Ising Emulation with Four-Body Interaction and All-to-All Connection

We propose and demonstrate a nonlinear optics approach to emulate Ising machines containing up to a million spins and with tailored two and four-body interactions with all-to-all connections. It uses a spatial light modulator to encode and control the spins in the form of the binary phase values of wavelets in coherent laser beams, and emulates the high-order interaction with frequency conversion in a nonlinear crystal at the Fourier plane. By adaptive feedback control, the system can be evolved into effective spin configurations that well approximate the ground states of Ising Hamiltonians with all-to-all connected many-body interactions. Our technique could serve as a new tool to probe complex, many-body physics and give rise to exciting applications in big data optimization, computing, and analytics.


Introduction
A wide range of modern applications across biology 1 , medicine 2 , finance 3 , and social networks 4 benefit from efficient processing and optimization of big data with complex structures and correlations. However, many such tasks are non-deterministic polynomial time hard (NP-hard), which could take existing supercomputers years to solve 5 . In this challenge, intense research efforts are underway to pave alternative approaches for computing and information processing. Among them, Ising machines have been shown to offer viable solutions for important NP-hard problems such as MAX-CUT 6 , protein folding 7 , and traveling salesman 8 , among others [9][10][11][12][13][14][15] . To this end, a variety of Ising machines have been demonstrated in effective spin systems of trapped atoms 16, 17 , polariton condensates 18 , superconducting circuits 19 , coupled oscillators 20-22 , nanophotonic waveguides [23][24][25] , randomly coupled lasers [26][27][28] , and time-multiplexed optical parametric oscillators 29,30 .
For the tasks of finding the ground states of many-body Hamiltonians, photonic systems enjoy the distinct advantages of high connectivity and speed [31][32][33][34][35][36] . For example, a fast coherent Ising machine can be realized in a looped optical parametric oscillator with temporally multiplexed pulses 37, 38 , albeit with limited spin numbers 37 or relying on photodetection and electronic feedback to emulate the spin-spin interaction 39, 40 . In contrast, a linear-optical Ising machine based on spatial light modulation was shown to subtend about 80,000 spins by coding them as the binary phases of pixels on a spatial light modulator (SLM) 41 . The far-field optical power of a modulated beam gives the expected energy of spin-spin interaction. The relatively simple setup yet high connectivity and scalability make this approach attractive to Ising machines with fully connected two-body interaction.
Yet, there are physical systems and numeric models whose dynamics cannot be fully captured by two-body interactions, and proper descriptions of multi-body interaction are required [42][43][44][45][46] . This poses a significant computational challenge, whose complexity and volume exceeds by far that of Ising problems with only two-body interaction, even for a moderate number of spins 47,48 . While a small class of many-body interaction can be decomposed onto a series of two-body interactions via some recursive or algebraic means [49][50][51] , they often subject to strict constraints 52,53 or require tedious error corrections 54,55 . For simulating complex systems and processing data with high-order correlation, suitable Ising machines remain desirable that support simultaneously high connectivity, multi-body interaction, and a large number of spins.
In this paper, we propose and experimentally demonstrate such an Ising machine hosting adjustable two-body interaction, four-body interaction, and all-to-all connections over a large number of spins. It realizes the spins as the binary phases of wavelets in a coherent laser beam, and implements effective multi-body interaction through nonlinear frequency conversion. Using SLM's (or equivalently, digital micromirror devices), one million spins are easily accessible. The fully connected two-body interaction is emulated with the optical power of the modulated light in the Fourier plane. The four-body interaction, also fully connected, is realized effectively by passing the modulated light through a lithium-niobate crystal in the Fourier plane for second harmonic (SH) generation. By simultaneously measuring the optical powers of the modulated light and its SH coupling into a fiber, complex Hamiltonians with all-to-all connected two-body and four-body interactions can be emulated over nearly one million spins. Through feedback control, the system can be evolved into the vicinity of the ground state of its Hamiltonian, exhibiting ferromagnetic, paramagnetic, and other novel nonlinear susceptibility phase transitions.
The present Ising emulator could pave a pathway to otherwise inaccessible territories of big data analytics and quantum simulation 45,46 . The high-order, many-body interaction can also serve as powerful activation functions for all optical machine learning 56,57 . For example, an immediate application is to use this machine as the q-state Potts model with two-body and four-body interactions on a square lattice 58 . Finally, while the current setup uses SH generation, even richer physics and controllability can be achieved by using other nonlinear optical processing like sum-frequency generation 59 , and four-wave mixing 60 , where other types of spin interaction and connection can be engineered.

Theoretical Analysis
The basic idea of the present Ising machine is illustrated in Fig. 1, which emulates chemical potential, two-body interaction, and four-body interaction over a large number of spins. Each spin is encoded as the binary phase of a pixel on a SLM. The total chemical potential energy is represented by the weighted sum of all spins. To realize the interactions, a coherent Gaussian pump beam is reflected off the SLM, focused using a Fourier lens to a nonlinear crystal for SH generation. The resulting beams at the original fundamental wavelength and the new SH wavelength are then separated at a dichroic mirror and captured by optical fibers. The pump power in the fiber is then measured to emulate the energy associated with spin-spin interaction, and that of the SH beam is to capture the four-body interaction among spins. Incorporating all three, a Hamiltonian describing the chemical potential, two-body interaction, and four-body interaction can be effectively constructed. To derive the effective Hamiltonian, we consider a Gaussian input pump beam of wavelength λ p , peak amplitude E 0 , and beam waist w p . It shines a SLM whose phase mask consists of pixels (m, n) centered around (x m , y n ), each giving 0 or π phase modulation. The transverse electric field immediately after the SLM is approximately 41 Here ξ mn = E 0 exp −(x 2 m + y 2 n )/w 2 p is the amplitude at pixel (m, n), Π is the rectangular function of width a, and σ mn = ±1 for the 0/π binary phase modulation.
The electric field is then transformed using a Fourier lens of focal length F and coupled into a periodic-poled lithium niobate (PPLN) crystal of length L, so that it reads at the center of the crystal (z = 0): Here, η mn = exp −2πi(xx m + yy n )/λ p F , κ p = (2πn p )/λ p , and n p is the refractive index of the pump in the PPLN crystal. For simplification, we introduce contracted notations ξ i and σ i , with ξ i=m+(n−1)N ≡ ξ mn , and σ i=m+(n−1)N ≡ σ mn , with i = 1, 2, ...N 2 to index the N × N spins (pixels). In our setup, only near-axis light is fiber coupled and measured, so that sinc(axπ/λ p F ), In the PPLN crystal, the pump creates its SH according to the following dynamics, which evolves from −L/2 to L/2. Here, κ h = (2πn h )/λ h is the wave number of the SH wave in the crystal with refractive index n h . ω p and ω h are the angular frequencies of the pump and SH waves, respectively. ∆κ = 2κ p − κ h − 2π/Λ is the phase mismatching, with Λ being the poling period. Equations (4)-(5) can in principle be solved by using split-step Fourier and adaptive step-size methods. However, the numeric solutions consume significant computational time and resources that increase exponentially with the spin number 61 , whose inefficiency calls for the present optical realization. Only under the conditions of phase matching (∆κ = 0), undepleted pump approximation, and negligibly small diffraction terms in Eqs. (4) and (5), the transverse electric-field of the output SH wave can be obtained analytically as At the crystal output, the pump and SH waves are each coupled into a single-mode fiber for detection, whose optical power is given by are the normalized back-propagated fiber modes of beam waist w p f and w h f for the pump and SH waves, respectively.
Substituting Eq. (3), (6) and (8) in Eq. (7), the detected power for the pump and SH waves is given in the form of and where J ij and J ijsr prescribe the strength of the two-body and four-body interactions, respectively.
As an example, Section 1 of the Supplementary Material presents the analytic results of J ij and J ijsr under the approximation of Eq. (6).
With Eq. (9) and (10), we can define a single parameter E to characterize the "energy" of the whole system, as where C = N 2 i=1 µ i σ i is the weighted sum of spins that represents their total chemical energy, with the local chemical potential µ i ∈ [−1, 1]. In this equation, α, β, and γ are free parameters defining the contribution of the chemical potential, two-body, and four-body interaction energy, respectively, to the total energy.
The total energy E can be minimized by optimizing the binary phase mask on the SLM through adaptive feedback control (see Fig. 4). This is equivalent to finding the ground-state solutions of the effective Hamiltonian whereĤ 1 is the Hamiltonian describing the chemical potential,Ĥ 2 for the two-body andĤ 4 for the four-body interaction, respectively, witĥ The two-body and four-body interactions can be tailored by modulating the input pump wave, varying the fiber optical modes, and modifying the phase matching conditions for the nonlinear  Fig. 2 (c) and (d), respectively. As shown, the many-body interactions can be tailored into complex forms by simple linear optics operations.

Experimental Setup
The experimental setup for the present nonlinear optical Ising machine is shown in Fig. 3. We use an optical pulse train at 1551.5 nm as the pump. Each pulse has 5 ps full-width at half-maximum 100D with sensors S132C and S130C). The measurement results are sent to a computer through a MATLAB interface for feedback control, which updates the phase mask on the SLM to find the ground state of the customized Ising Hamiltonian.

Results and discussions
As shown in Eq. (13), the chemical potential of each spin is flexibly defined by µ i and their collective contribution to the total energy is controlled by α. This provides the knob of studying the magnetization under a variety of local and global single-spin parameters. For this paper, however, we will focus on the many-body interaction and consider only µ i = 1 in all of the following results.
Meanwhile, we will leave fine tuning two and four-body interaction to our future work, but only control each's aggregated contribution to the total energy by varying β and γ, respectively.
To find the ground states of the total Hamiltonian, the SLM's initial phase mask is prepared in small clusters with randomly chosen 0 or π phases. The resulting pump and SH waves are coupled separately into single mode fibers, whose optical power is measured for feedback control. To minimize the total energy, we adaptively flip the spins within a randomly chosen cluster, following the standard Monte Carlo approach 63 . A flow chart of this procedure is shown in Fig. 4 Step 1: generate an initial random binary phase spin pattern on a SLM.
Step 2: define the total energy of the system. Step 3: define the range for s and τ and run t random trials for each.
Step    Figure 5 (a) shows how the optical pump power is increased, thus the decrease of total energy E to approach the ground state of the system. With α = 0, there is spontaneous symmetry breaking, as the system energy remains unchanged if all spins are flipped. As such, the feedback control will optimize the spins toward either positive or negative magnetic states with equal probability 41 . This is clear in Fig. 5 (b), where the magnetization trends both ways. For all initial phase trials, the absolute value of average magnetization can reached to ∼ 0.75.
In Fig. 6, we show the two sets of measurements for 400 × 400 and 800 × 800 spins, respectively, with purely four-body interaction, i.e., α = 0 and β = 0.  power evolution for γ = 1, for which the system ground states correspond to the minimum SH power. In both figures, different initial spins are optimized to give similar minimum SH power, which indicates the robustness of our optimization method. For γ = 1, the system evolves into a paramagnetic-like state that minimizes the SH power in the fiber. In Fig. 6(c), the values are close to the minimum detectable power of our optical sensor (∼5 nW). Because of a smaller pixel size in (c), the spin disorder is stronger to give lower SH power in (a). In opposite, (b) and (d) are for γ = −1, where the ground states are obtained at the highest SH power. In this case, the system exhibits a ferromagnetic-like behavior. For 400 × 400 and 800 × 800 spins, the optimization leads to similar maximum SH power despite different initial spin conditions. The convergence is slower for the latter case, as there are four times more spins to be optimized. Overall, our system can reliably and efficiently evolve into the vicinity of its ground state.
To further understand the optimization mechanism, we take images of the pump and SH beams at the crystal output by splitting them using flipping beam-splitters, as shown in Fig. 3.
Through 4f systems, the pump is imaged on a NIR-IR camera (FIND-R-SCOPE Model No. 85700 with pixel resolution of 17.6µm), and the SH on a CCD camera (Canon Rebel T6 with pixel pitch of 4.3µm). Figure 7 (a) and (b) show the pump and SH images, respectively, under different iteration numbers as they minimize the energy of 800×800 spin system with α = 0, β = 1, and γ = 1. As shown, both light beams become scattering and show speckle patterns as the optimization goes.
The resulting fiber-coupled SH and pump power is shown in Fig. 7 (c). Both decrease with the iteration numbers, dropping by orders of magnitude to minimum values after 150 iterations. Note that the final SH power level is very close to the detection level of the sensors, which prohibits further reduction via the present feedback control.
As seen in Eq. (12), the two-body and four-body interaction energies are dependent only on the relative alignment of the spins, but not on each's absolute orientation. Thus spontaneous symmetry breaking could occur during the optimization, leading to bifurcation 41 . To avoid this symmetry breaking, one could set a non-zero α to control the convergence direction of spin optimization. As an example, in Fig. 8, we compare the results with α = 1 and −1, both with β = −0.5 and γ = −1. As shown, α can indeed dictate the spin alignment to result in either positive or negative magnetization states. In both cases, starting from a rather randomized phase mask, the spins become relatively aligned to increase the SH and pump power, but magnetization can have positive or negative orientation as can be seen Fig. 8(b) and (d). The inset of Fig. 8(b,d) shows We last consider the cases where two-body and four-body interactions contribute oppositely to the total energy of the system. For instance, the two-body interaction can be attractive but four-body be repulsive, or vice versa. Such systems can be conveniently configured by defining the pre-factors β and γ in Eq. (12). The optimization will maximize one while minimizing the other. As an example, Fig. 9 plots the optimization trajectory for opposite two-body and four-body interactions. In Fig. 9(a), β = 1, and γ = −1, so that the system energy is maximized by reducing the pump power while increasing the SH power. Its opposite configuration is in Fig. 9(b), where

Conclusion
Using spatial phase modulation and second-harmonic generation, we have constructed Ising emulators with all-to-all connections and tailored chemical potential, two-body interaction, and fourbody interaction. Their ground-state solutions can be effectively and reliably approximated by adaptive feedback control, whose speed is currently limited by the processing time of the spatial light modulator. A significant speedup is achievable by using ferroelectric liquid-crystal based spatial light modulators or programmable plasmonic phase modulators 71,72 . At the present, the maximum number of accessible spins is about 1 million, and can be further increased by using modulators with more pixels or combining multiple modulators. While this study considers only second-harmonic generation, it can be straightforwardly extended to other nonlinear processes, such as sum-frequency generation, difference-frequency generation, four-wave mixing, for other interesting Ising machines [59][60][61][62]73 . Such nonlinear optical realizations of Ising machine could contribute as supplements for big data optimization and analyses that remain challenging for classical super-computers or forthcoming quantum machines but with a limited number of qubits.

Supplementary Materials SM1 Theoretical description for two-body and four-body interaction terms
To present the analytic results of J ij and J ijsr , Eq. (9) and (10) of the main text are described in this section. The detected pump power and SH power, each coupled into the single mode fiber, are given by is the electric field of the pump wave,  5) and the detected optical power of SH wave is ξ i ξ j ξ s ξ r ξ f is ξ f jr σ i σ j σ s σ r , (S.6) where ξ i={m+(N −1)n} = E 0 exp −(x 2 m + y 2 n )/w 2 p , (S.7) p F 2 (x m + x l ) 2 + (y n + y k ) 2 , (S.9) and A = (iω 2 h χ (2) L)/c 2 κ h . w p f and w h f are the beam waists of the normalized back-propagated fiber modes for pump and SH waves, respectively. It gives J ijsr σ i σ j σ s σ r , (S.11) where J ij = 2π(w p f ) 2 ξ i ξ j ξ f i ξ f j and J ijsr = 2π(w h f ) 2 A 2 ξ i ξ j ξ s ξ r ξ f is ξ f jr are the two-body, and fourbody interaction terms, respectively.
SM2 Results for the evolution of the pump and its SH with magnetization As an example, in Fig. S1, we compare the results with α = 1 and −1, both with β = −1 and