## Coherent Ising Machines with Optical Error Correction Circuits
Sam Reifenstein Satoshi Kako Farad Khoyratee Timoth´ ee Leleu Yoshihisa Yamamoto*
Sam Reifenstien, Satoshi Kako, Farad Khoyratee, Yoshihisa Yamamoto
PHI (Physics & Informatics) Laboratories, NTT Research Inc.
940 Stewart Drive, Sunnyvale, CA 94085, U.S.A.
Email Address:yoshihisa.yamamoto@ntt-research.com
Timoth´ ee Leleu
International Research Center for Neurointelligence, The University of Tokyo,
7-3-1 Hongo Bunkyo-ku, Tokyo 113-0033, JAPAN
Keywords: Coherent Ising machine, Chaotic solution search, Matrix-vector multiplication, Combinatorial optimization, Optical error correction
We propose a network of open-dissipative quantum oscillators with optical error correction circuits. In the proposed network, the squeezed/anti-squeezed vacuum states of the constituent optical parametric oscillators below the threshold establish quantum correlations through optical mutual coupling, while collective symmetry breaking is induced above the threshold as a decision-making process. This initial search process is followed by a chaotic solution search step facilitated by the optical error correction feedback. As an optical hardware technology, the proposed coherent Ising machine (CIM) has several unique features, such as programmable all-to-all Ising coupling in the optical domain, directional coupling ( J ij = J ji ) induced chaotic behavior, and low power operation at room temperature. We study the performance of the proposed CIMs and investigate how the performance scales with different problem sizes. The quantum theory of the proposed CIMs can be used as a heuristic algorithm and efficiently implemented on existing digital platforms. This particular algorithm is derived from the truncated Wigner stochastic differential equation. We find that the various CIMs discussed are effective at solving many problem types, however the optimal algorithm is different depending on the instance. We also find that the proposed optical implementations have the potential for low energy consumption when implemented optically on a thin film LiNbO 3 platform.
## 1 Introduction
Combinatorial optimization problems are ubiquitous in modern science, engineering, medicine, and business. Such problems are often NP-hard; hence, their runtime on classical digital computers is expected to scale exponentially. A representative example of NP-hard combinatorial optimization problems is the non-planar Ising model. [ 1 ] Special-purpose quantum hardware devices have been developed for finding solutions of Ising problems more efficiently compared tothan standard heuristic approaches. For example, a quantum annealing (QA) device exploits the adiabatic evolution of pure-state vectors using a timedependent Hamiltonian. [ 2 , 3 ] Another example is a coherent Ising machine (CIM), which exploits the quantumto-classical transition of mixed-state density operators in a quantum oscillator network. [ 4 , 5 , 6 , 7 ] Performance comparisons between QA devices and CIMs for various Ising models, such as complete, dense, and sparse graphs, have been reported. [ 8 ] Furthermore, theoretical performance comparisons between ideal gate-model quantum computers, implementing either Grover's search algorithm or the adiabatic quantum computing algorithm, and CIMs have been reported recently. [ 9 ] Although CIMs with all-to-all coupling among spins are highly effective, the use of an external FPGA circuit as well as an analog-todigital converter (ADC) and a digital-to-analog converter (DAC) not only results in considerable energy dissipation but also introduces a potential bottleneck for high-speed operation.
The standard linear coupling scheme of CIMs has been found to suffer from amplitude heterogeneity among the constituent quantum oscillators. Consequently, the Ising Hamiltonian is incorrectly mapped to the network loss, resulting in unsuccessful operation, especially in frustrated spin systems. [ 10 ] A novel error-correcting feedback scheme has been developed to resolve this problem [ 11 , 12 ] , which makes the solution accuracy of CIMs comparable to that of state-of-the-art heuristics such as break-out local search (BLS). [ 14 ] In this paper, we introduce a novel CIM architecture in which the error correction is implemented optically. In the proposed architecture, computationally intensive matrix-vector multiplication (MVM) and a nonlinear feedback function are implemented using phase-sensitive (degenerate) optical
parametric amplifiers, which are essentially the same device as the main-cavity optical parametric oscillator (OPO). This new CIM architecture can potentially be implemented monolithically in future photonic integrated circuits using thin-film LiNbO 3 platforms. LiNbO 3 platforms. [ 15 ]
A network of open dissipative quantum oscillators with optical error correction circuits is promising not only as a future hardware platform but also as a quantum-inspired algorithm because of its simple and efficient theoretical description. Numerical simulation of the time evolution of an N-qubit quantum system requires 2 N complex-number amplitudes. However, for a quantum oscillator network, various phasespace techniques of quantum optics have been developed over the last four decades. [ 18 , 19 , 20 ] The complete description of a network of quantum oscillators is now possible using N (or 2 N ) sets of stochastic differential equations (SDEs) based on positive-P, [ 21 ] truncated Wigner [ 22 , 23 , 24 ] or truncated Husimi [ 23 , 24 ] representations of the master equations. These SDEs can be used as heuristic algorithms on modern digital platforms. To completely described a network of low-Q quantum oscillators, a discrete map technique using a Gaussian quantum model is available, which is also computationally efficient. [ 25 ]
Similarly, a network of dissipation-less quantum oscillators with adiabatic Hamiltonian modulation is described using a set of N deterministic equations, which can also be used as a heuristic algorithm on modern digital platform. [ 27 , 28 , 29 ] Such heuristic algorithms are called simulated bifurcation machines (SBMs), [ 27 , 29 ] a variant of which will be studied in Section 6. Although the original SBM is inspired by dissipation-less adiabatic quantum computation, the version of SBM discussed in this paper (dSBM) is not a true unitary system, as dissipation is artificially added using inelastic walls to improve the performance of the algorithm. As both algorithms involve MVM as a computational bottleneck when simulated on a digital computer, we use the number of MVMs as the metric for performance comparison. We find that both types of systems have very similar performance in most cases, except graph types with great variation in vertex degree, where the SBM struggles consistently.
## 2 Semi-classical Model for Error Correction Feedback
In this section, we will describe several mutual coupling and error correction feedback schemes for CIMs. To simplify our argument, we consider a semi-classical deterministic picture. [ 10 ] The semi-classical model treated in this section is an approximate theory for the following fictitious machine. At an initial time t = 0, each signal pulse field is prepared in a vacuum state (Figure 1(a)), and each error pulse field is prepared in a weak coherent state (Figure 1(b)). When the pump fields p and p i are switched on at t ≥ 0, a vacuum field incident on the extraction beam splitter BS e from an open port is squeezed/anti-squeezed by a phase-sensitive amplifier (PSA) in this optical delay line (ODL) CIM, as shown in Figure 1(c). In other words, the vacuum fluctuation in the in-phase component ˜ X = 1 2 ( ˆ a +ˆ a † ) is deamplified by a factor of 1/G, while the vacuum fluctuation in the quadrature-phase component ˆ P = 1 2 i ( ˆ a -ˆ a † ) is amplified by a factor of G. Similarly, the vacuum fluctuations incident on the OPO pulse field owing to any linear loss of the cavity are all squeezed by the respective PSA. Moreover, the pump field and feedback injection field fluctuations along the in-phase component are also deamplified by the respective PSA (Figure (c)).
The truncated Wigner stochastic differential equation (W-SDE) for such a quantum-optic CIM with squeezed reservoirs has been derived and studied previously. [ 31 ] This particular CIM achieves the maximum quantum correlation among OPO pulse fields along the in-phase component as well as the maximum success probability, [ 31 ] because the quantum correlation among OPO pulse fields is formed by the mutual coupling of the vacuum fluctuations of OPO pulse fields without the injection of uncorrelated fresh reservoir noise in such a system. The following semi-classical model is considered as an approximate theory of the above-mentioned W-SDE in the limit of a large deamplification factor ( G 1). A full quantum description of a more realistic CIM with optical error correction circuits (without reservoir engineering) is given in Section 5. To overcome the problem of amplitude heterogeneity in the CIM [ 10 ] , the addition of an auxiliary variable for error detection and correction has been proposed. [ 11 , 12 ] This system has been
Figure 1: (a) Vacuum state and squeezed vacuum state. (b) Coherent state and squeezed coherent state. (c) Conventional CIM with vacuum noise injected from reservoirs, and a modified CIM with suppressed reservoir noise.
<details>
<summary>Image 1 Details</summary>
### Visual Description
## Quantum State Diagrams and Cavity Setup
### Overview
The image presents three diagrams illustrating quantum states and a cavity setup. Diagrams (a) and (b) depict squeezed vacuum and squeezed coherent states in phase space, respectively. Diagram (c) shows a cavity setup involving a PSA (Phase Sensitive Amplifier), beam splitters (BS), and signal/error pulses.
### Components/Axes
**Diagram (a): Squeezed Vacuum**
* **Axes:**
* X-axis labeled "X"
* Y-axis labeled "P"
* **Annotations:**
* "squeezed vacuum" (located at the top-right of the plot)
* "vacuum state" (located at the bottom-right of the plot, pointing to a dashed circle at the origin)
* **Data:** A vertically elongated, red-shaded region centered at the origin.
**Diagram (b): Squeezed Coherent State**
* **Axes:**
* X-axis labeled "X"
* Y-axis labeled "P"
* **Annotations:**
* "squeezed coherent state" (located at the top-right of the plot)
* "vacuum state" (located at the bottom-right of the plot, pointing to a dashed circle at the origin)
* **Data:** A horizontally elongated, red-shaded region centered on the X-axis, slightly offset from the origin.
**Diagram (c): Cavity Setup**
* **Components:**
* "Main Cavity": A large circle representing the cavity.
* "PSA": Phase Sensitive Amplifier, a rectangular box at the top of the cavity.
* "BS": Beam Splitters, two boxes at the bottom of the cavity.
* Input pulses: x1, x2, ..., xN, e1, e2, ..., eN entering the PSA.
* **Annotations:**
* "coherent state (squeezed coherent state)": Located above the PSA.
* "N signal pulses + N error pulses": Located to the right of the PSA, with an arrow indicating the direction of pulse propagation within the cavity.
* "vacuum fluctuation (squeezed vacuum)": Located below the main cavity.
* "coherent state (squeezed coherent state)": Located below the beam splitters.
* **Flow:** Pulses enter the PSA, circulate within the cavity (indicated by arrows), and interact with beam splitters.
### Detailed Analysis
**Diagram (a): Squeezed Vacuum**
* The squeezed vacuum state is represented by a probability distribution that is squeezed along the X-axis and elongated along the P-axis. This indicates reduced uncertainty in the X quadrature at the expense of increased uncertainty in the P quadrature. The dashed circle represents the vacuum state.
**Diagram (b): Squeezed Coherent State**
* The squeezed coherent state is similar to the squeezed vacuum state but is displaced from the origin along the X-axis. This indicates a non-zero mean value for the X quadrature.
**Diagram (c): Cavity Setup**
* The cavity setup involves a PSA that amplifies the signal pulses while suppressing the error pulses. The beam splitters are used to extract the coherent state from the cavity. The vacuum fluctuation represents the quantum noise within the cavity.
* The input pulses are labeled as x1, x2, ..., xN (signal pulses) and e1, e2, ..., eN (error pulses).
* The output from the beam splitters is labeled as f_i and x_i, e_i.
### Key Observations
* Diagrams (a) and (b) illustrate the concept of squeezing in quantum states, where uncertainty in one quadrature is reduced at the expense of increased uncertainty in the other.
* Diagram (c) shows a cavity setup designed to process quantum signals using a PSA and beam splitters.
### Interpretation
The image illustrates the manipulation of quantum states within a cavity. Squeezed states are used to reduce quantum noise and improve the sensitivity of quantum measurements. The cavity setup is designed to amplify signal pulses while suppressing error pulses, which is crucial for quantum communication and computation. The use of squeezed vacuum and squeezed coherent states suggests a strategy to enhance the signal-to-noise ratio in quantum information processing.
</details>
studied as a modification of the measurement feedback CIM. [ 31 ] The spin variable (signal pulse amplitude) x i and auxiliary variable (error pulse amplitude) e i obey the following deterministic equations: [ 11 ]
$$\frac { d x _ { i } } { d t } = - x _ { i } ^ { 3 } + \left ( p - 1 \right ) x _ { i } - e _ { i } \sum _ { j } \xi J _ { i j } x _ { j } ,$$
$$\frac { d e _ { i } } { d t } = - \beta e _ { i } \left ( x _ { i } ^ { 2 } - \alpha \right ) , & & ( 2 )$$
where J ij is the Ising coupling matrix, α , β and p are system parameters and ξ is a normalizing constant for J ij (see Appendix A for parameter selection). In many cases, we may modulate these parameters over time to achieve better performance (see Section 3 and Appendix C). To use this system as an Ising solver we consider the spin configuration σ i = sign( x i ) as a possible solution to the corresponding Ising problem. Even though noise is ignored in the above-mentioned equation, we can choose the initial x i amplitude randomly to create a diverse set of possible trajectories.
In this paper, we refer to this system of equations as CIM with chaotic amplitude control (CIM-CAC). The term 'chaotic' is used because CIM-CAC exhibits chaotic behavior (as discussed in Section 3). CIMCAC may refer to either the above-mentioned system of deterministic differential equations when integrated as a digital algorithm or an optical CIM that emulates the above-mentioned equations.
While studying the CIM-CAC equations, we have made the following modification:
$$z _ { i } = e _ { i } \sum _ { j } \xi J _ { i j } x _ { j } , & & ( 3 )$$
$$\frac { d x _ { i } } { d t } = - x _ { i } ^ { 3 } + \left ( p - 1 \right ) x _ { i } - z _ { i } ,$$
$$\frac { d e _ { i } } { d t } = - \beta e _ { i } \left ( z _ { i } ^ { 2 } - \alpha \right ) ,$$
which we refer to as CIM with chaotic feedback control (CIM-CFC). The only difference between Eqs. (2) and (5) is that the time evolution of the error variable e i monitors the feedback signal z i , rather than the internal amplitude x i . The dynamics of this new equation are very similar to those of CIM-CAC, which can be understood by observing that CIM-CAC and CIM-CFC have nearly identical fixed points. The motivation for studying CIM-CFC in addition to CIM-CAC is to gain a better understanding of
how these systems work. In addition, CIM-CFC may have slightly simpler dynamics, which simplifies its numerical integration.
The third system discussed in this paper has a very different equation:
$$z _ { i } = \sum _ { j } \xi J _ { i j } x _ { j } , & & ( 6 )$$
$$\frac { d x _ { i } } { d t } = - x _ { i } ^ { 3 } + \left ( p - 1 \right ) x _ { i } - f \left ( c z _ { i } \right ) - k \left ( z _ { i } - e _ { i } \right ) ,$$
$$\frac { d e _ { i } } { d t } = - \beta \left ( e _ { i } - z _ { i } \right ) .$$
The non-linear function f is a sigmoid-like function such as f ( z ) = tanh( z ), and p , k , c and β are system parameters (See Appendix A for parameter selection). The significance of this new feedback system is that the differential equation for the error signal e i is now linear in the 'mutual coupling signal' z i . In addition, z i is calculated simply as ∑ j ξJ ij x j without the additional factor e i as in Eq. (6). This means that the only nonlinear elements in this system are the gain saturation term -x 3 i and the nonlinear function f . For the results in this paper we will use f ( z ) = tanh( z ), however if a different function with the same properties is used the system will have similar behavior.
In the above-mentioned system, the two essential aspects of CIM-CAC and CIM-CFC are separated into two different terms. The term f ( cz i ) realizes mutual coupling while passively addressing the problem of amplitude heterogeneity, while the term k ( z i -e i ) introduces the error signal e i which helps to destabilize local minima. Therefore, we refer to this system as CIM with separated feedback control (CIM-SFC) in the remainder of this paper.
Another significant aspect of CIM-SFC (Eqs. (6),(7) and (8)) compared to CIM-CAC and CIM-CFC (Eqs. (1)-(5)) is that the auxiliary variables e i in CIM-SFC have a very different meaning. In CIM-CAC and CIM-CFC, e i is meant to be a strictly positive number that varies exponentially and modulates the mutual coupling term. In CIM-SFC, e i is instead a variable that stores sign information and takes the same range of values as the mutual coupling signal z i . The error signal e i can essentially be regarded as a low pass filter on z i , and the term k ( e i -z i ) can be regarded as a high pass filter on z i (in other words k ( e i -z i ) only registers sharp changes in z i ). The similarities and differences among CIM-SFC, CIM-CAC and CIM-CFC can be understood by observing the fixed points. In CIM-CAC and CIM-CFC, the fixed points are of the form: [ 11 ]
$$x _ { i } = \lambda _ { 1 } \sigma _ { i } ,$$
$$e _ { i } = \lambda _ { 2 } \frac { 1 } { h _ { i } \sigma _ { i } } , & & ( 1 0 )$$
$$h _ { i } = \sum _ { j } \xi J _ { i j } \sigma _ { j } . & & ( 1 1 )$$
with
Here, σ i is a spin configuration corresponding to a local minimum of the Ising Hamiltonian, and λ 1 and λ 2 are constants that depend on the system parameters. In CIM-SFC, the fixed points are generally very complicated and difficult to express explicitly. However, if we consider the limit of c 1, the fixed points will take the form:
$$x _ { i } = \lambda \sigma _ { i } , \tag* { ( 1 2 ) } x _ { i } = \lambda \sigma _ { i } ,$$
$$e _ { i } = \lambda h _ { i } , \quad ( 1 3 )$$
where λ is a number such that -λ 3 + ( p -1) λ = -1. Again, σ i is a spin configuration corresponding to a local minimum. This formula is only valid if f ( cz ) is an odd function that takes the value of +1 for cz 1 and -1 for cz 1. Therefore, it is important to choose an appropriate function f .
The important difference between the fixed points of these two types of systems is that in CIM-CAC and CIM-CFC, the error signal is
$$| e _ { i } | \varpropto \frac { 1 } { | h _ { i } | } ,$$
$$| e _ { i } | \varpropto | h _ { i } | \, .$$
In Section 5, we will see that this difference makes CIM-SFC more robust to quantum noise from reservoirs and pump sources. In the next section, we will investigate the similarities and differences among these three systems using numerical simulation.
## 3 Numerical Simulation of CIM-CAC, CIM-CFC and CIM-SFC
The originally proposed CIM architecture employs simple linear feedback without using an error detection/correction mechanism. In other words, the feedback term in Eq. (1) is simply ∑ j ξJ ij x j . [ 10 ] In this case, the Ising Hamiltonian cannot be properly mapped to the network loss owing to OPO amplitude heterogeneity, especially for frustrated spin systems, as shown in Figure 2 (a). Such a CIM often does not find a ground state of the Ising Hamiltonian; instead, it selects the lowest-energy eigenstate of the coupling (Jacobian) matrix [ J ij ].[10] This undesired operation is caused by a system's formation of heterogenous amplitudes.[10] We can address this problem partially by introducing a nonlinear filter function for the feedback pulse, such as tanh( ∑ j ξJ ij x j ). Thus, we can achieve the homogeneous OPO amplitudes, at least well above threshold, as shown in Figure 2 (b), and satisfy the proper mapping condition toward the end of the system's trajectory. However, such nonlinear filtering alone is not sufficiently powerful to prevent the machine from being trapped in numerous local minima. As the problem size N increases, NP-hard Ising problems are expected to have an exponentially increasing number of local minima; hence, a system that is easily trapped in these minima will be ineffective.
Figure 2: Trajectories of OPO amplitudes in CIMs with (a) linear feedback, (b) nonlinear filtering feedback, and (c) chaotic feedback control.
<details>
<summary>Image 2 Details</summary>
### Visual Description
## Chart Type: Multiple Time Series Plots Comparing CIM Configurations
### Overview
The image presents three time series plots comparing the pulse amplitude over time for different configurations of a CIM (presumably a computational or circuit model). The configurations are: (a) Conventional CIM, (b) CIM with nonlinear feedback, and (c) CIM with error correction (CFC). Each plot shows multiple lines, representing different pulse amplitudes evolving over a time period from 0 to 200.
### Components/Axes
**General for all three plots:**
* **X-axis:** "time", ranging from 0 to 200 in increments of 50.
* **Y-axis:** "pulse amplitude (x_i)", ranging from -2 to 2 in increments of 1.
* Each plot contains approximately 10 different colored lines, each representing a different pulse amplitude.
**Plot (a) Conventional CIM:**
* **Title:** "(a) Conventional CIM"
* The lines generally converge towards a stable state over time.
**Plot (b) CIM with nonlinear feedback:**
* **Title:** "(b) CIM with nonlinear feedback"
* The lines show a more stable behavior compared to (a), with less convergence.
**Plot (c) CIM with error correction (CFC):**
* **Title:** "(c) CIM with error correction (CFC)"
* The lines exhibit oscillatory behavior, indicating continuous correction and fluctuation.
### Detailed Analysis
**Plot (a) Conventional CIM:**
* Several lines start at different initial pulse amplitudes (between -1 and 1).
* Many lines converge towards the 0 amplitude level as time increases.
* Some lines initially increase or decrease before converging.
* One purple line starts near -0.5, decreases to approximately -1.5 around time 150, and then begins to increase slightly.
* One black line starts near 0.5 and increases to approximately 1.5.
**Plot (b) CIM with nonlinear feedback:**
* Lines start at various initial amplitudes, similar to plot (a).
* The lines tend to stabilize more quickly compared to plot (a).
* A dark blue line starts near 0.5 and increases steadily to approximately 1.5 by time 200.
* A green line starts near 0 and decreases slightly to approximately -0.25.
**Plot (c) CIM with error correction (CFC):**
* All lines exhibit oscillatory behavior, fluctuating between approximately -1 and 1.
* The oscillations appear to be somewhat synchronized, with peaks and troughs occurring at roughly the same time for different lines.
* The amplitude of the oscillations seems relatively consistent over the entire time period.
### Key Observations
* **Convergence vs. Oscillation:** The conventional CIM tends to converge to a stable state, while the CIM with error correction oscillates. The CIM with nonlinear feedback shows a more stable behavior than the conventional CIM.
* **Stability:** The nonlinear feedback seems to improve stability compared to the conventional CIM.
* **Error Correction:** The error correction mechanism introduces oscillations, preventing the system from settling into a fixed state.
### Interpretation
The plots demonstrate the impact of different feedback mechanisms on the behavior of a CIM. The conventional CIM tends to converge, suggesting it might be prone to settling into local minima. The nonlinear feedback improves stability, potentially avoiding these local minima. The error correction mechanism introduces oscillations, which could be a strategy to continuously explore the state space and avoid getting stuck. The choice of which configuration is "best" likely depends on the specific application and desired behavior of the CIM.
</details>
To destabilize the attractors caused by local minima and allow the machine to continue searching for a true ground state, we can introduce an error detection/correction variable expressed by Eq. (2) or (5). As shown in Figure 2 (c), the trajectory of a CIM with error correction variables will not reach equilibrium but continue to explore many states. Conversely, the systems in Figure 2 (a) and (b), which do not have an error correction variable ( e i ), will often converge rapidly on a fixed point corresponding to a high-energy excited state of the Ising Hamiltonian. Destabilizing the local minima will inevitably make the ground state unstable as well. Although this is undesirable, we can simply allow the machine to visit whereas in CIM-SFC, the error signal is
many local minima and then determine the one that has the lowest energy subsequently. Alternatively, we have found that by modulating the system parameters, the system will have a high probability of staying in a ground state toward the end of the trajectory (see Section 4 for further details).
The addition of another N degrees of freedom allows the machine to visit a local minimum, escape from it, and continue to explore nearby states; this is not possible with a conventional CIM algorithm. In this section, we will discuss the dynamics of the error correction schemes proposed in this paper: CIM-CFC and CIM-SFC.
Even though CIM-CFC and CIM-SFC are described by very different equations, the two systems were originally conceived through a similar concept. To understand why CIM-CFC and CIM-SFC are similar, we can consider these systems as follows. We introduce the 'mutual coupling signal' M i ( t ) = ∑ j ξJ ij x j ( t ) and the 'injection feedback signal' I i ( t ). Then, we can write both CIM-CFC and CIM-SFC in the form:
$$M _ { i } ( t ) = \sum _ { j } \xi J _ { i j } x _ { j } ( t ) ,$$
$$\frac { d x _ { i } } { d t } = - x _ { i } ^ { 3 } + ( p - 1 ) \, x _ { i } - I _ { i } ( t ) ,$$
where I i ( t ) depends on the time evolution of M i ( t ). Figure 3 shows how I i ( t ) (red) varies with respect to a mutual coupling field M i ( t ) (blue) for four different feedback schemes.
Figure 3: Mutual coupling field (blue) and injection feedback field (red) in four different feedback systems.
<details>
<summary>Image 3 Details</summary>
### Visual Description
## Diagram: Comparison of CIM Models
### Overview
The image presents a comparative analysis of four different Coupled-Inhibitory Mutator (CIM) models: Conventional CIM, CIM with nonlinear feedback, CIM-CFC, and CIM-SFC. For each model, the image provides the governing equation, a plot illustrating the behavior of the mutual coupling field and injection feedback field over time, and a brief description of the problem addressed by the model.
### Components/Axes
* **Columns:** The image is divided into four columns, each representing a different CIM model.
* Column 1: Conventional CIM
* Column 2: CIM with nonlinear feedback
* Column 3: CIM-CFC
* Column 4: CIM-SFC
* **Rows:** The image is divided into three rows for each model:
* Row 1: Equation defining the model
* Row 2: Plot of I(t) or M(t) vs. time
* Row 3: Problem addressed by the model
* **Plots:** Each plot has the following axes:
* Y-axis: Labeled "I(t) or M(t)", ranging from -4 to 4 with tick marks at -3, -2, -1, 0, 1, 2, 3, and 4.
* X-axis: Labeled "time", ranging from 0 to 10 with tick marks at 2, 4, 6, 8, and 10.
* **Legend:** Located within each plot, indicating:
* Blue dashed line: "mutual coupling field"
* Red solid line: "injection feedback field"
### Detailed Analysis
**1. Conventional CIM**
* **Equation:** I(t) = M(t)
* **Plot:**
* The blue dashed line (mutual coupling field) starts at approximately -1, increases linearly to 1 at time = 2, then increases linearly to 3 at time = 4. It remains constant at 3 until time = 6, then decreases linearly to -1 at time = 8, and remains constant at -1 until time = 10.
* The red solid line (injection feedback field) is identical to the blue dashed line.
* **Problem:** amplitude heterogeneity
**2. CIM with nonlinear feedback**
* **Equation:** I(t) = tanh(M(t))
* **Plot:**
* The blue dashed line (mutual coupling field) is identical to the mutual coupling field in the Conventional CIM plot.
* The red solid line (injection feedback field) starts at approximately -0.8, increases to approximately 0.8 at time = 2, then increases to approximately 1 at time = 4. It remains constant at 1 until time = 6, then decreases to approximately -0.8 at time = 8, and remains constant at -0.8 until time = 10.
* **Problem:** amplitude heterogeneity removed, trapped by local minima
**3. CIM-CFC**
* **Equation:**
* I(t) = e(t)M(t)
* de(t)/dt = -βe(t)(e(t)^2M(t)^2 - α)
* **Plot:**
* The blue dashed line (mutual coupling field) is identical to the mutual coupling field in the Conventional CIM plot.
* The red solid line (injection feedback field) starts at approximately -0.8, increases to approximately 1.5 at time = 2, then increases to approximately 2 at time = 4. It decreases to approximately -1 at time = 6, and remains constant at -1 until time = 10.
* **Problem:** amplitude heterogeneity removed, local minima destabilized
**4. CIM-SFC**
* **Equation:**
* I(t) = tanh(cM(t)) + k(M(t) - e(t))
* de(t)/dt = -β(e(t) - M(t))
* **Plot:**
* The blue dashed line (mutual coupling field) is identical to the mutual coupling field in the Conventional CIM plot.
* The red solid line (injection feedback field) starts at approximately -0.8, increases to approximately 1.5 at time = 2, then increases to approximately 1.5 at time = 4. It decreases to approximately -1 at time = 6, and remains constant at -1 until time = 10.
* **Problem:** amplitude heterogeneity removed, local minima destabilized
### Key Observations
* The "mutual coupling field" (blue dashed line) is identical across all four models.
* The "injection feedback field" (red solid line) varies significantly between the models, reflecting the different equations and feedback mechanisms.
* The Conventional CIM model has identical mutual coupling and injection feedback fields.
* The CIM with nonlinear feedback model compresses the amplitude of the injection feedback field due to the tanh function.
* The CIM-CFC and CIM-SFC models show a more complex behavior of the injection feedback field, with a sharper decrease after time = 6.
### Interpretation
The image illustrates how different mathematical formulations of the CIM model affect the behavior of the injection feedback field. The Conventional CIM model serves as a baseline, where the injection feedback field mirrors the mutual coupling field. The subsequent models introduce modifications to address specific problems, such as amplitude heterogeneity and the trapping of local minima. The nonlinear feedback model compresses the amplitude, while the CIM-CFC and CIM-SFC models destabilize local minima, potentially leading to more dynamic and adaptable system behavior. The equations provided offer a mathematical description of these behaviors, while the plots provide a visual representation of the system's response over time.
</details>
The similarity between CIM-CFC and CIM-SFC is as follows: if the mutual coupling field M i ( t ) remains constant for a certain period of time, then the injection feedback field I i ( t ) will converge on the value given by sign( M i ( t )). However, if M i ( t ) varies sharply, then I i ( t ) will deviate from its steady state values: +1 / -1. This small deviation is effective for triggering destabilization when the system is near a local minimum, which allows the machine to explore new spin configurations.
Although CIM-CFC and CIM-SFC were conceived on the basis of the same principle, the dynamics of the two systems seem to differ from each other. In particular, CIM-CFC (and CIM-CAC) nearly always features chaotic dynamics, as the trajectory is highly sensitive to the initial conditions. In the case of CIM-SFC, the trajectory will often immediately fall into a stable periodic orbit unless the parameters are dynamically modulated. At present, we do not have an exact theoretical reason for this difference in dynamics; this is purely an experimental observation. A more theoretical analysis in the case of CIMCAC can be found elsewhere. [11].
To demonstrate this difference, Figure 4 shows the correlation of pulse amplitudes between two initial conditions that are very close to each other. An initial condition for the pulse amplitude #1 (plotted on
the x-axis) is chosen from a zero-mean Gaussian with a standard deviation of 0.25, while the other initial condition for trajectory #2 (plotted in y-axis) is equal to that of trajectory #1 plus a small amount of noise (standard deviation 0.01).
Figure 4: Signal pulse amplitude correlations at different evolution time in CIM-SFC and CIM-CFC.
<details>
<summary>Image 4 Details</summary>
### Visual Description
## Scatter Plot Matrix: SFC vs. CFC Evolution Over Time
### Overview
The image presents a 2x4 matrix of scatter plots, comparing the evolution of two conditions, SFC and CFC, over time (T = 0, 10, 100, 4000). Each plot shows the relationship between two initial conditions, labeled as #1 and #2, with the x-axis representing condition #1 and the y-axis representing condition #2. The plots illustrate how the distribution of data points changes over time for each condition.
### Components/Axes
* **Rows:**
* Top Row: SFC (initial condition #2)
* Bottom Row: CFC (initial condition #2)
* **Columns:**
* Column 1: T = 0
* Column 2: T = 10
* Column 3: T = 100
* Column 4: T = 4000
* **Axes:**
* X-axis (horizontal): *x<sub>i</sub>* (initial condition #1). Scale ranges from -1.0 to 1.0, with tick marks at -1.0, -0.5, 0.0, 0.5, and 1.0.
* Y-axis (vertical): *x<sub>i</sub>* (initial condition #2). Scale ranges from -1.0 to 1.0, with tick marks at -1.0, -0.5, 0.0, 0.5, and 1.0.
* **Data Points:** Each plot contains multiple orange data points.
### Detailed Analysis
**SFC (Top Row):**
* **T = 0:** Data points form a relatively tight, positively sloped linear cluster.
* **T = 10:** The cluster remains linear and positively sloped, but the points are slightly more dispersed than at T = 0.
* **T = 100:** The linear relationship weakens, and the points become more scattered.
* **T = 4000:** The points are concentrated in the corners of the plot, forming four distinct clusters.
**CFC (Bottom Row):**
* **T = 0:** Similar to SFC at T = 0, the data points form a tight, positively sloped linear cluster.
* **T = 10:** The cluster remains linear and positively sloped, with a slight increase in dispersion compared to T = 0.
* **T = 100:** The points are significantly more scattered, forming a complex, non-linear pattern.
* **T = 4000:** The points are concentrated near the edges and corners of the plot, forming a more defined square-like shape with clusters in the corners.
### Key Observations
* Both SFC and CFC start with a strong positive correlation between initial conditions #1 and #2.
* As time progresses, the correlation weakens, and the data points become more dispersed.
* At T = 4000, both conditions exhibit clustering behavior, but the patterns are different. SFC clusters primarily in the corners, while CFC forms a more distributed pattern along the edges and corners.
### Interpretation
The plots illustrate the evolution of two systems (SFC and CFC) from an initial state of positive correlation to more complex and potentially chaotic states over time. The differences in the final distributions at T = 4000 suggest that the two systems evolve differently, possibly due to different underlying dynamics or parameters. The initial linear relationship indicates a strong dependency between the two initial conditions, which diminishes as the systems evolve independently. The clustering at later times suggests the emergence of stable states or attractors within the systems' phase space. The data suggests that while both systems start similarly, they diverge significantly over time, indicating different long-term behaviors.
</details>
Figure 4, shows the correlation of all 100 pulse amplitudes between the two initial conditions for a SherringtonKirkpatrick (SK) spin glass instance of problem size N = 100. In CIM-SFC (first row), the correlation remains even after 4000 time steps (round trips), which means that two initial conditions follow a nearly identical trajectory. However, in CIM-CFC (second row), we see that the xi variables become uncorrelated after around 100 time steps, even though the initial conditions of the two trajectories are very close. This indicates qualitatively that CIM-CFC is highly sensitive to the initial condition, whereas CIM-SFC is not.
This pattern tends to hold when different parameters and initial conditions are used. However, although CIM-SFC stays correlated in most cases, the two trajectories diverge under certain system parameters and initial conditions. This means that although CIM-SFC is less sensitive to the initial conditions compared to CIM-CFC, some chaotic dynamics likely occur during the search, especially when the parameters are modulated.
Another way to qualitatively observe the difference in dynamics is to simply observe the trajectories. Figure 5, shows examples of trajectories of both systems (10 out of 100 x i variables are shown) with fixed and linearly modulated system parameters. When the parameters are fixed, the difference between the two systems evident. CIM-SFC will rapidly become trapped in a stable periodic attractor, while CIMCFC will continue to search in an unpredictable manner. Therefore, the parameters are slowly modulated in CIM-SFC so that the system can find a ground state. CIM-CFC and CIM-CAC can find ground states with fixed parameters. However, we have found that modulation of the system parameters improves the performance of CIM-CFC and CIM-CAC considerably (see Appendix C for details).
In the lower left panel of Figure 5, the parameters c and p of CIM-SFC are linearly increased from low to high values ( p ranges from -1 to +1 and c ranges from 1 to 3). We can see that as the parameters change, the system may jump from one attractor to another and eventually end up in a fixed point/local minimum. By linearly increasing the parameters c and p from a low to high value in CIM-SFC, we are slowly transitioning the nonlinear term tanh( cz i ) from a 'soft spin' mode where the nonlinear coupling term has a continuous range of values between -1 and 1 to a 'discrete' mode where tanh( cz i ) will mostly take on the values +1 or -1. This transition is seems to be crucial for CIM-SFC to function properly.
For most fixed parameters, CIM-SFC rapidly approaches a periodic or fixed point attractor as shown in Figure 5; however, as mentioned earlier it is likely that for some specific values of c and p , CIM-SFC will feature chaotic dynamics similar to CIM-CFC. It has been shown [ 11 , 13 ] that chaotic dynamics are observed when solving hard optimization problems efficiently using a deterministic system. This trend is also observed in the simulated bifurcation machine [27, 29]. Whether or not CIM-SFC utilizes chaotic dynamics is beyond the scope of this paper. Whether CIM-SFC uses chaotic dynamics is beyond the
Figure 5: Signal pulse amplitude trajectories of CIM-SFC and CIM-CFC with fixed and modulated system parameters.
<details>
<summary>Image 5 Details</summary>
### Visual Description
## Time Series Plots: CIM-SFC vs. CIM-CFC
### Overview
The image presents four time series plots, arranged in a 2x2 grid. The plots compare two models, CIM-SFC and CIM-CFC, under two parameter conditions: "Fixed Parameters" (top row) and "Modulated Parameters" (bottom row). The y-axis represents the value of variables (x_i), and the x-axis represents time steps (T). Each plot shows 10 out of 100 variables.
### Components/Axes
* **Titles:**
* Top-Left: CIM-SFC (Fixed Parameters)
* Top-Right: CIM-CFC (Fixed Parameters)
* Bottom-Left: CIM-SFC (Modulated Parameters)
* Bottom-Right: CIM-CFC (Modulated Parameters)
* **Y-Axis:**
* Label (Left Side, Top Plots): "Fixed Parameters"
* Label (Left Side, Bottom Plots): "Modulated Parameters"
* Label (Right Side, All Plots): "x_i variables 10/100 shown"
* Scale: -1.0 to 1.0, with markers at -1.0, -0.5, 0.0, 0.5, and 1.0
* **X-Axis:**
* Label (All Plots): "T (time steps)"
* Scale (CIM-SFC Plots): 0 to 5000, with markers at 0, 1000, 2000, 3000, 4000, and 5000
* Scale (CIM-CFC Plots): 0 to 1000, with markers at 0, 200, 400, 600, 800, and 1000
* **Data Series:** Each plot contains multiple time series, each represented by a different color. There is no explicit legend.
### Detailed Analysis
**1. CIM-SFC (Fixed Parameters) - Top-Left Plot:**
* Trend: The time series exhibit oscillatory behavior with a consistent amplitude and frequency.
* Values: The oscillations range approximately from -0.8 to 0.8.
* Colors: Multiple colors are present, including blue, green, red, black, magenta, and cyan.
**2. CIM-CFC (Fixed Parameters) - Top-Right Plot:**
* Trend: The time series show irregular oscillations with varying amplitudes and frequencies.
* Values: The oscillations range approximately from -0.8 to 0.8.
* Colors: Multiple colors are present, including blue, green, red, black, magenta, and cyan.
**3. CIM-SFC (Modulated Parameters) - Bottom-Left Plot:**
* Trend: The time series start with irregular behavior, then converge towards a stable state. Some series approach 1.0, while others approach -1.0.
* Values: Initial values range from -1.0 to 1.0. Final values converge near -1.0 and 1.0.
* Colors: Multiple colors are present, including blue, green, red, black, magenta, and cyan.
**4. CIM-CFC (Modulated Parameters) - Bottom-Right Plot:**
* Trend: The time series exhibit initial fluctuations, then converge towards a stable state. Most series approach -1.0.
* Values: Initial values range from -1.0 to 1.0. Final values converge near -1.0.
* Colors: Multiple colors are present, including blue, green, red, black, magenta, and cyan.
### Key Observations
* **Fixed Parameters:** CIM-SFC shows regular oscillations, while CIM-CFC shows irregular oscillations.
* **Modulated Parameters:** Both models converge to stable states, but CIM-SFC shows convergence to both positive and negative values, while CIM-CFC primarily converges to negative values.
* **Time Scale:** CIM-SFC plots cover a longer time scale (0-5000) compared to CIM-CFC plots (0-1000).
### Interpretation
The plots compare the behavior of two models, CIM-SFC and CIM-CFC, under different parameter conditions. The "Fixed Parameters" plots show the inherent oscillatory dynamics of the models. The "Modulated Parameters" plots demonstrate how the models respond to changing conditions, with both models eventually reaching stable states. The difference in convergence behavior between CIM-SFC and CIM-CFC under modulated parameters suggests that the models have different sensitivities or responses to the modulation. The longer time scale of the CIM-SFC plots may indicate that this model takes longer to reach a stable state under modulated conditions.
</details>
scope of this paper. To answer this question, we need to further analyze how the parameters affect the dynamics of CIM-SFC and gain a deeper understanding of how CIM-SFC finds ground states.
## 4 Implementation of CIM with Optical Error Correction Circuits
Figure 6, together with Figure 1(c), shows a physical setup for CIM-CAC and CIM-CFC with optical error correction circuits. In our design, the main ring cavity stores both signal pulses with normalized amplitude, x i and error pulses with normalized amplitude e i , where i = 1 , 2 , · · · N . The signal pulses start from vacuum states | 0 〉 1 | 0 〉 2 · · · | 0 〉 N and are amplified (or deamplified) along the X-coordinate by a positive (or negative) pump rate p .
The error pulses start from a coherent state | α 〉 1 | α 〉 2 · · · | α 〉 N with α > 0 and are amplified (or deamplified) along the X-coordinate by the pump rate p ′ as described below. The squared amplitude of the error pulses is kept small ( e 2 i < 1) compared to the saturation level of the main cavity OPO. Thus the error pulses are controlled in a linear amplifier/deamplifier regime while the signal pulses are controlled in both a linear amplifier/deamplifier regime ( x 2 i < 1) and a nonlinear oscillator regime ( x 2 i > 1).
An extraction beamsplitter (BS e shown in Figure 1(c)) selects partial waves of the signal and error pulses that are amplified by a noise-free phase-sensitive amplifier (PSA 0 as shown in Figure 6(a)). PSA 0 amplifies the signal and error pulses to a classical level without introducing additional noise. The extracted amplitudes ˜ x i and ˜ e i suffer from the signal-to-noise ratio (SNR) degradation owing to the vacuum noise incident on BS e . However, they are amplified by a high-gain noise-free phase-sensitive amplifier PSA 0 to classical levels; hence, no further SNR degradation occurs even with large linear losses in the optical error correction circuits.
A small part of the PSA 0 output is sent to an optical homodyne detector that measures the extracted signal and error pulses with amplitudes ˜ x i and ˜ e i , respectively. The measurement error of the homodyne detection is determined solely by the reflectivity of BS e and the vacuum fluctuation incident on BS e (as described above). Figure 6(a) shows the output of the fan-out circuit at different time instances t = τ, 3 τ, 5 τ, ... separated by a signal pulse to signal pulse interval of 2 τ .
For instance, the signal pulse ( ˜ x j ) is first input into PSA j and then sent to optical delay line DL j with a delay time of (2 N -2 j + 1) τ . The phase-sensitive gain/loss of PSA j is set to √ G j = ξJ ij so that the amplified/deamplified signal pulse that arrives in front of the fan-in circuit at time t = 2 Nτ is equal to ξJ ij ˜ x j . Therefore the fan-in circuit will output a pulse with the desired amplitude of ∑ j ξJ ij ˜ x j . Suppose that PSA j has a phase-sensitive linear gain/loss of 10dB, then we can implement an arbitrary Ising coupling of range 10 -2 < | ξJ ij | < 1.
Figure 6: (a) Optical implementation of error correction circuits for CIM-CAC and CIM-CFC. (b) Pump pulse factory providing SHG pulses to the main cavity, and the error correction circuits. The pump pulse factory carries N 2 pulses spread over N optical cavities corresponding to the elements of the Ising coupling matrix J ij .
<details>
<summary>Image 6 Details</summary>
### Visual Description
## Diagram: Optical Circuit Schematics
### Overview
The image presents two optical circuit schematics, labeled (a) and (b), detailing different configurations of optical components for signal processing and frequency comb generation. Both diagrams illustrate the flow of optical signals through various elements such as beam splitters, phase-sensitive amplifiers, electro-optic modulators, and second harmonic generators.
### Components/Axes
#### Diagram (a)
* **Input:** Beam splitter (BS<sub>e</sub>) with inputs x̄<sub>j</sub> and ē<sub>i</sub>
* **Components:**
* PSA<sub>0</sub> (Phase-Sensitive Amplifier)
* Fan-Out: Splits the signal into N paths.
* x̄<sub>1(τ)</sub>, x̄<sub>2(3τ)</sub>, ..., x̄<sub>N(2N-1τ)</sub>: Represent different signal paths after the Fan-Out.
* PSA<sub>1</sub> + DL<sub>1</sub>, PSA<sub>2</sub> + DL<sub>2</sub>, ..., PSA<sub>N</sub> + DL<sub>N</sub>: Phase-Sensitive Amplifiers with Delay Lines.
* Fan-In: Combines the signals from the N paths.
* ∑ξ<sub>ij</sub>ē<sub>i</sub><sup>(m)</sup>, ∑ξ<sub>ij</sub>x̄<sub>j</sub>: Represent the summed signals.
* PSA<sub>e</sub>: Phase-Sensitive Amplifier.
* Homodyne Detection: Measures the signal.
* **Output:** BS<sub>i</sub>
* **Results:** x̄<sub>j</sub><sup>(m)</sup>, ē<sub>i</sub><sup>(m)</sup> and z̄<sub>i</sub><sup>(m)</sup>
#### Diagram (b)
* **Input:** Soliton frequency comb generator.
* **Components:**
* PSA<sub>p</sub>: Phase-Sensitive Amplifier.
* Fan-Out: Splits the signal into multiple paths.
* EOM: Electro-Optic Modulators (multiple).
* EOM<sub>1</sub>, EOM<sub>2</sub>, ..., EOM<sub>N</sub>: Specific Electro-Optic Modulators.
* PSA<sub>s</sub> Ring Cavity<sub>1</sub>, PSA<sub>s</sub> Ring Cavity<sub>2</sub>, ..., PSA<sub>s</sub> Ring Cavity<sub>N</sub>: Phase-Sensitive Amplifiers within Ring Cavities.
* SHG: Second Harmonic Generators (multiple).
* Fan-Out: Splits the signal into multiple paths.
* **Output:**
* P, P<sub>i</sub>, PSA<sub>e</sub>, PSA<sub>1</sub>, PSA<sub>2</sub>, ..., PSA<sub>N</sub>, PSA<sub>s</sub>, PSA<sub>0</sub>: Represent different output signals.
* **Main cavity + Error correction circuit:** Encloses the output signals.
### Detailed Analysis or ### Content Details
#### Diagram (a)
The diagram illustrates a signal processing scheme. An input signal is split into multiple paths, each path undergoing phase-sensitive amplification and delay. These paths are then recombined and further amplified before being measured via homodyne detection.
* The input signal is split by the beam splitter (BS<sub>e</sub>) into two components, x̄<sub>j</sub> and ē<sub>i</sub>.
* The Fan-Out splits the signal into N paths.
* Each path has a PSA and a delay line (DL).
* The signals are recombined in the Fan-In.
* The output is measured using homodyne detection.
#### Diagram (b)
This diagram depicts a soliton frequency comb generator coupled with an error correction circuit. The comb generator's output is processed through a series of electro-optic modulators and second harmonic generators. Some paths include phase-sensitive amplifiers within ring cavities. The final outputs are fed into a main cavity and error correction circuit.
* The soliton frequency comb generator provides the initial signal.
* The signal is split by the Fan-Out.
* Multiple paths include EOMs and SHGs.
* Some paths include PSAs within ring cavities.
* The final outputs are fed into the "Main cavity + Error correction circuit".
### Key Observations
* Diagram (a) focuses on signal processing using phase-sensitive amplification and delay lines.
* Diagram (b) focuses on frequency comb generation and error correction.
* Both diagrams use Fan-Out and Fan-In components to split and combine signals.
* Diagram (b) includes ring cavities and second harmonic generation, which are not present in diagram (a).
### Interpretation
The diagrams illustrate two different optical circuit configurations. Diagram (a) likely represents a scheme for manipulating and measuring optical signals, while diagram (b) represents a more complex system for generating and stabilizing a frequency comb. The use of phase-sensitive amplifiers, delay lines, electro-optic modulators, and second harmonic generators suggests advanced techniques for controlling the phase, frequency, and amplitude of optical signals. The error correction circuit in diagram (b) indicates a focus on achieving high stability and precision in the frequency comb generation process.
</details>
Next, the output of the fan-in circuit is input into another phase-sensitive amplifier PSA e that amplifies with a factor of √ G e = ˜ e i . This is achieved by modulating the pump power to PSA e based on the measurement result for ˜ e i . Finally, the output of PSA e is injected back into signal pulse ( x i ) of the main cavity via BS i (see Figure 1(c)). The extraction beamsplitter BS e outputs not only signal pulses but also error pulses that are used only for homodyne detection. Thus, we switch off the pump power to PSA 0 for the error pulses and deamplify the residual error pulses by PSA 1 , PSA 2 , ... , PSA N , PSA e . In this way we avoid any spurious injection of the error pulse back into the main cavity. The dynamics of the error pulse are governed solely by the pump power p i to the main cavity PSA, which is set to satisfy p i -1 = β ( α -˜ x i 2 ) or p i -1 = β ( α -˜ z i 2 ) .
One advantage of this optical implementation of CIM-CAC and CIM-CFC is that only one type of active device, a noise-free phase-sensitive (degenerate optical parametric) amplifier, and all the other elements are passive devices. This fact may allow for on-chip monolithic integration of the CIM system as well as low-energy dissipation in the computational unit, which will be discussed in Section 5.
A similar optical implementation of CIM-SFC is shown in Appendix F.
Figure 6(b) shows a pump pulse factory that provides the second harmonic generation (SHG) pulses to the main cavity PSA, post-amplifier PSA 0 , delay line amplifiers PSA 1 , PSA 2 , · · · PSA N and exit amplifier PSA e . The purpose of this pump pulse factory is to reduce the use of EOM modulators, which consume the most energy in the entire CIM. A soliton frequency comb generator produces a pulse train at a repetition frequency of 100 GHz and wavelength of 1.56 µ m wavelength. Before it is split into many branches, the pulse train is amplified by a pump amplifier PSA p . N storage ring cavities continuously produce the pump pulses for PSA 1 , PSA 2 , · PSA N in order to implement the MVM ∑ ξJ ij ˜ x j . For this
purpose, the pulses stored in the i-th ring cavity acquire the appropriate amplitudes to realize the gain √ G ij = ξJ ij . The time duration for using N EOM arrays is only one round trip of the ring cavity, i.e., N × 10 (psec). The out-coupling loss of the storage ring cavities is compensated for by the linear gain of the internal PSA s . The pump pulses for PSA p , PSA s , and PSA 0 have constant amplitudes and are hence driven directly by the PSA p output. The pump pulses P and P i for the signal and error pulses in the main cavity, as well as the exit PSA e , must be modulated during the entire computation time.
Another detail that needs to be accounted for when considering an optical implementation is the calculation of the Ising energy. In our digital simulation for generating the results presented in this paper, the Ising energy is calculated at every time step (round trip), and the smallest energy obtained is used as the result of the computation. This means that, in an optical implementation, we must measure the ˜ x i amplitude in every round trip and calculate the Ising energy using, for instance, an external ADC/FPGA circuit. This would defeat the purpose of using optics, as the digital circuit in the ADC/FPGA would then become a bottleneck in terms of time and energy consumption.
However, we have found that with proper parameter modulation as shown in Figure 5, it is possible to use only the final state of the system for the result and still have a high success probability. For the results on 800-spin Ising instances (SK model) presented in Section 6, we calculated how often a successful trajectory is in the ground spin configuration after the final time step. We found that, for CIM-SFC, in 100% of the 7401 successful trajectories the final spin configuration was in the ground state. In other words, if CIM-SFC visits the ground state at any point during the trajectory, then it will also be in the ground state at the end of the trajectory. Meanwhile, for CIM-CFC and CIM-CAC, this was true only in 75% of the time and 48% of the time, respectively. We believe that this difference among the three systems is a result of both the intrinsic dynamics and the parameters used.
This suggests that in a CIM with optical error correction, we can simply digitize the final measurement result of x i after many round trips to obtain the computational result, and still have a high success probability. In the case of CIM-CFC and CIM-CAC, it might be beneficial to read the spin configuration multiple times during the last few round trips, as the machine usually visits the ground state close to the end of the trajectory even if it does not stay there.
## 5 Quantum Noise Analysis and Energy Cost to Solution
As we propose implementation of these dynamical systems on analog optical devices, it is important to investigate the extent to which the noise from the physical systems (in this case, quantum noise from pump sources and external reservoirs) will degrade the performance. In this section, we present quantum models based on our optical implementation.
In our optical implementation for CIM-CAC, the real-number signal pulse amplitude µ i (in unit of photon amplitude) (in units of photon amplitude) obeys the following truncated Wigner SDE: [ 22 , 23 ]
$$\frac { d } { d t } \mu _ { i } = ( p - 1 ) \, \mu _ { i } - g ^ { 2 } \mu _ { i } ^ { 3 } + \tilde { \nu } _ { i } \sum _ { j } \xi J _ { i j } \tilde { \mu } _ { j } + n _ { i } ,$$
where the term pµ i represents the parametric linear gain and the term -µ i represents the linear loss rate; this includes the cavity background loss and extraction/injection beam splitter loss for mutual coupling and error correction. The nonlinear term -g 2 µ 3 i represents gain saturation (or back-conversion from signal to pump), where g is the saturation parameter. The saturation photon number is given by 1 /g 2 , which is equal to the average photon number of a solitary OPO at a pump rate of p = 2 (two times above the threshold). Furthermore, J ij is the ( i, j ) element of the N × N Ising coupling matrix, as described in Section 2. The time t is normalized by a linear loss rate; hence, the signal amplitude decays by a factor of 1 /e at time t = 1. In addition, ˜ µ j = µ j + ∆ µ j and ˜ ν i = ν i + ∆ ν i are the inferred amplitudes for the signal pulse and error pulse, respectively, and ∆ µ j and ∆ ν i represent the additional noise governed by vacuum fluctuations incident on the extraction beam splitter. They are characterized
by √ 1 -R B 4 R B w , where R B is the reflectivity of the extraction beam splitter and w is a zero-mean Gaussian random variable with a variance of one. Finally, n i is the noise injected from external reservoirs and pump sources. [ 22 , 23 ] It is characterized by the two time correlation functions 〈 n i ( t ) n i ( t ′ ) 〉 = ( 1 2 + g 2 µ 2 i ) δ ( t -t ′ ). We assume that the external reservoirs are in vacuum states and that the pump fields are in coherent states.
The real number error pulse amplitude ν i (in units of photon amplitude) is governed by
$$\frac { d } { d t } \nu _ { i } = \left ( p ^ { \prime } _ { i } - 1 \right ) \nu _ { i } + m _ { i } ,$$
where the correlation function for the noise term is given by 〈 m i ( t ) m i ( t ′ ) 〉 = 1 2 δ ( t -t ′ ). The pump rate p ′ i for the error pulse is determined by the inferred signal pulse amplitude ˜ x i = g ˜ µ i normalized by the saturation parameter,
$$p _ { i } ^ { \prime } - 1 = \beta \left ( \alpha - \tilde { x } _ { i } ^ { 2 } \right ) .$$
The error pulses start from coherent states | γ 〉 1 | γ 〉 2 · · · | γ 〉 N , for some positive real number 1 /g γ > 0. The absence of a gain saturation term in Eq. (17) implies that the error pulses are always pumped at below the threshold. Nevertheless, the error pulses represent exponentially varying amplitudes.
The parameter β governs the time constant for the error correction dynamics, and α is the squared target amplitude. This feedback model stabilizes the squared signal pulse amplitude ˜ x 2 i = g 2 ˜ µ 2 i to α through an exponentially varying error pulse amplitude e i = gν i . Eqs. (16) and (17) are rewritten for the normalized amplitudes x i and e i as
$$\frac { d } { d t } x _ { i } = \left ( p - 1 \right ) x _ { i } - x _ { i } ^ { 3 } + \tilde { e } _ { i } \sum _ { j } \xi J _ { i j } \tilde { x } _ { j } + g n _ { i } ,$$
$$\frac { d } { d t } e _ { i } = \left ( p ^ { \prime } _ { i } - 1 \right ) e _ { i } + g m _ { i } .$$
which are nearly identical to Eqs. (1) and (2) except for the noise terms.
CIM-CFC is also realized using the experimental setup shown in Figure 6. In this case, the relevant truncated Wigner SDE for the error pulse amplitude is still given by Eq. (17) or (20); however, the pump rate p ′ should be modified to
$$p _ { i } ^ { \prime } - 1 = \beta \left ( \alpha - \tilde { z } _ { i } ^ { 2 } \right ) .$$
$$\tilde { z } _ { i } = \sum _ { j } \xi J _ { i j } \tilde { x } _ { j } & & ( 2 2 )$$
with
Finally, CIM-SFC can also be realized using the experimental setup shown in Appendix F (Figure 17). In this case, Eqs. (19) and (20) should be modified as
$$\tilde { z } _ { i } = \sum _ { j } \xi J _ { i j } \tilde { x } _ { j } & & ( 2 3 )$$
$$\frac { d } { d t } x _ { i } = \left ( p - 1 \right ) x _ { i } - x _ { i } ^ { 3 } + k \left ( \tilde { e } _ { i } - \tilde { z } _ { i } \right ) + \tanh \left ( c \tilde { z } _ { i } \right ) + g n _ { i } ,$$
$$\frac { d } { d t } e _ { i } = - \beta \left ( e _ { i } - \tilde { z } _ { i } \right ) + g m _ { i } .$$
If we compare the semi-classical nonlinear dynamical models of CIM-CAC, CIM-CFC, and CIM-SFC, represented by Eqs. (1)-(8), with the quantum nonlinear dynamical models (truncated Wigner SDE), represented by Eqs. (19)-(25), we find that the main difference is the absence or presence of the vacuum noise and pump noise terms gn i and gm i , respectively. The other important difference is that ˜ x i and ˜ e i
Figure 7: Success probability P s vs. saturation parameter g 2 for CIM-SFC and CIM-CFC at N=100. The success probability is averaged over 100 SK instances. CIM-CAC is not shown; however, the result is nearly identical to that of CIM-CFC.
<details>
<summary>Image 7 Details</summary>
### Visual Description
## Chart Type: Line Graphs Comparing CIM-SFC and CIM-CFC
### Overview
The image presents two line graphs side-by-side, comparing the Average Success Probability (N=100) against the saturation parameter (σ²) for two different configurations: CIM-SFC (left) and CIM-CFC (right). Both graphs also show a "no noise" baseline as a horizontal dotted line. The x-axis (saturation parameter) is displayed on a logarithmic scale.
### Components/Axes
* **Titles:**
* Left Graph: CIM-SFC
* Right Graph: CIM-CFC
* **Y-Axis (Both Graphs):**
* Label: Average Success Probability (N=100)
* Scale: 0.0 to 1.0, with increments of 0.2 (0.0, 0.2, 0.4, 0.6, 0.8, 1.0)
* **X-Axis (Both Graphs):**
* Label: σ² (saturation parameter)
* Scale: Logarithmic scale from 10⁻⁸ to 10⁰ (10⁻⁸, 10⁻⁷, 10⁻⁶, 10⁻⁵, 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, 10⁰)
* **Legend (Both Graphs, Top-Left):**
* "no noise" - Represented by a gray dotted horizontal line.
### Detailed Analysis
**Left Graph: CIM-SFC**
* **Blue Line (Data Series):**
* Trend: The line remains relatively flat at approximately 0.55 from 10⁻⁸ to around 10⁻³, then sharply declines.
* Data Points:
* 10⁻⁸: ~0.55
* 10⁻⁶: ~0.55
* 10⁻⁴: ~0.55
* 10⁻³: ~0.57
* ~3 * 10⁻³: ~0.57
* ~6 * 10⁻³: ~0.55
* ~1.5 * 10⁻²: ~0.42
* ~3 * 10⁻²: ~0.09
* **Gray Dotted Line ("no noise"):**
* Horizontal line at approximately 0.55.
**Right Graph: CIM-CFC**
* **Blue Line (Data Series):**
* Trend: The line remains relatively flat at approximately 0.70 from 10⁻⁸ to around 10⁻⁴, then sharply declines.
* Data Points:
* 10⁻⁸: ~0.70
* 10⁻⁶: ~0.70
* 10⁻⁵: ~0.70
* 10⁻⁴: ~0.70
* ~3 * 10⁻⁴: ~0.70
* ~6 * 10⁻⁴: ~0.68
* ~1.5 * 10⁻³: ~0.55
* ~3 * 10⁻³: ~0.10
* ~6 * 10⁻³: ~0.01
* **Gray Dotted Line ("no noise"):**
* Horizontal line at approximately 0.70.
### Key Observations
* Both CIM-SFC and CIM-CFC show a stable average success probability at low saturation parameter values.
* The average success probability drops sharply for both configurations as the saturation parameter increases.
* CIM-CFC maintains a higher average success probability than CIM-SFC at lower saturation parameter values.
* The "no noise" baseline is higher for CIM-CFC than for CIM-SFC.
* The sharp decline in success probability occurs at a lower saturation parameter value for CIM-SFC compared to CIM-CFC.
### Interpretation
The graphs illustrate the impact of the saturation parameter (σ²) on the average success probability of two different configurations, CIM-SFC and CIM-CFC, with a fixed number of trials (N=100). The "no noise" baseline represents the ideal performance without any saturation effects.
The data suggests that CIM-CFC is more robust to the effects of the saturation parameter than CIM-SFC, as it maintains a higher success probability for a wider range of σ² values. The sharp decline in success probability indicates a critical threshold for the saturation parameter, beyond which the performance of both configurations degrades significantly. The difference in the "no noise" baselines suggests that CIM-CFC inherently performs better than CIM-SFC in the absence of saturation effects.
The trends indicate that controlling the saturation parameter is crucial for optimizing the performance of both CIM-SFC and CIM-CFC. The choice between the two configurations may depend on the expected range of saturation parameter values in a given application.
</details>
are inferred amplitudes with the vacuum noise contribution in the quantum model, whereas in the semiclassical model, the amplitudes x i and e i can be reproduced without additional noise.
Next, we will discuss the impact of quantum noise on the performance of CIM. As indicated in Eqs. (19)-(25), the relative magnitude of the quantum noise in the signal and error pulses is governed by the saturation parameter g . When g increases, the ratio between the normalized pulse amplitudes ( x i , e i ) and normalized quantum noise amplitudes ( gn i , gm i ) decreases. Therefore, the CIM performance is expected to degrade as g increases. However, as g increases, the OPO threshold pump power decreases (see Figure C1 in [31]), which suggests that the OPO energy cost to solution can be potentially reduced with increasing g .
Figure 7 shows the success probability P s for N = 100 Ising problems (SK model) plotted against the saturation parameter g 2 . The reflectivity of the extraction beam splitter R B is assumed to be R B = 0 . 1. The success probability P s is almost independent of the saturation parameter g 2 as long as g 2 10 -4 . However, when g 2 exceeds 10 -3 , the success probability drops rapidly owing to the decreased signal-toquantum noise ratio, as mentioned above.
Figure 8 shows the energy cost to solution for Ising problems (SK model) with N = 100 and N = 800, where we consider only the pump power to the main cavity PSA: E main = 2 ω (MVM) N ∆ t/g 2 , where MVM is the number of matrix-vector multiplication steps to solution and ∆ t is a round-trip time normalized by the signal lifetime ( ∼ 0 . 1).
Figure 8: Energy cost to solution in joules of CIM-SFC and CIM-CFC considering only pump power to main cavity PSA. The median ETS is plotted as a function of g 2 for N=100 and N=800 SK instances to show the optimal value of g 2 in each case.
<details>
<summary>Image 8 Details</summary>
### Visual Description
## Chart: Energy to Solution vs. Saturation Parameter
### Overview
The image presents two line charts comparing the energy to solution for two different computational methods, CIM-CFC and CIM-SFC, across varying saturation parameters. The chart on the left represents data for N=100, while the chart on the right represents data for N=800. Both charts have a logarithmic scale on both axes.
### Components/Axes
* **Titles:**
* Left Chart: "Energy To Solution (N=100)"
* Right Chart: "Energy To Solution (N=800)"
* **Y-axis (both charts):**
* Label: "OPO Energy To Solution (J)"
* Scale: Logarithmic, ranging from 10<sup>-11</sup> to 10<sup>-4</sup> on the left chart, and from 10<sup>-8</sup> to 10<sup>-2</sup> on the right chart.
* **X-axis (both charts):**
* Label: "σ<sup>2</sup> (saturation parameter)"
* Scale: Logarithmic, ranging from 10<sup>-8</sup> to 10<sup>-1</sup>.
* **Legend (both charts, located in the top-center):**
* Green line with circular markers: "CIM-CFC"
* Purple line with circular markers: "CIM-SFC"
### Detailed Analysis
**Left Chart (N=100):**
* **CIM-CFC (Green):**
* Trend: Initially decreases sharply, then flattens out, and finally increases at higher saturation parameter values.
* Data Points:
* At σ<sup>2</sup> ≈ 10<sup>-8</sup>, Energy ≈ 10<sup>-4</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-4</sup>, Energy ≈ 10<sup>-9</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-2</sup>, Energy ≈ 10<sup>-8</sup> J
* **CIM-SFC (Purple):**
* Trend: Decreases sharply, then flattens out at higher saturation parameter values.
* Data Points:
* At σ<sup>2</sup> ≈ 10<sup>-8</sup>, Energy ≈ 10<sup>-4</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-4</sup>, Energy ≈ 10<sup>-9</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-2</sup>, Energy ≈ 10<sup>-10</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-1</sup>, Energy ≈ 10<sup>-10</sup> J
**Right Chart (N=800):**
* **CIM-CFC (Green):**
* Trend: Decreases sharply, then flattens out at higher saturation parameter values.
* Data Points:
* At σ<sup>2</sup> ≈ 10<sup>-8</sup>, Energy ≈ 10<sup>-2.5</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-4</sup>, Energy ≈ 10<sup>-5</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-2</sup>, Energy ≈ 10<sup>-6</sup> J
* **CIM-SFC (Purple):**
* Trend: Decreases sharply, then increases sharply at higher saturation parameter values.
* Data Points:
* At σ<sup>2</sup> ≈ 10<sup>-8</sup>, Energy ≈ 10<sup>-3</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-4</sup>, Energy ≈ 10<sup>-5.5</sup> J
* At σ<sup>2</sup> ≈ 10<sup>-2</sup>, Energy ≈ 10<sup>-7</sup> J
### Key Observations
* For both N=100 and N=800, both CIM-CFC and CIM-SFC initially show a significant decrease in energy to solution as the saturation parameter increases.
* At N=100, CIM-SFC consistently requires less energy than CIM-CFC across the range of saturation parameters.
* At N=800, CIM-SFC requires less energy than CIM-CFC until σ<sup>2</sup> ≈ 10<sup>-2</sup>, where the energy required by CIM-SFC increases sharply.
* The energy scales are different between the two charts, with the N=800 chart showing higher energy values overall.
### Interpretation
The charts illustrate the relationship between the saturation parameter (σ<sup>2</sup>) and the energy required to reach a solution for two different computational methods (CIM-CFC and CIM-SFC) at two different problem sizes (N=100 and N=800). The initial decrease in energy with increasing saturation parameter suggests that a certain level of saturation is beneficial for both methods. However, at higher saturation levels, the behavior diverges. For N=100, CIM-SFC consistently outperforms CIM-CFC. For N=800, CIM-SFC initially outperforms CIM-CFC, but its performance degrades significantly at higher saturation levels. This suggests that the optimal saturation parameter and the relative performance of the two methods are dependent on the problem size (N). The sharp increase in energy for CIM-SFC at N=800 and high saturation levels could indicate a point where the method becomes unstable or inefficient.
</details>
In Figure 8, we can see that CIM-SFC is more robust to quantum noise compared to CIM-CFC, allowing us to potentially use a larger value of g 2 . This is to be expected owing to the different roles payed by the error variable e i in each system. In CIM-CFC, the feedback signal is calculated as
$$\tilde { z } _ { i } = \tilde { e } _ { i } \sum _ { j } \xi J _ { i j } \tilde { x } _ { j } ( t )$$
which is the main cause of performance degradation when the quantum noise is increased. This is because, if the coherent excitation of ˜ e i is large, then small errors in ∑ j J ij ξ ˜ x j ( t ) will be amplified, and
Table 1: Operational power of active photonic devices in 100 GHz CIM.
| Devices | Power consumption | Reference |
|-------------------------------------------|---------------------|-------------|
| Soliton frequency comb generator | 100 mW | [16] |
| Phase sensitive ampifier (PSA) 10 dB gain | 10 mW | [15] |
| Phase sensitive ampifier (PSA) 50 dB gain | 100 mW | [15] |
| EOM modulator | 400 mW | [17] |
conversely, if the coherent excitation of ∑ j ξJ ij ˜ x j ( t ) is large, then small errors in ˜ e i will be amplified. There are no such beat noise components in CIM-SFC. Therefore, CIM-SFC is more robust to quantum noise. Moreover, the nonlinear function tanh( c ˜ z i ) can help to suppress the quantum noise.
Although they are not shown, the results for CIM-CAC are nearly identical to those for CIM-CFC.
If we include the energy cost in the optical error correction circuit and pump pulse factory (as described in Figure 6), the energy cost is increased by several orders of magnitude, as shown in Figure 9. Here, we assume that the pump pulse energy for a small signal amplification ( ∼ 10 dB) in PSA 1 , PSA 2 , ... , PSA N and PSA e in the optical error correction circuit is 100 fJ/pulse, and that for a large signal amplification ( ∼ 50 dB) in PSA 0 is 1 pJ/pulse. These numbers correspond to the experimental values for a thin-film LiNO 3 ridge waveguide DOPO at a pump wavelength of 780 nm and a pump pulse duration of 100 fs[15]. The pump energy consumed in the optical error correction circuit is estimated as E correction = [( N +1) × 10 -13 +10 -12 ] N (MVM)( J ). The energy consumption in the pump pulse factory is attributed to three components: those of a 100-GHz soliton frequency comb generator, EOM modulators, and phasesensitive amplifiers (Figure 6(b)). The 100-GHz soliton frequency comb generator requires an input power of ∼ 100 mW. [ 16 ] The 100-GHz EOM modulators require an electrical input power of ∼ 400 mW each. [ 17 ] The energy cost per pulse for PSA p is ∼ 1 pJ, while those for N PSA s for the storage ring cavities are ∼ 100 fJ each. Note that N EOMs (EOM 1 , EOM 2 , . . . EOM N ) need to be operated only for one initial round-trip time, 10 -11 N (sec). The operational powers of the active devices in the 100-GHz CIM are summarized in Table 1. The energy cost in the pump pulse factory is E factory = [1 . 3 × 10 -11 N (MVM) + 4 × 10 -Table 2 summarizes the energy costs in three parts of the CIM.
Figure 9: Estimated energy cost to solution of optical and GPU implementations of CIM-CAC vs. problem size √ N . The energy cost to solution for the optical CIM is based on the results presented in Table 2. The energy cost of the GPU is based on the ˜ 200W power consumption of the Nvidia Tesla V100 GPU used.
<details>
<summary>Image 9 Details</summary>
### Visual Description
## Chart: Energy to Solution vs. Problem Size
### Overview
The image is a line chart comparing the energy required to solve a problem using two different implementations: an "Optical Implementation" and an "Nvidia Tesla V100". The chart plots energy (in Joules) on a logarithmic y-axis against the square root of the problem size on the x-axis.
### Components/Axes
* **Title:** None explicitly present in the image.
* **X-axis:** "problem size" with markers at √100, √200, √400, √500, √800, and √1000.
* **Y-axis:** "Energy to Solution (J)" with a logarithmic scale ranging from 10^-5 to 10^3.
* **Legend:** Located at the top-left of the chart.
* **Optical Implementation:** Represented by an orange line.
* **Nvidia Tesla V100:** Represented by a blue line.
### Detailed Analysis
* **Optical Implementation (Orange Line):** The energy consumption increases as the problem size increases.
* √100: Approximately 10^-5 J
* √200: Approximately 3 * 10^-4 J
* √400: Approximately 2 * 10^-3 J
* √500: Approximately 4 * 10^-3 J
* √800: Approximately 2 * 10^-2 J
* √1000: Approximately 2 * 10^-1 J
* **Nvidia Tesla V100 (Blue Line):** The energy consumption also increases as the problem size increases, but at a much higher rate than the optical implementation.
* √100: Approximately 2 * 10^-2 J
* √200: Approximately 3 * 10^-1 J
* √400: Approximately 6 J
* √500: Approximately 9 J
* √800: Approximately 60 J
* √1000: Approximately 150 J
### Key Observations
* The Nvidia Tesla V100 requires significantly more energy to solve the problem compared to the Optical Implementation across all problem sizes.
* The energy consumption of the Nvidia Tesla V100 increases more rapidly with problem size than the Optical Implementation.
* Both implementations show an increasing trend in energy consumption as the problem size increases.
### Interpretation
The chart demonstrates that for the problem being considered, the Optical Implementation is significantly more energy-efficient than the Nvidia Tesla V100. The difference in energy consumption becomes more pronounced as the problem size increases, suggesting that the Optical Implementation scales better in terms of energy efficiency. This implies that for larger problem sizes, the Optical Implementation would be the preferred choice if energy consumption is a critical factor. The logarithmic scale on the y-axis emphasizes the exponential difference in energy consumption between the two implementations.
</details>
Figure 9 shows the energy cost to solution if the CIM-CFC algorithm is implemented on GPU. The detailed description of this approach will be given in the next section. Even though the optical implementation of the error correction circuit and pump pulse factory as described in Figure 6 is technologically challenging, the energy cost can be decreased by several orders of magnitude compared to a modern GPU.
Table 2: Energy cost to solution in three subsystems in CIM. MVM: matrix-vector multiplication steps to solution, N : problem size, one round trip time: 10 -8 s, signal lifetime: 10 -7 s.
| Subsystem | Energy-to-solution |
|----------------------------------|----------------------------------------------------------------------------------------------------------------|
| Main cavity | 2 . 6 × 10 - 20 (MVM) N/g 2 |
| Optical error correction circuit | [( N +1)10 - 13 +10 - 12 ](MVM) N 10 - 13 N 2 (MVM) |
| Pump pulse factory | [ N × 10 - 13 +10 - 12 ](MVM) N +1 . 3 × 10 - 11 N (MVM)+4 × 10 - 12 N 2 10 - 13 N 2 (MVM) |
## 6 CIM - Inspired Heuristic Algorithms
## 6.1 Scaling performance of CIM-CAC, CIM-CFC, and CIM-SFC
To test whether the three classical nonlinear dynamics models given by Eqs. (1)-(8), are good Ising solvers, we can numerically integrate them on a digital platform. In this section, we will consider these CIMinspired algorithms when numerically integrated using an Euler step. In addition, to ensure numerical stability, we constrain the range of some variables, the details of which are presented in Appendix A.
The relevant performance metric is the time to solution or TTS (the number of integration time steps required to achieve a success rate of 99%). In particular, we study how the median TTS scales as a function of the problem size for randomly generated SK spin glass instances (the couplings are chosen randomly between +1 and -1). The median TTS is computed on the basis of a set of 100 randomly generated instances per problem size, and 3200 trajectories are used per instance to evaluate the TTS.
Figure 10 shows the median TTSs of the three CIM-inspired algorithms (CIM-CAC, CIM-CFC, and CIMSFC) are shown with respect to the problem size. The shaded regions represent 25th-75th percentiles. The linear behavior of the TTS with respect to √ N indicates that these algorithms have the same root exponential scaling of TTS that is also observed in physical CIMs with quantum noise from external reservoirs. [ 8 , 31 ] All three algorithms appear to have very similar scaling coefficients if the TTS is assumed to be of the form TTS ≈ A · B √ n . In addition to the similar scaling, all three algorithms show a similar spread (25th-75th percentile) in TTS, as indicated by the shaded region above. Although CIM-SFC may have a slightly larger spread in all cases, the spread does not appear to increase for larger problem sizes.
Figure 10: (left) TTSs of CIM-CAC, CIM-CFC and CIM-SFC vs. problem size √ N . The shaded regions represent the 25th-75th percentile TTS. (right) The 75th and 90th percentile TTS compared to the median TTS for all three systems as a function of the problem size.
<details>
<summary>Image 10 Details</summary>
### Visual Description
## Combined Line Chart: Algorithm Performance vs. Problem Size
### Overview
The image presents two line charts side-by-side, comparing the performance of different algorithms (CIM-CAC, CIM-CFC, CIM-SFC) against varying problem sizes. The left chart displays "time steps to solution" on a logarithmic scale, while the right chart shows the "ratio to median TTS" (Time To Solution), also on a logarithmic scale. The x-axis for both charts represents "problem size," expressed as the square root of the actual problem size.
### Components/Axes
**Left Chart:**
* **Title:** time steps to solution
* **Y-axis:** time steps to solution (logarithmic scale)
* Axis markers: 10^3, 10^4, 10^5, 10^6, 10^7
* **X-axis:** problem size (square root)
* Axis markers: √100, √200, √400, √500, √800, √1000
* **Legend (top-left):**
* Blue: CIM-CAC
* Green: CIM-CFC
* Red: CIM-SFC
* Shaded regions around each line represent some measure of variance or uncertainty.
**Right Chart:**
* **Title:** ratio to median TTS
* **Y-axis:** ratio to median TTS (logarithmic scale)
* Axis markers: 10^0, 10^1, 10^2
* **X-axis:** problem size (square root)
* Axis markers: √100, √200, √400, √500, √800, √1000
* **Legend (top-left):**
* Black (solid line): 75th percentile
* Black (dashed line): 90th percentile
### Detailed Analysis
**Left Chart (Time Steps to Solution):**
* **CIM-CAC (Blue):** The line slopes upward, indicating that the time steps to solution increase with problem size.
* √100: ~3000
* √200: ~20000
* √400: ~45000
* √500: ~60000
* √800: ~200000
* √1000: ~400000
* **CIM-CFC (Green):** The line slopes upward, indicating that the time steps to solution increase with problem size.
* √100: ~2500
* √200: ~10000
* √400: ~30000
* √500: ~40000
* √800: ~100000
* √1000: ~250000
* **CIM-SFC (Red):** The line slopes upward, indicating that the time steps to solution increase with problem size.
* √100: ~2000
* √200: ~8000
* √400: ~25000
* √500: ~35000
* √800: ~80000
* √1000: ~200000
**Right Chart (Ratio to Median TTS):**
* **75th percentile (Black, solid):** The line is relatively flat, suggesting the 75th percentile ratio to median TTS remains relatively constant across problem sizes.
* √100: ~1
* √200: ~1.5
* √400: ~1.5
* √500: ~1.3
* √800: ~1.5
* √1000: ~1.5
* **90th percentile (Black, dashed):** The line fluctuates more than the 75th percentile, but generally stays within a small range.
* √100: ~2.5
* √200: ~2
* √400: ~2.5
* √500: ~3
* √800: ~2
* √1000: ~2.5
### Key Observations
* In the left chart, all three algorithms (CIM-CAC, CIM-CFC, CIM-SFC) exhibit an increase in "time steps to solution" as the problem size increases.
* CIM-CAC consistently requires more time steps to solution than CIM-CFC and CIM-SFC.
* In the right chart, the 75th and 90th percentile ratios to median TTS are relatively stable across different problem sizes.
* The ratio to median TTS for CIM-SFC is higher than CIM-CAC and CIM-CFC at smaller problem sizes, but converges as the problem size increases.
### Interpretation
The left chart demonstrates the expected behavior that larger problems require more computational steps to solve. The relative positioning of the lines suggests that CIM-CAC is generally less efficient in terms of time steps than CIM-CFC and CIM-SFC for the problem set being analyzed.
The right chart provides insight into the variability of the algorithms' performance. The relatively stable 75th and 90th percentile lines indicate that the distribution of Time To Solution (TTS) is consistent across different problem sizes. The higher ratio to median TTS for CIM-SFC at smaller problem sizes suggests that it may have more variable performance in those scenarios, potentially due to factors like initialization or search space exploration. As the problem size increases, the algorithms' performance becomes more consistent relative to their median TTS.
</details>
## 6.2 Comparison with noisy mean field annealing (NMFA)
To show the importance of the auxiliary variable (error pulse) in CIM-SFC, we compared its performance to another CIM-inspired algorithm, namely noisy mean field annealing (NMFA). [ 26 ] NMFA also applies a hyperbolic tangent function to the mutual coupling term. However, it does not have an auxiliary variable and relies on (artificial) quantum noise to escape from local minima. Figure 11, compares the scal-
ing of NMFA to CIM-SFC with different values of the feedback parameter k . As k controls the strength of the destabilization force caused by the auxiliary variable, we can measure the importance of the term k ( z i -e i ) to the scaling behavior. When k = 0, CIM-SFC is nearly identical to NMFA. The fact that CIM-SFC with k = 0 shows slightly worse performance indicates that the noise included in NMFA likely has a small effect and may help destabilize the local minima (which can also be observed in Figure 7). The case k = 0 . 15 is shown as an intermediate case, and k = 0 . 2 is the (experimentally obtained) optimal value for k in CIM-SFC.
As can be seen, the addition of the error correction feedback term k ( z i -e i ) in Eq. (7) is effective in improving both the scaling and the spread of TTS for the SK instances. This implies that the 'correlated artificial noise' provided by the auxiliary variable is more effective in finding better solutions than the 'random quantum noise' from reservoirs.
Figure 11: (left) TTSs of CIM-SFC and NMFA vs. problem size √ N . The shaded regions represent 25-75 percentile TTS. The results for NMFA are from [12]. (right) TTS for N=400 SK instances for different values of k . The 75th and 90th percentiles are not shown for smaller values of k because they were too large to be computed.
<details>
<summary>Image 11 Details</summary>
### Visual Description
## Chart: Time Steps to Solution vs. Problem Size and k-value
### Overview
The image presents two line charts comparing the performance of different algorithms in terms of "time steps to solution." The left chart shows the relationship between "time steps to solution" and "problem size" for CIM-SFC with varying k-values (0.2, 0.15, 0.0) and NMFA. The right chart shows the relationship between "time steps to solution" and "k" for the 90th, 75th, median, and 25th percentiles. Both charts use a logarithmic scale for the y-axis ("time steps to solution").
### Components/Axes
**Left Chart:**
* **X-axis:** "problem size" with markers at approximately √100, √200, √400, √500, √800, and √1000.
* **Y-axis:** "time steps to solution" (logarithmic scale) with markers at 10^3, 10^4, 10^5, 10^6, 10^7, and 10^8.
* **Legend (top-left):**
* Red: CIM-SFC k=0.2
* Purple: CIM-SFC k=0.15
* Blue: CIM-SFC k=0.0
* Black: NMFA
**Right Chart:**
* **X-axis:** "k" with markers at 0.00, 0.05, 0.10, 0.15, 0.20, and 0.25.
* **Y-axis:** "time steps to solution" (logarithmic scale) with markers at 10^4, 10^5, 10^6, 10^7, and 10^8.
* **Legend (top-left):**
* Red dotted: 90th pct
* Red dashed: 75th pct
* Red solid: median
* Red dotted: 25th pct
### Detailed Analysis
**Left Chart:**
* **CIM-SFC k=0.2 (Red):** The line slopes upward.
* √100: ~4 x 10^3
* √200: ~9 x 10^3
* √400: ~3 x 10^4
* √500: ~5 x 10^4
* √800: ~1.5 x 10^5
* √1000: ~2 x 10^5
* **CIM-SFC k=0.15 (Purple):** The line slopes upward.
* √100: ~4 x 10^3
* √200: ~1.5 x 10^4
* √400: ~8 x 10^4
* √500: ~1 x 10^5
* √800: ~5 x 10^6
* √1000: ~8 x 10^6
* **CIM-SFC k=0.0 (Blue):** The line slopes upward.
* √100: ~4 x 10^3
* √200: ~2 x 10^4
* √400: ~1 x 10^5
* √500: ~2 x 10^5
* **NMFA (Black, dotted):** The line slopes upward.
* √100: ~4 x 10^3
* √200: ~2.5 x 10^4
* √400: ~2 x 10^5
* √500: ~4 x 10^5
* √800: ~1 x 10^7
**Right Chart:**
* **90th pct (Red, dotted):** The line decreases until k=0.2, then increases.
* 0.00: ~5 x 10^7
* 0.05: ~4 x 10^7
* 0.10: ~2 x 10^7
* 0.15: ~4 x 10^6
* 0.20: ~2 x 10^6
* 0.25: ~4 x 10^6
* **75th pct (Red, dashed):** The line decreases until k=0.2, then increases.
* 0.00: ~5 x 10^6
* 0.05: ~4 x 10^6
* 0.10: ~2 x 10^6
* 0.15: ~4 x 10^5
* 0.20: ~2 x 10^5
* 0.25: ~4 x 10^5
* **Median (Red, solid):** The line decreases until k=0.2, then increases.
* 0.00: ~5 x 10^5
* 0.05: ~4 x 10^5
* 0.10: ~2 x 10^5
* 0.15: ~4 x 10^4
* 0.20: ~2 x 10^4
* 0.25: ~4 x 10^4
* **25th pct (Red, dotted):** The line decreases until k=0.2, then increases.
* 0.00: ~2.5 x 10^4
* 0.05: ~2 x 10^4
* 0.10: ~1 x 10^4
* 0.15: ~2 x 10^3
* 0.20: ~1 x 10^3
* 0.25: ~2 x 10^3
### Key Observations
* In the left chart, for a given problem size, NMFA generally requires more time steps to solution than CIM-SFC with different k values.
* In the left chart, as the problem size increases, the time steps to solution also increase for all algorithms.
* In the right chart, the time steps to solution generally decrease as 'k' increases from 0.00 to 0.20, and then increase as 'k' increases from 0.20 to 0.25 for all percentiles.
* The right chart shows a clear trend of decreasing time steps to solution with increasing 'k' up to a point, suggesting an optimal 'k' value around 0.20.
### Interpretation
The charts suggest that CIM-SFC algorithms are more efficient than NMFA for the tested problem sizes. The optimal 'k' value for minimizing the time steps to solution appears to be around 0.20. The performance of all algorithms degrades as the problem size increases, which is expected. The percentile data in the right chart indicates the distribution of time steps to solution for different 'k' values, with the median representing the typical performance and the 25th and 75th percentiles indicating the range of variability. The 90th percentile shows the worst-case performance.
</details>
## 6.3 Comparison with discrete simulated bifurcation machine (dSBM)
We compared the performance of the CIM-inspired algorithms with that of another heuristic Ising solver, namely the discrete simulated bifurcation machine (dSBM). [ 27 , 29 , 28 ] Similar to CIM, dSBM also makes use of analog spins and continuous dynamics to solve combinatorial optimization problems.
We are aware that the authors of [29] seem to claim that dSBM is algorithmically superior to CIM-CAC by comparing the required number of MVMs to solution. Although the authors of [29] discussed the wall clock TTS of their implementations on many problem sets, when making the claim of algorithmic superiority, they only used the median TTS (in units of MVM) on SK instances for two problem sizes. In this section, we will provide a more detailed comparison of the three algorithms (CIM-CAC, CIM-CFC and CIM-SFC) with dSBM using MVM to solution (or equivalently, integration time steps to solution) as the performance metric. As mentioned before, this is a good comparison because all these algorithms will have MVM as the computational bottleneck when implemented on a digital platform. As discussed in Section 4, the computation of the Ising energy can be left until the end of the trajectory in most cases; thus, for this section, we will only consider the MVM involved in the computation of the mutual coupling term when calculating the MVM to solution.
The problem instance sets used in this section are:
1. A set of 100 randomly generated 800-spin SK instances (available upon request from the authors). This instance set contains fully connected instances with weights of +1 , -1.
2. The G-set instances that have been used as a benchmark for max-cut performance (available at https://web.stanford.edu/ yyye/yyye/Gset/). In this study, we consider 50 instances with a problem size of 800-2000. These instances have varying edge density and include weights of either +1 , 0 or +1 , 0 , -1.
3. Another set of 1000 randomly generated 800-spin and 1200-spin SK instance (available upon request from the authors) used to evaluate the worst-case performance.
To compare the performance on the 800-spin SK instances, the dSBM algorithm was also implemented on GPU. The parameters for dSBM were chosen on the basis of the parameters in [29] (see Appendix D).
Figure 12: Required number of MVMs to solution for CIM-CAC, CIM-CFC, and CIM-SFC vs. that for dSBM. The median TTS is indicated by the red lines and the 25th-75th percentiles are indicated by the shaded blue regions.
<details>
<summary>Image 12 Details</summary>
### Visual Description
## Scatter Plot Comparison: CIM vs dSBM
### Overview
The image presents three scatter plots comparing the performance of different CIM (Combinatorial Interaction Modeling) methods against a dSBM (degree-corrected stochastic block model) approach. Each plot compares a specific CIM variant (CAC, CFC, SFC) against dSBM in terms of "MVM to solution," likely referring to a measure of computational effort or resources required to reach a solution. The plots use a log-log scale, and each includes a dashed diagonal line, a horizontal line, and a vertical line, along with shaded regions.
### Components/Axes
* **Titles:**
* Left Plot: "CIM-CAC vs dSBM"
* Middle Plot: "CIM-CFC vs dSBM"
* Right Plot: "CIM-SFC vs dSBM"
* **X-axis (all plots):**
* Label: "CIM-[Variant] MVM to solution" (where [Variant] is CAC, CFC, or SFC)
* Scale: Logarithmic, ranging from 10^4 to 10^7
* Ticks: 10^4, 10^5, 10^6, 10^7
* **Y-axis (all plots):**
* Label: "dSBM MVM to solution"
* Scale: Logarithmic, ranging from 10^4 to 10^7
* Ticks: 10^4, 10^5, 10^6, 10^7
* **Data Points:** Blue dots representing individual data points.
* **Diagonal Line:** Dashed gray line representing the line of equality (x=y).
* **Horizontal Line:** Solid red line, indicating a specific dSBM MVM value. Appears to be around 2 * 10^5.
* **Vertical Line:** Solid red line, indicating a specific CIM MVM value. Appears to be around 2 * 10^5.
* **Shaded Regions:** Light blue shaded regions around the horizontal and vertical lines, indicating a range or tolerance around those values.
### Detailed Analysis
**Plot 1: CIM-CAC vs dSBM**
* **X-axis:** CIM-CAC MVM to solution
* **Y-axis:** dSBM MVM to solution
* **Data Points:** The blue data points are clustered around the intersection of the red lines, with a spread both above and below the diagonal line.
* **Trend:** The data points generally follow a positive correlation, but with significant scatter.
* **Horizontal Line:** Located at approximately 2 * 10^5 on the Y-axis.
* **Vertical Line:** Located at approximately 2 * 10^5 on the X-axis.
**Plot 2: CIM-CFC vs dSBM**
* **X-axis:** CIM-CFC MVM to solution
* **Y-axis:** dSBM MVM to solution
* **Data Points:** Similar to the first plot, the data points are clustered around the intersection of the red lines, with a spread both above and below the diagonal line.
* **Trend:** The data points generally follow a positive correlation, but with significant scatter.
* **Horizontal Line:** Located at approximately 2 * 10^5 on the Y-axis.
* **Vertical Line:** Located at approximately 2 * 10^5 on the X-axis.
**Plot 3: CIM-SFC vs dSBM**
* **X-axis:** CIM-SFC MVM to solution
* **Y-axis:** dSBM MVM to solution
* **Data Points:** Similar to the first two plots, the data points are clustered around the intersection of the red lines, with a spread both above and below the diagonal line.
* **Trend:** The data points generally follow a positive correlation, but with significant scatter.
* **Horizontal Line:** Located at approximately 2 * 10^5 on the Y-axis.
* **Vertical Line:** Located at approximately 2 * 10^5 on the X-axis.
### Key Observations
* All three plots show a similar pattern: a cluster of data points around the intersection of the red lines, with a general positive correlation between the CIM variant and dSBM.
* The scatter around the diagonal line suggests that the performance of the CIM variants and dSBM can vary significantly for individual instances.
* The red lines and shaded regions appear to highlight a specific performance range or target.
### Interpretation
The plots compare the "MVM to solution" (likely a measure of computational effort) for different CIM methods (CAC, CFC, SFC) against dSBM. The clustering of data points around the intersection of the red lines suggests that, on average, the CIM variants and dSBM have similar performance characteristics within the range defined by the red lines and shaded regions. However, the significant scatter indicates that the relative performance of these methods can vary considerably depending on the specific problem instance. The diagonal line represents the ideal scenario where both methods have equal "MVM to solution." Points above the line indicate that dSBM requires more effort, while points below the line indicate that the CIM variant requires more effort. The shaded regions likely represent a tolerance or acceptable range of performance.
</details>
In Figure 12compares the performance of the three algorithms (CIM-CAC, CIM-CFC, and CIM-SFC) instance by instance with that of dSBM on the 800-spin SK instance set. The ground-state energies used to evaluate the MVM to solution are the lowest energies found by the four algorithms. As all four algorithms found the same lowest energies, it is highly likely that these are true ground-state energies. The parameters for all four systems can be found in Appendix A. As shown in Figure 12, all four systems showed remarkably similar performance on the 800-spin instances when the parameters were optimized. It is important to note that with the parameters used in Figure 12, CIM-SFC did not find the ground state in one instance. However, if different parameters are used, CIM-SFC will find the ground state for this particular instance as well. Thus, although CIM-SFC can achieve high performance, it is highly sensitive to parameter selection.
The median TTS (in the units of MVM) of CIM-CFC, CIM-SFC and dSBM are nearly the same: around 2 × 10 5 . Furthermore, the spread in TTS of these three algorithms is similar. Although CIM-CAC shows slightly worse median TTS (by less than a factor of two), it is worth noting that the instances in which CIM-CAC performs better than dSBM tend to be the harder instances. This indicates that among the four algorithms, CIM-CAC may show slightly better worst-case performance. We investigate this further later in this section. This pattern can also be seen on the G-set.
Overall, all four algorithms show similar performance on the fully connected instance set, and it is impossible to determine which particular algorithm is the most effective one for this problem type. In addition to the similar median TTS and spread, there is a high level of correlation in the TTS among all four systems. This indicates either that instance difficulty is a universal property for all Ising heuristics or that there is something fundamentally common to the four algorithms. See Appendix E for a further discussion of the similarities and differences between these four systems.
Although CIM-SFC shows good performance on fully connected problem instances, it struggles on many G-set instances. In Appendix D, we discuss a partial reason for this failure; however, the full reason is yet to be understood. In the future, we expect to modify CIM-SFC or find better parameters so that it can solve all problem types; however, for now, we will just consider CIM-SFC as a fully connected (or densely connected) Ising solver and compare only CIM-CAC, CIM-CFC, and dSBM on the G-set. For some results of CIM- SFC on the G-set, see Appendix D.
All three algorithms (CIM-CAC, CIM-CFC, and dSBM) show fairly good performance on the G-set; however, we argue that CIM-CAC is the most consistently effective algorithm. CIM-CAC and dSBM
Q
Figure 13: (left) TTS on G-set graphs for CIM-CAC, CIM-CFC and dSBM. The best known cut values are found in [29]. The TTS for dSBM is from [29]. In this plot, the instances are separated into groups depending on the graph type and size. The dots above the plot indicate that the best known cut value was not found. (right) Histograms showing which algorithm realizes the slowest TTS and which algorithm has the fastest TTS. The column labeled 'none' indicates the two instances in which none of the three algorithms found the best known cut value. The parameters are chosen and optimized separately for each instance type. An instance-by-instance comparison and the parameters used are presented in Appendix B.
<details>
<summary>Image 13 Details</summary>
### Visual Description
## Algorithm Performance Comparison
### Overview
The image presents a comparison of three algorithms (CIM-CAC, CIM-CFC, and dSBM) based on their performance in solving instances of varying sizes (N = 800, N = 1000, N = 2000) and types (Random, Toroidal, Planar). The performance is measured by "MVM to solution" on a logarithmic scale. Additionally, the image includes bar charts indicating which algorithm is the "slowest" and "fastest" based on the number of instances.
### Components/Axes
**Main Scatter Plot:**
* **Y-axis:** "MVM to solution" (Matrix-Vector Multiplication to solution), logarithmic scale from 10^4 to 10^12.
* **X-axis:** Categorical axis representing different instance groups, divided by N values (800, 1000, 2000) and instance types (R, T, P) with corresponding values {1} or {1,-1}.
* N = 800:
* 1-5: R {1}
* 6-10: R {1,-1}
* 11-13: T {1,-1}
* 14-17: P {1}
* 18-21: P {1,-1}
* N = 1000:
* 42-47: R {1}
* 51-54: P {1}
* N = 2000:
* 22-26: R {1}
* 27-31: R {1,-1}
* 32-34: T {1,-1}
* 35-38: P {1}
* 39-42: P {1,-1}
* **Legend (top-left):**
* Blue: CIM-CAC
* Green: CIM-CFC
* Red: dSBM
**Bar Charts (top-right and bottom-right):**
* **Y-axis:** "Number of instances" (linear scale).
* **X-axis:** Algorithm names (CAC, CFC, dSBM, none).
* **Top Bar Chart:** "Which algorithm is slowest?"
* **Bottom Bar Chart:** "Which algorithm is fastest?"
**Legend (bottom):**
* R = Random
* T = Toroidal
* P = Planar
### Detailed Analysis
**Main Scatter Plot:**
* **CIM-CAC (Blue):** Generally shows lower MVM values compared to dSBM, especially for N=800 and N=1000. The values increase as N increases to 2000.
* **CIM-CFC (Green):** Similar to CIM-CAC, with relatively lower MVM values, but with some instances showing higher values. The values increase as N increases to 2000.
* **dSBM (Red):** Shows a wider range of MVM values, with many instances having significantly higher values than CIM-CAC and CIM-CFC, especially for N=2000.
**Specific Data Points (Approximate):**
* For N=800, R {1}, MVM values range approximately from 10^5 to 10^7 for all algorithms.
* For N=2000, P {1,-1}, dSBM reaches values close to 10^12, while CIM-CAC and CIM-CFC are around 10^7 to 10^9.
**Bar Charts:**
* **Slowest Algorithm:**
* CAC: ~7 instances
* CFC: ~6 instances
* dSBM: ~33 instances
* None: ~1 instance
* **Fastest Algorithm:**
* CAC: ~19 instances
* CFC: ~22 instances
* dSBM: ~7 instances
* None: ~2 instances
### Key Observations
* dSBM tends to require a higher number of MVM operations to reach a solution compared to CIM-CAC and CIM-CFC, especially for larger problem sizes (N=2000).
* CIM-CAC and CIM-CFC show more consistent and lower MVM values across different instance types and sizes.
* According to the bar charts, dSBM is most frequently the slowest algorithm, while CIM-CFC is most frequently the fastest.
### Interpretation
The data suggests that CIM-CAC and CIM-CFC are generally more efficient algorithms for solving these types of instances, as they require fewer matrix-vector multiplications to reach a solution. dSBM, while sometimes faster, is more often the slowest and exhibits a wider range of performance, indicating it may be more sensitive to the specific characteristics of the problem instance. The bar charts reinforce this conclusion, showing that dSBM is most often the slowest, while CIM-CFC is most often the fastest. The "none" category in both bar charts suggests that in some instances, none of the tested algorithms were particularly fast or slow compared to some baseline or expectation. The increase in MVM values as N increases to 2000 is expected, as larger problem sizes generally require more computational effort.
</details>
were able to find the best known cut values in 47 out of 50 instances, while CIM-CFC found the best known cut value in 45 out of 50 instances. It is worth noting that the simulation time used to calculate the TTS for dSBM [ 29 ] was much longer than that used in this study. Given the same simulation time, dSBM would most likely have solved only 45 out of 50 instances. As shown in Figure 13, CIM-CAC and CIM-CFC are faster (in units of MVM) than dSBM in most instances. More importantly, among the instances in which dSBM is faster, there are no cases where dSBM is significantly faster than CIM-CAC, other than G37, in which CIM-CAC did not find the best known cut value. Meanwhile, we found that CIM-CAC was more than an order of magnitude faster than dSBM in 13 out 50 instances. Therefore, we believe that CIM-CAC is a more reliable algorithm when considering many problem types.
The difference between CIM-CAC and CIM-CFC is subtle. This is to be expected, as the dynamics of the two systems are very similar. Although the performance of the two algorithms for G-set is nearly identical in most cases, for some of the harder instances, there are some cases in which CIM-CFC cannot find the best known cut value or CIM-CFC has a significantly longer TTS. This indicates that CIMCAC is fundamentally a more promising algorithm or that precise parameter selection for CIM-CFC is required.
As noted in Figure 12, the worst case performance of CIM-CAC may be slightly better than that of dSBM. To evaluate this further we created new sets of 1000 800-spin and 1200-spin SK instances. Figure 14 shows the number of instances solved as a function of the number of MVMs required to achieve a success probability of 99%. As can be seen in both cases, dSBM can solve the easier instances with fewer MVMs; however, for the hardest instances, CIM-CAC is faster. his can be understood by observing the intersection point of the two curves in Figure 14.
In nearly all cases, the best Ising energy found was the same for both the solvers when a similar number of MVMs were used (see Appendix A for the parameters). However, for two instances in the N=1200 set, the Ising energy found by CAC was not found by dSBM. This remained true even when 50,000 dSBM trajectories were used for these instances.
Our results suggest that dSBM may struggle considerably for some harder SK instances. However, we acknowledge that this could be a result of sub-optimal parameter selection for dSBM. The parameters used (see appendix A) were optimized manually to achieve a good median TTS; however, they may not
<details>
<summary>Image 14 Details</summary>
### Visual Description
## Chart Type: Comparative Line Graphs
### Overview
The image presents two line graphs comparing the performance of two algorithms, CIM-CAC and dSBM, in solving a set of instances. The graphs show the number of instances unsolved out of 1000 SK instances as a function of the number of MVMs (Million Vector Multiplications) required. The left graph corresponds to a dataset with N=800, while the right graph corresponds to N=1200. Both graphs use a log-log scale.
### Components/Axes
* **X-axis (horizontal):** "# MVMs" (Number of Million Vector Multiplications). Logarithmic scale, ranging from 10^3 to 10^8 for N=800, and 10^3 to 10^9 for N=1200.
* **Y-axis (vertical):** "# of instances unsolved out of 1000 SK instances". Logarithmic scale, ranging from 10^0 to 10^3.
* **Titles:** "N=800" (left graph), "N=1200" (right graph).
* **Legend (bottom-left of each graph):**
* Blue line: CIM-CAC
* Red line: dSBM
### Detailed Analysis
**Left Graph (N=800):**
* **CIM-CAC (Blue):** The line starts near 10^3 unsolved instances at 10^3 MVMs. It slopes downward, indicating that as the number of MVMs increases, the number of unsolved instances decreases. The line reaches approximately 1 unsolved instance around 10^7 MVMs.
* (10^3, 10^3)
* (10^4, ~300)
* (10^5, ~100)
* (10^6, ~20)
* (10^7, ~1)
* **dSBM (Red):** Similar to CIM-CAC, the line starts near 10^3 unsolved instances at 10^3 MVMs and slopes downward. However, it appears to decrease more slowly than CIM-CAC initially. The line reaches approximately 1 unsolved instance around 10^7 MVMs.
* (10^3, 10^3)
* (10^4, ~400)
* (10^5, ~100)
* (10^6, ~30)
* (10^7, ~2)
**Right Graph (N=1200):**
* **CIM-CAC (Blue):** The line starts near 10^3 unsolved instances at 10^3 MVMs. It slopes downward, indicating that as the number of MVMs increases, the number of unsolved instances decreases. The line reaches approximately 1 unsolved instance around 10^8 MVMs.
* (10^3, 10^3)
* (10^4, ~400)
* (10^5, ~100)
* (10^6, ~20)
* (10^7, ~2)
* (10^8, ~1)
* **dSBM (Red):** Similar to CIM-CAC, the line starts near 10^3 unsolved instances at 10^3 MVMs and slopes downward. However, it appears to decrease more slowly than CIM-CAC initially. The line reaches approximately 1 unsolved instance around 10^8 MVMs.
* (10^3, 10^3)
* (10^4, ~500)
* (10^5, ~100)
* (10^6, ~30)
* (10^7, ~5)
* (10^8, ~2)
### Key Observations
* Both algorithms, CIM-CAC and dSBM, show a decrease in the number of unsolved instances as the number of MVMs increases.
* For N=800, CIM-CAC appears to perform slightly better than dSBM, reaching lower unsolved instance counts for a given number of MVMs.
* For N=1200, the difference in performance between CIM-CAC and dSBM is less pronounced, especially at higher MVM counts.
* The N=1200 graphs show that more MVMs are required to solve the same number of instances compared to N=800.
### Interpretation
The graphs illustrate the trade-off between computational effort (MVMs) and the number of unsolved instances for two different algorithms. The downward slope of the lines indicates that increasing the computational effort generally leads to solving more instances. The comparison between CIM-CAC and dSBM suggests that CIM-CAC may be more efficient for smaller datasets (N=800), but the performance difference diminishes as the dataset size increases (N=1200). The shift of the curves to the right for N=1200 indicates that larger datasets require more computational resources to achieve the same level of problem-solving. The dashed lines likely represent some form of confidence interval or variance in the data.
</details>
# MVMs
# MVMs
Figure 14: Number of instances that remain unsolved (success probability under 99%) after a certain number of MVMs (time steps) for CIM-CAC and dSBM. The dotted lines represent the assumption of a log-normal distribution for the TTS on randomly generated SK instances. The instance sets used are 1000 randomly generated SK instances (different from the 100 instances used in Figure 12) of problem sizes N=800 (top) and N=1200 (bottom). The parameters for both systems can be found in Appendix A.
be the best parameters if one wants to solve the hardest instances. By contrast, for CIM-CAC, the optimal parameters for the median TTS appear to also perform well on the hardest instances.
To ensure that an Ising solver can find the true ground state of a given problem, the worst-case performance is very important. For this purpose, we believe that CIM-CAC is likely the more fundamentally superior algorithm, at least in the case of randomly generated SK instances. For the other CIM modifications, this is likely not true. In particular, for CIM-SFC, the worst-case performance is significantly worse than that of dSBM and CIM-CAC (as shown in Figure 10). In the future, it would be interesting to investigate the cause for this phenomenon and also to examine the spread of TTS for different problem types.
## 7 Conclusion
The new coherent Ising machines presented in this paper (CIM-CFC and CIM-SFC) have considerable potential as both digital heuristic algorithms and optically implemented physical devices. Rapid advances in the thin-film LiNbO 3 platform [ 15 , 16 , 17 ] as a photonic integrated circuit technology might enable the proposed optical CIM to surpass existing digital algorithms on a CMOS platform in terms of both speed and energy consumption.
The proposed CIM-inspired algorithms were shown to be fast and accurate Ising solvers even when implemented on an existing digital platform. In particular, we showed that their performance is very similar to that of other existing analog-system-based algorithms such as dSBM. This again brings up the question raised in [11] as to whether the simulation of analog spins on a digital computer can outperform a purely discrete heuristic algorithm. Finally, whether chaotic dynamics are necessary for a deterministic dynamical system to be a good Ising solver is left as an open issue for future research. [ 11 , 12 , 29 ]
## Appendix A: Optimization of simulation parameters
Here we summarize the simulation parameters used in our numerical experiments. The parameters are optimized empirically and thus do not necessarily reflect the true optimum values.
## Parameters used in Figure 5
| N step | CIM-SFC (upper left panel) | | CIM-CFC (upper right panel) |
|----------|------------------------------|--------|-------------------------------|
| | 500 | N step | 1000 |
| ∆ T | 0.4 | ∆ T | 0.4 |
| p | -1.0 | p | -1.0 |
| c | 1.0 | α | 1.0 |
| β | 0.3 | β | |
| k | 0.2 | | 0.2 |
| | CIM-SFC (lower left panel) | | CIM-CFC (lower right panel) |
| | | N step | 1000 |
| N step | 500 | ∆ T | 0.4 |
| ∆ T | 0.4 | T | 900 |
| p | -1.0 → 1.0 | r T | 100 |
| c | 1.0 → 3.0 | p | |
| β | 0.3 → 0.1 | p | -1.0 → 1.0 |
| | | α | 1.0 |
| k | 0.2 | β | 0.2 |
## Parameters used in Figure 10, 11, and 12
## CIM-CAC
In our simulation, the x i variables are restricted to the range [ -3 2 √ α, 3 2 √ α ] at each time step. The parameters p and α are modulated linearly from their starting to ending values during the T r time steps and are kept at the final value for an additional T p time steps. The initial value x i is set to a random value chosen from a zero-mean Gaussian distribution with a standard deviation of 10 -4 and e i = 1. Furthermore, 3200 trajectories are computed per instance to evaluate TTS. The actual parameters used for simulation are listed below:
| N step | 3200 |
|----------|------------|
| ∆ T | 0.125 |
| T r | 2880 |
| T p | 320 |
| p | -1.0 → 1.0 |
| α | 1.0 → 2.5 |
| β | 0.8 |
## CIM-CFC
In our simulation the x i variables are restricted to the range [ -1 . 5 , 1 . 5] and e i is restricted to the range [0 . 01 , ∞ ]. The parameter p is modulated linearly from its starting to ending values during the first T r time steps and kept at the final value for an additional T p time steps. The initial value x i is set to a random value chosen from a zero-mean Gaussian distribution with a standard deviation of 0 . 1 and e i = 1. Furthermore, 3200 trajectories are computed per instance to evaluate TTS. The actual parameters used for simulation are listed below:
| N step | 1000 |
|----------|------------|
| ∆ T | 0.4 |
| T r | 900 |
| T p | 100 |
| p | -1.0 → 1.0 |
| α | 1.0 |
| β | 0.2 |
## CIM-SFC
Restriction of x i and e i variables is not needed as this system is more numerically stable. The parameters p , c and β are modulated linearly from their starting to ending values during simulation. The initial value x i is set to a random value chosen from a zero-mean Gaussian distribution with standard deviation of 0 . 1 and e i = 0. 3200 trajectories are computed per instance to evaluate TTS. Actual parameters used for simulation are listed below:
| N step | 500 |
|----------|------------|
| ∆ T | 0.4 |
| p | -1.0 → 1.0 |
| c | 1.0 → 3.0 |
| β | 0.3 → 0.1 |
| k | 0.2 |
In addition to the above-mentioned parameters, it is important for the normalizing factor ξ for the mutual coupling term to be chosen as [31],
$$\xi = \sqrt { \frac { 2 N } { \sum J _ { i j } ^ { 2 } } } .$$
This choice is crucial for the successful performance of CIM-SFC but not for CIM-CAC and CIM-CFC.
Moreover, it is important to note that we used the same number of time steps for all the problem sizes in Figures 10 and 11. It is likely that the optimal number of time steps is smaller for smaller problem sizes, thus the scaling of TTS when the number of time steps is optimized separately for each problem size might be slightly worse than the reported scaling. However, we do not believe that this difference would be very significant. For the scaling of TTS for CIM-CAC when different parameters are chosen, see Appendix C.
## dSBM
For Figure 12, dSBM is implemented as described in [29]. The parameters used are
$$\frac { N \text { step} } { \Delta T } & = 2 0 0 0 \\ c & = 0 . 5$$
## Parameters used in Figure 14
The parameters for N=800 are the same as those for Figure 12. The parameters for N=1200 are listed below. The number of trajectories used for N=1200 was 3200 for most instances; however, to accurately evaluate the success probability, 10000-50000 trajectories were computed for the 10 hardest instances for both algorithms. Moreover, owing to the hardness in the case of N=1200, we are not very certain that the true ground state was found.
## CIM-CAC
| N step | 8000 |
|----------|------------|
| ∆ T | 0.125 |
| T r | 7200 |
| T p | 800 |
| p | -1.0 → 1.0 |
| α | 1.0 → 2.5 |
| β | 0.8 |
## dSBM
## Numerical Integration
An Euler step is used for integration in all the cases (except for dSBM). As described above we constrain the range of x i variables to ensure numerical stability. This is not necessary for performance but allows us to increase the integration time step by a factor of 2 or 3 without compromising the success probability. In Figure 15 we show the success probability of CIM-CAC with respect to the time step for both constrained and unconstrained systems.
The results in Section 5 for CIM-CFC do not use this numerical constraint the CIM in Section 5 is meant to be a physical machine, and a time step of 0.2 is used.
## Appendix B: Simulation Results for G-set
The results in Figure 13 for dSBM are taken directly from the GPU implementation of dSBM in [29]. The unit for TTS in Table ?? is time steps to solution, or equivalently, MVM to solution. In our simulation, 3200, 10000, or 32000 trajectories were generated to evaluate the TTS depending on the instance difficulty. The numbers in bold denote the best TTS among the three algorithms.
Figure 15: Success probability of CIM-CAC with respect to time step for both constrained and unconstrained systems. For the blue curve, the x i amplitudes are restricted to the range [ -1 . 5 √ α, 1 . 5 √ α ]
<details>
<summary>Image 15 Details</summary>
### Visual Description
## Chart: Effect of Time Step (CIM-CAC)
### Overview
The image is a line chart comparing the average success probability of "constrained" and "unconstrained" scenarios with varying time steps. The x-axis represents the time step, and the y-axis represents the average success probability.
### Components/Axes
* **Title:** Effect of Time Step (CIM-CAC)
* **X-axis:** time step
* Scale: 0.05, 0.1, 0.15, 0.2
* **Y-axis:** average success probability
* Scale: 0.00, 0.05, 0.10, 0.15, 0.20
* **Legend:** Located in the top-right corner.
* Blue: "contrained"
* Light Green: "unconstrained"
### Detailed Analysis
* **Constrained (Blue):** The line starts at approximately 0.07 at time step 0.05, gradually decreases to approximately 0.06 at time step 0.1, remains around 0.058 at time step 0.125, then decreases to approximately 0.042 at time step 0.15, and finally reaches 0 at time step 0.2.
* **Unconstrained (Light Green):** The line starts at approximately 0.06 at time step 0.05, then drops sharply to 0 at time step 0.1 and remains at 0 for the rest of the time steps.
### Key Observations
* The "constrained" scenario consistently shows a higher average success probability than the "unconstrained" scenario, except at time step 0.2 where both are 0.
* The "unconstrained" scenario experiences a sharp drop in success probability between time steps 0.05 and 0.1.
* Both scenarios show a decrease in average success probability as the time step increases, but the "constrained" scenario decreases more gradually.
### Interpretation
The chart suggests that constraining the system leads to a more stable and higher average success probability, especially at lower time steps. The "unconstrained" scenario is highly sensitive to the time step, with a rapid decline in success probability. As the time step increases, the success probability decreases for both scenarios, indicating that larger time steps may negatively impact the system's performance. The data demonstrates the importance of constraints in maintaining a higher success rate, particularly when the time step is small.
</details>
## CIM-CAC parameters for G-set
The variables are restricted as described in Appendix A, and the initial conditions are set in the same way. The following parameters are the same for all the G-set instances.
$$\frac { \alpha } { \beta } | \begin{array} { c | c } 1 . 0 & 3 . 0 \\ 0 . 3 \end{array}$$
The parameters p , ∆ T , and the number of time steps used in each phase are chosen by instance type as follows:
| N step | 4000 |
|----------|--------|
| ∆ T | 1.25 |
| c | 0.5 |
| Graph Type | Edge Weight | N | Instance # | p | N step | ∆ T | T r | T p |
|--------------|---------------|------|--------------|-------------|----------|-------|-------|-------|
| Random | { +1 } | 800 | 1-5 | -0.5 → 1.0 | 6666 | 0.075 | 6000 | 666 |
| Random | { +1, -1 } | 800 | 6-10 | -0.5 → 1.0 | 6666 | 0.075 | 6000 | 666 |
| Toroidal | { +1, -1 } | 800 | 11-13 | -4.0 | 5000 | 0.1 | 4500 | 500 |
| Planar | { +1 } | 800 | 14-17 | -1.0 | 20000 | 0.05 | 18000 | 2000 |
| Planar | { +1, -1 } | 800 | 18-21 | -1.0 | 20000 | 0.05 | 18000 | 2000 |
| Random | { +1 } | 1000 | 43-46 | -0.5 → 1.0 | 10000 | 0.1 | 9000 | 1000 |
| Planar | { +1 } | 1000 | 51-54 | -1.0 | 20000 | 0.05 | 18000 | 2000 |
| Random | { +1 } | 2000 | 22-26 | -0.5 → 1.0 | 20000 | 0.1 | 19000 | 1000 |
| Random | { +1, -1 } | 2000 | 27-31 | -0.5 → 1.0 | 20000 | 0.1 | 19000 | 1000 |
| Toroidal | { +1, -1 } | 2000 | 32-34 | -4.0 → -3.0 | 20000 | 0.1 | 19000 | 1000 |
| Planar | { +1 } | 2000 | 35-38 | -1.0 → -0.5 | 80000 | 0.05 | 78000 | 2000 |
| Planar | { +1 } | 2000 | 39-42 | -1.0 → -0.5 | 80000 | 0.05 | 78000 | 2000 |
## CIM-CFC parameters for G-set
The variables are restricted as described in Appendix A, and the initial conditions are set in the same way. The following parameters are the same for all the G-set instances.
$$\frac { \alpha } { \beta } \left | \begin{array} { c c } { 1 . 0 } \\ { 0 . 1 5 } \end{array} \right |$$
The parameters p , ∆ T , and the number of time steps used in each phase are chosen by instance type as follows:
| Graph Type | Edge Weight | N | Instance # | p | N step | ∆ T | T r | T p |
|--------------|---------------|------|--------------|-------------|----------|-------|-------|-------|
| Random | { +1 } | 800 | 1-5 | -1.0 → 1.0 | 4000 | 0.125 | 3600 | 400 |
| Random | { +1, -1 } | 800 | 6-10 | -1.0 → 1.0 | 2000 | 0.25 | 1800 | 200 |
| Toroidal | { +1, -1 } | 800 | 11-13 | -3.0 → -1.0 | 2000 | 0.25 | 1800 | 200 |
| Planar | { +1 } | 800 | 14-17 | -2.0 → 0.0 | 8000 | 0.125 | 7200 | 800 |
| Planar | { +1, -1 } | 800 | 18-21 | -2.0 → 0.0 | 4000 | 0.25 | 3600 | 400 |
| Random | { +1 } | 1000 | 43-46 | -1.0 → 1.0 | 5000 | 0.2 | 4500 | 500 |
| Planar | { +1 } | 1000 | 51-54 | -2.0 → 0.0 | 16000 | 0.125 | 15200 | 800 |
| Random | { +1 } | 2000 | 22-26 | -1.0 → 1.0 | 10000 | 0.2 | 9500 | 500 |
| Random | { +1, -1 } | 2000 | 27-31 | -1.0 → 1.0 | 10000 | 0.2 | 9500 | 500 |
| Toroidal | { +1, -1 } | 2000 | 32-34 | -3.0 → -1.0 | 40000 | 0.1 | 39000 | 1000 |
| Planar | { +1 } | 2000 | 35-38 | -2.0 → 0.0 | 80000 | 0.05 | 78000 | 2000 |
| Planar | { +1 } | 2000 | 39-42 | -2.0 → 0.0 | 40000 | 0.1 | 39000 | 1000 |
## Results on G-set
Success probability and TTS of CIM-CAC, CIM-CFC and dSBM on G-set graphs. The success probability for CIM-CAC and CIM-CFC is for a single trajectory. The results for dSBM are from the GPU implementation in [29].
| Instance | CIM-CAC TTS | CIM-CAC P s | CIM-CFC TTS | CIM-CFC P s | dSBM TTS | dSBM P s * |
|------------|--------------------|---------------|---------------|---------------|---------------------|--------------|
| G1 | 90805 | 0.286875 | 60078 | 0.264062 | 339332 | 0.987 |
| G2 | 920217 | 0.0328125 | 1330454 | 0.01375 | 2578124 | 0.82 |
| G3 | 169551 | 0.165625 | 249212 | 0.07125 | 400343 | 0.996 |
| G4 | 209086 | 0.136563 | 239389 | 0.0740625 | 361673 | 0.983 |
| G5 | 226881 | 0.126562 | 226449 | 0.078125 | 618221 | 0.972 |
| G6 | 188908 | 0.15 | 104083 | 0.0846875 | 190728 | 0.979 |
| G7 | 562413 | 0.053125 | 146490 | 0.0609375 | 201889 | 0.974 |
| G8 | 431029 | 0.06875 | 399118 | 0.0228125 | 358947 | 0.954 |
| G9 | 411604 | 0.071875 | 223836 | 0.0403125 | 1095704 | 0.867 |
| G10 | 1388073 | 0.021875 | 622470 | 0.0146875 | 1410031 | 0.407 |
| G11 | 337563 | 0.0659375 | 222079 | 0.040625 | 282524 | 0.98 |
| G12 | 224452 | 0.0975 | 78562 | 0.110625 | 407997 | 0.973 |
| G13 | 391011 | 0.0571875 | 373236 | 0.024375 | 800687 | 0.996 |
| G14 | 17291018 | 0.0053125 | 13080721 | 0.0028125 | 1469967245 | 0.005 |
| G15 | 521572 | 0.161875 | 462528 | 0.0765625 | 6782113 | 0.804 |
| G16 | 494303 | 0.17 | 737146 | 0.04875 | 9156329 | 0.992 |
| G17 | 3089154 | 0.029375 | 3256332 | 0.01125 | 33222397 | 0.283 |
| G18 | 556635 | 0.1525 | 507805 | 0.035625 | 14375986 | 0.074 |
| G19 | 1218307 | 0.0728125 | 179562 | 0.0975 | 417204 | 0.995 |
| G20 | 106718 | 0.578125 | 42428 | 0.352187 | 188349 | 0.98 |
| G21 | 3991180 | 0.0228125 | 574365 | 0.0315625 | 10080921 | 0.136 |
| G43 | 174031 | 0.2325 | 145625 | 0.14625 | 228908 | 0.992 |
| G44 | 244188 | 0.171875 | 257230 | 0.085625 | 263171 | 0.985 |
| G45 | 880856 | 0.0509375 | 970877 | 0.0234375 | 1754473 | 0.985 |
| G46 | 1528073 | 0.0296875 | 717957 | 0.0315625 | 610421 | 0.992 |
| G51 | 3732360 | 0.024375 | 2189016 | 0.0331 | 424989992 | 0.067 |
| G52 | 3122871 | 0.0290625 | 3820774 | 0.0191 | 92285270 | 0.213 |
| G53 | 17291018 | 0.0053125 | 11298922 | 0.0065 | 1676440469 | 0.043 |
| G54 | 294684837 | 0.0003125 | 1178886725 | 6.25e-05 | 49107077298 | 0.0006 |
| G22 | 2516544 | 0.0359375 | 2718081 | 0.0168 | 8401394 | 0.928 |
| G23 | N/A | 0.0 | N/A | 0.0 | N/A | 0 |
| G24 | 4785454 | 0.0190625 | 5733406 | 0.008 | 10585339 | 0.648 |
| G25 | 17291018 | 0.0053125 | 15327529 | 0.003 | 43414254 | 0.399 |
| G26 | 10117012 | 0.0090625 | 10941648 | 0.0042 | 10730290 | 0.643 |
| G27 | 835530 | 0.104375 | 649972 | 0.0684 | 832465 | 0.971 |
| G28 | 2389444 | 0.0378125 | 2413498 | 0.0189 | 1455914 | 0.952 |
| G29 | 4164219 | 0.021875 | 7404644 | 0.0062 | 5516820 | 0.737 |
| G30 | 42058344 | 0.0021875 | 32871041 | 0.0014 | 11002259 | 0.738 |
| G31 | 49075749 | 0.001875 | 92080375 | 0.0005 | 19923732 | 0.199 |
| G32 | 70802710 | 0.0013 | 1841975969 | 0.0001 | 150969342 | 0.093 |
| G33 | 306965291 | 0.0003 | N/A | 0.0 | 2204950868 | 0.005 |
| G34 | 20886506 | 0.0044 | 32801883 | 0.0056 | 84156149 | 0.231 |
| G35 | N/A | 0.0 | N/A | 0.0 | N/A | 0 |
| G36 | 368321503 | 0.0005 | 3683951938 | 0.0001 | 736790387782 | 0.0001 |
| G38 | 108264816 | 0.0017 | | 0.0025 | 1046295719 | 0.068 |
| | | | 147181162 | | | 0.107 |
| G39 | 4065368 1841975969 | 0.0443 0.0001 | 3497864 N/A | 0.0513 0.0 | 651087484 264354894 | 0.154 |
| G40 G41 | 17123320 | 0.0107 | 7916531 | 0.023 | 222414431 | 0.282 |
| G42 | 13969282 | 0.0131 | 18328423 | 0.01 | N/A | 0 |
*Note that P s for dSBM taken from [29] is the success probability for a batch of 160 trajectories. To calculate the P s for a single dSBM trajectory use 1 -(1 -P s ) 1 160 .
## Appendix C: Reasoning for parameter selection
The parameters are selected numerically for the most part; however, the choice of p , α , and β can be understood as follows. It is observed that the average residual energy visited by CIM-CAC during the search process can be roughly estimated by the formula. [ 11 ]
$$\Delta E _ { a v g } \approx K \frac { 1 - p } { \alpha \beta }$$
where K is a constant depending only on the problem type and size. This formula essentially predicts the effective sampling temperature of the system (although the distribution may not be an exact Boltzmann distribution). Based on this philosophy, we gradually reduce the 'system temperature' to produce an annealing effect. This is the motivation for increasing p and α . The different choices for the range of p on different G-set instances reflects the vastly different values for the constant K depending on the structure of the max-cut problem.
In a more general setting, the value of K can be predicted on the basis of the problem type; thus, the range for p and α can be chosen accordingly.
Although it has not been verified, a similar formula most likely holds for CIM-CFC; thus, the parameters for CIM-CFC are chosen in the same way.
## Optimal parameters with respect to problem size (CIM-CAC)
Figure 16: Performance of CIM-CAC with respect to problem size for different parameters. The fixed parameters are indicated in red while the parameter modulation is indicated in blue. The red and blue dotted lines are fits for the lower envelopes of the red and blue curves, respectively
<details>
<summary>Image 16 Details</summary>
### Visual Description
## Chart: Median TTS vs Problem Size
### Overview
The image is a line chart that plots the median Time To Solution (TTS) against the problem size. The y-axis (time steps to solution) is on a logarithmic scale, ranging from 10^4 to 10^8. The x-axis (problem size) displays the square root of the problem size, with values ranging from approximately 10 to 28 (corresponding to problem sizes of 100 to 800). Several lines represent different configurations or parameters, as indicated by the legend.
### Components/Axes
* **Title:** Median TTS vs Problem Size
* **X-axis:** problem size (values: √100, √200, √300, √400, √500, √800)
* **Y-axis:** time steps to solution (logarithmic scale, values: 10^4, 10^5, 10^6, 10^7, 10^8)
* **Legend (located on the right side of the chart):**
* Blue: T max = 400
* Dark Blue: T max = 200
* Medium Blue: T max = 100
* Light Blue: T max = 50
* Dark Red: fixed p=2.0 α=3.0
* Red: fixed p=1.5 α=2.5
* Light Red: fixed p=2.0 α=2.0
* Pink: fixed p=0.0 α=1.5
### Detailed Analysis
* **T max = 400 (Blue):** The line starts at approximately (√100, 2 * 10^5) and increases to approximately (√800, 7 * 10^6). The trend is upward.
* **T max = 200 (Dark Blue):** The line starts at approximately (√100, 1.5 * 10^5) and increases to approximately (√800, 6 * 10^6). The trend is upward.
* **T max = 100 (Medium Blue):** The line starts at approximately (√100, 8 * 10^4) and increases to approximately (√800, 4 * 10^6). The trend is upward.
* **T max = 50 (Light Blue):** The line starts at approximately (√100, 4 * 10^4) and increases to approximately (√800, 2 * 10^6). The trend is upward.
* **fixed p=2.0 α=3.0 (Dark Red):** The line starts at approximately (√100, 7 * 10^4) and increases to approximately (√500, 4 * 10^6). The trend is upward.
* **fixed p=1.5 α=2.5 (Red):** The line starts at approximately (√100, 4 * 10^4) and increases to approximately (√500, 2 * 10^6). The trend is upward.
* **fixed p=2.0 α=2.0 (Light Red):** The line starts at approximately (√100, 2 * 10^4) and increases to approximately (√500, 1 * 10^6). The trend is upward.
* **fixed p=0.0 α=1.5 (Pink):** The line starts at approximately (√100, 1.5 * 10^4) and increases to approximately (√300, 1.5 * 10^6). The trend is upward. This line is dashed.
### Key Observations
* All lines show an upward trend, indicating that the median time to solution increases as the problem size increases.
* The "T max" lines (blue shades) generally show a less steep increase compared to the "fixed" lines (red shades).
* The "fixed p=0.0 α=1.5" line (pink) has the steepest increase initially.
### Interpretation
The chart illustrates the relationship between problem size and the median time to solution (TTS) for different configurations. The upward trend in all lines suggests that larger problem sizes require more time steps to find a solution. The different slopes of the lines indicate that certain configurations (e.g., "fixed" parameters) are more sensitive to problem size increases than others (e.g., "T max" parameters). The "fixed p=0.0 α=1.5" configuration appears to scale very poorly with problem size, as indicated by its steep initial increase. The data suggests that tuning parameters can significantly impact the scalability of the solution.
</details>
Figure 16 shows the difference in scaling when the parameters are fixed (red shades) compared to when the parameters are modulated linearly (blue shades). In addition, the optimal annealing time (in other words, the optimal speed of modulation) will change with respect to the problem size, as we need long annealing times to get good results for large problem sizes. This pattern was also used when choosing parameters on the G-set.
Although Figure 16 is based on the Gaussian quantum model for CIM-CAC in MFB-CIM (see [31]), the difference in performance between this model and the noiseless model discussed in this paper is insignificant, as g 2 = 10 -4 was used. It should also be noted that a time step of 0.01 was used in Figure 16.
Therefore, the TTS in Figure 16 is an order of magnitude longer than the results presented in this study, where a time step of 0.125 was used.
## Appendix D: Results and Discussion for CIM-SFC on G-set
Based on our understanding of CIM-SFC, it is very important that the term tanh( cz i ) transitions from the 'soft spin' mode where cz i ≈ 0 and tanh( cz i ) ≈ cz i to the 'discrete spin' mode where | cz i | >> 0 and tanh( cz i ) ≈ sign( cz i ). Therefore, we use the normalizing factor ξ (as defined above) as this ensures that z i will on average be around √ 2 for a randomly chosen spin configuration, thus we can use the same value for c in all the cases and get similar results. However, this only works for instances such as SK instances where each node has equal connectivity; thus, we can expect z i to have roughly the same range of values for all i .
On some G-set instances, especially the planar graph instances, some nodes have a much larger degree; thus, cz i will be too large in some cases and too small in others regardless of the normalizing factor ξ used. This may be one of the reasons why CIM-SFC struggles on many G-set instances, especially planar graphs. This could also be the reason why dSBM struggles on planar graphs, as dSBM relies on the same normalizing factor to get good results. CIM-CAC and CIM-CFC do not need this normalizing factor, as they automatically compensate for different values of ∑ j J ij σ j , and this might be why they perform well on planar graphs.
Meanwhile, for toroidal graphs, the opposite is true, as ∑ j J ij σ j can only take on five different values for these graphs. This could mean that the transition from 'soft spin' to 'discrete spin' is rapid in the case of CIM-SFC; thus, we need to carefully tune the parameters to get good results on these graphs.
Although this observation regarding the analog/discrete transition may partially explain the poor results on the G-set, it is not a complete explanation. For example, CIM-SFC struggles on some random graphs (such as G9) that do not have the above-mentioned property, as each node has similar connectivity.
The results for CIM-SFC on the G-set as well as the parameters used are listed below (not all instances were tested).
## Results for CIM-SFC on the G-set
| Instance | TTS | P s |
|------------|----------|-----------|
| G1 | 28470 | 0.194 |
| G2 | 1531984 | 0.004 |
| G3 | 130388 | 0.046 |
| G4 | 435510 | 0.014 |
| G5 | 380685 | 0.016 |
| G6 | 112834 | 0.0202 |
| G7 | 163316 | 0.014 |
| G8 | 125360 | 0.0182 |
| G9 | 4604018 | 0.0005 |
| G10 | 395845 | 0.0058 |
| G11 | 195461 | 0.0572 |
| G12 | 128184 | 0.0859 |
| G13 | 311368 | 0.0363 |
| G43 | 1702012 | 0.0134375 |
| G44 | 1742812 | 0.013125 |
| G45 | 73671209 | 0.0003125 |
| G46 | 1927483 | 0.011875 |
For instances G14-G21 (800 node planar graphs) and G51-G54 (1000 node planar graphs), CIM-SFC shows a success probability of either zero or a very small nonzero value. Furthermore, 2000 node instances have not been tested.
## Parameters for CIM-SFC on G-set
## Common parameters
Parameters selected by problem type
| Graph Type | Edge Weight | N | Instance # | c | β | k | N step | ∆ T |
|--------------|---------------|------|--------------|-----------|------------|------|----------|-------|
| Random | { +1 } | 800 | 1-5 | 1.0 → 3.0 | 0.3 → 0.0 | 0.2 | 2666 | 0.15 |
| Random | { +1, -1 } | 800 | 6-10 | 1.0 → 3.0 | 0.3 → 0.0 | 0.2 | 500 | 0.4 |
| Toroidal | { +1, -1 } | 800 | 11-13 | 1.4 | 0.05 → 0.0 | 0.32 | 2500 | 0.4 |
| Random | { +1 } | 1000 | 43-46 | 1.4 → 4.2 | 0.2 → 0.0 | 0.2 | 5000 | 0.2 |
The parameters for CIM-SFC are chosen experimentally, and the understanding of how the parameters affect the performance and dynamics is limited. Once this system is studied more thoroughly, we will propose a more systematic method of choosing parameters so that good performance can be ensured on many different problem types.
## Appendix E: Similarities and Differences Between CIM and SBM Algorithms
Using continuous analog dynamics to solve discrete optimization problems is a somewhat new concept, and it is interesting to compare these different approaches. [ 13 , 11 , 29 ] In this appendix, we will briefly discuss some similarities and differences among the three CIM-inspired algorithms and the SBM algorithms.
All four systems discussed in Section 6, namely CIM-CAC, CIM-CAC, CIM-SFC, and dSBM, were originally inspired by the same fundamental principle:[10, 27]:
The function
$$H ( x ) = \sum _ { i } \left ( \frac { x _ { i } ^ { 2 } } { 4 } - \frac { 1 - p } { 2 } \right ) x _ { i } ^ { 2 } + c \sum _ { i } \sum _ { j } J _ { i j } x _ { i } x _ { j }$$
can be used as a continuous approximation of the Ising cost function.
In the original CIM algorithm, gradient descent is used to find the local minima of H , while H is deformed by increasing p . This system has two major drawbacks:[10]
1. local minima are stable;
2. incorrect mapping of the Ising problem to the cost function owing to amplitude heterogeneity.
All four algorithms discussed in Section 6 can be regarded as modifications of the original CIM algorithm, which aim to overcome these two flaws. [ 11 , 27 , 28 , 29 ] In all these algorithms, the first flaw is addressed by adding new degrees of freedom to the system; hence, there are now 2 N (instead of only N ) analog variables for N spins. In SBM, this is done by including both a position vector, x i , and a velocity/momentum vector y i , while in the modified CIM algorithms we add the auxiliary variable e i .
To address the second flaw, the creators of dSBM added discretization and 'inelastic walls', whereas in CIM- CFC and CIM-SFC, this discretization is not necessary. Using different mechanisms, all three algorithms ensure that the system only has fixed points at the local minima of the Ising Hamiltonian (during the end of the trajectory), something which is not true for the original CIM algorithm. Because these systems are fundamentally very similar, it should not be surprising that they achieve similar performance.
$$p \, - 1 . 0 \to 1 . 0$$
We also note that for dSBM to achieve good performance, it is necessary to use discretization and inelastic walls, which make the system discontinuous. This is particularly useful for implementation on a digital platform, which prefers discrete processes; however, when implementing these algorithms on an analog physical platform, this is not preferred. Meanwhile, in the case of CIM-CAC, CIM-CFC, and CIM-SFC, the system evolves continuously; thus, they are much more suitable for analog implementation, such as the optical CIM architecture proposed in this paper.
An interesting difference between the CIM and the original bifurcation machine [ 27 ] , which was referred to as aSBM in [29], is that aSBM is a completely unitary dissipation-less system. Because of this, aSBM relies on adiabatic evolution for computation (similar to quantum annealing), unlike the dissipative CIM and other Ising heuristics (such as simulated annealing or breakout local search [ 14 ] ) , which rely on some sort of dissipative relaxation. However, in [29], the new SBM algorithms deviate from this concept of adiabatic evolution by adding inelastic walls, thus making the new bifurcation machine a dissipative system in which information is lost over time. It would be interesting to try to understand whether dissipation is in fact necessary for a system to achieve the high performance of the algorithms discussed in this paper. For example, one could modify aSBM in a different way that addresses the problem of amplitude heterogeneity but retains the adiabatic nature. Whether this is possible is beyond the scope of this paper.
## Appendix F: Optical Implementation of CIM-SFC
Figure 17: Optical implementation of CIM-SFC.
<details>
<summary>Image 17 Details</summary>
### Visual Description
## Diagram: Quantum Reservoir Computing Architecture
### Overview
The image depicts a quantum reservoir computing architecture. It shows the flow of quantum information through a series of processing stages, including phase-sensitive amplifiers (PSAs), delay lines (DLs), fan-out and fan-in operations, and homodyne detection. The diagram illustrates how input signals are processed and transformed within the quantum reservoir to produce output results.
### Components/Axes
* **Input:** BS<sub>e</sub> (Beam Splitter input)
* **Output:** BS<sub>i</sub> (Beam Splitter output)
* **Processing Units:**
* PSA<sub>0</sub>: Phase-Sensitive Amplifier 0
* Fan-Out: Signal distribution unit
* PSA<sub>1</sub> + DL<sub>1</sub>: Phase-Sensitive Amplifier 1 + Delay Line 1
* PSA<sub>2</sub> + DL<sub>2</sub>: Phase-Sensitive Amplifier 2 + Delay Line 2
* PSA<sub>N</sub> + DL<sub>N</sub>: Phase-Sensitive Amplifier N + Delay Line N
* Fan-In: Signal aggregation unit
* PSA<sub>e</sub>: Phase-Sensitive Amplifier e
* DL<sub>e</sub>: Delay Line e
* **Signals:**
* x̃<sub>j</sub>: Input signal
* x̃<sub>1</sub>(τ): Signal after Fan-Out, first branch
* x̃<sub>2</sub>(3τ): Signal after Fan-Out, second branch
* x̃<sub>N</sub>(2N-1τ): Signal after Fan-Out, Nth branch
* ξ ∑ J<sub>ij</sub>x̃<sub>j</sub>(ω): Summation of weighted signals
* tanh(c z̃<sub>i</sub><sup>(m)</sup>) / z̃<sub>i</sub><sup>(m)</sup>: Non-linear transformation
* **Detection:**
* Homodyne detection (two instances)
* **Results:**
* Results x̃<sub>j</sub><sup>(m)</sup>: Output from the first homodyne detection
* Results z̃<sub>i</sub><sup>(m)</sup>: Output from the second homodyne detection
* **Weights:**
* J<sub>ij</sub>: Weights applied to the signals
### Detailed Analysis
1. **Input Stage:** The input signal x̃<sub>j</sub> enters the system through a beam splitter (BS<sub>e</sub>) and is processed by the first phase-sensitive amplifier (PSA<sub>0</sub>).
2. **Fan-Out:** The signal is then distributed by the Fan-Out unit into N branches. The signals in these branches are denoted as x̃<sub>1</sub>(τ), x̃<sub>2</sub>(3τ), ..., x̃<sub>N</sub>(2N-1τ).
3. **Processing Branches:** Each branch consists of a phase-sensitive amplifier (PSA<sub>i</sub>) and a delay line (DL<sub>i</sub>), where i ranges from 1 to N.
4. **Fan-In:** The signals from all N branches are aggregated by the Fan-In unit. The aggregated signal is represented as ξ ∑ J<sub>ij</sub>x̃<sub>j</sub>(ω), where J<sub>ij</sub> represents the weights applied to the signals.
5. **Output Stage:** The aggregated signal is then processed by another phase-sensitive amplifier (PSA<sub>e</sub>) and a delay line (DL<sub>e</sub>). A non-linear transformation, tanh(c z̃<sub>i</sub><sup>(m)</sup>) / z̃<sub>i</sub><sup>(m)</sup>, is applied before PSA<sub>e</sub>. The final output signal exits through a beam splitter (BS<sub>i</sub>).
6. **Homodyne Detection:** Homodyne detection is performed after PSA<sub>0</sub> and after the Fan-In unit. The results of these detections are x̃<sub>j</sub><sup>(m)</sup> and z̃<sub>i</sub><sup>(m)</sup>, respectively.
### Key Observations
* The architecture uses multiple phase-sensitive amplifiers and delay lines to process the input signal.
* The Fan-Out and Fan-In units distribute and aggregate the signals, respectively.
* Homodyne detection is used to measure the quantum states at different stages of the processing.
* The non-linear transformation tanh(c z̃<sub>i</sub><sup>(m)</sup>) / z̃<sub>i</sub><sup>(m)</sup> introduces non-linearity into the system.
### Interpretation
The diagram illustrates a quantum reservoir computing architecture, which leverages the complex dynamics of a quantum system to perform computational tasks. The input signal is processed through a series of non-linear transformations and feedback loops within the reservoir, allowing the system to learn and adapt to different input patterns. The homodyne detection provides a means to measure the quantum states and extract the computational results. The architecture is designed to exploit the unique properties of quantum mechanics, such as superposition and entanglement, to achieve computational advantages over classical systems. The weights J<sub>ij</sub> are crucial for training the reservoir to perform specific tasks. The delay lines introduce temporal dynamics, which are essential for processing time-dependent signals.
</details>
Figure 17 shows an optical implementation of CIM-SFC which is similar to that of CIM-CAC and CIMCFC shown in Figure 6. The feedback signal ˜ z i = ∑ j J ij ˜ x j is deamplified (rather than amplified) by PSA e with attenuation coefficient tanh( ˜ z i ( m ) ) ˜ z i ( m ) , where ˜ z i ( m ) is an optical homodyne measurement result of ˜ z i . This feedback signal is then injected back into signal pulse x i in the main cavity through BS i .
Part of the fan-in circuit output ˜ z i is delayed by a delay line DL e with delay time Nτ and combined with the error pulse e i (inside the main cavity). This implements the term e i -˜ z i in Eq. (8). The term -β ( e i -˜ z i ) on the right-hand side of Eq. (8) is implemented by a phase-sensitive amplifier PSA e of the main cavity. This is also a deamplification process. Finally, the error correction signal amplitude e i -˜ z i is coupled to the signal pulse x i inside the main cavity with a standard optical delay line.
## Conflict of Interest
The authors have no confict of interest, financial or otherwise.
## Acknowledgements
The authors wish to thank R. Hamerly, M. G. Suh, M. Jankowski, E. Ng and Y. Inui for their valuable discussions.
## References
- [1] F. Barahona, J. Phys. A 1982 , 15 , 3241.
- [2] 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, G. Rose, Nature 2011 , 473 , 194.
- [3] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, M. Troyer, Nat. Phys. 2014 , 10 , 218.
- [4] A. Marandi, Z. Wang, K. Takata, R. L. Byer, Y. Yamamoto, Nat. Photonics 2014 , 8 , 937.
- [5] T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, H. Takesue, Nat Photon. 2016 , 10 , 415.
- [6] P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, R. L. Byer, M. M. Fejer, H. Mabuchi, Y. Yamamoto, Science 2016 , 354 , 614.
- [7] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Takenouchi, K. Aihara, K. Kawarabayashi, K. Inoue, S. Utsunomiya, H. Takesue, Science 2016 , 354 , 603.
- [8] 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, Y. Yamamoto, Sci. Adv. 2019 , 5 , eaau0823.
- [9] K. Sankar, A. Scherer, S. Kako, S. Reifenstein, N. Ghadermarzy, W. B. Krayenhoff, Y. Inui, E. Ng, T. Onodera, P. Ronagh, Y. Yamamoto, arXiv:2105.03528, v1.
- [10] Z. Wang, A. Marandi, K. Wen, R. L. Byer, Y. Yamamoto, Phys. Rev. A 2013 , 88 , 063853.
- [11] T. Leleu, Y. Yamamoto, P. L. McMahon, K. Aihara, Phys. Rev. Lett. 2019 , 122 , 040607.
- [12] T. Leleu, F. Khoyratee, T. Levi, R. Hamerly, T. Kohno, K. Aihara, arXiv:2009.04084, v3
- [13] M. Ercsey-Ravasz, Z. Toroczkai, Nature Phys , 2011 , 966 970.
- [14] U. Benlic, J. K. Hao, Eng. Appl. Artif. Intell. 2013 , 26 , 1162.
- [15] C. Wang, C. Langrock, A. Marandi, M. Jankowski, M. Zhang, B. Desiatov, M. M. Fejer, M. Lonˇ car, Optica 2018 , 5 , 1438.
- [16] B. Stern, X. Ji, Y. Okawachi, A. L. Gaeta, M. Lipson, Nature 2018 , 562 , 401.
- [17] C. Wang, M. Zhang, X. Chen, M. Bertrand, A. Shams-Ansari, S. Chandrasekhar, P. Winzer, M. Lonˇ car, Nature 2018 , 562 , 101.
- [18] P. D. Drummond, C W Gardiner, J. Phys. A 1980 , 13 , 2353.
- [19] P. D. Drummond, C. W. Gardiner, D. F. Walls, Phys. Rev. A 1981 , 24 , 914.
- [20] D. F. Walls, G. J. Milburn, Quantum Optics (Springer Science and Business Media) 2007 .
- [21] K. Takata, A. Marandi, Y. Yamamoto, Phys. Rev. A 2015 , 92 , 043821.
- [22] D. Maruo, S. Utsunomiya, Y. Yamamoto, Phys. Scr. 2016 , 91 , 083010.
- [23] Y. Inui, Y. Yamamoto, Phys. Rev. A 2020 , 102 , 062419.
- [24] Y. Inui, Y. Yamamoto, Entropy 2021 , 23 , 624.
- [25] E. Ng, T. Onodera, S. Kako, P. L. McMahon, H. Mabuchi, Y. Yamamoto, arXiv:2103.05629, v1.
- [26] A. D. King, W. Bernoudy, J. King, A. J. Berkley, T. Lanting, arXiv:1806.08422, v1
- [27] H. Goto, K. Tatsumura, A. R. Dixon, Sci. Adv. 2019 , 5 , eaav2372.
- [28] K. Tatsumura, M. Yamasaki, H. Goto, Nat. Electron. 2021 , 4 , 208.
- [29] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, K. Tatsumura, Sci. Adv. 2021 , 7 , eabe795.
- [30] Y. Inui, Y. Yamamoto, arXiv:2009.10328.
- [31] S. Kako, T. Leleu, Y. Inui, F. Khoyratee, S. Reifenstein, Y. Yamamoto, Adv. Quantum Technol. 2020 , 3 , 2000045.