# Ising Machines for Diophantine Problems in Physics

Steven A. Abel  0000-0003-1213-907X<EMAIL_ADDRESS>Luca A. Nutricati
0000-0002-5045-5113<EMAIL_ADDRESS>Institute for Particle
Physics Phenomenology, Durham University, Durham DH1 3LE, UK  Department of
Mathematical Sciences, Durham University, Durham DH1 3LE, UK

###### Abstract

Diophantine problems arise frequently in physics, in for example anomaly
cancellation conditions, string consistency conditions and so forth. We
present methods to solve such problems to high order on annealers that are
based on the quadratic Ising Model. This is the intrinsic framework for both
quantum annealing and for common forms of classical simulated annealing. We
demonstrate the method on so-called Taxicab numbers (discovering some
apparently new ones), and on the realistic problem of anomaly cancellation in
$U(1)$ extensions of the Standard Model.

##  I Introduction

As well as being of intrinsic interest in number theory, Diophantine problems
play an important role in fundamental physics. The purpose of this paper is to
develop methods for solving the kinds of Diophantine problems that frequently
occur in particle physics, where for example they appear in anomaly
cancellation conditions as systems of cubic equations. They also appear in
index theorems, in consistency conditions in string theory, in problems
relevant for computing the effective potential in string theory and for
finding vacua with small cosmological constant, as well as in powerful
nonperturbative methods that are used in field theory, such as ’t Hooft
anomaly matching tHooft:1979rat , and a multitude of other applications.

Unfortunately Diophantine problems are also notoriously difficult to solve.
Sometimes they can be significantly simplified using for example Gröbner basis
methods, but typically such a problem is computationally hard. This is
certainly the case for the typical anomaly cancellation problem, which entails
the solving of a coupled set of cubic equations with an independent rational
variable appearing for every charge of every particle. In a system whose size
is comparable to that of say the Standard Model of particle physics, an
exhaustive scan becomes infeasible even when the domain of allowed charges is
restricted. (See for example Refs. Allanach:2018vjg ; Allanach:2019uuu ;
Allanach:2020zna ; Allanach:2021bfe .) Indeed determining precisely which
complexity class a problem falls into is itself an important question. (See
for example Ref. TUNG1987324 , and for discussions in the string theory
context Refs. Denef:2006ad ; Denef:2017cxt ; Halverson:2018cio ;
Halverson:2019vmd .)

For such problems, scaleable Ising hardware solvers, which form the basis of
both simulated and quantum annealers (as introduced in Ref. PhysRevE.58.5355 ;
doi:10.1126/science.1057726 ), could have huge impact, particularly as NP
problems can be formulated as Ising problems with only polynomial overhead
RevModPhys.80.1061 ; 10.3389/fphy.2014.00005 . (For some practical
applications of Ising machines and for recent reviews see Refs. Hauke_2020 ;
Mohseni ).

This paper will show how to encode Diophantine problems of the kind described
above onto such machines. We shall do this by focusing on two Diophantine
tasks.

The first is the purely number theoretic one of finding what we will refer to
generally as “Taxicab” numbers, namely those numbers that can be expressed in
more than one way as sums of equal powers. The most famous example is the
number of Hardy and Ramanujan’s eponymous taxi, ${\rm Ta}(2)=1729$. This is
the smallest of the following list of numbers, all of which are expressible as
the sum of two cubes in two different ways:

$\displaystyle 1729\leavevmode\nobreak\ =\leavevmode\nobreak\
1^{3}+12^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 9^{3}+10^{3},$
$\displaystyle 4104\leavevmode\nobreak\ =\leavevmode\nobreak\
9^{3}+15^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 16^{3}+2^{3},$
$\displaystyle 20683\leavevmode\nobreak\ =\leavevmode\nobreak\
24^{3}+19^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 10^{3}+27^{3},$
$\displaystyle 32832\leavevmode\nobreak\ =\leavevmode\nobreak\
32^{3}+4^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 18^{3}+30^{3},$
$\displaystyle...$ (1)

We shall use the notation $(k,m,n)$, to refer to such numbers, where $k$ is
the power, while $m$ and $n$ are the number of terms on each side. Thus
Ta$(2)\leavevmode\nobreak\ =\leavevmode\nobreak\ 1729$ is defined to be the
smallest $(3,2,2)$ number, while Fermat’s theorem is the statement that
$(k,1,2)$ numbers only exist for $k=2$. Here we will develop annealing methods
to determine the above list of $(3,2,2)$ numbers (where we consider all
numbers in the list to be of interest not just the smallest). We also test our
methods on several variations, namely $(4,3,3)$, $(3,1,5)$, $(3,1,7)$,
$(3,6,6)$, $(3,7,7)$, $(3,8,8)$. Examples of most of these are known, and can
be found in Refs. Guy:1994a ; Meyrignac ; Eulernet ; Weisstein ; piezas ,
although some were discovered only with the advent of high performance
computing and appeared relatively recently. However some, such as $(3,7,7)$
and $(3,8,8)$ numbers, do not seem to have been known before. (Indeed as we
shall see the latter represent solutions in a search space of size $\sim
10^{24}$.)

The second task that we will consider is the physical one of finding solutions
to the anomaly cancellation conditions of a typical extension of the Standard
Model of particle physics. In four space-time dimensions this is as mentioned
also a cubic (i.e. homogeneous, third order) Diophantine problem, in which the
integers correspond to the numerators of rational gauge charges. (In $2d$
dimensions the equations are instead order $d+1$). In this case as well, we
will find that a sufficiently well-crafted encoding onto an annealer allows
one to solve for anomaly cancellation in the systems considered in Refs.
Allanach:2019uuu ; Allanach:2020zna ; Allanach:2021bfe , at an already
relatively advanced level, and in a very short time, and certainly without
having to perform any kind of exhaustive scan.

Although we will not consider ’t Hooft anomaly matching explicitly in our
discussion, it is worth mentioning that this nonperturbative procedure is
morally a generalisation of the Ta$(2)$ Taxicab problem. That is, when two
particle theories are duals of each other, then by ’t Hooft’s argument they
must realise the same set of global anomalies in two different ways. Similarly
physical trialities and quadralities are generalisations of the Ta$(3)$ and
Ta$(4)$ Taxicab problems (which realise the same number in 3 and 4 different
ways respectively). Such systems exist but are comparatively rare
Gadde:2013lxa ; Franco:2016nwv ; Franco:2016tcm .

Constructing suitably efficient encodings for these sorts of problems requires
significant advancement. Firstly as we shall review we are interested in the
kinds of set-ups found in quantum annealers, in which problems are encoded in
the Hamiltonian of a quadratic Ising model, which must then be minimised to
solve the problem. Thus the crux of the matter in encoding a non-trivial
system of cubic- and higher-order Diophantine equations is to implement a
reduction procedure that can represent the complete system as a single loss-
function represented by a spin-Hamiltonian that is at most quadratic. Of
course Diophantine problems, in particular factorisation, have been considered
on Ising model annealers before dattani2014quantum ; tanburn:2015a ;
Jiang:2018b ; peng2019factoring ; warren2019factoring ; wang2020prime , along
with quadratic systems of polynomial equations chang2019quantum ;
chang2019least . However both these problems can be mapped into the
optimisation of an order four Hamiltonian which can in turn be reduced to a
quadratic Ising model suitable for a quantum annealer, with only two layers of
reduction. By contrast, here we will be considering problems of order much
higher than two. To accomplish this, we use a procedure to automate the
reduction of an arbitrary order Hamiltonian to a quadratic one. This
procedure, which iterates that first appearing in Ref. rosenberg1975reduction
and then more recently in Refs. tanburn:2015a ; dattani:2019a ; Hauke_2020 ;
Abel:2022lqr , is completely problem-independent and therefore potentially
applicable to any set of Diophantine equations. It can perform the many layers
of reduction required to reach a quadratic spin-Hamiltonian representation of
the high order problems we will be considering.

Furthermore, the integers in the Diophantine equations are naturally encoded
using a binary representation, which is what ultimately will enable the
procedure to be much more efficient than a systematic scan. However this in
turn leads to large and highly connected Ising model encodings of the
Diophantine system of equations. Thus to improve efficiency we introduce a
technique we refer to as solution-mining. This innovation begins by finding
one solution in the conventional manner. A random perturbation from this
solution then serves as the new starting point for the next run. From there
the system often tunnels to a new solution nearby, and when it does so this
serves as the next starting point, and so on. This approach appears to be
effective for problems, like anomaly cancellation, that contain many coupled
Diophantine equations, and it allows the domain of solutions to wander in the
parameter space. We shall describe these techniques in Section II.

The rest of the paper is organised as follows. In Section III we begin with a
warm-up problem which is to solve a single Diophantine polynomial equation
with two variables. We use this simple problem to discuss how annealing is
implemented on a quantum annealer, and to compare the performance of currently
available quantum annealers with classical simulated annealers. Our conclusion
here is that quantum annealers are only just beginning to catch up with
classical simulated annealing in terms of the complexity of problems that can
be embedded and solved, but when they do become competitive the potential
speed-ups could be substantial for this kind of problem. Then Sections IV and
V consider the Taxicab problems, and anomalies respectively.

##  II Methods for encoding Diophantine problems

The general principle behind an Ising machine is to solve problems by
reformulating their solution as the minimisation of a function $H$ of spin
variables $\sigma_{\ell}\leavevmode\nobreak\ =\leavevmode\nobreak\ \pm 1$,
where $\ell$ labels the spin sites. On a quantum annealer this so-called
problem-Hamiltonian, $H$, is forced (by the physical architecture of the
annealer) to be quadratic in the spins,

$H(\sigma_{\ell})\leavevmode\nobreak\ =\leavevmode\nobreak\
\sum_{\ell}h_{\ell}\sigma_{\ell}+\sum_{\ell m}J_{\ell
m}\sigma_{\ell}\sigma_{m}\leavevmode\nobreak\ ,$ (2)

where the spins would of course correspond to physical qubits. We can postpone
further discussion of the physical realisation of the system (which will be
described later) and focus on the embedded abstract problem which will apply
to any annealer of this kind.

Specialising to the present case, we are interested in solving a set of
polynomial Diophantine equations $f_{A}(t_{i})=0$, where
$t_{i}\in{\mathbb{Z}}$ are the would-be integer solutions to the problem of
interest. These integers will be encoded on the annealer in a binary format,
namely we will use the following encoding:

$t_{i}\leavevmode\nobreak\ =\leavevmode\nobreak\
\tau_{i,0}+2\tau_{i,1}+\dots+2^{\beta-1}\tau_{i,\beta-1}+s_{i}\leavevmode\nobreak\
.$ (3)

We use $\tau$ to denote the binary variables corresponding to a given spin,

$\tau_{i,k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{1}{2}(1+\sigma_{i,k})\leavevmode\nobreak\ ,$ (4)

with $\tau_{i,k}\in\\{0,1\\}$, and where we allow classical integer shifts,
$s_{i}\in\mathbb{Z}$. These shifts can for example be negative to allow the
domain to include negative integers or, as will be the case for solution-
mining, they can be adjusted iteratively to explore the search space.

We would like to use such an encoding to solve the Diophantine equations, and
this can in principle be done by finding the minimum of a loss-function
Hamiltonian,

$H_{\rm D}\leavevmode\nobreak\ \equiv\leavevmode\nobreak\
\sum_{A}(f_{A}(t_{i}(\sigma_{\ell})))^{2}\leavevmode\nobreak\ .$ (5)

In addition one may wish to add several constraints: for example the Taxicab
numbers are usually defined as the smallest numbers expressible in two
different ways. Such conditions can be included with a constraint Hamiltonian,
$H_{C}$, which might in the case of the Taxicab numbers be simply the numbers
themselves (since they are positive). Thus we begin with an idealised (i.e.
non-quadratic) system,

$\tilde{H}(\sigma_{\ell})\leavevmode\nobreak\ =\leavevmode\nobreak\
H_{D}(t_{i}(\sigma_{\ell}))+H_{C}(t_{i}(\sigma_{\ell}))\leavevmode\nobreak\ .$
(6)

Note that solutions to the Diophantine equations all have $H_{D}=0$, so that
constraints imposed by $H_{C}$ would independently select the preferred
solution. However the converse is generally not true: that is one should avoid
over-weighting the constraints $H_{C}$ such that competing minima appear that
have lower $H_{C}$ but $H_{D}\neq 0$. Of course in many cases the desired
solutions are very rare, so it is often much more efficient (or more precisely
not an NP-hard problem) to simply apply any desired constraints by post-
processing the solutions (e.g. for the Taxicab numbers, one could simply
select by hand the smallest number).

###  II.1 Reduction

For a Diophantine system containing order-$d$ polynomials in $t_{i}$, the raw
Hamiltonian $\tilde{H}$ in Eq. (6) is an order-$2d$ polynomial in the spins,
$\sigma_{\ell}$. For example the ${\rm Ta}(2)$ problem yields an order-6 spin-
polynomial. Therefore we now arrive at the task of transforming $\tilde{H}$
into an equivalent quadratic problem-Hamiltonian, $H$, that is guaranteed to
have the same minima, but which can be implemented on the annealer.

For this task we shall use the reduction method described in the Appendix of
Ref. Abel:2022lqr . This method works by introducing auxilliary spins111We
think it is more accurate to use ‘auxilliary’ to refer to both abstract spins
and later to qubits, rather than the quantum computing term, ‘ancillary’. to
represent pairs of spins in the original Hamiltonian of Eq. (6), and is one of
the many methods in the comprehensive survey of Ref. dattani:2019a . (We
should remark that there exist “qubit-saving” reduction methods that do not
require auxilliary spins, but these are more task specific and currently
appear to be restricted to reduction of terms in the Hamiltonian with products
of 3 or 4 spins tanburn:2015a .)

The method works as follows. We begin with the raw polynomial
$\tilde{H}(\sigma_{\ell})$ written as a function of binary variables using Eq.
(4). Suppose $\tilde{H}$ has terms involving products of two binary variables
$\tau_{1}$ and $\tau_{2}$. Now consider adding to the polynomial $\tilde{H}$ a
quadratic term that involves the binary variables together with a new
auxiliary variable $\tau_{12}$, which is of the form

$Q(\mbox{\small$\tau_{12};\tau_{1},\tau_{2}$})\leavevmode\nobreak\
=\leavevmode\nobreak\
\Lambda(\tau_{1}\tau_{2}-2\tau_{12}(\tau_{1}+\tau_{2})+3\tau_{12})\leavevmode\nobreak\
.$ (7)

Inspection shows that a sufficiently large and positive overall coupling
$\Lambda$ enforces $\tau_{12}\leavevmode\nobreak\ =\leavevmode\nobreak\
\tau_{1}\tau_{2}$. Importantly the minimum at this point has $Q=0$. Therefore
we may replace the product $\tau_{1}\tau_{2}$ with $\tau_{12}$ wherever it
appears within $\tilde{H}$, and the new Hamiltonian is guaranteed to have the
same set of minima as the original $\tilde{H}$. Therefore the process can be
iterated until one arrives at the desired problem-Hamiltonian which is
quadratic in spins, and which is schematically of the form

$\displaystyle H_{D}+H_{C}\leavevmode\nobreak\ $
$\displaystyle=\leavevmode\nobreak\
\tilde{H}(\mbox{\small$\tau_{1},\tau_{2},\ldots,\tau_{12},\tau_{13},\ldots,\tau_{12,34},\tau_{12,56}\ldots$})$
$\displaystyle+\sum_{i>j}Q(\mbox{\small$\tau_{ij};\tau_{i},\tau_{j}$})+\sum_{i<j,k<m}Q(\mbox{\small$\tau_{ij,km};\tau_{ij},\tau_{km}$})$
$\displaystyle\qquad\leavevmode\nobreak\ +\ldots$ (8)

with the constraints imposed by the $Q$ terms ensuring that this quadratic
Hamiltonian has the same minima as the original order-$2d$ polynomial.

We can check that the reduction works correctly with the example order-3
Hamiltonian

$\displaystyle\tilde{H}\leavevmode\nobreak\ $
$\displaystyle=\leavevmode\nobreak\
\sigma_{1}\sigma_{2}\sigma_{3}\leavevmode\nobreak\ $
$\displaystyle\equiv\leavevmode\nobreak\
8\tau_{1}\tau_{2}\tau_{3}-4\tau_{1}\tau_{2}-4\tau_{1}\tau_{3}-4\tau_{2}\tau_{3}$
$\displaystyle\qquad\qquad\qquad\qquad+2\tau_{1}+2\tau_{2}+2\tau_{3}\leavevmode\nobreak\
,$ (9)

where we drop the constant $-1$ in translating to the binary variables. This
Hamiltonian has 4 minima at
$\sigma_{1}\sigma_{2}\sigma_{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ -1$
(which corresponds in binary language to any one of the $\tau_{\ell}$ being
zero, or all of them), as opposed to the seven solutions to
$\tau_{1}\tau_{2}\tau_{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0$. As
described above we can reduce the trilinear term by trading $\tau_{1}\tau_{2}$
for an auxiliary binary $\tau_{12}$ and adding the Hamiltonian
$Q(\mbox{\small$\tau_{12};\tau_{1},\tau_{2}$})$. The quadratic problem-
Hamiltonian (in QUBO language) is then

$\displaystyle H\leavevmode\nobreak\ $ $\displaystyle=\leavevmode\nobreak\
Q(\mbox{\small$\tau_{12};\tau_{1},\tau_{2}$})\leavevmode\nobreak\
+\leavevmode\nobreak\ 8\tau_{12}\tau_{3}$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \qquad-4\tau_{12}-4\tau_{1}\tau_{3}-4\tau_{2}\tau_{3}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\qquad\qquad+2\tau_{1}+2\tau_{2}+2\tau_{3}\leavevmode\nobreak\ .$ (10)

It is easy to verify that provided $\Lambda>2$ the original 4 degenerate
solutions hold in the new combined Hamiltonian as required.

With the increasing complexity of the raw Hamiltonian $\tilde{H}$ and a
limited number of physical qubits at our disposal, we of course aim to find a
reduction procedure that minimises the number of auxiliary variables.
Therefore, the central question is how can we choose the smallest set of spin
pairs that correctly collapses all the higher order terms to quadratics? In
the case of cubic to quadratic, finding a spin-optimised procedure is
equivalent to the set cover problem which can in turn be cast as 0-1 ILP
babbush2013resource . Both set cover and 0-1 ILP are well known to be NP-
complete by analogy with vertex cover boros2002pseudo . Therefore,
generalising this to arbitrary order Hamiltonians would recast our spin-
optimised problem into a task which is at least equivalent to solving $k-2$
NP-complete problems, where $k$ is the order of the Hamiltonian and therefore
$k-2$ are the required layers of reduction. For this reason, we shall use a
different approach based on a simple greedy algorithm which works as follows:
at each reduction stage it finds the pair of binary variables
$\tau_{i}\tau_{j}$ that appears most often in the Hamiltonian; wherever the
pair appears, we replace $\tau_{i}\tau_{j}$ with the auxilliary logical spin
$\tau_{ij}$, and add the penalty term in Eq. (7). The quadratised Hamiltonian
is constructed by repeating these three steps iteratively. In the language of
set covering, this is equivalent to the greedy heuristic algorithm first
proposed in Ref. chvatal1979greedy . In Fig. 1 we have collected three plots
which show how the average number of required auxilliary variables grows as we
increase the number of cubic interactions, the rate of this growth in the
linear central region and the time required to perform the reduction as a
function of the number of cubic couplings.

(a)

(b)

(c)

Figure 1: In (a) we plot the average number of auxilliary spins required to
quadratise the Hamiltonian, versus the square root of the number of cubic
interactions for different numbers of total spins $n$. These curves exhibit a
linear behavior in the central region where the slope is given by the square
root of the total number of spins, as we can see in (b). In (c) we plot the
time required to quadratise the Hamiltonian as a function of the number of
cubic couplings. We clearly see the time increasing linearly with the size of
the problem.

These results are very similar to those obtained in Ref. babbush2013resource
where the optimal reduction is found by solving exactly the equivalent 0-1
ILP. Two important remarks are in order. First, as we can see in Fig. 1a, both
methods saturate approximately when the Hamiltonian contains all possible
cubic interactions with $n$ qubits, namely when
$N_{\text{3-couplins}}\approx\binom{n}{3}$. Second, we see that in either case
the number of auxilliary spins increases linearly with the square root of the
number of cubic interactions and the growth rate is given by the square root
of the total number of spins, as shown in Fig. 1b. Nevertheless, as we can see
in Fig 1a, in some regions our procedure uses a number of auxiliaries that is
larger than the value at saturation, especially in the $n\leavevmode\nobreak\
=\leavevmode\nobreak\ 12$ and $n\leavevmode\nobreak\ =\leavevmode\nobreak\ 13$
cases. This is of course due to the fact that we are not seeking an exact
solution of the spin-optimised reduction problem. Indeed, a local optimal
choice of our reduction algorithm does not necessarily lead to a global
minimum in terms of the number of auxilliary spins. More precisely, in the
language of the equivalent set covering problem, it has been shown in Ref.
chvatal1979greedy  that this greedy algorithm returns to an approximate
solution which cannot be bigger than $H(n)$ times the minimum one, where
$H(n)$ is the $n$-th harmonic number and $n$ the size of the set to be covered
(namely in our case the set of all higher order couplings). However, in the
problems treated below, the greedy algorithm we adopt returns quadratised
Hamiltonians with at most $\sim 300$ logical spins, far below the limit
imposed by for example the number of available qubits in the currently
accessible quantum annealers, making it unnecessary to solve the problem
exactly. Finally, we should remark that this procedure is straightforwardly
generalisable to Hamiltonians of arbitrary order, requiring a number of steps
which grows roughly linearly with the size of the problem, as we can see in
Fig. 1c and also in Ref. chvatal1979greedy . As expected, this is in contrast
with the exact method discussed in Ref. babbush2013resource  which shows an
exponentially increasing amount of time with the increasing complexity of the
Hamiltonian.

Reduced Hamiltonians can be represented using connected graphs in which nodes
correspond to spins and links to couplings. As an example, Fig. 2 is a
representation of the quadratised Hamiltonian associated to the first
Diophantine equation in Table 2.

Figure 2: Representation of $H\leavevmode\nobreak\ =\leavevmode\nobreak\
(x_{1}^{2}+x_{2}^{2}-149)^{2}$ as a quadratised Ising model. Nodes,
corresponding to spins, can be arranged in a circle, and they are
quadratically coupled by the links. (Diagrams for non-quadratised models would
contain junctions in the couplings). Weaker couplings are represented in light
grey, gradually getting darker for higher coupling strengths. Linear couplings
are coloured from dark blue (large negative terms) to dark red (large positive
terms)

.

###  II.2 Solution mining for improved performance

Given the increase in difficulty with bitnumber, $\beta$, we will utilise a
method for improving the performance. This method allows one to explore larger
regions of parameter space (i.e. larger integers) without increasing $\beta$,
yielding in turn solutions with larger values.

The method operates iteratively, by at each run constructing a brand new
Hamiltonian from the previously found solutions. At say the $k$-th iteration,
we minimise the Hamiltonian looking for solutions of the form

$t_{i}^{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\tau_{i,0}^{(k)}+2\tau_{i,1}^{(k)}+s_{i}^{k}\leavevmode\nobreak\ ,$ (11)

where $k\leavevmode\nobreak\ =\leavevmode\nobreak\ 0,...,N$ (with $N$ being
the total number of anneal runs), and where $s^{k}_{i}$ is a classical shift
that centres the new search, which is determined from a solution found in the
$(k-1)$-th run: if we designate the previous solution $\hat{t}_{i}^{k-1}$,
then the $\\{s_{i}^{k}\\}$ are chosen such that

$\displaystyle t_{i}^{k}\leavevmode\nobreak\ \in\leavevmode\nobreak\
[\hat{t}_{i}^{k-1}-1\,,\,\hat{t}_{i}^{k-1}+2]\quad\text{or}$ $\displaystyle
t_{i}^{k}\leavevmode\nobreak\ \in\leavevmode\nobreak\
[\hat{t}_{i}^{k-1}-2\,,\,\hat{t}_{i}^{k-1}+1]\leavevmode\nobreak\ ,$ (12)

based on a random choice.

This procedure finds new solutions by performing a kind of “random tunnelling”
from previously found solutions (hence the name “solution-mining”). It
generally operates well when there are many variables in the system and many
different equations, because in such systems the solutions can be relatively
close in each dimension of the search space (even though the total Hamming
distance could be very large due to the large number of dimensions). For the
two specific example problems we are discussing here, it is not a useful
enhancement for finding Taxicab numbers because there one is seeking the
smallest numbers, and (as we shall see) the solutions to the Diophantine
system are very widely spaced. However for solving anomaly equations the
method is a significant improvement. In such systems, new solutions to the
anomaly equations tend to appear with consistent frequency when the allowed
charge size is increased, and it is the sheer number of anomaly equations and
charges that makes the problem difficult.

It should be noted that there is no additional cost for solution-mining
because even though a brand new Hamiltonian must be constructed at each stage,
the embedding graph remains the same if the values of $\beta$ do not change.
This means we construct an entirely new Hamiltonian $\tilde{H}$, but do not
need to perform a new reduction of the solution. On a quantum annealer we
perform reverse annealing (to be explained below) in order to collect the
solution and construct the new Hamiltonian at each stage, which then simply
has to be translated into new couplings via the updated $\\{s_{i}^{k}\\}$
values.

##  III Quantum versus classical annealers for Diophantine equations

In this section we shall compare quantum and classical annealers to solve
Diophantine equations, using the encoding methods described above on a simple
warm-up problem. It has been already shown in Ref. chang2019quantum  that a
quantum annealer can successfully be used to solve second order systems of
polynomial equations. On the other hand, simulated annealing and variations
thereof have been extensively applied for similar purposes (see Ref
Aguiar2014DiophantineEA ; moradi2019solving ) along with other techniques such
as genetic algorithms, particle swarm optimisation and so on (see for example
abraham2001diophantine ; abraham2010particle ; abraham2013finding ;
perez2013numerical ). In the following we shall see that quantum annealing is
also successful in solving equations of order higher than two, retrieving most
of the results found using Particle Swarm Optimisation in Ref.
abraham2010particle  and Fuzzy Adaptive Simulated Annealing in Ref.
Aguiar2014DiophantineEA . We shall find that simulated annealing on the
quadratised Hamiltonian can be used to solve problems which are still hard to
embed in a quantum annealer due to technological limitations (restricted
number of qubits and limited connectivity). Such problems will realistically
be solvable in the near future with annealers that have higher connectivity
but still with at most quadratic interactions.

Let us first outline the central features of the quantum annealer that we will
need for this study (for a comprehensive review see Ref. Hauke_2020 ). The
quantum annealer is defined by a Hamiltonian of physical qubits of the form

$\displaystyle\mathcal{H}(s)$ $\displaystyle\leavevmode\nobreak\
=\leavevmode\nobreak\
A(s)\sum_{\ell}\sigma_{\ell,x}+\,B(s)H(\sigma_{\ell,z})\,,$ (13)

where $H$ is the problem-Hamiltonian in Eq. 2, $\sigma_{\ell,x}$ and
$\sigma_{\ell,z}$ are the Pauli matrices acting on the $\ell^{\rm th}$ qubit,
and $A(s)$, $B(s)$ are smooth functions such that $A(1)\leavevmode\nobreak\
=\leavevmode\nobreak\ B(0)\leavevmode\nobreak\ =\leavevmode\nobreak\ 0$ and
$A(0)\leavevmode\nobreak\ =\leavevmode\nobreak\ B(1)\leavevmode\nobreak\
=\leavevmode\nobreak\ 1$. Meanwhile the internal couplings $h_{\ell}$ and
$J_{\ell m}$ are fixed in each anneal session.

Since $\mathcal{H}(1)\leavevmode\nobreak\ =\leavevmode\nobreak\
H(\sigma_{\ell})$, the usual quantum annealing strategy is adiabatically to
adjust the pre-factors $A(s)$ and $B(s)$ using the parameter $s$, such that
the system ends up in a global minimum of $H$. The time dependence of $s(t)$
is defined by the user in a so-called anneal-schedule. Thus during a reverse
anneal for example one begins at $s\leavevmode\nobreak\ =\leavevmode\nobreak\
1$ and with the system set in any eigenstate of
$\bigotimes_{\ell}\sigma_{\ell,z}$. Then one allows the system to evolve
quantum-mechanically by bringing $s$ to small values using a piecewise-linear
function $s(t)$ of time $t$, which completes back at $s(t_{\rm
final})\leavevmode\nobreak\ =\leavevmode\nobreak\ 1$ when the measurement of
final spins is made.

A technological limitation is that the connectivity of the quantum annealing
device in terms of the allowed non-vanishing couplings $J_{\ell m}$ between
qubits is limited. Let us be specific to the architecture we will be using in
this work, namely D-Wave’s LantingAQC2017  Advantage_system4.1: this annealer
contains 5627 qubits, connected in a _Pegasus_ structure, but only has a total
of 40279 couplings between them.

The warm-up problem that we will use to compare this kind of annealer with
classical simulated annealing is simply to find solutions of a generic
Diophantine equation

$f(x_{1},...,x_{N})\leavevmode\nobreak\ =\leavevmode\nobreak\ 0\,,$ (14)

where $x_{1},...,x_{N}\in\mathbb{Z}$ and
$f:\mathbb{Z}^{N}\rightarrow\mathbb{Z}$ is a generic order $k$ polynomial
function.

In order to solve such an equation, we again square it to form the problem-
Hamiltonian:

$H(x_{1},\dots,x_{N})\leavevmode\nobreak\ \equiv\leavevmode\nobreak\
[f(x_{1},\dots,x_{N})]^{2}\,.$ (15)

The next step is to binary encode $x_{1},\dots,x_{N}$ as in Eq. (3), choosing
the values of $s_{i}$ and $\beta$ depending on the specific problem.

The above Hamiltonian is therefore an order-$2k$ polynomial in $\tau_{i,j}$.
Nevertheless, regardless of how high the order is, it can always be reduced to
a quadratic Ising Hamiltonian making use of the procedure described in Sec.
II.1. Using the D-Wave’s quantum annealer we find all the solutions listed in
Table 2 of Ref. abraham2010particle  (which have also been found using Fuzzy
Adaptive Simulated Annealing in Ref. Aguiar2014DiophantineEA ). We report some
of them in Table 1.

Equation | Solution | $N_{\text{aux}}$ | $\beta$
---|---|---|---
$x_{1}^{2}+x_{2}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 625$ | 15,20 | 35 | 5
$x_{1}^{3}+x_{2}^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1008$ | 2,10 | 33 | 4
$x_{1}^{4}+x_{2}^{4}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1921$ | 6,5 | 12 | 3
⋮ |  |  |
$x_{1}^{7}+x_{2}^{7}\leavevmode\nobreak\ =\leavevmode\nobreak\ 4799353$ | 9,4 | 26 | 4
⋮ |  |  |
$x_{1}^{15}+x_{2}^{15}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1088090731$ | 4,3 | 12 | 3

Table 1: A selection of Diophantine equations and corresponding solutions
found using the D-Wave’s quantum annealer. $N_{\text{aux}}$ is the number of
auxilliary qubits necessary to quadratise the Hamiltonian, and $\beta$ is the
number of qubits used to encode each variable (see Eq. (3)).

To test this method further, we move on to Diophantine equations with
increasingly many variables. As one would expect, the higher the number of
variables, the higher will be the average number of interactions per qubit.
This often means that when the connectivity of the problem exceeds the native
connections supported by the D-Wave Quantum Processor Unit (QPU), a single
binary variable in the quadratic optimization problem needs to be represented
by two (or more) qubits (called a ‘chain’) instead of one. This procedure,
known as embedding is carried out by an embedding algorithm, and should be
carefully monitored as it can lead to so-called broken-chains that have two or
more physical qubits in the same chain taking different values. This
ultimately limits the size of problems that can be solved on quantum
annealers, while performing classical annealing on the same Ising Hamiltonian
turns out to be successful in all the examples treated in Table 2, where we
list the equations solved using both classical (cyan) and quantum (blue)
annealing.

Equation | $N_{\text{aux}}$ | $\beta$ | Solutions
---|---|---|---
$x_{1}^{2}+x_{2}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 149$ | 20 | 4 | (7, 10)
$x_{1}^{2}+x_{2}^{2}+x_{3}^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ 244$ | 18 | 4 | (12, 8, 6)
$x_{1}^{2}+\dots+x_{4}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 295$ | 24 | 4 | (14, 7, 1, 7), (11, 5, 7, 10)
$x_{1}^{2}+\dots+x_{5}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 325$ | 30 | 4 | (4, 8, 8, 9, 10)
$x_{1}^{2}+\dots+x_{6}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 420$ | 36 | 4 | (1, 11, 9, 10, 6, 9)
$x_{1}^{2}+\dots+x_{7}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 450$ | 70 | 5 |  | (10, 4, 17, 2, 6, 1, 2)
---
(2, 7, 2, 2, 18, 8, 1)
…
$x_{1}^{2}+\dots+x_{8}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 590$ | 80 | 5 |  | (9, 11, 1, 6, 5, 14, 3, 11)
---
(1, 11, 8, 9, 4, 15, 1, 9)
…
$x_{1}^{2}+\dots+x_{9}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 720$ | 90 | 5 |  | (14, 2, 13, 8, 9, 6, 8, 9, 5)
---
(7, 16, 1, 1, 20, 2, 1, 2, 2)
…
$x_{1}^{2}+\dots+x_{10}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 956$ | 100 | 5 |  | (12, 2, 8, 20, 3, 11, 11, 8, 2, 5)
---
(3, 5, 7, 5, 6, 4, 2, 14, 14, 20)
…
$x_{1}^{2}+\dots+x_{11}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1502$ | 110 | 5 |  | (6, 2, 9, 22, 9, 5, 12, 23, 1, 6, 9)
---
(9, 7, 9, 23, 1, 5, 1, 23, 11, 9, 2)
…
$x_{1}^{2}+\dots+x_{12}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ 3842$ | 120 | 5 |  | (1, 21, 4, 6, 12, 2, 28, 17, 18, 29, 31, 1)
---
(17, 8, 17, 32, 28, 5, 19, 16, 17, 14, 16, 3)
…

Table 2: A list of Diophantine equations solved using both quantum (blue) and
classical (cyan) annealing.

This reproduces the results in Refs. abraham2010particle ;
Aguiar2014DiophantineEA  obtained using Particle Swarm Optimisation and Fuzzy
Adaptive Simulated Annealing respectively.

##  IV Ramanujan (TaxiCab) numbers

Having demonstrated that problems such as these can in principle be already
solved on a quantum annealer, we are now ready to move on to the more
complicated class of problems outlined in the introduction, beginning in this
section with Taxicab numbers. In line with the above discussion, classical
simulated annealing turns out to be superior currently for these problems, so
we will use that annealing method in this and the next section.

In general, finding Taxicab numbers is not a trivial task and for higher
Taxicabs, such as $\text{Ta}(7)$, $\text{Ta}(8)$ etc., only an upper bound is
known Boyer . Indeed (using the $(k,m,n)$ notation for these numbers), it is
interesting to note that no $(5,2,2)$ numbers have been found, despite
searches up to $10^{26}$ (see Ref. Guy:1994a ).

Let us show explicitly how we use the reduction technique of Subsection II.1
to construct an Ising Model Hamiltonian whose ground states are precisely the
Taxicab numbers we want to find. As a first example we focus on
$\text{Ta}(2)$, i.e. we want to find four non-negative integer numbers such
that

$a^{3}+b^{3}\leavevmode\nobreak\ =\leavevmode\nobreak\ c^{3}+d^{3}\,,\quad
a\neq\\{c,d\\}\,.$ (16)

We again use binary encoding (see Eq. (3)) with $\beta\leavevmode\nobreak\
=\leavevmode\nobreak\ 5$, $s_{i}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1$
(numbers from $1$ to $32$) and $t_{i}\in\\{a,b,c,d\\}$. To impose the equality
between the two sums of cubes we define the following Hamiltonian

$H_{D}\leavevmode\nobreak\ =\leavevmode\nobreak\
(a^{3}+b^{3}-c^{3}-d^{3})^{2}\,.$ (17)

However, this is not the end of the story as must also encode the constraint
$a\neq\\{c,d\\}$ to avoid all the trivial minima of the above Hamiltonian,
which occur when $a=c$ and $b=d$ or vice versa. In other words we want to
construct the $H_{C}$ Hamiltonian such that it has its global minimum when
$a\neq c$ and $a\neq d$. It is more straightforward to write such a constraint
Hamiltonian directly in terms of binary variables $\tau_{i,k}$, where
$i\in\\{a,b,c,d\\}$ and $k=0,...,\beta-1$. It is easy to see that the
Hamiltonian

$\displaystyle H_{C}\leavevmode\nobreak\ $
$\displaystyle\equiv\leavevmode\nobreak\ H_{\delta}(a,c)+H_{\delta}(a,d)$
$\displaystyle\qquad\equiv\prod_{k=0}^{\beta-1}\left(1-(\tau_{a,k}-\tau_{c,k})^{2}\right)$
$\displaystyle\qquad\qquad+\prod_{k=0}^{\beta-1}\left(1-(\tau_{a,k}-\tau_{d,k})^{2}\right)$
(18)

achieves this. Explicitly, one finds that

$H_{C}\leavevmode\nobreak\ =\leavevmode\nobreak\ \begin{cases}0,\text{ when
}a\neq c\text{ and }a\neq d\,,\\\ 1,\text{ when }a\leavevmode\nobreak\
=\leavevmode\nobreak\ c\text{ and }a\neq d\,,\\\ 1,\text{ when
}a\leavevmode\nobreak\ =\leavevmode\nobreak\ d\text{ and }a\neq c\,,\\\
2,\text{ when }a\leavevmode\nobreak\ =\leavevmode\nobreak\
c\leavevmode\nobreak\ =\leavevmode\nobreak\ d\,.\end{cases}$ (19)

The Hamiltonian we shall use is then the sum

$\tilde{H}\leavevmode\nobreak\ =\leavevmode\nobreak\ H_{D}+H_{C}\,.$ (20)

Written is terms of $\tau$’s, this Hamiltonian is a polynomial of order
$2\beta$ for $\beta\leavevmode\nobreak\ \geq\leavevmode\nobreak\ 3$. Again,
setting $\beta=5$ and using the technique described in Sec. II we can reduce
it to a quadratic Hamiltonian by adding $98$ auxilliary spins.

Figure 3: Representation of the Ising Hamiltonian corresponding to the Ta$(2)$
problem.

In Fig. 3 we represent the reduced Hamiltonian for $\beta=4$. We see that
stronger couplings are rare among the interactions, which mostly form a very
complex network of weaker couplings in the background. Classical annealing on
the reduced Hamiltonian yields all the solutions written in Eq. (1), namely
all the Taxicab numbers with $a,b,c,d\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 32$.

Let us now push this further and attempt to solve more complicated
generalisation of the Taxicab problem, $(k,m,n)$ where

$\displaystyle(k,m,n)\leavevmode\nobreak\ $
$\displaystyle\equiv\leavevmode\nobreak\ a_{1}^{k}\,+\,...\,+\,a_{m}^{k}$
$\displaystyle=\leavevmode\nobreak\ b_{1}^{k}\,+\,...\,+\,b_{n}^{k},$ (21)

where $\\{a_{1},...,a_{m}\\}\neq\\{b_{1},...,b_{n}\\}$. Beginning with
$(4,3,3)$ numbers,

$\displaystyle(4,3,3)\leavevmode\nobreak\ $
$\displaystyle=\leavevmode\nobreak\ a^{4}\,+\,b^{4}\,+\,c^{4}$
$\displaystyle=\leavevmode\nobreak\ d^{4}\,+\,e^{4}\,+\,f^{4}\,,$ (22)

we define the following Hamiltonians

$H_{D}\leavevmode\nobreak\ =\leavevmode\nobreak\
(a^{4}+b^{4}+c^{4}-d^{4}-e^{4}-f^{4})^{2}\,,$ (23)

and

$H_{C}\leavevmode\nobreak\ =\leavevmode\nobreak\
H_{\delta}(a,d)+H_{\delta}(a,e)+H_{\delta}(a,f)\,,$ (24)

to impose the equality in Eq. (22) and also to force $a\neq d,e,f$. The order
of the complete Hamiltonian, which is the sum of Eq. 24 and Eq. 23, is
$2\beta$ for $\beta\geq 4$. Again, it can be reduced to a quadratic one by
adding $66$ auxilliary variables in the case with $\beta=4$ and $154$ in the
case with $\beta=5$. Several anneal runs (each with $10000$ reads) with
$\beta=4,5$ yield the following results

$(4,3,3)$ | $a$ | $b$ | $c$ | $d$ | $e$ | $f$
---|---|---|---|---|---|---
2673 | 3 | 6 | 6 | 7 | 2 | 4
16562 | 9 | 1 | 10 | 11 | 6 | 5
28593 | 2 | 13 | 2 | 9 | 6 | 12
35378 | 13 | 4 | 9 | 11 | 12 | 1
43218 | 11 | 13 | 2 | 14 | 7 | 7
54977 | 4 | 8 | 15 | 9 | 14 | 10
195122 | 21 | 5 | 2 | 9 | 13 | 20
324818 | 14 | 9 | 23 | 21 | 2 | 19
619337 | 28 | 8 | 5 | 7 | 26 | 20
847602 | 1 | 25 | 26 | 29 | 19 | 10
1071713 | 12 | 32 | 7 | 28 | 26 | 3
1178898 | 29 | 11 | 26 | 1 | 32 | 19
1328498 | 29 | 9 | 28 | 23 | 32 | 3
Table 3: List of $(4,3,3)$ numbers found using $\beta\leavevmode\nobreak\
=\leavevmode\nobreak\ 3,4,5$ and $10000$ reads per anneal run.

To comment on the efficacy of the method: the search space is of order
$32^{6}\sim 10^{9}$, and yet these solutions are found after order $10^{5}$
reads.

For the remainder of this section we consider the $(3,n,m)$ numbers, where
$n,m\in\mathbb{N^{+}}$. For this purpose we define the following Hamiltonian

$\tilde{H}\leavevmode\nobreak\ =\leavevmode\nobreak\
\left(\sum_{i=1}^{n}a_{i}^{3}-\sum_{i=1}^{m}b_{i}^{3}\right)^{2}\,,$ (25)

where in this case we do not add any constraint Hamiltonian to enforce
$\\{a_{i}\\}\neq\\{b_{i}\\}$ when $n=m$, because here it is sufficient to
simply check at the end of each anneal run if the minimum is trivial or not.

Tables 4-5 list some of the solutions found for $n=1$ and $m=5,7$ with
$\beta=5,6$. Cases with $n=m\equiv N$, with $N=6,7,8$ are listed in Tables
6-8.

(3,1,5) | $a_{1}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$
---|---|---|---|---|---|---
729 | 9 | 1 | 3 | 4 | 5 | 8
1728 | 12 | 3 | 10 | 4 | 8 | 5
68921 | 41 | 3 | 17 | 21 | 28 | 32
125000 | 50 | 2 | 8 | 24 | 36 | 40
185193 | 57 | 16 | 17 | 30 | 40 | 44
216000 | 60 | 11 | 16 | 25 | 45 | 47
262144 | 64 | 9 | 18 | 31 | 44 | 52
Table 4: A list of $(3,1,5)$ numbers found using $\beta\leavevmode\nobreak\ =\leavevmode\nobreak\ 5,6$. The reduction needs $120$ and $216$ auxilliary spins respectively. (3,1,7) | $a_{1}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$ | $b_{6}$ | $b_{7}$
---|---|---|---|---|---|---|---|---
2744 | 14 | 2 | 3 | 5 | 7 | 8 | 9 | 10
13824 | 24 | 3 | 5 | 8 | 9 | 13 | 15 | 19
32768 | 32 | 1 | 6 | 15 | 16 | 17 | 20 | 23
148877 | 53 | 3 | 21 | 24 | 28 | 29 | 32 | 36
205379 | 59 | 5 | 12 | 13 | 18 | 23 | 43 | 47
238328 | 62 | 17 | 20 | 22 | 31 | 32 | 38 | 46
Table 5: A list of $(3,1,7)$ numbers found using $\beta=5,6$. The reduction needs $160$ and $288$ auxilliary spins respectively. $(3,6,6)$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$ | $b_{6}$
---|---|---|---|---|---|---|---|---|---|---|---|---
5012 | 2 | 4 | 5 | 7 | 12 | 14 | 3 | 6 | 8 | 9 | 11 | 13
7975 | 1 | 7 | 8 | 10 | 14 | 15 | 3 | 4 | 9 | 11 | 12 | 16
8309 | 1 | 5 | 7 | 10 | 14 | 16 | 4 | 6 | 9 | 12 | 13 | 15
41873 | 3 | 8 | 11 | 14 | 26 | 27 | 1 | 6 | 13 | 19 | 22 | 28
48438 | 9 | 13 | 17 | 20 | 22 | 28 | 1 | 4 | 15 | 18 | 23 | 30
51318 | 1 | 10 | 15 | 17 | 21 | 32 | 3 | 5 | 9 | 16 | 28 | 29
52359 | 5 | 6 | 7 | 11 | 26 | 32 | 3 | 14 | 15 | 20 | 24 | 29
78730 | 2 | 3 | 23 | 26 | 28 | 30 | 11 | 14 | 16 | 20 | 31 | 32
86400 | 3 | 9 | 21 | 24 | 31 | 32 | 4 | 5 | 26 | 27 | 28 | 30
Table 6: $(3,6,6)$ solutions with $\beta=4,5$. The reduction needs $120$ and $240$ auxilliary spins respectively. $(3,7,7)$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$ | $a_{7}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$ | $b_{6}$ | $b_{7}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
39256 | 3 | 8 | 9 | 14 | 17 | 22 | 27 | 2 | 5 | 11 | 16 | 19 | 21 | 26
45063 | 3 | 5 | 7 | 13 | 14 | 19 | 32 | 2 | 9 | 10 | 15 | 18 | 23 | 28
46411 | 7 | 9 | 14 | 17 | 18 | 23 | 27 | 3 | 5 | 8 | 10 | 11 | 22 | 32
52094 | 1 | 6 | 13 | 17 | 22 | 23 | 28 | 3 | 7 | 14 | 16 | 18 | 21 | 31
63224 | 7 | 9 | 12 | 18 | 20 | 24 | 32 | 2 | 3 | 13 | 14 | 19 | 29 | 30
73276 | 6 | 12 | 14 | 15 | 24 | 29 | 30 | 9 | 17 | 18 | 20 | 23 | 27 | 28
77687 | 2 | 9 | 17 | 21 | 24 | 28 | 30 | 4 | 5 | 15 | 16 | 26 | 27 | 32
Table 7: $(3,7,7)$ solutions with $\beta=5$. The reduction needs $280$ auxilliary spins. $(3,8,8)$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$ | $a_{7}$ | $a_{8}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$ | $b_{6}$ | $b_{7}$ | $b_{8}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
50139 | 1 | 3 | 6 | 10 | 12 | 20 | 23 | 30 | 2 | 5 | 9 | 13 | 17 | 19 | 25 | 27
73206 | 1 | 3 | 4 | 17 | 20 | 25 | 26 | 30 | 5 | 8 | 9 | 10 | 19 | 21 | 28 | 32
78202 | 3 | 4 | 17 | 18 | 19 | 24 | 27 | 30 | 1 | 2 | 9 | 16 | 20 | 22 | 28 | 32
85418 | 2 | 3 | 9 | 16 | 18 | 23 | 31 | 32 | 6 | 10 | 14 | 15 | 24 | 26 | 27 | 30
Table 8: $(3,8,8)$ solutions with $\beta=5$. The reduction needs $320$
auxilliary spins. Note that even the smallest $(3,8,8)$ number represents a
solution in a search space of size $32^{16}\sim 10^{24}$.

Note that all the above solutions are non-trivial $(3,n,m)$ numbers, in that
they are not sums of smaller solutions. Indeed, it may happen that a $(3,n,m)$
number is actually the sum of $(3,p,q)$ and $(3,k,l)$ numbers with
$p+k\leavevmode\nobreak\ =\leavevmode\nobreak\ n$ and $q+l\leavevmode\nobreak\
=\leavevmode\nobreak\ m$. In order to avoid such trivial solutions, we have
simply removed them by hand at the end of each anneal run.

##  V Anomaly cancellation in the Standard Model with an extra $U(1)$

Having road-tested our reduction methods on simple problems, we now turn to a
physical application, namely the anomaly cancellation conditions in the
Standard Model extended by an extra $U(1)$ gauge symmetry. This is one of the
simplest and most studied extensions of the Standard Model (see Ref.
Langacker:2009  for a review of $Z^{\prime}$ physics), and it has been the
target of numerous experimental searches zyla2020 .

The generalities of anomaly cancellation for such systems have been discussed
in Refs. Appelquist:2003 ; Batra:2005rh ; liu2012searching ; ismail2017axial ;
Allanach:2018vjg ; tang2018flavor ; Ellis:2018 ; Allanach:2019uuu ;
Costa:2019anomaly ; Smolkovivc2019anomaly ; Allanach:2020geometric ;
Costa:2020 ; Allanach:2020zna ; Allanach:2021bfe ; Podo:2022gyj . In this work
we will for concreteness specialise to the models studied in Ref.
Allanach:2020zna . Here, the main assumption is that the chiral fermions
appear in the usual 3 families of quarks and leptons, together with 3 right-
handed neutrinos. The charges under the additional $U(1)$ are labelled by
$\\{Q_{i},U_{i},D_{i},L_{i},E_{i},N_{i}\\}$, respectively, with
$i\in\\{1,2,3\\}$ indicating the generation number. Under this assumption the
anomaly cancellation condition yields the following set of Diophantine
equations for the charges:

$\displaystyle\sum_{i=1}^{3}\,(6Q_{i}+3U_{i}+3D_{i}+2L_{i}+E_{i}+N_{i})\leavevmode\nobreak\
=\leavevmode\nobreak\ 0\,,$ (26)
$\displaystyle\sum_{i=1}^{3}\,(3Q_{i}+L_{i})\leavevmode\nobreak\
=\leavevmode\nobreak\ 0\,,$ (27)
$\displaystyle\sum_{i=1}^{3}\,(2Q_{i}+U_{i}+D_{i})\leavevmode\nobreak\
=\leavevmode\nobreak\ 0\,,$ (28)
$\displaystyle\sum_{i=1}^{3}\,(Q_{i}+8U_{i}+2D_{i}+3L_{i}+6E_{i})\leavevmode\nobreak\
=\leavevmode\nobreak\ 0\,,$ (29)
$\displaystyle\sum_{i=1}^{3}\,(Q_{i}^{2}-2U_{i}^{2}+D_{i}^{2}-L_{i}^{2}+E_{i}^{2})\leavevmode\nobreak\
=\leavevmode\nobreak\ 0\,,$ (30)
$\displaystyle\sum_{i=1}^{3}\,(6Q_{i}^{3}+3U_{i}^{3}+3D_{i}^{3}+2L_{i}^{3}+E_{i}^{3}+N_{i}^{3})\leavevmode\nobreak\
=\leavevmode\nobreak\ 0\,.$ (31)

A general solution to the above equations has already been found analytically
in Ref. Allanach:2020zna . However, we shall demonstrate here that these
problems can be also tackled using Ising model annealing (in practice here we
use simulated annealing, but ultimately quantum annealers will be
practicable).

As for the Taxicab problem, we begin constructing the Hamiltonian by simply
squaring and summing the left hand side of all the above equations. We encode
all the variables involved as in Eq. (3) with
$t_{i}\in\\{Q_{i},U_{i},D_{i},L_{i},E_{i},N_{i}\\}$ and take $s_{i}=-1$ for
all the charges. We set $\beta=2$, thus looking for solutions with entries
from $-1$ to $2$. Note that although the number of bits we use to represent
each variable is relatively low, the number of possible configurations of
these $3\times 6\leavevmode\nobreak\ =\leavevmode\nobreak\ 18$ charges with
possible values in $[-1,2]$ is already quite high: $4^{18}\sim 10^{10}$. It is
worth mentioning that in this particular case a comprehensive scan can be
completed with far fewer attempts due to generation permutation symmetry in
the equations. Indeed, it is easy to see that the anomaly equations are
invariant under arbitrary permutations of $\\{A_{1},A_{2},A_{3}\\}$, where
$A\in\\{Q,U,D,L,E,N\\}$, giving an $(S_{3})^{6}$ permutational symmetry that
could be exploited if we were looking for solutions by exhaustive scanning
over all the different configurations.

Of course our goal here is to avoid using such tricks, but to instead find
solutions using annealing on the reduced Ising Hamiltonian. For $\beta=2$ the
reduction requires only $18$ auxiliaries. We have performed several anneal
runs with $10000$ reads obtaining an average of $60$ distinct solutions per
anneal run. In the following table we present a sample of three of them.

$Q_{1}$ | $Q_{2}$ | $Q_{3}$ | $U_{1}$ | $U_{2}$ | $U_{3}$ | $D_{1}$ | $D_{2}$ | $D_{3}$ | $L_{1}$ | $L_{2}$ | $L_{3}$ | $E_{1}$ | $E_{2}$ | $E_{3}$ | $N_{1}$ | $N_{2}$ | $N_{3}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 0 | -1 | -1 | 1 | 0 | 1 | -1 | 0 | -1 | 0 | 1 | -1 | 0 | 1 | 0 | 1 | -1
-1 | 1 | 0 | 0 | -1 | 1 | -1 | 1 | 0 | -1 | 0 | 1 | 0 | 1 | -1 | 1 | 0 | -1
1 | -1 | 0 | -1 | 1 | 0 | 1 | 0 | -1 | 1 | 0 | -1 | -1 | 1 | 0 | -1 | 0 | 1

Table 9: A sample of three solutions found using $\beta\leavevmode\nobreak\
=\leavevmode\nobreak\ 2$ and $10000$ reads in each anneal run.

One might expect higher values of $\beta$ to lead to new solutions with bigger
entries, along with those previously found. However, for this specific
problem, classical annealing turns out to be unfruitful for $\beta>2$. To
explain why it is useful to inspect how the energy gap $\Delta$ between the
ground state and the first excited state scales as a function of the size of
the problem. It can be shown (see Ref. chang2019least ) that

$\Delta\leavevmode\nobreak\ \sim\leavevmode\nobreak\
\mathcal{O}\left(\frac{2^{-n\mu}}{m\alpha}\right)\leavevmode\nobreak\ ,$ (32)

where $\alpha$ and $\mu$ are the number of additions and multiplications
respectively in the Hamiltonian written in terms of binary variables, $m$ is
the number of equations we want to solve and $n$ is the the effective
precision, which is the difference between the largest and smallest nonzero
absolute values representable among all the variables in the system.
Increasing $\beta$ makes all these parameter bigger, including $n$ and $\mu$,
causing an exponential shrinkage of the energy gap between the ground states
and the first excited states, which in turn considerably affects the
algorithm’s performance.

To improve our results and find solutions with bigger entries we use the
solution-mining method described in Sec. II.2. After $30$ anneal runs this
yields $153$ solutions with entries between $-13$ and $13$. Note that a
complete scan on all possible such configurations, even exploiting the
$(S_{3})^{6}$ permutational symmetry, would involve $\binom{13\times
2+1}{3}^{6}\sim 10^{20}$ trials, which is infeasible with conventional
computing methods. It should be noted that we do not make use of the
permutational symmetry and the Ising machine is in principle succeeding within
a search space of $26^{18}\sim 3\times 10^{25}$, although it is not yet clear
how exhaustive the method of small $\beta$ plus solution mining can eventually
be.

In Table 10 we present a sample of ten of the solutions found.

$Q_{1}$ | $Q_{2}$ | $Q_{3}$ | $U_{1}$ | $U_{2}$ | $U_{3}$ | $D_{1}$ | $D_{2}$ | $D_{3}$ | $L_{1}$ | $L_{2}$ | $L_{3}$ | $E_{1}$ | $E_{2}$ | $E_{3}$ | $N_{1}$ | $N_{2}$ | $N_{3}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
-1 | 0 | 1 | -1 | 0 | 1 | 1 | -1 | 0 | 1 | 0 | -1 | 0 | -1 | 1 | 1 | 0 | -1
0 | -2 | 2 | 1 | -1 | 2 | -2 | 0 | 0 | 0 | 1 | -1 | 0 | -1 | -1 | 0 | 1 | 1
3 | -1 | -2 | -1 | -2 | 3 | -4 | 2 | 2 | 0 | -3 | 3 | 2 | -2 | 0 | 2 | -3 | 1
3 | -2 | -1 | 1 | -3 | 3 | -4 | 3 | 0 | -1 | 0 | 1 | -1 | 0 | 0 | 3 | -3 | 1
-1 | 1 | 0 | -2 | -1 | 4 | -5 | 4 | 0 | -2 | -1 | 3 | 0 | 2 | -3 | 1 | -2 | 2
1 | -1 | 0 | 0 | -2 | 5 | -6 | 4 | -1 | -1 | 0 | 1 | 0 | -1 | -2 | 0 | -2 | 5
1 | 0 | 0 | -1 | -2 | 6 | -7 | 4 | -2 | -3 | 0 | 0 | 2 | 1 | -4 | 0 | 0 | 7
2 | -1 | -1 | 2 | -3 | 4 | -6 | 2 | 1 | 0 | 0 | 0 | -3 | 1 | -1 | -3 | -2 | 8
2 | -2 | -2 | 2 | 1 | 2 | -2 | -1 | 2 | 6 | -1 | 1 | -3 | -1 | -5 | -10 | -2 | 9
1 | -3 | 0 | 1 | 5 | 2 | -2 | 1 | -3 | 2 | 1 | 3 | -5 | -3 | -4 | -13 | 0 | 13

Table 10: A sample of ten solutions found using the solution-mining method.

The first of these solutions is equivalent to one found previously in Table 9.
This is because in the first anneal run the algorithm looks for solutions
centered around zero, i.e. with entries between $[-1,2]$. Then it starts
exploring the neighborhood of the solution found in the previous anneal run,
gradually finding solutions with larger entries.

##  VI Conclusions

Diophantine problems in physics and beyond are often computationally hard. In
this paper we have investigated the use of Ising machines for finding
solutions to such problems, and have shown that they can succeed within search
spaces that are vast. For example finding the (3,8,8) number 50139, the lowest
number we find that can be written as the sum of eight cubes in two different
ways, represents a search in a space of size $10^{24}$. Furthermore using
“solution mining” for the task of anomaly cancellation, we find a proven
ability to find solutions in search spaces of order $10^{26}$.

The methods described here are valid for any Ising machine, including both
quantum and simulated annealers. We found that currently available quantum
annealers can already solve a large variety of Diophantine problems at a
relatively advanced level, but they are not yet competitive with their
simulated counterparts. Nevertheless there are several reasons to believe that
ultimately quantum annealers will become dominant for such tasks. The first is
to do with the current obstacle to performing higher order tasks on a quantum
annealer which is the fact that with the increasing complexity of the
equations, finding a suitable embedding (i.e. an embedding such that the
connectivity and the number of auxilliary qubits required are within the
limits imposed by the quantum processing unit) is a non trivial task. Indeed
it is worth noting that deciding whether a graph $G$ can be embedded in a
graph $H$ is itself an NP-complete problem (when $H$ is arbitrary but as in
current annealers $G$ is built out of Chimera or Pegasus sub-graphs)
lobe2021minor . Thus one expects significant improvement in embedding as the
overall size of the annealer, and the connectivity of its sub-graphs is
continually increased.

The second reason to be optimistic about quantum annealers is in the potential
speed-up in the way that they find global minima. For example there are many
techniques open to quantum annealers such as diabatic annealing that have the
potential to avoid anneal times increasing exponentially with the difficulty
of the problem diabatic , an issue that is seen in both adiabatic quantum
annealing and in simulated annealing. These physical aspects of quantum
annealing which were also crucial in the quantum field theory tunneling
studies in Refs. Abel:2020ebj ; Abel:2020qzm ; Abel:2021fpn  cannot be
efficiently simulated classically. Indeed discrete problems generally favour
quantum annealers because in a sense these machines can operate by performing
a quantum gradient descent by tunnelling. By contrast any simulation method
requires a defined dynamical process for its evolution, in order to hop
between potential solutions, and this tends to become increasingly delicate
with the difficulty of the problem. Thus for certain problems quantum
tunneling could be an enormous advantage. For example, it has already been
shown that quantum annealing overcomes simulated annealing in a large variety
of cases (see for example Refs. santoro2002theory ; santoro2004 ;
denchev2016computational ). Particularly interesting are the results found in
Ref. denchev2016computational  showing quantum annealing outperforming
classical annealing by a factor of $10^{8}$ in finding the minimum of a
particular crafted problem with tall and narrow energy barriers separating
local minima. This is precisely the kind of configuration one expects when
solving Diophantine problems.

In summary, we believe that the methods presented here should be an effective
heuristic search method for the many discrete problems one encounters in
physics. In particular it will be interesting to employ them in the string
theory landscape context, and compare them to the genetic algorithmic and
machine learning techniques that have been studied in Refs. Abel:2014xta ;
He:2017 ; Ruehle:2017mzq ; Carifio:2017bov ; Bull:2018 ; Mutter:2018sra ;
Klaewer:2018sfl ; Brodie:2019dfx ; Bull:2019cij ; Cole:2019enn ;
Halverson:2019tkf ; Ruehle:2020jrk ; Krippendorf:2021uxu ; AbdusSalam:2020ywo
; Bena:2020xrh ; Larfors:2020ugo ; Constantin:2021for ; Abel:2021rrj ;
Loges:2021hvn .

Acknowledgements: We would like to thank Ben Allanach, Andrei Constantin, Juan
Criado, Jim Halverson, Thomas Harvey, Andre Lukas, Fabian Ruehle and Michael
Spannowsky for helpful discussions. S.A. is supported by the STFC under grant
ST/P001246/1.

## References

  * (1) G. ’t Hooft, “Naturalness, chiral symmetry, and spontaneous chiral symmetry breaking,” NATO Sci. Ser. B 59 (1980) 135–157.
  * (2) B. C. Allanach, J. Davighi, and S. Melville, “An Anomaly-free Atlas: charting the space of flavour-dependent gauged $U(1)$ extensions of the Standard Model,” JHEP 02 (2019) 082, [arXiv:1812.04602 [hep-ph]]. [Erratum: JHEP 08, 064 (2019)].
  * (3) B. C. Allanach, B. Gripaios, and J. Tooby-Smith, “Solving local anomaly equations in gauge-rank extensions of the Standard Model,” Phys. Rev. D 101 no. 7, (2020) 075015, [arXiv:1912.10022 [hep-th]].
  * (4) B. C. Allanach, B. Gripaios, and J. Tooby-Smith, “Anomaly cancellation with an extra gauge boson,” Phys. Rev. Lett. 125 no. 16, (2020) 161601, [arXiv:2006.03588 [hep-th]].
  * (5) B. C. Allanach, B. Gripaios, and J. Tooby-Smith, “Semisimple extensions of the Standard Model gauge algebra,” Phys. Rev. D 104 no. 3, (2021) 035035, [arXiv:2104.14555 [hep-th]].
  * (6) S. Tung, “Computational complexities of diophantine equations with parameters,” Journal of Algorithms 8 no. 3, (1987) 324–336. https://www.sciencedirect.com/science/article/pii/0196677487900137.
  * (7) F. Denef and M. R. Douglas, “Computational complexity of the landscape. I.,” Annals Phys. 322 (2007) 1096–1142, [arXiv:hep-th/0602072].
  * (8) F. Denef, M. R. Douglas, B. Greene, and C. Zukowski, “Computational complexity of the landscape II—Cosmological considerations,” Annals Phys. 392 (2018) 93–127, [arXiv:1706.06430 [hep-th]].
  * (9) J. Halverson and F. Ruehle, “Computational Complexity of Vacua and Near-Vacua in Field and String Theory,” Phys. Rev. D 99 no. 4, (2019) 046015, [arXiv:1809.08279 [hep-th]].
  * (10) J. Halverson, M. Plesser, F. Ruehle, and J. Tian, “Kähler Moduli Stabilization and the Propagation of Decidability,” Phys. Rev. D 101 no. 4, (2020) 046010, [arXiv:1911.07835 [hep-th]].
  * (11) T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse ising model,” Phys. Rev. E 58 (Nov, 1998) 5355–5363. https://link.aps.org/doi/10.1103/PhysRevE.58.5355.
  * (12) E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, “A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem,” Science 292 no. 5516, (2001) 472–475, [https://www.science.org/doi/pdf/10.1126/science.1057726]. https://www.science.org/doi/abs/10.1126/science.1057726.
  * (13) A. Das and B. K. Chakrabarti, “Colloquium: Quantum annealing and analog quantum computation,” Rev. Mod. Phys. 80 (Sep, 2008) 1061–1081. https://link.aps.org/doi/10.1103/RevModPhys.80.1061.
  * (14) A. Lucas, “Ising formulations of many np problems,” Frontiers in Physics 2 (2014) . https://www.frontiersin.org/article/10.3389/fphy.2014.00005.
  * (15) P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, “Perspectives of quantum annealing: methods and implementations,” Reports on Progress in Physics 83 no. 5, (May, 2020) 054401. https://doi.org/10.1088/1361-6633/ab85b8.
  * (16) N. Mohseni, P. McMahon, and T. Byrnes, “Ising machines as hardware solvers of combinatorial optimization problems.” Nature Reviews Physics (2022) . https://doi.org/10.1038/s42254-022-00440-8.
  * (17) R. K. Guy, “Sums of like powers. euler’s conjecture, §d1 in unsolved problems in number theory,” 1994.
  * (18) J.-C. Meyrignac, “Computing minimal equal sums of like powers.”  http://euler.free.fr.
  * (19) Eulernet.  http://euler.free.fr/database.txt.
  * (20) E. W. Weisstein, “Diophantine equation–nth powers (replace “n” by the power): From mathworld–a wolfram web resource.”  https://mathworld.wolfram.com/DiophantineEquationnthPowers.html.
  * (21) T. Piezas, “A collection of algebraic identities.”  https://sites.google.com/site/tpiezas/Home.
  * (22) A. Gadde, S. Gukov, and P. Putrov, “(0, 2) trialities,” JHEP 03 (2014) 076, [arXiv:1310.0818 [hep-th]].
  * (23) S. Franco, S. Lee, and R.-K. Seong, “Brane brick models and 2d (0, 2) triality,” JHEP 05 (2016) 020, [arXiv:1602.01834 [hep-th]].
  * (24) S. Franco, S. Lee, R.-K. Seong, and C. Vafa, “Quadrality for Supersymmetric Matrix Models,” JHEP 07 (2017) 053, [arXiv:1612.06859 [hep-th]].
  * (25) N. S. Dattani and N. Bryans, “Quantum factorization of 56153 with only 4 qubits,” arXiv preprint (2014) . https://arxiv.org/abs/1411.6758.
  * (26) R. Tanburn, E. Okada, and N. Dattani, “Reducing multi-qubit interactions in adiabatic quantum computation without adding auxiliary qubits. part 1: The "deduc-reduc" method and its application to quantum factorization of numbers,” 2015.  https://arxiv.org/abs/1508.04816.
  * (27) S. Jiang, K. A. Britt, A. J. McCaskey, T. S. Humble, and S. Kais, “Quantum annealing for prime factorization,” 2018.  https://arxiv.org/abs/1804.02733.
  * (28) W. Peng, B. Wang, F. Hu, Y. Wang, X. Fang, X. Chen, and C. Wang, “Factoring larger integers with fewer qubits via quantum annealing with optimized parameters,” SCIENCE CHINA Physics, Mechanics & Astronomy 62 no. 6, (2019) 1–8. https://doi.org/10.1007/s11433-018-9307-1.
  * (29) R. H. Warren, “Factoring on a quantum annealing computer,” Quantum Information & Computation 19 no. 3-4, (2019) 252–261.
  * (30) B. Wang, F. Hu, H. Yao, and C. Wang, “Prime factorization algorithm based on parameter optimization of Ising model,” Scientific reports 10 no. 1, (2020) 1–10. https://www.nature.com/articles/s41598-020-62802-5.
  * (31) C. C. Chang, A. Gambhir, T. S. Humble, and S. Sota, “Quantum annealing for systems of polynomial equations,” Scientific Reports 9 no. 1, (Jul, 2019) . https://doi.org/10.1038%2Fs41598-019-46729-0.
  * (32) T. H. Chang, T. C. Lux, and S. S. Tipirneni, “Least-squares solutions to polynomial systems of equations with quantum annealing,” Quantum Information Processing 18 no. 12, (2019) 1–17. https://doi.org/10.1007/s11128-019-2489-x.
  * (33) I. G. Rosenberg, “Reduction of bivalent maximization to the quadratic case.” Cahiers du Centre d’Etudes de Recherche Operationnelle 17 (1975) 71–74.
  * (34) N. Dattani, “Quadratization in discrete optimization and quantum mechanics,” 2019\.  https://arxiv.org/abs/1901.04405.
  * (35) S. Abel, J. C. Criado, and M. Spannowsky, “Completely Quantum Neural Networks,” [arXiv:2202.11727 [quant-ph]].
  * (36) R. Babbush, B. O'Gorman, and A. Aspuru-Guzik, “Resource efficient gadgets for compiling adiabatic quantum optimization problems,” Annalen der Physik 525 no. 10-11, (Sep, 2013) 877–888. https://doi.org/10.1002%2Fandp.201300120.
  * (37) E. Boros and P. L. Hammer, “Pseudo-Boolean optimization,” Discrete Applied Mathematics 123 no. 1-3, (2002) 155–225. https://doi.org/10.1016/S0166-218X(01)00341-9.
  * (38) V. Chvatal, “A greedy heuristic for the set-covering problem,” Mathematics of operations research 4 no. 3, (1979) 233–235. https://doi.org/10.1287/moor.4.3.233.
  * (39) H. Aguiar, “Diophantine Equations and Fuzzy Adaptive Simulated Annealing,”  2014\.
  * (40) M. Moradi, “Solving ill-conditioned linear equations using Simulated Annealing method,” Journal of Hyperstructures 7 no. 1, (2019) . http://www.jhs-uma.com/index.php/JHSMS/article/viewFile/327/120.
  * (41) S. Abraham and M. Sanglikar, “A Diophantine equation solver-a genetic algorithm application,” in Mathematical Colloquium journal, vol. 15.  2001\.
  * (42) S. Abraham, S. Sanyal, and M. Sanglikar, “Particle swarm optimization based diophantine equation solver,” 2010.  https://arxiv.org/abs/1003.2724.
  * (43) S. Abraham, S. Sanyal, and M. Sanglikar, “Finding numerical solutions of Diophantine equations using ant colony optimization,” Applied Mathematics and Computation 219 no. 24, (2013) 11376–11387. https://doi.org/10.1016/j.amc.2013.05.051.
  * (44) O. Pérez, I. Amaya, and R. Correa, “Numerical solution of certain exponential and non-linear Diophantine systems of equations by using a discrete particle swarm optimization algorithm,” Applied Mathematics and Computation 225 (2013) 737–746. https://doi.org/10.1016/j.amc.2013.10.007.
  * (45) T. Lanting, “The D-Wave 2000Q Processor,” 2017.  https://www.youtube.com/watch?v=_Y9sVY-XBfI&index=2&list=PLAlerseOylzZw8A-R4lyU3mqfYfQSggiz. presented at AQC 2017.
  * (46) C. Boyer, “New Upper Bounds for Taxicab and Cabtaxi Numbers.”  http://www.christianboyer.com/taxicab/.
  * (47) P. Langacker, “The physics of heavy ${Z}^{{}^{\prime}}$ gauge bosons,” Rev. Mod. Phys. 81 (Aug, 2009) 1199–1228. https://link.aps.org/doi/10.1103/RevModPhys.81.1199.
  * (48) P. Z. et al. (Particle Data Group), “Review of Particle Physics,” Progress of Theoretical and Experimental Physics 2020 no. 8, (08, 2020) , [https://academic.oup.com/ptep/article-pdf/2020/8/083C01/34673722/ptaa104.pdf]. https://doi.org/10.1093/ptep/ptaa104. 083C01.
  * (49) T. Appelquist, B. A. Dobrescu, and A. R. Hopper, “Nonexotic neutral gauge bosons,” Phys. Rev. D 68 (Aug, 2003) 035012. https://link.aps.org/doi/10.1103/PhysRevD.68.035012.
  * (50) P. Batra, B. A. Dobrescu, and D. Spivak, “Anomaly-free sets of fermions,” J. Math. Phys. 47 (2006) 082301, [arXiv:hep-ph/0510181].
  * (51) J.-Y. Liu, Y. Tang, and Y.-L. Wu, “Searching for $Z^{{}^{\prime}}$ Gauge Boson in an Anomaly-Free U(1)′ Gauge Family Model,” J. Phys. G 39 (2012) 055003, [arXiv:1108.5012 [hep-ph]].
  * (52) A. Ismail, W.-Y. Keung, K.-H. Tsao, and J. Unwin, “Axial vector $Z^{\prime}$ and anomaly cancellation,” Nucl. Phys. B 918 (2017) 220–244, [arXiv:1609.02188 [hep-ph]].
  * (53) Y. Tang and Y.-L. Wu, “Flavor non-universal gauge interactions and anomalies in B-meson decays,” Chin. Phys. C 42 no. 3, (2018) 033104, [arXiv:1705.05643 [hep-ph]]. [Erratum: Chin.Phys.C 44, 069101 (2020)].
  * (54) J. Ellis, M. Fairbairn, and P. Tunney, “Anomaly-free models for flavour anomalies,” The European Physical Journal C 78 no. 3, (2018) 1–9.
  * (55) D. B. Costa, B. A. Dobrescu, and P. J. Fox, “General solution to the ${U}(1)$ anomaly equations,” Phys. Rev. Lett. 123 (Oct, 2019) 151601. https://link.aps.org/doi/10.1103/PhysRevLett.123.151601.
  * (56) A. Smolkovič, M. Tammaro, and J. Zupan, “Anomaly free Froggatt-Nielsen models of flavor,” Journal of High Energy Physics 2019 no. 10, (2019) 1–50.
  * (57) B. C. Allanach, B. Gripaios, and J. Tooby-Smith, “Geometric General Solution to the $U(1)$ Anomaly Equations,” JHEP 05 (2020) 065, [arXiv:1912.04804 [hep-th]].
  * (58) D. B. Costa, “Anomaly-free ${U(1)}^{m}$ extensions of the Standard Model,” Phys. Rev. D 102 (Dec, 2020) 115006. https://link.aps.org/doi/10.1103/PhysRevD.102.115006.
  * (59) A. Podo and F. Revello, “Integer solutions to the anomaly equations for a class of chiral gauge theories,” [arXiv:2205.03428 [hep-th]].
  * (60) E. Lobe and A. Lutz, “Minor embedding in broken chimera and pegasus graphs is np-complete,” 2021.  https://arxiv.org/abs/2110.08325.
  * (61) E. Crosson and D. Lidar, “Prospects for quantum enhancement with diabatic quantum annealing,” Nat. Rev. Phys. 3 (Sep, 2021) 466––489. https://doi.org/10.1038/s42254-021-00313-6.
  * (62) S. Abel, N. Chancellor, and M. Spannowsky, “Quantum Computing for Quantum Tunnelling,” [arXiv:2003.07374 [hep-ph]].
  * (63) S. Abel and M. Spannowsky, “Observing the fate of the false vacuum with a quantum laboratory,” P. R. X. Quantum. 2 (2021) 010349, [arXiv:2006.06003 [hep-th]].
  * (64) S. Abel, A. Blance, and M. Spannowsky, “Quantum Optimisation of Complex Systems with a Quantum Annealer,” [arXiv:2105.13945 [quant-ph]].
  * (65) G. E. Santoro, R. Martonák, E. Tosatti, and R. Car, “Theory of Quantum Annealing of an Ising Spin Glass,” Science 295 no. 5564, (Mar, 2002) 2427–2430. https://doi.org/10.1126%2Fscience.1068774.
  * (66) R. Martoňák, G. E. Santoro, and E. Tosatti, “Quantum annealing of the traveling-salesman problem,” Phys. Rev. E 70 (Nov, 2004) 057701. https://link.aps.org/doi/10.1103/PhysRevE.70.057701.
  * (67) V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, “What is the computational value of finite-range tunneling?” Physical Review X 6 no. 3, (2016) 031015. https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.031015.
  * (68) S. Abel and J. Rizos, “Genetic Algorithms and the Search for Viable String Vacua,” JHEP 08 (2014) 010, [arXiv:1404.7359 [hep-th]].
  * (69) Y.-H. He, “Machine-learning the string landscape,” Physics Letters B 774 (2017) 564–568. https://www.sciencedirect.com/science/article/pii/S0370269317308365.
  * (70) F. Ruehle, “Evolving neural networks with genetic algorithms to study the String Landscape,” JHEP 08 (2017) 038, [arXiv:1706.07024 [hep-th]].
  * (71) J. Carifio, J. Halverson, D. Krioukov, and B. D. Nelson, “Machine Learning in the String Landscape,” JHEP 09 (2017) 157, [arXiv:1707.00655 [hep-th]].
  * (72) K. Bull, Y.-H. He, V. Jejjala, and C. Mishra, “Machine learning CICY threefolds,” Physics Letters B 785 (2018) 65–72. https://www.sciencedirect.com/science/article/pii/S0370269318306117.
  * (73) A. Mütter, E. Parr, and P. K. S. Vaudrevange, “Deep learning in the heterotic orbifold landscape,” Nucl. Phys. B 940 (2019) 113–129, [arXiv:1811.05993 [hep-th]].
  * (74) D. Klaewer and L. Schlechter, “Machine Learning Line Bundle Cohomologies of Hypersurfaces in Toric Varieties,” Phys. Lett. B 789 (2019) 438–443, [arXiv:1809.02547 [hep-th]].
  * (75) C. R. Brodie, A. Constantin, R. Deen, and A. Lukas, “Machine Learning Line Bundle Cohomology,” Fortsch. Phys. 68 no. 1, (2020) 1900087, [arXiv:1906.08730 [hep-th]].
  * (76) K. Bull, Y.-H. He, V. Jejjala, and C. Mishra, “Getting CICY High,” Phys. Lett. B 795 (2019) 700–706, [arXiv:1903.03113 [hep-th]].
  * (77) A. Cole, A. Schachner, and G. Shiu, “Searching the Landscape of Flux Vacua with Genetic Algorithms,” JHEP 11 (2019) 045, [arXiv:1907.10072 [hep-th]].
  * (78) J. Halverson, B. Nelson, and F. Ruehle, “Branes with Brains: Exploring String Vacua with Deep Reinforcement Learning,” JHEP 06 (2019) 003, [arXiv:1903.11616 [hep-th]].
  * (79) F. Ruehle, “Data science applications to string theory,” Phys. Rept. 839 (2020) 1–117.
  * (80) S. Krippendorf, R. Kroepsch, and M. Syvaeri, “Revealing systematics in phenomenologically viable flux vacua with reinforcement learning,” [arXiv:2107.04039 [hep-th]].
  * (81) S. AbdusSalam, S. Abel, M. Cicoli, F. Quevedo, and P. Shukla, “A systematic approach to Kähler moduli stabilisation,” JHEP 08 no. 08, (2020) 047, [arXiv:2005.11329 [hep-th]].
  * (82) I. Bena, J. Blåbäck, M. Graña, and S. Lüst, “The tadpole problem,” JHEP 11 (2021) 223, [arXiv:2010.10519 [hep-th]].
  * (83) M. Larfors and R. Schneider, “Explore and Exploit with Heterotic Line Bundle Models,” Fortsch. Phys. 68 no. 5, (2020) 2000034, [arXiv:2003.04817 [hep-th]].
  * (84) A. Constantin, T. R. Harvey, and A. Lukas, “Heterotic String Model Building with Monad Bundles and Reinforcement Learning,” [arXiv:2108.07316 [hep-th]].
  * (85) S. Abel, A. Constantin, T. R. Harvey, and A. Lukas, “Evolving Heterotic Gauge Backgrounds: Genetic Algorithms versus Reinforcement Learning,” Fortsch. Phys. 70 no. 5, (2022) 2200034, [arXiv:2110.14029 [hep-th]].
  * (86) G. J. Loges and G. Shiu, “Breeding Realistic D-Brane Models,” Fortsch. Phys. 70 no. 5, (2022) 2200038, [arXiv:2112.08391 [hep-th]].