uncertainties on the AstroSpectroSB1 solution – which also enter into our
joint fit – are underestimated. In calculating the astrometric solution, the
uncertainties of the epoch-level data are inflated by multiplying by a
constant, $c=\sqrt{\chi^{2}\nu^{-1}\left(1-\frac{2}{9\nu}\right)^{-3}}$, after
which $F_{2}=0$ (Halbwachs et al., 2022; Pourbaix et al., 2022).

##  Appendix C TESS light curves

The eleanor light curves of Gaia BH2 (TIC 207885877) are shown in Figure 20.
We extract these using a $15\times 15$ pixel target pixel file, with the
background estimated in a $31\times 31$ pixel region. We use the systematics-
corrected flux in the default aperture chosen by eleanor; i.e., the aperture
with the lowest combined differential photometric precision (CDPP) on one-hour
timescales. The expected flux contamination from background stars in the
aperture is of order 10%.

The source is clearly not strongly variable on short timescales. Much of the
apparent variability is noise; some variability on longer timescales may be
the result of imperfect background subtraction and removal of systematics. The
light curve from sector 38 is noisier than the one from sector 11, mainly
because sector 38 used $3\times$ shorter exposure times.

Power spectra for both sectors are shown in the right panels of Figure 20. The
red line shows the result of smoothing the power spectra with a Gaussian
kernel with $\sigma=5\,\mu\rm Hz$. Several potentially significant peaks are
evident in the smoothed power spectra. We tentatively identify the peak with
$\nu=61\,\mu\rm Hz$ as the global asteroseismic parameter $\nu_{\rm max}$,
since it is the strongest feature evident in data from both sectors. We
emphasize that the identification is tentative since the time baseline and
signal-to-noise ratio of the TESS data is still insufficient to unambiguously
constrain $\nu_{\rm max}$ in a source of this apparent magnitude (e.g. Stello
et al., 2022). The implied surface gravity is $\log g\approx\log
g_{\odot}+\log\left(\nu_{{\rm max}}/\nu_{{\rm
max,\odot}}\right)+\frac{1}{2}\log\left(T_{{\rm eff}}/T_{{\rm
eff},\odot}\right)\approx 2.68$, where $g$ is in cgs units and $\log
g_{\odot}\approx 4.43$, $\nu_{\rm max,\odot}=3090\,\rm\mu Hz$, and $T_{\rm
eff,\odot}=5777\,\rm K$.

Figure 20: TESS data for Gaia BH2. Left panels show normalized light curves
from the full frame images, with an exposure time of 1800s for sector 11 and
600s for sector 38. Right panels show power spectra, and red lines show power
spectra smoothed with a Gaussian kernel with $\sigma=5\,\mu\rm Hz$. Vertical
dashed line shows a tentative detection of $\nu_{\rm max}$.


Nuclei with enhanced Schiff moments in practical elements for atomic and
molecular EDM measurements

This is a resource paper concerning enhancement of observable time-reversal
breaking effects by nuclear structure, aimed at AMO experimentalists. It’s
intended to support a white paper [1] by providing some orientation on what
can be said about some particular isotopes. Any conclusions are qualitative,
and the reader should consult and cite the primary references and reviews
rather than this arXiv alone.

J.A. Behr, TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T2A3 Canada

<EMAIL_ADDRESS>
##  1 Qualitative info on enhanced atomic and molecular EDM’s from nuclear
phenonmena

Many nuclei have low-lying excited states with same spin as the ground state
and opposite parity. A number of these have recently been suggested in the AMO
literature as possibly exhibiting static or vibrational octupole deformation.
A goal is to identify practical candidates with enhanced Schiff moments in
particular elements where advanced AMO techniques could offer increased
sensitivity to boosted atomic and molecular EDM’s. Isotopes are sought in
addition to 225Ra that are long-lived enough to use away from isotope
production facilities [2, 3] or are producible from generators like 223FrAg
[4, 5]. There is a similar need to choose quantifiably large Schiff moments in
crystals with stable rare earth elements to optimize experiments searching for
CP-violating couplings to long-wavelength axionlike dark matter [6]. Some of
the new parity doublets suggested are supported by indirect experimental info,
e.g. BE(3) strength measured or calculated in nearby even-even nuclei.
Unfortunately, other nuclear structure phenomena that would not enhance Schiff
moments can account for closely spaced parity doublets. Whether or not the
ground state and higher-spin states built on it look like one of a pair of
static or vibrational octupole bands, often there is limited information about
the opposite-parity excited states. Considerable effort over three decades has
gone into identifying definitive experimental nuclear signatures for static
and vibrational octupole deformation for many nuclei [7, 8, 9, 10].

Here we examine informally to what extent a number of these suggested low-
lying parity doublets would need more experimental nuclear structure
information and/or microscopic calculations to identify and quantify possible
Schiff moment enhancement. Whenever we list experimental evidence against
static octupole deformation for a few interesting AMO cases, we try to qualify
to what extent collective octupole vibrations can have similar Schiff moment
enhancement.

Note that Ref. [1] considers AMO methods less element-specific, thus possibly
applicable to any Z$\geq$86 nucleus 111This paper avoids the word ‘actinide,’
sensibly restricted to ‘f-element’ by AMO physicists– who spend years on their
relatively complex atomic and molecular spectra– but used more broadly in the
nuclear structure literature. determined to have Schiff enhancements,
including quantum logic spectroscopy [11], frozen noble gas hosts [12], and
ultracold molecule clock states [13]. The possible parity doublet in 229Pa has
energy splitting 60$\pm$50 eV [14], and experiments to resolve its status are
considered [15].

###  1.1 TRV enhancements from nuclear structure: background

A continuum between static octupole deformation and collective octupole
vibrations. Schiff moment-induced atomic and molecular EDM’s are thought to be
similar in systems with static octupole deformation and in systems with ‘soft’
(i.e. virtual) excitation of octupole vibrations that are more common, though
those octupole vibrations must be collective to helpfully enhance the Schiff
moment. Engel, Friar, and Hayes (EFH) [16] discuss to what extent these
extremes are different aspects of formally similar phenomena (see other
reviews like Ref. [9] for related discussions of octupole vibrations and
static deformation), and showed that in macroscopic models a Schiff moment
enhancement from collective octupole vibrations also scales with the inverse
of the parity doublet’s energy splitting, with the static deformation
$\langle\beta_{3}\rangle^{2}$ replaced by the virtual deformation average
$\langle{\bar{\beta}_{3}}^{2}\rangle$. EFH point out additional physical
effects that must be calculated for the vibrators, which could enhance the
Schiff moment further. Flambaum and Zelenvinsky apply a similar analytical
calculation of collective octupole vibrations to cases of interest 223Rn and
223,225Ra, finding similar Schiff enhancement compared to previous
calculations assuming static octupole deformation [17]. Auerbach et al. [18]
with microscopic QRPA find similarly large Schiff moment enhancement from
collective vibrations in the limit when their energies are small: they find in
Z$\geq$86 nuclei (with N and Z separated from the region of known static
octupole deformation) that octupole vibrational configurations are only small
components of the complex nuclear wavefunctions, suppressing the Schiff moment
enhancement. Note that configurations are assumed purely octupole in some
analytical approximation formulae, though Ref. [19] calculates such a
suppression factor. This author is unaware of any calculations that consider
interference between static and collective vibrational effects, which
presumably would be included naturally in future microscopic calculations
[20].

Nuclear magnetic quadrupole moments The possible enhancements of Schiff
moments between two and five orders of magnitude from static octupole
deformation have produced much attention in designing AMO experiments.
Additionally, atomic and molecular EDM’s can be enhanced by other nuclear
effects by about an order of magnitude, potentially leading to exciting
discovery potential if the AMO experiment has sufficient sensitivity.
Contributions of the more ubiquitous nuclear magnetic quadrupole moment (MQM)
to atomic and molecular EDM’s are free from atomic shielding, so have been
reliably calculated [21, 22]. MQM’s provide order-of-magnitude enhancement for
any nucleus with quadrupole deformation and nuclear spin $\geq$ 1, and if the
atomic or molecular state has nonzero angular momentum. Experimental
difficulties from the nonzero electron spin needed can be mitigated by
choosing systems with small magnetic g-factor, allowing in principle long
coherence times, and its known existence in many stable nuclei is driving
dedicated experimental efforts [23, 1].

TRV matrix elements That all atomic EDM enhancements from nuclear structure
also require difficult-to-calculate matrix elements of parameterized TRV
hadronic interactions was made apparent in Ref. [24], with possible
microscopic theoretical improvements of the spatial wavefunctions necessary
highlighted in Ref. [20]. Some analytical formulae include an estimate of one
TRV matrix element from a one-body effective potential [19, 22]. There are no
known helpful experimental benchmarks for the most important operator
$\sigma\cdot p$ in heavy nuclei.

###  1.2 Observables

Only two observables are qualitatively discussed here– consult the reviews for
complete information.

Charge radii changes with neutrons Adding a neutron to a nucleus increases its
charge radius. For a given element, adding a neutron to an odd-N nucleus
usually increases the charge radius by less than adding a neutron to an even-N
nucleus– intuitively, this is reasonable when the extra neutron simply pairs
in spin with the odd valence neutron in the same orbital, rather than
populating another orbital or driving a change in collective properties.
Octupole deformation is a possible explanation for unusual isotopes with
inverted differences in neutron-induced change of charge radii (“IDNCR”) in
Rn,Fr,Ra [25, 26, 31], with measurements further delineating IDNCR at higher
neutron number [27, 28] and in Ac [29]. Other explanations than octupole
phenomena in lower-Z nuclei are considered in Ref. [30], which also mentions
the possibility that collective vibrational effects might similarly produce
IDNCR. For some isotopes, IDNCR is the only experimental information, yet
thought it can be consistent with static octupole deformation or collective
octupole vibrations, unfortunately can’t be definitive.

BE(3)’s, BE(1)’s A possible smoking gun is electric octupole transition
strength B(E3) between the opposite-parity bands, but AMO people should recall
the long-wavelength expansion works really well in nuclei (r/$\lambda$ $\sim$
1 MeV/197 MeV-fm $\sim$ 0.5%). So if an E1 is allowed, the E3’s start smaller
by over 8 orders of magnitude and are usually considered unmeasurable. Thus
measurements of E3 strength between the ground $0^{+}$ state and excited 3-
state in nearby even-even nuclei are possible and useful [32, 33]. Strong
B(E1)’s are expected between the opposite-parity bands of static octupole
deformation, while weak B(E1)’s are often considered to be consistent with
collective octupole vibrations: since there can be other reasons for either
extreme, B(E1)’s are not considered definitive unless there are thorough
theory considerations.

Odd-odd nuclei Note that nuclei with both odd N and odd Z can have octupole
deformation [34], possibly doubling the number of cases with Schiff
enhancement. However, the nuclear structure is so much more complicated than
odd-even nuclei that it’s much less likely the Schiff enhancement could be
definitively determined, either experimentally or theoretically.

##  2 Examples

Here are considerations of isotopes of AMO interest.

###  2.1 Rare earths

A number of rare earth elements have been laser-cooled and trapped, and
several isotopes of such elements with low-lying parity doublets are being
suggested for static octupole deformation [2, 3].

153Eu has been examined thoroughly. 152,154Eu are odd Z- odd N nuclei with
half-lives of 13 and 8 years that have structure and E1 strengths suggesting
octupole deformation, while they and stable 153Eu exhibit IDNCR consistent
with octupole deformation in the ground states [35]. Such evidence inspired
detailed experimental examination of excited states in 153Eu, and the
g-factors of the identified excited opposite-parity band are known to be about
3 times larger than the ground state band [36, 37]. Now, it’s natural for
perturbations such as Coriolis mixing to change properties of partner bands
differently– nearby configurations tend to have the natural parity of the
ground-state– thus certain properties don’t have to be identical.
Nevertheless, such different g-factors discourage an interpretation of this
doublet in 153Eu as exhibiting static octupole deformation, while relatively
large E1’s between the bands are hard to interpret in terms of octupole
vibrations [36, 37]. So 153Eu remains ambiguous, and needs further modelling
and calculation to determine whether its identified parity doublet has
octupole phenomena.

Other stable rare earths Several other stable rare earth isotopes have low-
lying parity doublets (161,163Dy, 165Er, 155Gd, 153Sm) [2, 3], while nearby
even-even nuclei with measured 3- states and BE(3) strength support the
possibility of collective vibrational octupole configurations in such nuclei
as well as in 153Eu [38]. A recent global calculation of particle-bound even Z
- even N nuclei predicts static octupole deformation in the ground states of
stable 146,148,150Nd and 150Sm [39], while a region in Z and N of radioactive
rare earth nuclei with much higher N than the line of stability has been
identified to have static octupole deformation, both experimentally and
theoretically [39, 9, 33]. These stable rare earth nuclei proposed for Schiff
moment enhancements need further experimental and/or microscopic theoretical
supporting evidence for static or collective vibrational octupole deformation
in their ground and excited states.

###  2.2  Z$\geq$86

AMO experiments have been considered for nuclei with supporting experimental
evidence for static octupole deformation and/or collective octupole vibrations
for ground and excited states of a parity doublet, including 225Ra, 223Fr,
225,227Ac, and 227,229Th (see Fig. 1).

Figure 1: Overly simplistic summary of Z$\geq$86 isotopes of possible AMO
interest that have experimental evidence suggesting extremes in the continuum
of static octupole deformation vs. collective octupole vibrations, either of
which can enhance TRV Schiff moments. Only the dominant configuration is
indicated– see text and references for any details. Darker blues have longer
half-lives (see nndc.bnl.gov).

Radium: The solid experimental evidence supporting static octupole deformation
in I=1/2 225Ra and I=3/2 223Ra is thoroughly documented in review articles [9,
8].

Francium: 223Fr is an illustrative case where experiment and theory suggest a
large enhancement of the Schiff moment, even though the nuclear wavefunctions
are thought to contain significant admixtures of static octupole deformation
and collective octupole vibrations. IDNCR exists for 223,225Fr [25, 40] (it
barely exists for 221Fr, a conclusion that does not change with a slightly
revised optical isotope shift [41]). Concerning excited states, it’s
instructive to note that 221Fr is a rather complex “transition” nucleus; e.g.,
its possible static octupole deformed excited parity-doublet band actually has
spins out of order in energy. This is understood as a case in neutron number N
between 219Fr, which has no special deformation, and 223Fr, which has “fairly
strong” evidence [42] from its structure and feeding from alpha- and beta-
decay for a static octupole deformation parity doublet, albeit with a
relatively large parity-doublet splitting of 161 keV. The experimental parity
doublet of 223Fr has been modelled including admixtures with static octupole
and collective octupole vibrations [42]; however, E1 transitions between the
223Fr bands were later found to be relatively weak, consistent with 223Fr
having collective octupole vibrations [43, 44, 45]. Ref. [19], which considers
theoretical admixtures of other configurations in all its calculations,
predicts the static octupole enhancement of EDM’s for 221Fr to be an order of
magnitude less than 223Fr and 225Ra. 225Fr has weak B(E1)’s inconsistent with
static octupole deformation but consistent with collective octupole vibrations
[45]. (Note Ref. [45] has a thorough pedagogical discussion of the B(E1)
phenomenology, including an updated figure of the B(E1)’s in Ref. [7].)

Actinium: Iπ=3/2- 227Ac has a low-lying 3/2+ excited state, reviewed in [7].
There is considerable evidence consistent with these being the lowest-lying
states of bands with static octupole deformation, including measurements of
g-factors. But there is also experimental evidence against static octupole
deformation for this parity doublet, including measured BE(1)’s an order of
magnitude smaller than other cases, and transfer reactions populating excited
states that can be explained by reflection-symmetric potentials [7]. 227Ac
could be a case with large collective vibrational octupole enhancement, and
EDM experiments involving actinium molecules [2] could also be done in 225Ac,
which has strong supporting evidence for static octupole deformation [7] and a
10-day half-life producible by a 229Th generator.

Thorium: Some EDM experimental techniques work best for nuclear spin I=1/2. In
addition to the well-known I=1/2 225Ra, 227Th likely has I=1/2 [46]. There is
strong experimental and theoretical support for static octupole deformation
with an identified parity double band for an I=1/2 ground state, and a large
Schiff moment enhancement [8, 9].

An explicit band of excited states is identified in 229Th as a vibrational
octupole partner of a ground state band, with energy splitting 146 keV [47,
48], and thus from the general expectations of Ref. [16] a full microscopic
calculation could show a large Schiff moment for 229Th. It is in this sense
that recent BE(1) measurements in 228Th supporting vibrational collective
octupole effects in 229Th are considered to support it as an enhanced atomic
or molecular EDM case [49], even though detailed consideration of the excited
states of 229Th has classified it as a ‘transition’ isotope with larger N than
nuclei with demonstrated static octupole deformation [50, 8].

Radon: IDNCR exists for 221,223Rn [25], but there is only one excited state
energy known for 221Rn and no excited state info for 223Rn. Favorable static
octupole arguments from early calculations [19] have been more recently
weakened by BE(3) measurements in nearby even-even nuclei [9], which instead
suggest collective vibrational octupole effects are possible [51]. To quantify
the Schiff enhancements would either way need detailed excited-state nuclear
structure measurements, though this has remained elusive in several
experimental attempts.

Summary All nuclei with reasonably definitive supporting experimental evidence
for static octupole deformation are radioactive, with half-lives on the order
of weeks or shorter. Longer-lived nuclei like 229Th and 227Ac, along with
stable rare earth isotopes, show evidence for collective octupole vibrations
that could also produce quantifiably large Schiff enhancements, if microscopic
calculations with sufficiently small uncertainty can be developed. Any
prospective parity doublet needs support from experimental observables to
define what fraction of the ground state and excited states have static or
collective vibrational octupole configurations. Microscopic theory
improvements of this aspect, as well as to quantify matrix elements of the
several TRV isospin- and range-dependent interactions [20], would enable a
move from the current exciting discovery phase to defining multiple
complementary experiments to extract different microscopic sources of TRV.

Acknowledgements P. Butler, J. Engel, N. Hutzler, A. Jayich, D. DeMille, A.
Vutha, W. Nazarewicz, A. Madden, G. Neyens, and G. Hackman all have had
helpful comments and/or orientation. This author is, of course, responsible
for any misinterpretations of the literature.

## References

  * [1] “Opportunities for Fundamental Physics Research with Radioactive Molecules,” N. Hutzler et al., white paper to appear soon
  * [2] “Electric dipole moments of atoms and molecules produced by enhanced nuclear Schiff moments,” V.V. Flambaum and V.A. Dzuba, Phys. Rev. A 101, 042504 (2020).
  * [3] “Enhanced nuclear Schiff moment in stable and metastable nuclei,” V.V. Flambaum and H. Feldmeier, Phys. Rev. C 101, 015502 (2020)
  * [4] “Francium dimer molecules: Atomic S-wave scattering lengths and E1 Matrix Elements,” J.A. Behr, D. DeMille, A. Teigelhöfer. TRIUMF Letter of Intent S2139LOI https://mis.triumf.ca/science/experiment/view/S2139LOI
  * [5] “Prospects for assembling ultracold radioactive molecules from laser-cooled atoms,” J. Klos, H. Li, E. Tiesinga, and S. Kotochigova, arXiv:2111.14896 (2021).
  * [6] “The Piezoaxionic Effect,” A. Arvanataki, A. Madden, K. Van Tilburg, arXiv:2122.11466 (2021).
  * [7] “Octupole Shapes in Nuclei,” Irshad Ahmad and Peter A. Butler, Annu. Rev. Nucl. Part. Sci. 1993 43 71
  * [8] “Intrinsic reflection asymmetry in atomic nuclei,” P.A. Butler and W. Nazarewicz, Rev. Mod. Phys. 68 349 (1996)
  * [9] “Octupole collectivity in nuclei,” P.A. Butler, J. Phys. G 43 073002 (2016)
  * [10] “Evolution of Octupole Deformation in Radium Nuclei from Coulomb Excitation of Radioactive 222Ra and 228Ra Beams,” P. A. Butler, L. P. Gaffney, P. Spagnoletti, K. Abrahams, M. Bowry, J. Cederkäll, G. de Angelis, H. De Witte, P. E. Garrett, A. Goldkuhle, C. Henrich, A. Illana, K. Johnston, D. T. Joss, J. M. Keatings, N. A. Kelly, M. Komorowska, J. Konki, T. Kröll, M. Lozano, B. S. Nara Singh, D. O’Donnell, J. Ojala, R. D. Page, L. G. Pedersen, C. Raison, P. Reiter, J. A. Rodriguez, D. Rosiak, S. Rothe, M. Scheck, M. Seidlitz, T. M. Shneidman, B. Siebeck, J. Sinclair, J. F. Smith, M. Stryjczyk, P. Van Duppen, S. Vinals, V. Virtanen, N. Warr, K. Wrzosek-Lipska, and M. Zielińska, Phys. Rev. Lett. 124, 042503 (2020)
  * [11] A. Jayich, in Ref. 1 Hutzler et al..
  * [12] “Orientation-dependent hyperfine structure of polar molecules in a rare-gas matrix: A scheme for measuring the elecron electric dipole moment,” A.C. Vutha, M. Horbatsch, and E.A. Hessels, Phys. Rev. A 98 032513 (2018); J.T. Singh, Hyp. Int. 240:29 (2019).
  * [13] “Electron Electric Dipole Moment Searches Using Clock Transitions in Ultracold Molecules,” Mohit Verma, Andrew M. Jayich, and Amar C. Vutha, Phys. Rev. Lett. 125 153201 (2020)
  * [14] “Electron capture decay of 58-min ${}^{229}_{92}$U and levels in ${}^{229}_{91}$Pa,” I. Ahmad, R. R. Chasman, J. P. Greene, F. G. Kondev, and S. Zhu, Phys. Rev. C 92 024313 (2015)
  * [15] J.T. Singh, in Ref. 1 Hutzler et al.
  * [16] “Nuclear octupole correlations and the enhancement of atomic time-reversal violation,” J. Engel, J. L. Friar, and A. C. Hayes, Phys. Rev. C 61, 035502 (2000)
  * [17] “Enhancement of nuclear Schiff moments and time-reversal violation in atoms due to soft nuclear octupole vibrations,” V.V. Flambaum and V.G. Zelevinsky, Phys. Rev. C 68 035502 (2003)
  * [18] “Nuclear Schiff moment in nuclei with soft octupole and quadrupole vibrations,” N. Auerbach, V. F. Dmitriev, V. V. Flambaum, A. Lisetskiy, R. A. Seńkov, and V. G. Zelevinsky, Phys. Rev. C 74, 025502 (2006);
  * [19] “Enhanced T-odd, P-odd electromagnetic moments in reflection asymmetric nuclei” V. Spevak, N. Auerbach, and V. V. Flambaum, Phys. Rev. C 56 1357 (1997).
  * [20] J. Engel, in Ref. 1 Hutzler et al.
  * [21] “Enhanced T-Nonserving Nuclear Moments,” W.C. Haxton and E.M. Henley, Phys. Rev. Lett. 51 1937 (1983)
  * [22] “Possibility of investigating P- and T-odd nuclear forces in atomic and molecular experiments,” O.P. Sushkov, V.V. Flambaum, and I.B. Khriplovich, Sov. Phys. JETP 60 873 (1984).
  * [23] “Polyatomic molecules as quantum sensors for fundamental physics,” N.R. Hutzler, Quantum Science and Technology 5 044011 (2020).
  * [24] “Correlating Schiff Moments in the Light Actinides with Octupole Moments,” Jacek Dobaczewski, Jonathan Engel, Markus Kortelainen, and Pierre Becker, Phys. Rev. Lett. 121, 232501 (2018)
  * [25] “Nuclear Radii and Unstable Isotopes,” E. W. Otten, in “Treatise on Heavy Ion Science,” edited by D. A. Bromley, Vol. 8 (Springer US, Boston, MA, 1989) pp. 598.
  * [26] “Mean square charge radii of radium isotopes and octupole deformation in the 220–228Ra region,” S.A. Ahmad et al., Nucl. Phys. A 483 284 (1988)
  * [27] “Laser spectroscopy of francium isotopes at the borders of the region of reflection asymmetry,” I. Budin$\check{\rm c}$evi$\acute{\rm c}$, J. Billowes, M. L. Bissell, T. E. Cocolios, R. P. de Groote, S. De Schepper, V. N. Fedosseev, K. T. Flanagan, S. Franchoo, R. F. Garcia Ruiz, H. Heylen, K. M. Lynch, B. A. Marsh, G. Neyens, T. J. Procter, R. E. Rossel, S. Rothe, I. Strashnov, H. H. Stroke, and K. D. A. Wendt Phys. Rev. C 90 014317 (2014)
  * [28] “Laser-spectroscopy studies of the nuclear structure of neutron-rich radium,” K.M. Lynch, S.G. Wilkins, J. Billowes, C.L. Binnersley, M.L. Bissell, K. Chrysalidis, T.E. Cocolios, T. Day Goodacre, R.P. de Groote, G.J. Farooq-Smith, D.V. Fedorov, V.N. Fedosseev, K.T. Flanagan, S. Franchoo, R.F. Garcia Ruiz, W. Gins, R. Heinke, Á. Koszorús, B.A. Marsh, P.L. Molkanov, P. Naubereit, G. Neyens, C.M. Ricketts, S. Rothe, C. Seiffert, M.D. Seliverstov, H.H. Stroke, D. Studer, A.R. Vernon, K.D.A. Wendt, and X.F. Yang, Phys. Rev. C 97 024309 (2018)
  * [29] “Search for octupole-deformed actinium isotopes using resonance ionization spectroscopy,” E. Verstraelen, A. Teigelhöfer, W. Ryssens, F. Ames, A. Barzakh, M. Bender, R. Ferrer, S. Goriely, P.-H. Heenen, M. Huyse, P. Kunz, J. Lassen, V. Manea, S. Raeder, and P. Van Duppen, Phys. Rev. C 100, 044321 (2019)
  * [30] “Charge radii in covariant density functional theory: A global view,” U. C. Perera, A. V. Afanasjev, and P. Ring, Phys. Rev. C 104 064313 (2021)
  * [31] “The breaking of intrinsic reflection symmetry in nuclear ground states,” G.A. Leander, R.K. Sheline, P. Möller, P. Olanders, I. Ragnarsson, and J. Sierk, Nucl. Phys. A388 (1982) 452
  * [32] “Studies of pear-shaped nuclei using accelerated radioactive beams,” L.P. Gaffney et al., Nature 497, 199 (2013).
  * [33] “Pear-shaped atomic nuclei,” P.A. Butler, Proc. R. Soc. A 476: 20200202.
  * [34] “Octupole deformation in odd-odd nuclei,” R. K. Sheline, Phys. Rev. C 37 423 (1988).
  * [35] “Octupole Deformation in the Odd-Odd Nucleus 1S2Eu,” R.K. Sheline and P.C. Sood, Progress of Theoretical Physics, 81 1057 (1989).
  * [36] “‘Lifetimes of states in the opposite-parity bands of 153 Eu: Recoil-distance measurements following Coulomb excitation,” J.F. Smith, M.W. Simon, R.W. Ibbotson, P.A. Butler, A. Aprahamian, A.M. Bruce, D. Cline, M. Devlin, G. D. Jones, P. M. Jones, and C. Y. Wu , Phys. Rev. C 58 3171 (1998).
  * [37] “Opposite-parity bands in 153Eu,” C. J. Pearson, W. R. Phillips, J. L. Durell, B. J. Varley, W. J. Vermeer, W. Urban, and M. K. Khan, Phys. Rev. C 49, R1239 (1994)
  * [38] “Reduced Electric-Octupole Transition Probabilities, B(E3; 0${}_{1}^{+}\rightarrow 3_{1}^{-}$)— An Update,” T. Kibedi and R.H. Spear, Atomic Data and Nuclear Data Tables 80 35 (2002).
  * [39] “Landscape of pear-shaped even-even nuclei,” Yuchen Cao, S. E. Agbemava, A. V. Afanasjev, W. Nazarewicz, and E. Olsen Phys. Rev. C 102, 024311 (2020).
  * [40] “Isotope Shifts, Spins and Hyperfine Structures of 118,146Cs and of Some Francium Isotopes,” A. Coc, C. Thibault, F. Touchard, H.T. Duong, P. Juncar, S. Liberman, J. Pinard, M. Carre, J. Lerme, J.L. Vialle, S. Buttgenbach, A.C. Mueller, A. Pesnelle, and The ISOLDE Collaboration, Nucl. Phys. A468 (1987) 1.
  * [41] “Efficient Collection of 221Fr into a Vapor Cell Magneto-optical Trap,” Z.-T. Lu, K. L. Corwin, K. R. Vogel, C. E. Wieman, T. P. Dinneen, J. Maddi, and H. Gould, Phys. Rev. Lett. 79 994 (1997).
  * [42] “Experimental and theoretical study of the nuclear structure of 223Fr,” R. K. Sheline, C. F. Liang, P. Paris, J. Kvasil, and D. Nosek, Phys. Rev. C 51, 1708 (1995)
  * [43] “The nuclear structure of 223Fr,” W. Kurcewicz, G. Løvhøiden, T.F. Thorsteinsen, M.J.G. Borge, D.G. Burke, M. Cronqvist, H. Gabelmann, H. Gietz, P. Hill, N. Kaffrell, R. A. Naumann, K. Nybø, G. Nyman, J. Rogowski, and the ISOLDE Collaboration, Nucl. Phys. A539 451 (1992)
  * [44] A critical correction to a wrong tabulation of a 223Fr E1 strength in Kurcewicz 1992 is in footnote 9 p. 136 of Burke 1997– the true E1 strength is 2 orders smaller. Note Sheline 1995 uses Kurcewicz 1992 without the correction.
  * [45] “Search for stable octupole deformation in 225Fr,” D.G. Burke, W. Kurcewicz, G. Løvhøiden, M.J.G. Borge, M. Cronqvist, H. Gabelmann, H. Gietz, P. Hill, N. Kaffrell, S. Mattsson, R.A. Naumann, K. Nybø, G. Nymang, J. Rogowski, G.L. Struble, T.F. Thorsteinsen, ISOLDE Collaboration, Nucl. Phys. A 612 (1997) 91.
  * [46] “The ﬁrst experimental evidence for the (M1 + E2) mixed character of the 9.2 keV transition in 227Th,” A. Kovalík, A.Kh. Inoyatov, L.L. Perevoshchikov, M. Ryšavý, D.V. Filosofov, P. Alexa, and J. Kvasil, Physics Letters B820 136593 (2021)
  * [47] “Nuclear Structure of 229Th,” E. Ruchowska, W. A. Plóciennik, J. Zylicz, H. Mach, J. Kvasil, A. Algora, N. Amzal, T. Bäck, M. G. Borge, R. Boutami, P. A. Butler, J. Cederkäll, B. Cederwall, B. Fogelberg, L. M. Fraile, H. O. U. Fynbo, E. Hagebo, P. Hoff, H. Gausemel, A. Jungclaus, R. Kaczarowski, A. Kerek, W. Kurcewicz, K. Lagergren, E. Nacher, B. Rubio, A. Syntfeld, O. Tengblad, A. A. Wasilewski, and L. Weissman, Phys. Rev. C 73 044326 (2006)
  * [48] “The nuclear structure of 229Th,” K. Gulda, W. Kurcewicz, A.J. Aas, M.J.G. Borge, D.G. Burke, B. Fogelberg, I.S. Grant, E. Hagebø b , N. Kaffrell, J. Kvasil, G. Løvhøiden, H. Mach, A. Mackova, T. Martinez, G. Nyman, B. Rubio, J.L. Tain, O. Tengblad, T.F. Thorsteinsen, ISOLDE collaboration, Nucl. Phys. A 703 45 (2002)
  * [49] “Direct measurement of the intrinsic electric dipole moment in pear-shaped thorium-228” M. M. R. Chishti, D. O’Donnell, G. Battaglia, M. Bowry, D. A. Jaroszynski, B. S. Nara Singh, M. Scheck, P. Spagnoletti, and J. F. Smith, Nature Physics 16 853 (2020)
  * [50] “Intrinsic Reflection Asymmetry in Odd-A Nuclei,” G.A. Leander and R.K. Sheline, Nucl. Phys. A413 375 (1984).
  * [51] “The observation of vibrating pear-shapes in radon nuclei,” P.A. Butler, L.P. Gaffney, P. Spagnoletti, J. Konki, M. Scheck, J.F. Smith, K. Abrahams, M. Bowry, J. Cederkäll, T. Chupp, G. de Angelis, H. De Witte, P.E. Garrett, A. Goldkuhle, C. Henrich, A. Illana, K. Johnston, D.T. Joss, J.M. Keatings, N.A. Kelly, M. Komorowska, T. Kröll, M. Lozano, B.S. Nara Singh, D. O’Donnell, J. Ojala, R.D. Page, L.G. Pedersen, C. Raison, P. Reiter, J.A. Rodriguez, D. Rosiak, S. Rothe, T.M. Shneidman, B. Siebeck, M. Seidlitz, J. Sinclair, M. Stryjczyk, P. Van Duppen, S. Vinals, V. Virtanen, N. Warr, K. Wrzosek-Lipska, and M. Zielinska, Nat. Commun. 10, 2473 (2019)


# Classical vs Quantum Annealing and Manifold Reduction in Soft-Spin
Minimizers of Ising Hamiltonians

James S. Cummins1, Hayder Salman2,3 and Natalia G. Berloff
<EMAIL_ADDRESS>Department of Applied Mathematics and Theoretical
Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United
Kingdom
2 School of Mathematics, University of East Anglia, Norwich Research Park,
Norwich, NR4 7TJ, UK, 3 Centre for Photonics and Quantum Science, University
of East Anglia, Norwich Research Park, Norwich, NR4 7TJ, UK

###### Abstract

We investigate the minimization of the Ising Hamiltonians, comparing the
dynamics of semi-classical soft-spin models with quantum annealing. We
systematically analyze how the energy landscape for the circulant couplings of
a Möbius graph evolves with increased annealing parameters. Our findings
indicate that these semi-classical models face challenges due to a widening
dimensionality landscape. To counteract this issue, we introduce the ‘manifold
reduction’ method, which restricts the soft-spin amplitudes to a defined phase
space region. Concurrently, quantum annealing demonstrates a natural
capability to navigate the Ising Hamiltonian’s energy landscape due to its
operation within the comprehensive Hilbert space. Our study indicates that
physics-inspired or physics-enhanced optimizers will likely benefit from a
blend of classical and quantum annealing techniques.

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 field $\psi_{i}=r_{i}\exp[i\theta_{i}]$ that has the
amplitude $r_{i}$ (referred to as ‘soft mode’) and phase $\theta_{i}$
(discrete e.g. $\theta_{i}=\\{0,\pi\\}$ or continuous ‘spin’) degrees of
freedom.

Optical parametric oscillator (OPOs) based coherent Ising machines (CIM) [2,
3, 4, 5], lasers [6, 7, 8], spatial light modulators [9], lattices of
polariton [10, 11], and photon condensates [12], Microsoft analogue iterative
machine [13] and Toshiba simulated bifurcation machine [14] can all minimize
the classical hard-spin Ising Hamiltonians (HSIH) $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
$s_{i}=(\cos\theta_{i},\sin\theta_{i})$) using soft-spin bifurcation dynamics
via the Aharonov-Hopf bifurcation [15]. 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.
Dimensionality annealing was proposed to bring advantages to physical
optimizers by bridging minima separated by significant energy barriers [16].
In parallel to these methods, quantum annealing is another approach to
minimize the HSIH. Despite numerous studies contrasting classical and quantum
methods, the limitations of currently available hardware and limitations of
simulating quantum systems classically has led to conflicting conclusions as
to whether a quantum advantage can potentially be realized using quantum
annealing [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].

Here, we analyse and contrast the annealing of semi-classical 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. 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 [25]. 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 annealing finding the global minima.

The ONNs based on laser operation are non-Hermitian systems that tend to
minimise losses on their route to coherence. The losses can be written as an
‘energy’ (‘cost’) function to be minimized. For instance, CIM’s energy
landscape to be minimised (in a classical limit) is

${E}=\frac{C}{4}\sum_{i=1}^{N}(p(t)-x_{i}^{2})^{2}-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}J_{ij}x_{i}x_{j},$
(1)

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 a
gradient descent as

${\dot{x}_{i}}=-\frac{\partial{E}}{\partial
x_{i}}=C(p(t)x_{i}-x_{i}^{3})+\sum_{i\neq j}J_{ij}x_{j}.$ (2)

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 network (HT) ${\dot{x}_{i}}=p(t)x_{i}+\sum_{i\neq j}J_{ij}x_{j}$
also used [26] for Ising Hamiltonian minimization [27]. Another approach uses
the second-order resonance to project the XY onto the Ising dynamics [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, what are the distinguishing features
between quantum and classical annealing?’. Answers often rely on the coupling
matrix $\mathbf{J}$. An instructive problem 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 the issues of statistical convergence [25]. We comprehensively explore
ONN dynamics on circulant coupling matrices, revealing evolution patterns
dictated by the eigenvalues of $\mathbf{J}$ and the spectral gap. We identify
certain unoptimizable regions using ONNs in classical annealing limits and
suggest improvements via constraining soft amplitudes to their mean value.

Cyclic graphs with $N$ nodes are characterized by an $N$x$N$ circulant
coupling matrix $\mathbf{J}$, constructed through cyclical permutations of any
$N-$vector. These graphs inherently have vertex permutation symmetry,
signifying boundary periodicity and uniform neighbourhoods. 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)]$ [29, 30, 31].

We consider the Ising Hamiltonian minimization on a particular form of a
cyclic graph – a Möbius ladder graph with even number of vertices $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,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$ 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$. When
$N/2$ is odd, $S_{0}$ is always the ground state with energy $E(J)=-(J+2)N/2$.
When $N/2$ is even, the $S_{0}$ configuration has energy $E_{0}(J)=(J-2)N/2$
and $S_{1}$ has energy $E_{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_{\rm crit}<J$]. The eigenvalues for Möbius ladders 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 Supp.
Inf. for the detailed derivation of the spectra).

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. (1)(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 realised 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. (1)(b,c), in agreement with the dynamical simulations
presented in Figs. (2)(b). 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. (1)(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}$ from the steady state
of equation on the dynamics of $X_{L}$ from Eq. (2), 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 (1)(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_{1}$). The bifurcation likelihood
for this state is near zero, as illustrated by Fig.(1)(e). 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}$, shown in Fig.(1)(d). Yet, increasing amplitudes bring
increased height of the energy barriers while preventing system transitions,
see Fig.(1)(f). Supp. Inf. Fig.(3) provides further context, depicting
attraction basins for diverse $p$ within a fixed $J_{\rm e}<J<J_{\rm crit}$.

Figure 1: (a-c) Schematic representation of the states realised by soft-spin
models of Eq. (1) on the Möbius ladder graphs with varying cross-ring
couplings (shown in red). (a),(b) and (c) depict states that map onto $S_{0}$,
$S_{1}$ Ising states for $N=8$ and $S_{1}$ state for $N=12$, respectively. (c)
uses a different node arrangement that illustrates the graph relationship with
the topology of the Möbius strip. The same colours are used to show equal
intensities; the larger sizes correspond to larger intensities. (d) Regions of
different global minima of Eq. (1): $E_{1}$ in the blue region and $E_{0}$ in
the pink region, in $J-p$ space for $N=8,C=1$. Two critical values of $J$ are
shown as solid black lines. Between these lines, $S_{0}$ is expected as the
hard-spin Ising model global minimizer. Thin lines show the contours
$E_{1}=E_{0}$ for $C=1,1.2,1.5,2,4$ in that region. (e) Success probability of
reaching $E_{0}$ (labelled as $SP_{0}$) and $E_{1}$ (labelled as $SP_{1}$)
states of the soft-spin energy in Eq. (1) from a point ${\bf x}$ with randomly
chosen components $x_{i}$ in $[-1,1]$ for different values of $p$ and
$J=0.4,N=8,C=1$. For larger values of $p$ a third state of higher energy
appears with the success probability $SP_{2}$; when projected on spins
$s_{i}=x_{i}/|x_{i}|$ this state corresponds to $S_{1}$. (f) The height of the
minimum energy barrier between $E_{1}$ and $E_{0}$ calculated as the energy
difference between $E_{1}$ and energy of the nearest saddle point is shown as
a black dashed line for $J=0.4,C=1$. The difference between $E_{0}$ and
$E_{1}$ is shown in red.

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 [32]; (see Supp. Inf). 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. (1)(f).

These 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 recognised before [33, 34], 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)=\epsilon(1-x_{i}^{2})$, where $\epsilon$ 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$

$x_{i}\rightarrow(1-\delta)x_{i}+\delta Rx_{i}/|x_{i}|.$ (3)

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 (2)e shows the probability of finding the ground state of the Ising
Hamiltonian using the HT networks [26], Eq.(2) (CIM-I), Eq.(2) with individual
pumping adjustments $p\rightarrow p_{i}$ according to
$\dot{p_{i}}=\epsilon(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(\epsilon\,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 utilising 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 the soft-spin nonlinear models.

We consider the transverse field Ising model with the circulant coupling
matrix $\mathbf{J}$ given by

$\displaystyle\hat{H}$
$\displaystyle=-\frac{1}{2}\sum_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{N,N}J_{ij}\hat{S}_{i}^{z}\hat{S}_{j}^{z}-\sum_{i=1}^{N}h_{i}\hat{S}_{i}^{z}-\gamma(t)\sum_{i=1}^{N}\hat{S}_{i}^{x}\,,$
(4) $\displaystyle\hat{S}_{i}^{\alpha}$
$\displaystyle=\mathds{1}\otimes\mathds{1}\otimes\cdots\mathds{1}\otimes\hat{S}^{\alpha}\underbrace{\otimes\mathds{1}\otimes\cdots\otimes\mathds{1}\otimes\mathds{1}}_{i-1\text{
terms }},\,\,\alpha=x,y,z\,,$ (5)

where $\hat{S}^{\alpha}$ are the spin-$1/2$ Pauli matrices, and $\mathds{1}$
is the $2\times 2$ identity matrix. The first term, $\hat{H}_{0}$, is diagonal
and represents the classical Ising Hamiltonian $H_{I}$. We will also consider
the effect of a symmetry-breaking longitudinal magnetic field given by the
second term, $\hat{H}_{1}$, although initially we will exclude this term by
setting $h_{i}=0$. The transverse field, $\hat{H}_{2}$, results in a non-
diagonal Hamiltonian operator and gives rise to the quantum Ising model. We
set the time-dependent annealing term to $\gamma(t)=B/\sqrt{t}$ for some
constant $B$ [35].

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
$|\Psi\rangle=\sum_{\xi}C_{\xi}|\xi\rangle$ with
$\sum_{\xi}\left|C_{\xi}\right|^{2}=1$, where the $C_{\xi}^{\prime}s$ are
complex numbers and the basis element
$|\xi\rangle\equiv\left|\xi_{1}\cdots\xi_{N}\right\rangle=\left|\xi_{N}\right\rangle\otimes\cdots\otimes\left|\xi_{1}\right\rangle,\quad\xi_{k}=\\{\lvert\downarrow\rangle,\,\,\lvert\uparrow\rangle\\}\,,\quad
k=1,\cdots,N.$

We begin with an initial state that is a ground state of the Hamiltonian
$\hat{H}_{2}$. The initial state at time $t_{i}$ can then be expressed as
$|\Psi(t_{i})\rangle=\lvert\psi_{\rightarrow}\rangle\otimes\cdots\otimes\lvert\psi_{\rightarrow}\rangle$
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 Supp. Inf. for details). As $t\rightarrow\infty$,
$\gamma(t)\rightarrow 0$ and the contribution of $\hat{H}_{2}$ decays to bring
about the target Hamiltonian. Provided $\gamma(t)$ is varied adiabatically, we
can expect that the system evolves, 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 Hamiltonian $\hat{H}_{0}$ given by
$P_{{}_{\rm GS}}=|\langle\phi_{{}_{\rm
GS}}\,|\,\mathopen{}\Psi(t)\rangle|^{2}$. For $J<J_{\rm crit}$, the $S_{0}$
ground state has a two-fold degeneracy. Therefore, we can expect that the
probability of finding one of the ground states is $P_{{}_{\rm GS}}=1/2$. In
Fig. 2f, we present numerical simulations of the time evolution of the success
probability for finding the ground state of $N=8$ spin system for $J=0.35$.
For comparison, we have also included the results for simulated annealing [36]
and classical annealing by evolving a master equation [35] (see Supp. Inf. for
details). In the former, we allow transition probabilities for single-spin
flips only, whereas in the latter, we allow for all spin flips. As can be seen
from Fig. 2f, only simulated annealing does not converge to the correct ground
states of the system. In contrast, quantum and the classical annealing recover
the correct ground states even though our value of $J$ lies within the range
$J_{e}<J<J_{\rm crit}$.

Figure 2: (a) Evolution of $N=8$ soft spins for $J=0.35$ and (b) $J=0.55$
according to Eq. (2). In each case, the ground state is recovered. The
amplitudes connected by the frustrated edges are lower than in the rest of the
system and are shown in red. In all runs, $C=1$, $p_{0}=J-2$,
$\epsilon=0.003$, $\Delta t=0.1$, and each $x_{i}(0)$ is chosen randomly from
a uniform distribution in the range $[-0.001,0.001]$; (c) Similar evolution
for probability amplitude of $\lvert\uparrow\rangle$ state in quantum
annealing simulation with $J=0.35$ and (d) $J=0.6$. The amplitudes colored in
red correspond to the frustrated spins as in (b); (e) Ground state probability
for HT, CIM-I, CIM-II, CIM-III, and quantum annealing (QA) for the Möbius
ladder graph with $N=8$. For CIM-III, for each value of $J$ the optimum value
of $0<\delta<1$ is chosen based on a set of preliminary runs in which $\delta$
is varied. Two thousand runs are used to calculate the probability of finding
the ground state $P_{\rm GS}$ for each value of $J$. For QA, $B=5$ and $\Delta
t=0.1$. Inset to (e) shows the same plots but for $N=100$; (f) Ground state
probability of target Hamiltonian with $J=0.35$ and $B=5$ for quantum
annealing (QA), single-spin simulated annealing (SA), classical annealing (CA)
and corresponding probabilities expected for adiabatic simulated (SA-ad) and
adiabatic quantum (QA-ad) annealing. Insets show Bloch vector for single-spin
and magnitude of Bloch vector $|{\mathbf{u}}|$; (g) similar results as (f) but
with symmetry-breaking terms added to Hamiltonian and Bloch sphere shows
typical trajectories of two neighbouring spins.

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 Supp. Inf. for details). In the
inset of Fig. (2)f, we show the evolution of the Bloch vector with time
evaluated for one of the spins (other spins show similar behaviour). 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 the emergence of Bell states: one of four maximally entangled
states. A definite state emerges only upon measurement, which then
subsequently collapses the corresponding wavefunction to one specific
configuration.

For direct 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.005\lvert\xi\rangle_{S_{1}}$
where $\lvert\xi\rangle_{S_{0}}$ and $\lvert\xi\rangle_{S_{1}}$ correspond to
the $S_{0}$ or $S_{1}$ state, respectively. This enforces the evolution
towards a specific ground state without the final state being in a quantum
superposition. The resulting Bloch vector for $J=0.35$ is shown in the inset
of Fig. 2g and now indicates evolution that ends at the surface of the Bloch
sphere, reaching either the $\lvert\uparrow\rangle$ or
$\lvert\downarrow\rangle$ state. This evolution 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 in the $J_{\rm e}<J<J_{\rm crit}$ range; (see Fig. 2e)
and demonstrates that correlations play a key role in facilitating the system
to converge to the true ground state. Although, its performance appears to
degrade near $J_{\rm crit}$, this is simply a consequence of keeping $B$
fixed. In general $B$ will need to be adjusted as we approach $J_{\rm crit}$
to account for the vanishing energy gap of the spectrum of the coupling matrix
${\mathbf{J}}$.

The corresponding time-dependent probability of finding each spin, $i$, in the
$\lvert\uparrow\rangle$ state is presented in Fig. 2c for $J=0.35$ (and for
$J=0.6$ in Fig. 2d). We see that even though we have strongly perturbed the
system towards the $S_{1}$ state for $J=0.35$, 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 summary, we analysed the optimization of Ising Hamiltonians, contrasting
the classical dynamics of semi-classical soft-spin models with quantum
annealing. We discussed semi-classical models’ challenges 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 optimising 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.

##  I Acknowledgements

J.S.C. acknowledges the PhD support from the EPSRC, N.G.B acknowledges the
support from the Julian Schwinger Foundation Grant No. JSF-19-02-0005, HORIZON
EUROPE/UKRI CL4-2021-DIGITAL-02-16, and Weizmann-UK Make Connection grant
142568.

##  II Supplementary Information

##  III Circulant graphs in physical systems

An all-optical scalable ONN was recently proposed for cyclic graphs; spatial
light modulators (SLMs) are used to discretize the optical field where each
pixel defines a different pulse amplitude [37]. 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
[38]). The condensates arranged in a circle interact with the nearest
neighbours, but the interactions beyond this decay exponentially fast.
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 synthesised [39]. 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 [40].

##  IV Eigenvectors and Eigenvalues of the Möbius Ladder coupling matrix

To find the eigenvalues and eigenvectors of $N\times N$ matrix $\mathbf{J}$,
for even $n$ we use the roots of unity, so the solutions of $\omega^{N}=1$
that are $\omega_{k}=\exp(i2\pi k/N)$ for $k=0,...,N-1$. The corresponding
eigenvectors are $(1,\omega_{k},\omega_{k}^{2},...,\omega_{k}^{N-1})$ which
can be verified by direct computation. Then, from the first row of
$\mathbf{J}$ we form the polynomial
$f(\omega)=-\omega-J\omega^{N/2}-\omega^{N-1}$ and evaluate it at the unit
roots $\omega_{k}=\exp[i2\pi k/N]$ to obtain the eigenvalues
$\lambda_{k}=f(\omega_{k})=-2\cos(2\pi k/N)-J(-1)^{k}$ with the corresponding
eigenvectors $v_{k}=(1,\omega_{k},\omega_{k}^{2},...,\omega_{k}^{N-1}).$ The
largest of $f(\omega_{k})$ is either $\lambda_{N/2}=2-J$ or $\lambda_{N/2\pm
1}=J+2\cos(2\pi/N)$ depending on whether $J<J_{\rm e}$ or $J>J_{\rm e}$ with
$J_{\rm e}\equiv=1-\cos(2\pi/N)$ is the value where these two eigenvalues
cross. The corresponding real-valued, and mutually orthogonal eigenvectors
$\mu_{k}$ can be formed from $v(\omega_{k})$ as
$\mu_{k}=Re[v(\omega_{k})]+Im[v(\omega_{k})]$ [30]. For two largest
eigenvalues the eigenvectors are $\mu_{N/2}=(1,-1,1,-1,...,-1)$ and
$\mu_{N/2\pm 1}=(1,\pm\cos(2\pi/N)\pm\sin(2\pi/N),...,\pm\cos(2\pi
k/N)\pm\sin(2\pi k/N),...,\pm\cos(2\pi(N-1)/N)\pm\sin(2\pi(N-1)/N))$. If $N/2$
is even, then $\mu_{N/2\pm 1}$ have the components with two zero values at the
positions separated by $N/2-1$ sign alternating components.

We illustrate this construction for the Möbius ladder coupling matrix
$\mathbf{J}$ with $N=8$ considered in the main text

${\bf J}=\left(\begin{array}[]{cccccccc}0&-1&0&0&-J&0&0&-1\\\
-1&0&-1&0&0&-J&0&0\\\ 0&-1&0&-1&0&0&-J&0\\\ 0&0&-1&0&-1&0&0&-J\\\
-J&0&0&-1&0&-1&0&0\\\ 0&-J&0&0&-1&0&-1&0\\\ 0&0&-J&0&0&-1&0&-1\\\
-1&0&0&-J&0&0&-1&0\\\ \end{array}\right).$ (6)

The eigenvalues are $\lambda_{0}=f(\omega_{0})=f(1)=-2-J,$
$\lambda_{1}=f(\omega_{1})=-\sqrt{2}+J,$ $\lambda_{2}=f(\omega_{2})=-J,$
$\lambda_{3}=f(\omega_{3})=\sqrt{2}+J,$ $\lambda_{4}=f(\omega_{4})=2-J,$
$\lambda_{5}=f(\omega_{5})=\sqrt{2}+J,$ $\lambda_{6}=f(\omega_{6})=-J,$ and
$\lambda_{7}=f(\omega_{7})=-\sqrt{2}+J.$

Figure 3: Eigenvalues of the Möbius ladder graph for $N=8$ as functions of
$J$. $J_{\rm e}$ is the value of $J$ where the red dashed line shows the two
largest eigenvalues cross. $J_{\rm crit}$ shows where their energies equal
$E_{1}=E_{0}$. The ground state corresponds to $S_{0}$ for $J<J_{\rm crit}$
and $S_{1}$ for $J>J_{\rm crit}$.

The eigenvector that corresponds to $\lambda_{3}$ is
$\mu_{3}=(1,-1,1,-1,1,-1,1,-1)$ and the eigenvector that corresponds to, say,
$\lambda_{4}$ is $\mu_{4}=(1,-\sqrt{2},1,0,-1,\sqrt{2},-1,0)$. The soft spin
system, therefore, follows $(+,-,+,+,-,+,-,-)$ or $(+,-,+,-,-,+,-,+)$
direction at the onset of the pitch-fork bifurcation when $J>J_{\rm e}$ while
$\lambda_{4}$ becomes the dominant eigenvalue of matrix $\mathbf{J}$. Figure 3
illustrates how these eigenvalues vary as a function of $J$.

##  V Basins of attraction

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 the 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) of
main text) 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.

Figure 4: Basins of attraction of the soft spin energy in Eq. (2) as defined
in the main text. We take $N=8$, $J=0.4$, $C=1$, various laser powers $p$,
20,000 randomly distributed $x_{i}$ in $[-1,1]$ and show which minimum is
reached via gradient descent using Newton’s method. To characterise points,
the average magnetization $m=\sum_{i}x_{i}/N$ (vertical axis) and the
correlations along the circle between $x_{i}$ defined as $X_{\rm
corr}=\sum_{i}(x_{i}-m)(x_{i+1}-m)/\sum_{i}(x_{i}-m)^{2}$ (horizontal axis),
are used. For small $p$, the basin of attraction is dominated by $S_{1}$ state
as any initial state descents to $E_{1}$. As $p$ grows, the ratio of the
volume of the basins of attraction of the $S_{1}$ state to the volume of the
basins of attraction of the $S_{0}$ approaches $4$. At the critical value of
$p\approx-0.08715$, both $S_{1}$ and $S_{0}$ states have the same energy, and
after that $S_{0}$ state becomes the ground state: this is indicated by the
switch between ground (blue) and excited (red) states.

##  VI Critical points

The main text gives the energy of the CIM ${E}$ in Eq. (1). We can determine
the critial (minima and the saddle) points by setting ${\partial E}/\partial
x_{i}=0$ for all $i=1,...,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.

Figure 5: Critical points of the CIM energy (defined by Eq. (1) of the main
text) for $N=8,C=1$ and different values of laser power $p$. Minima are shown
in green, saddles with one, two, three and 4+ unstable directions are shown as
red, light blue, blue and black points, respectively. The $S_{0}$ state is
furthest from the origin, which becomes the global minimum for
$p>p_{c}=-0.0872$. The minima for $p=2$ are $S_{0}$, $S_{1}$,
$(-,-,+,-,+,-,-,+)$, $(+,-,-,+,+,-,+,-)$ and $(+,-,-,+,+,-,-,+)$ in increasing
energy.

##  VII Solution of the Time-Dependent Schrödinger Equation

The wavefunction is evolved according to the time-dependent Schrödinger
equation (with $\hbar=1$) given by

$\displaystyle i\frac{\mbox{d}}{\mbox{d}t}\,\lvert\Psi(t)\rangle$
$\displaystyle=\hat{H}(t)\,\lvert\Psi(t)\rangle\,,$ (7)
$\displaystyle\hat{{H}}(t_{i})\,|\Psi(t_{i})\rangle$
$\displaystyle=\epsilon_{gs}\,|\Psi(t_{i})\rangle\,,$ (8)

where $\epsilon_{gs}$ denotes the ground state energy of the system at the
initial time.

To evolve the time-dependent Hamiltonian given by Eq.(7), we use a second
order accurate Strang time-splitting method where $\hat{H}_{0}$ is evolved for
half a time-step $\Delta t$ followed by $\hat{H}_{1}$ for a full time-step and
then $\hat{H}_{0}$ for another half a time-step. The resulting time-
integration scheme can then be written as

$\displaystyle\lvert\Psi(t_{n+1})\rangle=$
$\displaystyle\exp\left(-i\frac{\Delta
t}{2}\hat{H}_{D}\right)\left(-\frac{i}{2}\int_{t_{n}}^{t_{n+1}}\hat{H}_{2}(\tilde{t})\mbox{d}\tilde{t}\right)$
$\displaystyle\times$ $\displaystyle\exp\left(-i\frac{\Delta
t}{2}\hat{H}_{D}\right)\lvert\Psi(t_{n})\rangle\,,$ (9)

where $\hat{H}_{D}=\hat{H}_{0}+\hat{H}_{1}$ is the diagonal part of the
Hamiltonian operator.

By placing the Hamiltonian operator $\hat{H}_{2}(t)$ containing the time-
dependent term in the middle of the split-step algorithm, we ensure that we
have a symmetric time-splitting method. The time integral appearing in Eq. (9)
was evaluated analytically. In our simulations, we set $\Delta t\equiv
t_{n+1}-t_{n}=0.1$. The simulations were performed in MATLAB. The exponentials
of the diagonal and non-diagonal Hamiltonian can then be readily evaluated
using the `expm` function [41].

##  VIII Computation of Bloch Vectors in Quantum Annealing Simulations

The single-spin reduced density matrix $\hat{\rho}_{1,k}$ is obtained by
taking the partial trace of the $2^{N}\times 2^{N}$ density matrix
$\hat{\rho}$ over the Hilbert space of the other $N-1$ spins. For the $k$’th
spin, this is defined as

$\displaystyle\hat{\rho}_{1,k}(t)=\text{Tr}_{\\{N\backslash
k\\}}\hat{\rho}\,,$ (10)

where $\\{N\backslash k\\}$ denotes the $N$ spin Hilbert space excluding the
$k$’th spin. The single-spin density matrix can then be parameterised as

$\displaystyle\hat{\rho}_{1,k}$
$\displaystyle=\frac{1}{2}\left(\mathds{1}+\mathbf{u}_{k}\cdot\hat{\mathbf{S}}\right)=\frac{1}{2}\left(\begin{array}[]{cc}1+w_{k}&u_{k}-iv_{k}\\\
u_{k}+iv_{k}&1-w_{k}\end{array}\right)\,,$ (13)

where ${\mathbf{u}}_{k}=(u_{k},v_{k},w_{k})$ defines the corresponding Bloch
vector and $\hat{\mathbf{S}}=(\hat{S}^{x},\hat{S}^{y},\hat{S}^{z})$ correspond
to the vector of spin-1/2 Pauli matrices. For pure states the single-spin
reduced density matrix has rank 1, with a magnitude of the Bloch vector
$|\mathbf{u}_{k}|=1$. The surface of the Bloch sphere, therefore, represents
all the possible pure states whereas the interior of the sphere corresponds to
mixed states. The collapse of the Bloch vector towards the origin implies that
the state represents a maximally entangled Bell-like state.

##  IX Master Equation for Classical/Simulated Annealing

To model simulated annealing, we follow the method described in [35] and
introduce the Master equation for the transition probability $P_{i}(t)$ for
each spin configuration.

$\displaystyle\frac{\mbox{d}P_{i}(t)}{\mbox{d}t}=\sum_{j=1}^{2^{N}}A_{ij}(t)P_{j}(t)$
(14)

The $2^{N}\times 2^{N}$ matrix $A_{ij}(t)$ describes the transition rates. The
master equation can be written in the form

$\displaystyle\frac{\mbox{d}P_{i}(t)}{\mbox{d}t}$ $\displaystyle=\sum_{i\neq
j}A_{ij}(t)P_{j}(t)+A_{ii}(t)P_{i}(t)$ (15) $\displaystyle=\sum_{i\neq
j}\left(A_{ij}(t)P_{j}(t)-A_{ji}(t)P_{i}(t)\right)$ (16)

where we have made use of the conservation of probability given by

$\displaystyle\frac{\mbox{d}}{\mbox{d}t}\sum_{j}P_{j}(t)=\sum_{i,j}(A_{ji}P_{i})=0\,,$
(17)

to arrive at the final equality. Since the normalisation condition must hold
for any probabilities $P_{j}$, it follows that

$\displaystyle\sum_{j}A_{ji}=0\,,\qquad\text{or}\qquad A_{ii}=-\sum_{i\neq
j}A_{ji}\,.$ (18)

Using Eq.(18) to represent the diagonal terms of the Master equation ensures
that a numerical integration of this equation continues to conserve the
normalisation of the probabilities.

The precise form of the transition probabilities is problem specific although
it is common to use the Boltzmann distribution. In our work, we follow [35]
and use the Bose-Einstein distribution such that

$\displaystyle
A_{ij}(t)=\begin{cases}\left\\{1+\exp\left[\dfrac{\left(E_{i}-E_{j}\right)}{T(t)}\right]\right\\}^{-1}\,\,,&(\text{
single-spin flip) }\\\ -\sum_{k\neq i}A_{ki}\,\,,&(i=j)\\\
0\,,&\text{(otherwise) }\,.\end{cases}$

The form given above that is used for our simulated annealing simulations
means that entries of $A_{ij}$ are non-zero only for transitions corresponding
to single-spin flips. The annealing is performed by varying the temperature
$T(t)$ with time. To maintain consistency with our quantum annealing
simulations, we have varied the temperature according to $T(t)=D/\sqrt{t}$,
where $D$ is a free parameter which we set to $D=5$.

For our classical annealing (CA) simulations, we did not zero out any of the
transition probabilities in order to infer how collective transitions of spins
at each time-step, as opposed to only single-spin transitions, affects the
performance of classical algorithms. 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. It is also useful to make the observation that
quantum and classical annealing can be closely related to one another if one
formulates quantum annealing in imaginary time following a Wick rotation. It
then follows that the evolution of the $N$-spin wavefunction $|\Psi(t)\rangle$
is given by

$\displaystyle\frac{\mbox{d}}{\mbox{d}t}|\Psi(t)\rangle=\left(\mu(t)-\hat{{H}}(t)\right)|\Psi(t)\rangle\,\,.$
(19)

Here, $\mu(t)$ plays the role of a Lagrange multiplier which ensures that the
normalisation of the wavefunction is conserved in analogy with the
modification introduced above to the diagonal term of the Master equation.
Therefore, by comparing quantum annealing, simulated annealing, and classical
annealing, we can distinguish between the effects of retaining all-spin
transitions from the difference of evolving our equations in real and
imaginary time.

##  X Results of Quantum Annealing Simulations

Figure 6: Time-evolution of ground state probability of target Hamiltonian
with $J=0$ and $B=5$ for quantum annealing (QA), single-spin simulated
annealing (SA), classical annealing (CA) and corresponding probabilities
expected for adiabatic simulated (SA-ad) and adiabatic quantum (QA-ad)
annealing. (a) and (b) correspond to simulation without symmetry breaking term
which corresponds to a doubly degenerate ground state. Each figure corresponds
to the projection of the probability density onto each one of the ground
states; (c) and (d) correspond to a simulation with a symmetry breaking term
added which lifts the degeneracy and leads to a unique ground state.

Here, we include additional simulation results to complement the results
presented in the main text. In Figs. 2f and 2g of the main text, we presented
results for the time evolution of the ground state probability for $J=0.35$
and $B=5$. In Fig. 4, we also present results for $J=0$ for the probability
amplitude of the eigenstates with the two energy minima to discern the role of
grounstate degeneracy and how this is lifted by the symmetry breaking terms.
In Figs. 6a and 6b we 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, all
models including 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. 6c and 6d we present simulations for the case where a symmetry
breaking term is included as in the main text. Here, we used
$0.05\lvert\xi\rangle_{S_{0}}+0.05\lvert\xi\rangle_{S_{1}}$ to introduce the
symmetry breaking term. In contrast to the previous case, the behaviour 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.

Figure 7: Time-evolution of probability amplitude of $\lvert\uparrow\rangle$
state (left panels) and magnitude of Bloch vectors (right panels) in quantum
annealing computation with $J=0.35$ and $B=5$. The dashed lines correspond to
spins that at the end of the annealing align along the
$\lvert\downarrow\rangle$ state whereas solid lines correspond to spins that
align with the $\lvert\uparrow\rangle$ state. The red and blue colored lines
correspond to the color of the spins shown in Fig. 1b of main text.

These results demonstrate the importance of the symmetry breaking terms and
how they affect the time-evolution of the ground state probabilities. In order
to perform a more systematic study of the impact of including the symmetry
breaking terms, in Fig. 7 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 magntude 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 pure 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. 6c and 6d. 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.

## References

  * Kasprzak  _et al._ [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. Keeling, F. Marchetti, M. Szymańska, R. André, J. Staehli, _et al._ , Nature 443, 409 (2006).
  * Yamamoto  _et al._ [2017] Y. Yamamoto, K. Aihara, T. Leleu, K.-i. Kawarabayashi, S. Kako, M. Fejer, K. Inoue, and H. Takesue, npj Quantum Information 3, 1 (2017).
  * Inagaki  _et al._ [2016] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, _et al._ , Science 354, 603 (2016).
  * McMahon  _et al._ [2016] P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, _et al._ , Science 354, 614 (2016).
  * Honjo  _et al._ [2021] T. Honjo, T. Sonobe, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Kazama, K. Enbutsu, T. Umeki, R. Kasahara, K. ichi Kawarabayashi, and H. Takesue, Science Advances 7, eabh0952 (2021), https://www.science.org/doi/pdf/10.1126/sciadv.abh0952 .
  * Pal  _et al._ [2020] V. Pal, S. Mahler, C. Tradonsky, A. A. Friesem, and N. Davidson, Physical Review Research 2, 033008 (2020).
  * Babaeian  _et al._ [2019] M. Babaeian, D. T. Nguyen, V. Demir, M. Akbulut, P.-A. Blanche, Y. Kaneda, S. Guha, M. A. Neifeld, and N. Peyghambarian, Nature communications 10, 1 (2019).
  * Parto  _et al._ [2020] M. Parto, W. Hayenga, A. Marandi, D. N. Christodoulides, and M. Khajavikhan, Nature materials 19, 725 (2020).
  * Pierangeli  _et al._ [2019] D. Pierangeli, G. Marcucci, and C. Conti, Physical review letters 122, 213902 (2019).
  * 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 (2017).
  * Kalinin  _et al._ [2020] K. P. Kalinin, A. Amo, J. Bloch, and N. G. Berloff, Nanophotonics 9, 4127 (2020).
  * Vretenar  _et al._ [2021] M. Vretenar, B. Kassenberg, S. Bissesar, C. Toebes, and J. Klaers, Physical Review Research 3, 023167 (2021).
  * Kalinin  _et al._ [2023] K. Kalinin, G. Mourgias-Alexandris, H. Ballani, N. G. Berloff, J. H. Clegg, D. Cletheroe, C. Gkantsidis, I. Haller, V. Lyutsarev, F. Parmigiani, L. Pickup, _et al._ , arXiv preprint arXiv:2304.12594 (2023).
  * Goto  _et al._ [2021] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, and K. Tatsumura, Science Advances 7, eabe7953 (2021).
  * Syed and Berloff [2022] M. Syed and N. G. Berloff, “Physics-enhanced bifurcation optimisers: All you need is a canonical complex network,” (2022).
  * Calvanese Strinati and Conti [2022] M. Calvanese Strinati and C. Conti, Nature Communications 13, 7248 (2022).
  * Liu  _et al._ [2015] C.-W. Liu, A. Polkovnikov, and A. W. Sandvik, Phys. Rev. Lett. 114, 147203 (2015).
  * King  _et al._ [2022] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Nature Physics 18, 1324 (2022).
  * Heim  _et al._ [2015] B. Heim, T. F. Rønnow, S. V. Isakov, and M. Troyer, Science 348, 215 (2015), https://www.science.org/doi/pdf/10.1126/science.aaa4170 .
  * King  _et al._ [2023] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, Nature 617, 61 (2023).
  * Albash and Lidar [2018a] T. Albash and D. A. Lidar, Phys. Rev. X 8, 031016 (2018a).
  * Yan and Sinitsyn [2022] B. Yan and N. A. Sinitsyn, Nature Communications 13, 2212 (2022).
  * Albash and Lidar [2018b] T. Albash and D. A. Lidar, Phys. Rev. X 8, 031016 (2018b).
  * Hamerly  _et al._ [2019] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, K. Enbutsu, T. Umeki, R. Kasahara, S. Utsunomiya, S. Kako, K. ichi Kawarabayashi, R. L. Byer, M. M. Fejer, H. Mabuchi, D. Englund, E. Rieffel, H. Takesue, and Y. Yamamoto, Science Advances 5, eaau0823 (2019), https://www.science.org/doi/pdf/10.1126/sciadv.aau0823 .
  * Marsh  _et al._ [2023] B. P. Marsh, R. M. Kroeze, S. Ganguli, S. Gopalakrishnan, J. Keeling, and B. L. Lev, “Entanglement and replica symmetry breaking in a driven-dissipative quantum spin glass,” (2023), arXiv:2307.10176 [quant-ph] .
  * Hopfield [1982] J. J. Hopfield, Proceedings of the national academy of sciences 79, 2554 (1982).
  * Hopfield and Tank [1985] J. J. Hopfield and D. W. Tank, Biological cybernetics 52, 141 (1985).
  * Kalinin and Berloff [2018a] K. P. Kalinin and N. G. Berloff, arXiv preprint arXiv:1806.11457 (2018a).
  * Kalman and White [2001] D. Kalman and J. E. White, The American Mathematical Monthly 108, 821 (2001).
  * Zhang and Yang [2007] H. Zhang and Y. Yang, International journal of quantum chemistry 107, 330 (2007).
  * Gancio and Rubido [2022] J. Gancio and N. Rubido, Chaos, Solitons & Fractals 158, 112001 (2022).
  * Yamamura  _et al._ [2023] A. Yamamura, H. Mabuchi, and S. Ganguli, “Geometric landscape annealing as an optimization principle underlying the coherent ising machine,” (2023), arXiv:2309.08119 [cond-mat.dis-nn] .
  * Kalinin and Berloff [2018b] K. P. Kalinin and N. G. Berloff, New Journal of Physics 20, 113023 (2018b).
  * Leleu  _et al._ [2020] T. Leleu, F. Khoyratee, T. Levi, R. Hamerly, T. Kohno, and K. Aihara, arXiv e-prints , arXiv (2020).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998).
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983), https://www.science.org/doi/pdf/10.1126/science.220.4598.671 .
  * Calvanese Strinati  _et al._ [2021] M. Calvanese Strinati, D. Pierangeli, and C. Conti, Phys. Rev. Appl. 16, 054022 (2021).
  * Ohadi  _et al._ [2016] H. Ohadi, R. Gregory, T. Freegarde, Y. Rubo, A. Kavokin, N. Berloff, and P. Lagoudakis, Physical Review X 6, 031032 (2016).
  * Ayoub and Psaltis [2021] A. B. Ayoub and D. Psaltis, Scientific Reports 11, 1 (2021).
  * Qiang  _et al._ [2016] X. Qiang, T. Loke, A. Montanaro, K. Aungskunsiri, X. Zhou, J. L. O’Brien, J. B. Wang, and J. C. Matthews, Nature communications 7, 1 (2016).
  * Norambuena  _et al._ [2020] A. Norambuena, D. Tancara, and R. Coto, European Journal of Physics 41, 045404 (2020).