Demonstration of chip-based coupled degenerate optical parametric oscillators for realizing a nanophotonic spin-glass

The need for solving optimization problems is prevalent in various physical applications, including neuroscience, network design, biological systems, socio-economics, and chemical reactions. Many of these are classified as non-deterministic polynomial-time hard and thus become intractable to solve as the system scales to a large number of elements. Recent research advances in photonics have sparked interest in using a network of coupled degenerate optical parametric oscillators (DOPOs) to effectively find the ground state of the Ising Hamiltonian, which can be used to solve other combinatorial optimization problems through polynomial-time mapping. Here, using the nanophotonic silicon-nitride platform, we demonstrate a spatial-multiplexed DOPO system using continuous-wave pumping. We experimentally demonstrate the generation and coupling of two microresonator-based DOPOs on a single chip. Through a reconfigurable phase link, we achieve both in-phase and out-of-phase operation, which can be deterministically achieved at a fast regeneration speed of 400 kHz with a large phase tolerance.

Simulated GVD of the silicon nitride (SiN) microresonator for the fundamental transverse electric (TE) mode for two different waveguide cross sections of 730×1050 nm and 730×1100 nm.
We characterize the interference between the DOPO using a setup similar to Fig. 2(c) in the main text. The SiN chip is pumped using two frequency nondegenerate pump lasers that are offset by ±1 free spectral range (FSR) from the degeneracy point. The pump wavelengths are 1546.1 nm and 1554.2 nm. One of the pumps is modulated using an acousto-optic modulator with 250-ns pulses at a repetition rate of 2 MHz. Supplementary  Figure 2(a) shows the measured temporal interference signal along with each individual DOPO, and we observe constructive interference between the two signals. In addition, we observe fast modulations in the interference signal and DOPO 1 . To understand the modulations, we measure the beatnote between an external continuous wave laser and two different DOPO states [ Supplementary Figure 2 (b)]. Modulations in the time trace are observed for an DOPO in state 1, and we observe sideband in the radio frequency (RF) spectrum that is 54 MHz from DC. We attribute these sidebands to nondegenerate frequency components that are generated in addition to the DOPO signal within a single microresonator cavity resonance. These frequency components are generated due to degenerate four-wave mixing (FWM) from a single pump from the fact that the pumps lie in the anomalous GVD regime.While it is possible to suppress the modulations via pump-cavity detuning, which can be seen in the relatively flat temporal profile for DOPO 2 which is in state 2, operation in the normal GVD regime allows for a stronger and more stable DOPO signal.

Supplementary Note 2 -Characterization of coupling strength
We characterize the ratio between the DOPO 2 power and the power that is coupled from DOPO 1 by measuring the power in the bus waveguide that is coupled to the microresonator corresponding to DOPO 2 when DOPO 2 is oscillating and not oscillating. The device under investigation here is the same as that from the main text. For the case where DOPO 2 is not oscillating, we use the integrated platinum heaters to detune the cavity resonance. Supplementary Figure 3 shows the measured optical spectra for the two cases. We observe a power difference of 26.4 dB, which corresponds to a ratio between the coupled field from DOPO 1 and the DOPO 2 field of 0.0479.  Figure 3: DOPO coupling ratio. Measured output spectra from output port of DOPO 2 when DOPO 2 is (blue) and is not (magenta) oscillating. DOPO 1 is always oscillating for both cases. We measure the difference between the DOPO 2 power and the power coupled from DOPO 1 to be 26.4 dB.

Supplementary Note 3 -Characterization of uncoupled DOPO
We investigate the behavior of two uncoupled DOPO's generated from the same pump. The device is shown in Supplementary Figure 4. The structure is similar to that shown in Fig.2(a) in the main text but the coupling waveguide is not included here. The microresonators have an FSR of 500 GHz and a cross section of 730×1050 nm. The two pumps are coupled in and are split 50/50 using a multimode interference (MMI) splitter and sent to DOPO 1 (bottom) and DOPO 2 (top). The pumps are offset by ±2 FSR's from the degeneracy point. The experimental setup for DOPO generation and the interference measurement is the same as that shown in Supplementary Figure 2(c). The high wavelength pump is modulated using an acousto-optic modulator with 310ns pulses at a repetition rate of 400 kHz in order to ensure that the DOPO is extinguished after each generation such that each successive DOPO initiation is uncorrelated. Supplementary Figure 5 shows the generated spectrum for DOPO 1 and DOPO 2 .
Supplementary Figure 6 shows the time domain interference measurement (green) along with the temporal measurement of each individual DOPO signal (red and blue). In contrast to the coupled DOPO measurement shown in Fig. 4 in the main text, we observe that there is no fixed phase relationship between the two generated DOPO's indicating that each of the OPO's are oscillating independently and generating uncorrelated bi-phase states.

µm
Supplementary Figure 4: Microscope image of uncoupled DOPO device. An multi-mode inteference (MMI) splitter is used to route the pumps to DOPO 1 (bottom) and DOPO 2 (top). Contrary to the device shown in Fig. 2(b) in the main text, the coupling waveguide is omitted.

Supplementary Note 4 -Theoretical modeling of coupled DOPO system
The propagation of a slowly varying light field envelope B (m) (z, t) in a one-dimensional dispersive medium exhibiting instantaneous Kerr nonlinearity is described by a generalized nonlinear Schrödinger equation which can be expressed in the frequency domain as where Ω is the frequency variable,B (m) (z, Ω) = F B (m) (z, t) with F indicating the Fourier transform, α 0 is a propagation loss coefficient, β 1 is the reciprocal of the light group velocity v g at a reference frequency ω 0 , β k (k ≥ 2) are dispersion coefficients associated with the Taylor expansion of the propagation constant β(ω) at ω 0 , and γ is the nonlinear coefficient. Note that denotes convolution andB(z, Ω) =B * (z, −Ω), with * indicating the complex conjugate. In order to describe the evolution of an optical field in a microresonator, Supplementary Equation 1 must be complemented by the cavity boundary conditioñ where δ 0 = 2pπ − ω 0 t R is the linear phase detuning experienced by the carrier field with frequency ω 0 relative to the nearest (pth) cavity resonance, L and t R = β 1 L are the cavity length and roundtrip time, respectively, θ is the input coupling coefficient, andẼ in (Ω) is the input pump field. It is convenient to introduce a change of variablẽ B (m) (z, Ω) =Ẽ (m) (z, Ω)e iβ1Ωz so that Supplementary Equations 1 and 2 are rewritten as This change of variable is equivalent to defining a shifting time reference frame, or τ = t − β 1 z such that τ and Ω form a pair of conjugate variables. Note that the mathematical form of the convolution term is preserved.
Let us now consider a specific pump field of the formẼ in (Ω) = A in [δ(Ω−Ω + )+δ(Ω−Ω − )], where Ω + = 2nπ·FSR and Ω − = −2nπ · FSR. Introducing this pump form into Supplementary Equation 4 and carrying out the usual procedures [2], we obtain the Lugiato-Lefever equation [3] for the dual-pump case Here t = mt R is a slow evolutionary time variable which keeps track of the number of roundtrips completed by the intracavity field. The pumps are symmetrically located with respect to the frequency ω 0 and each of them is an integer multiple (n) of FSR away from ω 0 . Therefore, in the absence of dispersion, the two pumps and the degenerate signal are detuned to the same extent relative to their respective nearest resonances, that is, by δ 0 . It has been shown that Supplementary Equation 5 allows modeling of the dynamics of Kerr optical parametric oscillators under certain circumstances [4], however predicts for our current experimental system that there is no stable oscillation of the degenerate signal. A major limitation of Supplementary Equation 5 lies in the fact that the pumps and the signal are constrained to the same detuning δ 0 . It completely rules out our ability to independently and individually tune the two pumps. To incorporate this additional experimental degree of freedom into our model, we are led to consider a more general case such that the pump becomes Assuming the fieldẼ (m) in Supplementary Equation 3 varies slowly over a length-scale of one roundtrip, Supplementary Equation 3 can be integrated over the cavity length L by the first-order Euler's method. We substitute the result of the integration step into Supplementary Equation 4 and obtaiñ Generally, |δΩ + | = |δΩ − | and the spectral positions of the pumps are no longer symmetric about ω 0 . The symmetry can be restored by introducing a shifted frequency Ω such that the pump term in the shifted frequency is where δΩ = δΩ + − δΩ − . Substituting the shifted frequency into Supplementary Equation 9 results iñ Let us pay a close attention to the term e iβ1Ω L . Since the pumps are no longer separated from each other by an integer number of the cavity FSR, we should be cautious with the definition of the frequency Ω . More specifically, the frequency grid with a step size of FSR is not appropriate. Instead, the frequency axis should be discretized with adjacent points spaced by As an example, we see that the qth frequency component from the central cavity mode will experience an additional frequency shift [2qπ · FSR + qδΩ/(2n)]β 1 L = 2qπ + qδΩ/(2n · FSR), since β 1 L = 1/FSR. In other words, we have In view of the continuous variable Ω , Supplementary Equation 14 can equivalently be written as As in Supplementary Equation 5, we define the slow evolutionary time t = mt R so that roundtrip-to-roundtrip variation of the field envelope can be approximated as a first-order derivative, that is,

Substituting Supplementary Equation 19 into Supplementary Equation 18
and taking the inverse Fourier transform, we obtain the following mean-field equation For most practical situations, | | FSR so that | β k+1 | |β k | and corrections to the dispersion terms can be neglected. Let us additionally define Here, α is half the fractional roundtrip loss, δ 0 is the effective phase detuning, δ 1 is the mode-dependent detuning, and Ω 0 corresponds to the pump offset frequency from the degeneracy point. By substituting these new variables into Supplementary Equation 20, we obtain the following modified Lugiato-Lefever equation Finally, we extend Supplementary Equation 25 to emulate a system of N coupled DOPO's. To this end, we introduce the complex coupling coefficient κ jk = |κ jk |e iφ jk , which describes a linear coupling from the k th to j th DOPO's with strength |κ jk | and phase φ jk . The emulation of a system of N coupled DOPO's amounts to solving the following coupled Lugiato-Lefever equations with j, k = 1, 2, 3, ..., N . Our experimental system of two DOPO's with unidirectional coupling can be simulated by setting N = 2, and κ 12 = 0 and κ 21 = κ with κ being the experimentally measured coupling coefficient. The coupling phases of 0 and π correspond to the cases of in-phase and out-of-phase couplings, respectively. Under these conditions, dropping the primes in Supplementary Equation 26 leads to the coupled equations in the main text.

Supplementary Note 5 -Characterization of DOPO with low pump power
We have characterized a single SiN microresonator to determine the oscillation threshold of the system. The SiN microresonator has a cross section of 730×1050 nm and an FSR of 500 GHz. The quality factor is 1.6 million. Supplementary Figure 7 shows the generated spectra for three different combined pump powers. We achieve oscillation with combined pump powers as low as 50 mW. The data point is included in Fig. 6 in the main text.

Supplementary Note 6 -Characterization of system scalability
In order to allow for simultaneous parametric oscillation, the resonances for DOPO 1 and DOPO 2 must be matched. To characterize the deviation of the resonances for the two microresonators in our coupled-DOPO system, we have performed measurements of the absolute wavelength of the lower wavelength pump for the 10 devices that are on one chip. Supplementary Figure 8 shows the low wavelength pump resonance for DOPO 1 and DOPO 2 (top) and the difference between the two resonance frequencies (bottom). Here, each die is 1×1 mm so the 10 devices span over 1 cm. We calculate an average difference between the two frequencies of 4.68 GHz. We attribute the variation in the resonance frequencies between the two microresonators within the same die to the resolution of electron beam lithography. Based on numerical modeling of our waveguide geometry, a 1-nm variation in the microresonator width corresponds to a resonance shift of 13 GHz. We see that all but one die has a frequency difference well within this range, confirming that we have variation within 1-nm. In addition to the variation within a single die, we also observe a decrease in resonance frequency of 120 GHz over the 10 dies. We believe that this is due to the nonuniformity of the SiN film thickness due to the deposition process. There are several ways to address this issue. First, our silicon nitride films are deposited in a two-step process, between the two steps, we can twist the wafer orientation and change the wafer position to compensate the film thickness variations [5]. Second, a post chemical mechanical planarization (CMP) process can be used to reduce the film thickness variations. Many parameters can be tuned in the CMP machine such as backside pressure, rotation speed and load pressure, down load pressure, such that different film removal rates can be achieved at different locations. Previous demonstrations have shown reduction of the root-mean-squared roughness of the top surface of SiN from 0.38 nm to 0.08 nm [6]. Third, a wafer map can be measured using measurement tools such as Filmetrics, which can provide nm or even sub-nm film thickness accuracy. With the wafer map, the devices' widths can be designed to compensate the relative height difference across the wafer. With this accuracy in film thickness the deviation in FSR will be largely due to the electron beam resolution. Previous work has shown that it is possible to fabricate structure below 1 nm in size [7], which is consistent with our characterization based on our measured resonance variations.  Next, we calculate the expected free-spectral range (FSR) mismatch based on a 1 nm difference in microresonator width. The resonance frequency satisfied the condition 2πm = n(ω) ω c L. We assume that DOPO 2 is tuned thermally to match DOPO 1 . The FSR's of DOPO 1 and DOPO 2 can be expressed as FSR 1 = c [n1(ω)+ω , respectively, where n is the effective index, L is the circumference of the resonator, and ∆n is the thermal index change. In order for the two resonances to match at ω = ω 1 , the condition n 1 (ω 1 ) = n 2 (ω 1 ) + ∆n must be satisfied. Thus, the FSR equations at ω = ω 1 simplify to . We use a finite element mode solver to model the wavelength dependent effective index of our device. For our coupled-DOPO device in the main text, the waveguide cross section is 730×1050 nm. For a 1 nm change in waveguide width, we calculate a 43.5 MHz change in FSR. This FSR difference is sufficient for our current system which has a resonance linewidth of 300 MHz. However, for higher Q-factor microresonators, the FSR variation must be significantly reduced. To achieve this, we consider a wider width microresonator with a cross section of 730×3000 nm. Remarkably, the FSR difference is significantly reduced to 2.9 MHz for a 1-nm change in width. This allows for alignment of the DOPO's within the 19 MHz linewidth for a Q = 10 7 . We have verified through numerical modeling that the system indeed oscillates with 1 mW of pump power (Supplementary Figure 9). The GVD still satisfies the conditions required for degenerate parametric oscillation [1], and we observe two distinct DOPO phases for different initiations that are offset by π. This approach provides a path towards scaling our system to a large number of coupled-DOPO's.
In addition, we characterize the thermal tuning properties of our microresonators. Here, we use two separate microresonators, Ring 1 and Ring 2 , with similar FSR's and characterize the frequency shift using integrated heaters for a resonance near 1560 nm. To better characterize the thermal frequency shift, we choose two microresonators with a relatively large resonance offset. Supplementary Figure 10 shows the transmission measurement of the two rings. We characterize the resonance frequency shift ∆f m as the electrical power to the integrated heater for Ring 1 is changed. When the heater is turned off, the difference in resonance wavelengths near 1560 nm is 0.4 nm. We measure the transmission for several electrical power values and extract a thermal frequency coefficient of 0.42 GHz/mW. For the DOPO system described in the main text, we apply 10.9 mW of electrical power to the heater. Based on this coefficient, the frequency shift ∆f m = 4.6 GHz, which agrees well with the frequency difference shown in Supplementary Figure 8.

Supplementary Note 7 -Numerical demonstration of an integrated photonic Ising machine
The aim of this section is to establish that our proposed system indeed minimizes the Ising Hamiltonian and behaves as a photonic Ising machine. We employ the coupled Lugiato-Lefever equations (Supplementary Equation 26) to tackle the Max-Cut problem of an undirected graph with 100 nodes. The graph we consider is the configuration known as Möbius ladder, the Max-Cut problem of which has been used as a bench-mark problem for evaluating the performance of physical photonic Ising machines [8]. We encode this problem by enforcing κ j,j+1 = κ j,j−1 = κ j,j+ N 2 = −|κ|, where N = 100, and |κ| = 0.048 · θ ≈ 9.1 × 10 −5 is the strength of the coupling used in our experiment. The periodic boundary condition E N +1 = E 1 and E 0 = E N is also enforced. The lack of directions in the graph is satisfied by the symmetry constraint κ jk = κ kj . By numerically solving the 100 coupled Lugiato-Lefever equations and tracking the evolution of the signal phases of all 100 DOPO's, we indeed confirm that our system attempts to minimize the Ising Hamiltonian. Typical evolutions of the signal powers and phases are shown in Supplementary Figure 11(a) and (b). In addition, Supplementary Figure 11(c) shows the evolution of the associated Ising Hamiltonian, which clearly reveals that the system finds the ground state with the associated minimum energy of -146 (see the argument below). In this specific trial, we use δ 0 = 0.0031 and δ 1 = 7.5 × 10 −16 s. By ramping the detuning at a rate of 10 −6 α per roundtrip (see Supplementary Figure 12), we achieve an annealing time T ann = 0.2 μs with a 35 % ground state success probability. We calculate a time to solution T sol = log(0.01) log(1−P ) = 2.1 μs [9,10]. This indicates that the spatial-multiplexed DOPO system potentially offers substantial acceleration as compared to the fiber-based DOPO system (T sol = 2.3 ms) [10] and the photonic recurrent Ising sampler (T sol < 1 ms) [11] for solving cubic Max-Cut problems with N = 100, and compares favorably to the state-of-theart digital and electronic systems [12].
Supplementary Figure 11: Simulated evolution of the 100 coupled DOPO's. Evolutions of the (a) signal powers, (b) signal phases, and (c) the associated Ising Hamiltonian.
The highly symmetric nature of the Möbius ladder configuration enables us to deduce the maximum cut as follows. Let us first consider a one-dimensional ring of N = 2n (n is an arbitrary natural number such that N is an even number) nodes with the nearest-neighbor coupling of equal weighting. The i th node can have either s i = +1 (spin up) or s i = −1 (spin down). This ring configuration satisfies the periodic boundary condition s N +1 = s 1 . The Max-Cut problem of this graph can be implemented by assigning J i,i+1 = J i+1,i = −1, and 0 otherwise, with the resulting Hamiltonian H = N i=1 s i s i+1 . It is easy to see that the energy is minimized when we assign spins of alternating signs along the ring. The associated minimum energy is H min = −N . Note that there is an inherent degeneracy arising from an arbitrary choice of the spin in the first randomly chosen node. We ignore this degeneracy. The one-dimensional ring can be extended to the corresponding Möbius-ladder configuration by adding an edge between each node and the node half-way across the ring, that is, J i,i+N/2 = −1, for all i = 1, 2, ..., N . Now two possible scenarios arise and we must consider them separately. (1) If N = 2n is not a multiple of 4, the same spin configuration as that of the corresponding ring minimizes the system energy with the new minimum energy of H min = −3N/2. This is because N/2 is an odd number so that s i s i+N/2 = (−1) N/2 (s i ) 2 = −1. Therefore, each new connection reduces energy by one unit. (2) If N = 2n is a multiple of 4, the spin configuration of the ring is no longer the optimal solution. This is because now N/2 is an even number such that s i s i+N/2 = +1. Moreover, this penalty in energy gets multiplied by N/2 which is the number of cross-connecting edges. We therefore search for a more optimal configuration by introducing at least one domain wall. We readily see that the configuration with a pair of domain walls that are separated by N/2 nodes is the most energy-favorable solution. These domain walls only add a fixed energy penalty of +2, independent of N , while all the conflicts are resolved for the cross-connections. Therefore, the minimum energy is H min = − 3N 2 − 2 + 2 = − 3N 2 + 4. Our example simulation case of N = 100 nodes is an element of the second scenario and we expect the minimum energy to be H min = −146, which is consistent with what our simulation predicts [see Supplementary Figure 11(c)] and with [8]. We also see that the solution is degenerate with a multiplicity of N/2, due to the underlying symmetry of the graph itself and the absence of external magnetic field. The absolute location of the domain wall pair does not matter, as long as they are on opposite ends. This degeneracy is also demonstrated in our simulation which shows no preference for any one specific configuration of all possible optimal solutions.

Supplementary Note 8 -Theoretical analysis of DOPO network
In this section, we theoretically investigate the operation of a large-scale network of Kerr DOPO's as a photonic Ising machine. We first treat a single isolated DOPO. Although we have derived a modified Lugiato-Lefever equation (Supplementary Equation 25) to model our DOPO, it is not easy to reveal the behavior of a DOPO directly from Supplementary Equation 25 in its current form. We therefore simplify our analysis by only considering the evolutions of the two pumps and the frequency-degenerate signal. We also limit ourselves to the degenerate detuning case by setting δ 1 = 0, and only retain the second-order dispersion coefficient, i.e. β k = 0 for k ≥ 3. We first normalize Supplementary Equation 25 according to the conventions presented in [13], and introduce the ansatz E(t, τ ) = E S (t) + E + (t)e −iΩ0τ + E − (t)e iΩ0τ into the normalized Lugiato-Lefever equation, where E S and E ± are slowly varying Fourier amplitudes of the signal and two pumps, respectively. After simplification, we obtain the following set of coupled ordinary differential equations where S, ∆ and d 2 are the dimensionless pump strength, detuning and GVD coefficient, respectively. In the small signal limit, we can find the steady-state pump fields by solving Supplementary Equations 27 and 28. Let us assume we have found the steady-state pump fields E + = E − = √ P e iφp , with power P and phase φ p . Next, we separate Supplementary Equation 29 into its real and imaginary parts E S = s R + is I and obtain d dt The eigenvalues of the system matrix above can be found to be with the corresponding eigenvectors The signal is amplified when the eigenvalue λ + exceeds 0, that is, when the gain exceeds the normalized loss of −1.  Having rigorously established the binary phase operation of our DOPO, we next look at a network of coupled DOPOs with arbitrary connections. We introduce the index i to denote the i th DOPO. We assume that the coupling strength is uniform throughout the network, which is consistent with our simulation of the unweighted max-cut problem presented in the previous section. We also neglect the coupling of the pumps which are assumed to be negligible. In this case, all the DOPO's have identical pumps. Then, we have At threshold, we can also assume that all DOPO's have identical magnitudes and are only allowed to vary in phase by π. It can be readily seen that the DOPO's will orient their spins such that the overall system loss is minimized. In the ideal situation, the coupling phases will be either 0 or π, in which case the loss is minimized if e iψ ik E S,k = E S,i and the equation is reduced to Note how the loss has been reduced by N |κ|, where N is the number of DOPO's that are connected to the ith DOPO. We immediately deduce two decisive factors which influence the ground state probability: (1) the rate at which the pumps are tuned into resonance, and (2) the coupling strength. In particular, the former is tied to the fact that the system must be allowed to reside in the ideal detuning range where only the ground state experiences gain, before it crosses the oscillation threshold for the next excited state. Finally, the above analysis can be generalized to non-ideal situations where the coupling phases are not exactly 0 or π. In this case, the effective loss becomes 1 − |κ| k =i cos(ψ ik ) and the system becomes less effective in distinguishing the ground state from the next excited state. In particular, when ψ ik = π/2 + mπ, where m is an arbitrary integer, the DOPO network will operate essentially in a random manner, which is consistent with our experimental result of the frustrated case.
We numerically verify the significance of the rate at which the pumps are tuned into resonance from below the oscillation threshold of the DOPO network. The Lugiato-Lefever equations (Supplementary Equation 26) are employed in this investigation and we defer the further discussion of Supplementary Equations 27-29 to a later study. We perform our simulation as follows. We initially set the pump detunings from the resonances (δ 0,initial = −1.6 · α) such that the system is below the threshold for the collective oscillation. While numerically integrating the Lugiato-Lefever equations, we ramp up the detunings at a constant rate to a final value (δ 0,final = −1.45 · α) slightly above the threshold. We investigate several different ramping speeds. For each ramping speed, we perform 100 independent trials in order to estimate the success probability of the global minimum search. The result is shown in Supplementary Figure 12, where the detuning rate (per roundtrip) is normalized to the total cavity loss α. For comparison, we have run 100 independent trials for the scenario in which the detuning is fixed at the final value (i.e. δ 0 = δ 0,final ) throughout the simulation, the result of which is marked as a horizontal black dashed line in the same figure. The ramping speed is clearly a decisive factor in the overall performance of our Ising machine. In particular, we observe from the direct comparison with the fixed detuning case that the system must be tuned at a sufficiently slow rate in order to improve the performance. Also on the same axis, we superimpose a plot of time to solution for the different detuning rates. As expected, for a problem of fixed geometry and size, the time to solution remains relatively constant regardless of the ramping speed. Slight fluctuations are likely to be statistical in nature (see, for example, that the ground state probability is slightly lower for a rate of 10 −5 than 10 −4 ) and should diminish with an increased number of samples. While in principle, it is possible to simulate at an even slower rate, the direct integration of the Lugiato-Lefever equations is computationally too wasteful to pursue such a study. We will report elsewhere a more comprehensive investigation of our Ising machine based on a more efficient truncated Fourier-mode analysis which was briefly introduced above (Supplementary Equations 27-29).  Inset shows the zoom-in of the crossing region. We extract a loss of 0.45 dB per intersection.

Supplementary Note 9 -Characterization of 4-DOPO system
We preform measurements on a 4-DOPO system as shown in Supplementary Figure 13(a). MMI splitters in series are used to route the pumps to four SiN microresonators. Unidirectional coupling is implemented between each adjacent DOPO by using a coupling waveguide between the bus waveguides. To maintain a single plane, the coupling waveguide between DOPO 4 and DOPO 1 crosses the bus waveguide for DOPO 2 , DOPO 3 , and the coupling waveguide between DOPO 1 and DOPO 2 . Supplementary Figure 13(b) shows the generated spectrum for one of the DOPO's.

Supplementary Note 10 -Characterization of crossing loss
As the number of microresonators increase in our system, in order to allow for coupling between the DOPO's in a single device layer, we require crossings between waveguides. To verify that we have low loss for the intersections between waveguides, we design devices with different number of crossings to characterize the loss (Supplementary Figure 14). We measure 0.45 dB loss per crossing. This can be further optimized by using a taper design at the crossing point.