Ising Hamiltonian Minimization: Gain-Based Computing with Manifold Reduction of Soft-Spins vs Quantum Annealing

Pursuing enhanced computing speed and power efficiency has led to exploring alternatives to traditional electronic systems in solving complex tasks. Optical Neural Networks (ONNs) promise unprecedented parallelism, potentially superior speeds, and reduced power consumption. ONNs encode neural weights as phase shifts or changes in light intensity, with activation functions instantiated via nonlinear optical materials or components, or via a strong hybridization to matter excitations [1]. They offer the potential to operate in the terahertz range, vastly surpassing the gigahertz frequencies of conventional electronic systems that can be exploited in machine learning and combinatorial optimization. The common feature of ONNs is to utilize a network of optical oscillators dynamically described by a coupled system of soft-spin models on complex-valued fields $\psi_{i}=r_{i}\exp[i\theta_{i}]$ that have amplitude $r_{i}$ (referred to as ‘soft mode’) and phase $\theta_{i}$ (discrete e.g. $\theta_{i}\in\{0,\pi\}$ or continuous ‘spin’) degrees of freedom. Each spin in the network can be associated with the quadrature of the optical complex-valued fields, thereby in the classical limit reducing the system to a model of real soft-spins given by $r_{i}\cos\theta_{i}$ which we analyze hereafter.

Optical parametric oscillator based coherent Ising machines (CIMs) [2, 3, 4, 5], lasers [6, 7, 8], spatial light modulators (SLMs) [9], lattices of polariton [10, 11] and photon condensates [12], Microsoft’s analog iterative machine [13], and Toshiba’s simulated bifurcation machine [14] can all minimize the classical hard-spin Ising Hamiltonians $H_{\rm I}=-\sum_{i,j}J_{ij}s_{i}s_{j}$ with $s_{i}=\pm 1$ for a coupling matrix $\mathbf{J}$, and other spin Hamiltonians (e.g. XY Hamiltonians $H_{\rm XY}=-\sum_{i,j}J_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}$ with $\mathbf{s}_{i}=(\cos\varphi_{i},\sin\varphi_{i})$) using soft-spin bifurcation dynamics via the Aharonov-Hopf bifurcation [15]. This principle of operation leads to an exciting new paradigm known as “gain-based computing”. The concept behind gain-based computing is that computational problems can be encoded in the gain and loss rates of driven-dissipative systems, which are then driven through a symmetry-breaking transition (bifurcation), selecting a mode that minimizes losses. Such soft-spin models exploit enhanced dimensionality, marked by small energy barriers during amplitude bifurcation, but also complicate the energy landscape with numerous local minima. In parallel to these methods, quantum annealing is another approach to minimize the hard-spin Ising Hamiltonian. Despite numerous studies contrasting classical and quantum methods, the limitations of currently available hardware and the limitations of simulating quantum systems classically have led to contrasting conclusions as to whether a quantum advantage can potentially be realized using quantum annealing [16, 17, 18, 19, 20, 21, 22, 23] and in particular, how quantum annealers such as D-Wave perform in comparison with CIM [24]. In the latter, the connectivity of the coupling matrix was assumed to be a key factor in performance differences between these machines [24].

An all-optical scalable ONN was recently proposed for cyclic graphs; SLMs are used to discretize the optical field where each pixel defines a different pulse amplitude [25]. The SLM with $M_{x}\times M_{y}$ pixels is set up with a transmission function $\tilde{J}_{k}$ which multiplies the Fourier transform of the amplitudes at each round trip. The SLM couples the fields with coupling matrix $J_{ij}\equiv\tilde{J}_{j-i+1}$, which corresponds to a circulant graph. An alternative setup allows for any general coupling matrix. However, there is an $N=M_{x}$ limit to the number of pulses. Circulant graphs such as Möbius ladders therefore lend themselves well to optical solvers, where $N=M_{x}\times M_{y}\sim 10^{6}$ spins can be defined.

The couplings are often geometrical in polariton condensates, photon condensates, and laser cavities (e.g. the sign and amplitude of the coupling strength depend on the distance between condensates and outflow wavenumber [26]). The condensates arranged in a circle interact with the nearest neighbors, but the interactions beyond this decay exponentially. Previously, various ways of establishing long-range interactions in polariton-based XY-Ising machines were discussed. An easier way to achieve the couplings between remote sites is to use digital micro-mirror devices (DMDs) to direct the light across the ring. DMDs were shown to perform complex (amplitude and phase) modulation. By splitting the complex field into real and imaginary parts and using the time modulation scheme of the DMD, a complex signal could be synthesized [27]. Reflecting the entire ring of condensates on itself with a radial displacement implements a 3-regular cyclic graph. Cyclic graphs are known to be computationally intractable for classical computers for sampling probability distributions of quantum walks [28].

Using ONNs for optimization has shown promise, yet key questions remain: ‘what are suitable benchmarks for optical machines, how to guide annealing to aid optimization, what are the ONN energy landscape dynamics during annealing to ensure the optimal state is achieved, and what are the distinguishing features between quantum and classical annealing?’. Answers often rely on the coupling matrix $\mathbf{J}$. An instructive problem encoded in $\mathbf{J}$ should be technologically feasible, have controllable couplings, possess non-trivial structures resistant to simple local perturbations, and be mathematically tractable. Moreover, it is better to have deterministic rather than random couplings to avoid issues of statistical convergence [29].

Here, we analyse and contrast gain-based computing for soft-spin Ising models (SSIM) with quantum annealing for circulant coupling matrices, which allow complete control of frustration, energy gaps, and the structure of critical points. Furthermore, the potential to realize them in future optical systems [25] make them more suitable to consider over previously reported benchmarks [30, 31, 32]. A highlighted challenge for SSIM annealing lies in the opposing relationship between local and global minima when mapping the Ising Hamiltonian to the energy of the soft-spin system [29]. Notably, we demonstrate that quantum annealing within the whole Hilbert space of the hard-spin system navigates this challenge. Additionally, we suggest that ‘manifold reduction’, aligning amplitudes to the mean, is needed to augment the likelihood of SSIM finding the global minima.

ONNs based on laser operation leverage quantum-inspired principles such as coherence, interference, and parallelism. They are dissipative systems that tend to minimize losses on their route to coherence. The losses can be written as an ‘energy’ (‘cost’) function to be minimized. For instance, in the classical limit, CIM’s energy landscape to be minimized is

where $x_{i}$ are quadratures of the OPOs, $p(t)$ describes the effective laser pumping power (injection minus linear losses), and $C$ corresponds to the strength of saturable nonlinearity. As $p(t)$ grows from a large negative to large positive $p(t)=p_{\infty}$ values, $E$ anneals from the dominant convex first term on the right-hand side of Eq. (1) that is minimized at $x_{i}=0$ for all $i$, to the minimum of the second term which is the scaled target Ising Hamiltonian with $x_{i}=\pm\sqrt{p_{\infty}}$. The temporal change of $p(t)$, therefore, is the annealing parameter combined with gradient descent as

The operation of CIM, therefore, relies on the gradient descent of an annealed energy landscape. All ONN soft spin optimizers exploit this central principle, while the details of the nonlinearity or the gradient dynamics can vary from platform to platform [15]. In particular, CIM dynamics is an example of the Hopfield-Tank (HT) network ${\dot{x}}_{i}=p(t)x_{i}+\sum_{i=1}^{N}J_{ij}x_{j}$ also used for Ising Hamiltonian minimization [33, 34]. Another approach uses second-order resonance to project the XY onto the Ising dynamics [35]. In the next section, we describe the principles of gain-based computing and contrast it with simulated and quantum annealing.

Gain-based computing is a computational paradigm where problems are encoded in the gain and loss rates of driven-dissipative systems, as illustrated in Fig. (1)(a). These systems undergo a symmetry-breaking transition when various physical modes are excited from the vacuum state. As these modes grow, the loss function evolves until a coherent state that minimizes losses emerges. The mode that achieves the minimum of the loss function is amplified, as shown in Fig. (1)(a). Gain-based computing leverages soft-spin models, which provide enhanced dimensionality and small energy barriers during amplitude bifurcation. Although these models create a complex energy landscape with numerous local minima, making optimization challenging, they are also rich in potential solutions.

Simulated annealing (SA), on the other hand, is a classical optimization technique; see Fig. (1)(b). SA probabilistically explores the solution space by simulating the cooling of a material to reach a state of minimum energy. It uses thermal fluctuations to escape local minima, as the system trajectory shown in blue indicates, with the probability of accepting worse solutions decreasing over time. This simulates a cooling process that gradually refines the search towards the global minimum. Implemented on classical computing systems using stochastic algorithms, SA explores the energy landscape by thermal fluctuations, with a gradual reduction in temperature controlling the balance between exploration and exploitation. The performance of simulated annealing is influenced by the cooling schedule, which determines how the temperature is reduced over time, as well as the specific parameters of the algorithm.

In contrast, quantum annealing (QA) is a quantum computation method used to find the ground state of a system’s energy; see Fig. (1)(c). QA operates by evolving the system from an initial Hamiltonian, which is usually simple and convex, to the target Hamiltonian that encodes the optimization problem. This evolution relies on the principles of quantum mechanics, specifically quantum tunneling, to explore the energy landscape. Quantum annealing utilizes quantum fluctuations to escape local minima and tunnel through energy barriers, potentially leading to faster convergence to the global minimum. This approach can be advantageous in navigating complex landscapes with high barriers between local minima. In Fig. (1)(c), the varying energy landscape is shown as the annealing from the initial convex Hamiltonian to the target Hamiltonian takes place in time. The system starts at the ground state of the initial Hamiltonian and remains in the ground state if annealing is sufficiently slow. The blue line shows the system state at each moment.

Cyclic graphs with $N$ nodes are characterized by the circulant coupling matrix $\mathbf{J}\in\mathbb{R}^{N\times N}$, constructed through cyclical permutations of an $N$-vector. These graphs inherently have vertex permutation symmetry, signifying boundary periodicity and uniform neighborhoods. The structure of a circulant matrix is contained in any row, and its eigenvalues and eigenvectors can be analytically derived using the $N$ roots of unity of a polynomial, where the row components of the matrix act as coefficients: $\lambda_{n}=\sum_{j=1}^{N}J_{1,j}\cos[\frac{2\pi n}{N}(j-1)]$ [36, 37, 38]. We consider the minimization of the Ising Hamiltonian on a particular form of cyclic graph – Möbius ladder graphs with tunable hardness. These have even $N$ such that the $i$-th vertex has two edges connecting it to vertices $i\pm 1$ with antiferromagnetic coupling with strength $J_{i,i\pm 1}=-1$ (circle couplings), and an additional antiferromagnetic coupling with vertex $i+N/2$ with strength $J_{i,i+N/2}=-J$ for $J>0$ (cross-circle couplings). We denote by $S_{0}$ the state where the spins alternate along the ring so that $s_{i}s_{i\pm 1}=-1$ for all $i$ (Fig. (2)(a)), and by $S_{1}$ the state where the spins alternate everywhere except at two positions on the opposite sides of the ring: $s_{i}s_{i+1}=-1$ for all $i\neq i_{0}$ and $s_{i_{0}}s_{i_{0}+1}=s_{i_{0}+N/2}s_{i_{0}-1+N/2}=1$ (Fig. (2)(b)). When $N/2$ is odd, $S_{0}$ is always the ground state with energy $H_{\rm I}(J)=-(J+2)N/2$. When $N/2$ is even, the $S_{0}$ configuration has energy $H_{\rm I}^{(0)}(J)=(J-2)N/2$ and $S_{1}$ has $H_{\rm I}^{(1)}(J)=4-(J+2)N/2$. Therefore, $S_{0}$ [$S_{1}$] is the global minimum (while $S_{1}$ [$S_{0}$] is the excited state) if $J<J_{\rm crit}\equiv 4/N$ [$J>J_{\rm crit}$]. The eigenvalues of the coupling matrix $\bf{J}$ for the Möbius ladder with $J_{1,j}\in\{-1,0,-J\}$ are $\lambda_{n}=-2\cos(2\pi n/N)-J(-1)^{n}$. Equating the two largest eigenvalues $2\cos(2\pi/N)+J$ and $2-J$ gives the value of $J=J_{\rm e}=1-\cos(2\pi/N)$ at which the leading eigenvectors change. When $J_{\rm e}<J<J_{\rm crit}$ the eigenvalues for $S_{0}$ are less than that for $S_{1}$, despite $S_{0}$ being the lower energy state (see Appendix A for the derivation of the spectra). This is in contrast to computationally simple problem instances, in which the ground state minimizer is located at the hypercube corner of the projected eigenvector corresponding to the largest eigenvalue [39].

SSIM in Eq. (1) has real amplitudes $x_{i}$. As the laser pumping $p(t)$ increases from negative values, the minimizers ${\bf x}^{*}$ of Eq. (1) and minima of $E$ change. We associate Ising spins with $x_{i}$ via $s_{i}=x_{i}/|x_{i}|$. We expect that the soft-spin energy state $E_{0}$ that corresponds to the hard-spin Ising state $S_{0}$, and depicted in Fig. (2)(a), is symmetric in amplitudes as all spins experience the same frustration of the cross-circle coupling, so all amplitudes have the same modulus $|x_{i}|=X.$ From Eq. (2), $X$ satisfies $X=\sqrt{p(t)+(2-J)/C}$, with the corresponding soft-spin energy $E_{0}=(J-2)N(2-J+2Cp)/4.$ This state can be realized from a vacuum state when $p(t)$ exceeds $(J-2)/C.$ The soft-spin energy state $E_{1}$ corresponding to $S_{1}$, when two side edges are frustrated, is asymmetric in amplitudes. This asymmetry is shown schematically in Figs. (2)(b) and (c), in agreement with the dynamical simulations presented in Fig. (3). This occurs because the lower energy is achieved if the amplitudes connected by the frustrated edges $|x_{i}|=X_{L}$ are lower than in the rest of the system. For $N=8$ in Fig. (2)(b), there are two types of amplitudes: 4 nodes with $\pm|X_{L}|$ and 4 with amplitudes $|x_{i}|=X_{B}$, where $X_{B}=(1-J-Cp)X_{L}+CX_{L}^{3}$ as obtained from the steady states of equation Eq. (2) governing the dynamics of $X_{L}$, while the steady-state on the evolution of $X_{L}$ gives $(p+1+J)X_{B}+X_{L}=X_{B}^{3}$. By solving the polynomial equation for $X_{L}$, we can compute $E_{1}$ across any $p$, $J$, $N$, and $C$. This allows us to discern regions in this parameter space where the global minimum aligns with either $E_{0}$ or $E_{1}$ and confirm if these states correspond with the hard-spin Ising Hamiltonian’s global minimum. Figure (2)(d) depicts distinct regions in the $J-p$ parameter space. Within the $J_{\rm e}<J<J_{\rm crit}$ interval $S_{0}$ emerges as the hard-spin Ising model’s lowest energy state. For the soft-spin model, however, only the region shown in pink corresponds to this state ($E_{0}$). Figure (2)(d) shows that for values $J_{\rm e}<J<J_{\rm crit}$, as laser power $p$ rises, the $E_{0}$ state becomes the energy minimum for the soft-spin model, aligning with the hard-spin Ising Hamiltonian’s $S_{0}$. Yet, the success probability of converging to the true ground state does not increase beyond $0.2$ as shown in Fig. (2)(e). This is a consequence of increasing amplitudes that bring the increased height of the energy barriers that prevent the system from transitioning to the state $S_{0}$ (see Fig. (2)(f)). Figure (4) depicts the basins of attraction for various $p$ and fixed $J_{\rm e}<J<J_{\rm crit}.$ The basins of attraction are defined as the sets of points randomly distributed on $[-1,1]$ from which gradient descent leads to different minima. At the threshold of large negative $p$, the basin of attraction of $E_{1}$, which is the ground state of $E$ as given by Eq. (1), dominates. As $p$ increases, the basin of attraction of the excited state $E_{0}$ increases while at small positive values of $p$, $E_{0}$ becomes the ground state. With a further increase of $p$, other states with even higher energy appear.

The space structure of soft spin models can be further understood by considering the critical points of their energy landscape for different annealing parameter values [40]. We can determine the critical (minima and the saddle) points by setting ${\partial E}/\partial x_{i}=0$ for all $i=1,\ldots,N$ and classify them using the Hessian matrix. The number of critical points grows exponentially fast with $p$; however, not in terms of energy and the distance from the state $x_{i}=0$ $\forall$ $i$ as Fig. (5) illustrates.

The $S_{1}$ state is always further away from the origin than other critical points ($S_{0}$ and saddle points). At the same time, the transition between minima $E_{0}$ and $E_{1}$ is possible only through a saddle point whose energy relative to $E_{1}$ and $E_{0}$ defines the height of the energy barriers; see Fig. (2)(f).

The aforementioned considerations suggest that amplitude heterogeneities have a severely detrimental effect on the optimization process in some regions of parameter space as they allow the soft-spin energy landscape to find and follow its ground state, which is quite different from the ground state of the hard-spin Ising Hamiltonian. This problem has been recognized before [41, 42], but in the context of the final state, so various feedback schemes were suggested to bring all amplitudes to the same value, say $\pm 1$, at the end of the simulations. This could be achieved, for instance, by changing the laser intensity individually for each spin as $\dot{p_{i}}(t)=\varepsilon(1-x_{i}^{2})$, where $\varepsilon$ is a small constant parameter. However, as our results on the simple circulant graphs illustrate, this feedback does not change the most essential part of the dynamics during the pitchfork bifurcation from the vacuum state. Moreover, this feedback becomes important only for amplitudes sufficiently close to $\pm 1$ when the barriers between states are already too high.

Instead, we suggest introducing feedback restricting the soft spin energy landscape to keep them close to the average value. This restriction can be achieved by modifying the signal intensities bringing them towards the average mass per particle defined by the squared radius of the quadrature $R({\bf x})\equiv\sum_{i=1}x_{i}^{2}/N$ as

If $\delta=0$, then no adjustment is made. If $\delta=1$, then all amplitudes are set to the same (average) value. For $0<\delta<1$, $1/\delta$ determines the proportion of the effective space for the restricted evolution.

Figure (6) shows the probability of finding the ground state of the Ising Hamiltonian using the HT networks: Eq. (2) (CIM-I), Eq. (2) with individual pumping adjustments $p\rightarrow p_{i}$ according to $\dot{p_{i}}=\varepsilon(1-x_{i}^{2})$ (CIM-II), and Eq. (2) with manifold reduction by Eq. (3) (CIM-III). For CIM-I and CIM-III, we set $p(t)=(1-p_{0})\tanh(\varepsilon\,t)+p_{0}$. CIM-III shows a significant improvement in finding the ground state compared with other models. Thus, in soft-spin models, the imperative to constrain the manifold implies that dimensional annealing should be tailored according to the energy landscape’s characteristics. Quantum annealing, on the other hand, harnesses dimensional annealing within an extended Hilbert space. By utilizing only linear dynamics at the expense of operating within this higher-dimensional phase space, it can effectively navigate energy barriers. Next, we study the quantum evolution on the Ising energy landscape of circulant coupling matrices in order to contrast its performance with soft-spin nonlinear models.

We consider the transverse field Ising model given by

where $\hat{S}^{\alpha}$ are the spin-$1/2$ Pauli matrices, $\mathds{1}$ is the $2\times 2$ identity matrix, and $\otimes$ denotes a tensor product. The first term, $\hat{H}_{0}$, is diagonal and corresponds to the operator representation of the classical Ising Hamiltonian $H_{\rm I}$; the second term is a symmetry-breaking longitudinal magnetic field; the third term is a transverse field that results in a non-diagonal Hamiltonian operator and gives rise to the quantum Ising model. We will take the annealing term to be of the form $\gamma(t)=B/\sqrt{t+t_{0}}$ for some constant $B$ [43] and set $t_{0}=0.5$. Our quantum system is made up of $N$ spin-$1/2$ subsystems each having a basis ${\lvert\downarrow\rangle,\lvert\uparrow\rangle}$. A general state $\lvert\Psi\rangle$ of the $N$-spin system can then be written as

where the $C_{\xi}$’s are complex numbers and the basis element

for $k=1,\cdots,N$. We begin with an initial state, which is the ground state of the transverse field Hamiltonian. The initial state at time $t_{i}$ can then be expressed as

where for each subsystem $\lvert\psi_{\rightarrow}\rangle=(\lvert\uparrow\rangle+\lvert\downarrow\rangle)/\sqrt{2}$. The wavefunction is then evolved according to the time-dependent Schrödinger equation (see Appendix B for details) [44]. As $t\rightarrow\infty$, $\gamma(t)\rightarrow 0$ and the contribution of the last term decays to bring about the target Hamiltonian. Provided $\gamma(t)$ is varied adiabatically, the state evolves whilst remaining in the true ground state of the system, and settles into the target Hamiltonian’s desired ground state at sufficiently long times.

To determine the probability of finding the ground state, we compute the projection of $\lvert\Psi(t)\rangle$ onto the ground state $\lvert\phi_{{}_{\rm GS}}\rangle$ of the classical Hamiltonian, $\hat{H}_{0}$, given by $P_{{}_{\rm GS}}=|\langle\phi_{{}_{\rm GS}}\,|\,\mathopen{}\Psi(t)\rangle|^{2}$. In Fig. (7), we present numerical simulations of the time evolution of the success probability for finding the ground state of an $N=8$ spin system with $J=0$ and $B=5$. For comparison, we have also included the results for simulated annealing [45] and classical annealing by evolving a master equation [43] (see Appendix D for details). In the former, we allow transition probabilities for single-spin flips only, whereas, in the latter, we allow for all spin flips to reveal the importance of spin correlations on the success probability of finding the ground states. Such collective transitions can be important when topological constraints associated with particular spin configurations can render certain single-spin transitions ineffective at escaping local energy minima.

For $J<J_{\rm crit}$, the $S_{0}$ ground state has a two-fold degeneracy. Therefore, in the absence of a symmetry breaking term, we can expect that the probability of finding one of the ground states is $P_{{}_{\rm GS}}=1/2$. Figures (7)(a) and (b) present simulations for the case where no symmetry breaking term is included so that the system contains two degenerate energy minima. Time evolution of the probabilities for finding the system in one of the two degenerate ground states show an equal probability for the system to be found in either one of these states. Moreover, the results are relatively similar regardless of which numerical method is considered. Therefore, quantum annealing, simulated annealing and classical annealing show a similar performance in tracking the ground states as indicated by the curves representing adiabatic evolution of the system.

In Figs. (7)(c) and (d), we present simulations for the case where a symmetry breaking term is included. To introduce the symmetry breaking term, we used $0.05\lvert\xi\rangle_{S_{0}}+0.05\lvert\xi\rangle_{S_{1}}$, where $\lvert\xi\rangle_{S_{0}}$ and $\lvert\xi\rangle_{S_{1}}$ correspond to the $S_{0}$ and $S_{1}$ states, respectively. In contrast to the previous case, the behavior of the different models is now markedly different. In particular, we observe that both quantum annealing and classical annealing correctly evolve with the true ground state as indicated by the curves corresponding to the adiabatic evolution. Moreover, due to the symmetry breaking, there is one unique ground state that the two methods can follow. In contrast, simulated annealing is not always successful at tracking the true ground state even for this case where $J=0$. We found a success probability of only 67%, whereas the remaining probability is associated with the system converging to what is now a metastable state. These results demonstrate the importance of the symmetry breaking terms and how they affect the time-evolution of the ground state probabilities.

In Figs. (8)(a) and (8)(b), we present results with and without symmetry breaking terms for the case with $J=0.35$ and $B=5$ that corresponds to a hard region of the parameter space for the soft-spin models. As can be seen, now the success probability for simulated annealing degrades even in the absence of symmetry breaking terms. In contrast, classical annealing and quantum annealing continue to perform well. Although the convergence of simulated annealing can be enhanced for slower annealing rates, in general, the success probabilities are lower than the other algorithms we have investigated over a range of annealing schedules (see also Ref. [43]).

To compare the quantum annealing and semi-classical soft-spin simulations, we computed the single-spin reduced density matrix $\hat{\rho}_{1,i}$ from the pure state $\lvert\Psi(t)\rangle$. In general, the single-spin density matrix will correspond to entangled states. This is illustrated by recovering the Bloch vector from $\hat{\rho}_{1,i}(t)$ (see Appendix C for details). In the inset of Fig. (8)(a), we show the evolution of the Bloch vector with time evaluated for one of the spins (other spins show similar behavior) for a simulation with $J=0.35$ and $B=5$ in the absence of a symmetry breaking term. We see that the spin is initially aligned along the equator (consistent with the form of $\lvert\rightarrow\rangle$) but shrinks towards the origin as the state evolves. The departure of the Bloch vector from the surface of the Bloch sphere is indicative of quantum entanglement while its dynamics towards the origin signals a spin state that is maximally entangled with the rest of the system. A definite state emerges only upon measurement, which then subsequently collapses the corresponding wavefunction to one specific configuration.

These results demonstrate the striking differences between the states of a fully quantum mechanical description, and a semi-classical description considered earlier. To facilitate comparison with the deterministic semi-classical simulations, we removed the ground state degeneracy in our quantum annealing simulations by introducing the symmetry-breaking term $\hat{H}_{1}$ in Eq. (4). We set $h_{i}$ to correspond to $0.05\lvert\xi\rangle_{S_{0}}+0.05\lvert\xi\rangle_{S_{1}}$. This enforces the evolution towards a specific ground state as can be seen by contrasting the success probability for the ground states presented in Figs. (8)(a) and (8)(b). The resulting Bloch vector is shown in the inset of Fig. (8)(b), and now indicates evolution that ends at the surface of the Bloch sphere, reaching either the $\lvert\uparrow\rangle$ or $\lvert\downarrow\rangle$ state which represents a final state that is not in quantum superposition. Since the Bloch vector does not remain on the surface of the Bloch sphere during the evolution, it clearly demonstrates that though individual spins converge towards a non-correlated value, their evolution bears the imprint of inter-spin correlations. Unlike the semi-classical models, our quantum annealing algorithm consistently identifies the correct ground state for a wide range of parameters in the interval $J_{\rm e}<J<J_{\rm crit}$; (see Fig. (6)) and demonstrates that correlations play a key role in facilitating the system to converge to the true ground state. However, its performance appears to degrade near $J_{\rm crit}$. In contrast, the CIM gain-based algorithm is less sensitive near $J_{\rm crit}$ and indicates an advantage of gain-based computing on such a Möbius ladder graph. The corresponding time-dependent probability of finding each spin, $i$, in the $\lvert\uparrow\rangle$ state is presented in Fig. (8)(c) for $J=0.35$ (and for $J=0.6$ in Fig. (8)(d)). As can be seen from the initial evolution of the single-spin probability amplitude, we strongly perturb the system towards state $S_{1}$ through the form of the symmetry breaking terms used. Despite this, the results emphasize the quantum annealing algorithm’s capacity to find the correct ground state during gradual $\gamma(t)$ quenches, leveraging the quantum system’s expanded phase space.

In order to perform a more systematic study of the impact of including the symmetry breaking terms on the results presented, in Fig. (9) we present results for simulations performed by using a symmetry breaking term of the form $0.05\lvert\xi\rangle_{S_{0}}+h_{1}\lvert\xi\rangle_{S_{1}}$, where $h_{1}$ is varied from $0.005$ to $0.1$. For each value of $h_{1}$, we have evaluated the time evolution of the probability amplitudes $\langle\uparrow\,|\,\mathopen{}\hat{\rho}_{1,k}|\,\uparrow\rangle$ for each spin $k$, as well as the time evolution of the magnitude of the corresponding Bloch vectors $|{\mathbf{u}}_{k}|$. The results demonstrate that as $h_{1}$ is increased, the initial evolution of the probability amplitudes is to align the spins towards the $S_{1}$ state which is caused by the increasing contribution of the symmetry breaking term. However, as the system navigates the energy landscape, quantum correlations develop as indicated by the decreasing amplitude of the Bloch vectors of the individual spins. This emerging quantum entanglement of the spins prevents the system from becoming stuck in local energy minima and subsequently allows the spins to readjust in order to track the true ground state. Subsequently, the system converges to the true ground state that is well described by a product state as the magnitude of the Bloch vectors converge to unity. We note that during the evolution, the maximal entanglement occurs at the time when the projection of some of the spins flips to the opposite direction. This time also coincides with the time where simulated annealing fails to track the correct ground state in comparison to quantum annealing as reflected in the results of Figs. (7) and (8). We, therefore, conclude that quantum correlations play a key role in allowing quantum annealing to outperform other methods in this region of the parameter space of the Möbius circulant graph.

In summary, we analyzed the optimization of Ising Hamiltonians, contrasting the classical dynamics of semi-classical soft-spin models with quantum annealing. We discussed the challenges that arise with using semi-classical models, which are due to a broadening dimensionality landscape, especially when the models’ global minimum maps to the Ising Hamiltonian’s excited state. A solution, termed ‘manifold reduction’, is presented, constraining the soft-spin amplitudes and restricting the dimensionality landscape. On the other hand, we showed that quantum annealing inherently can traverse the Ising Hamiltonian’s energy terrain, operating within an extensive Hilbert space. The findings highlight the importance of understanding the influence of dimensionality and the energy landscape overall on optimizing physical systems. Furthermore, they demonstrate how extensions of semi-classical models to include quantum effects has the potential to assist the annealing in navigating the system towards the true ground state.

J.S.C. acknowledges the PhD support from the UKRI EPSRC DTP EP/T517847/1, N.G.B acknowledges the support from the Julian Schwinger Foundation Grant No. JSF-19-02-0005, HORIZON EIC-2022-PATHFINDERCHALLENGES-01 HEISINGBERG Project 101114978, and Weizmann-UK Make Connection grant 142568.