Control of amplitude homogeneity in coherent Ising machines with artificial Zeeman terms

A coherent Ising machine (CIM) is an open-dissipative Ising solver using optical pulses generated from a degenerate optical parametric oscillator as analog magnetizations. When solving real-world optimization problems with CIM, this solver has two difficulties: mutual coupling induced amplitude inhomogeneity and absence of natural way to implement Zeeman terms. For the approximate Gaussian formulation of CIMs with amplitude control feedback, we add artificial Zeemam terms using the target amplitude information. Here we show, for 16-spin CIM with Zeeman terms, the amplitude control increases the performance, particularly when Zeeman terms are competing against mutual coupling coefficients. Coherent Ising machines represent an optical approach to unconventional computing that fail if inhomogeneity of the analog spins’ amplitude becomes too great. Here, an approach to amplitude control is shown to improve their performance even when Zeeman terms are included.

A s the Moore's law stagnates, special purpose accelarators to solve particular problems are expected to mitigate increasing time and energy costs for computation. Finding ground states of non-planar Ising models has been known to be NP hard in computational complexity theory, and has also drawn attention since various combinatorial optimization problems can be mapped into Ising Hamiltonian. There have been various attempts to use CMOS accelerators or physical systems for solving Ising problems [1][2][3][4][5][6][7][8][9] . A coherent Ising machine (CIM) [10][11][12][13][14][15][16][17][18] uses the optical pulses in degenerate optical parametric oscillators (DOPOs) as binary spins to find the ground state of Ising problem In CIM, optical pulses are mutually coupled by dissipative circuits rather than unitary gates to realizeJ rr 0 coupling. There have been various attempts to use it for solving combinatorial optimization problems, including traveling salesman problem 19 , lead optimization in drug discovery 20 , multiple-input multipleoutput optimization for wireless communication 21 , and compressed sensing for medical imaging 22 . In the attempts to use CIMs for real world problems, there have been two difficulties, the inhomogeneity of amplitudes 23 and the lack of natural Zeeman terms, which are denoted ash r in Eq. (1). Particularly, in a previous attempt 20 , implementation of Zeeman terms with mean absolute amplitudes didn't work if the additional parameter to represent the strength of Zeeman terms deviated from the optimal value. On the other hand, the correction of pulse-amplitudes' inhomogeneity using the feedback technique has been suggested 24 , and used for solving random weight problems [25][26][27] . In this paper, we extend the model of amplitude controlled CIM and propose the implementation scheme of Zeeman terms.

Results and discussion
Coherent Ising machine. A coherent Ising machine with measurement feedback coupling has a ring cavity structure shown in Fig. 1a. Ising spins are represented by the signal pulses generated by χ (2) degenerate optical parametric amplification in periodically poled LiNbO 3 (PPLN) waveguide device. Small parts of generated signal pulses are extracted by outlet coupler and their amplitudes are measured by optical homodyne detection. With the measured amplitudes (we denote asμ r for r-th pulse), the coherent amplitude of feedback pulse to the r-th signal pulse is calculated in the digital circuit (FPGA) to the form of matrix vector product ∑ r 0J rr 0μ r 0 . By using intensity and phase modulation, the amplitudes of the feedback pulses are set to the calculated values. The feedback pulses are injected from the injection coupler.
We present the quantum theoretical model of measurementfeedback CIM, that are obtained by assuming the round-trip-time (Δt) is much smaller than the linear dissipation time. When we normalize the linear dissipation time to 1, The master equation of r-th DOPO is written as ∂ρ ∂t Here,â r is the annihilation operator of the r-th signal mode. The first term on the right hand side (R.H.S.) represents the linear dissipation of the signal mode. The second and third terms show the effect of χ (2) nonlinear interaction between signal and pump modes shown in Fig. 1b. When the r-th pump mode is coherently excited, the r-th signal mode has the parametric gain represented by the HamiltonianĤ r ¼ i_ 2 pðâ y2 r Àâ 2 r Þ. Inversely, the signal mode photons are converted into the pump mode photon by two photon absorption (see Supplementary Note 1). The amplitude μ of a solitary DOPO follow the equation of motion dμ dt ¼ À ∂VðμÞ ∂μ , where V(μ) is the effective potential VðμÞ ¼ ð1 À pÞ μ 2 2 þ g 2 4 μ 4 23 . As shown in Fig. 1c, below the oscillation threshold (p < 1), the effective potential has a minimum at μ = 0. On the other hand, it has two minima above the threshold (p > 1). When we assume that the DOPO has the spin value σ ¼ μ jμj , it transits from fluctuating random value for p < 1 to the stable value representing the low energy phase of spin.
To realize coupling of optical pulses, as shown in Fig. 1a, the measurement feedback CIM has two couplers added to the ring cavity, which are called outlet coupler and injection coupler. We assume two couplers have the same reflectance R B (we assume R B ≪ 1) which is related to the coupling induced linear dissipation rate j via j = R B Δt. By the outlet coupler and homodyne detection, the master equation has the additional terms which represent measurement-induced linear loss and state-reduction 28 : ∂ρ ∂t S:R: The first term of R.H.S. represents the additional loss by adding outlet coupler, and the second term represents the impact of homodyne measurement. W r are real random numbers representing the quantum vacuum noise incident from the open port of the outlet coupler. By homodyne measurement, one of the real numbers is selected from normally distributed vacuum noise, and measured values satisfy W r ðtÞW r 0 ðt 0 Þ ¼ δ rr 0 δðt À t 0 Þ under ensemble averaging. Related to the injection coupler, the master Fig. 1 Model of Coherent Ising machine (CIM). a CIM with measurement and feedback (cited from ref. 16 , SHG, ADC, FPGA, DAC, IM/PM, PPLN, and LO represent second harmonic generation, analog-to-digital converter, field programmable gate array, digital-to-analog converter, intensity modulator/phase modulator, periodically poled lithium niobate, and local oscillator, respectively). b χ (2) interaction of optical modes via virtual excitation of valence band electronic state. c Effective potential V(μ) of mean amplitude μ for p = 0.5 and 2 (g 2 = 10 −3 ). equation has the following terms 28,29 : ∂ρ ∂t F:B: The first term of R.H.S. represents the additional loss by injection coupler, and the second term represents the coherent feedback injection. The total master equation of measurement feedback CIM is written as ∂ρ ∂t ¼ ∑ r .
Gaussian approximation of positive-P model. When the two photon absorption coefficient g 2 is small, the direct calculation of quantum master equation using photon number states requires huge computational cost (Above-threshold photon number has the order of g −2 for single DOPO.). Instead, we use the positive-P distribution function 30 , which is the generalization of Glauber-Sudarshan diagonal P function 31,32 . To model measurement-induced state reduction in positive-P simulation, we consider the following form instead of Eq. (3), Here,X r ¼â r þâ y r ffiffi 2 p and ΔX r ¼X r À hX r i. The outlet-coupling loss, measurement induced mean amplitude shift 33 , and measurement induced fluctuation reduction 33 are introduced in the form of Gaussian measurement theory 34,35 .
The master equation composed of Eqs. (2), (4) and (5) is exactly converted to the positive-Pc-number stochastic differential equations (CSDEs). As shown in Supplementary Note 2, we used Gaussian approximation by neglecting higher order fluctuation products, to the positive-P CSDEs (shown as Supplementary Eqs. (S10) and (S11)). By adding amplitude control 24,26 and Zeeman terms, we obtained following equations of Gaussian approximated positive-P (GAPP) model, Here, μ r is the mean amplitude, n r and m r represent variances of quantum fluctuations of r-th optical pulse. The auxiliary parameter e r is introduced to reduce amplitude inhomogeneity 24,26 by modulating the mutual coupling term as shown in Eq. (7). τ(t) is the target value of normalized squared amplitudes 24,26 . β is the strength to fix squared amplitudes to the target value.
Performance without amplitude control and Zeeman terms. Here, we start with simple examples with no amplitude control (β = 0 and e r | t=0 = 1) and no Zeeman termsh r ¼ 0. First, we consider antiferromagnetically coupled two DOPOs The ground state energy (E 0 ) is −1 and in the ground state configuration two spins have the opposite sign. We consider the time development starting from μ r = m r = n r = 0. The parametric excitation p depends on time as pðtÞ ¼ 1 þ tanh tþ2 10 (Fig. 2a, Supplementary Note 3). The time step was set to Δt = 2 × 10 −3 . Single run of time development for g 2 = 10 −3 and j = 2 is shown in Fig. 2b, c. With antiferromagneticJ rr 0 matrix, the two DOPOs have the oppositely signed mean amplitudes μ r (Fig. 2b). In Fig. 2c, the variances hΔX 2 r i ¼ m r þ n r þ 1 2 are always smaller than 1 both for r = 1, 2. For two DOPOs, we calculate the success probability P s . If energy calculated from Eq. (1) and spin values σ r ¼ μ r jμ r j at the final time step of the run is equal to the exact ground state energy (E 0 = − 1), the run was counted as success. In Fig. 2d, we show the g 2 dependent success probability obtained from GAPP. For each run, we simulate the time development until t = 4 for the coupling coefficient j = 2. The results from photon number expansion of effective density matrix (EDM) (See Methods and Supplementary Note 4) is shown with black circles. Although the results from GAPP deviate from those of EDM for large g 2 due to the use of Eq. (5) and the Gaussian approximation to neglect higher-order fluctuation products, the deviation from EDM model is smaller than those of Gaussian approximated truncated Wigner (GATW) model 26 (See Methods and Supplementary Note 5) shown with red broken line. Positive-P Gaussian simulation is the better approximation to the density matrix simulation than GATW. This is because we need to truncate higher order terms of Fokker-Planck equation to derive CSDEs in GATW. With the blue line, we show the GAPP results without measurement induced state-reduction, that used ð ∂ρ ∂t Þ S:R: ! j 2 ∑ r ð½â r ;ρâ y r þ h:c:Þ instead of Eq. (5). From this comparison measurement induced state-reduction enhances the performance of measurement feedback CIM.
Next, we considered 3-spin random system, where each independent component ofJ rr 0 ðr > r 0 Þ, i.e.,J 12 ;J 23 ;J 31 are randomly chosen from 21 discrete values (−1, −0.9, ⋯ , 0.9, 1) 26,36 . We prepared up to 10 7 instances and solved each instance only once (10 7 instances for GAPP, GATW and 10 5 instances for EDM). The ground state energy for each instance was obtained by brute-force search. The mean success probabilities P s are shown in Fig. 2e. GAPP is also the better approximation of density matrix simulation than GATW. When g 2 is small, the impact of measurement induced state-reducton was smaller than 2-spin system.
Setting target amplitude. Next, we introduce the model with amplitude control 24 and zero Zeeman terms,h r ¼ 0. We use the following target function for normalized squared amplitude (Supplementary Note 6), This is the approximated form of squared amplitude in solitary DOPO (In Supplementary Note 6, we compare it with the exact steady-state values in ref. 37 .). To see the effectiveness of amplitude control with Eq. (12), we consider the one-dimensional ring with antiferromagnetic nearest-neighbor and next-nearest-neighbor coupling 38 , whereJ rr 0 ¼ À1 if r À r 0 ± 1; ± 2ðmodNÞ, and J rr 0 ¼ 0 for other r À r 0 . The N = 16 ring and the ground state are shown in Fig. 3a. For β = 10, 1, 0.1 and β = 0 (no amplitude control), success probabilities are plotted for various j in Fig. 3b. We used g 2 = 10 −3 in these simulations. Each run was simulated until t = 10 and we used the final measured valueμ r to calculate the energy in Eq. (1), via σ r ¼μ r jμ r j . When the CIM has the ground state energy (E 0 = − 16), the run was counted as success. In Fig. 3b, we obtained higher peak success probability when the strength of amplitude control β is large. Particularly when we used β = 10, the peak value was more than five times larger than the conventional CIM (β = 0) without amplitude control. In Fig. 3c, d, the time development of mean amplitudes is shown for j = 2 and β = 0, 10. The amplitude control with the target function in Eq. (12) solved the problem of amplitude inhomogeneity, and realized the larger success probability even for the frustrated lattice.
Zeeman biased model. We consider a simple model, the anti-ferromagnetically and fully connected four DOPOs (J rr 0 ¼ À1ðr ≠ r 0 Þ;J rr ¼ 0) with non-zero Zeeman terms,h r ≠ 0. The ground state has two up-spins and two down-spins with zero Zeeman term 11 , but with uniform strong bias it has a ferromagnetic order (Fig. 4a). We set the strength of Zeeman bias as the parameter,h r ¼ ζ. In Fig. 4b, we show the ζ dependent success probability for β = 10 and 0. We fixed the coupling and nonlinear saturation to j = 2, g 2 = 10 −3 . With amplitude control (β = 10), the success probability is increased from that with no amplitude control (β = 0). However, even with β = 10, the success probability has dropped at two spots, around ζ = 1 and around ζ = 3. This is understood in terms of the ground state transition for varying ζ. As shown in Fig. 4c, when ζ ∈ (1, 3) the ground state has three up-spins and one down-spin. There is a dip in success probability when ζ is slightly smaller than 1, where there is only small difference between the ground state energy (E 0 ) with two up-spins and two down-spins and the first excited state energy (E 1 ) with three up-spins and one down-spin. The first excitation state is excited by the fluctuations because of small energy difference. On the other hand, at ζ = 1 both states (two up-spins and two down-spins, three up-spins and one down-spin) are ground states and are counted as success. Therefore, the success probability was almost one for both β = 10, 0. In Fig. 4d, we show the probability for the system in the first excited state (we call P(E 1 )), and the energy gap between ground state (E 0 ) and first excitation state (E 1 ). When E 1 − E 0 is close to zero, the probability to find the first excitation state P(E 1 ) increases.  Complete graph with random weights and random Zeeman terms. Next, we consider the system of N = 16 complete graph with a random weight matrix. Each component of J rr 0 ðr > r 0 Þ andh r is randomly chosen from 21 discrete values (−1, − 0.9, ⋯ , 0.9, 1). We show three examples of random instances in Supplementary Note 7. In the numerical simulation, we calculated mean success probability P s by generating 10 5 problem instances and performing the simulation only once for each instance. The success was judged at t = 20 by comparison between the CIM energy calculated from σ r ¼μ r =jμ r j and the ground state energy (E 0 ) obtained from the brute-force search. For g 2 = 10 −3 , we show the j-dependent success probability for β = 10, 1, 0.1, 0 (Fig. 5a). The peak success probability was the largest when we use the strongest amplitude control β = 10. However, as different from the results in Fig. 3b, β = 0.1 had slightly better peak success probability than β = 1. In this case, the β dependence has two peaks as shown in Fig. 5b.
Next, we compare our approach which uses target amplitude to realize Zeeman terms (Eq. (7)) with two models. The first one is conventional CIM 14,15,17 , that has no amplitude control and no Zeeman terms (β ¼ 0;h r ¼ 0). Each component ofJ rr 0 ðr > r 0 Þ is chosen from 21 discrete uniform values (−1, −0.9, ⋯ , 1). The second one is the model using mean field to represent Zeeman terms 20,39 , dμ r dt In Fig. 5c, we compare our realization of Zeeman terms, the realization using mean absolute amplitude 20,39 (both have β = 10 amplitude control), and conventional CIM with no amplitude control and no Zeeman terms. For 16-spin randomJ rr 0 matrix, the simulation with our method had the slightly better performance than the one with mean absolute amplitude. Under amplitude control, two realizations of Zeeman terms had the  . c Comparison of j-dependence for β = 10 coherent Ising machine (CIM) and previously considered models, conventional CIM which doesn't have amplitude control and Zeeman terms, and β = 10 amplitude controlled CIM with Zeeman terms provided by Eq. (13). d Dependence on the relative size of Zeeman terms ζ (j = 2). e The success probability (P s ) and the probability to find first excitation (P(E 1 )) of one random instance (I1). f The excitation energy of the first excitation state from the ground state (E 1 − E 0 ) in I1-instance. better peak performance than the conventional CIM. The comparison with auxiliary spin method to realize Zeeman terms 7,21 is shown in Supplementary Note 8.
In Fig. 5d, we calculated the mean success probability when we change the relative magnitudes of Zeeman terms, by multiplying ζ toh r , i.e.,h r =ζ 2 ½À1; À0:9; Á Á Á ; 0:9; 1. For certain values of ζ, particularly when ζ is an integer value, there are small peaks of performance. As we saw in Fig. 4, if the ground states degenerate accidentally, the success probability increases because all degenerate states with ground state energy E 0 are counted as success. The small peaks of performance are the results of the fact that degeneration occurs more likely for integer ζ than for fractional ζ. The characterstics for large ζ explicitly depended on whether we used the amplitude control or not. With amplitude control the mean success probability increased for large ζ. On the other hand without amplitude control (β = 0) the success probability decreased. In Fig. 5e, f, we show the ζ dependent success probability P s , and characteristics related to the first excitation state with energy E 1 for one random instance (I1), whoseJ rr 0 matrix andh r vector are shown in Supplementary Note 7. As shown in Fig. 5e, the dip of the success probability P s is closely related to the increase of probability in finding the first excitation state P(E 1 ), which increases when the energy gap (E 1 − E 0 ) is close to zero (Fig. 5f).

Conclusion
We show the numerical simulation results of measurement feedback CIM modeled as GAPP model, with correction of amplitude inhomogeneity whose target amplitude is also used for implementing Zeeman terms. For N = 16 frustrated system with nearest-neighbor and next-nearest-neighbor coupling and N = 16 complete graphs with random weights and Zeeman terms, the suggested CIM with amplitude control has the improved performance than the CIM without amplitude control. Particularly as the magnitudes of Zeeman termsh r increased with fixedJ rr 0 , the success probability increased with amplitude control, whereas the success probability decreased without amplitude control. Studying the dependence on saturation coefficient g 2 is also important (Supplementary Note 9), since large g 2 is expected to reduce energy to solution and GAPP is more accurate than GATW model 26 in such a region. The method in the paper will be effective in the application for compressed sensing 22 . Recently, the experimental measurement feedback CIM with N~100, 000 spins was achieved 40 . GAPP will be applied to such a large N system. For current experimental CIM, where the round-trip time Δt is not much smaller than the photon lifetime, the discrete model 41 can be more accurate method to simulate experimental system.

Methods
Effective density matrix (EDM) equation. To simulate measurement feedback CIM with density matrix equation, we add the terms to represent effective increase of single mode fluctuation that occurs when the fluctuation (∝W r ) in Eq. (4) is ensemble-averaged (related to the last term of Eq. (8) in ref. 28  . The total density matrixρ is written as direct products of the density matrices of each DOPO (we useρ ðrÞ to represent r-th DOPO), since DOPOs are coupled via mean amplitudes hâ r þâ y r i and real random numbers W r . For r-th DOPO, density matrix can be expanded using photon number states jN r iðN r ¼ 0; 1; Á Á ÁÞ, asρ ðrÞ ¼ ∑ N r N 0 r ρ ðrÞ  Fig. 2d, e, each run started from the vacuum state where only ρ ðrÞ 0;0 ¼ 1 components are non-zero and a run was counted as success when the energy calculated from σ r ¼ μ r jμ r j and Eq. (1) is identical to the ground state energy from the brute-force search. Here, μ r ¼ hâ r þâ y r i 2 were calculated as 1 2 ∑ N r ffiffiffiffiffiffiffiffiffiffiffiffiffiffi N r þ 1 p ðρ ðrÞ Þ N r þ1;N r þ ρ ðrÞ N r ;N r þ1 Þ.

Equations of Gaussian approximated truncated Wigner (GATW) simulation.
In GATW simulation, we used following two equations to represent r-th DOPO.
dμ r dt ¼ Àð1 À p þ jÞμ r À g 2 μ 3 r þ j ∑ r 0J rr 0μ r 0 þ ffi ffi j p V r À 1 2 W r ; ð15Þ dV r dt ¼ À2ð1 À p þ jÞV r À 6g 2 μ 2 r V r À 2j V r À Here, μ r ¼ hX r i= ffiffi ffi 2 p are the mean amplitudes, and V r ¼ hΔX 2 r i are the inphase variances. The derivations of these equations are provided in Supplementary Note 5. In Fig. 2d, e, the time development started from the vacuum state, μ r ¼ 0; V r ¼ 1 2 , and a run was counted as success when the energy calculated from σ r ¼