## All-optical scalable spatial coherent Ising machine
Marcello Calvanese Strinati, 1, ∗ Davide Pierangeli, 1, 2 and Claudio Conti 1, 2, 3
1 Physics Department, Sapienza University of Rome, 00185 Rome, Italy
2 Institute for Complex Systems, National Research Council (ISC-CNR), 00185 Rome, Italy
3 Centro Ricerche Enrico Fermi (CREF), Via Panisperna 89a, 00184 Rome, Italy
(Dated: November 15, 2021)
Networks of optical oscillators simulating coupled Ising spins have been recently proposed as a heuristic platform to solve hard optimization problems. These networks, called coherent Ising machines (CIMs), exploit the fact that the collective nonlinear dynamics of coupled oscillators can drive the system close to the global minimum of the classical Ising Hamiltonian, encoded in the coupling matrix of the network. To date, realizations of large-scale CIMs have been demonstrated using hybrid optical-electronic setups, where optical oscillators simulating different spins are subject to electronic feedback mechanisms emulating their mutual interaction. While the optical evolution ensures an ultrafast computation, the electronic coupling represents a bottleneck that causes the computational time to severely depend on the system size. Here, we propose an all-optical scalable CIM with fully-programmable coupling. Our setup consists of an optical parametric amplifier with a spatial light modulator (SLM) within the parametric cavity. The spin variables are encoded in the binary phases of the optical wavefront of the signal beam at different spatial points, defined by the pixels of the SLM. We first discuss how different coupling topologies can be achieved by different configurations of the SLM, and then benchmark our setup with a numerical simulation that mimics the dynamics of the proposed machine. In our proposal, both the spin dynamics and the coupling are fully performed in parallel, paving the way towards the realization of size-independent ultrafast optical hardware for large-scale computation purposes.
## I. INTRODUCTION
Solving large-scale optimization problems is extremely useful to several different fields of modern science, with applications ranging from biology to finance and social science [1-5]. These problems often belong to the nondeterministic polynomial (NP-hard) computational complexity class [6]: Finding the optimal solution requires computational resources that scale exponentially with the size of the system, making these problems intractable using conventional computer architectures. A tremendous amount of interest has been recently attracted by the development of unconventional computational methods (heuristic solvers) to solve probabilistically, but efficiently, large-scale optimization problems. A key observation behind these heuristic methods is that optimization problems can be mapped onto specific classical Ising models efficiently [7], i.e., in a polynomial (P) time. Solving the specific optimization problem then translates into the NP-hard problem of finding the ground state (GS) of the corresponding Ising Hamiltonian [8].
In recent years, several physical systems have been demonstrated to evolve according to the classical Ising Hamiltonian, therefore providing valuable ad hoc platforms to solve the Ising model for large-scale optimization purposes. Remarkable examples include two-component Bose-Einstein condensates [9, 10], superconducting circuits [11], trapped ions [12, 13], digital computers [14-17], electrical oscillators [18], optoelectronical oscillators [19], polariton condensates [20-22], laser networks [23, 24],
∗ marcello.calvanesestrinati@gmail.com
and coupled optical parametric oscillators (OPOs) [2538], which are the focus of this work. These networks, called coherent Ising machines (CIMs) exploit the fact that, when driven above the oscillation threshold, a second-order phase transition takes place [39, 40]: In the long-time limit, the phase of each OPO takes values 0 or π with respect to the reference phase enforced by the pump. Because of the bistable nature of its phase, an OPO is suitable to simulate a classical Ising spin, and systems of coupled OPOs in proper conditions can simulate the dynamics of coupled Ising spins [29, 41, 42].
Nowadays, major issues in realizing CIMs for realistic optimization problems involve, on one hand, the physical conditions (e.g., temperature) in which the machine has to operate, and on the other hand, the scalability and the connectivity that these systems can implement. In this respect, photonic systems offer a versatile platform to realize large-scale CIMs with general connectivity, while working at room temperature and being constructed from off-the-shelf components. An implementation of an alloptical CIM with few spins using time-multiplexed OPOs has been reported in Refs. [26, 27], where different OPOs are different temporal pulses within a nonlinear cavity, and optical delay lines are used to couple different OPOs. This approach allows the realization of arbitrary coupling topology, but cannot be scaled up to a large number of spins. An all-optical CIM with a largenumber of spins was reported in Ref. [28], implementing the one-dimensional nearest-neighbor coupling via a Mach-Zehnder interferometer. While this other approach allows the implementation of several spins, the coupling topology was limited to nearest-neighbor coupling. Large-scale CIMs with arbitrary coupling topology
FIG. 1. (a) All-optical spatial CIM. The setup consists of the following: (i) χ (2) NLM pumped by a laser at wavelength λ p non-resonant with the cavity mirrors (M), which amplifies the signal at λ s = 2 λ p , (ii) Coupling mechanism, (iii) Extraction of the signal by a beam splitter (BS) with reflection and transmission coefficients R out and T out = √ 1 -R 2 out , respectively, and detection (D). The coupling is implemented using different configurations of the SLM depending on the implemented graph. (b) Circulant graphs: A first lens (L1) transforms the signal in Fourier space, the SLM multiplies the FT of the field by ˜ Q k , and a second lens (L2) transforms the modulated field back to real space. (c) General graphs, based on the vector-matrix multiplication scheme. The SLM works in real space, multiplying the incoming field by Q ij . The modulated signal is focused along the y -axis by a cylindrical lens (CL1), defocused, and rotated by 90 ◦ on the xy -plane by another cylindrical lens (CL2) with 45 ◦ rotated optical axis. The pixels in panels (b) and (c) depict the OPOs encoded on the xy -plane: An amplitude A j,τ is defined as a single pixel in panel (b) , whereas as a column of M y pixels in panel (c) .
<details>
<summary>Image 1 Details</summary>
### Visual Description
## Optical System Diagram: Experimental Setup
### Overview
The image depicts a three-part technical diagram of an optical system. Part (a) shows a schematic layout with labeled components and light paths, while parts (b) and (c) provide detailed views of specific subsystems with spatial light modulators (SLMs) and coupling elements.
### Components/Axes
**Part (a):**
- **Labels**: Pump, D, BS, Coupling, M (repeated), χ² NLM
- **Axes**: x (rightward), y (upward), z (out of page)
- **Light Paths**:
- Green line (Pump beam)
- Red lines (signal paths)
- **Key Components**:
- BS: Beam Splitter
- NLM: Nonlinear Mixer (χ²)
- M: Mirrors (three instances)
- D: Detector
**Parts (b) and (c):**
- **SLM**: Spatial Light Modulator (labeled with ~Q_k and Q_ij heatmaps)
- **L1/L2**: Lenses
- **CL1/CL2**: Coupling Lenses
- **Heatmaps**: Color-coded intensity/phase patterns (no numerical values visible)
### Detailed Analysis
**Part (a):**
- Pump beam (green) enters from left, interacts with Beam Splitter (BS), then splits into multiple paths.
- Red paths represent signal beams interacting with:
- Nonlinear Mixer (NLM) for χ² frequency conversion
- Mirrors (M) for beam steering
- Coupling element for beam combination
- Detector (D) positioned at lower-left for signal collection.
**Parts (b) and (c):**
- SLMs modulate light spatially:
- (b) Shows ~Q_k pattern (checkerboard-like)
- (c) Shows Q_ij pattern (gradient stripes)
- Lenses (L1/L2) focus beams before/after SLM
- Coupling elements (CL1/CL2) combine beams with different modulation patterns
### Key Observations
1. **Symmetry**: Mirror configurations (M) in (a) suggest balanced beam steering.
2. **Modulation Diversity**: SLM patterns differ between (b) and (c), implying distinct experimental conditions.
3. **Color Coding**: Green = pump, Red = signal paths, Heatmaps = modulation states.
### Interpretation
This system demonstrates:
1. **Nonlinear Optical Processing**: The χ² NLM enables frequency mixing, critical for applications like parametric amplification.
2. **Adaptive Optics**: SLMs with programmable patterns (Q_k, Q_ij) suggest dynamic control over beam properties.
3. **Beam Management**: Mirrors and couplers optimize path alignment and power distribution.
The absence of numerical data in heatmaps suggests this is a conceptual schematic rather than experimental results. The system appears designed for investigating nonlinear light-matter interactions with programmable spatial modulation.
</details>
have been demonstrated using time-multiplexed OPOs in hybrid optical-electronic systems [32, 43], or optoelectronic oscillators [19]: Optical signals evolve in time subject to electronic measure and feedback mechanisms to emulate the spin-spin interaction. The optical nature of the setup ensures ultrafast computation. However, the presence of the electronic feedback inherently represents a 'bottleneck' that introduces an additional computational time that scales quadratically with the system size [44]. The realization of a scalable fully optical CIM without the hindrance of electronic components is thus highly desirable.
A step towards the realization of an all-optical scalable CIM has been recently made by exploiting spatial degrees of freedom of light [34, 45]. This approach relies on the usage of spatial light modulators (SLMs) to encode a spin variable σ j into the binary phase modulation of the optical field shining the j -th pixel of the SLM. An electronic feedback mechanism was also present. To date, a proposal of an all-optical scalable CIM implementing an arbitrary coupling topology is still missing. In this paper, we propose and theoretically validate an all-optical scalable spatial CIM with fully-programmable coupling. We first propose two different configurations of the SLM to realize two different classes of coupling, and then estimate the computational performance of our machine by means of a numerical simulation that closely mimics the temporal dynamics of coupled OPOs within a parametric cavity. We find that our proposed setup converges close to the minimum of the Ising Hamiltonian after a computational time that is orders of magnitude smaller compared to existing hybrid electro-optical setups.
This paper is organized as follows. In Sec. II, we discuss the scheme of our proposed setup. In Sec. III, we show our numerical results, and draw our conclusions in Sec. IV.
## II. SCHEME OF THE PROPOSED SETUP
The scheme of our machine is shown in Fig. 1 a . We consider an optical parametric cavity with a second-order nonlinear medium ( χ (2) NLM) of length L , shone by a pump laser (green beam) at wavelength λ p . We take z as the direction of propagation of light. The pump beam within the NLM parametrically amplifies degenerate pairs of signal and idler photons (red beam) at wavelength λ s = 2 λ p via spontaneous parametric down conversion. The spatial configuration of the signal wavefront on the xy -plane encodes the amplitudes { A j,τ } of the N OPO fields ( j = 1 , . . . , N ) at round trip number τ . These amplitudes are in general complex, i.e., the phase of the optical field on the xy -plane can take any value. However, the presence of the χ (2) NLM forces the optical phases to be either 0 or π with respect to the pump phase, making the amplitudes { A j,τ } effectively real, as detailed below.
Different configurations of the SLM with M x × M y pixels realize different couplings between different OPOs, as shown in panels (b) and (c) . In panel (b) , the SLM is placed at the focal plane of a first lens (L1). The discretization of the field in real (and thus momentum) space is enforced by the SLM: Different pixels define different OPO amplitudes, thereby defining N = M x × M y different OPOs. Since in commercial SLMs one typically has M x , M y ∼ 10 3 , this scheme allows to define N ∼ 10 6 OPOs. The SLM acts as a programmable matrix with a transmission function ˜ Q k , which multiplies at each round trip the Fourier transform (FT) of the incoming OPO amplitudes ˜ A k,τ = FT[ { A j,τ } ]. A second lens (L2) with same focal length as L1 transforms the modulated fields ˜ Q k ˜ A k,τ back to real space, yielding the inverse Fourier transform (IFT) of ˜ Q k ˜ A k,τ . By convolution theorem, the resulting field on each pixel after the FT-SLM-IFT sequence is A ′ j,τ = ∑ i Q i -j +1 A i,τ , where Q j is the IFT of ˜ Q k . The output takes the form of a cou-
pled field A ′ j,τ = ∑ i Q ij A i,τ , where the coupling matrix Q has entries Q ij ≡ Q i -j +1 . This matrix represent a rotationally invariant (or circulant) graph [46], where all nodes are equivalent. Notable examples are the nearestneighbor Ising chain and the M¨ obious ladder [26-28, 33].
While the coupling scheme in panel (b) gives the possibility to encode a large number of OPOs, granting a straightforward experimental implementation, it allows the implementation of a limited class of graph. To overcome this issue, we propose in panel (c) a different scheme for a general coupling matrix Q . This setup is based on the vector-matrix multiplication scheme [47]: The different OPOs are arranged on the xy -plane as different column vectors, such that the signal wavefront shines all M y pixels on a given column of the SLM with a uniform field. The SLM multiples in real space the vectorized signal with amplitude A j,τ by Q ij , such that the amplitude at the point ( i, j ) of the field after the SLM is Q ij A j,τ . Acylindrical lens (CL1) focuses the signal wavefront onto a single column along the y -axis, whose amplitude at point i is given by A ′ i,τ = ∑ j Q ij A j,τ . Propagation in free space defocuses the signal on the xy -plane, obtaining a vectorized signal arranged as a row matrix. Subsequent rotation by 90 ◦ of the field by a second cylindrical lens (CL2) recovers the structure of the signal as column vectors. Now, each column encodes the amplitude A ′ j,τ = ∑ i Q ij A i,τ , i.e., a coupled field with general Q ij . As such, this scheme implements any coupling matrix, but with the drawback that the OPOs need to be redundantly defined over M y pixels, limiting to N = M x the number of OPOs in the system.
We stress that the two coupling schemes presented here are fully optical and process all interactions in parallel, without need of electronic feedbacks. Since in our scheme the propagation of the OPOs within the cavity also occurs in parallel, our scheme realizes a size-independent large-scale spatial CIM [44], with critical advantages in terms of scaling and computational time compared to the existing hybrid optical-electronic devices.
In our setup, the binary nature of the phase of each OPO is enforced by the χ (2) NLM in Fig. 1 a . We inject into the NLM a pump field with spatially uniform wavefront that, at each pixel j on the xy -plane, mixes inside the NLM with the signal field. The subsequent dynamics along the z -axis, independently at each point, follows the second-order nonlinear wave equation [48] for the degenerate signal field { A j,τ } and pump { B j,τ } :
<!-- formula-not-decoded -->
where the star denotes complex conjugation. Here, κ = 2 πχ (2) /λ s n 2 , where χ (2) and n = n ( λ s ) are the nonlinear coefficient and the index of refraction of the NLM, respectively, and λ s = 2 π/k s . We assume perfect phase matching 2 k s = k p , where k s ( p ) is the wave vector of the signal (pump), which is achieved by tuning the temperature of the NLM to equalize the index of refraction for the signal and pump. To show phase-dependent amplification, we
FIG. 2. Histogram of time evolution of the (a) real and (b) imaginary part of the N = 112 OPO amplitudes (in units of A 0 , see text), computed using Eq. (3), for short times. As evident, the cavity dynamics enhances the real part of the OPO amplitudes and suppresses their imaginary part.
<details>
<summary>Image 2 Details</summary>
### Visual Description
## Line Plots: Real and Imaginary Components of A_j,τ/A_0 vs. τ
### Overview
The image contains two line plots (panels a and b) showing the real (Re) and imaginary (Im) components of a complex quantity A_j,τ normalized by A_0, plotted against a parameter τ. Both panels use color gradients to encode an additional parameter (likely j), with distinct scales and τ ranges.
---
### Components/Axes
- **Panel (a): Real Component**
- **X-axis (τ)**: Ranges from 0 to 200 (linear scale).
- **Y-axis (Re[A_j,τ]/A_0)**: Scaled by 10⁻³, ranging from -4×10⁻³ to +4×10⁻³.
- **Color Legend**: Vertical bar on the right, ranging from 1 (dark blue) to 10 (yellow), indicating parameter j.
- **Lines**: 10 distinct curves (one per j value), colored according to the legend.
- **Panel (b): Imaginary Component**
- **X-axis (τ)**: Ranges from 0 to 40 (linear scale).
- **Y-axis (Im[A_j,τ]/A_0)**: Scaled by 10⁻⁴, ranging from -5×10⁻⁴ to +5×10⁻⁴.
- **Color Legend**: Vertical bar on the right, ranging from 1 (dark blue) to 100 (yellow), indicating parameter j.
- **Lines**: 10 distinct curves (one per j value), colored according to the legend.
---
### Detailed Analysis
#### Panel (a): Real Component
- **Trends**:
- At τ ≈ 0, all lines start near 0 and are tightly clustered.
- As τ increases, lines diverge vertically, indicating growing variability in Re[A_j,τ]/A_0.
- By τ ≈ 200, the spread reaches ±4×10⁻³, with higher j values (yellow lines) showing larger amplitudes.
- **Key Data Points**:
- For j = 1 (dark blue), Re[A_j,τ]/A_0 ≈ 0.001 at τ = 100.
- For j = 10 (yellow), Re[A_j,τ]/A_0 ≈ 0.004 at τ = 200.
#### Panel (b): Imaginary Component
- **Trends**:
- At τ ≈ 0, lines start near 0 but spread rapidly, reaching ±5×10⁻⁴ by τ = 10.
- As τ increases, lines converge toward 0, suggesting damping or stabilization.
- Higher j values (yellow lines) exhibit larger initial amplitudes but decay faster.
- **Key Data Points**:
- For j = 1 (dark blue), Im[A_j,τ]/A_0 ≈ -0.0002 at τ = 20.
- For j = 100 (yellow), Im[A_j,τ]/A_0 ≈ -0.0005 at τ = 10.
---
### Key Observations
1. **Divergence in Panel (a)**: The real component grows more variable with increasing τ, with higher j values amplifying this effect.
2. **Convergence in Panel (b)**: The imaginary component decays toward zero as τ increases, with higher j values showing faster damping.
3. **Scale Differences**: Panel (a) uses a smaller y-axis scale (10⁻³ vs. 10⁻⁴) but a narrower j range (1–10 vs. 1–100 in panel b).
4. **Color Consistency**: Lines in both panels match their respective color legends (e.g., dark blue = j=1, yellow = j=10/100).
---
### Interpretation
- **Dynamic Behavior**: The real component (panel a) suggests a system where variability or oscillations grow with τ, potentially indicating instability or resonance effects. The imaginary component (panel b) implies damping or energy dissipation, as higher j values decay faster.
- **Parameter Sensitivity**: The imaginary component is more sensitive to j (color scale up to 100), while the real component is dominated by lower j values (scale up to 10). This could reflect different physical or mathematical roles for j in the two components.
- **Temporal Scaling**: The shorter τ range in panel (b) (0–40 vs. 0–200) may indicate that the imaginary component’s behavior is transient, while the real component evolves over a longer timescale.
---
### Notes on Data Extraction
- All axis labels, legends, and scaling factors were transcribed directly from the image.
- Color-to-value mappings were cross-verified between legends and line colors.
- Spatial grounding confirmed legends are positioned on the right, with labels on axes and panels.
</details>
rewrite by omitting j and τ in Eq. (1) A ( z ) = u ( z ) e iφ ( z ) and B ( z ) = u p ( z ) e iφ p ( z ) , where u ( u p ) and φ ( φ p ) are the amplitude and phase of the signal (pump), respectively. This allows to separate the dynamics of the relative phase θ = φ p -2 φ and of the amplitudes u and u p in Eq. (1) as
<!-- formula-not-decoded -->
From Eq. (2), one can see that the evolution of θ has two fixed points (modulo 2 π ): θ = 0, i.e., φ -φ p / 2 = 0 , π and θ = π , i.e., φ -φ p / 2 = ± π/ 2, corresponding to two distinct regimes: (i) Parametric amplification, where energy is transferred from the pump to the field, and (ii) Up-conversion, where energy is converted from the signal to the pump. Focusing on the first case, which is our case of interest, θ flows towards θ = 0, fixing φ to be either 0 or π with respect to φ p / 2, thereby manifesting phase-dependent amplification. In terms of the original variables, taking φ p = 0 (real pump), the evolution along z amplifies the real part of the fields Re[ A j,τ ] and suppresses their imaginary parts Im[ A j,τ ].
## III. NUMERICAL RESULTS
We now discuss the experimental realizability of our setup in Fig. 1. We follow Refs. [41, 49] using realistic experimental values of the system parameters. Each OPO amplitude A j,τ evolves into A j,τ +1 after a round trip by undergoing (i) Parametric amplification within the NLM [Eq. (1)], (ii) Coupling, and (iii) Measurement and losses. We then capture the OPO dynamics by the following map:
<!-- formula-not-decoded -->
In Eq. (3), NLM[ A l,τ ] represents the l -th OPO amplitude at the exit of the NLM (i.e., for z = L ), which is computed by integrating Eq. (1) for z ∈ [0 , L ] with initial conditions A j,τ ( z = 0) = A j,τ -1 and B j,τ ( z = 0) = B 0 ,
<latexi sh
1\_b
64="3FdYy2w fS
HLp
0k
X
>A
B
c
VD
N
J
U
vq
g
Z
u
M
oQ
r
T
j
K
+
Wz
E
7n
/
R
C
I
m
G
O
P
<latexi sh
1\_b
64="HUIE5m9MoNr
2TBc0
+uD
yw
>A
V
LS
FJ
vqk
g
Z
Q
j
K
d
Wz
7n
/
R
X
Y8
C
p
f
G
O
P
FIG. 3. Ising energy from the OPO phases during the round trips for (ML1),(ML2) M¨ obius ladder, (K1),(K2) complete graph, (ER1),(ER2) Erd˝ os-R´ enyi graph, and (BA1),(BA2) Barab´ asi-Albert graph. Upper panels refer to N = 112, while lower panels to N = 224. Different lines refer to different initial conditions. Horizontal black dashed lines mark the GS energy found by diagonalizing J for ML, and the minimal energy found by Metropolis annealing for K, ER, and BA. The insets in the upper panels depict the different connectivities in the four cases (black dots are nodes, and red and green lines depict negative and positive edge weights, respectively), where N = 16 is used for illustration purposes only.
<details>
<summary>Image 3 Details</summary>
### Visual Description
## Line Graphs with Network Diagrams: E_Ising vs τ Across Models
### Overview
The image contains eight panels arranged in a 2x4 grid. Each panel includes:
1. A line graph showing E_Ising (y-axis) vs τ (x-axis) on a logarithmic scale
2. A network diagram illustrating different graph structures
3. A legend with colored lines corresponding to data series
4. Model labels (ML1, K1, ER1, BA1, ML2, K2, ER2, BA2)
### Components/Axes
**Graph Elements:**
- **Y-axis (E_Ising):** Ranges from -60 to 0 in 5-unit increments
- **X-axis (τ):** Logarithmic scale from 10⁰ to 10⁴
- **Legend:** Positioned in top-right of each graph, contains 5-7 colored lines (red, blue, green, orange, purple, black)
- **Model Labels:** (ML1), (K1), (ER1), (BA1), (ML2), (K2), (ER2), (BA2) positioned above graphs
**Diagram Elements:**
- **Network Structures:**
- (ML1)/(K1): Wheel graph (central hub + peripheral nodes)
- (ER1)/(BA1): Random graph with varying connectivity
- (ML2)/(K2): Grid-like lattice structure
- (ER2)/(BA2): Scale-free network with power-law degree distribution
### Detailed Analysis
**Graph Trends:**
1. **(ML1):**
- Red line: Starts at -5, drops sharply to -30 by τ=10²
- Blue line: Begins at -10, declines gradually to -25 by τ=10³
- Green line: Flat at -15 until τ=10³, then drops to -30
2. **(K1):**
- Orange line: Sharp decline from 0 to -20 by τ=10¹
- Purple line: Gradual decrease from -5 to -25 by τ=10³
3. **(ER1):**
- Black line: Steady decline from 0 to -30 by τ=10²
- Red line: Slow decrease from -5 to -20 by τ=10³
4. **(BA1):**
- Blue line: Rapid drop from 0 to -25 by τ=10¹
- Green line: Moderate decline from -10 to -35 by τ=10³
5. **(ML2):**
- Orange line: Steep decline from 0 to -40 by τ=10²
- Purple line: Gradual decrease from -5 to -30 by τ=10³
6. **(K2):**
- Red line: Sharp drop from 0 to -35 by τ=10¹
- Blue line: Moderate decline from -10 to -35 by τ=10³
7. **(ER2):**
- Green line: Steady decrease from 0 to -40 by τ=10²
- Orange line: Gradual decline from -5 to -30 by τ=10³
8. **(BA2):**
- Purple line: Rapid drop from 0 to -45 by τ=10¹
- Black line: Moderate decline from -10 to -40 by τ=10³
### Key Observations
1. **Universal Decay Pattern:** All models show E_Ising decreasing with increasing τ
2. **Network Structure Correlation:**
- Scale-free networks (BA) exhibit fastest initial decay
- Grid structures (K) show intermediate decay rates
- Random networks (ER) demonstrate slower decay
3. **Model Variations:**
- ML models (ML1/ML2) show more gradual declines than kinetic models (K1/K2)
- BA2 (N=224) exhibits most pronounced decay compared to BA1 (N=112)
### Interpretation
The data suggests that network topology significantly influences the decay rate of E_Ising:
- Scale-free networks (BA) demonstrate fastest decay, likely due to hub nodes accelerating information propagation
- Grid structures (K) show intermediate behavior, consistent with regular connectivity patterns
- Random networks (ER) exhibit slowest decay, possibly due to less efficient information flow
- Machine learning models (ML) appear to optimize decay rates compared to traditional kinetic models (K)
The logarithmic τ scale indicates that decay processes follow power-law dynamics across all models. The consistent downward trend across all panels suggests a fundamental relationship between network structure and system stability in this Ising-like model.
</details>
where B 0 is the uniform pump amplitude at the entrance of the NLM, which provides the gain. There are two sources of loss: The SLM transmission function Q , and measurement and intrinsic loss encoded in R out . The balance between gain and losses during the initial round trips defines the value B 0 , th of the oscillation threshold: B 0 , th = -log[ R out ρ ( Q )] /κL , where ρ ( Q ) is the spectral radius of Q .
We shine the NLM with a pump laser at λ p = 532nm with spot radius w 0 = 100 µ m. The physical properties of the NLM are encoded in χ (2) , λ s , and n (determining κ ), and L . Here, we use χ (2) = 10 -11 m / V, λ s = 1064nm, n = 2 ( κ 1 . 48 × 10 -5 V -1 ), and L = 10cm, which defines the characteristic length scale in our simulations. To integrate Eq. (1), we use as characteristic field amplitude A 0 6 . 77 × 10 3 V / m, which yields the numerical rescaled nonlinear constant ˜ κ := κz 0 A 0 = 10 -2 , and thus integrate Eq. (1) for z ∈ [0 , 1] (in units of L ). At τ = 0, the signal consists of white noise with amplitude 10 -3 A 0 .
The SLM transmission function Q is written by including self-interaction and off-diagonal coupling terms J , which is a real symmetric matrix: Q = a 1 + b J , where 1 is the identity matrix. Since the SLM provides phase and amplitude modulation, Q ij is in general a complex number with | Q ij | < 1. The lossy nature of the SLM implies ρ ( Q ) < 1. To benchmark our proof-of-principle machine, we simulate N = 112 coupled OPOs and solve the MAX-CUT problem for four undirected graphs for which the optimization problem belongs to different classes of computational complexity (P and NP-hard) [41, 50]: The M¨ obius ladder (ML) [51], which is a circulant P-graph realized with the scheme in Fig. 1 b , and three NP-hard graphs, specifically the random Erd˝ os-R´ enyi (ER) and scale-free Barab´ asi-Albert (BA) graphs [52] with approximately 20% edge density, and the random complete (K) graph [53], realized with the scheme in Fig. 1 c . For
TABLE I. Values of the pump amplitude B 0 at the entrance of the NLM used in the simulation, for the data in Fig. 3.
| | ML | K | ER | BA |
|---------|----------------|-----------------|----------------|----------------|
| N = 112 | 1 . 2 B 0 , th | 1 . 33 B 0 , th | 1 . 3 B 0 , th | 1 . 3 B 0 , th |
| N = 224 | 1 . 2 B 0 , th | 1 . 36 B 0 , th | 1 . 3 B 0 , th | 1 . 3 B 0 , th |
the ML, the nonzero entries of J are J i,i +1 = α and J i,i + N/ 2 = α , with negative α . For the ER and BA graphs, J ij = 0 , ± β , where the three values are randomly chosen from the appropriate probability distribution to yield the chosen edge density [52], while for the K graph, J ij = ± γ , where the sign is randomly chosen with equal probability. We set a = 0 . 96, b = 0 . 04, α = -0 . 2, and β = 0 . 05, and γ = 0 . 03, for which ρ ( Q ) 0 . 98 in all cases. In this setup, by setting T out = √ 0 . 1, we obtain threshold pump powers P th ∼ 200 mW.
We first show in Fig. 2 the histogram of the time evolution of (a) Re[ A j,τ ], and (b) Im[ A j,τ ], for short times, specifically for ML. As evident, Re[ A j,τ ] is exponentially amplified, whereas Im[ A j,τ ] is instead suppressed. This confirms that the system is correctly in the phase-dependent amplification regime. Next, we simulate the OPO dynamics for the chosen graphs and compute the Ising energy from OPO phase configuration as E Ising ( τ ) = -(1 / 2) ∑ ij J ij σ i ( τ ) σ j ( τ ), where σ i ( τ ) = sgn( A i,τ ) and sgn( · ) is the sign function. We compare the energy for two different values of N , namely N = 112 and N = 224. The result is shown in Fig. 3. Different panels refer to different graphs and different N as in the labels, using the values of the pump amplitude B 0 at the entrance of the NLM given in Table I. Different solid lines refer to different initial conditions of the fields. The black horizontal lines mark the minimal Ising energy: The P
nature of the problem with the ML allows to find this value exactly by diagonizaling J , selecting the eigenvector with maximal eigenvalue [41]. Instead, for the NP-hard cases of K, ER, and BA, we estimate the minimal value by a Metropolis annealing algorithm. The fact that the steady-state value of E Ising coincides with the exact minimal value for ML, while sometimes it does not for K, ER, and BA, reflects the different computational complexity of the optimization [41, 50].
The key result in Fig. 3 is that our system approaches a steady state with energy close to the minimum energy of the corresponding Ising Hamiltonian after about 10 3 round trips, for both values of N . Since in our scheme all OPOs evolve in time in parallel, our device allows to use short cavities (typically D ∼ 1 m long) that in turn ensures short round trips times τ RT = 2 nD/c ∼ 10 ns, independent of N . We then estimate the total computational time as approximately 10 µ s. As such, the parallelization of the dynamics allows to envision orders of magnitude shorter computational time compared to existing CIM realizations [32, 33, 54, 55].
## IV. CONCLUSIONS
In conclusion, we propose a fully-optical scalable spatial CIM implementing different connectivities. The binary nature of the signal phase at different points on the wavefront is enforced by the NLM within the parametric cavity, and the spatial discretization of the wavefront is defined by the SLM within the cavity implementing the spin-spin optical coupling. The number of spins that our system encodes critically depends on the specific configuration of the SLM. To implement a fully-programmable
- [1] J. J. Hopfield, 'Neural networks and physical systems with emergent collective computational abilities,' PNAS 79 , 2554-2558 (1982).
- [2] M. Gilli, D. Maringer, and E. Schumann, Numerical Methods and Optimization in Finance (Elsevier Science, 2011).
- [3] Q. Zhang, D. Deng, W. Dai, J. Li, and X. Jin, 'Optimization of culture conditions for differentiation of melon based on artificial neural network and genetic algorithm,' Sci. Rep. 10 , 3524 (2020).
- [4] A. Degasperi, D. Fey, and B. N. Kholodenko, 'Performance of objective functions and optimisation procedures for parameter estimation in system biology models,' npj Syst. Biol. Appl. 3 , 20 (2017).
- [5] M. Ohzeki, S. Okada, M. Terabe, and S. Taguchi, 'Optimization of neural networks via finite-value quantum fluctuations,' Sci. Rep. 8 , 9950 (2018).
- [6] R. M. Karp, 'Reducibility among combinatorial problems,' in Complexity of Computer Computations (1972) pp. 85-103.
- [7] A. Lucas, 'Ising formulations of many NP problems,' Frontiers in Physics 2 , 5 (2014).
CIM, we propose a setup based on the vector-matrix multiplication scheme, where the SLM works in real space of the field. This scheme allows the implementation of any graph, however limiting the number of spins to N ∼ 10 3 due to the redundant encoding of the spin variables on the SLM pixels. We then propose an alternative coupling scheme with the SLM working in momentum space of the field, which allows the implementation of a limited class of graphs but it can host N ∼ 10 6 spins.
The all-optical nature of our machine presents a step towards the realization of large-scale scalable CIMs. First, in our proposal, both the OPO dynamics and their mutual coupling are fully parallelized, which makes the reach for the optimal solution size independent. Second, the parallel encoding of all OPOs allows for short cavity lengths, and thus drastically smaller computational time, compared to state-of-the art realizations with hybrid electronic-optical setups, where the OPOs are arranged as a temporal sequence of pulses and cavity lengths of approximately 1 km are needed. Another important advance of our proposal compared to existing realizations is that no measurement is performed during the time evolution to realize the mutual interaction. This makes our machine suitable to study fundamental features of coupled OPOs beyond optimization, like the emergence of robust macroscopic quantum entanglement [56, 57].
## ACKNOWLEDGEMENTS
We acknowledge funding from Sapienza Ricerca, PRIN PELM (20177PSCKT), QuantERA ERA-NET Co-fund (Grant No. 731473, Project QUOMPLEX), H2020 PhoQus Project (Grant No. 820392).
- [8] F. Barahona, 'On the computational complexity of Ising spin glass models,' J. Phys. A 15 , 3241-3253 (1982).
- [9] T. Byrnes, K. Yan, and Y. Yamamoto, 'Accelerated optimization problem search using Bose-Einstein condensation,' New J. Phys. 13 , 113025 (2011).
- [10] T. Byrnes, S. Koyama, K. Yan, and Y. Yamamoto, 'Neural networks using two-component Bose-Einstein condensates,' Sci. Rep. 3 , 2531 (2013).
- [11] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, 'Quantum annealing with manufactured spins,' Nature 437 , 194-198 (2011).
- [12] K. Kim, M.-S. Chang, S. Korenblit, R Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan, and C. Monroe, 'Quantum simulation of frustrated Ising spins with trapped ions,' Nature 465 , 590-593 (2010).
- [13] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and
J. J. Bollinger, 'Engineered two-dimensional Ising interactions in a trapped-ion quantum simulator with hundreds of spins,' Nature 484 , 489-492 (2012).
- [14] A. D. King, W. Bernoudy, J. King, A. J Berkley, and T. Lanting, 'Emulating the coherent Ising machine with a mean-field algorithm,' arXiv:1806.08422 (2018).
- [15] E. S. Tiunov, A. E. Ulanov, and A. I. Lvovsky, 'Annealing by simulating the coherent Ising machine,' Opt. Express 27 , 10288-10295 (2019).
- [16] H. Goto, K. Tatsumura, and A. R. Dixon, 'Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems,' Sci. Adv. 5 , eaav2372 (2019).
- [17] K. Tatsumura, M. Yamasaki, and H. Goto, 'Scaling out Ising machines using a multi-chip architecture for simulated bifurcation,' Nat. Electron. 4 , 208-217 (2021).
- [18] J. Chou, S. Bramhavar, S. Ghosh, and W. Herzog, 'Analog coupled oscillator based weighted Ising machine,' Sci. Rep. 9 , 14786 (2019).
- [19] F. B¨ ohm, G. Verschaffelt, and G. Van der Sande, 'A poor man's coherent Ising machine based on optoelectronic feedback systems for solving optimization problems,' Nat. Commun. 10 , 3538 (2019).
- [20] N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. T¨ opfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis, 'Realizing the classical XY Hamiltonian in polariton simulators,' Nat. Mat. 16 , 1120-1126 (2017).
- [21] K. P. Kalinin and N. G. Berloff, 'Simulating Ising and n -state planar Potts models and external fields with nonequilibrium condensates,' Phys. Rev. Lett. 121 , 235302 (2018).
- [22] K. P. Kalinin and N. G. Berloff, 'Global optimization of spin Hamiltonians with gain-dissipative systems,' Sci. Rep. 8 , 17791 (2018).
- [23] S. Utsunomiya, K. Takata, and Y. Yamamoto, 'Mapping of Ising models onto injection-locked laser systems,' Opt. Express 19 , 18091-18108 (2011).
- [24] C. Tradonsky, I. Gershenzon, V. Pal, R. Chriki, A. A. Friesem, O. Raz, and N. Davidson, 'Rapid laser solver for the phase retrieval problem,' Sci. Adv. 5 , 10 (2019).
- [25] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, 'Coherent Ising machine based on degenerate optical parametric oscillators,' Phys. Rev. A 88 , 063853 (2013).
- [26] A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, 'Network of time-multiplexed optical parametric oscillators as a coherent Ising machine,' Nat. Photonics 8 , 937 (2014).
- [27] K. Takata, A. Marandi, R. Hamerly, Y. Haribara, D. Maruo, S. Tamate, H. Sakaguchi, S. Utsonomiya, and Y. Yamamoto, 'A 16-bit coherent Ising machine for onedimensional ring and cubic graph problems,' Sci. Rep. 6 , 34089 (2016).
- [28] T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, 'Large-scale Ising spin network based on degenerate optical parametric oscillators,' Nat. Photonics 10 , 415-419 (2016).
- [29] R. Hamerly, K. Inaba, T. Inagaki, H. Takesue, Y. Yamamoto, and H. Mabuchi, 'Topological defect formation in 1D and 2D spin chains realized by network of optical parametric oscillators,' Int. J. Mod. Phys. B 30 , 1630014 (2016).
- [30] W. R. Clements, J. J. Renema, Y. H. Wen, H. M. Chrzanowski, W. S. Kolthammer, and I. A. Walms-
ley, 'Gaussian optical Ising machines,' Phys. Rev. A 96 , 043850 (2017).
- [31] T. Wang and J. Roychowdhury, 'Oscillator-based Ising machine,' arXiv:1709.08102 (2017).
- [32] Takahiro Inagaki, Yoshitaka Haribara, Koji Igarashi, Tomohiro Sonobe, Shuhei Tamate, Toshimori Honjo, Alireza Marandi, Peter L. McMahon, Takeshi Umeki, Koji Enbutsu, Osamu Tadanaga, Hirokazu Takenouchi, Kazuyuki Aihara, Ken-ichi Kawarabayashi, Kyo Inoue, Shoko Utsunomiya, and Hiroki Takesue, 'A coherent Ising machine for 2000-node optimization problems,' Science 354 , 603-606 (2016).
- [33] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, K. Enbutsu, T. Umeki, R. Kasahara, S. Utsunomiya, S. Kako, K. Kawarabayashi, R. L. Byer, M. M. Fejer, H. Mabuchi, D. Englund, E. Rieffel, H. Takesue, and Y. Yamamoto, 'Experimental investigation of performance differences between coherent Ising machines and a quantum annealer,' Sci. Adv. 5 , eaau0823 (2019).
- [34] D. Pierangeli, G. Marcucci, and C. Conti, 'Large-scale photonic Ising machine by spatial light modulation,' Phys. Rev. Lett. 122 , 213902 (2019).
- [35] L. Bello, M. Calvanese Strinati, E. G. Dalla Torre, and A. Pe'er, 'Persistent coherent beating in coupled parametric oscillators,' Phys. Rev. Lett. 123 , 083901 (2019).
- [36] T. Wang and J. Roychowdhury, 'Oim: Oscillator-based Ising machines for solving combinatorial optimisation problems,' in Unconventional Computation and Natural Computation (Springer International Publishing, Cham, 2019) pp. 232-256.
- [37] Y. Okawachi, M. Yu, J. K. Jang, X. Ji, Y. Zhao, B. Y. Kim, M. Lipson, and A. L. Gaeta, 'Demonstration of chip-based coupled degenerate optical parametric oscillators for realizing a nanophotonic spin-glass,' Nat. Commun. 11 , 4119 (2020).
- [38] D. Pierangeli, G. Marcucci, and C. Conti, 'Adiabatic evolution on a spatial-photonic Ising machine,' Optica 7 , 1535-1543 (2020).
- [39] E. Goto, 'The parametron, a digital computing element which utilizes parametric oscillation,' Proc. IRE 47 , 1304 (1959).
- [40] J. W. F. Woo and R. Landauer, 'Fluctuations in a parametrically excited subharmonic oscillator,' IEEE J. Quantum Electron QE-7 , 435 (1971).
- [41] M. Calvanese Strinati, L. Bello, E. G. Dalla Torre, and A. Pe'er, 'Can nonlinear parametric oscillators solve random Ising models?' Phys. Rev. Lett. 126 , 143901 (2021).
- [42] M. Erementchouk, A. Shukla, and P. Mazumder, 'Computational capabilities of nonlinear oscillator networks,' arXiv:2105.07591 (2021).
- [43] Yoshitaka Haribara, Shoko Utsunomiya, and Yoshihisa Yamamoto, 'A coherent Ising machine for max-cut problems: Performance evaluation against semidefinite programming and simulated annealing,' in Principles and Methods of Quantum Information Technologies , edited by Yoshihisa Yamamoto and Kouichi Semba (Springer Japan, Tokyo, 2016) pp. 251-262.
- [44] D. Pierangeli, M. Rafayelyan, C. Conti, and S. Gigan, 'Scalable spin-glass optical simulator,' Phys. Rev. Applied 15 , 034087 (2021).
- [45] S. Kumar, H. Zhang, and Y.-P. Huang, 'Large-scale Ising emulation with four body interaction and all-to-all connections,' Comm. Phys. 3 , 108 (2020).
- [46] P. J. Davis, Circulant Matrices , Chelsea Publishing Series (Chelsea, 1994).
- [47] J. Spall, X. Guo, T. D. Barrett, and A. I. Lvovsky, 'Fully reconfigurable coherent optical vector-matrix multiplication,' Opt. Lett. 45 , 5752-5755 (2020).
- [48] R.W. Boyd, Nonlinear Optics (Elsevier Science, 2008).
- [49] M. Calvanese Strinati, I. Aharonovich, S. Ben-Ami, E. G. Dalla Torre, L. Bello, and A. Pe'er, 'Coherent dynamics in frustrated coupled parametric oscillators,' New J. Phys. 22 , 085005 (2020).
- [50] K. P. Kalinin and N. G. Berloff, 'Complexity continuum within Ising formulation of NP problems,' arXiv:2008.00466 (2020).
- [51] R. K. Guy and F. Harary, 'On the M¨ obius ladders,' Canad. Math. Bull. 10 , 493-496 (1967).
- [52] R. Albert and A.-L. Barab´ asi, 'Statistical mechanics of complex networks,' Rev. Mod. Phys. 74 , 47-97 (2002).
- [53] D. Gries and F. B. Schneider, A Logical Approach to Discrete Math (Springer-Verlag, 1993).
- [54] Y. Haribara, H. Ishikawa, S. Utsunomiya, K. Aihara, and Y. Yamamoto, 'Performance evaluation of coherent Ising machines against classical neural networks,' Quantum Sci. Technol 2 , 044002 (2017).
- [55] Y. Yamamoto, K. Aihara, T. Leleu, K. Kawarabayashi, S. Kako, M. Fejer, K. Inoue, and H. Takesue, 'Coherent Ising machines-optical neural networks operating at the quantum limit,' njp Quantum Information 3 , 49 (2017).
- [56] S. Kiesewetter and P. D. Drummond, 'Weighted phasespace simulations of feedback coherent Ising machines,' arXiv:2105.04190 (2021).
- [57] Z.-Y. Zhou, C. Gneiting, J. Q. You, and F. Nori, 'Generating and detecting entangled cat states in dissipatively coupled degenerate optical parametric oscillators,' Phys. Rev. A 104 , 013715 (2021).