# Integrability in the chiral model of magic angles

Simon Becker<EMAIL_ADDRESS>ETH Zurich, Institute for Mathematical
Research, Rämistrasse 101, 8092 Zurich, Switzerland  ,  Tristan Humbert
<EMAIL_ADDRESS>ENS Paris, Département de Mathématiques et
Applications, Rue d’Ulm, Paris, France  and  Maciej Zworski
<EMAIL_ADDRESS>Department of Mathematics, University of California,
Berkeley, CA 94720, USA.

###### Abstract.

Magic angles in the chiral model of twisted bilayer graphene are parameters
for which the chiral version of the Bistritzer–MacDonald Hamiltonian exhibits
a flat band at energy zero. We compute the sums over powers of (complex) magic
angles and use that to show that the set of magic angles is infinite. We also
provide a new proof of the existence of the first real magic angle, showing
also that the corresponding flat band has minimal multiplicity for the
simplest possible choice of potentials satisfying all symmetries. These
results indicate (though do not prove) a hidden integrability of the chiral
model.

##  1\. Introduction and statement of results

When two sheets of graphene are stacked on top of each other and twisted, it
has been observed that at certain angles, coined the _magic angles_ , the
composite system becomes superconducting. In this article, we study the chiral
limit of the Bistritzer-MacDonald Hamiltonian [BiMa11, CGG22, Wa*22]

$H(\alpha)=\begin{pmatrix}0&D(\alpha)^{*}\\\ D(\alpha)&0\end{pmatrix}\text{
with }D(\alpha)=\begin{pmatrix}D_{\bar{z}}&\alpha U(z)\\\ \alpha
U(-z)&D_{\bar{z}}\end{pmatrix}$

where the parameter $\alpha$ is proportional to the inverse relative twisting
angle. After a simple rescaling, the potential is a smooth and periodic
function satisfying

$U(z+a_{\ell})=\bar{\omega}U(z),\quad U(\omega z)=\omega U(z),\text{ and
}U(\bar{z})=\overline{U(z)},$ (1.1)

where $\omega=e^{2\pi i/3}$ and $a_{\ell}=\frac{4}{3}\pi i\omega^{\ell}.$ The
simplest example of such a potential and our canonical choice of $U$ is

$U_{0}(z)=\sum_{k=0}^{2}\omega^{k}e^{\frac{1}{2}(z\bar{\omega}^{k}-\bar{z}\omega^{k})}.$
(1.2)

Even though the potential $U(z)$ is only periodic with respect to $\Gamma=4\pi
i(\omega{\mathbb{Z}}\oplus\omega^{2}{\mathbb{Z}})$ the first property implies
that the matrix potential, and thus $D(\alpha)$, commutes with the translation
operator

$\mathscr{L}_{a}w(z):=\begin{pmatrix}\omega^{a_{1}+a_{2}}&0\\\
0&1\end{pmatrix}w(z+a),\ \ \ a\in\tfrac{1}{3}\Gamma,$ (1.3)

where $w\in\mathbb{C}^{2}$ and $a=\tfrac{4}{3}\pi i(\omega
a_{1}+\omega^{2}a_{2}),\ \ a_{j}\in\mathbb{Z}.$ We note that if $\Gamma^{*}$
is the dual (reciprocal) lattice of $\Gamma$, then $3\Gamma^{*}$ is the dual
lattice of $\frac{1}{3}\Gamma$.

Figure 1.  Plots of the first 5 non-negative eigenvalues of
$H_{\mathbf{k}}(0.3)$ acting on $\mathscr{H}$ (see (1.4) and (1.5)), as
function of $k=(k_{1}\omega^{2}-k_{2}\omega)/\sqrt{3}$ in in a fundamental
cell of $3\Gamma^{*}$, parametrized by $(k_{1},k_{2})$ $|k_{j}|<\frac{3}{2}$.
See also [BHZ22, Figure 4] for more information and comparison with band
structure of other models.

When moving to functions with values in
$\mathbb{C}^{4}=\mathbb{C}^{2}\times\mathbb{C}^{2}$ (on which $H(\alpha)$
acts) we extend the action of $\mathscr{L}_{\mathbf{a}}$ to an action on each
$\mathbb{C}^{2}$ component. We then consider the Floquet spectrum of

$H_{k}(\alpha)=\begin{pmatrix}0&D(\alpha)^{*}-\bar{k}\\\
D(\alpha)-k&0\end{pmatrix}\text{ with }k\in 3\Gamma^{*},$ (1.4)

defined by $(H_{k}(\alpha)-E_{j}(\alpha,k))w_{j}(\alpha,k)=0$, where
eigenvalues of positive energy are labelled with $j\geq 1$ in ascending order,
as a self-adjoint operator on

$\mathscr{H}:=\\{\mathbf{v}\in
L^{2}(\mathbb{C}/\Gamma):\mathscr{L}_{\mathbf{a}}\mathbf{v}=\mathbf{v},\ \
\mathbf{a}\in\tfrac{1}{3}\Gamma\\},$ (1.5)

with the domain given by $\mathscr{H}\cap H^{1}(\mathbb{C}/\Gamma)$ such that

$\operatorname{Spec}_{L^{2}(\mathbb{C};\mathbb{C}^{4})}(H(\alpha))=\bigcup_{k\in{\mathbb{C}}}\operatorname{Spec}_{\mathscr{H}}(H_{k}(\alpha)).$

This Hamiltonian is an effective one-particle model which exhibits perfectly
flat bands at magic angles. This appearance of perfectly flat bands in the
chiral limit was explained by Tarnopolsky, Kruchkov and Vishwanath [TKV19]
with the help of Jacobi theta functions111As was pointed out to us by Alex
Sobolev a similar argument appeared in the work of Dubrovin and Novikov
[DuNo80] who studied magnetic Hamiltonians on tori.. An equivalent spectral
theoretic characterization of magic angles was then provided in [Be*22]: if we
define the following compact Birman-Schwinger operator

$T_{k}=(2D_{\bar{z}}-k)^{-1}\begin{pmatrix}0&U(z)\\\ U(-z)&0\end{pmatrix}.$
(1.6)

then (see [Be*22, Theorem 2] we have the following equivalence [BHZ22, §2.3])

$0\in\bigcap_{k\in{\mathbb{C}}}\operatorname{Spec}_{\mathscr{H}}(H_{k}(\alpha))\
\Longleftrightarrow\
\left\\{\begin{array}[]{l}\alpha^{-1}\in\operatorname{Spec}_{\mathscr{H}}(T_{k_{0}})\\\
\text{ for some
}k_{0}\in{\mathbb{C}}\setminus(3\Gamma^{*}-\\{0,i\\}),\end{array}\right.$
(1.7)

where $\mathscr{H}$ is defined in (1.5). In other words, the spectrum of
$T_{k_{0}}$ is independent of
$k_{0}\in{\mathbb{C}}\setminus(3\Gamma^{*}-\\{0,i\\})$ and characterizes the
values of $\alpha\in{\mathbb{C}}$ at which the Hamiltonian exhibits a flat
band at zero energy. Since the parameter $\alpha$ is inherently connected with
the twisting angle, we shall refer to $\alpha$’s at which (1.7) occurs as
_magic_ and denote their set by $\mathcal{A}\subset\mathbb{C}$.

Figure 2.  The set $\mathcal{A}$ of magic $\alpha$’s for which (1.7) holds,
that is, the first band is flat. The positive elements of $\mathcal{A}$ are
the reciprocals of the “physically relevant” positive angles. Potential (1.2)
is responsible for the regularity of the set which seems to indicate hidden
integrability. For more general potentials the distribution is more
complicated – see https://math.berkeley.edu/~zworski/multi.mp4 for
$U_{\theta}(z)=(\cos^{2}\theta)U(z)+(\sin^{2}\theta)\sum_{k=0}^{2}\omega^{k}e^{\bar{z}\omega^{k}-z\bar{\omega}^{k}}$
which satisfies the required symmetries (1.1). The animation also indicates
changing multiplicities.

The analysis of magic angles is therefore reduced to a spectral theory problem
involving a single compact non self-adjoint operator. Since even non-trivial
non self-adjoint compact operators do not necessarily have non-zero
eigenvalues, the existence of a parameter $\alpha$ at which the Hamiltonian
exhibits a flat band at zero energy is non-trivial. In [Be*22] the existence
of such a complex parameter $\alpha\in{\mathbb{C}}\setminus\\{0\\}$ was first
concluded by showing that
$\operatorname{tr}_{\mathscr{H}}(T_{k}^{4})=8{\pi}/{\sqrt{3}}$ which implied
existence of a non-zero eigenvalue222In [Be*22] we considered the trace on
$L^{2}(\mathbb{C}/\Gamma;\mathbb{C}^{2})$ which gave this answer multiplied by
$9$.. This result was improved by a computer-assisted proof [WaLu21] in which
Watson and Luskin used the complex-analytic characterization of magic angles
from [TKV19] to prove existence of the first _real_ magic angle and obtained
explicit bounds on its position.

In this article, we exhibit a general form of traces of powers of $T_{k}$.
This suggests a hidden integrability of the Hamiltonian $H(\alpha)$ for
potentials satisfying (1.1), as all traces exhibits special arithmetic
properties. With our current techniques, we do not have explicit control on
the full set of traces which would imply a complete understanding of all magic
angles. These are already visible in the regular but evasive structure of the
set of of magic $\alpha$, $\mathcal{A}\subset\mathbb{C}$ – see Figure 2.

######  Theorem 1.

For $\ell\geq 2$ and $U=U_{0}$ with $U_{0}$ as in (1.2)

$\operatorname{tr}(T_{k}^{2\ell})=\sum_{\alpha\in\mathcal{A}}\alpha^{-2\ell}=\frac{\pi}{\sqrt{3}}\,q_{\ell}\
\text{ with }q_{\ell}\in\mathbb{Q}.$ (1.8)

In addition, we are able to express the rational numbers
$q_{\ell}\in\mathbb{Q}$ in terms of a finite sum involving residues of
rational functions which is fully presented in Theorem 4. A generalization of
Theorem 1 which extends this result to more general potentials $U$ (1.1) is
presented in Theorem 5. As we show in §6, it is already possible to conclude
directly from Theorem 1 that

######  Theorem 2.

Let $U=U_{0}$ with $U_{0}$ as in (1.2). There exist infinitely many magic
$\alpha$’s, that is,

$|\mathcal{A}|=\infty.$

This theorem will follow from the more general Theorem 6 and the observation
that by the aforementioned explicit computation
$\operatorname{tr}_{\mathscr{H}}(T_{k}^{4})=8{\pi}/{\sqrt{3}}$ for $U=U_{0}$
there is at lest one complex magic angle. We then focus on real magic angles.
Since the operator $T_{k}^{2}$ is Hilbert-Schmidt, we can use the regularized
determinant to study real magic $\alpha.$ Compared with the initial approach
proposed in [TKV19], this approach has two advantages. Unlike the series
expansion in [TKV19, WaLu21], the regularized determinant is an entire
function with explicit error bounds in terms of the Hilbert-Schmidt norm. In
addition, the Taylor coefficients of the determinant are polynomials of traces
as in Theorem 1. This leads to

######  Theorem 3.

The chiral Hamiltonian with $U=U_{0}$ and $U_{0}$ as in (1.2), exhibits a flat
band of multiplicity 2 at a real magic $\alpha_{*}\in(0.583,0.589)$, which is
minimal, in the sense that the Hamiltonian does not possess a flat band for
any $\alpha$ satisfying $|\alpha|<|\alpha_{*}|$, that is,

$|\mathcal{A}\cap(0.583,0.589)|=1,\ \ \ \mathcal{A}\cap
D_{\mathbb{C}}(0,\alpha^{*})=\emptyset,$

where the counting $|\bullet|$ respects multiplicities. In particular, the
flat bands of multiplicity 2 are uniformly gapped from all other bands.

Remark. Compared with results in [WaLu21] which require floating-point
arithmetic, our proof of existence relies only on exact symbolic computations,
the exact evaluation of residues to compute traces of powers of $T_{k}$ and
the summation of finitely many matrix entries to estimate the Hilbert-Schmidt
norm.

##  2\. Preliminaries

From now on, we consider a potential $U\in
C^{\infty}(\mathbb{C}/\Gamma;\mathbb{C})$ satisfying the first two symmetries
of (1.1). The last symmetry $U(\bar{z})=\overline{U(z)}$ will only be needed
in Corollary 5 to ensure that all traces are real. We recall that an
orthonormal basis of $L^{2}(\mathbb{C}/\Gamma;\mathbb{C})$ is given by setting

$e_{\nu}(z):=e^{\frac{i}{2}(\bar{\nu}z+\nu
z)}/\sqrt{\mathrm{Vol}(\mathbb{C}/\Gamma)},\quad\nu\in\Gamma^{*}=\frac{1}{\sqrt{3}}\big{(}\mathbb{Z}+\omega\mathbb{Z}\big{)}.$

We can express the potential $U$ in this basis. A straightforward calculation
gives

######  Proposition 2.1.

Let $a=\frac{4\pi i}{3}(\omega a_{1}+\omega^{2}a_{2})\in\Gamma/3$,
$a_{j}\in\mathbb{Z}$. Then $u\in L^{2}(\mathbb{C}/\Gamma)$ satisfies

$u(z+a)=\bar{\omega}^{(a_{1}+a_{2})}u(z),\quad u(\omega z)={\omega}u(z),$

if and only if

$u(z)=\sum_{n\in\mathbb{Z}^{2}}c_{n}e_{k_{n}}$ (2.1)

where for $n\in{\mathbb{Z}}^{2}$

$c_{n}=(u,e_{k_{n}})_{L^{2}(\mathbb{C}/\Gamma)},\quad
k_{n}={\frac{\omega^{2}(2+3n_{1})-\omega(2-3n_{2})}{\sqrt{3}}},$

satisfies

$c_{n}=\omega
c_{(-n_{2})(n_{1}-n_{2}-1)}=\omega^{2}c_{(n_{2}-n_{1}+1)(-n_{1})}.$

If in addition $\overline{u(\bar{z})}=u(z)$ then

$\overline{c_{n}}=c_{(-n_{2})(-n_{1})}=\omega
c_{n_{1}(n_{1}-n_{2}-1)}=\omega^{2}c_{(n_{2}-n_{1}-1)n_{2}}.$

Our aim is to obtain trace formulae for the powers of the compact operator
$T_{k}$ defined by in (1.6). Since odd powers of $T_{k}$ have only off-
diagonal components, it is clear that the traces of odd powers vanish. Thus,
it is sufficient to compute the traces of powers of

$T_{k}^{2}=\begin{pmatrix}(2D_{\bar{z}}-k)^{-1}U(z)(2D_{\bar{z}}-k)^{-1}U(-z)&0\\\
0&(2D_{\bar{z}}-k)^{-1}U(-z)(2D_{\bar{z}}-k)^{-1}U(z)\end{pmatrix}.$ (2.2)

The invariance of the trace under cyclic permutations shows that it is
sufficient to compute traces of powers of

$\mathcal{A}_{k}:=(2D_{\bar{z}}-k)^{-1}U(z)(2D_{\bar{z}}-k)^{-1}U(-z):L^{2}(\mathbb{C}/\Gamma;\mathbb{C})\to
L^{2}(\mathbb{C}/\Gamma;\mathbb{C}),\quad k\notin\Gamma^{*}.$ (2.3)

We shall study traces of powers of $\mathcal{A}_{k}$ on smaller $L^{2}$
spaces, which we define below, for $(p_{1},p_{2})\in\mathbb{Z}_{3}^{2}$, by

$\begin{split}L^{2}_{(p_{1},p_{2})}(\mathbb{C}/\Gamma;\mathbb{C})&:=\Big{\\{}u\in
L^{2}:u(z+2i(\omega
a_{1}+a_{2}\omega^{2}))=e^{i(a_{1}p_{1}+a_{2}p_{2})}u(z);\\\ &\qquad\ \
a_{j}\in\tfrac{2\pi}{3}\mathbb{Z}\Big{\\}}\end{split}$ (2.4)

whose ${\mathbb{C}}^{2}$-valued analogues are defined, using (1.3), as

$\begin{split}L^{2}_{(p_{1},p_{2})}(\mathbb{C}/\Gamma;\mathbb{C}^{2})&:=\Big{\\{}u\in
L^{2}:\mathscr{L}_{a}u(z)=e^{i(a_{1}p_{1}+a_{2}p_{2})}u(z);a_{j}\in\tfrac{2\pi}{3}\mathbb{Z}\Big{\\}}.\end{split}$
(2.5)

We remark that the operator $(2D_{\bar{z}}-k)^{-1}$ acts diagonally on the
Fourier basis and thus preserves the $L^{2}_{(p_{1},p_{2})}$ spaces. On the
other hand, multiplication by $U(\pm z)$ does not preserve the space but one
has by the translational symmetry defined in (1.1)

$U(\pm z):L^{2}_{(p_{1},p_{2})}\to L^{2}_{(p_{1}\mp 1,p_{2}\mp
1)},\quad(p_{1},p_{2})\in\mathbb{Z}_{3}^{2}.$

In total, we have

$L^{2}_{(p_{1},p_{2})}\xrightarrow{U(-z)}L^{2}_{(p_{1}+1,p_{2}+1)}\xrightarrow{(2D_{\bar{z}}-k)^{-1}}L^{2}_{(p_{1}+1,p_{2}+1)}\xrightarrow{U(z)}L^{2}_{(p_{1},p_{2})}\xrightarrow{(2D_{\bar{z}}-k)^{-1}}L^{2}_{(p_{1},p_{2})}.$

This shows that we can restrict the operator $\mathcal{A}_{k}$ to the
subspaces $L^{2}_{(p_{1},p_{2})}$. From now on, we will denote by $A_{k}$ the
restriction of $\mathcal{A}_{k}$ to $L^{2}_{(1,1)}$. We then define the
unitary multiplication operator

$\begin{gathered}U_{(p_{1},p_{2})}:L^{2}_{(0,0)}(\mathbb{C}/\Gamma,\mathbb{C}^{2})\to
L^{2}_{(p_{1},p_{2})}(\mathbb{C}/\Gamma,\mathbb{C}^{2}),\\\
U_{(p_{1},p_{2})}v(z):=e^{\frac{i}{2}(z\bar{p}+\bar{z}p)}v(z),\ \
p=\tfrac{1}{\sqrt{3}}(\bar{\omega}p_{1}-\omega p_{2}),\ \ \
p_{j}\in\mathbb{Z}_{3},\\\
U_{(p_{1},p_{2})}T_{k}U^{*}_{(p_{1},p_{2})}=T_{k-p},\ \ \
k\notin\Gamma^{*}.\end{gathered}$

The $k$-independence of the spectrum of $T_{k}$ implies then

$\begin{split}\operatorname{Spec}_{L^{2}_{(0,0)}(\mathbb{C}/\Gamma,\mathbb{C}^{2})}(T_{k}^{2})\setminus\\{0\\}&=\operatorname{Spec}_{L^{2}_{(1,1)}(\mathbb{C}/\Gamma,\mathbb{C}^{2})}(T_{k}^{2})\setminus\\{0\\}\\\
&=\operatorname{Spec}_{L^{2}_{(1,1)}(\mathbb{C}/\Gamma,\mathbb{C})}(A_{k})\setminus\\{0\\},\end{split}$
(2.6)

where $k\in D(0,r)\setminus\\{0\\}$, and the last equality is meant in the
sense of sets: multiplicities of elements in the top row are twice the
multiplicities of elements in the bottom row.
We also note that $A_{k}$ is defined for $k\in D(0,r)$ since
$D_{\bar{z}}^{-1}$ is defined on $L^{2}_{(p,p)}$, $p\not\equiv 0\\!\\!\\!\mod
3$. Since $\mathbb{C}\ni k\to\mathcal{A}_{k}|_{L^{2}_{(1,1)}}$ is an analytic
family of operators with compact resolvents and the spectrum is independent of
${k}\in{\mathbb{C}}\setminus 3\Gamma^{*}$, it follows that
$\operatorname{Spec}(A_{k})=\operatorname{Spec}(A_{0})$ [Ka80, Theorem
$1.10$]. From (2.6) we obtain, as sets,

$\operatorname{Spec}_{L^{2}_{(0,0)}}(T_{p}^{2})\setminus\\{0\\}=\operatorname{Spec}_{L^{2}_{(1,1)}}(A_{k})\setminus\\{0\\},\
\ \ p\in D(0,r)\setminus\\{0\\},\ \ k\in D(0,r),$ (2.7)

with multiplicities on the left, twice the multiplicities on the right. Since
$k=0$ is included in the set of possible $k$ for $A_{k}$. Indeed, the set of
possible values of $k$ is
$\mathbb{C}\setminus\big{(}(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i)\big{)}$. We
conclude together with [BHZ22, Theorem $2$] that

$\operatorname{dim}\ker_{\mathscr{H}}(D(\alpha))=\operatorname{dim}\ker_{L^{2}_{(1,1)}}(A_{0}-\alpha^{-2}).$
(2.8)

We end this preliminary section by stating and proving the main three
properties we will use for our calculation.

######  Lemma 2.2.

Consider a potential $U\in C^{\infty}(\mathbb{C}/\Gamma;\mathbb{C})$
satisfying the first two symmetries of (1.1) with a finite number of non zero
Fourier mode in its decomposition (2.1). Define the operator $A_{k}$ for
$k\notin(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i)$, where
$i:=-(\omega^{2}-\omega)/\sqrt{3}$, to be the restriction of $\mathcal{A}_{k}$
defined in (2.3) on the space $L^{2}_{(1,1)}$. For $\ell\geq 2$, one has:

  * •

The trace is constant in $k$

$\operatorname{tr}(A_{k}^{\ell})=\tau_{\ell}\text{ independent of
}k\in\mathbb{C}\setminus(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i).$ (2.9)

  * •

The function $\mathbb{C}\setminus(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i)\ni
k\mapsto\langle A_{k}^{\ell}e_{m},e_{m}\rangle_{L^{2}}$ is a finite sum of
rational fractions on the complex plane $\mathbb{C}$ with degree equal to
$-2\ell$ and with $($a finite number of$)$ poles contained in
$(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i)$.

  * •

For any $\gamma\in\Gamma^{*}$ and for any
$k\notin(3\Gamma^{*}-i)\cup(3\Gamma^{*}-2i)$, we have

$\langle A_{k}^{\ell}e_{3\gamma+i},e_{3\gamma+i}\rangle_{L^{2}}=\langle
A_{k-3\gamma}^{\ell}e_{i},e_{i}\rangle_{L^{2}}.$

###### Proof.

The first point is a consequence of the independence of the spectrum of
$T_{k}$ in $k$ (see [BHZ22, §2.3]) as well as the relation (2.7).
For the last two points, we prove by induction that $k\mapsto
A_{k}^{\ell}e_{3\gamma+i}$, where $\gamma\in\Gamma^{*}$, is of the form

$A_{k}^{\ell}e_{3\gamma+i}=\sum_{\nu\in F}R_{\nu+3\gamma}(k)e_{\nu+3\gamma},$
(2.10)

where $F\subset(3\Gamma^{*}+i)$ is a finite set and $R_{\nu}(k)$ is a sum of
rational fraction of degree $-2\ell$ with poles located on
$(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i)$. Moreover, we will prove that the one has
the relation $R_{\nu+3\gamma}(k)=R_{\nu}(k-3\gamma).$
The result is clear for $\ell=0$. Suppose the result true for $\ell$, let’s
prove it holds for $\ell+1$. The main observation is that multiplication by
$U(\pm z)$ acts as a shift on the Fourier basis. The multiplication by $U(-z)$
sends $e_{\nu}$ to a linear combination of $e_{\ell}$ for
$\ell\in(3\Gamma^{*}+2i)$. Then applying $(D(0)-k)^{-1}$ multiplies the
coefficient of $e_{\ell}$ by $(\ell-k)^{-1}$. Multiplying by $U(z)$ gives back
a linear combination of $e_{\nu}$ with $\nu\in(3\Gamma^{*}+i)$. Finally,
applying $(D(0)-k)^{-1}$ multiplies the coefficient of $e_{\nu}$ by
$(\nu-k)^{-1}$. This means that, using the induction hypothesis (2.10),

$A_{k}^{\ell+1}e_{3\gamma+i}=\sum_{\nu\in F}R_{\nu+3\gamma}(k)\sum_{\eta\in
L}\sum_{\beta\in
L}\frac{a_{\eta}}{k-(\nu+\beta-\eta+3\gamma)}\frac{a_{\beta}}{k-(\nu+\beta+3\gamma)}e_{\nu+\beta-\eta+3\gamma},$

where $L\subset 3\Gamma^{*}+i$ is a finite subset that depends only on $U$ and
$a_{\bullet}$ are constants. Thus, it is clear from this formula that the
induction carries on to $\ell+1$. This concludes the proof of the Lemma. ∎

##  3\. Trace computations

We prove the following result.

######  Theorem 4.

Let $A_{k}:L^{2}_{(1,1)}\to L^{2}_{(1,1)}$ be a meromorphic family of Hilbert-
Schmidt operators defined for $k\notin(3\Gamma^{*}-i)\cup(3\Gamma^{*}+i)$. We
suppose that $A_{k}$ satisfies the three properties stated in Lemma 2.2. Then
one has, for any $\ell\geq 2$,

$\begin{split}\tau_{\ell}&=\frac{2i\pi\omega}{3\sqrt{3}}\sum_{n\in\mathbb{Z}}n\left[\sum_{m\in\mathbb{Z}}\mathrm{Res}\big{(}\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}},\sqrt{3}(m\omega^{2}-n\omega)+i\big{)}\right.\\\
&\ \ \ \ \ \ \ \ \ \ \ \
\left.+\sum_{m\in\mathbb{Z}}\mathrm{Res}\big{(}\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}},\sqrt{3}(m\omega^{2}-n\omega)+2i\big{)}\right],\end{split}$
(3.1)

where the all infinite sums are in fact finite.

###### Proof.

We want to give a semi-explicit formula for $\tau_{\ell}$ in terms of the
residue of the rational fraction (second point in Lemma 2.2)

$k\in\mathbb{C}\setminus\big{(}3\Gamma^{*}-i\big{)}\cup\big{(}3\Gamma^{*}+i\big{)}\mapsto\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}.$

We first start by writing, using that $A_{k}^{\ell}$ is trace-class for
$\ell\geq 2$,

$\tau_{\ell}=\sum_{\gamma\in\Gamma}\langle
A_{k}^{\ell}e_{3\gamma+i},e_{3\gamma+i}\rangle_{L^{2}}.$

We start with the relation which follows directly from (2.9),

$3\tau_{\ell}=\int_{0}^{3}\operatorname{tr}\big{(}A^{\ell}_{t\omega^{2}/\sqrt{3}}\big{)}dt=\int_{0}^{3}\sum_{n\in\mathbb{Z}^{2}}\langle
A_{t\omega^{2}/\sqrt{3}}^{\ell}e_{3\gamma_{n}+i},e_{3\gamma_{n}+i}\rangle_{L^{2}}dt.$

Here, we wrote
$\gamma_{n}:=\frac{1}{\sqrt{3}}(n_{1}\omega^{2}-n_{2}\omega)\in\Gamma^{*}$. We
now use the third property stated in Lemma 2.2 to write

$3\tau_{\ell}=\int_{0}^{3}\sum_{n\in\mathbb{Z}^{2}}\langle
A_{t\omega^{2}/\sqrt{3}-3\gamma_{n}}^{\ell}e_{i},e_{i}\rangle_{L^{2}}dt.$

The second property in Lemma 2.2 implies that

$\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}={\mathcal{O}}\big{(}k^{-2\ell}\big{)}.$
(3.2)

Since we assume that $\ell\geq 2$, this justifies the exchange of integration
and summation such that

$3\tau_{\ell}=\sum_{n\in\mathbb{Z}^{2}}\int_{0}^{3}\langle
A^{\ell}_{(t+3n_{1})\omega^{2}/\sqrt{3}-3n_{2}\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}dt.$

We make the change of variable $s=t+3n_{1}$ and sum in $n_{1}$ to get

$3\tau_{\ell}=\sum_{n\in\mathbb{Z}}\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3n\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}dt.$
(3.3)

We now consider

$\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3n\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}\
dt-\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3(n+1)\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}dt,\
\ n\in\mathbb{Z}.$

This is equal to the limit of the integral over a parallelogram $\Gamma_{n,R}$
with sides $\frac{1}{\sqrt{3}}[-R\omega^{2}-3n\omega,R\omega^{2}-3n\omega]$,
$\frac{1}{\sqrt{3}}[R\omega^{2}-3n\omega,R\omega^{2}-3(n+1)\omega]$,
$\frac{1}{\sqrt{3}}[R\omega^{2}-3(n+1)\omega,-R\omega^{2}-3(n+1)\omega]$ and
$\frac{1}{\sqrt{3}}[-R\omega^{2}-3(n+1)\omega,-R\omega^{2}-3n\omega]$. Here,
we used (3.2) to prove that the integral over the small parallel sides tends
to $0$. In particular, because there is only a finite number of poles, we see
that, for $|n|$ large enough, one has

$\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3n\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}dt-\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3(n+1)\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}dt=0.$

Using formula (3.3) as well as a partial summation, this allows us to rewrite
the full trace as a telescopic sum

$3\tau_{\ell}=\sum_{n\in\mathbb{Z}}n\left[\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3n\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}\
dt-\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3(n+1)\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}dt\right].$
(3.4)

The residue theorem shows that for $n\in\mathbb{Z}$ and $R$ large enough,

$\begin{split}\int_{\Gamma_{n,R}}\langle
A^{\ell}_{z}e_{i},e_{i}\rangle_{L^{2}}dz&=\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3n\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}}\frac{\omega^{2}}{\sqrt{3}}dt\\\
&\ \ \ \ \ -\int_{\mathbb{R}}\langle
A^{\ell}_{t\omega^{2}/\sqrt{3}-3(n+1)\omega/\sqrt{3}}e_{i},e_{i}\rangle_{L^{2}})\frac{\omega^{2}}{\sqrt{3}}dt.\end{split}$
(3.5)

Applying the residue theorem and using (3.4) gives (3.1). ∎

The consequences of this formula are summarized in the following theorem:

######  Theorem 5.

Consider a potential $U\in C^{\infty}(\mathbb{C}/\Gamma;\mathbb{C})$
satisfying the first two symmetries of (1.1) with finitely many non-zero
Fourier modes $c_{n}\in\mathbb{Q}(\omega/\sqrt{3})$ appearing in the
decomposition (2.1). Then for any $\ell\geq 2$, one has
$\tau_{\ell}\in\pi\mathbb{Q}(\omega/\sqrt{3})$. If $U$ also has the third
symmetry of (1.1) then the traces are real and thus
$\tau_{\ell}\in\pi\mathbb{Q}/\sqrt{3}$. In particular, for all potentials
satisfying all three symmetries in (1.1), including $U=U_{0}$ defined in
(1.2), one has

$\forall\ell\geq
2,\quad\operatorname{tr}(T_{k}^{2\ell})=\sum_{\alpha\in\mathcal{A}(U)}\alpha^{-2\ell}=\frac{\pi}{\sqrt{3}}q_{\ell},\;\;\;q_{\ell}\in\mathbb{Q},$

where $\mathcal{A}(U)$ is the set of magic angles counting multiplicity for a
potential $U.$

###### Proof.

Under the hypothesis of the corollary, the function $k\mapsto\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}$ is a rational fraction with
coefficients in $\mathbb{Q}(\omega/\sqrt{3})$. This ring is actually a field
as $\omega/\sqrt{3}$ is algebraic on $\mathbb{Q}$. Now, taking partial
fraction expansion of $k\mapsto\langle A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}$
in $\mathbb{Q}(\omega/\sqrt{3})(X)$ (the space of rational fractions with
coefficients living in $\mathbb{Q}(\omega/\sqrt{3})(X)$) gives coefficients in
$\mathbb{Q}(\omega/\sqrt{3})$. In particular, the residues of $k\mapsto\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}$ live in the field
$\mathbb{Q}(\omega/\sqrt{3})(X)$. But the uniqueness of the partial fraction
expansion now gives that these are also the residues of $k\mapsto\langle
A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}$ in $\mathbb{C}(X)$. Using the trace
formula stated in Theorem 4, this yields

$\forall\ell\geq 2,\quad\tau_{\ell}\in\pi\mathbb{Q}(\omega/\sqrt{3}).$

If we add the last symmetry of (1.1), the trace is real so that

$\forall\ell\geq
2,\quad\tau_{\ell}\in\pi\mathbb{Q}(\omega/\sqrt{3})\cap\mathbb{R}\Rightarrow\tau_{\ell}=\frac{\pi}{\sqrt{3}}q_{\ell},\;\;\;q_{\ell}\in\mathbb{Q}.$

∎

This rationality condition suffices to prove that there is an infinite number
of magic angles as long as there exists at least one magic angle.

######  Theorem 6.

Under the assumptions and with the same notation as in Theorem 5 one has the
implication

$|\mathcal{A}(U)|>0\Rightarrow|\mathcal{A}(U)|=+\infty.$

In particular, the set of magic angles for our canonical potential $U_{0}$
defined in (1.2) is infinite.
Let $N\geq 0$, for a tuple $a=(a_{n})_{\\{n;\|n\|_{\infty}\leq N\\}}$, define
$U_{a}$ to be the potential defined by (2.1). Then the above implication holds
for a generic $($in the sense of Baire$)$ set of coefficients
$a=(a_{n})_{\\{n;\|n\|_{\infty}\leq N\\}}\in\mathbb{C}^{(2N+1)^{2}}$ that
contains $(\mathbb{Q}(\omega/\sqrt{3}))^{(2N+1)^{2}}.$

###### Proof.

We start by observing that since $\pi$ is transcendental on $\mathbb{Q}$, it
is also transcendental in $\mathbb{Q}(\omega/\sqrt{3}).$ Now, assume by
contradiction, that there exist only finitely many eigenvalues
$\lambda_{i}\in\mathbb{C}$ for $i=1,..,N$ of $A_{k}^{2}$. Then we define the
$n$-th symmetric polynomial

$e_{n}(\lambda_{1},\ldots,\lambda_{N})=\sum_{1\leq
j_{1}<j_{2}<\cdots<j_{n}\leq N}\lambda_{j_{1}}\dotsm\lambda_{j_{n}}.$

Newton identities show that this polynomial can be expressed as

$e_{n}(\lambda_{1},\ldots,\lambda_{N})=(-1)^{n}\sum_{m_{1}+2m_{2}+\cdots+nm_{n}=n\atop
m_{1}\geq 0,\ldots,m_{n}\geq
0}\prod_{i=1}^{n}\frac{(-\operatorname{tr}A_{k}^{2i})^{m_{i}}}{m_{i}!i^{m_{i}}}$
(3.6)

where $e_{n}=0$ for $n>N.$ Theorem 1 shows that

$\prod_{i=1}^{n}(\operatorname{tr}A_{k}^{{2i}})^{m_{i}}\in\mathbb{Q}\left(\frac{\omega}{\sqrt{3}}\right)\pi^{m_{1}\cdots
m_{n}}.$

The power $m_{1}\cdots m_{n}$ from sequences allowed in (3.6) is maximized by
the unique choice $m=(n,0,\ldots,0)$. The Newton identities for $n>N$ then
imply that the transcendental number $\pi$ is a root of a polynomial with
coefficients in $\mathbb{Q}\left({\omega}/{\sqrt{3}}\right).$ But then all
these coefficients vanish, this is equivalent to the fact that the spectrum is
empty (because of the determinant function, see (4.4)). For our particular
choice of potential $U_{0}$, the fact that $\operatorname{tr}A^{2}_{k}=0$
contradicts [Be*22, Theorem $3$] so the set of magic angles is non empty, and
thus infinite.
Now, let $a=(a_{n})_{\\{n;\|n\|_{\infty}\leq N\\}}\in\mathbb{C}^{(2N+1)^{2}}$
and assume that $\mathcal{A}(U_{a})\neq\emptyset$. Then, we can find an open
neighbourhood of $a$, $\Omega_{a}\ni a$, such that for coefficients
$b=(b_{n})_{\\{n;\|n\|_{\infty}\leq N\\}}\in\Omega_{a}$ we have
$\mathcal{A}(U_{b})\neq\emptyset$. Take $q=(q_{n})_{\\{n;\|n\|_{\infty}\leq
N\\}}\in(\mathbb{Q}(\omega/\sqrt{3}))^{(2N+1)^{2}}\cap\Omega_{a}$ for which we
then have $|\mathcal{A}(U_{q})|=\infty.$ Continuity of eigenvalues of $T_{k}$
as the potential $U$ changes shows that the
$V_{m,a}:=\\{b\in\Omega_{a}:|\mathcal{A}(U_{b})|\geq m\\}$ is open and dense
in $\Omega_{a}$. Hence, the set coefficients for which
$0<|\mathcal{A}_{b}|<\infty$ is given by
$\bigcup_{m\in\mathbb{N}}\bigcup_{q\in(\mathbb{Q}+i\mathbb{Q})^{2N+1}}\Omega_{q}\setminus
V_{m,q}$. It is then meagre and does not contain
$(\mathbb{Q}(\omega/\sqrt{3}))^{(2N+1)^{2}}.$ ∎

##  4\. Fredholm determinants and the first magic angle

In this section, we explain how to compute the first few traces from our
formula and show the existence of a simple real magic angle, i.e. prove
Theorem 3. From now on, our choice of potential is given by $U=U_{0}$ defined
in 1.2. Here we recall some facts from [Be*22, BHZ22] needed in this paper.

###  4.1. Fourier coordinates

For our numerics, it is convenient to use rectangular coordinates $z=2i(\omega
y_{1}+\omega y_{2})$, see [Be*22, §3.3] for details. In these coordinates, we
may introduce

$\begin{gathered}\mathcal{D}_{k}:={\omega^{2}(D_{y_{1}}+k_{1})-\omega(D_{y_{2}}+k_{2})},\\\
\mathcal{V}(y):=\sqrt{3}(e^{-i(y_{1}+y_{2})}+\omega
e^{i(2y_{1}-y_{2})}+\omega^{2}e^{i(-y_{1}+2y_{2})}),\end{gathered}$ (4.1)

with periodic periodic boundary conditions (for $y\mapsto y+2\pi\mathbf{n}$,
$\mathbf{n}\in{\mathbb{Z}}^{2}$). In the following, we shall write
$\mathcal{V}_{\pm}(y):=\mathcal{V}(\pm y).$ The operator $\mathcal{A}_{k}$,
defined in (2.3), reads in the new coordinates

$\begin{gathered}\mathcal{D}_{k}^{-1}\mathcal{V}_{+}\mathcal{D}_{k}^{-1}\mathcal{V}_{-}:L^{2}(\mathbb{C}/2\pi(\mathbb{Z}+i\mathbb{Z});\mathbb{C})\to
L^{2}(\mathbb{C}/2\pi(\mathbb{Z}+i\mathbb{Z});\mathbb{C}).\end{gathered}$

On the Fourier transform side we introduce the equivalent of operators (4.1)

$\begin{gathered}\widehat{\mathcal{D}}_{k}:={\omega^{2}(D+k_{1})-\omega(D+k_{2})},\text{
with }D=\operatorname{diag}(\ell)_{\ell\in{\mathbb{Z}}}\\\
\widehat{\mathcal{V}}_{\pm}(y):=\sqrt{3}\left(J^{\pm}\otimes J^{\pm}+\omega
J^{\mp 2}\otimes J^{\pm}+\omega^{2}J^{\pm}\otimes J^{\mp
2}\right),\end{gathered}$ (4.2)

where $J$ is the right-shift $J((a_{n})_{n})=(a_{n+1})_{n}$ – see [Be*22,
(3.17)]. The spaces $L^{2}_{(p_{1},p_{2})}(\mathbb{C}/\Gamma;\mathbb{C})$,
introduced in (2.4), correspond to

$\begin{gathered}\ell_{(p_{1},p_{2})}^{2}:=\\{f\in\ell^{2}(\mathbb{Z}^{2}):\forall\,n\notin\left(3\mathbb{Z}+p_{1}\right)\times\left(3\mathbb{Z}+p_{2}\right),\>\>f_{n}=0\\}.\end{gathered}$

As in [Be*22, §3.3], we introduce auxiliary operators $J^{p,q}:=J^{p}\otimes
J^{q},\>\>p,q\in\mathbb{Z}$. For a diagonal matrix
$\Lambda=(\Lambda_{i,j})_{i,j\in\mathbb{Z}}$ acting on
$\ell^{2}(\mathbb{Z}^{2})$, we define a new diagonal matrix

$\Lambda_{p,q}:=(\Lambda_{i+p,j+q})_{i,j\in\mathbb{Z}}.$

We recall the following properties [Be*22, (3.24)]

$J^{p,q}\Lambda
J^{p^{\prime},q^{\prime}}=\Lambda_{p,q}J^{p+p^{\prime},q+q^{\prime}}=J^{p+p^{\prime},q+q^{\prime}}\Lambda_{-p^{\prime},-q^{\prime}}.$
(4.3)

Denoting the inverse of $\widehat{\mathcal{D}}_{k}^{-1}$ by

$\Lambda=\Lambda_{k}:=\widehat{\mathcal{D}}_{k}^{-1},\ \ \
\Lambda_{m,n}=\frac{1}{\omega^{2}(m+k_{1})-\omega(n+k_{2})},\ \ \
(k_{1},k_{2})\notin\mathbb{Z}^{2},$

we see that $\mathcal{A}_{k}$ reads in the new Fourier coordinates

$\frac{1}{3}\widehat{\mathcal{A}}_{k}=\Lambda\Lambda_{1,1}+\omega\Lambda\Lambda_{1,-2}+\omega^{2}\Lambda\Lambda_{-2,1}+\omega\Lambda\Lambda_{1,1}J^{3,0}+\omega^{2}\Lambda\Lambda_{1,1}J^{0,3}$
$+\omega\Lambda\Lambda_{-2,1}J^{-3,0}+\omega^{2}\Lambda\Lambda_{1,-2}J^{0,-3}+\Lambda\Lambda_{-2,1}J^{-3,3}+\Lambda\Lambda_{1,-2}J^{3,-3}$

with $\widehat{A}_{k}$ the analogous restriction to $\ell^{2}_{(1,1)}.$

###  4.2. Fredholm determinants

We start by defining the regularized Fredholm determinant

$\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})=\prod_{\lambda\in\operatorname{Spec}(A_{k})}E_{1}(\alpha^{2}\lambda)\text{
with }E_{1}(z)=(1-z)e^{z}$ (4.4)

where the product respects multiplicities. We find from (2.6) that
$\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})=0\Leftrightarrow\alpha^{-1}\in\operatorname{Spec}(T_{k})\setminus\\{0\\}.$
The symmetry of the spectrum of $\widehat{A}_{k}$,
$\operatorname{Spec}(\widehat{A}_{k})=\overline{\operatorname{Spec}(\widehat{A}_{k})},$
implies that
$\alpha\mapsto\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})$ is real-
valued on the real axis. To show existence and simplicity of magic angles, in
the representation, we therefore use the following Lemma which provides ab
initio bounds on the Fredholm determinants and its derivatives.

######  Lemma 4.1.

The determinant
${\mathbb{C}}\ni\alpha\mapsto\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})$
in (4.4) is an entire function, independent of $k\in{\mathbb{C}}$, which for
any $n,m\in\mathbb{N}_{0}$ satisfies

$\begin{split}\left\lvert\partial_{\alpha}^{m}\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})-\partial_{\alpha}^{m}\sum_{j=0}^{n}\mu_{j}\frac{(-\alpha^{2})^{j}}{j!}\right\rvert&\leq\sum_{j=n+1}^{\infty}\partial_{|\alpha|}^{m}\Bigg{(}\frac{\sqrt{e}\inf_{k\in{\mathbb{C}}}\|\widehat{A}_{k}\|_{2}|\alpha|^{2}}{\sqrt{j}}\Bigg{)}^{j}\end{split}$

with $\|A_{0}\|_{2}\leq 2$, where

${\mu}_{j}:=\operatorname{det}\begin{pmatrix}0&j-1&0&\cdots&0\\\
\sigma_{2}&0&j-2&\cdots&0\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\
\sigma_{j-1}&\sigma_{j-2}&\cdots&0&1\\\
\sigma_{j}&\sigma_{j-1}&\sigma_{j-2}&\cdots&0\end{pmatrix},\text{ with
}\sigma_{j}=\operatorname{tr}\widehat{A}_{k}^{j}.$ (4.5)

###### Proof.

The expression (4.4) is well-defined since $\widehat{A}_{k}$ is a Hilbert-
Schmidt operator and the Taylor coefficients $\mu_{j}$ are for example stated
in [Si77, (6.13)]. Indeed, since $|E_{1}|\leq e^{\frac{|z|^{2}}{2}}$ for
$z\in{\mathbb{C}}$ and
$\sum_{\lambda\in\operatorname{Spec}(\widehat{A}_{k})}|\lambda|^{2}\leq\|\widehat{A}_{k}\|^{2}_{2},$
we conclude that

$|\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})|\leq\exp\left({\frac{|\alpha|^{4}\|\widehat{A}_{k}\|^{2}_{2}}{2}}\right).$
(4.6)

Cauchy estimates for the entire function $f(z):=\det_{2}(1+z\hat{A}_{k})$ show
by using the growth bound (4.6)

$|\mu_{j}|\leq\frac{j!}{|\alpha|^{2j}}\exp\left({\frac{|\alpha|^{4}\|\widehat{A}_{k}\|^{2}_{2}}{2}}\right)$

which is optimized at $|\alpha|^{2}=\frac{\sqrt{j}}{\|\widehat{A}_{k}\|_{2}}$,
such that

$|\mu_{j}|\leq\frac{\|\widehat{A}_{k}\|^{j}_{2}e^{j/2}j!}{j^{j/2}}.$

The Taylor coefficients $\mu_{j}$ are then given by the Plemelj-Smithies
formula [Si77] stated in (4.5). Since they only depend on traces $\sigma_{j}$
which are independent of $k$, it follows that the regularized Fredholm
determinant is an entire function independent of $k.$ Hence, it suffices to
study the determinant for $k=0.$

If we write $\widehat{A}_{0}=(\widehat{A}_{0}(n))_{n\in\mathbb{Z}^{2}}$ and
let $P_{m}$ be the projection onto $(3\\{-m,-m+1,...,m\\}+1)^{2},$ then

$\begin{split}\|\widehat{A}_{0}\|_{2}&\leq\|P_{M}\widehat{A}_{0}\|_{2}+\|(\operatorname{id}-P_{M})\widehat{A}_{0}\|_{2}.\end{split}$
(4.7)

The first term constitutes the Hilbert-Schmidt norm of a finite matrix which
can be explicitly computed from the matrix elements using symbolic
calculations, indeed

$\|P_{M}\widehat{A}_{0}\|_{2}=\sqrt{\operatorname{tr}(P_{M}\widehat{A}_{0}\widehat{A}_{0}^{*}P_{M})}\leq
5\text{ for }M=760.$

To estimate the second term, we may use that the operator norm of
$\mathcal{V}_{\pm}$ satisfies $\|\mathcal{V}_{\pm}\|=3\sqrt{3}$, therefore one
has

$\begin{split}\|(\operatorname{id}-P_{M})\widehat{A}_{0}\|_{2}&\leq
9\|(\operatorname{id}-P_{M})(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(1,1)}\to\ell^{2}_{(1,1)}}\|_{4}\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(2,2)}\to\ell^{2}_{(2,2)}}\|_{4}.\end{split}$
(4.8)

We recall that by definition

$\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(1,1)}\to\ell^{2}_{(1,1)}}\|_{4}=\Bigg{(}\sum_{m\in(3{\mathbb{Z}}+1)^{2}}|\omega^{2}m_{1}-\omega
m_{2}|^{-4}\Bigg{)}^{1/4}.$

A simple change of variables shows that
$\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(1,1)}\to\ell^{2}_{(1,1)}}\|_{4}=\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(2,2)}\to\ell^{2}_{(2,2)}}\|_{4}.$
Then, a direct computation shows that in terms of

$g(m)=\frac{3((m_{1}+1)^{2}+(m_{2}+1)^{2}+(m_{1}+m_{2})^{2})}{2}-2$

we have

$\begin{split}\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(2,2)}\to\ell^{2}_{(2,2)}}\|_{4}=\frac{1}{\sqrt{3}}\Bigg{(}\sum_{m\in{\mathbb{Z}}^{2}}\frac{1}{g(m)^{2}}\Bigg{)}^{1/4}.\end{split}$

While an explicit computation shows using exact symbolic calculations

$\sum_{|m|_{\infty}\leq 6}\frac{1}{g(m)^{2}}\leq\frac{24}{7}$ (4.9)

Then, we may use for $|m|_{\infty}>6$ that $g(m)\geq|m|^{2}+5^{2}$, such that
we can estimate the remainder

$\begin{split}&\sum_{|m|_{\infty}\geq
7}\frac{1}{g(m)^{2}}\leq\int_{6}^{\infty}\frac{2\pi r}{(r^{2}+5^{2})^{2}}\
dr=\frac{\pi}{61}\\\
&\Rightarrow\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(2,2)}\to\ell^{2}_{(2,2)}}\|_{4}\leq\left(\frac{8}{21}+\frac{\pi}{549}\right)^{1/4}.\end{split}$
(4.10)

Inserting this estimate into (4.8), we find along the lines of
(LABEL:eq:remainder)

$\begin{split}\|(\operatorname{id}-P_{M})\widehat{A}_{0}\|_{2}&\leq\frac{213}{10}\|(\operatorname{id}-P_{M})(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(1,1)}\to\ell^{2}_{(1,1)}}\|_{4}\\\
&\leq\frac{213}{10}\frac{1}{\sqrt{3}}\Bigg{(}\int_{760}^{\infty}\frac{2\pi
r}{(r^{2}+759^{2})^{2}}\ dr\Bigg{)}^{1/4}<\frac{1}{2},\end{split}$ (4.11)

which shows that $\|\widehat{A}_{0}\|_{2}<\frac{11}{2}.$ ∎

Using the preceding error estimate with the explicit traces in Table 1, we
conclude the existence of a first real magic angle in the next Proposition.
The Proposition also completes the proof of Theorem 3. Indeed, (2.8) implies
together with [BHZ22, Theorem 2] the existence of a 0 gap between the two flat
bands of the Hamiltonian and the remaining bands.

######  Proposition 4.2.

There exists a simple real eigenvalue $\frac{1}{\alpha_{*}^{2}}$ to the
operator $\widehat{A}_{k}$, independent of $k\in{\mathbb{C}}$, with
$\alpha_{*}\in(0.583,0.589)$ such that
$(\frac{1}{\alpha_{*}^{2}},\infty)\subset{\mathbb{R}}\setminus\operatorname{Spec}(\widehat{A}_{k}).$

###### Proof.

To see that this is the first real magic angle, we first notice that the
operator norm of $\widehat{A}_{0}$ is bounded by

$\|\widehat{A}_{0}\|\leq(3\sqrt{3})^{2}\|(\mathscr{D}_{0}^{-1})_{\ell^{2}_{(1,1)}\to\ell^{2}_{(1,1)}}\|^{2}=9.$

This estimate shows that $\alpha\in\mathbb{R}^{+}$ with
$1/(\alpha^{2})\in\operatorname{Spec}(\widehat{A}_{0})$ satisfies
$\alpha\geq\frac{1}{3}.$ A finite number of traces as explicitly computed in
Table 1 are then relevant to prove the existence of a magic angle.

$p$ |  $\sigma_{p}\frac{\sqrt{3}}{\pi}$
---|---
1 |  $2/3$
2 |  $4$
3 |  $96/7\approx\textit{13.71}$
4 |  $40$

$p$ |  $\sigma_{p}\frac{\sqrt{3}}{\pi}$
---|---
5 |  ${28680}/{247}\approx\textit{116.14}$
6 |  $2206080/6517\approx\textit{338.51}$
7 |  $1957475168/1983163\approx\textit{987.05}$
8 |  $39948260880/13882141\approx\textit{2877.67}$

Table 1. First eight exact traces of $A^{p}_{k}$,
$\sigma_{p}=\operatorname{tr}(A_{k}^{p})$, with floating point approximation,
where ${\color[rgb]{0,0,1}\sigma_{1}}:=\lim_{n\to\infty}\sum_{|i|\leq
n}\langle A_{k}e_{i},e_{i}\rangle$ is not absolutely summable as $A_{k}$ is
not of trace-class, computed using Theorem 4 in the version stated as Theorem
7 in the appendix. One sees that the ratio of $\sigma_{p}/\sigma_{p-1}\approx
1/0.5857^{2}=2.91507,$ for $p$ large, where $0.5857$ is the first magic angle.

For $\nu\in\mathbb{R}^{+}$ we find

$r_{i}\leq\bigg{(}\frac{2\nu}{\alpha}\bigg{)}^{i}\frac{\big{(}\frac{\nu}{\sqrt{N}}\big{)}^{N-i}}{1-\frac{\nu}{\sqrt{N}}}\text{
for }r_{0}:=\sum_{k=N}^{\infty}\bigg{(}\frac{\nu}{\sqrt{k}}\bigg{)}^{k}\text{
and
}r_{1}:=\sum_{k=N}^{\infty}\frac{2k}{\alpha}\bigg{(}\frac{\nu}{\sqrt{k}}\bigg{)}^{k}.$

Evaluating the bound for $N=17$ and
$\nu=\sqrt{e}\|\widehat{A}_{0}\|_{2}\alpha^{2}$, as in the error bound stated
in Lemma 4.1, with upper bound $\|\widehat{A}_{0}\|_{2}=5.5$, we obtain for
$\alpha=\frac{3}{5}$ that $r_{0}\leq\frac{1}{50}$ and $r_{1}\leq\frac{1}{2}$.
The existence of a root follows from studying

$f(\alpha):=\sum_{k=0}^{16}\mu_{k}\frac{(-\alpha^{2})^{k}}{k!},\ \ \
\sup_{\alpha\in(1/3,\beta)}f^{\prime}(\alpha)\leq
g(\beta=0.6):=\sum_{k=2}^{20}a_{k}(\beta),$

where the summation starts at $k=2$ since $\mu_{1}=0$, with

$
a_{k}(\beta)=\begin{cases}2\mu_{k}\frac{(-1)^{k}\big{(}\tfrac{1}{3}\big{)}^{2k-1}}{(k-1)!},&\quad\text{
if }\mu_{k}(-1)^{k}<0\\\
2\mu_{k}\frac{(-1)^{k}(\frac{3}{5})^{2k-1}}{(k-1)!},&\quad\text{ if
}\mu_{k}(-1)^{k} \geq 0.\end{cases}$

One then checks (using computations involving integers only)

$\begin{split}&f(0.583)>\frac{1}{40},\qquad f(0.589)<-\frac{1}{40},\text{ and
}g(\frac{3}{5})<-\frac{7}{10}.\end{split}$

We conclude that there is $\alpha_{*}\in(0.583,0.589)$ such that
$\operatorname{det}_{2}(1-\alpha_{*}^{2}\widehat{A}_{k})=0$ and
$\partial_{\alpha}|_{\alpha=\alpha_{*}}\operatorname{det}_{2}(1-\alpha^{2}\widehat{A}_{k})<0.$
The non-existence of any other $\alpha\in(\tfrac{1}{3},\alpha_{*})$ at which
the determinant vanishes follows from the monotonicity of $f.$ ∎

Appendix: Trace formula in Fourier coordinates.

In this section we give an auxiliary version of Theorem 4 that we used for our
computer assisted computation of traces. Using the relation (4.3), the
diagonal part of $\mathcal{A}_{k}^{\ell}$ is of the form

$\begin{gathered}((\mathcal{A}_{k}^{\ell})_{ii})_{i\in{\mathbb{Z}}}=3^{\ell}\sum_{\pi\in\Theta_{\ell}}\omega^{m_{\pi}}\prod_{i=1}^{\ell}\Lambda_{\tilde{\alpha}_{i},\tilde{\beta}_{i}}\Lambda_{{\tilde{\gamma}_{i}},{\tilde{\delta}_{i}}},\\\
\pi:=\left[(\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1}),(\alpha_{2},\beta_{2}),...,(\gamma_{\ell},\delta_{\ell})\right],\end{gathered}$
(A.1)

where

$\begin{gathered}\tilde{\alpha}_{i}=\sum_{j=1}^{i-1}\alpha_{j}+\gamma_{j}\quad\tilde{\beta}_{i}=\sum_{j=1}^{i-1}\beta_{j}+\delta_{j},\quad{\tilde{\gamma}_{i}}=\alpha_{i}+\sum_{j=1}^{i-1}\alpha_{j}+\gamma_{j},\\\
\tilde{\delta}_{i}=\beta_{i}+\sum_{j=1}^{i-1}\beta_{j}+\delta_{j},\quad
m_{\pi}:=\tfrac{2}{3}\sum_{i=1}^{\ell}(\gamma_{i}+\beta_{i}).\end{gathered}$
(A.2)

In (A.1), the sum is over elements of the finite set

$\begin{split}\Theta_{\ell}&:=\Bigg{\\{}\pi=\left[(\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1}),(\alpha_{2},\beta_{2}),...,(\gamma_{\ell},\delta_{\ell})\right],\>\sum_{j=1}^{\ell}\alpha_{j}+\gamma_{j}=\sum_{j=1}^{\ell}\beta_{j}+\delta_{j}=0,\>\\\
&\quad(\alpha_{i},\beta_{i})\in\\{(1,1),(-2,1),(1,-2)\\},(\gamma_{i},\delta_{i})\in\\{(-1,-1),(2,-1),(-1,2)\\}\Bigg{\\}}.\\\
&\end{split}$ (A.3)

Using (A.1), the diagonal part of ${A}_{k}^{\ell}$, is of the form

$3^{\ell}\sum_{\pi\in\Theta_{\ell}}\omega^{m_{\pi}}\prod_{i=1}^{\ell}\Lambda^{\prime}_{\alpha_{i},\beta_{i}}\Lambda^{\prime}_{\gamma_{i},\delta_{i}},\
\ \
\pi=\left[(\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1}),(\alpha_{2},\beta_{2})...,(\gamma_{\ell},\delta_{\ell})\right],$

where $\Lambda^{\prime}$ corresponds to the matrix where we only kept the
coefficients $(n,m)$ where
$(n,m)\in\left(3\mathbb{Z}+1\right)\times\left(3\mathbb{Z}+1\right)$ i.e

$\Lambda^{\prime}_{m,n}=\frac{1}{\omega^{2}(3m+1+k_{1})-\omega(3n+1+k_{2})}.$

Theorem 4 then reduces to

######  Theorem 7.

Let ${\ell}\geq 2$ and $\Theta_{\ell}$ be as in (A.3) with coefficients
$\tilde{\alpha}_{i},..,\tilde{\delta}_{i},m_{\pi}$ as in (A.2). Then the
traces are given by

${\operatorname{tr}}\left({A}_{k}^{\ell}\right)=-\frac{2i\omega\pi}{3}\sum_{\pi\in\Theta_{\ell}}\sum_{(\eta_{i},\varepsilon_{i})\in\\{({\tilde{\alpha}_{i}},{\tilde{\beta}_{i}}),({\tilde{\gamma}_{i}},{\tilde{\delta}_{i}}),1\leqslant
i\leqslant
l\\}}{\operatorname{Res}}(f_{\pi},-\gamma_{(\eta_{i},\varepsilon_{i})})\varepsilon_{i},$

where with $\gamma_{(a,b)}=\omega^{2}a-\omega b$

$f_{\pi}(k):=3^{\ell}\omega^{m_{\pi}}\prod_{i=1}^{\ell}\frac{1}{(k+\gamma_{({\tilde{\alpha}_{i}},{\tilde{\beta}_{i}})}+\mu)(k+\gamma_{({\tilde{\gamma}_{i}},{\tilde{\delta}_{i}})}+\mu)},\quad\mu:=\omega^{2}-\omega.$

###### Proof.

This is just a re-writing of formula of Theorem 4 in these rectangular
coordinates. Indeed, the $(0,0)$-th entry of the matrix $A_{k}^{\ell}$ is, in
these notation, equal to $\sum_{\pi\in\Theta_{\ell}}f_{\pi}(k)$. Because we
work in a the Hilbert space $L^{2}(\mathbb{C}/\Gamma,\mathbb{C})$, this entry
is also equal to $\langle A_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}}$. Now, the
poles of this function are exactly described by
$\gamma_{({\tilde{\alpha}_{i}},{\tilde{\beta}_{i}})}+\mu$ and
$\gamma_{({\tilde{\gamma}_{i}},{\tilde{\delta}_{i}})}+\mu$ (this is a
consequence of formula A.1) . Note however that in these coordinates, the
poles get rescaled by $\sqrt{3}$, this is why $\mu=-i\sqrt{3}$ replaces $-i$.
On the level of residues, this explains why a $\sqrt{3}$ does not appear in
this formula. Finally, we remark that in this decomposition,
$\gamma_{({\tilde{\alpha}_{i}},{\tilde{\beta}_{i}})}+\mu\in\sqrt{3}\big{(}3\Gamma^{*}-i\big{)}$
and
$\gamma_{({\tilde{\gamma}_{i}},{\tilde{\delta}_{i}})}+\mu\in\sqrt{3}\big{(}3\Gamma^{*}-2i\big{)}$
thus corresponding to the splitting appearing in the formula of stated in
Theorem 4. ∎

Acknowledgements. We would like to thank the two anonymous referees for their
careful reading of the paper and many helpful comments. TH and MZ gratefully
acknowledge partial support by the National Science Foundation under the grant
DMS-1901462 and by the Simons Foundation Targeted Grant Award No. 896630.

Conflict of interests. The authors have no competing interests to declare that
are relevant to the content of this article.

## References

  * [Be*21] S. Becker, M. Embree, J. Wittsten and M. Zworski, Spectral characterization of magic angles in twisted bilayer graphene, Phys. Rev. B 103, 165113, (2021).
  * [Be*22] S. Becker, M. Embree, J. Wittsten and M. Zworski, Mathematics of magic angles in a model of twisted bilayer graphene, Probab. Math. Phys. 3 (2022), 69–103.
  * [BHZ22] S. Becker, T. Humbert and M. Zworski, Fine structure of flat bands in a chiral model of magic angles, arXiv:2208.01628.
  * [BiMa11] R. Bistritzer and A. MacDonald, Moiré bands in twisted double-layer graphene. PNAS, 108, 12233–12237, (2011).
  * [CGG22] E. Cancès, L. Garrigue, D. Gontier, _A simple derivation of moiré-scale continuous models for twisted bilayer graphene._ arXiv:2206.05685.
  * [DuNo80] B.A. Dubrovin and S.P. Novikov, Ground states in a periodic field. Magnetic Bloch functions and vector bundles. Soviet Math. Dokl. 22, 1, 240–244, (1980).
  * [Ka80] T. Kato, _Perturbation Theory for Linear Operators_ , Corrected second edition, Springer, Berlin, (1980).
  * [SGG12] P. San–Jose, J. González, and F. Guinea, _Non-Abelian gauge potentials in graphene bilayers,_ Phys. Rev. Lett. 108, 216802 (2012).
  * [Si77] B. Simon, _Notes on infinite determinants of Hilbert space operators_ , Adv. in Math. 24 (1977), 244-273.
  * [TKV19] G. Tarnopolsky, A.J. Kruchkov and A. Vishwanath, _Origin of magic angles in twisted bilayer graphene_ , Phys. Rev. Lett. 122, 106405, (2019).
  * [Wa*22] A. B. Watson, T. Kong, A. H. MacDonald, and M. Luskin _Bistritzer-MacDonald dynamics in twisted bilayer graphene_ , arXiv:2207.13767.
  * [WaLu21] A. Watson and M. Luskin, Existence of the first magic angle for the chiral model of bilayer graphene, J. Math. Phys. 62, 091502 (2021).


# Hybrid Classical-Quantum Simulation of MaxCut using QAOA-in-QAOA ††thanks: *
Equal Contribution

Aniello Esposito*  EMEA Research Lab
Hewlett Packard Enterprise
Basel, Switzerland
<EMAIL_ADDRESS>Tamuz Danzig  Classiq Technologies
Tel Aviv, Israel
<EMAIL_ADDRESS>
###### Abstract

The Quantum approximate optimization algorithm (QAOA) is a leading hybrid
classical-quantum algorithm for solving complex combinatorial optimization
problems. QAOA-in-QAOA (QAOA2) uses a divide-and-conquer heuristic to solve
large-scale Maximum Cut (MaxCut) problems, where many sub-graph problems can
be solved in parallel. In this work, an implementation of the QAOA2 method for
the scalable solution of the MaxCut problem is presented, based on the Classiq
platform. The framework is executed on an HPE-Cray EX supercomputer by means
of the Message Passing Interface (MPI) and the SLURM workload manager. The
limits of the Goemans-Williamson (GW) algorithm as a purely classical
alternative to QAOA are investigated to understand if QAOA2 could benefit from
solving certain sub-graphs classically. Results from large-scale simulations
of up to 33 qubits are presented, showing the advantage of QAOA in certain
cases and the efficiency of the implementation, as well as the adequacy of the
workflow in the preparation of real quantum devices. For the considered
graphs, the best choice for the sub-graphs does not significantly improve
results and is still outperformed by GW.

###### Index Terms:

QAOA, MaxCut, supercomputing, hybrid classical-quantum

##  I Introduction

Quantum computation (QC) is an active area of research that has attracted key
players in academia and industry, primarily due to its potential to solve
specific problems much more efficiently that their classical counterparts [1,
2]. However, current Noisy Intermediate-Scale Quantum (NISQ) devices feature a
modest number of qubits and useful compute time is limited due to the
decoherence [3, 4]. This has driven the development of methods to partition
quantum circuits whose components are executed on smaller quantum devices with
assistance of classical processing [5, 6, 7, 8, 9, 10, 11, 12, 13]. These
methods can be roughly subdivided in circuit knitting and entanglement forging
and are naturally implemented as hybrid classical-quantum workflows. While
quantum devices are expected to perform well for certain problems involving a
small amount of data but high complexity, current classical supercomputers
continue to serve researchers in areas with data-intensive workloads relying
on software tools and infrastructure developed over decades. Therefore, in
conjunction with a number of possibly tightly integrated NISQ devices,
classical supercomputers would be valuable candidates for the execution of
hybrid workflows while sustaining the classical portion of the computation and
communication of potentially large data volumes. However, while considerable
resources have been invested in quantum circuit simulators to accompany the
development of quantum hardware, the investigation of high performance
computing (HPC) systems featuring quantum hardware has though driven
programming environments[14, 15, 16, 17] but other operational aspects such as
workflows[18, 19, 20, 21, 22, 23, 24, 25] have not received a comparable
attention.

This work focuses on the hybrid workflow, which is investigated by means of
the SLURM[26] workload manager on a HPE-Cray EX supercomputer which executes
quantum circuits on simulated devices. The case study consists of a MaxCut
problem, which has gained popularity in the QC community due to its natural
mapping to qubit measurement results. The solution can be approximated by the
QAOA [27] on quantum devices for general quantum circuit execution[28] or
conversely formulated as quadratic unconstraing binary optimization (QUBO)
problem and solved with quantum annealers[29]. The most popular classical
alternative for approximate solutions is given by the GW algorithm.
Partitioning schemes for the MaxCut problem and QAOA have been extensively
investigated, either targeting the generated quantum circuit[30, 31] or the
original classical graph[32], i.e. QAOA2, using a divide-and-conquer
heuristic. The latter approach is considered in this work and implemented with
the support of the Classiq platform [33], being able to synthesize more
optimized quantum circuits compared to a manual construction. Circuits are
executed with the aer simulator [34] using distributed memory parallelism over
MPI. This workflow naturally maps to the Multiple-Program-Multiple-Data (MPMD)
or heterogeneous jobs paradigm of SLURM. While QAOA2 considers only QAOA for
the sub-graphs, the SLURM features allows the allocation of a mixture of
quantum and classical resources giving a choice between QAOA or a classical
solution of the sub-graph problem, for instance, with GW to maximize the
overall MaxCut. To that end, both methods are applied to a set of graphs to
assess the suitability of such a run-time decision mechanism. The present
implementation also provides a testbed to further refine such a mechanism and
possibly include machine learning (ML) methods similar to what has been
considered in [35].

The paper starts with a consideration of related work, followed by details and
remarks about the methods. Comparisons between the methods for moderate graph
sizes are presented before the final method is applied to larger graphs.
Finally, concluding remarks and an outlook are given.

##  II Related Work

Hybrid classical-quantum workflows have been investigated[18, 19, 20, 21, 22,
23, 24, 25] considering cloud environments and more tightly integrated
infrastructures. A number of expectations for components such as workload
managers or circuit depths have emerged with related challenges. In this work,
the ingredients of the workflow are kept simple, using a state-of-the-art
workload manager and circuit simulator in order to focus on the suitability of
a supercomputer by means of a well-known example application, where the
circuit depth is optimized by means of the Classiq platform. Partitioning of
the MaxCut problem with circuit knitting of the resulting QAOA formulation has
also been investigated[30, 31]. However, while this approach allows, in
principle, the use of multiple quantum devices simultaneously, it discards the
possibility to employ classical methods for a portion of the problem where it
could be more advantageous. In this work, the QAOA2 method is considered as a
more versatile example of hybrid application, where the GW method [36] is used
for the classical solution of a sub-graph if more advantageous.

ML methods have been considered by Moussa et al. [35] who created a classifier
to predict the suitability of a graph instance for eitherQAOA or GW. They
achieved 96% prediction accuracy, though not at the qubit counts considered in
this work. Their approach could be applied to a QAOA2 context as a next step.
The present work provides a testbed to train and test such selection
mechanisms. Furthermore, with a large dataset of QAOA results, a neural
network can be trained to predict initial parameters for subsequent QAOA
simulations or computations on real quantum hardware. This could improve the
number of iterations in the hybrid scheme of QAOA while preserving the
accuracy [37].

Finally, it is worth mentioning that exact classical methods to solve the
MaxCut have been extensively investigated[38] but the node counts in typical
applications are still limited compared to GW. Also probabilistic methods from
statistical physics such as simulated annealing [39] have been considered for
MaxCut.

##  III Methods

###  III-A The MaxCut problem

A graph $G=(V,E)$ is considered, with a set of nodes $V={1,2,...,N}$, and a
set of edges $(i,j)\in E$ between the nodes $i$ and $j$ with weights
$w_{ij}=w_{ji}$. The goal is to divide the nodes into two groups in such a way
that maximizes the sum of weights of edges that touch nodes from the two
different groups (cut) rather than the same group. The problem represents a
NP-hard combinatorial problem.

The problem Hamiltonian for a given graph is:

$H_{C}=\frac{1}{2}\sum_{(i,j)\in E}w_{ij}\left(1-Z_{i}Z_{j}\right)\quad,$ (1)

where each qubit represents a node, and $Z_{i}$ is the Pauli operator that
acts on qubit $i$. If the two qubits of the edge $(i,j)$ are measured
differently, it adds to the overall sum. If they are measured similarly, there
is no contribution to the sum.

###  III-B Quantum Approximate Optimization Algorithm

The method belongs to the family of Variational Quantum Algorithm (VQA) and
makes use of a specific parameterized quantum circuit structure (ansatz) [27]
composed of $p$-layers

$\ket{\psi_{p}(\vec{\beta},\vec{\gamma})}=\left[\prod_{l=1}^{p}e^{-i\beta_{l}H_{M}}e^{-i\gamma_{l}H_{C}}\right]\ket{+}^{\otimes
n},$ (2)

where $H_{C}$ is called the problem Hamiltonian, which is defined by the
specific problem to be solved, $H_{M}$ is the mixer Hamiltonian that takes the
same structure for any given problem, $n$ is the number of qubits in the
quantum circuit, and $\vec{\gamma}$ and $\vec{\beta}$ are the variational real
parameters that correspond to $H_{C}$ and $H_{M}$, respectively. The objective
function is given by

$F_{p}(\vec{\beta},\vec{\gamma})=\bra{\psi_{p}(\vec{\beta},\vec{\gamma})}H_{C}\ket{\psi_{p}(\vec{\beta},\vec{\gamma})}\quad,$
(3)

which needs to be maximized (or minimized) by changing $\vec{\gamma}$ and
$\vec{\beta}$ values in each iteration by a classical optimizer. The number of
shots for the circuit simulation used in this work is $4096$. The rigid
structure of the Ansatz can be considered as a time-discretization of
adiabatic quantum computation [27, 40, 41], and as $p\rightarrow\infty$, QAOA
obtain the optimal solution [27]. The MaxCut problem has become the de-facto
standard benchmark for QAOA due to its simple formulation and suitability for
binary encoding [42, 43]. Several works claim that QAOA has the potential to
achieve quantum advantage. Farhi and Harrow argue that QAOA can reach quantum
advantage based on complexity theoretic assumptions [44]. Guerreschi claims
that "QAOA for MaxCut requires hundreds of qubits for quantum speed-up" [45].
The number of qubits indeed grows in that direction. Crooks shows that when
the number of layers is more than $8$, QAOA achieves better results than GW
[40]. Shaydulin et al. found that QAOA gives the best empirical scaling of any
known classical algorithm for the low auto-correlation binary sequence problem
[46]. Bravyi et al. proposed a non-local variation of QAOA, called recursive-
QAOA (RQAOA), which numerically outperforms standard QAOA [47]. This method
can also be leveraged using QAOA2 to get a good global solution for very large
problems in a reasonable time. The VQA family of algorithms are especially
suitable for the NISQ era because they use shallow quantum circuits with the
assistance of classical computers [48, 49]. In this work, once the QAOA
circuit is executed, the bit string corresponding to the highest amplitude in
the resulting state vector is chosen as a solution for sake of simplicity. A
more appropriate approach would have been to consider a number of highest
amplitudes and chose the bit string yielding the highest cut among them.

###  III-C QAOA-in-QAOA

The method proposed by Zhou et al. [32] revisits the MaxCut problem via the
divide-and-conquer heuristic by seeking the solutions of sub-graphs in
parallel and then merging these solutions to obtain the global solution. It
consists of the following steps.

  1. 1.

Specify the number of qubits $n$ in our quantum computers or simulators, the
number of layers in the quantum circuit, and the number of iterations.

  2. 2.

For the dividing procedure, the graph $G$ is partitioned into sub-graphs in
which the number of nodes does not exceed a specified number of qubits $n$
from the previous step. The greedy modularity method from the NetworkX Library
is used, which maximizes the modularity of the graph. If a sub-graph has more
nodes than $n$, the sub-graph is divided into fewer sub-graphs, recursively.

  3. 3.

All sub-graphs are solved with QAOA in parallel over different (simulated)
quantum devices.

  4. 4.

In the merging process, the overall solutions of all sub-graphs are
considered. In order to obtain a better total solution, a new graph is
generated as follows:

    1. (a)

Each sub-graph is represented by a node.

    2. (b)

Each edge that connects between two sub-graphs and is also part of the cut, we
multiply its weight with $-1$. Each edge that is not part of the cut, we don’t
change.

    3. (c)

Take the sum on all edges between each two sub-graphs to add an edge with a
single weight to the new graph.

  5. 5.

The resulting graph is solved with QAOA. If a node in the new graph is $-1$,
all the nodes in the sub-graph represented by this node are flipped. If the
new graph has more nodes than in the quantum computer, the process is repeated
recursively.

A total of $\sim\sum_{k=1}^{a}N/n^{k}=N(n^{a}-1)/(n^{a}(n-1))$ graphs have to
be considered for $a$ levels. The condition $\lceil N/n^{a}\rceil\approx n$,
i.e. that the lowest level graph is computed with $n$ qubits, yields
$a\approx\lceil log_{n}N\rceil-1$. If sufficient quantum devices are available
to be used in parallel, the complexity reduces to $\sim O(\log_{n}N)$. The
original work shows experimentally that for different types of graphs, QAOA2
achieves competitive or even better results compared to GW for 2000 nodes. On
the theoretical side, the algorithm is proven to be better than a random
solution MaxCut, and a lower bound to the approximation ratio of QAOA2 is
given, which depends on the partition method and the local solver.

###  III-D Goemans-Williamson Algorithm

It is the best-known classical approximation algorithm for the MaxCut problem,
with an approximation ratio of $0.878$. The computationally most intensive
part is the optimization of a semi-definite rogramming (SDP) problem. Once the
SDP is solved, a slicing to determine the node values is applied $30$ times,
and the average value of the cut is taken. This provides a more suitable way
to compare to QAOA, which is not repeated and averaged in the present work.
For dense problems, the time to solution and memory consumption of GW grow
like $O(N^{6.5})$ and $O(N^{4})$[50], respectively, where $N$ is the number of
nodes in the graph. The cvxpy and networkx Python packages have been used for
the solution of the SDP problem and the graphs, respectively. The cvxpy
disposes of a partial OpenMP shared memory parallelization via the splitting
conic solver (SCS) in addition to the underlying threaded BLAS/LAPACK
libraries but the scalability is rather limited and therefore only one OpenMP
thread is used.

###  III-E Circuit Synthesis

The Classiq platform [33] allows a user to build a high-level functional
quantum model. The synthesis engine takes this model along with optimization
preferences and global constraints to produce an optimized quantum circuit
while considering many different implementations. The synthesis engine can
optimize over circuit depth, number of qubits, two-qubit gates, or create an
optimized implementation into a specific hardware architecture, etc. In the
present case, the description of a high-level combinatorial optimization
problem is converted into an optimized gate-level quantum circuit.

###  III-F Workflow Execution

The MPMD paradigm in SLURM can be used to compute the solutions of the MaxCut
problems of the sub-graphs generated by QAOA2, either classically by the GW
method or quantum mechanically (simulated) by the QAOA algorithm. This
distinction is only made for the first partitioning of the original graph in
QAOA2 for the sake of simplicity in this preliminary investigation. A
schematic view of the procedure is given in Fig. 2, where a coordinator could
inspect the sub-graphs and calculate the most appropriate resource allocation
in advance. Furthermore, in order to save idle time caused by a possibly
scarce usage of a quantum device, the heterogeneous jobs paradigm of SLURM
could be employed, as shown in Fig. 1.

Figure 1: Heterogeneous jobs for the reduction of idle time of a quantum
device. Before the first heterogeneous job finishes, a second one can already
start using the quantum device. In this scenario, the quantum device is always
accessed exclusively by a single user.

Figure 2: Distribution scheme for the solution of the sub-graphs generated by
QAOA2, either classically or quantum mechanically (simulated). A coordinator
executed on a dedicated MPI rank handles the partitioning and collection of
results. Usually, the consumption of classical and quantum resources does not
start at the same time. However this can be achieved by splitting,
checkpointing, and restarting the classical part appropriately.

An HPE-Cray EX supercomputer featuring compute nodes with two AMD EPYC 7763
64-Core processors and 512 GB of memory is used for the experiments.
Distributed memory parallelization is achieved through the mpi4py[51] package.

##  IV Results

To identify an advantage of QAOA over GW for the solution of the sub-graphs in
QAOA2, a grid search over a number of circuit layers $p\in\\{3,4,\ldots,8\\}$
from Eq. 2 and values of the initial change to the variables for the COBYLA
optimizer $\rm rhobeg\in\\{0.1,0.2,0.3,0.4,0.5\\}$ has been considered, where
the employed Erdős–Rényi graphs are generated with the networkx Python
package. The number of iterations is chosen to be linearly dependent on $p$
and ranges from $30$ to $100$ steps. The grid search is applied to every graph
in a set with node counts ranging from $15$ to $25$ and edge probabilities
from $0.1$ to $0.5$. A graph instance with uniform edges and one with edge
weights randomly chosen in $[0,1]$ is created for every node count and edge
probability. The proportions of cases in which QAOA yields strictly larger
MaxCut values than GW are recorded, as well as the limiting case where QAOA
reaches $[95,100)\%$ of the GW value. When QAOA is strictly better than GW, a
score for the corresponding data point in the search grid is increased by $1$
and finally normalized. All results are summarized in Fig. 3.

(a) Proportions of cases in which QAOA is strictly better than GW for
unweighted (left) and weighted (right) graphs.

(b) Proportions of cases in which the QAOA value is $[95,100)\%$ of GW for
unweighted (left) and weighted (right) graphs.

(c) Proportions of cases in the grid search for which QAOA is strictly better
than GW for unweighted (left) and weighted (right) graphs.

Figure 3: Results from a grid search over different number of circuit layers
and optimization parameters. Weighted and unweighted graphs with different
node counts $15$ to $25$ and edge probabilities $0.1$ to $0.5$ have been
considered, i.e. 30 grid points.

This creates a simple, yet instructive, knowledge base about which type of
parameterization of QAOA is more suitable for a type of graph or whether a
classical solution is better overall. This knowledge can in turn be used to
optimally process a set of sub-graphs resulting from a step in QAOA2. The grid
search reveals that the QAOA has a partial advantage for graphs with small
edge connection probabilities and that the most successful parameter
combination is $(\rm rhobeg=0.5,p=6)$. For higher $p$-layers, QAOA is expected
to reach better result using more iterations or better initial parameters. In
a naive approach, all sub-graphs with small edge probabilities would be
simulated with QAOA using the most successful parameterization or classically
otherwise. A more advanced approach could use ML to train a predictive model
and test the performance using the workflow described in Sec. III-F. The same
analysis is repeated for node counts $30$ to $33$ and edge probabilities $0.1$
and $0.2$ and summarized in Table I.

Table I: Analogous data to Fig. 3 (a) and (b) for node counts $30$ to $33$ and edge probabilities $0.1$ and $0.2$. The top table shows the proportions of cases in which QAOA2 is strictly better than GW and the bottom table shows the proportions of cases in which the QAOA2 value is $[95,100)\%$ of GW. Number of | Weighting | Edge Probability
---|---|---
nodes |  | 0.1 | 0.2
30 | yes | 0.1 | 0.1
| no | 0.167 | 0
31 | yes | 0.267 | 0.033
| no | 0 | 0.067
32 | yes | 0.1 | 0.033
| no | 0.1 | 0
33 | yes | 0.033 | 0.033
| no | 0.167 | 0.033
30 | yes | 0.133 | 0.2
| no | 0.33 | 0.1
31 | yes | 0.1 | 0.1
| no | 0.2 | 0.033
32 | yes | 0.167 | 0.067
| no | 0.167 | 0.133
33 | yes | 0.067 | 0.167
| no | 0.2 | 0.067

The occurrences of QAOA being strictly better than GW are less frequent and no
distinct point in the parameter grid can be identified as for the lower node
counts. Though, a high rhobeg or a high number of layers seem more successful.
The simulation of QAOA for $33$ qubits takes approximately 10 minutes on 512
compute nodes for $p=8$. The QAOA2 is finally applied to a set of larger
unweighted graphs with node counts $\\{500,1000,1500,2000,2500\\}$ and edge
probability $0.1$. The sub-graphs from the first partition are analyzed with
the same parameter grid search from before, and the QAOA solution with the
highest MaxCut value is stored along with a GW solution for that sub-graph,
where the best result among the two is recorded. In case of further iterations
in the QAOA2 method, the classical solution is chosen. The results are
summarized in Fig. 4 together with the GW solution of the original graph and a
random partition solution with approximation.maxcut from the NetworkX
algorithms suite.

Figure 4: MaxCut value for unweighted graphs with different node counts and
edge probability $0.1$. The sub-graphs from the first partition of the QAOA2
method are solved either all with QAOA (blue), or all classically with GW
(green), or with the best among the two (yellow). The result from GW applied
to the original graph is reported as well (black) together with a random
partition solution (red) from the NetworkX algorithms suite. Data is relative
to the QAOA solution.

The GW method applied to the full graph is superior compared to the other
schemes up to $2000$ nodes, where abnormal terminations are encountered, and
diminishes steadily compared to QAOA2 for larger node counts. The abnormal
termination of GW beyond $2000$ nodes seems to be related to the triplet
representation in the Eigen library used by cvxpy and not yet to memory
consumption. Choosing the best solution among QAOA and GW for the sub-graphs
yields slightly better results compared to solving all sub-graphs either with
QAOA or GW in QAOA2. Finally, all methods are better than a random cut. The
overhead incurred by the coordination of the various sub-graph solutions is
minimal and overall an almost ideal scaling is achieved.

##  V Conclusion

The QAOA method has been compared to the GW algorithm for the MaxCut problem
of graphs with node counts up to $33$ and equivalent number of simulated
qubits. For certain parameters, the QAOA can outperform its classical
counterpart. This knowledge base motivates the selection of a most appropriate
method for the solution of sub-graphs in the QAOA2 method. Experimental
results show that the best choice for a sub-graph does not yield significantly
better results compared to a purely classical or (simulated) quantum
mechanical solution and is still substantially worse than the GW method for
the entire graph. However, the GW advantage diminishes for larger node counts,
which leads to the assumption that QAOA2 could have an effective advantage at
some point. This motivates the investigation of other graph types and
partitions including more statistics. Also, considering a larger number of
amplitudes in the resulting state vectors is expected to significantly improve
the QAOA results. The presented simulation infrastructure based on SLURM,
which constitutes another research interest, can be further exploited for this
analysis. In a future work, the extension to heterogeneous SLURM jobs as well
as the use of ML methods for the selection of the best method will be
investigated.

## Acknowledgment

We want to thank Nadav Yoran and Adam Goldfeld from Classiq Technologies as
well as Alfio Lazzaro, Jess Jones, and David Brayford from HPE for fruitful
discussions.

## References

  * [1] M. A. Nielsen and I. L. Chuang, _Quantum computation and quantum information_. Cambridge university press, 2010.
  * [2] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell _et al._ , “Quantum supremacy using a programmable superconducting processor,” _Nature_ , vol. 574, no. 7779, pp. 505–510, 2019.
  * [3] K. Sharma, S. Khatri, M. Cerezo, and P. J. Coles, “Noise resilience of variational quantum compiling,” _New Journal of Physics_ , vol. 22, no. 4, p. 043006, 2020.
  * [4] S. Stein, N. Wiebe, Y. Ding, P. Bo, K. Kowalski, N. Baker, J. Ang, and A. Li, “Eqc: ensembled quantum computing for variational quantum algorithms,” in _Proceedings of the 49th Annual International Symposium on Computer Architecture_ , 2022, pp. 59–71.
  * [5] P. de Schoulepnikoff, O. Kiss, S. Vallecorsa, G. Carleo, and M. Grossi, “Hybrid ground-state quantum algorithms based on neural schr$\backslash$" odinger forging,” _arXiv preprint arXiv:2307.02633_ , 2023.
  * [6] P. Huembeli, G. Carleo, and A. Mezzacapo, “Entanglement forging with generative neural network models,” _arXiv preprint arXiv:2205.00933_ , 2022.
  * [7] G. Gentinetta, F. Metz, and G. Carleo, “Overhead-constrained circuit knitting for variational quantum dynamics,” _arXiv preprint arXiv:2309.07857_ , 2023.
  * [8] L. Brenner, C. Piveteau, and D. Sutter, “Optimal wire cutting with classical communication,” _arXiv preprint arXiv:2302.03366_ , 2023.
  * [9] C. Piveteau and D. Sutter, “Circuit knitting with classical communication,” _IEEE Transactions on Information Theory_ , 2023.
  * [10] K. Mitarai and K. Fujii, “Overhead for simulating a non-local channel with local channels by quasiprobability sampling,” _Quantum_ , vol. 5, p. 388, 2021.
  * [11] ——, “Constructing a virtual two-qubit gate by sampling single-qubit operations,” _New Journal of Physics_ , vol. 23, no. 2, p. 023021, 2021.
  * [12] T. Peng, A. W. Harrow, M. Ozols, and X. Wu, “Simulating large quantum circuits on a small quantum computer,” _Physical review letters_ , vol. 125, no. 15, p. 150504, 2020.
  * [13] S. Bravyi, G. Smith, and J. A. Smolin, “Trading classical and quantum computational resources,” _Physical Review X_ , vol. 6, no. 2, p. 021043, 2016.
  * [14] J.-S. Kim, A. McCaskey, B. Heim, M. Modani, S. Stanwyck, and T. Costa, “Cuda quantum: The platform for integrated quantum-classical computing,” in _2023 60th ACM/IEEE Design Automation Conference (DAC)_. IEEE, 2023, pp. 1–4.
  * [15] “Qaptiva (atos),” https://atos.net/en/solutions/high-performance-computing-hpc/quantum-computing-qaptiva.
  * [16] X.-C. Wu, P. Khalate, A. Schmitz, S. Premaratne, K. Rasch, S. Daraeizadeh, R. Kotlyar, S. Ren, J. Paykin, F. Rose _et al._ , “Intel quantum sdk version 1.0: Extended c++ compiler, runtime and quantum hardware simulators for hybrid quantum-classical applications,” _Bulletin of the American Physical Society_ , 2023.
  * [17] J. K. Lee, O. T. Brown, M. Bull, M. Ruefenacht, J. Doerfert, M. Klemm, and M. Schulz, “Quantum task offloading with the openmp api,” _arXiv preprint arXiv:2311.03210_ , 2023.
  * [18] A. Esposito, J. R. Jones, S. Cabaniols, and D. Brayford, “A hybrid classical-quantum hpc workload,” in _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , vol. 2. IEEE, 2023, pp. 117–121.
  * [19] A. Elsharkawy, X.-T. M. To, P. Seitz, Y. Chen, Y. Stade, M. Geiger, Q. Huang, X. Guo, M. A. Ansari, M. Ruefenacht _et al._ , “Challenges in hpcqc integration,” in _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , vol. 2. IEEE, 2023, pp. 405–406.
  * [20] M. N. Farooqi and M. Ruefenacht, “Exploring hybrid classical-quantum compute systems through simulation,” in _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , vol. 2. IEEE, 2023, pp. 127–133.
  * [21] B. Mete, M. Schulz, and M. Ruefenacht, “Predicting the optimizability for workflow decisions,” in _2022 IEEE/ACM Third International Workshop on Quantum Computing Software (QCS)_. IEEE, 2022, pp. 68–74.
  * [22] M. Schulz, L. Schulz, M. Ruefenacht, and R. Wille, “Towards the munich quantum software stack: Enabling efficient access and tool support for quantum computers,” in _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , vol. 2. IEEE, 2023, pp. 399–400.
  * [23] M. Ruefenacht, B. G. Taketani, P. Lähteenmäki, V. Bergholm, D. Kranzlmüller, L. Schulz, and M. Schulz, “Bringing quantum acceleration to supercomputers,” 2022.
  * [24] M. Schulz, M. Ruefenacht, D. Kranzlmüller, and L. B. Schulz, “Accelerating hpc with quantum computing: It is a software challenge too,” _Computing in Science & Engineering_, vol. 24, no. 4, pp. 60–64, 2022.
  * [25] R. F. da Silva, R. M. Badia, V. Bala, D. Bard, P.-T. Bremer, I. Buckley, S. Caino-Lores, K. Chard, C. Goble, S. Jha _et al._ , “Workflows community summit 2022: A roadmap revolution,” _arXiv preprint arXiv:2304.00019_ , 2023.
  * [26] A. B. Yoo, M. A. Jette, and M. Grondona, “Slurm: Simple linux utility for resource management,” in _Workshop on job scheduling strategies for parallel processing_. Springer, 2003, pp. 44–60.
  * [27] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” _arXiv preprint arXiv:1411.4028_ , 2014.
  * [28] R. Shaydulin and M. Pistoia, “Qaoa with $\mathbb{N}\cdot p\geq 200$,” _arXiv preprint arXiv:2303.02064_ , 2023.
  * [29] W. van der Schoot, D. Leermakers, R. Wezeman, N. Neumann, and F. Phillipson, “Evaluating the q-score of quantum annealers,” in _2022 IEEE International Conference on Quantum Software (QSW)_. IEEE, 2022, pp. 9–16.
  * [30] A. Lowe, M. Medvidović, A. Hayes, L. J. O’Riordan, T. R. Bromley, J. M. Arrazola, and N. Killoran, “Fast quantum circuit cutting with randomized measurements,” _Quantum_ , vol. 7, p. 934, 2023.
  * [31] M. Bechtold, J. Barzen, F. Leymann, A. Mandl, J. Obst, F. Truger, and B. Weder, “Investigating the effect of circuit cutting in qaoa for the maxcut problem on nisq devices,” _arXiv preprint arXiv:2302.01792_ , 2023.
  * [32] Z. Zhou, Y. Du, X. Tian, and D. Tao, “Qaoa-in-qaoa: solving large-scale maxcut problems on small quantum machines.”
  * [33] Classiq Technologies Ltd., “Classiq platform.” [Online]. Available: https://platform.classiq.io/
  * [34] J. Doi and H. Horii, “Cache blocking technique to large scale quantum computing simulation on supercomputers,” in _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_. IEEE, 2020, pp. 212–222.
  * [35] C. Moussa, H. Calandra, and V. Dunjko, “To quantum or not to quantum: towards algorithm selection in near-term quantum optimization.”
  * [36] M. X. Goemans and D. P. Williamson, “Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming,” _Journal of the ACM (JACM)_ , vol. 42, no. 6, pp. 1115–1145, 1995.
  * [37] O. Amosy, T. Danzig, E. Porat, G. Chechik, and A. Makmal, “Iterative-free quantum approximate optimization algorithm using neural networks,” 8 2022. [Online]. Available: https://arxiv.org/abs/2208.09888v1
  * [38] D. Rehfeldt, T. Koch, and Y. Shinano, “Faster exact solution of sparse maxcut and qubo problems,” _Mathematical Programming Computation_ , pp. 1–26, 2023.
  * [39] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi, “Optimization by simulated annealing,” _science_ , vol. 220, no. 4598, pp. 671–680, 1983.
  * [40] G. E. Crooks, “Performance of the quantum approximate optimization algorithm on the maximum cut problem,” _arXiv preprint arXiv:1811.08419_ , 2018.
  * [41] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, “Quantum computation by adiabatic evolution,” _arXiv preprint quant-ph/0001106_ , 2000.
  * [42] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, “Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices,” _Physical Review X_ , vol. 10, no. 2, p. 021067, 2020.
  * [43] M. Willsch, D. Willsch, F. Jin, H. De Raedt, and K. Michielsen, “Benchmarking the quantum approximate optimization algorithm,” _Quantum Information Processing_ , vol. 19, pp. 1–24, 2020.
  * [44] E. Farhi and A. W. Harrow, “Quantum supremacy through the quantum approximate optimization algorithm,” _arXiv preprint arXiv:1602.07674_ , 2016.
  * [45] G. G. Guerreschi and A. Y. Matsuura, “Qaoa for max-cut requires hundreds of qubits for quantum speed-up,” _Scientific reports_ , vol. 9, no. 1, pp. 1–7, 2019.
  * [46] R. Shaydulin, C. Li, S. Chakrabarti, M. DeCross, D. Herman, N. Kumar, J. Larson, D. Lykov, P. Minssen, Y. Sun, Y. Alexeev, J. M. Dreiling, J. P. Gaebler, T. M. Gatterman, J. A. Gerber, K. Gilmore, D. Gresh, N. Hewitt, C. V. Horst, S. Hu, J. Johansen, M. Matheny, T. Mengle, M. Mills, S. A. Moses, B. Neyenhuis, P. Siegfried, R. Yalovetzky, and M. Pistoia, “Evidence of scaling advantage for the quantum approximate optimization algorithm on a classically intractable problem,” 8 2023. [Online]. Available: https://arxiv.org/abs/2308.02342v1
  * [47] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, “Obstacles to variational quantum optimization from symmetry protection,” _Physical Review Letters_ , vol. 125, no. 26, p. 260505, 2020.
  * [48] J. Preskill, “Quantum computing in the nisq era and beyond,” _Quantum_ , vol. 2, p. 79, 2018.
  * [49] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke _et al._ , “Noisy intermediate-scale quantum algorithms,” _Reviews of Modern Physics_ , vol. 94, no. 1, p. 015004, 2022.
  * [50] R. Y. Zhang and J. Lavaei, “Sparse semidefinite programs with near-linear time complexity,” _Proceedings of the IEEE Conference on Decision and Control_ , vol. 2018-December, pp. 1624–1631, 7 2018.
  * [51] L. Dalcin and Y.-L. L. Fang, “mpi4py: Status update after 12 years of development,” _Computing in Science & Engineering_, vol. 23, no. 4, pp. 47–54, 2021.