# Solving the homogeneous Bethe-Salpeter equation with a quantum annealer

Filippo Fornetti  Dipartimento di Fisica e Geologia, Università di Perugia,
and INFN, Sezione di Perugia, Via Alessandro Pascoli, 06123 Perugia, Italy
Alex Gnech  Department of Physics, Old Dominion University, Norfolk, Virginia
23529, USA and Theory Center, Jefferson Laboratory, Newport News, Virginia
23606, USA  Tobias Frederico  Instituto Tecnológico de Aeronáutica, 12228-900
São José dos Campos, Brazil  Francesco Pederiva  Dipartimento di Fisica,
Università di Trento, and INFN-TIFPA, Via Sommarive 14, Povo 38123, Trento,
Italy  Matteo Rinaldi  INFN, Sezione di Perugia, Via Alessandro Pascoli, 06123
Perugia, Italy and Dipartimento di Fisica e Geologia, Università di Perugia,
Italy  Alessandro Roggero  Dipartimento di Fisica, Università di Trento, and
INFN-TIFPA, Via Sommarive 14, Povo 38123, Trento, Italy  Giovanni Salmè  INFN,
Sezione di Roma, P.le A. Moro 2, 00185 Rome, Italy  Sergio Scopetta
<EMAIL_ADDRESS>Dipartimento di Fisica e Geologia, Università di
Perugia, and INFN, Sezione di Perugia, Via Alessandro Pascoli, 06123 Perugia,
Italy  Michele Viviani  INFN, Sezione di Pisa, Largo Pontecorvo 3, 56100,
Pisa, Italy

###### Abstract

The homogeneous Bethe-Salpeter equation (hBSE), describing a bound system in a
genuinely relativistic quantum-field theory framework, was solved for the
first time by using a D-Wave quantum annealer. After applying standard
techniques of discretization, the hBSE, in ladder approximation, can be
formally transformed in a generalized eigenvalue problem (GEVP), with two
square matrices: one symmetric and the other non symmetric. The latter matrix
poses the challenge of obtaining a suitable formal approach for investigating
the non symmetric GEVP by means of a quantum annealer, i.e to recast it as a
quadratic unconstrained binary optimization problem. A broad numerical
analysis of the proposed algorithms, applied to matrices of dimension up to
64, was carried out by using both the proprietary simulated-annealing package
and the D-Wave Advantage 4.1 system. The numerical results very nicely compare
with those obtained with standard classical algorithms, and also show
interesting scalability features.

##  I Introduction

The past decade has seen an exponential growth of research activity and
interest in quantum computing in an impressive number of fields (see, e.g.,
Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] for recent reviews which could offer
an admittedly partial idea of the state-of-the-art). The driving force behind
all of this has been the astonishing advances in new quantum-hardware
implementation which occurs with close frequency, as shown by the long list of
major achievements one can find in literature and media (see, e.g. Ref. [12]
for a an introduction of the main quantum-computing platforms, as well as
algorithms). However, one should recall that an easy path toward a universal,
fault-tolerant quantum computer is currently hindered by noise, decoherence
and scalability [1].

In general, one can perform calculations that take advantage of either digital
gate-based quantum computing or analog quantum annealing (see, e.g., Ref. [13]
for a detailed discussion of the two approaches). The first choice is more
widely adopted, since it is meant as the quantum analog of a classical
general-purpose computer, and it has the possibility to address a broader set
of problems (see, e.g., Ref. [14], for recent experimental results obtained
using an IBM quantum machine).The second possibility, the quantum annealing,
is mainly devoted to optimization problems, that can be properly encoded into
an Ising Hamiltonian. In this case, one searches for the ground-state of the
system by exploiting quantum fluctuations (see, e.g., Refs. [15, 16, 17, 18,
19] for an introduction to the foundations of quantum annealing). Quantum
annealing, that has its basis in the adiabatic theorem, have been exploited in
several practical applications (see, e.g., Refs.[20, 21, 22, 23, 24] and
references quoted therein for a non exhaustive list of quantum-annealing
applications) outperforming the results obtained with simple classical
heuristics, like simulated annealing. In this work we exploit this second
possibility to solve the non symmetric Generalized EigenValue Problem (GEVP).
In particular we use the D-wave Advantage quantum annealer (QA) [25, 26, 27]
which has a network of more than 5000 qubits, with 15-way connectivity.

The application of the quantum annealer for finding the solution of the GEVP
for symmetric matrices it is very well known in literature [28, 29, 30, 31,
32]. In this work, we aim instead to find the minimal/maximum right eigenvalue
(and corresponding eigenvector) of a non-symmetric GEVP by the use of a
quantum annealer. Notice that the proposed technique could also be used for
the left eigenvectors with simple modifications. In our specific case we focus
on the non symmetric GEVP stemming from the discretization method adopted for
solving the homogeneous Bethe-Salpeter equation (hBSE) [33, 34] in ladder
approximation, directly in the four-dimensional (4D) Minkowski space (see also
Refs. [35, 36] for a detailed introduction to the topic, and the brief resumé
in Appendix A). In particular we focus on the numerical solutions of the hBSE
for two massive scalar fields interacting through the ladder exchange of
another massive scalar field (see Ref. [37] for the results obtained with a
classical computation).

The non symmetric GEVP is given in general by the expression

$\displaystyle A~{}{\bf v}_{i}=\lambda_{i}~{}B~{}{\bf v}_{i}$ (1)

where the $n\times n$ matrix $A$ is a real and non symmetric, the $n\times n$
matrix $B$ is non singular and symmetric and, ${\bf v}_{i}$ is the eigenvector
corresponding to the i-th eigenvalue $\lambda_{i}$. It is important to
underline that this form is relevant for many areas, e.g., for the analysis of
mechanical systems [38, 39], in fluid mechanics [40], for investigating
criticality in nuclear reactors [41] or the onset of oscillatory instability
in a heated cavity [42], for lattice QCD [43], in industry [44], etc.

At the present stage, we reduce the GEVP into a standard eigenvalue problem by
inverting the matrix $B$ through a classical algorithm, before passing it to
the quantum annealer to compute the real maximal/minimal eigenvalue and
eigenvector. This is a crucial first step of the investigation since it
permits to asses the viability of this algorithm before being generalized.

The problem is then divided in two steps. The first step is to translate the
GEVP into a quadratic unconstrained binary optimization111One aims to get a
quadratic polynomial over binary variables. (QUBO), constructing the suitable
objective function (OF) in terms of binary variables. The second step, carried
out through the proprietary software, is to map the QUBO problem onto such
Ising-model, embodied by the qubit architecture (see Ref. [45] for the first
implementation of a QUBO problem on a D-wave QA and Refs. [25, 26, 27] for the
D-Wave Advantage system). The quantum annealear is then used to search for the
global minimum of the corresponding Ising Hamiltonian implemented on the
physical quantum processing units (QPUs).

The paper is organized as follows. In Sec. II, the physics problem from where
the non symmetric GEVP stemmed from is briefly illustrated. In Sec. III, the
OF, which our GEVP is mapped onto, and its expression in terms of the binary
variables are presented in detail. In Sec. IV, the results of our numerical
investigation on the annealing simulator, SA, as well as the D-wave Advantage
System QA are discussed at some length and compared. Finally, in Sec. V, our
conclusions are drawn with a focus on the perspectives. Appendices are devoted
to supply more details on the physics case and the construction of the OF.

##  II The physics problem

As previously stated, our aim is to solve the homogeneous Bethe-Salpeter
equation [33] of a system composed by two massive scalars bound through the
ladder exchange of a third massive scalar, directly in Minkowski space. A
variational method, tailored for a QA, is applied in order to obtain the
ground-state of the system. This non-perturbative quantum-field-theory problem
was solved in Ref. [37] by using the Nakanishi integral representation [46]
(NIR) of the Bethe-Salpeter (BS) amplitude, i.e. the key quantity for
describing the ground-state of the system and eventually calculating relevant
observables. Following Ref. [37], one can formally transform the hBSE in
Minkowski space into a GEVP adopting a standard discretization method. In
particular, a polynomial ortho-normal basis, composed by the Cartesian product
of Laguerre and Gegenbauer polynomials, was adopted. It should be pointed out
that solutions obtained directly in Minkowski space make a more direct
comparison between the theoretical outcomes and the experimental results
possible, without resorting to Euclidean solutions (see, e.g., Ref. [47, 48,
49] for an introduction to the Dyson-Schwinger plus hBSE framework, in
Euclidean space). Let us recall that Euclidean solutions could be directly
adopted for describing physical observables, once the assumptions of the
famous theorems by Osterwalder and Schrader [50, 51], stating necessary and
sufficient conditions for formally bridging Euclidean and Minkowskian quantum-
field theory, are fulfilled and the set of unavoidable approximations is
actually under control.

From a bird’s eye view, the two-body BS amplitude allows to reconstruct the
residue of the four-leg Green’s function at the bound-state mass pole.
Moreover, one can show that such an amplitude is the solution of a homogeneous
integral equation, i.e. the hBSE. Notice that, in relativistic quantum-field
theory, the hBSE has the same role in the bound-state description that the
Schrödinger equation has in non relativistic quantum-mechanics.

In conclusion, one is able to formally transform the Minkowskian hBSE into the
non symmetric GEVP of the form in Eq.(1). The eigenvalue is proportional to
the inverse of the square coupling constant present in the interaction vertex
(produced by the underlying field theory). Since the coupling constant is real
for physical reasons, real eigenvalues must be sought. It should be pointed
out that the eigenvalues and corresponding eigenvectors (eigenpair) in a non
symmetric GEVP can be real or complex. In case they are complex they appear in
conjugated pairs. Moreover, the eigenvectors ${\bf v}_{i}$ do not have the
property to be orthogonal. As specified above, in our case we are interested
only in the real eigenvalues, and particularly the one corresponding to the
largest, positive real $\lambda_{i}$, since it is the inverse of the minimal
coupling constant that allows the existence of a bound system with a given
mass [37].

##  III QA modeling of a non-symmetric generalized eigenvalue problem

In order to perform numerical calculations on a QA, the actual problem has to
be expressed as a QUBO model or, equivalently, to find an Ising-model form for
the problem under scrutiny (see, e.g., Ref. [52] for a tutorial on formulating
and using QUBO models).

For a symmetric GEVP, where the eigenvalues are real, there are already
several investigations with specific application to different problems (see,
e.g., Refs. [53, 54, 55, 29, 28, 30]). In all these applications, the
Rayleigh-Ritz variational principle is adopted to transform the eigenvalue
problem in a QUBO form. For the non symmetric GEVP, involving also complex
eigenpairs, dedicated efforts are still needed. Recalling that QUBO problems
are by definition strictly symmetric, one should be able to replace the non
symmetric original problem with a symmetric one, which allows to evaluate the
set of real eigenpairs we are interested in. This will be achieved by
transforming the initial GEVP into a symmetric one with an equivalent
spectrum. Thus, the same algorithm used for applying the Rayleigh-Ritz
variational principle, with minimal changes, can be adopted. In performing
such a transformation, the minimization problem acquires a quadratic
dependence on the parameter playing the role of eigenvalue. This means that
each eigenvalue is a minimum. This on one side can be usefully exploited for
addressing the entire spectrum on the other side require to guide the search
toward the relevant region where the eigenvalue has to be found.

In this Section we go through the formalism we adopted for the non-symmetric
case with some detail.

###  III.1 Formalism

Since the physical case under consideration generates a non-symmetric matrix
$A$ in Eq. (1), the Rayleigh-Ritz variational principle is not directly
exploitable. A first trivial attempt to overcome this problem might be to
single out the symmetric part of the matrix $A$, and solve the corresponding
symmetric GEVP, recalling that $B$ by itself is symmetric. Unfortunately, this
formal manipulation does not return the original global minimum. This is
illustrated in Appendix C, where an inequality between the eigenvalues of a
non-symmetric matrix and the ones of its symmetric part is discussed in the
simpler case of a standard eigenvalue problem. Therefore, the eigenvalue
problem $\Bigl{[}(A+A^{T})/2\Bigr{]}~{}{\bf v}_{i}=\lambda_{i}~{}B~{}{\bf
v}_{i}$ cannot be used for determining the minimal/maximal eigenvalue of Eq.
(1). The same issue is met also for other methods of symmetrization, like
$A^{T}A$ (see Appendix C for more discussion and examples).

A possible way to overcome this problem was discussed in Ref. [38], where the
symmetrization of the original problem is recast as the product of
$A-\tilde{\lambda}~{}B$ with its transpose, i.e.

$\displaystyle{S}=\Bigl{[}A^{T}-\tilde{\lambda}~{}B^{T}\Bigr{]}\Bigl{[}A-\tilde{\lambda}~{}B\Bigr{]}\,,$
(2)

where $\tilde{\lambda}$ can be complex (see below).

The resulting OF, which now contains a symmetric product, reads

$\displaystyle f(A,B,{\bf v},\tilde{\lambda})={\bf
v}^{T}~{}\Bigl{[}A^{T}-\tilde{\lambda}~{}B^{T}\Bigr{]}\Bigl{[}A-\tilde{\lambda}~{}B\Bigr{]}~{}{\bf
v}=$ $\displaystyle={\bf
v}^{T}~{}\Bigl{[}A^{T}A-2\tilde{\lambda}M+\tilde{\lambda}^{2}B^{T}B\Bigr{]}~{}{\bf
v}\geq 0~{}~{},$ (3)

where ${\bf v}$ is chosen to be a real normalized vector and $M$ is a
symmetric matrix, defined as:

$\displaystyle M={(B^{T}A+A^{T}B)\over 2}~{}~{}.$ (4)

In what follows, it should be kept in mind that our goal is to find a real
eigenpair of the analyzed GEVP, although the OF $f(A,B,{\bf
v},\tilde{\lambda})$ vanishes also for complex eigenpairs.

The striking difference between the global-minimum search of $f(A,B,{\bf
v},\tilde{\lambda})$ and the standard eigenvalue problem is the quadratic
dependence upon the variable $\tilde{\lambda}$. Calling $\lambda$ the solution
of the quadratic form in Eq. (3) for a given real vector ${\bf v}$, one gets

$\displaystyle\lambda({\bf v})=\lambda^{R}({\bf v})\pm i\lambda^{I}({\bf v})$
(5)

where the real part is

$\displaystyle\lambda^{R}({\bf v})={\bar{M}({\bf v})\over|B~{}{\bf
v}|^{2}}~{}.$ (6)

with $\bar{M}({\bf v})={\bf v}^{T}~{}M~{}{\bf v}$. The imaginary part is
written as follows

$\displaystyle\lambda^{I}({\bf v})={\sqrt{{\cal I}({\bf v})}~{}\over|B~{}{\bf
v}|^{2}}$ (7)

where

$\displaystyle{\cal I}({\bf v})=||A~{}{\bf v}||^{2}~{}||B~{}{\bf
v}||^{2}-\bar{M}^{2}({\bf v})~{}~{},$ (8)

with $||A~{}{\bf v}||^{2}={\bf v}^{T}~{}A^{T}A~{}{\bf v}$ and equivalently for
$||B~{}{\bf v}||^{2}$. Notice that ${\cal I}({\bf v})$ is a positive quantity,
since

$\displaystyle{\cal I}({\bf v})=||A~{}{\bf v}||^{2}~{}||B~{}{\bf
v}||^{2}-\left[{\bf v}^{T}~{}{B^{T}A+A^{T}B\over 2}{\bf v}\right]^{2}$
$\displaystyle\geq||A~{}{\bf v}||^{2}~{}||B~{}{\bf v}||^{2}-\left[{||{\bf
v}^{T}B^{T}A{\bf v}||+||{\bf v}^{T}A^{T}B{\bf v}||\over 2}\right]^{2}$
$\displaystyle\geq||A~{}{\bf v}||^{2}~{}||B~{}{\bf v}||^{2}-||{\bf
v}^{T}B^{T}A{\bf v}||^{2}\geq 0~{}.$ (9)

Therefore, our task is to set up an algorithm for finding a real vector ${\bf
v}$ that minimizes both the OF and the imaginary part of the eigenvalue,
$\lambda^{I}({\bf v})$, so that the set $\\{{\bf v},\lambda^{R}({\bf v})\\}$
is the searched eigenpair. Notice that a real vector ${\bf v}$ that minimizes
${\cal I}({\bf v})$, also minimizes the OF $f(A,B,{\bf
v},\tilde{\lambda}=\lambda^{R})$. However, the opposite is not true, since the
OF vanishes also for complex eigenpairs.

In Ref. [38], the search was carried out by classically minimizing ${\cal
I}({\bf v})$ through a gradient-descent method. Unfortunately, ${\cal I}({\bf
v})$ has a quartic dependence on ${\bf v}$, and therefore the corresponding
minimization cannot straightforwardly be translated into a QUBO problem. The
quartic problem can be recast into a quadratic one by adding suitable
penalties, at the price of a substantial increase of the dimension of the QUBO
problem (see, e.g., Ref. [56, 52]). However this will require to substantially
increase the matrix size to pass to the annealer. In order to overcome the
difficulties, we approached the problem in a different way.

First of all, given the exploratory nature of our investigation, instead of
directly using the OF in Eq. (3), we introduced a hybrid algorithm, namely, we
transformed the initial GEVP into a standard one, leaving the generalized
problem for future investigations. We exploit the symmetry and non singularity
of the matrix $B$, decomposing it by using the standard $LDL^{T}$
factorization, with $L$ a lower triangular matrix, having diagonal elements
equal to $1$, and $D$ a non singular diagonal matrix (see, e.g., the LAPACK
library). Hence, the original GEVP, Eq. (1), becomes

$\displaystyle C~{}L^{T}{\bf v}_{i}=\lambda_{i}~{}L^{T}{\bf v}_{i}~{},$ (10)

where

$\displaystyle C=D^{-1}L^{-1}~{}A~{}[L^{T}]^{-1}~{}~{},$ (11)

is non symmetric, and the real vector $L^{T}{\bf v}_{i}$ is no longer
normalized. Since the global-minimum search is not affected by multiplicative
factors, we can use a normalized vector formally defined by

$\displaystyle{\bf w}=L^{T}{\bf v}/||L^{T}{\bf v}||~{}.$ (12)

Then, one can address the QUBO problem by replacing the initial OF with an
equivalent one, given by

$\displaystyle f(C,{\bf w},\tilde{\lambda})={\bf
w}^{T}~{}\Bigl{[}C^{T}-\tilde{\lambda}\Bigr{]}\Bigl{[}C-\tilde{\lambda}\Bigr{]}~{}{\bf
w}~{}.$ (13)

The real part of the solution of the quadratic form in $\tilde{\lambda}$, for
a given ${\bf w}$, is

$\displaystyle\lambda^{R}({\bf w})={1\over|{\bf w}|^{2}}~{}{\bf
w}^{T}{C^{T}+C\over 2}{\bf w}~{},$ (14)

while the imaginary part reads

$\displaystyle\lambda^{I}({\bf w})=\sqrt{{|C{\bf w}|^{2}\over|{\bf
w}|^{2}}-\Bigl{[}\lambda^{R}({\bf w})\Bigr{]}^{2}}~{}~{}.$ (15)

Since we have chosen ${\bf w}^{T}{\bf w}=1$, it is understood that in order to
get the normalized eigenvector of the original GEVP, one has to apply the
inverse transformation to ${\bf w}$ and normalize the result, i.e.

$\displaystyle{\bf v}=[L^{T}]^{-1}~{}{\bf w}/||[L^{T}]^{-1}~{}{\bf w}||~{}.$
(16)

###  III.2 Algorithm implementation

The OF in Eq. (13) is explicitly written as follows:

$\displaystyle f(C,{\bf w},\tilde{\lambda})={\bf w}^{T}{\cal
S}(\tilde{\lambda}){\bf w}=\sum_{i,j}w_{i}w_{j}{\cal S}_{i,j}\,,$ (17)

where ${\cal S}$ is an $n\times n$ matrix ($n$ is dictated by the
discretization method applied to the hBSE) given by

$\displaystyle{\cal
S}_{i,j}(\tilde{\lambda})=\tilde{\lambda}^{2}\delta_{i,j}-\tilde{\lambda}(C_{i,j}+C_{j,i})+\sum_{\ell}C_{\ell,i}C_{\ell,j}\,,$
(18)

and ${\bf w}\equiv\\{w_{i}\\}$ is a normalized vector. By using the binary
representation introduced in Appendix B, the components $w_{i}\in[-1,1]$ are
expressed in terms of the binary basis as follows

$\displaystyle
w_{i}=-q_{i,b}+\sum_{\ell=1}^{b-1}\frac{q_{i,\ell}}{2^{\ell}}~{},$ (19)

where $b$ is the number of bits, $q_{i,\ell}=0,1$ and the bit $q_{i,b}$
carries the sign, i.e. for $w_{i}\geq 0$ ($w_{i}<0$) one assigns $q_{i,b}=0$
($q_{i,b}=1$). Notice that the binary expression in Eq. (19) belongs to the
interval $[-1,1)$.

The string $\\{-1,1/2,1/4,~{}\dots~{},1/2^{b-1}\\}$ is the so-called precision
vector, ${\bf p}$, since the last term $1/2^{b-1}$ controls the precision at
which $w_{i}\in[-1,1)$ is approximated. Within our notation ${\bf p}$ is a
column vector with dimension $b$. Using the precision vector ${\bf p}$, it is
possible to construct a rectangular matrix, $P^{T}$, with $n$ rows and
$n\times b$ columns, that allows to transform a real vector ${\bf w}$ into its
binary expression, and eventually $f(C,{\bf w},\tilde{\lambda})$ in a form
suitable for a QUBO evaluation by using a QA. In a compact form, one writes

$\displaystyle{\bf w}=P^{T}~{}{\bf x}~{},$ (20)

where

$\displaystyle P^{T}=\begin{pmatrix}~{}{\bf p}^{T}&0&\dots&0\\\ 0&~{}{\bf
p}^{T}&\dots&0\\\ \vdots&\vdots&\ddots&0\\\ 0&0&\dots&{\bf
p}^{T}\end{pmatrix}~{},$ (21)

and

$\displaystyle{\bf x}=\begin{pmatrix}{\bf q}_{1;b}\\\ {\bf q}_{2;b}\\\
\vdots\\\ {\bf q}_{n;b}\end{pmatrix}~{},$ (22)

with ${\bf q}_{i;b}\equiv\\{q_{i,1},q_{i,2},\dots,q_{i,b}\\}$. Inserting Eq.
(20) in Eq. (17) one obtains (see Appendix B for details)

$\displaystyle f(C,{\bf w},\tilde{\lambda})={\bf x}^{T}~{}P{\cal
S}(\tilde{\lambda})P^{T}~{}{\bf x}~{}.$ (23)

One easily recognizes that ${\bf x}$ is a column vector containing $n\times b$
bits, leading to the expected QUBO form of our problem, and $P{\cal
S}(\tilde{\lambda})P^{T}$ is a square symmetric matrix, with dimensions
$(n\times b)~{}\times~{}(n\times b)$.

Starting with some given real value of $\tilde{\lambda}$ (see the next
subsection), that guides the algorithm towards a real-valued eigenpair, the SA
or the QA returns a binary vector ${\bf x}$ for each annealing cycle $N_{A}$.
Among the $N_{A}$ vectors ${\bf x}$ we select the one that gives the minimal
value of $f(C,{\bf w},\tilde{\lambda})$. Such a vector allows to get ${\bf
w}$, necessary for calculating $\lambda^{R}({\bf w})$ and $\lambda^{I}({\bf
w})$, using Eqs. (14) and (15), respectively. To accurately determine ${\bf
w}$, which minimizes at the same time the number of bits required in the QA, a
two-step search is implemented: the first step is the guess phase (GP), and
the second one is the iterative gradient-descent phase (DP) (see, e.g., Refs.
[28, 29] and [30] for the symmetric case).

###  III.3 Initial guess for the non-symmetric case

Differently from what happens for the symmetric case [28, 29, 30, 32], the
non-symmetric one leads to a non linear dependence on the scalar parameter
$\tilde{\lambda}$, and therefore the search needs a particular care. It should
be emphasized that the OF is minimized by the entire set of real eigenpairs.
Hence, in order to point to the largest, positive real eigenvalue, we need to
add further information to our algorithm. We can guide the search of the
minimum by using the Gershgorin circle theorem (see, e.g., Ref. [57, 58], for
an introduction). Noteworthy, it provides bounds for each eigenvalue of a
given complex matrix, allowing the search range to be narrowed.

The theorem yields

$\displaystyle|\lambda-a_{ii}|=\sum_{j\neq i}\left|{a_{ij}x_{j}(\lambda)\over
x_{i}(\lambda)}\right|\leq\sum_{j\neq i}\left|a_{ij}\right|=R_{G}(a_{ii})~{},$
(24)

where $\lambda$ is an eigenvalue of the given matrix, with elements $a_{ij}$,
and $x_{i(j)}(\lambda)$ are components of the corresponding eigenvector, such
that $|x_{j}/x_{i}|\leq 1$, $\forall~{}j\neq i$. Hence, an eigenvalue belongs
to a proper disc, with radius $R_{G}(a_{ii})$ and center $a_{ii}$. The same is
valid by summing over the columns, once the transpose is considered. It is
also worth mentioning a second Gershgorin theorem [58]. It proves that if a
disc is disconnected from all the others, then it contains one and only one
eigenvalue, which necessarily is real (the complex ones are conjugated).

Inspired by the Gershgorin theorem, a suitable permutation can be applied in
order to rearrange the matrix $C$ in Eq. (11), so that the diagonal elements
are in decreasing order, i.e. $c_{11}\geq c_{22}\geq\dots\geq c_{(n\times
b)(n\times b)}$. The initial value of $\tilde{\lambda}$ is chosen equal to
$c_{11}$. The matrix $P{\cal S}(\tilde{\lambda})P^{T}$ is passed to the QA or
the SA. For each of the $N^{GP}_{A}$ annealing cycles the binary vectors ${\bf
x}_{\alpha}$, and the corresponding values of the OF is returned. With
$\alpha$ we label the annealing cycle. Following the Gershgorin constraint, we
analyze the set

$\\{{\bf x}_{\alpha},f(C,{\bf
w}_{\alpha},\tilde{\lambda}=c_{11}),\lambda^{R}({\bf
w}_{\alpha}),\lambda^{I}({\bf w}_{\alpha})\\}~{},$

and eliminate the solutions such that

$|\lambda^{R}({\bf w}_{\alpha})-c_{11}|>R_{G}(c_{11})~{}.$

Among the surviving solutions, the one that satisfies the condition

$\displaystyle{f}_{best}^{GP}=\min_{{\bf w}_{\alpha}}f(C,{\bf
w}_{\alpha},\tilde{\lambda}=\lambda^{R}({\bf w}_{\alpha}))$ (25)

is retained, and $\lambda^{R}(\bf{w}_{\alpha})$ is calculated. The
minimization search is repeated three times, $i=1,2,3$, testing if the
condition ${f}^{GP}_{best;i}<{f}^{GP}_{best;i-1}$ is satisfied. If there is an
improvement, the next run receives
$\tilde{\lambda}=\lambda^{R}(\bf{w}_{\alpha_{best;i}})$, otherwise the
iteration stops and the GP is closed. Empirically, we found this procedure is
more effective in finding a good initial eigenvector than just increasing the
number of annealing samples in a single iteration. After completing the
minimization search, the best pair, i.e. $\lambda^{\bf R}(\bf{w}_{{GP}})$ and
the corresponding real vector ${\bf w}_{{GP}}$, that satisfies the condition
in Eq. (25), is passed to the gradient-descent phase. The algorithm of the GP
is schematically illustrated in Tab. 1, and the numerical results that
illustrate the behavior of the algorithm are discussed in Section IV.1.

Table 1:  The guess-phase search of the minimum for the OF in Eq. (17), obtained from Eq. (3) after applying a classical $LDL^{T}$ factorization. An assigned number of bits $b$ is used for the transformation to the QUBO form, Eq. (23) (see text).  | Algorithm for the guess-phase, with a given number of bits $b$
---|---
1: |  Assign the input value $\tilde{\lambda}=c_{11}$, suggested by the Gershgorin theorem
2: |  While ${f}^{GP}_{best;i}<{f}^{GP}_{best;i-1}$, for $i=1,2,3$
3: |  Minimize the QUBO form, Eq. (23), with $N_{A}^{GP}$ annealing cycles
4: |  Look for acceptable $\lambda^{R}({\bf w}_{\alpha_{i}})$ falling in the Gershgorin disc
5: |  Among the acceptable solutions, check if ${f}^{GP}_{best;i}<{f}^{GP}_{best;i-1}$ and continue the loop, replacing $\tilde{\lambda}$ by $\lambda^{R}({\bf w}_{\alpha_{best;i}})$
6: | Put $\lambda^{R}({\bf w}_{{GP}})=\lambda^{R}({\bf w}_{\alpha_{best;i}})$ and ${\bf w}_{{GP}}={\bf w}_{\alpha_{best;i}}$ and pass those quantity to the gradient-descent phase

###  III.4 Gradient-descent phase

The algorithm proceeds through a second phase: the gradient-descent phase,
which is based on the Hessian of the OF in Eq. (17). Let us call
${\bm{\delta}}(z)\equiv\\{\delta_{i}(z)\\}$ the variation of the vector ${\bf
w}^{z}$, i.e. ${\bm{\delta}}(z)={\bf w}^{z}-{\bf w}^{z-1}$, where $z$ labels
the iterations or zoom steps of the gradient-descent phase. Note that ${\bf
w}^{0}={\bf w}_{{GP}}$. Then, the following binary expression for the
${\bm{\delta}}(z)$ components holds

$\displaystyle\delta_{i}(z)={1\over
2^{z}}~{}\Biggl{(}-q_{i,b}+\sum_{\ell=1}^{b-1}{q_{\ell,b}\over
2^{\ell}}\Biggr{)}={1\over 2^{z}}~{}{\bf p}^{T}\cdot{\bf q}_{i;b}~{},$ (26)

where the prefactor $1/2^{z}$ decreases at each iteration by increasing $z$
(cf. Refs. [28, 29] and Ref. [30]). The prefactor $1/2^{z}$ allows to control
the refinement of our search, increasing the precision of the gradient-descent
method. It should be pointed out that by increasing the number of iterations,
i.e. zoom steps, one could decrease the initial number of bits $b$ in the
expansion of the real vectors involved in the minimization, without worsening
the accuracy of the algorithm (see also Refs. [28, 29] for the symmetric
GEVP).

By expanding $f(C,{\bf w}^{z},\tilde{\lambda})$ around ${\bf w}^{z-1}$ up to
the second order, one gets (cf. Eqs. (17) and (18))

$\displaystyle\Delta f(C,{\bf
w}^{z},\tilde{\lambda})=2{\bm{\delta}}^{T}(z){\cal S}(\tilde{\lambda}){\bf
w}^{z-1}+{\bm{\delta}}^{T}(z){\cal S}(\tilde{\lambda}){\bm{\delta}}(z)~{}.$
(27)

The term linear in ${\bm{\delta}}(z)$ can be recast in a QUBO form by i)
recalling that $q_{\ell}=q^{2}_{\ell}$, being $q_{\ell}=0,1$, and ii)
transforming the vector ${\cal S}(\tilde{\lambda}){\bf w}^{z-1}$ in a
$(n\times b)\times(n\times b)$ diagonal matrix, where the vector components
are the diagonal elements (see Ref. [28]). Then, the OF suitable for the
gradient-descent phase is obtained. Specifically, recalling that one can
divide or multiply by a constant an OF without affecting the minimization
procedure, one has

$\displaystyle\widehat{f}(C,{\bf\delta}(z),\tilde{\lambda})={\bm{\delta}}^{T}(z)~{}{\cal
Q}(z,\tilde{\lambda})~{}{\bm{\delta}}(z)$ (28)

where the square matrix ${\cal Q}$ is given by

$\displaystyle\Bigl{[}{\cal
Q}(z,\tilde{\lambda})\Bigr{]}_{ij}=2\delta_{ij}~{}\Bigl{[}{\cal
S}(\tilde{\lambda})~{}{\bf w}^{z-1}\Bigr{]}_{i}+{\cal
S}(\tilde{\lambda})_{ij}~{},$ (29)

with the elements of the symmetric matrix ${\cal S}(\tilde{\lambda})$ written
in Eq. (18). The transformation to the binary expression of the vectors
${\bm{\delta}}(z)$ follows the rule given in Eq. (20).

The gradient-descent phase proceeds through iteration on $z$, with
$z=1,\dots,N_{z}$, and aims to minimize the OF in Eq. (28), starting with the
value of $\tilde{\lambda}=\lambda^{R}({\bf w}_{{GP}})$ and the corresponding
vector ${\bf w}_{{GP}}$. At each zoom step $z$, an inner loop (labelled with
$i$) is opened. Within this loop, the QA, or the SA, returns an ensemble of
$\alpha_{i}=1,2,\dots,N^{DP}_{A}$ qubits states, from which we select the one
with minimal energy, i.e.

$\displaystyle\widehat{f}^{DP}_{best;i}(z)=\min_{{\bm{\delta}}_{\alpha_{i}}(z)}\widehat{f}(C,{\bm{\delta}}_{\alpha_{i}}(z),\tilde{\lambda}=\lambda^{R}({\bf
w}^{z}_{\alpha_{i}}))\,,$ (30)

where $i$ is the inner-loop index for a fixed $z$ and ${\bf
w}^{z}_{\alpha_{i}}={\bm{\delta}}_{\alpha_{i}}(z)+{\bf w}^{z}_{best,{i-1}}$.
For the value $\alpha_{i}=\alpha_{best}$ in correspondence of
$\widehat{f}^{DP}_{best;i}(z)$, we calculate the values $\lambda^{R(I)}({\bf
w}^{z}_{\alpha_{best}};i)$, that is input for the next iteration starting with
the new eigenvector ${\bf w}^{z}_{best,i}={\bf w}^{z}_{\alpha_{best},i}$.
When, for a fixed $z$, the condition
$\widehat{f}^{DP}_{best;i+1}(z)\geq\widehat{f}^{DP}_{best;i}(z)$ is found,
then the internal loop on $i$ is terminated and the value of $z$ updated,
namely $z\rightarrow z+1$. The final best result obtained for a given $z$ is
used as starting point in the $z+1$ iteration. The procedure is repeated until
the desired precision is reached ($z^{max}$).

Finally, one obtains the eigenvector of the GEVP in Eq. (1) by inverting the
relation in Eq. (12), i.e.

$\displaystyle{\bf v}_{best}=\Bigl{[}L^{T}\Bigr{]}^{-1}{\bf
w}_{best}^{z_{max}}/\left|\left|\Bigl{[}L^{T}\Bigr{]}^{-1}{\bf
w}_{best}^{z_{max}}\right|\right|~{},$ (31)

where $z_{max}$ indicates the final run of the iterative gradient-descent
method and ${\bf w}^{z_{max}}_{best}$ is the eigenvector corresponding to the
final eigenvalue

$\displaystyle\lambda_{best}=\lambda^{R}({\bf w}_{best}^{z})\,.$ (32)

It should be recalled that the Hessian is always positive, so that a
minimization path is ensured. The gradient-descent algorithm is sketched in
Tab. 2.

Table 2: The gradient-descent phase search of the global minimum on both SA and QA. The gradient-descent phase is based on the Hessian of the OF in Eq. (3) (see text). The same guess-phase number of bits $b$ is used.  | Algorithm for the gradient-descent phase
---|---
1: |  Use $\tilde{\lambda}=\lambda^{R}({\bf w}_{{GP}})$, from the guess phase
2: |  While $1/2^{z}\geq\epsilon^{\prime}$ ($z\geq 1$)
3: |  Transform the obj. function in Eq. (28) to QUBO form with $1/2^{z}$ in ${\bf p}$
4: |  While $\widehat{f}_{best;z}^{DP}(i)<\widehat{f}_{best;z}^{DP}(i-1)$ for i =1,2,3,…
5: |  Minimize the QUBO form in Eq. (28), with $N^{DP}_{A}$ annealing cycles
6: |  Search the minimal value among the $N^{DP}_{A}$ energies returned by the annealer ($\alpha_{best}$)
7: |  Replace $\tilde{\lambda}$ with $\lambda^{R}_{z}({\bf{w}}_{\alpha_{best}};i)$
8: |  Pass to $z+1$ with $\tilde{\lambda}=\lambda^{R}_{z}({\bf w}_{\alpha_{best}};N_{z})$

By repeating the entire minimization procedures, i.e guess plus gradient-
descendent phases, given the non-deterministic nature of the annealing
process, errors and fluctuations of the final results, $\\{\lambda_{best},{\bf
v}_{best}\\}$, are generated. In order to study such uncertainties, the entire
algorithm was run by $N_{run}$ times, keeping fixed the initial matrices $A$
and $B$. This provides an estimate of the uncertainties associated to the
annealing process, and how these uncertainties propagate through the
algorithm. We performed this statistical analysis using both the noisy
annealer simulator, provided by D-Wave Systems software-package, with large
$N_{run}$, and the QA with a substantially smaller $N_{run}$, given the
limitation on running-time at our disposal. In what follows, the results
obtained with the simulator are indicated by inserting SA at the beginning of
the figure captions, while the results obtained with the D-wave Advantage
System 4.1 are labelled by QA.

##  IV Results

In order to study the reliability of the algorithm, we used the matrices $A$
and $B$ obtained from the discretization of the hBSE describing a system of
two massive scalars, with mass equal to $m$ and strong binding energy
$E_{B}/m=1.0$. They interact through the ladder-exchange of a massive scalar,
with mass $\mu/m=0.15$ (see Eq. (37) in Appendix A, and for more details Ref.
[37]). In our investigation we classically evaluate the $LDL^{T}$
factorization, exploiting the non singularity of the symmetric matrix $B$, and
then we solve the standard eigenvalue problem given in Eq. (13). We leave to
future works the generalized case. Let us recall that we focus only on the
search of the largest real eigenpair. It should be recalled that, in our
approach, the largest real eigenvalue corresponds to the lowest coupling
constant, which is able to bind the two massive scalars. Finally, we give a
first insight on the scalability of the proposed algorithm by considering
matrix dimensions $n_{M}=4,8,12,16,24,32$.

In an effort to better optimize the available quantum-computing resources, a
challenging issue is to guide the choice of the input parameters of the
proposed algorithm. They are given by: i) the number of bits $b$, needed for
the binary representation of the real vectors involved in the problem, and ii)
the number of annealing cycles, i.e. $N^{GP}_{A}$, for the guess phase, and
$N^{DP}_{A}$, for the gradient-descent one. This important optimization task,
which is time consuming, was carried out by using the SA, and the
corresponding results are discussed in the first part of this Section. In the
second part, we illustrate the outcomes obtained by using the QA.

In particular, to perform the SA analysis, a matrix with dimension $n_{M}=32$
was adopted and the entire algorithm was run for $N_{run}=500$ times, in order
to gather enough statistics and determine a reasonable set of input parameters
to be used with the QA. For the remaining part of this work, the following
symbols are adopted: i) $\bar{\lambda}_{best}(\bar{\lambda}_{best}^{I})$ is
the average of $\lambda^{R}(\lambda^{I})$ obtained from Eq. (14) (Eq. (15))
over the $N_{run}$ runs of the algorithm for $z^{max}$; ii) $\overline{||{\bf
v}_{true}-{\bf v}_{best}||}$ indicates the mean Euclidean distance, where
${\bf v}_{true}$ is the real eigenvector of the original GEVP, and ${\bf
v}_{best}$ computed for $z^{max}$ is given in Eq. (31).

###  IV.1  Running the guess-phase on the SA

Selecting a reliable initial guess in the first phase of the algorithm,
sketched in Tab. 1, is a crucial step for obtaining an accurate result in the
subsequent gradient-descent phase. As already mentioned, all the real-valued
eigenpairs are acceptable for the annealer, and therefore one has to guide the
search toward the target solution, i.e. the real eigenpair with the largest
real eigenvalue. In addition, one has to carefully consider another source of
numerical challenges: the (necessarily) finite number of bits $b$. In fact,
one has to assign $b$ for representing the components of the normalized vector
entering the definition of the QUBO OF in Eq. (23), and eventually construct
the matrix $P{\cal S}P^{T}$. This obvious limitation can generate spurious
minima, and direct the path toward a wrong target. An example is shown in Fig.
1, where the case of the three largest real eigenvalues,
$\lambda_{1}>\lambda_{2}>\lambda_{3}$, is studied. For bit numbers $b=3,4,6,8$
and running the entire algorithm only one time (i.e. $N_{run}=1$), the values
of $\lambda^{R}({\bf w}_{\alpha})$, Eq. (14), obtained in the guess-phase from
each of the $N^{GP}_{A}=2000$ annealing cycles on the SA, are plotted as a
function of the corresponding annealer energies divided by $[\lambda^{R}]^{2}$
(this is suggested by the quadratic dependence of the OF on $\tilde{\lambda}$
and the purpose of assigning a normalization to different outcomes). Notice
that the annealer energies correspond to the values of Eq. (23) as returned by
the annealer in the QUBO form i.e. with the binary vectors not normalized.

Figure 1: (Color online) [SA] Scatter plot of the $N_{A}^{GP}=2000$ guess-
phase values $\lambda^{R}_{\alpha}$, Eq.(14), vs. the corresponding annealer
energies divided by $[\lambda^{R}]^{2}$, for the bit numbers $b=3,4,6,8$ used
to get the QUBO matrix $P{\cal S}P^{T}$, Eq. (23), with dimension $32\times
b$. Dashed lines: the three largest real eigenvalues,
$\lambda_{1}>\lambda_{2}>\lambda_{3}$ of the matrix $C$, Eq. (11). Solid line:
$|c_{11}-R_{G}(c_{11})|$ with $R_{G}(c_{11})=\sum_{j\neq 1}|c_{1j}|$, i.e. the
lower limit of the Gershgorin disc, with center $c_{11}=0.190$.

In Fig. 1, the dashed horizontal lines represent the three largest real
eigenvalues (calculated classically) of the $32\times 32$ matrix $C$ in Eq.
(11). The solid line is the lower extremum of the Gershgorin disc, with center
$c_{11}=0.190$. Such an extremum is the lowest value admissible for
$\lambda_{1}$. In the various panel, one observes several regions where the
results appear to clusterize despite the initial guess
$\tilde{\lambda}=c_{11}$ is very close to the actual eigenvalue
$\lambda_{1}=0.188225$. Once a more detailed analysis is performed, it can be
shown that the clusterization is present already at $b=3$, although it is not
observable in the figure. Moreover, some clusters approach the minima
generated by the other two eigenpairs, with eigenvalues $\lambda_{2}$ and
$\lambda_{3}$, while others are spurious minima, that are generated by the
limited precision of the assigned binary representation, i.e. $3$ bits. By
increasing $b$, the clusters migrate close to the actual eigenvalues, with
smaller and smaller energies (notice the different $x$-axis scale for the four
panels). At the same time, a larger amount of points are approaching
$\lambda_{1}$, given the increasing precision, although they spread over a
wider range of energy as expected. This is due to the increase of states
degeneration when the number of $b$ increases, since the returned bit strings
are neither normalized nor orthogonal.

Interestingly, $\lambda_{1}$ is approached from below by increasing $b$. This
behavior can be understood considering the bias due to the spectrum of the
matrix under scrutiny. In fact, the decreasing ordering of the real
eigenvalues, that follows from the chosen ordering of the matrix elements
along the main diagonal, can generate the observed pattern, since the
annealing process should clearly spend a non negligible amount of time to
explore the minima corresponding to the part of the spectrum below
$\lambda_{1}$. It turns out that the feature persists, even after running many
times the code.

Increasing $b$ allows to refine the solution found by the annealer. On the
other hand, a larger value of $b$ implies a larger dimension of the matrix
representing the OF. Unfortunately, such matrix dimension is constrained by
the actual topology of the qubits network during the process of minor-
embedding (see below). This limits the precision of our solution in the
initial guess phase. Fortunately, as also discussed in Refs. [28, 29], a
reasonable compromise can be reached by properly increasing the number of
iterated annealing cycles and adopting a small number of bits, as investigated
in the following gradient-descent phase.

In order to gain insights on the mapping of the OF in Eq. (13) and the one
given in terms of the binary strings, it is useful to show an analysis
equivalent to the one in Figs. 1, but focusing on the values of the OF in Eq.
(13), where normalized vectors are used. In particular, the cases with $b=3$
and $b=8$ are shown in Fig. 2, where the squares are the solutions selected
after applying the Gershgorin criterium with $R_{G}(c_{11})$. The
clusterization is more evident than in Fig. 1 and presents a different overall
behavior compared to the results directly obtained from the annealing process.
In particular, the accuracy of the OF mapping onto the binary basis improves
by increasing the number of bits, and the annealer spend more and more time
close to the actual minima.

Summarizing, two main consequences can be deduced from the calculations shown
in Figs. 1 and 2. First, in the guess phase, it is not possible to select the
final value of $\lambda^{R}$ only on the basis of the annealer energies, as in
the symmetric case (see, e.g., Refs. [28, 29]), but we need the Gershgorin
criterium, to overcome the issue of the quadratic dependence on
$\tilde{\lambda}$ in the OF. Second, we note that by exploiting the mentioned
quadratic dependence, the annealer becomes a powerful tool to explore a wide
portion of the spectrum of our eigenvalue problem. In principle, it is
possible to identify also the other eigenpairs, by studying the clusterization
of the returned energies. This appealing feature will be investigated
elsewhere.

Figure 2:  (Color online). [SA] The same as the panels $b=3$ and $b=8$, in
Fig. 1, but with the values of $f(C,{\bf
w}_{\alpha},\lambda^{R})/(\lambda^{R})^{2}$, Eq. (13), on the $x$-axis,
instead of the annealing energies (see Eq. (23)). Recall that here the vector
${\bf w}_{\alpha}$ is normalized. Red squares are the solutions selected after
applying the $R_{G}(c_{11})$ cut, and blue crosses are the solutions outside
the Gershgorin disc.

In Fig. 3, the GP success-rate, i.e. the number of times the algorithm is able
to find a solution falling inside the Gershgorin disc, is shown as a function
of the GP annealing cycles, $N_{A}^{GP}$, for $b=3$ and running the code
$N_{run}=500$ times for the $32\times 32$ matrix. A number of $\sim 200$
annealing cycles seems to guarantee a $100\%$ probability of having at least
one point within the Gershgorin disc, by using the SA.

Figure 3:  [SA] Guess-phase success-rate vs. $N_{A}^{GP}$, for a $32\times 32$
matrix in Eq. (13). The solid line was obtained considering $N_{run}=500$
independent samples with $b=3$ (see Tab. 1 for the outline of the algorithm in
the guess-phase).

Since it is crucial to receive as reliable input as possible from the guess
phase, so that a better final outcome could be reached, we attempted to find a
tighter bound, if any, for $\tilde{\lambda}$ to be used in the initial run of
the gradient-descent phase. In order to implement this, we arbitrarily reduced
the Gershgorin radius, $R_{G}(c_{11})$, to a chosen value $r$. In Fig. 4, the
GP success-rate for the above $32\times 32$ matrix is analyzed by i) taking
$b=3$ and ii) running $N_{run}=500$ times the GP algorithm with annealing
cycles $N_{A}^{GP}=200$. As expected, the GP success-rate decreases by
reducing $r$, and we were not able to find any solution below $r/R_{G}<0.4$.
In conclusion, for the matrices considered in this work, the most reasonable
and practical choice is $r=R_{G}(c_{11})$.

Figure 4: [SA] Guess-phase success-rate vs. an arbitrary reduction of the
Gershgorin radius (represented by the ratio $r/R_{G}$), for a $32\times 32$
matrix in Eq. (13). The solid line was obtained by considering $N_{run}=500$
independent samples, with $b=3$ and $N_{A}^{GP}=200$ (see Tab. 1 for the
outline of the algorithm in the guess phase).

###  IV.2 Running the gradient-descent phase on the SA

The main goal of the gradient-descent phase is the improvement of the accuracy
in determining the searched eigenpair. We pursued this scope by decreasing the
prefactor $1/2^{z}$, with $z\geq 1$, in Eq. (26). In this subsection we focus
on the impact of the input parameters by increasing this prefactor. Notice
that for all the analysis carried out in this section, we run $N_{run}=500$
independent simulations.

In Fig. 5, the relative precision
$|1-\overline{\lambda}_{best}(z)/\lambda_{true}|$ (cf. Eq. (32)) dependence on
the number of bits, $b\geq 2$, is presented for various zoom steps $z$ (in
what follows, $z=0$ indicates the guess phase and the values $z\geq 1$ belong
to the gradient-descent phase). The dots in the figure are the mean, over
$500$ independent runs, of the relative precision obtained by using
$N_{A}^{GP}=200$, and $N_{A}^{DP}=20$. An overall decreasing pattern can be
observed for increasing $b$ and fixed $z$, that can be largely ascribed to a
better identification of the minimum during the guess phase, $z=0$ (see
Section IV.1). Interestingly, the numerical outcomes shown in Fig. 5 suggest
that it is possible to establish a fair trade-off between the number of bits
and the number of zoom steps in order to achieve increasingly accurate
results, as already observed for the symmetric GEVP (see, e.g., Refs. [28,
29]). One should recall that the lower the number of bits, the lower the
dimension of the QUBO matrix, and thus the less difficult the eigenpair search
becomes.

Figure 5: (Color online)[SA] Relative precision
$|1-\bar{\lambda}_{best}(z)/\lambda_{true}|$, at fixed zoom step, vs. the
number of bits $b$. In the calculations, we used $N_{A}^{GP}=200$,
$N_{A}^{DP}=20$. The dots represent the mean of the relative precision
obtained from $N_{run}=500$ independent samples.

In Fig. 6, the relative precision at fixed zoom step was studied by changing
the number of annealing cycles during the gradient-descent phase,
$N_{A}^{DP}$, keeping fixed $N_{A}^{GP}=200$ and $b=3$. The figure shows that
there is basically no improvement by increasing $N^{DP}_{A}$. In the gradient-
descent phase for a given $z$, the algorithm requires a minimal number of
annealing cycles to find the best minima for ${\bm{\delta}}(z)$. This is one
of the advantage of using such an algorithm.

Figure 6: (Color online)[SA] The relative precision
$|1-\bar{\lambda}_{best}(z)/\lambda_{true}|$, at fixed zoom step, vs.
$N_{A}^{DP}$, the number of annealing cycles in the gradient-descent phase. In
the calculations, $N_{A}^{GP}=200$ and $b=3$ were used. The dots represent the
mean of the relative precision performed over $N_{run}=500$.

In conclusion, the campaign of numerical calculations performed on the SA has
suggested to adopt the following input parameters for the studies of our
algorithm on the QA: $N^{GP}_{A}=200$, $N^{DP}_{A}=20$, and $b=2$ or $b=3$,
but with a suitable number of zoom steps. This preliminary analysis is
critical to fine-tune the strategy when the QA comes into play and a limited
run-time is available. It also sheds light on the different behavior of QA and
SA, with obvious practical implications for the analysis of large matrices.

###  IV.3  Analysis of the QA results

We perform the experimental study of the algorithm by using the Advantage
system 4.1 quantum annealer, provided by D-Wave Systems. As a first step, we
adopted $N^{GP}_{A}=200$ and $N^{DP}_{A}=20$ in order to study the whole set
of matrices whose relevant eigenpairs we aim to determine. As for the number
of bits, we used $b=3$ for $n_{M}<16$ and $b=2$ for $n_{M}\geq 16$, in order
to avoid limitations on the matrix dimensions that can be processed by the QA.
In fact, the current QA has a given upper-bound on the number of logical
qubits to be mapped onto sets of physical qubits, organized with an assigned
topology222This process is called minor-embedding or simply embedding. Recall
that the matrix dimension in the QUBO problem is $n_{M}\times b$, and
therefore $n_{M}\times b$ logical qubits are fully connected with weights
given by the values of the matrix elements. For the Advantage 4.1 system [26,
27], based on the Pegasus topology of the physical qubits, a maximal size of a
fully-connected set of logical qubits (clique) is $177$, but it can be reached
under very peculiar conditions, not fully achievable in our actual problem, if
accurate results have to be obtained.

In Tab. 3 we show the summary of the QA results for the problem in Eq. (1),
with dimensions $n_{M}=4,8,12,16,24,32$. As target we selected a nominal
relative precision on the eigenvalue of $\epsilon^{\prime}=10^{-3}$ that in
the algorithm corresponded to a value of $z=9$.

In the fifth column, we present the mean value $\bar{\lambda}_{best}$ obtained
running the algorithm $N_{run}$ times. The agreement between the true
eigenvalue, $\lambda_{true}$ in the fourth column, and the mean value
$\bar{\lambda}_{best}$ is smaller than the nominal precision requested
$\epsilon^{\prime}=10^{-3}$. More importantly, all the results obtained in our
runs are better than the requested precision. This can be observed by
comparing the $68\%$ confidence interval over the $N_{run}$ runs with the
nominal precision $\epsilon^{\prime}=10^{-3}$.

In Tab. 3 we also show: i) the ratio
$\bar{\lambda}^{I}_{best}/\bar{\lambda}_{best}$, and ii) the mean Euclidean
distance $\overline{||{\bf v}_{true}-{\bf v}_{best}||}$. It is worth noting
that the ratio and the mean Euclidean distance have almost the same magnitude,
and slightly increases with the matrix dimension $n_{M}$. In view of this, a
remarkable degree of reliability of the proposed algorithm can be inferred.

$n_{M}$ | $b$ | $N_{run}$ | $\lambda_{true}$ | $~{}~{}\bar{\lambda}_{best}$ | $\bar{\lambda}^{I}_{best}/\bar{\lambda}_{best}$ | $\overline{||{\bf v}_{true}-{\bf v}_{best}||}$
---|---|---|---|---|---|---
4 | 3 | 80 | $0.188026$ | $0.188012^{+1\cdot 10^{-5}}_{-7\cdot 10^{-6}}$ | $0.00024^{+2\cdot 10^{-5}}_{-5\cdot 10^{-6}}$ | $0.00024^{+3\cdot 10^{-5}}_{-2\cdot 10^{-5}}$
8 | 3 | 80 | $0.188204$ | $0.18820^{+2\cdot 10^{-5}}_{-2\cdot 10^{-5}}$ | $0.0003^{+1\cdot 10^{-4}}_{-1\cdot 10^{-4}}$ | $0.0003^{+1\cdot 10^{-4}}_{-1\cdot 10^{-4}}$
12 | 3 | 80 | $0.188203$ | $0.18821^{+2\cdot 10^{-5}}_{-2\cdot 10^{-5}}$ | $0.0005^{+1\cdot 10^{-4}}_{-1\cdot 10^{-4}}$ | $0.0006^{+2\cdot 10^{-4}}_{-1\cdot 10^{-4}}$
16 | 2 | 80 | $0.188203$ | $0.18820^{+4\cdot 10^{-5}}_{-3\cdot 10^{-5}}$ | $0.0009^{+1\cdot 10^{-4}}_{-1\cdot 10^{-4}}$ | $0.0011^{+2\cdot 10^{-4}}_{-2\cdot 10^{-4}}$
24 | 2 | 80 | $0.188225$ | $0.18822^{+5\cdot 10^{-5}}_{-4\cdot 10^{-5}}$ | $0.0013^{+1\cdot 10^{-4}}_{-3\cdot 10^{-4}}$ | $0.0015^{+3\cdot 10^{-4}}_{-2\cdot 10^{-4}}$
32 | 2 | 200 | $0.188225$ | $0.18823^{+4\cdot 10^{-5}}_{-3\cdot 10^{-5}}$ | $0.0016^{+2\cdot 10^{-4}}_{-2\cdot 10^{-4}}$ | $0.0018^{+4\cdot 10^{-}4}_{-3\cdot 10^{-}4}$
Table 3: [QA] Results obtained by using the Advantage 4.1 quantum annealer,
provided by D-Wave Systems. In the Table, $n_{M}$ is the dimension of the
matrices $A$, $B$ and $C$ involved in the calculations (see Eqs. (1) and (11),
respectively); $b$ is the number of bits adopted for the matrix representation
in the QUBO problem (see Eq. (23)); $N_{run}$ is the number of runs of the
entire algorithm (guess phase + gradient-descent phase);
$\bar{\lambda}_{best}$ is the average of
$\lambda_{best}(z_{fin})=\lambda^{R}({\bf w}_{\alpha_{fin}}^{z_{fin}})$ (cf.
Eq. (14) and the definition below Eq. (31)), over the outcomes returned after
running the entire algorithm $N_{run}$ times; in
$\bar{\lambda}^{I}_{best}/\bar{\lambda}_{best}$ the numerator
$\bar{\lambda}^{I}_{best}$ is obtained as $\bar{\lambda}_{best}$, but
following Eq. (15); $\overline{||{\bf v}_{true}-{\bf v}_{best}||}$ is the mean
Euclidean distance between the true eigenvector and the one from the QA
(recall that the vectors ${\bf v}_{i}$ are normalized). In addition to the
mean values, the $68\%$ confidence interval is shown. All the results are
obtained using $z=9$ zoom steps, corresponding to an upper-bound for the
relative precision equal to $\epsilon^{\prime}=10^{-3}$.

In Fig. 7 we plot the two dimensional distribution of the results obtained
running the code for $N_{run}=200$ in the $n_{M}=32$ case for $z=6$ and $z=9$
as function of $(1-\lambda_{best}/\lambda_{true})$, and the Euclidean distance
$||{\bf v}_{true}-{\bf v}_{best}||$, where $\lambda_{best}$ and ${\bf
v}_{best}$ are the solution found by the algorithm (notice the difference with
$\overline{\lambda}_{best}$). It is worth noticing that the relevant area
substantially shrinks when passing from the distribution of the results for
$z=6$ to $z=9$, as shown by the changes in the $x-y$ scales, but it is not
possible to identify a direct correlation between the two quantities. However,
the results seems to accumulate on a circular region that does not change
shape increasing $z$ indicating that the correlation does not depend on the
specific zoom phase. We observed also that the projection on the $x$ axis, the
distribution of $(1-\lambda_{best}/\lambda_{true})$, approaches a Gaussian
distribution centered in zero with a smaller and smaller width as $z$
increases, while the projection on the $y$ axis, the Euclidean distance,
starts from zero and extend to positive values, still with a Gaussian fall-
off. Plainly, these distributions are very important for a better
understanding of the quantum-hardware behavior and should be carefully
investigated in a future work, when a much longer QA running-time will be
available for exploring larger $N_{run}$ and lower $\epsilon^{\prime}$.

Figure 7:  (Color online) [QA] Two dimensional distribution of the results for
the $n_{M}=32$ case, after running $N_{run}=200$ times the entire algorithm as
a function of the Euclidean distance $||{\bf v}_{true}-{\bf v}_{best}||$ vs.
the precision $(1-\lambda_{best}/\lambda_{true})$. Upper panel: $z=6$ results.
Lower panel: $z=9$ results.

Finally, an interesting comparison between the outcomes obtained for matrix
dimensions $n_{M}=4$ and $32$ from SA and QA, is shown in Figs. 8, 9, and 10,
11, for both the eigenvalue and the eigenvector, respectively. For matrix
dimensions $n_{M}=4$, the analysis has been carried out by using $b=3$ and
pushing the number of zoom step beyond $z=30$ (recall that in the calculations
shown in previous figures the maximal value is $z=9$, leading to a nominal
precision $\epsilon^{\prime}=10^{-3}$). For the case $n_{M}=32$, the number of
bits was reduced to $b=2$, as already explained, and the zoom steps exceeded
$z=20$.

In Figs. 8 and 9, the mean relative precision
$|1-\bar{\lambda}_{best}(z)/\lambda_{true}|$ for fixed $z$ is shown for
$n_{M}=4$ and $n_{M}=32$, respectively. Red squares and arrows represent the
mean and the worst result obtained over $N_{run}=1000$ simulations,
respectively. Green triangles and arrows are the same as before, but running
the entire algorithm on QA by $N_{run}=10$ times for $n_{M}=4$ and by
$N_{run}=8$ times for $n_{M}=32$. We reduced the number of runs used for this
study compared to the one in Table 3, because of the larger running-time
requested to go beyond $z=9$. Finally, blue dots are the results obtained
using the exact, classical solver for the QUBO problem at each zoom step.
These blue dots are not available for $n_{M}=32$, since it should be necessary
to poll $2^{64}$ states and infer the suitable distributions (i.e. the well-
known Feynman’s argument in favor of a quantum computer).

As clearly seen in Figs. 8 and 9, the actual accuracy improves substantially
as $z$ increases for both SA and QA. In general, it remains below the nominal
precision $1/2^{z}$, given by the dashed lines. For $n_{M}=4$, the results
obtained with SA and QA are below $1/2^{z}$ up to $z=25$, reaching a plateau
with a rough value of $\sim 10^{-8}$. For $n_{M}=32$, the SA reach a plateau
roughly around $\sim 10^{-6}$ while the QA shows a behavior close to the one
seen for $n_{M}=4$. From the comparison between $n_{M}=4$ and $n_{M}=32$ SA
results, one could ascribe the different behavior to a lower SA’s performance
in managing large matrices, rather than an issue with the algorithm. In fact,
the results obtained with the QA just follows the expected behavior. One could
understand heuristically the presence of a plateau by considering that the
actual zero of the OF $f(C,{\bf w},\tilde{\lambda})$ is numerically about
$10^{-16}$, and therefore a precision for $\tilde{\lambda}$ of the order of
$\sqrt{f}$, i.e. $\sim 10^{-8}$, is expected. Properly rewriting Eq. (13), one
gets

$\displaystyle{f(C,{\bf w},\tilde{\lambda})\over|{\bf
w}|^{2}}=\Bigl{[}\tilde{\lambda}-\lambda^{R}({\bf
w})\Bigr{]}^{2}+\Bigl{[}\lambda^{I}({\bf w})\Bigr{]}^{2}~{},$ (33)

that for real $\tilde{\lambda}$ leads to the above estimate.

The Euclidean distance between the target eigenvector and the results
obtained, at each zoom step, with the SA and the QA, is shown in Figs. 10 and
11, for $n_{M}=4$ and $n_{M}=32$, respectively. The symbols are the same as in
Figs. 8 and 9, while the error bars represent the $99.9\%$ confidence level.
As expected, the Euclidean distance decreases by increasing $z$, until it
reaches a plateau following a pattern similar to the eigenvalue relative
precision. Again, for $n_{M}=32$, the QA performs better than the SA.

Summarizing, the proposed algorithm opens an actual window on the study of the
spectrum of non symmetric matrices by using a QA, that nicely seems
outperforms the software providing the simulated annealing. Clearly, the main
issue is given by the size of the matrices to be used, but this will be surely
overcome by the planned improvements of the QPU topology, so that larger and
larger matrices will be embedded.

Figure 8:  (Color online) [SA - QA] Mean relative precision,
$|1-\bar{\lambda}_{best}(z)/\lambda_{true}|$, vs. the zoom step in the
gradient-descent phase, for $n_{M}=4$ case with $b=3$. The dashed line is the
nominal precision at each zoom step, equal to $1/2^{z}$. Red squares and
arrows represent the mean and the worst result respectively, obtained by
averaging over $N_{run}=1000$ runs on the SA. Green triangles and arrows are
the same as before obtained with the QA, but averaging over $N_{run}=10$ runs.
Blue dots are the results obtained using the exact solver for the QUBO problem
at each zoom step.  Figure 9:  (Color online) [SA - QA] The same as Fig. 8,
but for $n_{M}=32$ and $b=2$. The QA results are obtained by averaging over
$N_{run}=8$ runs (see text).  Figure 10:  (Color online) [SA - QA] Mean
Euclidean distance between the target eigenvector and the true one, at each
zoom step, for the $n_{M}=4$ case with $b=3$. Green triangles: results from
the QA. Red squares: outcomes from the SA. Error bars represent the $99\%$
confidence level over 10 (1000) runs on the QA (SA). Figure 11:  (Color
online) [SA - QA] The same as Fig. 10, but for $n_{M}=32$ and $b=2$. The QA
results are averaged over $N_{run}=8$ runs.

###  IV.4 Scalability on the QA

A first study of the scalability of the proposed hybrid algorithm running on
the QA has been carried out by focusing on the total time spent in the
annealing cycles. In our case, the total annealing-time is theoretically given
by $3\times N^{GP}_{A}\times t+z^{max}\times i^{max}\times(N^{DP}_{A}\times
t)$ (cf. Tabs. 1 and 2) where $t\sim 20\mu$s is the annealing time selected,
$z^{max}$ is the maximal value of zoom steps adopted in the gradient-descent
phase, e.g. for the results shown in the third column of Tab. 3 one has
$z^{max}=9$. Moreover, $i^{max}$ is the maximal number of search cycles for a
given $z$ (see step 4 in Tab. 2). In Tab. 4, $T[ms]$ is the total annealing-
time averaged on $N_{run}$ runs of the two-phase algorithm. One notices that a
slightly increasing behavior is present when the matrix dimension,
$n_{M}\times b$, increases, at fixed $b$. But, further studies with larger
matrices will be necessary for fully assessing a linear increasing of the
total annealing time on the QA.

Another parameter to be considered when studying the scalability of the
algorithm is the number of physical qubits, $N_{qubits}$, used by the annealer
to map the logical qubits of the QUBO problem on the topology of the hardware
(Pegasus in this case), through the embedding process. In general, since each
physical qubit is not connected to all the other, the mapping is not a one-to-
one correspondence, but it is a one-to-many, i.e. one logical qubit is
represented by a set of physical qubits, all in principle forced to have the
same value. This formal step is elaborated by the proprietary software D-Wave
Ocean, before a problem is submitted to the quantum annealer. The software is
heuristic and a different embedding, with a different number of physical
qubits, is used at each QA run. However, the exact number of physical qubits
adopted for mapping the QUBO square matrix is provided by the software, so
that the distribution of $N_{qubits}$ over different runs can be eventually
obtained. The mean value and the $68\%$ confidence level of the number of this
distribution is then showed in the last column of Table 4. It should be noted
that the number of involved physical qubits is roughly given by the matrix
dimension times the so-called chain length, i.e. the number of qubits needed
for representing a single logical qubit. In our analysis, the chain length
linearly grows with the matrices dimension as shown in Ref. [26].333 In Ref.
[26], for Advantage 4.1, it is shown that one can reach at most
$N_{clique}=177$ fully-connected logical bits, with a chain-length equal to
17, while for $N_{clique}=647$ the chain length is about 7, not too far from
what we found. It should be recalled the critical role of the accuracy one
aims to, as well as the condition number of the involved matrices. In our
case, the matrix elements, being the outcomes of an actual physical problem,
are within an interval several orders of magnitude wide. The overall memory
required by the algorithm grows therefore quadratically with the dimension of
the original matrix.

$n_{M}\times b$ | $N_{run}$ | $T[ms]$ | $N_{qubits}$
---|---|---|---
12(4x3) | 80 | $18.5^{+0.9}_{-0.7}$ | $24.1^{+0.9}_{-1.1}$
24(8x3) | 80 | $21.7^{+1.3}_{-0.9}$ | $81.3^{+2.1}_{-2.3}$
32(16x2) | 80 | $23.0^{+1.4}_{-1.0}$ | $140.9^{+1.9}_{-1.9}$
36(12x3) | 80 | $23.0^{+1.4}_{-1.3}$ | $177.6^{+6.4}_{-6.6}$
48(24x2) | 80 | $23.0^{+1.0}_{-1.0}$ | $306.6^{+11.8}_{-15.6}$
64(32x2) | 200 | $23.6^{+1.2}_{-1.2}$ | $529.8^{+30.6}_{-33.8}$
Table 4: [QA] Total annealing time, $T[ms]$, and total number of physical
qubits $N_{qbits}$, averaged on $N_{run}$ runs, for the results shown in Tab.
3.

##  V Conclusions and Perspectives

A hybrid algorithm, suitable for a quantum annealer, was implemented to
evaluate the largest real eigenvalue and corresponding eigenvector of a
generalized eigenvalue problem involving a non symmetric matrix. This
numerical problem stems from the discretization of the homogeneous Bethe-
Salpeter equation describing a bound state of two massive scalars, that
interact by exchanging a massive scalar (see, e.g., Ref. [37]). In the current
initial stage of the study, the non singularity of the symmetric matrix was
exploited, so that a classical $LDL^{T}$ factorization of the symmetric matrix
was used in order to deal with a simpler QUBO problem.

The numerical results were obtained by running our two-phase alogorithm both
on an Advantage 4.1 quantum annealer, provided by the D-Wave Systems, and a
simulator, based on the proprietary software Ocean. We first tested our
algorithm, based on an OF suggested in Ref. [38] and supplemented by the
result of the valuable Gershgorin circle theorem [58], on the SA. This first
investigation led to establish a practical set of input parameters, given by:
i) the number of bits for expressing the real components of the involved
vectors in the binary basis, and ii) the number of annealing cycles in both
the guess phase and the gradient-descent one. Notably, the studies carried out
on the SA have confirmed that the trade-off between the bit number and the
zoom-step number, that controls the nominal precision, is favorable as already
found in the symmetric case [28, 29]. Then, a minimal bit number can be
chosen, so that the matrix dimension in the QUBO problem does not exceed the
current limitations of the quantum hardware. After strengthening our numerical
experience on the SA, we performed a numerical campaign by running $N_{run}$
times our algorithm on the D-Wave QA, obtaining very encouraging results. We
successfully approached the target eigenpair, extending our studies up to a
matrix dimension $n_{M}=32$, with a corresponding QUBO-matrix dimension equal
to $64$. We have shown that the algorithm used in combination with the QA is
able to compute the eigenpair corresponding to the largest eigenvalue with
$100\%$ reliabilty and with improvable precision up to $\sim 10^{-8}$. As to
the scalability of the algorithm, Tab. 4 yields a promising slightly linear
increasing for growing matrix dimension, but more studies have to be performed
with larger and larger $n_{M}$, before drawing definite conclusions.

With an eye to the future, the next challenge is to improve the algorithm in
order to address the generalized eigenvalue problem in its full glory, i.e.
without exploiting the non singularity of the symmetric matrix.

###### Acknowledgements.

The Authors gratefully thank Q@TN - Quantum Science and Technology in Trento,
sponsored by Università degli Studi di Trento, FBK, INFN, CNR and CINECA, for
providing access and run-time to the CINECA Quantum Computing facility (during
the period 2021-2023), where all the numerical results shown in the present
work were obtained by using a D-Wave Systems quantum annealer. A.G.
acknowledges support from the DOE Topical Collaboration “Nuclear Theory for
New Physics,” award No. DE-SC0023663 and Jefferson Lab that is supported by
the U.S. Department of Energy, Office of Nuclear Science, under Contracts No.
DE-AC05-06OR23177. During part of the development of this work A.G. was a
postdoc at ECT∗. T. F. thanks the financial support from CNPq (Conselho
Nacional de Desenvolvimento Científico e Tecnológico) grant 306834/2022-7,
CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), Finance
Code 001, FAPESP (grant 2019/07767-1) and Instituto Nacional de Ciência e
Tecnologia - Física Nuclear e Aplicações Proc. No. 464898/2014-5.

##  Appendix A The homogeneous Bethe-Salpeter equation

In this Appendix, some details about the physical case, which suggested the
non symmetric GEVP to be investigated with the QA, are briefly illustrated.

The simplest hBSE is the one that allows to dynamically describe a bound
system composed by two massive scalars interacting through the exchange of a
massive scalar (the case with a simple ladder-exchange of a massless scalar
between massless scalars is known as the Wick-Cutkowsky model [59, 60]). In 4D
Minkowski space, it reads (see Ref. [37] for more details)

$\displaystyle\Phi_{b}(k,p)=G_{0}^{(1)}(k,p)~{}G_{0}^{(2)}(k,p)~{}\int\frac{d^{4}k^{\prime}}{(2\pi)^{4}}i~{}{\cal
K}(k,k^{\prime},p)\Phi_{b}(k^{\prime},p)~{},$ (34)

where $\Phi_{b}(k,p)$ is the BS amplitude of a two-body bound system,
$p=p_{1}+p_{2}$ the total momentum of the system, with square mass
$M^{2}=p^{2}$, $k=(p_{1}-p_{2})/2$ the relative momentum, and $i~{}{\cal K}$
the interaction kernel, that contains all the irreducible diagrams [33]. In
Eq.(34) , and $G_{0}^{(i)}(k,p)$ is the free scalar propagator, given by

$\displaystyle G_{0}^{(i)}(k,p)=\frac{i}{(\frac{p}{2}\pm
k)^{2}-m^{2}+i\epsilon}~{}~{},$ (35)

with $m$ the mass of the scalar constituents.

To proceed for obtaining actual numerical solutions is helpful to adopt the
NIR [46, 36] of the BS amplitude, so that its analytic structure is made
explicit. Within the NIR framework, the BS amplitude is written as a proper
folding of i) a non singular weight function that depends upon real variables
and appears in the numerator (for the system under consideration, one variable
is compact and the other is non compact), and ii) a denominator that contains
the analytic structure. One writes

$\displaystyle\Phi_{b}(k,p)=~{}i~{}\int_{-1}^{1}dz^{\prime}\int_{0}^{\infty}d\gamma^{\prime}~{}{g_{b}(\gamma^{\prime},z^{\prime};\kappa^{2})\over\left[\gamma^{\prime}+\kappa^{2}-{k}^{2}-p\cdot
kz^{\prime}-i\epsilon\right]^{2+n}}$ (36)

where $g_{b}(\gamma^{\prime},z^{\prime};\kappa^{2})$ is the Nakanishi weight
function (NWF), $n\geq 1$ is a suitable power that can be chosen with some
degree of arbitrariness [36] (one should properly redefine the NWF) and
$\kappa^{2}$ is given by

$\displaystyle\kappa^{2}=m^{2}-{M^{2}\over 4}~{},$ (37)

that is a measure of the binding, defined by $B=2m-M\geq 0$. Inserting such a
relation in Eq. (37), one gets $\kappa^{2}>0$ for bound states. Interestingly,
the dependence upon $z^{\prime}$ of
$g_{b}(\gamma^{\prime},z^{\prime};\kappa^{2})$ is even as expected by the
symmetry property of the BS amplitude for the two-scalar system and
$g_{b}(\gamma^{\prime},z^{\prime}=\pm~{}1;\kappa^{2})=0$, as illustrated in
Ref. [37]. It is crucial to emphasize that once the NWF is known, then the BS
amplitude can be reconstructed via Eq. (36) and the evaluation physical
observables can be carried out.

If the 4D kernel , $i{\cal K}$ is explicitly known, one can perform the
analytic integration of both sides of Eq. (34). To this end, instead of using
the Minkowskian four-momentum $k$ in Cartesian coordinates it is very useful
to introduce their light-front (LF) combination, that amount to
$k\equiv\\{k^{\pm}=k^{0}\pm k^{3},{\bf k}_{\perp}\\}$ (for interested readers,
Ref. [61] yields details on the Dirac proposal of the light-front dynamics,
and in general on the relativistic Hamiltonian description of a dynamical
system). Hence, one gets an integral equation for the NWF, avoiding the
difficulties related to the indefinite metric of the Minkowski space. In
particular, the integral equation reads

$\displaystyle\int_{0}^{\infty}d\gamma^{\prime}~{}\frac{g_{b}(\gamma^{\prime},z;\kappa^{2})}{[\gamma^{\prime}+\gamma+z^{2}m^{2}+(1-z^{2})\kappa^{2}-i\epsilon]^{2}}=\int_{0}^{\infty}d\gamma^{\prime}\int_{-1}^{1}dz^{\prime}\;V^{LF}_{b}(\gamma,z;\gamma^{\prime},z^{\prime})g_{b}(\gamma^{\prime},z^{\prime};\kappa^{2})\
.$ (38)

where the new kernel $V^{LF}_{b}$, so-called Nakanishi kernel, is related to
the 4D BS kernel, $i{\cal K}$, in Eq. (34), as follows

$\displaystyle
V^{LF}_{b}(\gamma,z;\gamma^{\prime},z^{\prime})=~{}ip^{+}~{}\int_{-\infty}^{\infty}{dk^{-}\over
2\pi}~{}G_{0}^{(1)}(k,p)~{}G_{0}^{(2)}(k,p)$
$\displaystyle\times\int\frac{d^{4}k^{\prime}}{(2\pi)^{4}}\frac{i{\cal
K}(k,k^{\prime},p)}{\left[{k^{\prime}}^{2}+p\cdot
k^{\prime}z^{\prime}-\gamma^{\prime}-\kappa^{2}+i\epsilon\right]^{3}}~{},$
(39)

with $p^{+}=M$, once the intrinsic frame where ${\bf p}_{\perp}=0$ is chosen
(one adopts the rest frame for carrying out the calculation, given our
interest on the intrinsic properties of the system).

Two comments are in order. Given the dependence on $k^{2}$ in both scalar
propagators and interaction kernel, one expects double poles in the Cartesian
variable $k^{0}$. If one uses LF variables, one remains with single poles in
$k^{-}$ and $k^{+}$, obtaining a simplification in the formal treatment. To
numerically evaluate the integral equation (38), one has to specify the
structure of the interaction kernel. A quantitative analysis was carried out
by using the exchange of one massive scalar particle, that generates the
ladder approximation of the hBSE, by iteration. This means that the integral
equation automatically produces an infinite number of scalar exchanges, so
that the bound-system pole of the four-leg Green’s function can be
established. The ladder kernel is given by

$\displaystyle i{\cal
K}^{(Ld)}(k,k^{\prime})={i(-ig)^{2}\over(k-k^{\prime})^{2}-\mu^{2}+i\epsilon}~{},$
(40)

where $g$ is the coupling constant and $\mu$ the mass of the exchanged scalar.
A quantitative analysis of Eq. (38) with the kernel in Eq. (40) was carried
out in Ref. [37] by using a complete orthonormal basis for expanding the NWF
$g_{b}(\gamma,z;\kappa^{2})$. In particular, one writes

$\displaystyle
g^{(Ld)}_{b}(\gamma,z;\kappa^{2})=\sum_{\ell=0}^{N_{z}}\sum_{j=0}^{N_{\gamma}}~{}{\cal
C}_{\ell j}(\kappa^{2})~{}G_{\ell}(z)~{}{\cal L}_{j}(\gamma)\ ,$ (41)

where i) the functions $G_{\ell}(z)$ are given in terms of even Gegenbauer
polynomials, $C^{(5/2)}_{2\ell}(z)$, by

$\displaystyle
G_{\ell}(z)=4~{}(1-z^{2})~{}\Gamma(5/2)~{}\sqrt{{(2\ell+5/2)~{}(2\ell)!\over\pi\Gamma(2\ell+5)}}~{}C^{(5/2)}_{2\ell}(z)\
,$ (42)

and ii) the functions ${\cal L}_{j}(\gamma)$ are expressed in terms of the
Laguerre polynomials, $L_{j}(a\gamma)$, by

$\displaystyle{\cal
L}_{j}(\gamma)=\sqrt{a}~{}L_{j}(a\gamma)~{}e^{-a\gamma/2}~{}~{}.$ (43)

The following orthonormality conditions are fulfilled

$\displaystyle\int^{1}_{-1}dz~{}G_{\ell}(z)~{}G_{n}(z)=\delta_{\ell
n}~{}~{},\quad\int_{0}^{\infty}~{}d\gamma~{}{\cal L}_{j}(\gamma)~{}{\cal
L}_{\ell}(\gamma)=a\int_{0}^{\infty}~{}d\gamma~{}e^{-a\gamma}~{}L_{j}(a\gamma)~{}L_{\ell}(a\gamma)=~{}\delta_{j\ell}~{}.$
(44)

By inserting the expansion of the NWF and approximating the kernel with only
the one-scalar exchange (i.e. the ladder approximation), Eq. (38) becomes

$\displaystyle\sum_{\ell^{\prime}j^{\prime}}\int_{0}^{\infty}d\gamma\int_{-1}^{1}dz\int_{0}^{\infty}d\gamma^{\prime}\int_{-1}^{1}dz^{\prime}~{}G_{\ell}(z){\cal
L}_{j}(\gamma)~{}{\cal
B}(\gamma,z;\gamma^{\prime},z^{\prime};\kappa^{2})~{}G_{\ell^{\prime}}(z^{\prime}){\cal
L}_{j^{\prime}}(\gamma^{\prime})~{}{\cal
C}_{\ell^{\prime}j^{\prime}}(\kappa^{2})$
$\displaystyle=\alpha~{}\sum_{\ell^{\prime}j^{\prime}}\int_{0}^{\infty}d\gamma\int_{-1}^{1}dz\int_{0}^{\infty}d\gamma^{\prime}\int_{-1}^{1}dz^{\prime}~{}G_{\ell}(z){\cal
L}_{j}(\gamma)~{}{\cal
V}^{Ld}_{b}(\gamma,z;\gamma^{\prime},z^{\prime})~{}G_{\ell^{\prime}}(z^{\prime}){\cal
L}_{j^{\prime}}(\gamma^{\prime})~{}{\cal
C}_{\ell^{\prime}j^{\prime}}(\kappa^{2}).$ (45)

where $\alpha=g^{2}/(32\pi^{2})$, ${\cal
B}(\gamma,z;\gamma^{\prime},z^{\prime};\kappa^{2})$ is a function symmetric
under the exchange $\gamma\to\gamma^{\prime}$ and $z\to z^{\prime}$ given by

$\displaystyle{\cal
B}(\gamma,z;\gamma^{\prime},z^{\prime};\kappa^{2})={\delta(z-z^{\prime})\over[\gamma^{\prime}+\gamma+z^{2}m^{2}+(1-z^{2})\kappa^{2}-i\epsilon]^{2}}~{},$
(46)

and ${\cal V}^{Ld}_{b}(\gamma,z;\gamma^{\prime},z^{\prime})$ is a function non
symmetric in the above exchange and is obtained from the following integral
(see Ref. [37])

$\displaystyle{\cal
V}^{Ld}_{b}(\gamma,z;\gamma^{\prime},z^{\prime})=~{}32\pi~{}p^{+}~{}\int_{-\infty}^{\infty}{dk^{-}\over
2\pi}~{}~{}\int\frac{d^{4}k^{\prime}}{(2\pi)^{4}}\frac{G_{0}^{(1)}(k,p)~{}G_{0}^{(2)}(k,p)}{\left[(k-k^{\prime})^{2}-\mu^{2}+i\epsilon\right]~{}\left[{k^{\prime}}^{2}+p\cdot
k^{\prime}z^{\prime}-\gamma^{\prime}-\kappa^{2}+i\epsilon\right]^{3}}~{}.$
(47)

By applying the above steps, one is able to formally transform the 4D hBSE
into a GEVP with the following structure

$\displaystyle{A}~{}{\bf v}_{i}=\lambda_{i}~{}{B}~{}{\bf v}_{i}~{},\quad{\rm
with~{}}i=1,2\dots~{},n~{},$ (48)

where ${A}$ and ${B}$ are real square matrices, calculated from rhs and lhs of
Eq. (45) respectively. Moreover, the matrix $A$ is non symmetric, while $B$ is
symmetric, and they have dimension $n=N_{z}\times N_{\gamma}$. The eigenvector
${\bf v}_{i}$ , corresponding to the i-th eigenvalue
$\lambda_{i}=1/\alpha_{i}$, contains the coefficients ${\cal C}_{\ell
j}(\kappa^{2})$ of the expansion of the NWF (cf. Eq. (41)). It should be
pointed out that the eigenvalues and eigenvectors can be real and complex
conjugated, but we are interested only in real eigenvalues, since we aim to to
describe a physical system (recall that the eigenvalue is the inverse of the
square coupling constant, in the ladder approximation we adopted).

##  Appendix B The binary basis

This Appendix illustrates how to transform a QUBO problem in a form suitable
for the evaluation with a quantum annealer. The needed step is given by the
transformation of each component of the involved real vectors into their
expression in terms of a binary string, by using $b$ binary variables. For
instance, a generic real number, $a$, such that $|a|\leq 1$, can be expressed
in terms of a binary basis with $b$ bits as follows (see also Ref. [28, 29]
for details)

$\displaystyle a=-q_{b}+q_{1}~{}{1\over 2}+q_{2}~{}{1\over
2^{2}}+q_{3}~{}{1\over 2^{3}}+~{}\dots~{}+q_{b-1}~{}{1\over 2^{b-1}}~{},$

where $a\in[-1,1)$ (see below), $q_{i}=0,1$ and the $b$-th bit carries the
sign, i.e. for $a\geq 0$ ( $a<0$) one uses $q_{b}=0$ ($q_{b}=1$). In a compact
form, one can rewrite $a$ as a scalar product between two vectors with $b$
components, viz.

$\displaystyle a={\bf p}^{T}\cdot{\bf q}_{b}~{},$ (50)

where ${\bf q}_{b}\equiv\\{q_{1},q_{2},\dots,q_{b}\\}$ is a binary string (a
column vector in our notation) and ${\bf
p}^{T}\equiv\\{-1,1/2,1/4,~{}\dots~{},1/2^{b-1}\\}$ is a row vector,
representing the transpose of the so-called precision vector, since
$1/2^{b-1}$ controls the precision to which one wants to approximate the
number $a$. Notice that $1$ can be reached only asymptotically, since

$\displaystyle\lim_{b\to\infty}\sum_{i=1}^{b-1}{1\over 2^{i}}={1\over
2}~{}\lim_{b\to\infty}\left[{1-2^{-(b-1)}\over 1/2}\right]\to 1~{},$ (51)

while $-1$ corresponds to the binary vector $1,0,\dots~{}\dots,0$, truncated
at any order. Therefore, one has to consider the right-open interval
$a\in[-1,1)$.

One quickly generalizes the compact expression in Eq. (50) to a vector ${\bf
v}$ with real components and dimension $n$, since one has

$\displaystyle\begin{pmatrix}v_{1}\\\ v_{2}\\\ \vdots\\\
v_{n}\end{pmatrix}=\begin{pmatrix}~{}{\bf p}^{T}&0&\dots&0\\\ 0&~{}{\bf
p}^{T}&\dots&0\\\ \vdots&\vdots&\ddots&0\\\ 0&0&\dots&{\bf
p}^{T}\end{pmatrix}~{}\begin{pmatrix}{\bf x}_{1b}\\\ {\bf x}_{2b}\\\ \vdots\\\
{\bf x}_{nb}\end{pmatrix}$ $\displaystyle=P^{T}~{}{\bf X}~{},$ (52)

where ${\bf x}_{ib}$ is the $i$-th binary vector with dimension $b$. The
vector ${\bf X}\equiv\\{{\bf x}_{i}\\}$ in the last line has dimension:
$n\times b$, and the rectangular matrix ${\bf P}^{T}$ has $n$ rows and
$n\times b$ columns. In the last line ${\bf X}\equiv\\{{\bf x}_{i}\\}$.

Hence, for a $n\times n$ matrix $C$, and a given vector ${\bf v}$ with real
components, one can write

$\displaystyle{\bf v}^{T}~{}C{\bf
v}=\begin{pmatrix}v_{1}&v_{2}&\dots&v_{n}\end{pmatrix}~{}C~{}\begin{pmatrix}v_{1}\\\
v_{2}\\\ \vdots\\\ v_{n}\end{pmatrix}=\begin{pmatrix}{\bf x}^{T}_{1b}&{\bf
x}^{T}_{2b}&\dots&{\bf x}^{T}_{nb}\end{pmatrix}\begin{pmatrix}~{}{\bf
p}&0&\dots&0\\\ 0&~{}{\bf p}&\dots&0\\\ \vdots&\vdots&\ddots&0\\\
0&0&\dots&{\bf p}\end{pmatrix}~{}C~{}\begin{pmatrix}~{}{\bf
p}^{T}&0&\dots&0\\\ 0&~{}{\bf p}^{T}&\dots&0\\\ \vdots&\vdots&\ddots&0\\\
0&0&\dots&{\bf p}^{T}\end{pmatrix}~{}\begin{pmatrix}{\bf x}_{1b}\\\ {\bf
x}_{2b}\\\ \vdots\\\ {\bf x}_{nb}\end{pmatrix}$ $\displaystyle={\bf
X}^{T}~{}PCP^{T}~{}{\bf X}~{},$ (53)

recalling that ${\bf p}$ is a column vector.

##  Appendix C Inequalities for non-symmetric-matrix eigenvalues

In this Appendix, some inequalities between the eigenvalues of a $n\times n$
non symmetric matrix and the ones of its symmetric part are briefly reviewed,
following the results of Ref. [62]. Indeed, in Ref. [62], the most general
case of a linear transformation in a $n$-dimensional unitary space was
considered, while we limit ourselves to our specific needs.

Let us introduce the symmetric combination of a real, non symmetric matrix $A$
and its transpose $A^{T}$

$\displaystyle S={A+A^{T}\over 2}~{}.$ (54)

Recall that $Tr\\{S\\}=\sum_{i}\\{\lambda_{i}(S)\\}=Tr\\{A\\}=\sum_{i}\Re
e\\{\lambda_{i}(A)\\}$ for real $A$, where $\lambda_{i}(S)$ (real) and
$\lambda_{i}(A)$ (both real and, possibly, complex conjugated) are the
eigenvalues of $S$ and $A$, with $i=1,~{}\dots~{},n$, respectively.

From Theorem II of Ref. [62], one can deduce that the maximal eigenvalue of
$S$ is an upper-bound of the maximal eigenvalue of $A$, as sketched in what
follows. It turns out that

$\displaystyle\sum_{i=1}^{q}\lambda_{i}(S)\geq\sum_{i=1}^{q}\Re
e\\{\lambda_{i}(A)\\}~{},~{}~{}~{}~{}{\rm with}~{}~{}q\leq n~{}~{},$ (55)

where the two sets of eigenvalues are ordered in a decreasing way, i.e.
$\lambda_{1}(S)\geq\lambda_{2}(S)\geq\dots~{}\dots\geq\lambda_{n}(S)$ and the
same for the real parts of $\lambda_{i}(A)$. If $q=1$, then one has the
following inequality

$\displaystyle\lambda_{max}(S)\geq\Re e\\{\lambda_{max}(A)\\}$ (56)

The minimal eigenvalue of $S$ is lower-bound of the minimal eigenvalue of $A$.
Recalling that $Tr\\{A\\}=Tr\\{S\\}$ one gets

$\displaystyle\sum_{i=1}^{n-1}\Re e\\{\lambda_{i}(A)\\}+\Re
e\\{\lambda_{min}(A)\\}=\sum_{i=1}^{n-1}\lambda_{i}(S)+\lambda_{\min}(S)~{},$

and using Eq. (55) with $q=n-1$ one eventually has

$\displaystyle\Re
e\\{\lambda_{min}(A)\\}=\Biggl{[}\sum_{i=1}^{n-1}\lambda_{i}(S)-\sum_{i=1}^{n-1}\Re
e\\{\lambda_{i}(A)\\}\Biggr{]}$
$\displaystyle+\lambda_{\min}(S)\geq\lambda_{\min}(S)~{}.$ (58)

It should be pointed out that the antisymmetric combination ${\cal
A}=(A-A^{T})/2$ leads to bounds on the imaginary part of the eigenvalues
$\lambda(A)$.

For the sake of completeness, one can mention that if i) the largest
eigenvalue of the symmetric matrix $S$ is positive and ii) $\lambda_{max}(A)$
is real and positive, then $\lambda_{max}=\lambda_{max}(S)$, after applying
Eq. (56). In fact, one has

$\|S\|_{p}\leq\|A\|_{p}=\|A^{T}\|_{p}\;,$ (59)

for any Schatten $p-$norm, that becomes a constraint between the spectral
norms for $p=\infty$, i.e.

$\begin{split}\|S\|_{\infty}&=\max\left[|\lambda_{min}(S)|,|\lambda_{max}(S)|\right]\\\
&\leq\max\left[|\lambda_{min}(A)|,|\lambda_{max}(A)|\right]=\|A\|_{\infty}\end{split}$
(60)

From the assumed assumptions, one can writes

$\displaystyle\|S\|_{\infty}=\lambda_{max}(S)\leq|\lambda_{max}(A)|=\|A\|_{\infty}$
(61)

and after taking into account Eqs. (56) one writes

$\begin{split}\Re
e&\\{\lambda_{max}(A)\\}^{2}\leq\lambda_{max}(S)^{2}\leq|\lambda_{max}(A)|^{2}\\\
&\quad\quad=\Re e\\{\lambda_{max}(A)\\}^{2}+\Im m\\{\lambda_{max}(A)\\}^{2}\\\
\end{split}$ (62)

This means that if the largest eigenvalue of $A$ is real then the largest
eigenvalue of $S$ will be equal to it. Unfortunately the assumptions are quite
restrictive and are not fulfilled by our actual matrices.

For the symmetric combination $AA^{T}$, one has an analogous result. For the
sake of simplicity (and since we are interested in), we focus on real
eigenvalues of the non symmetric matrix $A$. One writes

$\displaystyle A~{}{\bf x}_{\lambda(A)}=\lambda(A)~{}{\bf
x}_{\lambda}~{}~{},\quad{\bf x}^{T}_{\lambda(A)}A^{T}=\lambda(A)~{}{\bf
x}^{T}_{\lambda(A)}$ (63)

and

$\displaystyle{\bf x}^{T}_{\lambda(A)}~{}A^{T}A~{}{\bf
x}_{\lambda(A)}=\lambda^{2}(A)~{}{\bf x}^{T}_{\lambda(A)}{\bf
x}_{\lambda(A)}~{}.$ (64)

Let us now consider the eigenvalue problem for the symmetric matrix
$P=A^{T}A$. One has real eigenvalues and eigenvectors generating an
orthonormal basis. To establish formalism, one writes

$\displaystyle P~{}{\bf x}_{\lambda_{i}(P)}=\lambda_{i}(P){\bf
x}_{\lambda_{i}(P)}~{}.$ (65)

The expansion of an eigenvector of $A$ on the orthonormal basis ${\bf
x}_{\lambda_{i}(P)}$ yields

${\bf x}_{\lambda(A)}=\sum_{i}c_{i}~{}{\bf x}_{\lambda_{i}(P)}$

that leads to

$\displaystyle{\bf x}^{T}_{\lambda(A)}~{}A^{T}A~{}{\bf
x}_{\lambda(A)}=\Bigl{[}\lambda(A)\Bigr{]}^{2}$
$\displaystyle=\sum_{i}\lambda_{i}(P)c^{2}_{i}|{\bf x}^{T}_{\lambda(A)}{\bf
x}_{\lambda_{i}(P)}|^{2}$ $\displaystyle\geq\lambda_{min}(P)~{}{\bf
x}^{T}_{\lambda(A)}{\bf x}_{\lambda(A)}~{},$ (66)

and then

$\displaystyle\lambda_{max}(P)\geq\Bigl{[}\lambda(A)\Bigr{]}^{2}\geq\lambda_{min}(P)~{},$
(67)

where an obvious inequality involving $\lambda_{max}(P)$ has been used for the
lhs. In conclusion, again the minimal(maximal) eigenvalue of the naively
symmetrized matrix is a lower(upper) bound, and a variational method can
return $\lambda_{min}(P)\Bigl{(}\lambda_{max}(P)\Bigr{)}$ and not
$\lambda(A)$.

Another possible approach to be mentioned is based on the following symmetric
form

$\widetilde{A}=\begin{pmatrix}0&A\\\ A^{T}&0\end{pmatrix}~{}.$ (68)

Since the non symmetric matrices $A$ and its transpose has the same set of
eigenvalues but not the right, or left, eigenvectors one should solve the
following symmetric GEPV

$\displaystyle\begin{pmatrix}0&A\\\ A^{T}&0\end{pmatrix}~{}\begin{pmatrix}{\bf
w}_{i}\\\ {\bf v}_{i}\end{pmatrix}=\lambda_{i}~{}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}~{}~{}\begin{pmatrix}{\bf w}_{i}\\\ {\bf v}_{i}\end{pmatrix}$
(69)

where $A{\bf v}_{i}=\lambda_{i}~{}{\bf v}_{i}$ and $A^{T}{\bf
w}_{i}=\lambda_{i}~{}{\bf w}_{i}$. Clearly, this approach doubles the
dimension of the initial problem and therefore is confronted with the limits
of the matrix dimension tractable by using a quantum annealer.

## References

  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018), arXiv:1801.00862 [quant-ph] .
  * Preskill [2022] J. Preskill, The Physics of Quantum Information (2022) arXiv:2208.08064 [quant-ph] .
  * Bauer  _et al._ [2023] C. W. Bauer _et al._ , Quantum Simulation for High-Energy Physics, PRX Quantum 4, 027001 (2023), arXiv:2204.03381 [quant-ph] .
  * Funcke  _et al._ [2023] L. Funcke, T. Hartung, K. Jansen, and S. Kühn, Review on quantum computing for lattice field theory, in  _Proceedings of The 39th International Symposium on Lattice Field Theory — PoS(LATTICE2022)_ (Sissa Medialab, 2023).
  * Delgado  _et al._ [2022] A. Delgado, K. E. Hamilton, P. Date, J.-R. Vlimant, D. Magano, Y. Omar, P. Bargassa, A. Francis, A. Gianelle, L. Sestini, D. Lucchesi, D. Zuliani, D. Nicotra, J. de Vries, D. Dibenedetto, M. L. Martinez, E. Rodrigues, C. V. Sierra, S. Vallecorsa, J. Thaler, C. Bravo-Prieto, su Yeon Chang, J. Lazar, C. A. Argüelles, and J. J. M. de Lejarza, Quantum computing for data analysis in high energy physics (2022), arXiv:2203.08805 [physics.data-an] .
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Quantum chemistry in the age of quantum computing, Chemical reviews 119, 10856 (2019).
  * Bauer  _et al._ [2020] B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Quantum algorithms for quantum chemistry and quantum materials science, Chemical Reviews 120, 12685 (2020).
  * Emani  _et al._ [2021] P. S. Emani, J. Warrell, A. Anticevic, S. Bekiranov, M. Gandal, M. J. McConnell, G. Sapiro, A. Aspuru-Guzik, J. T. Baker, M. Bastiani, _et al._ , Quantum computing at the frontiers of biological sciences, Nature Methods 18, 701 (2021).
  * Flöther [2023] F. F. Flöther, The state of quantum computing applications in health and medicine, arXiv preprint arXiv:2301.09106 (2023).
  * Bova, Francesco  _et al._ [2021] Bova, Francesco, Goldfarb, Avi, and Melko, Roger G., Commercial applications of quantum computing, EPJ Quantum Technol. 8, 2 (2021).
  * Herman  _et al._ [2022] D. Herman, C. Googin, X. Liu, A. Galda, I. Safro, Y. Sun, M. Pistoia, and Y. Alexeev, A survey of quantum computing for finance, arXiv preprint arXiv:2201.02773 (2022).
  * Cheng  _et al._ [2023] B. Cheng _et al._ , Noisy intermediate-scale quantum computers, Front. Phys. (Beijing) 18, 21308 (2023), arXiv:2303.04061 [quant-ph] .
  * Alexeev  _et al._ [2021] Y. Alexeev _et al._ , Quantum Computer Systems for Scientific Discovery, PRX Quantum 2, 017001 (2021), arXiv:1912.07577 [quant-ph] .
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Finnila  _et al._ [1994] A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll, Quantum annealing: A new method for minimizing multidimensional functions, Chemical physics letters 219, 343 (1994).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse ising model, Physical Review E 58, 5355 (1998).
  * Morita and Nishimori [2008] S. Morita and H. Nishimori, Mathematical foundation of quantum annealing, Journal of Mathematical Physics 49 (2008).
  * Santoro and Tosatti [2006] G. E. Santoro and E. Tosatti, Optimization using quantum mechanics: quantum annealing through adiabatic evolution, Journal of Physics A: Mathematical and General 39, R393 (2006).
  * Das and Chakrabarti [2008] A. Das and B. K. Chakrabarti, Colloquium: Quantum annealing and analog quantum computation, Rev. Mod. Phys. 80, 1061 (2008), arXiv:0801.2193 [quant-ph] .
  * Baldassi and Zecchina [2018] C. Baldassi and R. Zecchina, Efficiency of quantum vs. classical annealing in nonconvex learning problems, Proceedings of the National Academy of Sciences 115, 1457 (2018).
  * Hauke  _et al._ [2020] 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, 054401 (2020).
  * Abel  _et al._ [2022] S. Abel, A. Blance, and M. Spannowsky, Quantum optimization of complex systems with a quantum annealer, Phys. Rev. A 106, 042607 (2022), arXiv:2105.13945 [quant-ph] .
  * King  _et al._ [2023] A. D. King _et al._ , Quantum critical dynamics in a 5,000-qubit programmable spin glass, Nature 617, 61 (2023), arXiv:2207.13800 [quant-ph] .
  * Rajak  _et al._ [2023] A. Rajak, S. Suzuki, A. Dutta, and B. K. Chakrabarti, Quantum annealing: An overview, Philosophical Transactions of the Royal Society A 381, 20210417 (2023).
  * McGeoch and Farré [2020] C. McGeoch and P. Farré,  _The D-wave Advantage system: An overview_ , Tech. Rep. 14-1049A-A (D-Wave Systems Inc., Burnaby, BC, Canada, 2020).
  * McGeoch and Farré [2021] C. McGeoch and P. Farré,  _The advantage system: Performance update_ , Tech. Rep. 14-1054A-A (D-Wave Systems Inc., Burnaby, BC, Canada, 2021).
  * McGeoch and Farré [2022] C. McGeoch and P. Farré,  _Advantage processor overview_ , Tech. Rep. 14-1058A-A (D-Wave systems, 2022).
  * Krakoff  _et al._ [2021] B. Krakoff, S. M. Mniszewski, and C. F. A. Negre, A QUBO algorithm to compute eigenvectors of symmetric matrices, (2021), arXiv:2104.11311 .
  * Krakoff  _et al._ [2022] B. Krakoff, S. M. Mniszewski, and C. F. Negre, Controlled precision qubo-based algorithm to compute eigenvectors of symmetric matrices, Plos one 17, e0267954 (2022).
  * Illa and Savage [2022] M. Illa and M. J. Savage, Basic elements for simulations of standard-model physics with quantum annealers: Multigrid and clock states, Phys. Rev. A 106, 052605 (2022), arXiv:2202.12340 [quant-ph] .
  * Tosti Balducci  _et al._ [2022] G. Tosti Balducci, B. Chen, M. Möller, M. Gerritsma, and R. De Breuker, Review and perspectives in quantum computing for partial differential equations in structural mechanics, Frontiers in Mechanical Engineering , 75 (2022).
  * Stockley and Briggs [2023] A. Stockley and K. Briggs, Optimizing antenna beamforming with quantum computing, in  _2023 17th European Conference on Antennas and Propagation (EuCAP)_ (IEEE, 2023) pp. 1–5.
  * Salpeter and Bethe [1951] E. E. Salpeter and H. A. Bethe, A Relativistic Equation for Bound-State Problems, Phys. Rev. 84, 1232 (1951).
  * Gell-Mann and Low [1951] M. Gell-Mann and F. Low, Bound states in quantum field theory, Phys. Rev. 84, 350 (1951).
  * Itzykson and Zuber [2012] C. Itzykson and J.-B. Zuber,  _Quantum field theory_ (Courier Corporation, 2012).
  * Nakanishi [1971] N. Nakanishi,  _Graph Theory and Feynman Integrals_ (Gordon and Breach, New York, 1971).
  * Frederico  _et al._ [2014] T. Frederico, G. Salmè, and M. Viviani, Quantitative studies of the homogeneous Bethe-Salpeter equation in Minkowski space, Phys. Rev. D 89, 016010 (2014).
  * Alliney  _et al._ [1992] S. Alliney, F. Laudiero, and M. Savoia, A variational technique for the computation of the vibration frequencies of mechanical systems governed by nonsymmetric matrices, Applied mathematical modelling 16, 148 (1992).
  * Kim  _et al._ [2007] N.-I. Kim, D. K. Shin, and M.-Y. Kim, Improved flexural–torsional stability analysis of thin-walled composite beam and exact stiffness matrix, International Journal of Mechanical Sciences 49, 950 (2007).
  * Li and Vuik [2004] C. Li and C. Vuik, Eigenvalue analysis of the simple preconditioning for incompressible flow, Numerical Linear Algebra with Applications 11, 511 (2004).
  * Scheben and Graham [2011] F. Scheben and I. G. Graham, Iterative methods for neutron transport eigenvalue problems, SIAM Journal on Scientific Computing 33, 2785 (2011).
  * Burroughs  _et al._ [2004] E. A. Burroughs, L. A. Romero, R. B. Lehoucq, and A. G. Salinger, Linear stability of flow in a differentially heated cavity via large-scale eigenvalue calculations, International Journal of Numerical Methods for Heat & Fluid Flow 14, 803 (2004).
  * Orginos and Richards [2015] K. Orginos and D. Richards, Improved methods for the study of hadronic physics from lattice QCD, J. Phys. G 42, 034011 (2015).
  * Yarkoni  _et al._ [2022] S. Yarkoni, E. Raponi, T. Bäck, and S. Schmitt, Quantum annealing for industry applications: introduction and review, Rept. Prog. Phys. 85, 104001 (2022), arXiv:2112.07491 [quant-ph] .
  * Johnson  _et al._ [2011] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, _et al._ , Quantum annealing with manufactured spins, Nature 473, 194 (2011).
  * Nakanishi [1963] N. Nakanishi, Partial-Wave Bethe-Salpeter Equation, Phys. Rev. 130, 1230 (1963).
  * Roberts and Williams [1994] C. D. Roberts and A. G. Williams, Dyson-Schwinger equations and their application to hadronic physics, Prog. Part. Nucl. Phys. 33, 477 (1994), arXiv:hep-ph/9403224 [hep-ph] .
  * Maris and Roberts [2003] P. Maris and C. D. Roberts, Dyson-Schwinger equations: A Tool for hadron physics, Int. J. Mod. Phys. E12, 297 (2003), arXiv:nucl-th/0301049 [nucl-th] .
  * Bashir  _et al._ [2012] A. Bashir, L. Chang, I. C. Cloët, B. El-Bennich, Y.-X. Liu, C. D. Roberts, and P. C. Tandy, Collective perspective on advances in Dyson-Schwinger Equation QCD, Commun. Theor. Phys. 58, 79 (2012), arXiv:1201.3366 [nucl-th] .
  * Osterwalder and Schrader [1973] K. Osterwalder and R. Schrader, Axioms for Euclidean Green’s functions, Commun. Math. Phys. 31, 83 (1973).
  * Osterwalder and Schrader [1975] K. Osterwalder and R. Schrader, Axioms for Euclidean Green’s Functions. 2., Commun. Math. Phys. 42, 281 (1975).
  * Glover  _et al._ [2018] F. Glover, G. Kochenberger, and Y. Du, A tutorial on formulating and using qubo models, arXiv preprint arXiv:1811.11538 (2018).
  * Teplukhin  _et al._ [2019] A. Teplukhin, B. K. Kendrick, and D. Babikov, Calculation of molecular vibrational spectra on a quantum annealer, Journal of chemical theory and computation 15, 4555 (2019).
  * Pastorello and Blanzieri [2019] D. Pastorello and E. Blanzieri, Quantum annealing learning search for solving qubo problems, Quantum Information Processing 18, 1 (2019).
  * Teplukhin  _et al._ [2020] A. Teplukhin, B. K. Kendrick, and D. Babikov, Solving complex eigenvalue problems on a quantum annealer with applications to quantum scattering resonances, Phys. Chem. Chem. Phys. 22, 26136 (2020).
  * Chang  _et al._ [2019] C. C. Chang, A. Gambhir, T. S. Humble, and S. Sota, Quantum annealing for systems of polynomial equations, Sci. Rep. 9, 10258 (2019), arXiv:1812.06917 [quant-ph] .
  * Wolkowicz and Styan [1980] H. Wolkowicz and G. P. Styan, Bounds for eigenvalues using traces, Linear algebra and its applications 29, 471 (1980).
  * Varga [2010] R. S. Varga,  _Geršgorin and his circles_ , Vol. 36 (Springer Science & Business Media, 2010).
  * Wick [1954] G. C. Wick, Properties of Bethe-Salpeter Wave Functions, Phys. Rev. 96, 1124 (1954).
  * Cutkosky [1954] R. E. Cutkosky, Solutions of a Bethe-Salpeter equations, Phys. Rev. 96, 1135 (1954).
  * Dirac [1949] P. A. M. Dirac, Forms of Relativistic Dynamics, Rev. Mod. Phys. 21, 392 (1949).
  * Fan [1950] K. Fan, On a theorem of $\rm{Weyl}$ concerning eigenvalues of linear transformation: $\rm{II}$, Proc. Nat. Acad. Sci. 36, 31 (1950).