# Quantized Vortices Mediated Annealing

Airat Kamaletdinov1 and Natalia G. Berloff1,2<EMAIL_ADDRESS>1Skolkovo Institute of Science and Technology, Bolshoy Boulevard 30, bld.1,
Moscow, 121205, Russian Federation  2Department of Applied Mathematics and
Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, United
Kingdom

###### Abstract

The complex matter-field lattice model is a ubiquitous and universal physics
model that directly links to many universal spin models. However, finding the
ground state of such a model for the most general interactions between the
lattice sites is an NP-hard problem. The knowledge of the ground state is
crucially essential for the range of applications from quantum annealing to
predicting new materials and understanding the behaviour of spin glasses. Here
we show how letting the quantised quasi-3D vortex loops pass through the
system during the dissipative quench dynamics helps it find its ground state.
Such quantized vortices mediated annealing (QVMA) is achieved by introducing
and dynamically changing selected couplings between the lattice sites that
either affect the local dimensionality of the medium or change the local
mobility of the defects. Analysing the performance of QVMA on strongly glassy
systems, we demonstrate that QVMA represents a highly efficient mode of
physics enhanced computation.

Introduction. Many different physical implementations and models for quantum
and classical analog computing have been proposed in recent years. The
physical implementations include superconducting qubits, neutral atoms in
optical lattices, trapped ions, nuclear magnetic resonances, quantum dots,
photons, nuclear spins, optical computing systems Georgescu  _et al._ (2014).
Each physical implementation assumes a particular model of computation ranging
from the quantum circuit model characteristic of a universal quantum
computation Nielsen and Chuang (2002), adiabatic quantum annealing Farhi  _et
al._ (2001), topological computing based on the braiding of a certain type of
quasiparticles called anyons Wilczek (1982), quantum random walk Aharonov  _et
al._ (1993), quantum Turing machines Bernstein and Vazirani (1997),
measurement based quantum computing, permutational computing Jordan  _et al._
(2010), the one clean qubit model Moussa  _et al._ (2010), gain-dissipative
optical systems Marandi  _et al._ (2014); Nixon  _et al._ (2013); Berloff  _et
al._ (2017); Cai  _et al._ (2019). The variability of models of quantum
computation provides alternative frameworks for both discovering new quantum
algorithms and devising new physical implementations of quantum computers or
simulators. In particular, the adiabatic model seems particularly promising
for implementation in superconducting systems Kaminsky  _et al._ (2004), while
loss minimisation of gain-dissipative systems is suitable for optical
implementation Marandi  _et al._ (2014); Nixon  _et al._ (2013); Berloff  _et
al._ (2017). Simulators can be classified into digital when they use discrete
quantum-gate operations and analog when they implement a spin Hamiltonian in
an analog fashion Buluta and Nori (2009). The former is hugely demanding
technology to design and build. It requires single quantum particles to be
localized and controlled for long enough to apply hundreds or thousands of
precise manipulations to them and then measured individually with high
fidelity. The latter type of simulators is considered the most promising type
of simulators achievable on the shorter time scale using a wide variety of
experimental platforms. A particular type of computationally intractable
problem that such simulator can address consists of finding a global minimum
of universal classical spin models such as Ising, XY, Potts, etc., with the
most general coupling matrix. The universality means that any minimization
problem of interest can be mapped into any of these universal classical spin
models with at most polynomial overhead on the number of variables De las
Cuevas and Cubitt (2016). For instance, the minimization of spin Hamiltonians
is performed in quantum annealing machines using adiabatic mode of quantum
computation (e.g. D-wave) and in optical computing using the gain-dissipative
principle of computation (e.g. Coherent Ising Machines). Such physics enhanced
computation is an emergent field of research that promises to facilitate the
developments in machine learning, neural networks and large-scale
optimization.

In this Letter, we propose and justify a novel model of physics enhanced
computation based on quantized vortices that aims at finding the ground state
of the Hamiltonian of the universal complex matter field lattice model
Svistunov  _et al._ (2015)

$H=-\frac{1}{2}\sum_{i,j}J_{ij}(\psi_{i}^{*}\psi_{j}+c.c)+\frac{g}{2}\sum_{i}(1-|\psi_{i}|^{2})^{2},$
(1)

where $\psi_{i}$ represents the wavefunction associated with the $i-$th
lattice site, $g>0$ is the strength of the nonlinear repulsive self-
interactions within the $i-$th site, and $J_{ij}$ is the strength of the
coupling between the lattice sites $i$ and $j$. A large variety of the
classical and quantum complex-field systems on a lattice from ultra-cold BECs
in an optical lattice to superconducting qubits, planar magnets, 2D solids,
Josephson-junction arrays, photon and polariton systems, superfluid and
superconductor films, etc. realise this Hamiltonian Domb (2000). The
correspondence of this model with the XY model can be seen by observing that
the second term on the right-hand side represents the penalty for the
deviation of the occupation $|\psi_{i}|^{2}$ from one. For sufficiently strong
nonlinear self-interactions, therefore, the ground state (GS) configuration
favours $\psi_{i}=\exp(i\theta_{i})$ and reduces to the XY model $H=H_{XY}.$
The XY Hamiltonian for $N$ spins denoted as a vector on a unit circle ${\bf
s}_{i}=(\cos\phi_{i},\sin\phi_{i})$ takes form
$H_{XY}=-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}J_{ij}{\bf s}_{i}\cdot{\bf
s_{j}}=-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}J_{ij}\cos(\phi_{i}-\phi_{j}).$

The time evolution of the system obeys the Hamiltonian dynamics
$i\dot{\psi_{i}}=\partial H/\partial\psi_{i}^{*}$ and for Eq. (1) becomes

$i\frac{d\psi_{i}}{dt}=-\sum_{j,j\neq
i}J_{ij}\psi_{j}+g(1-|\psi_{i}|^{2})\psi_{i}.$ (2)

If the particular case of $J_{ij}=J>0$ representing the nearest neighbor (NN)
interaction on a square grid, Eq. (2) reduces to the discrete Gross-Pitaevskii
model Kevrekidis  _et al._ (2001) where the Laplacian operator is replaced
with its discretized version. Low-density Bose–Einstein condensates (BECs)
embedded into optical lattices also obey a lattice model with $J_{ij}=J$ (NN
interactions of the Josephson type) when the tight-binding approximation is
employed Brazhnyi and Konotop (2004). The sign and magnitude of the
interactions can be engineered to take a range of values Struck  _et al._
(2011) and beyond NN Davis  _et al._ (2019). From the point of global
minimisation, we are interested in the spin-glass couplings that include long-
range interactions (or even all-to-all) so we consider the sign and magnitude
varying $J_{ij}$. The systems governed by anisotropic and potentially
frustrating couplings are ubiquitous from magnetic nanoparticle ensembles to
molecular clusters bonded by electric dipolar and quadrupolar interactions
Gallina and Pastor (2020). Finding GS of these systems is a highly non-trivial
task and sophisticated algorithms borrowed from computer science are often
employed.

In this Letter, we address the question of bringing the system to the ground
state by experimentally accessible means and, therefore, to the GS of the
lattice model having two crucial applications in mind.

Firstly, the quantum annealing protocol starts with preparing the system in a
GS of an ”easy” (often referred to as the ”driver”) Hamiltonian Kadowaki and
Nishimori (1998); Farhi  _et al._ (2001). The ground state of the driver
Hamiltonian is assumed to be unique and easily prepared. The drivers are
routinely taken as the ground state of the system with ferromagnetic couplings
$J_{ij}=J>0$ Hen and Sarandy (2016). However, the connectivity of the driver
Hamiltonian is not arbitrary. The embedding of practical optimization problems
on experimental quantum annealers requires the construction of driver
Hamiltonians that commute with the constraints of the problem Hen and
Spedalieri (2016). The resulting connectivity (adjacency) matrix is likely to
contain lower-dimensional subsystems (1D and 2D), preventing the GS
achievement of the driver Hamiltonian. We show that the nonlinearity,
dissipation and quasi-3D nonplanar adiabatic connections can resolve this
problem while driving the system to the true GS.

Secondly and more importantly, we show that the analysis of local mobility of
quantised vortices in the physical system can be used to develop an annealing
scheme (that does not have to be adiabatic) to drive the spin glass to a lower
energy state and even to GS with dramatically improved success probabilities
in comparison with unassisted quenched dynamics or with annealing from the
ferromagnetic state.

Classical lattice model meets Khalatnikov’s two-fluid model. Equation (2)
describes the time dynamics of a conservative Hamiltonian system: a
condensate, however, has no preferred reference frame. In the actual
experiments, the reference frame is set by noncondensed or reservoir particles
that have to enter the theory. For instance, Khalatnikov’s theory of mutual
friction between the superfluid and normal fluid at finite temperatures
involves three independent kinetic coefficients to describe the energy
transfer between the components Khalatnikov (1965). According to this theory,
the integrated superfluid momentum equation is
$\phi_{t}+\frac{1}{2}u^{2}+\mu=-\zeta_{3}\nabla\cdot(\rho{\bf w}),$ where
$\phi$ is the phase of the field, ${\bf u}=\nabla\phi$ is the superfluid
velocity, $\rho$ is the superfluid density, $\mu$ is the chemical potential,
${\bf w}$ is the velocity of the normal fluid relative to the superfluid and
$\zeta_{3}>0$ is the third kinetic coefficient of Khalatnikov’s theory. Using
the equation of mass conservation $\rho_{t}+\nabla\cdot(\rho{\bf u})=0$ and
assuming that the normal component is stationary in our reference frame (${\bf
w}=-{\bf u}$) we obtain

$\frac{\partial\phi}{\partial
t}+\frac{1}{2}u^{2}+\mu=-\zeta_{3}\frac{\partial\rho}{\partial t}.$ (3)

To modify Eq. (2) according to Khalatnikov’s theory, we need to consider
instead

$\displaystyle i\frac{d\psi_{i}}{dt}=$ $\displaystyle-\sum_{j,j\neq
i}J_{ij}\psi_{j}+g\;(1-|\psi_{i}|^{2})\psi_{i}$ (4) $\displaystyle-
iK\psi_{i}\sum_{j,j\neq i}J_{ij}(\psi_{j}^{*}\psi_{i}-\psi_{j}\psi_{i}^{*}),$

where $K=\zeta_{3}$. To see the correspondence with Eq. (3) we use the
Madelung transformation $\psi_{i}=\sqrt{\rho_{i}}\exp[i\phi_{i}]$ in Eq. (4)
so that the imaginary part yields $\dot{\rho_{i}}=2\sum
J_{ij}\sqrt{\rho_{i}\rho_{j}}\sin\phi_{ij}$ and the real part becomes
$\dot{\phi_{i}}-\sum
J_{ij}\sqrt{\frac{\rho_{j}}{\rho_{i}}}\cos\phi_{ij}+g\;(1-\rho_{i})=-K\dot{\rho_{i}},$
where we denoted $\phi_{ij}=\phi_{i}-\phi_{j}$. If the continuous wavefunction
is discretized on a square grid with the lattice constant $\Delta$ and
$J_{ij}=1/2\Delta^{2}$ NN interactions, then $-\sum
J_{ij}\cos\phi_{ij}\approx-\frac{1}{2\Delta^{2}}(4-\frac{1}{2}(\phi_{i}-\phi_{i\pm
1})^{2}-\frac{1}{2}(\phi_{i}-\phi_{i\pm N})^{2})$ $\approx
const+\frac{1}{2}(\nabla\phi)^{2}\approx const+\frac{1}{2}u^{2}.$ Associating
$\mu=const+g(1-\rho_{i})$ and $\rho_{i}\approx\rho_{j}$ completes the
justification of the Khalatnikov’s correction to Eq. (2). We note that Eq. (4)
for $J_{ij}=J>0$ and NN interactions only (and the time-dependent $K$) were
used to study the non-equilibrium quenches De Wijn  _et al._ (2013); Tarkhov
(2020) in laser-induced melting of the charge-density-wave order in solid Zong
_et al._ (2019). The steady state is achieved when $\sum
J_{ij}\sqrt{\rho_{i}\rho_{j}}\sin\phi_{ij}=0$ which corresponds to the minimum
of the XY model if $\rho_{i}=1$ as the penalty in the Hamiltonian $H$
dictates. Now we consider the quenching dynamics of the system described by
Eq. 4 and determine the conditions that prevent the system from achieving its
GS.

Low-dimensional topologies of coupling matrix. The formation of domains and
topological structures on the route to GS is well understood for the classical
lattice models with $J_{ij}=J>0$ in the context of BEC formation (quench
dynamics) from a strongly degenerate gas of weakly interacting bosons Berloff
and Svistunov (2002); Weiler  _et al._ (2008); Sun  _et al._ (2017). In
toroidal 1D geometry, a global and persistent phase winding can spontaneously
form Das  _et al._ (2012). In 2D, the establishment of the long-range order
and a genuine condensate in GS is replaced by the persistence of the isolated
vortices according to the Berezinskii–Kosterlitz–Thouless (BKT) transition
scenario Berezinskii (1972); Kosterlitz and Thouless (1973). We conclude that
the existence of fully or partially disjoined toroidal 1D or planar subsets in
the adjacency matrix ${\bf J}$ prevents the system from achieving its GS. In
toroidal 1D configurations, the global phase windings characterize the excited
states of the system, while topological protection prohibits further
transition to the GS. Strongly converging vortex-antivortex pairs (VAPs) in 2D
can eventually annihilate, bringing the system to the vortex-free GS. However,
the larger is effective separation of vortices, the less likely they interact,
the more likely the 2D subsystem remains in the excited state. The detailed
analysis of these two scenarios is presented in the Supplementary Information
(SI), where we also studied their strong dependence on nonlinearity $g$ and
the rate of dissipation $K$. Both of these parameters can vary from system to
system and often can be controlled within one system (e.g. using Feshbach
resonance in ultracold atomic BECs to change $g$ Bauer  _et al._ (2009) or the
cooling rate to control $K$).) The effective vortex core size depends on $g$
as $\xi\sim(J/g)^{1/2}$ which sets a competition between two opposing effects.
Stronger nonlinearity leads to stronger interactions and, therefore, more
efficient vortex annihilation. On the other hand, sufficiently large $g$
reduces the scale at which isolated vortices can still interact with each
other, which leads to freezing of plane vortices; see Fig. 1(c,d). Moreover,
for each size of the 2D subsystem, there is an optimal $g$ that maximizes the
success probability for finding GS; see Fig. 1(the inset of d).

Figure 1:  Ferromagnetic couplings $J_{ij}=1$. (a - c) The time snapshots of
the evolution of Eq. (4) on a 2D square $100\times 100$ lattice with NN planar
interactions from (a) the randomly distributed initial state with
$|\psi_{i}|=1$ and uniformly distributed $\phi_{i}$ (shown as black arrows) to
(c) the final steady-state with the persistence of isolated vortices. (d)
Average Hamiltonian time evolution for 100 different initial conditions and
various strengths of nonlinearity $g$. Label $P$ [$NP$] implies planar NN
[non-planar] couplings discussed in the main text. Here $g_{1}=10^{2}$,
$g_{2}=10^{3}$, $g_{3}=2\cdot 10^{3}$, $g_{4}=6\cdot 10^{3}$. (inset of d) The
success probability of finding the ground (vortex-free) state during the
evolution for $K=100$ as a function of $g$. Shading represents the variability
of the result for different runs.

The distinctive feature of the quench dynamics of 3D systems is the formation
of vortex filaments emerging from the local symmetry breaking Berloff and
Svistunov (2002). The closed filaments – vortex rings– are the torus-shaped
regions of phases rotating by $2\pi$ around a closed contours. When the system
evolves according to Eq. (4), the vortex rings decay bringing the system to
the true GS without any topological defects persisting at the final state.
This suggests that we can help the system achieve its GS by temporarily
introducing non-planar (quasi-3D) couplings that prevent the appearance of the
global winding and turn isolated VAPs into vortex rings to enable them to
shrink annihilate. To describe this procedure, we consider the graph
configuration of the interacting spins $G=\big{(}U,V\big{)}$, as shown in Fig.
2. The nodes of $G$ are characterised by $\psi_{i}\in U$ while
$(\psi_{i},\psi_{j})\in V$ represent the edges; see Fig. 2(a). ${\bf J}$
becomes the adjacency matrix of $G$. We define the coupling scheme
$S^{k}=\big{(}U^{k},V^{k}\big{)}$ of the individual node $\psi_{k}$ as the
subgraph of G containing all nodes $\psi_{j}\in U^{k}$ adjacent to $\psi_{k}$
and all edges $(\psi_{k},\psi_{j})\in V^{k}$ between them; see Fig. 2(a). We
say that graph $G$ and the corresponding ${\bf J}$ are generated by the scheme
${\bf S}$ if each individual node of $G$ can be described by ${\bf S}$.

Figure 2:  (a) Representation of interactions of individual node $\psi_{i}$
according to ${\bf S}$. Red [blue] edges/cells represent ferromagnetic [anti-
ferromagnetic] interactions with the central cell. (b) An illustration of one
non-planar coupling used in this paper (brown cells represent the additional
non-planar interactions). (c) Schematic representation of interacting VAPs on
a quasi-3D lattice with non-planar couplings.

In 2D square geometry and NN interactions, $G$ is planar. To implement a
quasi-3D geometry we introduce additional non-planar interactions shown in
Fig. 2(b). With these couplings, $G$ becomes non-planar and VAPs are replaced
with vortex filaments, as schematically illustrated in Fig. 2(c). Under the
action of non-energy preserving dissipation in Eq. (4), the radius $R(t)$ of
3D vortex rings obeys $R^{\prime}(t)=-\frac{1}{2}a(K)\log(8R(t))/R(t),$ where
$a(K)$ is the rate of the vortex ring dissipation Groszek  _et al._ (2018);
Berloff and Youd (2007). An average radius of the ring arising during the
quench with linear dimensionless domain size $L/\xi$ is $R_{0}=L/8\xi$. To
estimate the characteristic time for VAPs to pass through the ferromagnetic
quasi-3D system, we place VAP with separation $2R_{0}$ and numerically
calculate $a(K)$ from Eq. (4) for different values of $K$ and $g$ as
$a(K)=R_{0}^{2}/{\log(8R_{0})T},$ where $T$ is the dimensionless time it takes
to completely annihilate VAP. Here, we used the result that the square of the
vortex line length decays linearly with time Berloff and Youd (2007). We
determined the phenomenological form of $a(K)$ as $a(K)=1.35K+50.5$ (the fit
of this formula to the numerical data is shown in Fig. S.5). The VAP
dissipation rate $a(K)$ does not depend on $g$, which supports our assumption
that VAP motion in quasi-3D geometry resembles the motion of the vortex rings
in 3D.

This suggests that to help the system to reach the ground state, we can
introduce additional non-planar couplings $B_{ij}$ (see Fig. 2(b)) by
replacing planar $J_{ij}$ with $J_{ij}+\eta(t)B_{ij},$ where $\eta(t)$ is a
monotonically decaying function of $t$, chosen here as the Fermi function
$\eta(t)=(\exp[5(t-t_{0}+1)]+1)^{-1}.$ For sufficiently long attenuation times
$t_{0}$ ($\sim T$), all VAPs annihilate and bring the system to the GS with
100% probability; see Fig. 1(d) for the characteristic time evolution of $H$.

We fulfilled the first task on the proposed route to find the GS of the
general coupling matrix that consists of using quench dynamics while changing
the system’s topology to allow the topological defects to dissipate and
annihilate while propagating through the system. However, for the general
(sign and magnitude changing) form of adjacency matrix ${\bf J}$ the potential
energy landscape becomes rough so that topological structures are likely to
become trapped in one of many local minima (traps), manifesting massive
complexity of the problem of finding GS for spin glasses. However, not every
glassy system traps excitations. For these systems, the vortex phase
imprinting (or quench with vortices spontaneously generated) brings the system
to the GS, as we illustrate in the SI. So next, we develop an annealing scheme
that allows excitations to escape the traps. We consider the local mobility of
VAPs in the neighborhood of each node. If VAPs are mobile, no changes to
couplings are necessary during the annealing phase. If VAPs get trapped around
some $\psi_{i}$, the task is to find the minimum set of edges connected to
$\phi_{i}$ that allow VAPs to pass when changed. To fully characterize all
possible states, we focus on physically relevant couplings that are generated
by schemes ${\bf S}$ shown in the inset to Fig. 3(a) and $J_{ij}=\pm J$. That
gives us $2^{16}+2^{8}+2^{8}$ possible coupling configurations.

Figure 3: (a) Probability distribution for VAPs being trapped in different
spin glass lattice models. Mobility of VAPs is analysed by numerical
integration of Eq. (4) with $K=100$, $g=500$ on $40\times 40$ doubly-periodic
2D square lattices using all possible $2^{16}+2^{8}+2^{8}$ couplings ${\bf J}$
generated by schemes shown in the inset and using $50$ random VAP placements
for each scheme. (b) Energy states found by three different annealing schemes
starting from 80 random initial configurations and numerical integration of
Eq. 4. ”$A$” indicates the original glassy problem to be solved.
$(X\rightarrow A)$ shows the distribution of the found energy states using the
quenching process (if $X=A$), using annealing from the ferromagnetic state (if
$X=F$), and using the QVMA scheme (if $X=B$). Red [blue] cells in annealing
schemes represent ferromagnetic [anti-ferromagnetic] interactions with the
central cell.

The coupling matrices ${\bf J}$ have been constructed by the following rule:

$\displaystyle\hat{R}(Y)S^{(X,Y)}\hat{R}(X)=\begin{bmatrix}0&&1\\\ &\udots&\\\
1&&0\end{bmatrix}^{Y}\overbrace{\begin{gathered}\vspace{-0.1cm}\includegraphics[height=31.2982pt]{LUScheme.pdf}\end{gathered}}^{S}\;\begin{bmatrix}0&&1\\\
&\udots&\\\ 1&&0\end{bmatrix}^{X},$ (6)
$\displaystyle\hat{R}(P)=\begin{bmatrix}0&0&1\\\ 0&\hat{I}_{3\times 3}&0\\\
1&0&0\end{bmatrix}^{F(P)}\quad F(P)=\frac{P-P\;\rm{mod}\;2}{2}.$

Here $S^{(X,Y)}$ is a coupling scheme connecting $\psi{(X,Y)}\in G$ with its
neighbours with coordinates denoted by $X,Y\in\mathbb{N}$ (Eq. (5) uses one
example of such $S$).

Firstly, we need to determine the mobility of VAPs for different schemes. To
do this, we prepare each $40$ doubly-periodic lattice in the steady-state with
some low energy $H_{0}$ (this can be done by integrating Eq. (2) in imaginary
time or using the gain-dissipative networks Kalinin and Berloff (2018)). We
inject a VAP at random positions $(X_{i},Y_{i})$ by multiplying $\psi(X,Y)$ by
$v(X-X_{0},Y-Y_{0})v^{*}(X-X_{1},Y-Y_{1})$ where
$v(X,Y)=[X+iY]/[X^{2}+Y^{2}+\xi]$ and
$(X_{1}-X_{0})^{2}+(Y_{1}-Y_{0})^{2}=10^{2}$ Berloff (2004). We numerically
integrate Eq. (4) until the steady state is found with energy $H_{i}$. If
$H_{i}<H_{0}$, we redefine as $H_{0}=H_{i}$ for the rest of the calculations.
For each random VAP placement, VAPs either become trapped ($H_{i}>H_{0}$ for
some $i$ ) or not ($H_{i}\leq H_{0}$ for all $i$). Figure 3(a) depicts the
probability of VAPs to be trapped somewhere on the lattice as the function of
$\sum_{i,j}J_{ij}/N^{2}D$ remarkably illustrating high probability (around
60%) of VAPs being trapped for any configuration away from either
predominantly ferromagnetic or antiferromagnetic couplings.

Secondly, to demonstrate our quantized vortex mediated annealing (QVMA)
approach, we selected the set $A$ of worst coupling schemes denoted by $A_{i}$
with $D=16$ that trap VAPs with 100% probability. We denote by
$B=\\{\\{B_{j}\\}_{i}\\}$ the sets of coupling schemes where VAPs are free to
move while each $B_{j}\in\\{B_{j}\\}_{i}$ differs from $A_{i}\in A$ by a
single coupling strength. In total there were $|A|=3545$ and $|B|=5418$
coupling schemes with this property. We study the quench dynamics according to
Eq. (4) starting from $\\{\\{B_{j}\\}_{i}\\}$ and slowly (on the timescale of
the vortex annihilation $T$) transform couplings from one such set (when
several sets satisfy this property we will take the one with the minimum
energy configuration reached at the end of the annealing). Such annealing
transformation is realised by $J_{ij}=J^{B/F}_{ij}$ if $t<T_{i}$,
$J_{ij}=(1-(t-T_{i})/\Delta T)J^{B/F}_{ij}+(t-T_{i})J^{A}_{ij}/\Delta T$ if
$T_{i}\leq t\leq T_{f}$, and $J_{ij}=J^{A}_{ij}$ if $t>T_{f}.$ Here, $\Delta
T=T_{f}-T_{i}$, $T_{i}$ and $T_{f}$ represent the beginning and the end moment
of annealing: $T_{i}\approx T$ is the time required to reach a steady-state,
and $T_{f}\approx 15T$. $J_{ij}^{A[B/F]}$ are the couplings given by scheme
$A_{i}$ [$B_{i}$ or all ferromagnetic couplings.] For such set $A$, we
calculated the fraction of matrices for which the lowest energy state within
0.5% of the corresponding spectral width was found using 80 random initial
states per each matrix and three methods: (1) QVMA ($B\rightarrow A$); (2)
quench dynamics without annealing ($A\rightarrow A$) and (3) from the
ferromagnetic state with annealing ($F\rightarrow A$). QVMA ($B\rightarrow A$)
achieved minimization of 76% of all $3545$ matrices, compared with 13% for
$F\rightarrow A$ and 44% for $A\rightarrow A$. The success probabilities of
finding various energy states from 80 different initial configurations are
depicted in Fig. 3(b) for three different schemes $A_{i}$ demonstrating a
striking improvement of QVMA over other methods.

In summary, we argued that two main characteristics of the coupling graph of
the universal complex matter field lattice systems prevent it from obtaining
GS during the dissipative quench dynamics: (i) the presence of low (one- or
two-) dimensional subgraphs that leads to either global winding or isolated
noninteracting and, therefore, non-dissipative topological defects and (ii)
the existence of trapping regions for mobile excitations. We have shown that
both characteristics can be efficiently annealed from the modified coupling
graphs with favourable characteristics bringing the original system to the
ground state with high probability. Such modified coupling graphs satisfy the
following properties: they are (a) non-planar (to facilitate quasi-3D vortex
evolution); (b) do not trap vortex excitations; (c) have a minimum difference
in couplings with the original graph. This is the first time the engineered
mobility of quantized vortices was used as a mediator in the system’s
transition to the GS. Our work suggests that bringing the graph topology,
quench dynamics, and topological defects together may open new routes in
global physically enhanced optimization.

##  I Supplementary Information: Quantized Vortices Mediated Annealing

###  I.1 Role of nonlinearity and dissipation in one- and two-dimensional
graph topologies

Here we consider 1D and 2D periodic lattices of size $N=100$ and $N\times
N=100\times 100$ respectively, with purely ferromagnetic nearest neighbour
interactions. We define the winding number $W$ of phase configuration as the
total number of times that phases rotate by $2\pi$ in positive (counter-clock)
or negative (clockwise) directions. We numerically integrate Eq. (M.4) of the
main text starting with $|\psi_{i}|=1$ and randomly distributed phases of
$\psi_{i}$. The quenching dynamics in the one-dimensional periodic lattices
can result in the ground state with zero phase differences between the
neighboring sites so that $W=0$ or in the formation of the global phase
winding $W\neq 0$ as shown in Fig. 4(c). Phase configurations with $W\neq 0$
characterize the excited states of the system, and topological protection
prohibits their further transition to the ground state. We studied the
numerical evolution of one thousand 1D spin-lattices (each of size $N=100$)
with various $g$ and $K$ system parameters. Figure 4(a) shows that the
formation of the excited states with irreducible nonzero $W$ value strongly
depends on nonlinearity of the self-interactions $g$ and the dissipation level
$K$ of the system. Attainment of the final steady-state is estimated using the
order parameter

$O(t)=\bigg{|}\frac{1}{N}\sum_{i=1}^{N}\psi_{i}(t)\bigg{|}^{2}.$ (7)

The final steady-state is characterized by the stabilization of $O(t)\leq 1$
value in time. $O(t)=1$ corresponds to the ground state. As shown in Fig. 4
(a), the formation of stable excited states with $W\neq 0$ decreases with the
increase of $g$ in 1D systems. The local symmetry breaking occurring in 2D
lattices leads to the formation of plane vortices from domains with
independent choices of broken symmetry.

Strongly converging vortex-antivortex pairs can eventually annihilate,
bringing the system to the vortex-free ground state. The rate of the system
transition to the ground state strongly depends on their separation distances
which are characterized by the vortex core sizes $\xi$ that are set by the
balance between the interactions between the sites that scales as $J/\xi^{2}$)
and the potential energy of self-interactions characterised by $g$. Therefore,
$\xi\sim(J/g)^{1/2}$.

There is a competition between two opposing effects: stronger nonlinearity
leads to stronger interactions and therefore more efficient vortex
annihilation. On the other hand, sufficiently large $g$ values reduce the
scale at which isolated vortices can still interact with each other. For
sufficiently large $g$ values and fixed lattice size this results in the final
state with the presence of ”frozen” plane vortices, as shown in Fig. M.1 (c)
of the main text.

We verify this by starting from the set of 100 randomly distributed initial
conditions with $|\psi_{i}|\equiv 1$, and numerically evolving Eq. (M.4) (from
the main text) up to the stabilization of $H(t)$ dependence in time. Averaging
$H(t)$ dynamics over all 100 initial conditions, we provide $\langle
H(t)\rangle$ distribution shown in Fig. M.1 (d). Comparing the final state
Hamiltonian with the expected value of the ground state
$H_{GS}=-\sum_{ij}J_{ij}$, we conclude whether the ”frozen” vortices remain in
the system by the end of the quenching process or not. Figure M.1 (d inset)
depicts the obtained success probability distribution for achieving the
vortex-free ground state for different magnitudes of nonlinear self-
interactions $g$ in the $100\times 100$ system.

Note that the periodic boundary conditions we used endow the lattice with a
toroidal topological structure that has zero Euler characteristic $\chi$.
According to the Poincaré–Hopf theorem, the sum of the singularity indices
should satisfy the expression $\chi=\sum I_{i}=0$, where $I_{i}={\rm sign}W$.
Therefore, local violation of the $U(1)$ symmetry leads to the formation of an
equal number of vortices and anti-vortices, which makes the annihilation of
all vortex-anti-vortex pairs fundamentally achievable.

Figure 4: (a) Probability distribution of Winding Numbers obtained by
numerical integration of Eq. (M.4) (of the main text) on 1D lattice of size
$N=100$ with NN interactions and various $g$ and $K$ parameters with the set
of $1000$ randomly distributed initial states ($|\psi_{i}|=1$ and random
phases $\phi_{i}\in[0,2\pi]$). (b) The phase representations of 10
horizontally oriented initial states. (c) The phase representations of
different final states (depicted are only 10 of 1000).

###  I.2 Time evolution of energy during the quench process in glassy systems
with and without annealing

The energy time dynamics with and without annealing shows different stages of
the transformation. When spin-glass system undergoes the quench process, the
energy monotonically decreases in time and achieves the excited state as the
localised excited states dissipate and get trapped by the potential energy
landscape. When QVMA is implemented, the energy of the system may increase
when the annealing phase starts, but as the trapped localised states are
allowed to escape the system achieves a lower state at the end of the
annealing process. Figure 5 schematically illustrates this process.

Figure 5: Schematics representation of the quench dynamics of the average
energy of the system. Time evolution of the system marked as $A\rightarrow A$
does not involve annealing procedure. The other one (marked as $B\rightarrow
A$) starts with a coupling scheme $B$ that differs from $A$ by one coupling
strength while admitting free propagation of localised excitations through the
system. The adiabatic transformation of $B$ into $A$ takes place between
$t=T_{i}$ and $t=T_{f}$ as discussed in the main text. Here $T_{i}\approx T$
denotes the time required for the system to pass to the steady state and
$T_{f}\approx 15\cdot T$.

###  I.3 Effect of vortex imprinting on transition to a lower energy state

Here we show that by introducing the vortex-antivortex pair we can bring a
glassy system into a lower energy state. Such imprinting can be realized
experimentally, for example, using a laser beam White  _et al._ (2014). Such
imprinting is accompanied by an energy injection, which may be enough to
overcome the local potential barrier and bring the system into a lower energy
state. To demonstrate this, we numerically integrate Eq. (M.4) from the main
text with $K=100$, $g=500$, and $J_{ij}$ as shown in the inset to Fig. 6(f)
starting with $\psi_{i}(t=0)=1$. During the time evolution the system settles
into an excited steady state with energy $H_{1}$. At time $T_{0}$ we introduce
vortex-antivortex pair as discussed in the main text. This leads to a spike in
the energy followed by a slow relaxation to the true ground state with energy
$H_{0}$ as Fig. 6 illustrates.

Figure 6:  Time-evolution of the energy (see Eq. (M.1) (from the main text) of
spin-glass of size $N\times N=60\times 60$ governed by Eq. (M.4) (of the main
text) with couplings given in the inset (f) from uniformly distributed initial
state $\psi_{i}\equiv 1$ and $K=100,g=500$. The consecutive time snapshots of
phases following the vortex pair injection at $t=T_{0}$ are given in (a - e).

###  I.4 Vortex trapping probability for other coupling rules

Clearly, Fig. M.3 of the main text does not exhaust all possible types of
couplings. The rule we have chosen for constructing matrix $J_{ij}$, described
by Eq. (M.5) of the main text, is only one of many possible rules. However,
the general characteristics of the vortex trapping probability we discussed in
the main text are quite robust: the trapping probability saturates at about
60%-70% away from the couplings that are predominantly ferro- or
antiferromagnetic as Fig. 7(a) illustrates. The coupling strengths represented
in Fig. 7(b) give another rule for constructing matrix $J_{ij}$. Here, the
edges of same color have the same ferromagnetic or anti-ferromagnetic coupling
strengths. For each of $2^{8}$ possible choices of $J_{ij}\in\\{-1,1\\}$
satisfying described rule, the system evolution from 50 different initial
conditions with a randomly located single vortex-antivortex pairs separated by
distance $20$ on $40\times 40$ lattice sites was investigated by the numerical
integration of Eq. (M.4) with $K=100$ and $g=500$.

Figure 7: (a) The probabilities of vortex trapping in spin-glasses of size
$N\times N=40\times 40$ obtained by numerical numerical evolution of Eq. (M.4)
of the main text with the structure of interactions shown in (b) and
$K=100,g=500$. Here, $D=8$ is the degree of each node, different colors
correspond to different ferromagnetic or anti-ferromagnetic couplings, and
edges with the same color have the same coupling strength. Probabilities has
been calculated via $50$ vortex-antivortex pairs (VAPs) imprinting for each of
$2^{8}$ possible $J_{ij}$ structures. Detailed description of how this
distribution was obtained is given in the last part of the main text.

###  I.5 Evaluation of the vortex dissipation

Figure 8: Vortex-antivortex pairs characteristic dissipation rate $a(K)$
calculated for different values of $g$ on the (a) ferromagnetic quasi-3D
lattice (generated by scheme represented in Fig. M.2(b) of the main text) and
(b) glassy quasi-3D lattice (generated by scheme depicted in the (b inset)).
The red dotted line represents the approximation fit of the $a(K)$ dependence
cited in the main text which was obtained by numerical evaluation of the
characteristic time $T$ required for injected VAPs governed by Eq. (M.4) to
converge and annihilate.

To quantify the dynamics of vortices on spin glasses, we evaluate the VAPs
decay rate dependence $a(K)$ for glassy spin systems using the procedure
presented in the main text for ferromagnetic couplings; see Fig. 8(a). The
scheme with additional long-range non-planar interactions represented in Fig.
8(b inset) endow the lattice with quasi-3D geometry. The lattice generated by
this scheme does not contain any trapping regions for vortices. Therefore, we
expect that nothing will prevent the vortices from passing through the
lattice, and the system has to reach its global minimum eventually. Using the
expression $a(K)=R_{0}^{2}/{\log(8R_{0})T}$ derived in the main text, we
estimate the characteristic time for the vortices to pass through the glassy
system placing VAPs at a distance of $2R_{0}$ from each other. Here $R_{0}$
denotes the average radius of the vortex ring arising in the purely
ferromagnetic 3D system. Evolving Eq. (M.4) of the main text in time we
calculate the time $T$ required for the vortices to annihilate. Finally, we
calculate $a(K)$ for several $g$ values; see Fig. 8(b). We conclude that the
vortex decay rate dependence on dissipation level remains close to linear for
glassy quasi-3D systems in the absence of traps. However, the evolution
proceeds much slower in comparison with the case of purely ferromagnetic
systems considered previously and shown in Fig. 8(a).

###  I.6 Effect of the distance between trapping regions on the vortex
dynamics

We studied the influence of the separation between the vortex traps on the
capacity of vortices to avoid them. The simplest vortex trap can be realized
by the insertion of the frustrated square loop into the lattice. We will call
such frustrated or non-frustrated square loops represented in Fig. 9(b inset)
plaquettes. Twelve internal and external coupling strengths of the inserted
plaquette are depicted in Fig. 9(b inset), and coupling strength can take
either $+1$ or $-1$ value. We denote by $M$ the total sum of the internal and
external plaquette coupling strengths. Starting from the ground state of the
purely ferromagnetic system with an imprinted VAP, we insert two plaquettes of
the same structure at a distance $r$ from each other between vortices. VAPs
propagate through the system and enter the plaquettes region. There they can
either be trapped or pass freely. Vortex trapping rate can be estimated by
comparing the final state with the system’s global minimum. We evolve the Eq.
(M.4) with $K=200$ and $g=500$ numerically with all possible $2^{12}$
configurations of inserted plaquettes and different separations $r$ between
them. The probabilities of vortices to be trapped by the inserted plaquettes
are shown in Fig. 9(a). The smaller is the distance between the plaquettes the
higher is the escape probability.

Figure 9: (a) The probabilities of VAPs to get trapped on two plaquettes
placed at a distance $r$ from each other, between the VAPS and along the
trajectory of their motion. The plaquettes are the frustrated or non-
frustrated square loops with the structure shown in (b inset). Colors depict
twelve internal and external coupling strengths of the inserted plaquette, and
each color can take either $+1$ or $-1$ value. Here $M$ is the total sum of
the internal and external plaquette coupling strengths. For short separations
$r$ the attraction of VAPs assists in the escape probability. (b) The
dependence of the system’s dissipation level $K$ on the trapping probabilities
for plaquettes separated by $r=16$ and $g=500$.

## References

  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Nielsen and Chuang (2002) M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” (2002).
  * Farhi  _et al._ (2001) E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 (2001).
  * Wilczek (1982) F. Wilczek, Physical review letters 49, 957 (1982).
  * Aharonov  _et al._ (1993) Y. Aharonov, L. Davidovich, and N. Zagury, Physical Review A 48, 1687 (1993).
  * Bernstein and Vazirani (1997) E. Bernstein and U. Vazirani, SIAM Journal on computing 26, 1411 (1997).
  * Jordan  _et al._ (2010) S. P. Jordan, D. Gosset, and P. J. Love, Physical Review A 81, 032331 (2010).
  * Moussa  _et al._ (2010) O. Moussa, C. A. Ryan, D. G. Cory, and R. Laflamme, Physical review letters 104, 160501 (2010).
  * Marandi  _et al._ (2014) A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, Nature Photonics 8, 937 (2014).
  * Nixon  _et al._ (2013) M. Nixon, E. Ronen, A. A. Friesem, and N. Davidson, Physical review letters 110, 184102 (2013).
  * Berloff  _et al._ (2017) N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. Töpfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis, Nature materials 16, 1120 (2017).
  * Cai  _et al._ (2019) F. Cai, S. Kumar, T. Van Vaerenbergh, R. Liu, C. Li, S. Yu, Q. Xia, J. J. Yang, R. Beausoleil, W. Lu, _et al._ , arXiv preprint arXiv:1903.11194 (2019).
  * Kaminsky  _et al._ (2004) W. M. Kaminsky, S. Lloyd, and T. P. Orlando, arXiv preprint quant-ph/0403090 (2004).
  * Buluta and Nori (2009) I. Buluta and F. Nori, Science 326, 108 (2009).
  * De las Cuevas and Cubitt (2016) G. De las Cuevas and T. S. Cubitt, Science 351, 1180 (2016).
  * Svistunov  _et al._ (2015) B. V. Svistunov, E. S. Babaev, and N. V. Prokof’ev,  _Superfluid states of matter_ (Crc Press, 2015).
  * Domb (2000) C. Domb,  _Phase transitions and critical phenomena_ (Elsevier, 2000).
  * Kevrekidis  _et al._ (2001) P. Kevrekidis, K. Rasmussen, and A. Bishop, International Journal of Modern Physics B 15, 2833 (2001).
  * Brazhnyi and Konotop (2004) V. Brazhnyi and V. Konotop, Modern Physics Letters B 18, 627 (2004).
  * Struck  _et al._ (2011) J. Struck, C. Ölschläger, R. Le Targat, P. Soltan-Panahi, A. Eckardt, M. Lewenstein, P. Windpassinger, and K. Sengstock, Science 333, 996 (2011).
  * Davis  _et al._ (2019) E. J. Davis, G. Bentsen, L. Homeier, T. Li, and M. H. Schleier-Smith, Physical review letters 122, 010405 (2019).
  * Gallina and Pastor (2020) D. Gallina and G. Pastor, Physical Review X 10, 021068 (2020).
  * Kadowaki and Nishimori (1998) T. Kadowaki and H. Nishimori, Physical Review E 58, 5355 (1998).
  * Hen and Sarandy (2016) I. Hen and M. S. Sarandy, Physical Review A 93, 062312 (2016).
  * Hen and Spedalieri (2016) I. Hen and F. M. Spedalieri, Physical Review Applied 5, 034007 (2016).
  * Khalatnikov (1965) I. Khalatnikov, Inc., New York (1965).
  * De Wijn  _et al._ (2013) A. De Wijn, B. Hess, and B. Fine, Journal of Physics A: Mathematical and Theoretical 46, 254012 (2013).
  * Tarkhov (2020) A. Tarkhov, (2020).
  * Zong  _et al._ (2019) A. Zong, A. Kogar, Y.-Q. Bie, T. Rohwer, C. Lee, E. Baldini, E. Ergeçen, M. B. Yilmaz, B. Freelon, E. J. Sie, _et al._ , Nature Physics 15, 27 (2019).
  * Berloff and Svistunov (2002) N. G. Berloff and B. V. Svistunov, Physical Review A 66, 013603 (2002).
  * Weiler  _et al._ (2008) C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
  * Sun  _et al._ (2017) Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, Physical review letters 118, 016602 (2017).
  * Das  _et al._ (2012) A. Das, J. Sabbatini, and W. H. Zurek, Scientific reports 2, 1 (2012).
  * Berezinskii (1972) V. Berezinskii, Sov. Phys. JETP 34, 610 (1972).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Bauer  _et al._ (2009) D. M. Bauer, M. Lettner, C. Vo, G. Rempe, and S. Dürr, Nature Physics 5, 339 (2009).
  * Groszek  _et al._ (2018) A. J. Groszek, D. M. Paganin, K. Helmerson, and T. P. Simula, Physical Review A 97 (2018), 10.1103/physreva.97.023617.
  * Berloff and Youd (2007) N. G. Berloff and A. J. Youd, Physical review letters 99, 145301 (2007).
  * Kalinin and Berloff (2018) K. P. Kalinin and N. G. Berloff, New Journal of Physics 20, 113023 (2018).
  * Berloff (2004) N. G. Berloff, Journal of Physics A: Mathematical and General 37, 1617 (2004).
  * White  _et al._ (2014) A. C. White, B. P. Anderson, and V. S. Bagnato, Proceedings of the National Academy of Sciences 111, 4719 (2014).

*