Network of TimeMultiplexed Optical Parametric Oscillators as a Coherent Ising Machine
Finding the ground states of the Ising Hamiltonian [1] maps to various combinatorial optimization problems in biology, medicine, wireless communications, artificial intelligence, and social network. So far no efficient classical and quantum algorithm is known for these problems, and intensive research is focused on creating physical systems  Ising machines  capable of finding the absolute or approximate ground states of the Ising Hamiltonian [2, 3, 4, 5, 6]. Here we report a novel Ising machine using a network of degenerate optical parametric oscillators (OPOs). Spins are represented with abovethreshold binary phases of the OPOs and the Ising couplings are realized by mutual injections [7]. The network is implemented in a single OPO ring cavity with multiple trains of femtosecond pulses and configurable mutual couplings, and operates at room temperature. We programed the smallest nondeterministic polynomial time (NP) hard Ising problem on the machine, and in 1000 runs of the machine no computational error was detected.
Many important combinatorial optimization problems belong to the NPhard or NPcomplete classes, and solving them efficiently has been beyond the reach of classical and quantum computers. Currently, simulated annealing [2] and various approximate algorithms [8, 9, 10] are among the widely used methods in classical digital computers. Recently, quantum annealing [3, 4, 6] and adiabatic quantum computation [5] have been studied intensively for solving those NP problems. Even though relatively large machines are recently implemented to perform quantum annealing [6], their comprehensive potentials are yet to be explored [12, 11], and the quest for a novel computing machine to solve largescale NPhard problems continues.
Recently, degenerate OPOs pumped by femtosecond laser pulses have been shown useful for a wide range of classical and quantum applications. While operated above their oscillation threshold, they are used in generation of broadband midinfrared frequency combs which are intrinsically phase and frequency locked to the pump [13]. Below threshold, they are employed for realization of quantum networks using wavelength division multiplexing [14]. Intriguingly, these OPOs demonstrate nonequilibrium phase transition at threshold, which has been previously exploited for optical quantum random number generation [15]. Here we explore the computational capability of this phase transition in a coherent network of timedivisionmultiplexed femtosecond OPOs.
The Hamiltonian of an Ising model with spins and without an external field is given by [16]:
(1) 
with being the coupling between the and spins, and and representing the projections of the spins, the eigenvalues of which are either +1 or 1. To implement an artificial Ising spin system, elements with binary degrees of freedom and configurable couplings are required. This has been demonstrated experimentally in cold ion or atom system [17, 18] with practical difficulties in scalability. Artificial Ising spin systems based on a network of injection locked lasers are also proposed [20, 19].
Degenerate OPOs are particularly suitable to represent Ising spins because of binary phase operation above threshold [21]. Fig. 1a illustrates the operation of a degenerate OPO for below and above oscillation threshold in the inphase and quadraturephase coordinates. Below the threshold, the signal field is a squeezed vacuum state [22], and by increasing the pump field it undergoes spontaneous symmetry breaking around the oscillation threshold [23]. This nonequilibrium phase transition leads to oscillation in one of the two phase states ( or ) [21], which are exploited to represent an Ising spin (). The couplings between the spins () are realized by mutual injections of the signal fields of the and OPOs. The resulting OPO network has a phasestate dependent photon decay rate corresponding to the energy landscape of the original Ising Hamiltonian [7], i.e. different spin configurations (phase states) have different total photon losses in the network, which is illustrated in Fig. 1b.
The search for the ground states is conducted by gradually raising the gain via increasing the pump field, as illustrated in Fig. 1b. As the parametric gain approaches the lowest possible loss, the network goes through the OPO phase transition, and as a result it is expected to oscillate in one of the exact or approximate ground states. Further increase in the pump pins the gain to the threshold value and thus erroneous oscillation in the excited states continues to be suppressed. Although the below or abovethreshold behavior of an OPO network can be efficiently simulated on a classical computer, no efficient algorithm is yet known for comprehensive simulation of the phase transition point.
We use a numerical model based on cnumber Langevin equations, which are compatible with fundamental quantum mechanical master equations or FokkerPlanck equations. We simulate OPO networks of up to 20000 OPOs to solve some of the NPhard maximum cut (MAXCUT) problems. For these problems, the OPO network outperforms one of the wellknown approximate algorithms in terms of the computation time and the accuracy of the results (see Supplementary Information). However, despite these promising numerical results, no mathematical proof exists to guarantee superior performance of the OPO network for all Ising problems, and the possibility of failure for some instances exists.
For the experiment, we exploit the novel scheme of “timedivisionmultiplexing”, in which a single ring resonator accommodates many Ising spins simultaneously, which are free from mismatch and phase decoherence noise. Figure 2a shows the experimental setup comprising a 4OPO system. The ring resonator has a roundtrip time () of four times the pump pulse repetition period (), as shown in Fig. 2b, and therefore the pulses represent four temporally separated independent OPOs. The couplings between these independent but identical OPOs are realized using three pairs of output and input couplers in the resonator paths. In each pair, 4% of the intracavity power is outcoupled and precisely delayed in a freespace delay line; each delay is an integer multiple of the repetition period, and 4% of the outcoupled light is injected back to the resonator. This corresponds to 4% of field coupling between the OPOs. Figure 2b depicts how these delay lines provide all the possible couplings () among the four temporally separated OPOs.
To measure the phase states of the 4OPO system, the output is sent to a onebit delay Michelson interferometer as shown in Fig. 2a. The delay difference of the arms is locked to the repetition period of the system (). Therefore, the interferometer measures the differential phase between adjacent OPOs, i.e. the output is the interference of OPO and OPO, and changes with the time slot of measurement. We use fast and slow detectors at the output of the interferometer and in Supplementary Information we show the expected outputs for different phase states.
The 4OPO system is pumped with a femtosecond mode locked laser with a central wavelength of 1045 nm, and a repetition period of 4 ns. It operates at degeneracy and the detailed characterization is presented in Supplementary Information. The pump field is gradually increased from below to about two times above the oscillation threshold over a time of 180 s, while the cavity photon lifetime is 60 ns.
First, we present the measurements using a slow detector in Fig. 3. When the couplings are blocked, the OPOs are expected to operate independently, and therefore a random uniformly distributed selection of the phase states are expected each time the OPO is restarted. A sample of the interferometer output for this case is shown in Fig. 3a. The output toggles between the three distinct intensity levels (0, , and ). Assuming all the states are equally probable to occur, corresponds to 12 phase states and and correspond to 2 different phase states (Supplementary Information). Indeed the frequency of these events are 6:1:1, as shown in Fig. 3a.
Having inphase (outofphase) couplings among the adjacent OPOs corresponds to ferromagnetic (antiferromagnetic) couplings in an Ising spin ring with the expected outcome of aligned (antialigned) spins. In the experiment, when the shortest delay line is set to inphase (outofphase) coupling, the OPOs operate at the same (alternating) phase states. We show this behavior in Fig. 3b where the phase of delay 1 is swept and the other delays are blocked. The sweeping speed is slow compared with the 1 kHz restarting frequency of the OPOs. The output of the slow detector is always at when phase of delay 1 is around zero (). This corresponds to or phase states confirming that inphase coupling between OPO and OPO leads to all OPOs oscillating in the same phase state. When the phase of delay 1 is around 180, i.e. the field of OPO is shifted by and injected to OPO, the output intensity remains at level, meaning that the OPOs are either in or . As indicated in the plot, the coupled OPOs tolerate a coupling phase deviation of at least 30 around these two points of inphase or outofphase couplings. This regenerative behavior is due to the phase sensitive nature of the parametric gain in the degenerate OPO and makes the system immune to the phase noise in the environment. Similar behavior is observed for the other delay lines and the detailed results are presented in Supplementary Information.
We focus on implementation of an NPhard MAXCUT problem, which is mapped to an Ising problem [24] and corresponds to a frustrated Ising spin system [17]. The problem is to find a subset of vertices in a cubic graph, in which all vertices have three edges, such that the number of edges between the subset and its complementary subset is maximized. For a cubic graph with 4 vertices, which has 6 edges, any subset choice with two vertices is an answer of the MAXCUT problem (or MAX2SAT problem), and any other subset choice is not. To program this problem on the OPO Ising machine, all the mutual couplings between the OPOs are set to out of phase, i.e. . The total photon decay rate of the network is lowest for the answers of the MAXCUT problem [7] corresponding to two OPOs in and two OPOs in , hence upon each gradual pumping the 4OPO network is expected to oscillate in such a phase state.
The phases of the three delays are locked to using the residual pump interference signal at the transmission port of the input couplers, and a sample of the interferometer output is shown in Fig. 3c. It matches the expected outcome of being either or with the rate of twice the rate of . However, since corresponds to both answer and nonanswer states, it is necessary to use a fast detector to confirm that the nonanswer states do not occur. Other configurations for the coupling phases (corresponding to other Ising problems) are also investigated using the slow detector and the results follow the expected outcomes (see Supplementary Information).
When a fast detector is used at the interferometer output and all the couplings are blocked, after the 4OPO system is turned on, one of the four possible pulse patterns (Supplementary Information) is detected as shown in Fig. 4a. In these patterns, the pulses are separated by 4 ns –the repetition period of the pump– and each pulse has either a low or a high intensity corresponding to destructive or constructive interference of the consecutive OPOs, respectively. The lowlevel pulses are nonzero because of the diffraction mismatch of the interferometer arms. For the case of no coupling, the 4OPO system is restarted 1000 times, and the histogram of the eight phase states is shown in Fig. 4b. For each restart, the pulse train is detected and the corresponding phase state is inferred. The result indicates uniform distribution of the phase states confirming that four temporally separated OPOs are operating independently in the same resonator.
To verify that the OPO network is capable of solving the NPhard MAXCUT problem, the histogram of the OPO states is measured for the case of all the delays locked to the phase of for 1000 trials. The result is depicted in Fig. 4c. No excited phase state is detected, and the distribution of the answer states is close to uniform. The error rate of computation is less than which is limited by the length of measurements.
In conclusion, we demonstrate a coherent network of timedivisionmultiplexed degenerate OPOs that is capable of solving the NPhard MAXCUT problem for =4. The practical scalability enabled by multipulse operation in a long ring cavity, simple mapping of the OPO bistable phase states to the Ising spin, and intriguing quantum phase transition at the threshold can make the OPO network suitable for obtaining approximate solutions for largescale NPhard Ising problems with reasonable accuracy and speed. Time division multiplexing enables increasing the size of the network by increasing the cavity roundtrip time (), and implementing delay lines, with each delay corresponding to one coupling at a given time slot (see Supplementary Information for a practical implementation). To program a largescale machine to arbitrary Ising problems, intensity modulators are needed to turn on and off each delay as well as phase modulators to choose either inphase or outofphase couplings. Exploiting optical fiber technologies and planar light wave circuits can enable implementation of compact large Ising machines.
Acknowledgements
The authors would like to thank S. E. Harris, H. Mabuchi, M. Armen, S. Utsunomiya, S. Tamate, K. Yan, and Y. Haribara, for their useful discussions, and K. Ingold, C. W. Rudy, C. Langrock, and K. Urbanek for their experimental supports. The work is supported by the FIRST Quantum Information Processing project.
Contributions
A.M. and Y.Y. conceived the idea and designed the experiment. A.M. and K.T. carried out the experiment. Z.W. performed the numerical simulations. Y. Y. and R.L.B. guided the work. A.M. wrote the manuscript with input from all authors.
Supplementary Information
1 Principle of Operation
The proposed computational concept is in sharp contrast with the classical and quantum annealing techniques. The cartoon of Fig. S1 illustrates the challenge in classical and quantum annealing to go from a metastable excited state to the ground state through multiple thermal hopping or quantum mechanical tunneling through many metastable states. Classical simulated annealing employs a downward vertical search, in which the temperature is repeatedly decreased and increased until the ground state is found. Quantum annealing exerts a horizontal search in the energy landscape with quantum tunneling. Therefore, with these methods, the computational time of finding the ground state increases with the increase in the number of metastable excited states or local minima. In contrast, the OPObased Ising machine searches for the ground state in an upward direction. The total energy, i.e. the ordinate in Fig. S1, is now replaced by the network loss. The ground state (optimum solution) has a minimum loss as shown in the figure. If we put a parametric gain (G) into such a network and increase it gradually, the first touch to the network loss happens at the ground state, which results in the singlemode oscillation of the ground state spin configuration. At the pump rate above this threshold point, the parametric gain is clamped at the same value of due to nonlinear gain saturation, so that all the other modes including local minima stay under the oscillation threshold. If we use the terminology of “negative temperature” to represent the parametric gain, the mentioned upward search corresponds to the heating process from for zero gain toward for high gain. In this sense, the OPO machine is a “heating machine” while the classical simulated annealing is a “cooling machine.”
2 Theoretical Modeling and Numerical Test of OPO Network
To simulate a network of degenerate OPOs, we can start with the quantum mechanical FokkerPlanck equation (QFPE) for a single OPO using the generalized Prepresentation [25]:
(S1) 
Here and are the normalized eigenvalues for the offdiagonal coherentstate expansion, , of the density matrix, and are the signal and pump photon decay rates, is the oscillation field amplitude at a normalized pump rate , is the normalized time, is the pump field amplitude, is the threshold pump amplitude. The average amplitudes of inphase and quadraturephase components of the signal wave are obtained by
(S2) 
Equation 2 can be cast into the cnumber Langevin equation (CLGE) for the inphase and quadraturephase components of the signal field via the KramersMoyal expansion [26]:
(S3) 
where and , and and are two independent Gaussian noise processes which represent the incident vacuum fluctuations at signal frequency and pump frequency .
The equivalence of the QFPE (2) and the CLGE (2) can be confirmed by comparing the squeezing and antisqueezing characteristics of the two quadrature components using
(S4) 
for the QFPE and
(S5) 
for the CLGE. As shown in Fig. 1 in [7], the degrees of squeezing and antisqueezing obtained by the two methods completely agree at a pump rate across the threshold ().
The QFPE and CLGE for an injectionlocked laser oscillator [27] can be extended to the mutually coupled degenerate OPOs using equation 2 or 2. The resulting CLGEs for a network of degenerate OPOs are given by
(S6) 
Here and are the normalized amplitudes of two quadrature components of the th OPO which corresponds to and in Eq. 2.
Performance of the proposed network of degenerate OPOs as an Ising machine is tested against the NPhard MAXCUT problems on cubic graphs for to and on random graphs for to . For a graph with vertices, the CLGEs are solved by the DormandPrince method as the differential equation solver [28], in which adaptive integration strength is introduced by evaluating the local truncation error.
Figure S2 shows the normalized buildup time when the OPO network reaches the steady state oscillation conditions after an abovethreshold pump rate () is turned on versus the graph order . We have numerically tested all graphs, for instance a total number of 510489 graphs are studied for =20. Most of the buildup time (up to 99% of all graphs) is independent of the graph order and is on the order of 100, as shown in Fig. S2. Only slight increase is observed for the worst case as shown in green triangles. Therefore, an actual computational time is determined mainly by the success probability for obtaining a ground state. Since the proposed OPO network is a stochastic machine driven by quantum noise, the success probability is always smaller than one.
Order  
4  6  8  10  12  14  16  20  
0.93  1.00  0.41  0.54  0.52  0.37  0.33  0.11  
1.05  1.30  1.30  1.30  1.00  0.85  0.85  0.82  
1.00  1.00  0.70  0.74  1.00  1.00  1.00  0.74  
Table S1 summarizes the performance of the OPO network in solving MAXCUT problems on cubic graphs. Here, denotes the worstcase success probabilities at a fixed pump rate and coupling coefficient , denotes the optimal pump rate for each worstcase instance, at which the optimum success probability is achieved under the same coupling coefficient . The success probability at the optimum pump rate for the worst instance is independent of the graph order and ranges from .
Graph  

G1  800  19176  12083  0.9591  0.9516  498.82 
G6  800  19176  2656  0.9559  0.9506  471.06 
G11  800  1600  629  0.9384  0.9254  406.24 
G14  800  4694  3191  0.9367  0.9274  498.26 
G18  800  4694  1166  0.9308  0.9223  430.24 
G22  2000  19990  14136  0.9349  0.9277  768.34 
G27  2000  19990  4141  0.9321  0.9270  780.18 
G32  2000  4000  1567  0.9328  0.9260  467.42 
G35  2000  11778  8014  0.9264  0.9202  602.34 
G39  2000  11778  2877  0.9214  0.9152  539.9 
G43  1000  9990  7032  0.9373  0.9309  542.92 
G48  3000  6000  6000  0.9463  0.9292  762.34 
G51  1000  5909  4006  0.9333  0.9242  491.68 
G55  5000  12498  11039  0.9070  0.9009  903.46 
G57  5000  10000  3885  0.9305  0.9259  648.72 
G59  5000  29570  7312  0.9114  0.9074  583.54 
G60  7000  17148  15222  0.9037  0.8995  918.98 
G62  7000  14000  5431  0.9295  0.9256  719.74 
G64  7000  41459  10466  0.9129  0.9092  666.16 
G65  8000  16000  6206  0.9284  0.9252  757.28 
G66  9000  18000  7077  0.9285  0.9251  786.9 
G67  10000  20000  7744  0.9285  0.9260  732.14 
G70  10000  9999  9863  0.9433  0.9379  687.28 
G72  10000  20000  7809  0.9284  0.9256  748.3 
G77  14000  28000  11046  0.9281  0.9256  842.24 
G81  20000  40000  15656  0.9268  0.9250  980.54 
The performance of the OPO network in solving the MAXCUT problems has also been examined on 71 benchmark instances of the socalled Gset graphs when and . These instances are randomly constructed by a machineindependent graph generator written by G. Rinaldi, with the number of vertices ranging from 800 to 20000, edge density from 0.02% to 6%, and geometry from random, almost planar to toroidal. The outcomes of running the OPO network 100 times for sample Gset graphs are summarized in Table. S2. It can be easily seen that both the best and average outputs of the OPO network are about 2  6% better than the 0.878performance guarantee of the celebrated GoemansWilliamson algorithm based on semidefinite programming (SDP) [29]. Since the differences between the best and the average values are within 1% for most of the instances, reasonable performance is expected for the OPO network even in a single run, which makes the OPO network favorable for applications when response time is the utmost priority. In addition, there is further room to improve the performance, for example by applying local improvement to the raw outcomes of the OPO network and operating the OPO network under optimum pump rate of and coupling strength of .
The average computational time of the OPO network in solving the MAXCUT problem on the worstcase cubic and Gset graphs is displayed in Fig. S3. The growth of the computation time fits well to a sublinear function . Computational complexities of bestknown algorithms for solving the SDP in the GoemansWilliamson algorithm are also plotted. If a graph with vertices and nodes is regular, the SDP can be approximately solved in almost linear time as using the matrix multiplicative weights method [30], where represents the accuracy of the obtained solution. This behavior is shown by the red solid line in Fig. S3. However, slower algorithms are required for general graphs. If the edge weights of the graph are all nonnegative, the fastest algorithm runs in time based on a Lagrangian relaxationbased method [31]. This computational time is plotted by the pink solid line. For graphs with both positive and negative edge weights, the SDP is commonly solved using the interiorpoint method which scales as [32]. This general computational time for SDP is shown by the blue dashed line. Since the OPO network is applicable to all types of graphs, the sublinear scaling of the computation time gives it a huge advantage over the SDP algorithm in solving largescale instances. For instance, the OPO network outputs a solution for the graph G81 with 20000 vertices and 40000 edges in a normalized time of 980, which corresponds to the actual time of sec, since the photon lifetime in the largescale fiberbased OPO network , which can handle this size of problem, is sec (see the last section in this Supplementary Information). This computational time is compared to the time required to run the SDP (interiorpoint method) using a 1.7 GHz core i7 machine of about sec, which is seven orders of magnitude larger than sec for the OPO network.
3 Measurement of Phase States
For the 4OPO system, the onebit delay interferometer is used for measuring the phase states. Table S3 shows the phase states and corresponding output pulse trains in the first and second columns. Since the complementary phase states result in the same output, for the 16 possible phase states 8 different pulse trains can occur. In the measurements with a fast detector, no time reference is used, therefore only four different pulse patterns can be detected. If a slow detector is used, the output will be proportional to the number of ones in the pulse train. Therefore, there will be three distinct output average intensity levels namely , , and , as shown in the third column of the table.
Interferometer  
State  Pulse Train  Slow Detector 
[1111]  
[0110]  
[0011]  
[1010]  
[1001]  
[0000]  
[0101]  
[1100]  
[1100]  
[0101]  
[0000]  
[1001]  
[1010]  
[0011]  
[0110]  
[1111]  

4 Details of Experimental Setup
The OPO design shares similarities with the experiment reported in [33]. The ring resonator of the OPO illustrated in Fig. 1b, has a round trip time of 16 ns (a perimeter of 4.8 m). The setup has two more flat mirrors than the schematic (corresponding to a folded bow tie configuration). All the flat mirrors, except M are gold coated with enhancement dielectric coatings at 2 m. One of the flat gold mirrors is placed on a translation stage with piezoelectric actuator (PZT). The dielectric mirror (M) has a coating which is antireflective at the pump wavelength with less than 0.2% reflection, and is highly reflective (99%) from 1.8 m to 2.4 m. The curved mirrors (M and M) have 50mm radius of curvature and are unprotected gold coated mirrors. The angle of incidence on these mirrors is 4, which is chosen to compensate the astigmatism introduced by the Brewstercut nonlinear crystal, and results in 1 mm of cavity stability range for the spacing between the curved mirrors. The signal beam has a waist radius of 8.3 m (1/e intensity) at the center of the crystal.
The 1mm long, Brewstercut, MgO:PPLN crystal has a poling period of 31.254 m, which is designed to provide degenerate parametric gain for a pump at 1035 nm with type 0 phase matching (ee+e) at 373 K temperature. The crystal operates at room temperature in the OPO, and even though the phase matching condition is not optimal for the pump (centered at 1045 nm), degenerate operation is achieved by length tuning of the cavity.
Input and output coupling of the signal are achieved with 2m thick nitrocellulose pellicles to avoid dispersion in the cavity and etalon effects. In Fig. 1b, the three pairs of “OC” and “IC” are uncoated pellicles (with % power reflection ), the “OC” for the main output is coated (with 15% power reflection at 2090 nm). The beam splitter in the interferometer (BS) is the same coated pellicle. For stabilizing the OPO cavity, another uncoated pellicle is used as an output coupler in the resonator (not shown in the schematic).
The pump is a freerunning modelocked Ybdoped fiber laser (Menlo Systems Orange) producing 80 fs pulses centered at 1045 nm with a repetition rate of 250 MHz, and maximum average power of 1 W. The filter is a long pass filter at 1850 nm on a Ge substrate to eliminate the pump and transmit the signal.
Gradual pumping is achieved by the chopper, as it causes a rise time (1090% power) of 180 s for introducing the pump. The cavity photon lifetime for the signal is estimated to be 60 ns, and the network is pumped 2.2 times above threshold.
4.1 Servo Loops
Five feedback servo controllers are used to stabilize the length of the cavity, the phase of the delay lines, and the armlength difference of the interferometer. All the controllers are based on “ditherandlock” scheme, where a slight modulation (less than 10 nm amplitude at a frequency between 5 and 20 kHz) is applied to a fast PZT, and the error signal is generated electronically by mixing the detector output and the modulated signal [13]. Identical electronic circuits are used with a controller 3dB bandwidth of 10 Hz.
For the delay lines, the interference of the pump at the other port of the input couplers are used as the input of the controller, and the controller locks the length to achieve destructive interference on the detector, which results in constructive interference on the other port that enters the cavity. The arm length difference of the interferometer is also locked similarly. No phase stabilization is required for the path from the OPO to the interferometer since all the OPO pulses experience the same path and phase change.
4.2 FreeRunning Pump
Here we show that the servo controllers used in the experiment suffice for implementation of the Ising machine and no stabilization on the pump is required. Slow changes (within the response time of the controller) in the carrierenvelop offset frequency () and repetition rate () of the pump do not affect the operation of the Ising machine. Smooth changes in of the pump is intrinsically transferred to the signal since signal pulses are generated from pump pulses. However, the effects of changes in on the Ising machine require taking into account the intrinsic phase locking of the degenerate OPO as well as the role of servo controllers.
The primary task of the servo controller of the OPO is to maximize the output power by matching the roundtrip phase in the resonator to the pulse to pulse phase slip of the pump (). The pulse to pulse phase slip is related to by:
(S7) 
Assuming the carrier fields for the pump and the signal pulses are defined as: , and , respectively, the phasesensitive gain dictates [13]:
(S8) 
Therefore, if the carrier phase of the pump changes by from one pulse to next, for a single OPO (), the pulse to pulse phase slip of the signal follows that by:
(S9) 
This means that the phase slip of the signal pulses is locked to the phase slip of the pump with a factor of half, which consequently means the of the pump and signal are locked [13]; the servo loop provides feedback to the cavity to follow this phase slip and maximize the output power. A similar behavior also happens for a doublyresonant OPO operating away from degeneracy, with a different ratio between the s of the pump and signal [34].
For OPOs in the cavity (i.e. ), when all OPOs are in the same phase state, the pulse to pulse phase slip in the pump transfers to the OPO to OPO phase slip by a factor of half. Changing phase state from one OPO to another simply means adding to the phase slip. When a delay line is locked to the top of the interference fringe of the pump pulses, the phase change provided by the delay line at the pump wavelength compensates the pulse to pulse phase slip of pump, i.e. . At the signal wavelength, because this phase change in the delay line is half (), which means that the servo controller compensates the OPO to OPO phase slip resulting from the .
Locking a delay line to top of the fringe of pump pulses corresponds to having either or phase change for the signal. This is also true for the interferometer. In the experiments, for all the servo controllers, we were able to precisely tune the length from one fringe to the next. This gave us the ability to try different configurations and find the desired coupling phases.
5 OPO Characterization
When no out coupler is used in the cavity, the OPO has a threshold of 6 mW of pump average power. With all the output and input couplers in the cavity, the threshold reaches 135 mW. Oscillation at degeneracy and away from degeneracy can be achieved depending on the cavity length. The OPO is pumped with 290 mW and the main output of the OPO has 15 mW of average power at degeneracy centered at 2090 nm with the spectrum shown in Fig. S4a. The interferometric autocorrelation of the signal pulses is shown in Fig. S4b suggesting a pulse length of 85 fs. The spatial profile of the output beam is very close to Gaussian as shown in Fig. S4c with a radius of 1 mm (at intensity). The average power of the signal in the delay lines are 2 mW, and the intracavity power is estimated to be 100 mW.
6 Extended SlowDetector Results
In Fig. S5 we show the results obtained using a slow detector for different combinations of couplings. Fig. S5 ac are obtained by scanning the phase of one delay line while the other delay lines are blocked. Delay 1 and 3 have similar effects, because they couple adjacent OPOs but in different directions (Fig. 2b). As shown in Fig. S5 a and c, inphase coupling by these delays results in the same phase state for all OPOs, and consequently highintensity interferometer output (); and outofphase coupling results in alternating phase states and consequently lowintensity interferometer output ().
In Fig. S5b we show the interferometer output while the phase of delay 2 is scanned. When the coupling is inphase, OPO 1 and 3 have the same phase state, and OPO 2 and 4 oscillate in the same phase state. However, these two pairs can either be the same or different, and therefore the output would be either or 0, as shown in Fig. S5b around phase of zero. On the other hand, outofphase coupling of delay 2 results in constant output as shown in the same plot. Regenerative behavior of the OPO and its insensitivity to a wide range of phase change in the couplings are observed in these three plots.
When the network is configured to the MAXCUT problem, we scanned the phase of the delays one by one, and the results are shown in Fig. S5df. Different delay phase configurations and the expected outcomes are shown in Table S4, where the last row corresponds to the MAXCUT problem with all antiferromagnetic couplings, and for each of the other rows one of the delay phases is different. For each plot in Fig. S5df, the center of the plot, where the phase of the scanned delay is , corresponds to the antiferromagnetic MAXCUT problem. The outputs follow the expected outcome.
Phase of Couplings  Expected Phase  
[D1, D2, D3]  States  Slow Detector 
100% in  

100% in  

100% in  

66.7% in ,  
33.3% in  

7 Scalability
Time division multiplexing facilitates scalability of the OPO network, and it benefits from intrinsically identical nodes in the network. It is also worth noting that for an Ising problem with sites, the number of possible couplings () is . However, in the timedivisionmultiplexed network only delay lines are required for realizing these couplings, which means that the physical size of the machine scales linearly with .
A network of OPOs can be realized in a single ring resonator with a roundtrip time of ( is the pulsetopulse interval), and constructing delay lines. Schematic of fiberbased implementation of such a machine is illustrated in Fig. S6. To avoid effects of nonlinearities and dispersion in optical fibers, picosecond pump pulses can be used in a long resonator and long delay lines comprising optical fiber components. As an example, for a pump with 10GHz repetition frequency ( ps), a resonator with 200 m of optical fiber results in 10000 temporally separated OPOs. An expected photon lifetime of such a fiberbased OPO network is about s, which promises a reasonably fast computational time for a MAXCUT problem with . The main challenge is stabilizing the phases of all these fiber links. Development of extremely lownoise phasestabilized long (100km) optical fibers [35] promises overcoming this challenge using the existing technologies. Moreover, the regenerative behavior of the degenerate OPO (as shown in Fig. 3b) suggests that the OPObased Ising machine can tolerate relatively large phase noise in the couplings.
Similar to the 4OPO network, each delay line provides a delay of equal to an integer multiple of the repetition period (), and is responsible for multiple of the Ising coupling terms in the form of . In one delay line, each of these couplings happens at one time slot, and one can use electrooptic phase and amplitude modulators (EOM) to synchronously switch the delay on and off depending on whether the corresponding coupling term is zero or not. This can be extended to synchronously controlling the phases and strengths of the couplings through the delay lines, and hence programming any arbitrary Ising problem on the machine.
References
 [1] F. Barahona, “On the computational complexity of Ising spin glass models,” J. Phys. A: Math. Gen. 15, 3241 (1982).
 [2] S. Kirkpatrick, D. Gelatt Jr., and M. P. Vecchi, “Optimization by simulated annealing,” Science 220, 4598 (1983).
 [3] T. Kadowaki, and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E 58, 5355 (1998).
 [4] G. E. Santoro, R. Martoak, E. Tosatti, and R. Car, “Theory of quantum annealing of an Ising spin glass,” Science 295, 2427 (2002).
 [5] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, “A quantum adiabatic evolution algorithm applied to random instances of an NPcomplete problem,” Science 292, 472 (2001).
 [6] M. W. Johnson et al, “Quantum annealing with manufactured spins,” Nature 473, 194 (2011).
 [7] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, “Coherent Ising machine based on degenerate optical parametric oscillators,” Phys. Rev. A. 88, 063853 (2013).
 [8] C. H. Papadimitriou, and K, Steiglitz, Combinatorial optimization: algorithms and complexity (Courier Dover Publications, 1998).
 [9] D. S. Hochbaum, Approximation algorithms for NPhard problems (PWS Publishing Co., 1996).
 [10] M. X. Goemans, and D. P. Williams, “Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming,” J. of the ACM 42, 1115 (1995).
 [11] N. G. Dickson et al, “Thermally assisted quantum annealing of a 16qubit problem,” Nature Communications 4, 1903 (2013).
 [12] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, “Defining and detecting quantum speedup,” Science 1252319 (2014).
 [13] A. Marandi, N. C. Leindecker, V. Pervak, R. L. Byer, and K. L. Vodopyanov, “Coherence properties of a broadband femtosecond midIR optical parametric oscillator operating at degeneracy,” Opt. Express 20, 7255 (2012).
 [14] J. Roslund, R. Medeiros de Araújo, S. Jiang, C. Fabre, and N. Treps, “Wavelengthmultiplexed quantum networks with ultrafast frequency combs,” Nature Photon. 8, 109 (2014).
 [15] A. Marandi, N. C. Leindecker, K. L. Vodopyanov, and R. L. Byer, “Alloptical quantum random bit generation from intrinsically binary phase of parametric oscillators,” Opt. Express 20, 19322 (2012).
 [16] E. Ising (1925) Beitrag zur Theorie des Ferromagnetismus. Z Phys 31:253258.
 [17] K. Kim, M.S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.D. Lin, L.M. Duan, and C. Monroe “Quantum simulation of frustrated Ising spins with trapped ions,” Nature 465, 590 (2010).
 [18] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner “Quantum simulation of antiferromagnetic spin chains in an optical lattice,” Nature 472, 307 (2011).
 [19] S. Utsunomiya, K. Takata, and Y. Yamamoto, “Mapping of Ising models onto injectionlocked laser systems,” Opt. Express 19, 18091 (2011).
 [20] K. Takata, S. Utsunomiya, and Y. Yamamoto, “Transient time of an Ising machine based on injectionlocked laser network,” New J. Phys. 14, 013052 (2012).
 [21] C. D. Nabors, S. T. Yang, T. Day, and R. L. Byer, “Coherence properties of a doublyresonant monolithic optical parametric oscillator,” J. Opt. Soc. Am. B 7, 815 (1990).
 [22] L. Wu, H. J. Kimble, J. L. Hall, and H. Wu, “Generation of squeezed states by parametric down conversion,” Phys. Rev. Lett. 57, 2520 (1986).
 [23] R. Graham, “Unified thermodynamics of dissipative structures and coherence in nonlinear optics,” in Coherence and Quantum Optics 851–872, (Springer, 1973).
 [24] A. Galluccio, M. Loebl, and J. Vondrk. “New algorithm for the Ising problem: Partition function for finite lattice graphs,” Phys. Rev. Lett. 84, 5924 (2000).
 [25] P. D. Drummond, K. J. McNeil, and D. F. Walls, “Nonequilibrium transitions in sub/second harmonic generation. II. Quantum theory,” Opt. Acta 28, 211 (1981).
 [26] H. Risken, The FokkerPlanck Equation: Methods of Solution and Applications, 2nd ed. (SpringerVerlag, Berlin, 1989).
 [27] H. A. Haus, and Y. Yamamoto, “Quantum noise of an injectionlocked laser oscillator,” Phys. Rev. A 29, 1261 (1984).
 [28] J. R. Dormand, and P. J. Prince, “A family of embedded RungeKutta formulae,” J. Comp. Appl. Math. 6, 19 (1980).
 [29] M. X. Goemans and D. P. Williams, “Improved approximation algorithms for maxcut and satisfiability problems using semidefinite programming,” J. of the ACM 42, 11151145 (1995).
 [30] S. Arora and S. Kale, “A combinatorial, primaldual approach to semidefinite programs,” in Proceedings of the thirtyninth annual ACM symposium on Theory of computing (2007).
 [31] P. Klein and H.I. Lu, “Efficient approximation algorithms for semidefinite programs arising from max cut and coloring,” in Proceedings of the twentyeighth annual ACM symposium on Theory of computing, 338947 (1996).
 [32] F. Alizadeh, “Interior point methods in semidefinite programming with applications to combinatorial optimization,” SIAM J. Optimiz. 5 (1995).
 [33] C. W. Rudy, A. Marandi, K. A. Ingold, S. J. Wolf, K. L. Vodopyanov, R. L. Byer, L. Yang, P. Wan, and J. Liu, “Sub50 fs pulses around 2070 nm from a synchronouslypumped, degenerate OPO,” Opt. Express 20, 27589 (2012).
 [34] K. F. Lee, J. Jiang, C. Mohr, J. Bethge, M. E. Fermann, N. Leindecker, K. L. Vodopyanov, P. G. Schunemann, and I. Hartl, “Carrier envelope offset frequency of a doubly resonant, nondegenerate, midinfrared GaAs optical parametric oscillator,” Opt. Lett. 38, 1191 (2013).
 [35] G. Grosche, O. Terra, K. Predehl, R. Holzwarth, B. Lipphardt, F. Vogt, U. Sterr, and H. Schnatz, “Optical frequency transfer via 146 km fiber link with relative accuracy,” Opt. Lett. 34, 2270 (2009).