Minor embedding with Stuart-Landau oscillator networks

We theoretically implement a strategy from quantum computation architectures to simulate Stuart-Landau oscillator dynamics in all-to-all connected networks, also referred to as complete graphs. The technique builds upon the triad structure minor embedding which expands dense graphs of interconnected elements into sparse ones which can potentially be realized in future on-chip solid state technologies with tunable edge weights. As a case study, we reveal that the minor embedding procedure allows simulating the XY model on complete graphs, thus bypassing a severe geometric constraint.


I. INTRODUCTION
Solving dense graph combinatorial optimization problems has been studied thoroughly in the field of quantum annealing, where qubits and their relative coupling strengths represent graph vertices and edge weights respectively. To solve arbitrary graph problems using preexisting qubit architectures, dense graphs can be embedded into new representative graphs using gauge field constraints [1,2], or connected ferromagnetic subgraphs using minor embedding techniques [3][4][5][6][7][8]. In the latter, a standardized minor embedding to the triad structure is commonly used, allowing any all-to-all connected (dense) graph to be mapped to a planar qubit architecture [3,9,10].
In this study, we investigate the feasibility of using minor embedding to a triad structure [3,10] to simulate the dynamics in dense networks of interacting classical oscillators instead of qubits. Our study is motivated by the recent developments in designing analogue computing strategies based on optical oscillatory networks to heuristically solve complex graph problems across various platforms such as: Exciton-polariton condensates [11][12][13][14][15][16], photon condensates [17], non-degenerate optical parametric oscillators [18][19][20], photon down-conversion oscillators [21], and coupled microlaser arrays [22,23]. Even digital computer algorithms designed to efficiently solve the dynamics of oscillatory networks have displayed impressive prowess in combinatorial optimization [24,25].
Optical platforms for unconventional analogue computation [21] have many desirable properties such as photonic parallelism, low cross-talk, ultrafast timescales (i.e., high sampling rate), and low power consumption using passive elements. However, engineering high degrees of graph connectivity-which has been achieved with good control in cold atom ensembles in optical cavities [26]presents still a major hurdle in many of the abovementioned photonic platforms. As an example, in planar * S.L.Harrison@soton.ac.uk microcavity exciton-polariton condensates the coupling strength between neighboring condensates decays exponentially with their spatial separation distance [27,28], making the coupling beyond nearest neighbors usually negligible. That, plus the need to control each intercondensate coupling strength makes it nigh on impossible to solve an all-to-all connected graph beyond a handful of condensate vertices. A scheme to achieve all-to-all coupling between polariton condensate modes was theoretically proposed using a fast time-modulated nonresonant laser [29]. We show that these connectivity problems can, in principle, be overcome using the minor embedding technique to the triad structure, which realizes a more experimentally friendly, on chip, strategy to simulate dynamics and emergent behaviours in all-to-all connected networks of nonlinear optical oscillators.
We focus on two measures, and compare them between non-embedded and corresponding embedded networks, to quantify the applicability of our technique. Firstly, calculate and characterize the emergence of coherence in disordered ferromagnetically connected Stuart-Landau oscillatory networks [30] with increasing coupling strength. Second, we investigate the efficiency of embedded oscillatory networks to approximate low energy solutions of the classical XY Hamiltonian through dynamical annealing [31,32].

II. THE STUART-LANDAU MODEL
The Stuart-Landau model describes a plethora of oscillatory systems and is formally derived from the normal form of an Andronov-Hopf bifurcation. We apply this universal model to describe the dynamics of our dissipatively coupled oscillators at the graph vertices written: Here, ψ n ∈ C denotes the nth oscillator, ω n its intrinsic frequency, and J n,m is the coupling strength between oscillator n and m corresponding to the weighted edge connecting graph vertices n and m. We can rewrite the Stuart-Landau model in polar coordinates with ψ n = ρ n e iθn arriving at: If the amplitudes of the oscillators are set as equal ρ n = ρ 0 , and with zero detuning ω n = 0, then Equation 3 describes the generalized Kuramoto model [33,34]. The fixed points of a Kuramoto network are described by the minima of the following Lyapunov potential [30], which is the same as the well known XY Hamiltonian. Therefore, the phases θ n in a Kuramoto network experience a gradient descent towards local minima of the XY Hamiltonianθ n = −∂L/∂θ n . In this respect, the phasors s n = [cos (θ n ), sin (θ n )] T play the role of two-dimensional classical spins which experience phase space flow towards these minima. This is in a similar spirit to Ising machines which undergo quantum or classical annealing into the ground state of some target Ising Hamiltonian [35]. In this sense, J n,m > 0 are said to be ferromagnetic links (in-phase coupling) and J n,m < 0 are antiferromagnetic links (anti-phase coupling). Surprisingly, even with coupling induced amplitude inhomogeneity ρ n = ρ m Stuart-Landau networks have displayed impressive results in approximating the global minima of the XY Hamiltonian in dense networks [32]. In Sec. IV B we will characterize the quality of using a minor embedding technique to simulate densely connected Stuart-Landau networks by comparing their performance in approximating low energy solutions of the XY Hamiltonian against corresponding non-embedded networks.

A. Triad Graph
In the process of creating a triad graph, we must first consider the undirected complete graph (all-to-all connected graph) of N vertices K N , with vertex set: V (K N ) and edge set: J(K N ). Each vertex V n is assigned an index n and each edge symmetrically connecting two distinct vertices V n and V m is denoted J n,m = J m,n . Through the process of minor embedding, K N is mapped to the triad graph K emb N [shown for K 5 in Figure 1(a,b)] by expanding each vertex V n ∈ V (K N ) to a chain of uniform coupled vertices of length N − 1, with intra-chain edge weights set to, Each vertex of the chain is adjacent to a single vertex of a another chain with inter-chain edge weights, J inter n,m = J n,m , ∀ J n,m ∈ J(K N ), (black edges). (6) Additional description can be found in [36]. Physically, J c > 0 gives precedence to in-phase locking between oscillators within each chain, which-in the context of the XY Hamiltonians-can be regarded as a ferromagnetic (FM) type coupling. The aim of FM coupling is to minimize the deviation between oscillators across each chain, in order to better simulate the dynamics of the complete graph K N [8]. In other words, the red vertex in the complete graph [red oscillator ψ 1 in Figure 1(a)] is represented by the average amplitude and phase of the oscillators in the red chain of the triad [see grey solid box in Figure 1(b)]. Notice that the number of black edges in the complete graph is the same as in the triad graph.
In this study, we will work close to the bifurcation threshold of the system separating the attenuated state ρ n = 0 from the oscillatory state ρ n = 0, and are thus interested only in the phases of the oscillators whose average in each colored chain is written, We will refer to this phase averaging as unembedding the triad graph [37,38]. That is, the averaged phasesθ n of the triad graph have been unembedded back to the vertices of the original complete graph.
Where the triad structure in Figure 1(b) is the hardware layout for investigating the minor embedding of complete graphs, it can be easier to picture the mapping process as the mathematical layout shown in We also consider the addition of an intra-chain edge between the first and last vertex of each chain [ Figure 1(d)], thus "looping" them, such that the looped graph has a uniform degree of connectivity. Although this is more difficult to realize experimentally in most platforms, we will compare results between unlooped and looped chains at minimising the XY Hamiltonian, and elucidate on how this more "symmetric" connectivity affects the performance of the triad graph.

B. XY Energy Extraction
As mentioned above, each vertex of the complete graph K N is associated with a complex valued number ψ n containing information on the state the nth oscillator. We define the state vector of the graph as ψ = [ψ 1 , ψ 2 , . . . ψ N ] T and its corresponding phase vector θ = [θ 1 , θ 2 , . . . θ N ] T . The energy of the graph K N is calculated from the XY Hamiltonian written, On the other hand, the energy of the triad graph can follow two methods. The first one defines the energy of the unembedded triad in the same way as Equation 8 but with the average phases across the chainsθ n , We refer to this as the unembedded energy of the triad. The other method directly uses the relative phase ∆θ inter n,m between all pairs of oscillators still "embedded" in the triad graph connected by a black inter-chain edge as depicted by the grey dashed boxes in Figure 1(a,b). We will refer to this as the embedded energy of the triad graph written, Note that when J c = 0 then minimising H emb XY is trivial since the graph forms just a set of N individual oscillator pairs connected by J n,m .

A. Coherence Properties
Here, we numerically investigate and compare the phase coherence properties (i.e., the ability to synchronize in-phase) in the oscillator network dynamics between the complete graph and the triad graph. We consider a complete FM graph where each oscillator is coupled equally to all the others with coupling strength J n,m = J > 0. The frequencies ω n are randomly chosen from a normal distribution with density g(ω) of mean ω and standard deviation σ. Naturally, we can always go into a rotating reference frame with frequencyω and therefore we can setω = 0 throughout our study without any loss of generality. We will use dimensionless units for all parameters and variables and fix σ = 1 in this section. We point out that the oscillator coupling matrix J = (J n,m ) ∈ R M ×M has always at least one positive eigenvalue for all graphs considered throughout the paper which means that the trivial ρ n = 0 solution is never stable. This means that the phases θ n (t) are well defined at all times when calculating the dynamics of Equation 1.
To understand the emergent coherence properties of the network it is useful to define a phase order parameter, commonly used in analysis of such networks [30], to capture the degree of phase coherence between the oscillators. For the complete graph it is written, If all the oscillators have the same phase θ n = θ m then r complete = 1, and in the limit of infinitely many uniformrandomly distributed phases on the interval [0, 2π) one has r complete = 0. We have checked that the special case of r complete = 0 where half of the oscillators are θ n = 0 and other half θ n = π (i.e., anti-phase ordering) does not appear in our results. We numerically integrate Equation and calculate the average coherence r complete at the final time t = T over 160 random realizations of ω n and initial conditions (i.e., Monte Carlo sampling). We then repeat our calculation over a range of coupling strengths J [see Figure 2(a)] observing a gradual transition from an incoherent state to a coherent state with increasing coupling strength. This is reminiscent of the coherence bifurcation in Kuramoto networks [30,34]. There, all phase oscillators are in an incoherent state r complete = 0 below some critical coupling strength J crit defined in the limit N → ∞. As J is increased through and above J crit , the system reaches a partially synchronized state r complete > 0, where oscillators at the centre of g(0) are synchronized while those at the tails of the distribution remain in an incoherent state such that the system is split in two dynamical groups [39]. As J is increased further, more oscillators join the synchronized group until the entire system becomes coherent [33,40], as we observe with the Stuart-Landau model. Finite size effects can also be clearly observed in Figure 2(a) when changing the number of oscillators N . Notably, larger coupling strengths are needed to synchronize smaller networks. We also point out that the presence of finite coherence values as J → 0 arises from these finite size effects which gradually decrease as N becomes larger.
We now move onto the triad graph. Similar to Equation 11, we can define coherence order parameters for the unembedded (averaged) phases of the triad graph, and for the average phase coherence within each chain of the triad, For more details, please see Figure 8 in section A. Just like for the complete graph, we average the coherences of the triad over 160 random samples of the dynamics denoted by . .
Considering two different orders of the scaling parameter J c = {1, 10}, we investigate how the triad graph's coherence properties relate to the complete graph (going to lower and larger orders of J c did not qualitatively change the findings). Note that the frequencies ω n are randomly drawn from g(ω) for all oscillators in the triad which represents experimental reality (i.e., we do not associate a single random frequency across each chain). As expected, when J c is small the inter-and intra-chain coherences in Figure 2(b,c) show poor coherence with weak dependence on J. When J c is large we instead observe a fast interchain coherence transition in Figure 2(b). This indicates that the triad graph managed to represent the embedded complete graph dynamics within J/J c < 0.1 which gives a figure of merit for the design requirements of possible triad graph platforms. Interestingly, in Figure 2(b) for large J the coherence is smaller in large graphs. This can be attributed to the fact that the chains within the triad graph need themselves to be coherent, r intra ≈ 1, in order to represent the embedded complete graph. But longer chains struggle more to settle on a phase and achieve good coherence as can be seen in Figure 2(c), such that the larger graphs in general require a stronger FM coupling strength in order to build coherence.
As both J and J c increase, both inter and intra-chain coherences converge to unity and all oscillators synchronize. In Figure 2(d-i) we show example simulations of varying coupling strength where the vertex colors represent the steady state oscillator phases. The effect of looping the chains like shown in Figure 1(d) is found to have little effect on the coherence properties of the system (see Figure 9).
To better understand the quality of the triad graph as a simulator for the phase dynamics of the complete graph we calculate the root mean square error (RMSE) between the phases in the complete graph θ n and their corresponding representations in the triad graphθ n as a function of J c . For brevity, we use the same set of randomly sampled detunings ω n with strength σ, and same random initial conditions ψ n (t = 0) across all calculations in order to quantify the similarity between the complete and the triad graph. We define a gauge-invariant

RMSE measure as
where θ nm = θ n − θ m for the complete graph phases and θ nm =θ n −θ m for the representative triad graph phases. The RMSE is calculated only for the last 200 timesteps of each simulation in order to skip any transient effects that might be strongly dissimilar at early times between the complete and the triad graph.
In Figure 3(a-d), we plot the RMSE with N = 5, 10, 15 and 20 respectively, which show that the RMSE converges towards a minimum with increasing J c , where we have marked with black circles the values of J c where the gradient of ln(η) surpasses 0.05.
The corresponding J c values where ∇ ln(η) = 0.05 are collated in Figure 3(e), showing that the optimum value of J c increases with the number of oscillators and, in-terestingly, appears to be independent of the strength of detuning σ. We apply a linear fit to this data, shown by the black dashed line, with J c = 0.34N + 2.39 which gives an estimation of what value of intra-chain coupling strength in the triad J c is needed in order to simulate the dynamics of the complete graph.

B. XY Energy Minimization
In this section we investigate the feasibility in using the minor embedding technique on complete graphs of Stuart-Landau oscillators to optimize the XY Hamiltonian. We compare the unembedded and embedded complete graph XY energies extracted from the triad structure [Equation 9,10] to the ground state of the XY Hamiltonian, as calculated using the basin hopping optimization method [41]. Here we measure the performance of the Stuart-Landau system through the normalized difference in the extracted XY energy from solving the dynamics of the triad graph of Stuart-Landau oscillators E SL compared to the complete graph XY energy solved using the basin hopping method E BH , The basin hopping algorithm is time consuming but a highly accurate optimization algorithm that serves as a good reference for the energies found by the solving the dynamics of triad embedded Stuart-Landau networks. Previous studies on complete graphs (i.e., no embedding techniques) have shown very good performance from numerically solving the dynamics of Stuart-Landau networks as compared to the basin hopping algorithm [32]. For completeness, we provide evidence in section B that the Stuart-Landau network indeed outperforms commercial optimizers, and also direct integration of the Kuramoto model, in finding low energy solutions of the XY Hamiltonian. Throughout this section, the performance is obtained by averaging over 160 unique complete graphs with weights randomly selected from a uniform distribution J n,m ∈ [−1, 1] (there is no qualitative difference in using a normal distribution of same variance). To find the optimum embedding parameter J c , we first consider a graph size of N = 10 in Figure 4. In Figure 4(a) we show the system performance for σ = 0 and scanning across the embedding parameter J c . The error of the extracted XY energies reduces as J c is increased, where the unembedded XY energy converges to zero faster than the embedded energy, with a reduction in error for both methods when the triad chains are looped. This result is to be expected, as larger values of J c reduce the distribution of oscillator amplitude and phase across each triad chain, achieving a better representation of the complete graph in corroboration with Figure 2(b). This is seen in Figure 4 Results in (a,c) are averaged over 160 random graphs and initial conditions and, in the case of (c), for random oscillator energies ωn.
10 moving paths, with each path representing a different triad chain.
Interestingly, looping the triad chains in the considered case of N = 10 oscillators is found to reduce the phase variation within each chain (see Figure 9), achieving lower error (see Figure 4). However, for larger N the error starts increasing (see Figure 5 in the Appendix) which we discuss in more detail at a later stage. We point out that the negative error for the embedded XY energy stems from the fact that min[H emb XY ] < min[H XY ] when J c → 0. In this limit, the triad graph breaks into N pair-coupled oscillators (only black edges remain) and the sum of their XY energies is trivially minimized by simply setting the relative phase in each pair to ∆θ n,m = 0, π for positive and negative couplings, respectively.
In Figure 4(c) we show the performance of Stuart-Landau networks in dynamically finding a good XY energy when including random oscillator frequencies ω n , normally distributed with standard deviation σ. Here, we focus on unlooped triad graphs with J c = 10 and J c = 20. As an additional reference, we also plot the performance of the Stuart-Landau oscillators in the original complete K 10 graph (black line). As expected, the fre- quencies ω n lead to desynchronization and fast increase in error with increasing σ, and the network steady state is lost with no well-defined phase relation that can enter into the XY Hamiltonian. This increase in error in the Stuart-Landau networks occurs for σ 1 10 which is in qualitative agreement with coherence transition in Figure 2(b).
Lastly, we show the performance of the triad graph as a function of both the problem (i.e., network) size N and J c in Figure 5. As before, the error converges to zero with increasing J c with better performance for the unembedded XY energy. As expected, larger networks struggle to find the correct (global) XY energy minimum and instead stabilize into the growing number of local minima. Interstingly, in Figure 5(b), when the chains are looped the unembedded XY energy error plateaus such that beyond J c = 10, the embedded energy does not noticeably change and the benefit of looping the chains is only apparent for N ≤ 10. This puzzling behaviour implies that looping the chains in the triad graph has generated a new family of stable attractors which do not correlate with the minima of the XY Hamiltonian. Such attractors could belong to twisted states in the chains which appear in sparse networks [42], which we explore in the following section.

C. Twisted States
Twisted states in oscillatory networks are characterized by a winding number, , denoting integer multiples of 2π phase winding about a network of FM coupled oscillators [42]. Typically, these states appear in looped networks (where the first and last oscillators in a chain are connected) with sparse connectivity up to the first few nearest neighbors. Although an interesting nontriv-ial state, this family of attractors is however detrimental to the performance of the looped triad graph to find good XY energy minima. For the triad to best represent the complete graph, the FM chains in the triad should have all sites following same phase dynamics (or as similar as possible) in order to represent the dynamics of a single oscillator in the complete graph. Appearance of a vortex phase gradient in the looped triad chains is clearly detrimental to this requirement.
Networks of N identical Kuramoto oscillators have recently been studied with unit FM coupling to µ(N − 1) nearest neighbors [42], where µ > 0.75 guarantees that all oscillators will converge to a synchronized in-phase configuration with = 0 . Below this threshold however, the loop of oscillators exhibit a constant relative phase difference between nearest neighbors of 2π /N , where = 0. These solutions are known as twisted states.
To explore the occurrence of twisted states in the Stuart-Landau model, we create networks of N identical oscillators (σ = 0) in loops (just like a single looped chain of a triad graph) with unitary FM coupling and plot the distribution of different winding numbers that appear in the system over 1000 unique realizations as a function of nearest neighbor connectivity and varying N in Figure 6(a-d). Remarkably, we find that twisted states appear frequently in our simulations. For larger N , the spread in is greater, but reduces quickly as the number of nearest neighbor connections increases. When N = 20 for example, we observe | | = 0, 1, 2 for first nearest neighbor coupling, but this range decreases to | | = 0, 1 when second nearest neighbor coupling is introduced, as depicted in Figure 6(e-i) where the oscillators (circles) have unit FM coupling following the black graph edges and their steady-state phases are shown by their color. In this analysis, | | = 2 occurred in 0.1% of the cases for N = 20 with first nearest neighbor connectivity, so it is very unlikely to observe a phase winding greater than 2π in a smaller chain. This is the likely reason that looped-chain triad graphs do not converge to zero energy error in Figure 5(b), but instead plateau to a larger error that grows in value with increasing N . As non-zero winding numbers in Stuart-Landau chains grows in probability with N , we emphasize that it is important to sample the system over many realizations starting with different random initial conditions in order to increase the chances of finding an XY energy closest to the XY ground state.

D. Dynamic Pumping
The presence of random couplings induces an amplitude inhomogeneity in the Stuart-Landau network ρ n = ρ m , meaning its fixed points will, in general, not coincide exactly to the minima of the XY model. As a possible improvement to achieve amplitude homogeneity [32,43], we investigate the addition of a dynamic pumping mechanism that feeds back the amplitude of each oscillator at each step in numerical integration in order to adjust the gain of each node respectively. To explore this feedback mechanism, we adjust Equation 1 slightly to follow: Here, P n (t) is the pumping rate of oscillator ψ n where initially all oscillators are injected equally with P n (0) = 0 for all n. In terms of exciton-polariton condensates [11,16], the pumps represent an arrangement of spatially modulated nonresonant excitation beams. After the initial pumping in the first time step of integration, we apply a feedback mechanism to balance the occupation of each oscillator, following: where controls the rate of change of P n , ρ t is the target amplitude of the oscillators. We numerically integrate Equations 16 and 17 for 100 unique randomly-connected complete graphs and corresponding triad graphs with varying N , J c = 20, and for = 0 and = 0.04. The resulting performance, and amplitude dynamics for a single instance of the network is shown in Figure 7. We do indeed see a considerable drop in error for smaller networks when the feedback is present. However, for larger networks the error is dominated by the approximative nature of the embedding procedure and the role of feedback becomes less important. It could be possible to design a more complex feedback procedure which not only eliminates amplitude inhomogeneity but also helps balance the phases in each chain of the triad graph, but this is beyond the scope of the current study. Note, due to the large intra-chain couplings J c in the triad graph, the oscillator amplitudes reach saturation an order of magnitude sooner than for the equivalent complete graph.

V. CONCLUSIONS
We have demonstrated that the dynamics of an random all-to-all coupled Stuart-Landau oscillator network can be approximated using a minor embedding technique, regularly applied in the design of quantum computing platforms. Here, a complete (dense) graph is embedded into a sparse triad graph defined by a single embedding parameter J c . We show that the steady-state phases in the embedded Stuart-Landau oscillator network correspond to good approximation of the optimal solutions in the corresponding graph XY Hamiltonian, achieving good performance by simply adjusting its embedding parameter. The results are compared against the standard complete graph Stuart-Landau oscillator dynamics, and the basin hopping method.
The convergence of the minor embedded graph dynamics to that of the complete graph offers up the triad structure as a potential testbed for mapping out dense graph problems to continuous-phase coupled oscillator systems where fully-controllable all-to-all couplings are not practicable (such as polariton condensates, photonic condensates, and coupled laser arrays). This opens perspectives on designing analogue computing hardware aimed at heuristically solving dense graphs problems (such as optimising the XY Hamiltonian) across a wide range of platforms in a similar spirit to quantum computing platforms. Optical systems might hold a particularly strong promise in this regards since their GHz operation rates can provide rapid sampling of the objective function energy landscape that can serve as good trial points for more sophisticated optimizers.
Appendix A: Coherence

Definitions
In the main text we extract the coherence across all oscillators in the complete graph r complete , the coherence across the average phases of the chains r inter and the average coherence across each chain r intra . All three parameters are defined alongside a schematic in Figure 8. Here, θ n is the phase of oscillator n in the complete graph, θ intra n n is the phase of the n th oscillator in triad chain n andθ n is the average phase over the N − 1 oscillators in chain n.

Looping Triad Chains
In addition to the unlooped triad representation of the uniform ferromagnetic (FM) complete graph in the main text, we also consider the effect of looping the triad chains  on the coherence of the system under the same conditions as the unlooped case [9]. When the minor embedded chains are looped, all oscillators are symmetrically coupled and thus the triad graph coherence is higher for the same coupling strength compared to the unlooped equivalent, as the edge effects of the unlooped chain coherence are eliminated. For J c = 1, r inter and r intra remain at a low coherence as the encoded complete graph weights are not dominant in the system. For J c = 10, there is a coherence build up with increasing J, indicating that the triad coherence dynamics represent the complete graph for large J c . Here we will benchmark the performance of dense Stuart-Landau (SL) networks against commercially available global optimizers. We have applied the Basin Hopping (BH) algorithm from the Python SciPy library as a benchmark in the main manuscript. However, we point out that we are benchmarking dense graphs meaning that a N = 50 vertex graph has N (N − 1) = 2450 weighted edges to be optimized. The reason we limit ourselves to N = 50 vertices in this study is because running the BH algorithm with 1000 iterations (for good convergence) for larger graphs on a single 2.6 GHz Intel Sandybridge processor (components of the IRIDIS 4 supercomputer at University of Southampton) exceeded 100 hours.
As a result, to investigate graphs of N > 50, instead of using the BH algorithm we have decided to use the speedier and commercially available global optimizer GlobalSearch [44] (GS) from the Global Optimization Toolbox [45] of Matlab TM . This algorithm uses a scatter search method to generate feasible trial points which are then evaluated using a chosen optimizer method to find multiple local minima. The GS algorithm then iteratively analyses points that converge using a score function which updates on-the-fly and rejects those points that are unlikely to improve the best minimum found so far. We have decided to use two well known optimization methods for the GS algorithm. (1) The trust region method since we can easily compute both the gradient and the Hessian matrix of the XY Hamiltonian (i.e., the objective function) making it quite fast and accurate and thus an appropriate choice. (2) The sequential quadratic programming (SQP) gradient descent method since it also benefits from knowing the gradient of the objective function. Our search region is bounded on θ n ∈ [−4, 4] which is taken larger than the periodic range [−π, π] in order to more efficiently find minima that might be close to values around θ n = ±π. All other options of the GS algorithm were set to default as they did not considerably improve the efficiency of the optimizer.
Just like in the main manuscript, we quantify the performance between the SL model and the GS algorithm using the ratio of their energy difference, The factor 1/2 is added so that if E SL = −E GS then the difference is exactly unity. Our results are presented in 10(a) where we have averaged over 1000 random dense graphs with weights sampled from the interval J n,m ∈ [−1, 1] going up to N = 200 spins. Amazingly, even after supplying knowledge of the gradient and the Hessian of the objective function (H XY in main manuscript) the SL model starts outperforming the GS algorithm around N ≈ 30 spins. We also show in 10(b) the time taken to iterate the SL model using an explicit Runge-Kutta (4,5) method until it converged to a steady state (red line), and the time it took the trust region (blue line) and gradient descent (cyan line) methods in the GS algorithm to provide a solution. The results evidence the considerably better efficiency in iterating the SL model than applying the GS algorithm in finding a low energy value to the XY Hamiltonian.

Stuart-Landau versus Kuramoto networks
We additionally investigate the performance of the SL network at minimising the XY Hamiltonian in comparison to the fixed point solutions of the Kuramoto (KM) model. As mentioned in the manuscript, relating the steady states of the SL model with the minima of the XY Hamiltonian is a heuristic approach since, in general, oscillatory systems like lasers and driven-dissipative condensates have freely evolving amplitudes and therefore it demands investigation on how well such heuristic systems will perform in approximating the XY ground state energy. To backup this point, we show in 10(c) the average energy difference between the KM and the SL model in minimising the XY Hamiltonian over 1000 random dense graphs. Amazingly, around N ≈ 30 spins (oscillators) the SL model starts outperforming the KM model. This means that the SL oscillators are much more efficient in exploring their state space during their transient growth phase from the initial "vacuum state" |ψ n (t = 0)| 0 (in the context of quantum annealing) rather than the KM oscillators which are always locked onto the unit circle. Therefore, SL oscillators can clearly outperform KM oscillators at minimising the XY Hamiltonian of densely connected oscillator networks.