The collective behavior of numerous natural and social systems—ranging from phase transitions in matter to neural network dynamics and financial market volatility—can be described by a universal spin model that originated in physical science to study magnetism.1 The Hamiltonian of a spin model is given by Hr=−∑ijJijSi·Sj, where Jij is the interaction strength between the i’th and the j’th spin. The spin S is a unit vector, and Jij is a function of spin-spin displacement (rij). The spatial profile of J(rij) plays a crucial role in determining the system’s physical properties. Short-range interactions regulate critical behavior and lead to distinct phase transitions,2 whereas long-range interactions are likely to induce correlated states for quantum entanglement. Particularly, tuning the interaction strength among several neighboring terms in J(rij) can give rise to exotic magnetic complexity,3^,4 replica symmetry breaking,5^,6 or unusual topological effects.7^–9 Therefore, exploring spin systems of a general J(rij) has been actively pursued both theoretically and experimentally, for the potential to unravel novel effects and to solve combinatorial optimization problems.10

Optical platforms offer abundant opportunities for simulating spin models via diverse light-matter interactions, such as nonlinear optical effects,11^,12 spontaneous parametric down-conversion,13^,14 lasing,15^–17 and exciton-polaritons.18^,19 A typical scheme of the optical spin model is supported by an array of resonant structures, such as coupled waveguides20^,21 or microcavities,18^,22 where the J(rij) is realized by the field overlapping between adjacent optical modes through the leakage of evanescent fields. This strategy stands for a group of experimental efforts that aims to construct spin interactions from real-space field overlapping. However, implementing a general J(rij) necessitates intricate structure design and fabrication, which in turn restricts the functionality and scalability of the spin system. In addition to real-space approaches, theoretical analysis of spin models can be significantly enriched by examining the momentum spectrum of J(rij).23 For instance, the critical behavior of a short-range interaction corresponds to a second momentum (k2) in the Fourier transform of J(rij).24 Experimentally, the Fourier transform serves as a powerful tool in various spatial light calculations, including image processing,25 spatial differentiation,26^–28 and convolution.29^,30

Lately, a spatial photonic Ising machine (SPIM) has been reported as a promising optical architecture for solving large-scale spin models.31 Leveraging the parallel processing of free space light, the optical circuit of the SPIM is well-suited for linear matrix computation with a low power consumption and high operation speed independent of the spin model scalability.32 The proposed SPIM utilizes the central spot of the momentum space corresponding to a specific all-to-all spin interaction. Based on the primitive version, several upgraded SPIMs have utilized gauge transformation33^–35 and multiplexing techniques36^–39 to generalize the spin interaction functions. In general, the methods stemming from the central momentum of light leave the entire momentum space unexplored. Because of that, the generalization of spin Hamiltonian requires extra computational resources, which can be achieved through additional optical measurements, an increased number of wavelengths, or a larger spatial area of the spatial light modulator (SLM) chip. The required resources rely on the complexity of spin Hamiltonian. For example, for nearest-neighboring spin interaction with full-rank Hamiltonians, the required additional resource for computing the full Hamiltonian grows linearly with the number of spins.

In this work, we introduce a momentum-space-modulated spin Hamiltonian, which allows us to implement arbitrary rij-dependent (displacement-dependent) interactions J(rij) in a simple spatial light platform. As examples, we perform two groups of experiments to show the ability of our optical simulator. The first experiment applies a J1-J2-J3 Ising spin model by tuning the interaction strength ratio among the nearest neighbor (NN), next-to-NN, and 3rd NN terms, demonstrating distinct magnetic ground states that occur in iron chalcogenides. The second experiment performs optical annealing for the XY spin model with NN interaction, exhibiting a Berezinskii-Kosterlitz-Thouless (BKT) dynamics that is governed by vortex proliferation. These experiments exhibit high calculation precision, with application potentials to solve general spin models using a simple spatial light architecture.

We establish a general real-and-momentum space correspondence of spin Hamiltonians (see Note 1 in the Supplementary Material for details) Hk=−∑ijJ(rij)Si·Sj=−∬V(k)I(k)dkxdky.(1)

Here, the sum is taken over all spin pairs, including self-interactions. The modulation function V(k)∼∑ijJ(rij)eik·rij corresponds to the Fourier transform of J(rij), and I(k) is the normalized momentum-space intensity of light. The unit vector S is mapped to the phase of the optical field E=exp(iφ) via the relationship of S=(cosφ,sinφ). It results in an XY model for φ∈(−π,π] and an Ising model for φ∈{0,π}. Equation (1) reveals a general correspondence between a two-body spin Hamiltonian and a V(k)-modulated diffraction of light. The insertion of V(k) in the momentum space is equivalent to applying a corresponding spin interaction function J(rij) for the real-space spins. Specifically, for V(k)=1, Eq. (1) reduces to the Plancherel theorem for zero spin interaction. For V(k)=δ(k), it maps to a specific all-to-all spin interaction. Therefore, in this work, different spin Hamiltonians, such as short-range interaction, long-range interaction, and arbitrary J(rij), can be feasibly implemented here by applying different V(k) in the momentum space.

The experiment scheme is illustrated in Fig. 1(d). A plane wave laser impinges onto a phase-only SLM, generating a square array of optical fields with uniform amplitude and arbitrary phases. The phase of light is mapped to the in-plane spin vector in the same way as we have mentioned before [Figs. 1(a) and 1(b)]. The lens transforms the real-space light into the Fourier spectrum I(k), which is captured by a camera (see Note 2 and Fig. S1 in the Supplementary Material). The spin Hamiltonian Hk is experimentally obtained by measuring the V(k)-weighted light intensity. Alongside the optically computed Hamiltonian, the feedback is assisted by the computer system, which generates random phases via a Markov chain,40 and accepts new phase configurations by comparing the successively observed Hk according to the Metropolis algorithm.41 This iterative process ensures that light evolves towards smaller Hk.

The J1-J2-J3 model plays a pivotal role in revealing the intricate magnetism complexity of iron chalcogenides.3^,42^,43 The Hamiltonian of J1-J2-J3 model is given by Hr=−∑NNJ1Si·Sj−∑2NNJ2Si·Sj−∑3NNJ3Si·Sj. Here, the first sum is taken over all NN spin pairs, the second sum over next-to-NN spin pairs, and so on. For simplicity, we set J1=−1, and we define two parameters to describe the interaction ratios, R1=J2/J1, and R2=J3/J1, akin to the notation in Ref. 4. The theoretically calculated phase diagram for the J1-J2-J3 model at T=0 is depicted in Fig. 2(a).4 Here, T is the temperature normalized by J1/kB. There are four types of states separated by linear boundaries determined by the interaction ratios R1 and R2. Among them, the antiferromagnetic ground state is the most common one that can be achieved without considering J2 and J3. As R1 and R2 increase, distinct phases emerge, such as the double-stripe pattern (4×4) and staggered-dimer pattern (4×2) [Fig. 2(b)]. These patterns are important ground-state candidates for the iron chalcogenides.

We apply the J1-J2-J3 model to our simulator to solve the ground states, starting with the case of (R1,R2)=(0.5,0.9). The interaction function J(rij) and the Fourier transform V(k) are calculated and shown in Fig. 3(a), without applying periodic boundary conditions. A 10×10 random spin configuration is initially encoded onto the SLM, as depicted in Fig. 3(b). Meanwhile, the diffraction of light in the momentum space is captured by the camera. Because of the binary phase encoding for the Ising model, the diffraction in the momentum space is centrosymmetric with I(k)=I(−k). To reach the ground state without being trapped in a local minimum, we have divided the annealing process by 10 temperatures from T=1.8 to T=0 [Fig. 3(e)]. At every temperature, we run 2000 iterative steps to update the phase configuration to reach thermal equilibrium (see Note 4 and Figs. S8 and S9 in the Supplementary Material). During this process, the experimentally obtained Hk is depicted in Fig. 3(e). For high-T cases, the Hk is overall decreasing as the iteration continues, but it also exhibits strong fluctuation due to the simulated thermal effect from the Metropolis algorithm. At around T∼1, a minimal Hk is reached in good agreement with the theoretically predicted critical temperature in Fig. 3(f). The annealing process is terminated after 2×104 iterative steps and a double-stripe pattern (4×4) is observed in Fig. 3(c), which agrees with the theoretical prediction in Fig. 2(b). The momentum space intensity I(k)sol.. The ground state is shown in the lower panel of Fig. 3(c), where we see four diffraction spots. During the experiment, a linear relationship between Hk and Hr is observed to verify Eq. (1), as depicted in Fig. 3(d).

Next, we systematically vary the interaction ratios R1 and R2 to uncover the complete phase diagram in Fig. 2. Particularly, the antiferromagnetic (AF) state [Fig. 4(a)], single-stripe [Fig. 4(d)], and double-stripe [Fig. 4(c)] are all in perfect agreement with the theoretical predictions. For (R1,R2)=(0.5,0.3), the solved state is slightly higher than the predicted ground state due to the competitive interaction of the staggered-stripe pattern [Fig. 4(b)]. The newly observed state for (R1,R2)=(0.5,0.9) is a mixture of two cases, which is different from that in Fig. 3(c). This phenomenon arises because of the ground state degeneracy because random flips during optical annealing can stochastically choose one spin distribution over another, leading to a different or mixture of ground state configurations. From the experimental results, we can see that as V(k) modulates the flow of light, it does not fix the spatial profile of the final I(k). Therefore, there is no target image as proposed in Ref. 31 because the Hamiltonian is defined by a real-valued number from Eq. (1), rather than an intensity distribution. More importantly, a single target image cannot map to all ground states of a spin system, especially for a frustrated spin system that has many degenerated modes.

We extend the use of our simulator to solve an XY model by extending the spins as quasi-continuous variables. A well-known phenomenon from the XY model is the BKT phase transition, which is a critical phenomenon characterized by the number of vortices.44 For simplicity, we consider a ferromagnetic NN interaction that is reduced from the J1-J2-J3 model by setting J1=1 and J2 =J3=0. Practically, an SLM can only generate discrete phases. Therefore, we utilize a q-state clock model45^,46 to approach the XY model by dividing a continuous 2π phase into q levels. In our experiments, we encode an array (20×20) of random phases onto the SLM with q=8. Typical phase distributions at different temperatures are shown in Fig. 5(a). It can be seen that the high-T system is strongly affected by thermal noise, and the phases change abruptly between neighboring spins. When the system is gradually cooling down, random phases tend to be collinear with each other by the NN spin interaction. The topological phase transition is characterized by the evolution of vortex number (Nv) as a function of T, as depicted in Fig. 5(b). The Nv crossover occurs around T=1, agreeing with the numerical calculation (critical temperature Tc∼0.9). Meanwhile, we capture Hr as a function of T, as presented in Fig. 5(c). The Hamiltonian of the XY system remains at a temperature of T≈0.5, indicating the optical noise level of q=8 model. In the q-clock model, the larger the value of q, the more susceptible the system is to noise. This is because, in a random process, a larger q corresponds to a smaller phase shift for each spin pixel (2π/q), resulting in a smaller variation in relative intensity in momentum space. Notably, this noise can be reduced by using a different photodetector.

At last, we perform a quenching experiment to showcase the out-of-equilibrium dynamics of the XY spin system. To do that, we start with a random phase distribution and set T=0. The quenching dynamics evolve through a coarsening process revealed in the Hk curve presented in Fig. 5(e)47 (see Note 6 and Figs. S11–S14 in the Supplementary Material). The final spin configuration is trapped in a typical vortex-pair state as shown in Fig. 5(d).

In summary, we have demonstrated an optical simulator to find the ground states of spin models with distinct displacement-dependent spin interactions by exploiting the momentum-space modulation of light. Although spin encoding relies on real-space optical phases, the spin interaction is governed by the modulated function V(k) in momentum space. Remarkably, our approach allows us to extract different Hamiltonians from a single optical measurement. Compared with previous SPIMs, which inevitably increase the computing steps for complex spin interactions or complicate the optical system, our system maintains a constant computing time consumption for different spin Hamiltonian from a single Fourier transform without sacrificing the scalability of primitive SPIM.

Currently, we have utilized linear optics to demonstrate two important physical phenomena, including the complex ground states from the J1-J2-J3 model, and the BKT transition from the XY model. It has been shown that a quadratic spin interaction can be realized using nonlinear optical effects,48^,49 which can be combined with our moment-space modulation method to implement a J1-J2-J3-K model3 to show the magnetism complexity of iron-based superconductors. Moreover, the implementation of an XY model is a critical step to generalize optical simulators with continuous spin variables adapting a wider group of combinatorial optimization problems. Leveraging the versatile multiplexing technique of optics, our system can be potentially utilized to realize asymmetric or nonreciprocal interactions, such as the Dzyaloshinskii–Moriya50^,51 and Haldane interactions to simulate nontrivial topological effects. For example, we can introduce asymmetric or random spin locations onto the SLMs. This allows for the implementation of asymmetric or arbitrary spin-spin interactions Jij, thereby broadening the functionality of our SPIM (Note 8 and Fig. S16 in the Supplementary Material). Moreover, the architecture of our spin model simulator is compatible with intracavity optical systems.52^–54 In this sense, V(k) can be performed by SLMs or high-performance metasurfaces55^–58 that are placed in the momentum plane of an intracavity system with coupled lasers or exciton-polaritons. As proof-of-demonstrations, we utilized a simple annealing algorithm in this work. In the near future, optimization algorithms adapting to the SPIMs59 and field programmable gate arrays60 can be applied to enhance the computation performance of our simulator for better application purposes.

Juan Feng is currently a PhD student at Shanghai Jiao Tong University (SJTU). She received her BS degree from Dalian University of Technology in 2019. Her current research interests include photonic Ising machines, spin-orbit interactions of light, and topological photonics.

Zengya Li is currently a PhD student at SJTU. He received his BS degree from the School of Physics and Microelectronics, Zhengzhou University, in 2020. His current research interests include nonlinear optics and nanophotonics.

Luqi Yuan is now a professor in the School of Physics and Astronomy at SJTU. He received his PhD in physics from Texas A&M University and was a postdoctoral scholar at Stanford University. With a focus on quantum optics and photonics, his research interests encompass a wide range of theoretical physics topics in topological photonics, non-Hermitian optics, synthetic dimensions, and hybrid quantum systems.

Erez Hasman is the Schlesinger chaired professor at the Technion – Israel Institute of Technology. He is a pioneer in the fields of optical metasurfaces, geometric-phase optical devices, and spin-optics. His current research interests include atomic-scale photonics, quantum materials, and quantum sensing and devices.

Bo Wang is an associate professor in the School of Physics and Astronomy at SJTU. He received his PhD from Peking University in 2017, and he did his 4-year postdoctoral research at the Technion – Israel Institute of Technology. His research interests include light manipulation at nanoscales, spin-orbit interactions of light, and photonic Ising machines.

Xianfeng Chen is a distinguished professor in the School of Physics and Astronomy at SJTU. He received his PhD in physics at SJTU in 1999. His current research focuses on nonlinear photonics, integrated photonics, biophotonics, and quantum information. In the past years, over 400 journal papers have been published in leading refereed journals and received more than 10,000 citations.