MnLargeSymbols’164 MnLargeSymbols’171

# Topological Gap Opening without Symmetry Breaking from Dynamical Quantum
Correlations

F. Paoletti  International School for Advanced Studies (SISSA), via Bonomea
265, 34136 Trieste, Italy  L. Fanfarillo  CNR-ISC, Via dei Taurini, Rome,
Italy  International School for Advanced Studies (SISSA), via Bonomea 265,
34136 Trieste, Italy  M. Capone  International School for Advanced Studies
(SISSA), via Bonomea 265, 34136 Trieste, Italy  CNR-IOM, Istituto Officina dei
Materiali, Consiglio Nazionale delle Ricerche, Via Bonomea 265, 34136 Trieste,
Italy  A. Amaricci  CNR-IOM, Istituto Officina dei Materiali, Consiglio
Nazionale delle Ricerche, Via Bonomea 265, 34136 Trieste, Italy

###### Abstract

Topological phase transitions are typically associated with the formation of
gapless states. Spontaneous symmetry breaking can lead to a gap opening
thereby obliterating the topological nature of the system. Here we highlight a
completely different destiny for a topological transition in presence of
interaction. Solving a Bernevig-Hughes-Zhang model with local interaction, we
show that dynamical quantum fluctuations can lead to the opening of a gap
without any symmetry breaking. As we vary the interaction and the bare mass of
the model, the continuous gapless topological transition turns into a first-
order one, associated with the presence of massive Dirac fermion at the
transition point showing a Gross-Neveu critical behaviour near the quantum
critical endpoint. We identify the gap opening as a condensed matter analog of
the Coleman-Weinberg mechanism of mass generation.

The discovery of symmetry protected topological phases of matter [1, 2, 3, 4,
5, 6, 7, 8] has enriched the landscape of phase transitions beyond the
conventional Landau paradigm of the symmetry-breaking [9, 10]. In presence of
a given symmetry the possible electronic band structures of an insulator can
be divided into distinct equivalence classes, which can only be connected
through the continuous closure of the energy gap from both sides of the
transition through a topological quantum phase transition (TQPT). The
corresponding formation of symmetry-protected massless Dirac fermions at the
transition is a distinctive feature of the topological insulators [11, 12, 13,
14, 15].

The presence of interactions can change the above scenario and a gapped state
can appear also at the transition. The standard phenomenology for this to
occur requires a spontaneous symmetry breaking (SSB) [4, 16, 17, 18, 19, 20].
Indeed, breaking a continuous symmetry opens a gap in the energy spectrum or,
equivalently, gives a finite mass for the Dirac electrons which is understood
in terms of the Anderson-Higgs mechanism. Clearly, a SSB can lead to break any
of the symmetry protecting the topological state, thus leaving behind a
topologically trivial long-range ordered phase. A similar scenario can be
described within a static mean-field picture in the channel where SSB takes
place.

A great deal of attention has been recently drawn in different fields [21, 22,
23] to novel possible mechanisms of spontaneous mass generation which preserve
the symmetry, beyond the conventional SSB description. Here we show that such
process describes the gap opening for Dirac electrons at the boundary of a
topological insulator. More concretely, we address the question whether or not
electron-electron interactions can drive the formation of a spontaneous mass
for the otherwise gapless electrons at a topological transition. The lack of a
gap closing at the TQPT is expected to change the character of the transition
which becomes necessarily discontinuous despite the symmetries protecting the
topological phase are preserved.

For the sake of definiteness we consider a two-dimensional Bernevig-Hughes-
Zhang (BHZ) model augmented via the inclusion of local electron-electron
interactions that preserve some symmetries of the model. Without interactions,
this model features a TQPT through the formation of a gapless state. The
control parameter of the transition is the energy splitting between two
electronic orbitals playing the role of a mass term. As a consequence the
difference in the occupation of the orbitals, or orbital polarization, is
expected to assume a prominent role. Since the orbital symmetry is broken by
the mass term, the concept of SSB does not apply to the TQPT. Within a mean-
field theory, the interactions simply dress the mass term, shifting the
topological transition without changing its nature with respect to the non-
interacting limit.

In this work we go beyond mean-field using a variational approach including
quantum fluctuations not only of the orbital polarization, but also in the
other particle-hole channels. We demonstrate a new scenario in which a gap
opens at the TQPT without breaking any of the symmetries of the model. We show
explicitly that the quantum fluctuations in the different channels make the
TQPT discontinuous for sufficiently large interactions [24, 25]. The first-
order line ends in a critical endpoint, where we show a Gross-Neveu quantum
critical behavior as a function of the relevant coupling strength [26]. As we
shall discuss in the following, the mechanism we revealed is reminiscent of
the Coleman-Weinberg (CW) theory of mass generation [27, 28].

We solve an interacting BHZ model on a square lattice [1, 29, 2, 30, 24]

$\mathcal{H}=\sum_{\mathbf{k}}\psi^{\dagger}_{\mathbf{k}}H^{0}_{\mathbf{k}}\psi_{\mathbf{k}}+\sum_{\bf
i}\mathcal{H}^{\mathrm{int}}_{\bf i}$ (1)

where
$\psi_{\mathbf{k}}=[c_{\mathbf{k}1\uparrow},c_{\mathbf{k}2\uparrow},c_{\mathbf{k}1\downarrow},c_{\mathbf{k}2\downarrow}]^{T}$
and the operators $c_{\mathbf{k}\alpha\sigma}$ annihilate an electron with
momentum $\mathbf{k}$, orbital $\alpha=1,2$ and spin
$\sigma=\uparrow,\downarrow$. If we define
$\Gamma_{\mu\nu}=\sigma_{\mu}\otimes\tau_{\nu}$ with $\sigma_{\mu}$ and
$\tau_{\nu}$ the Pauli matrices, respectively, in the spin and orbital
subspaces the single-particle Hamiltonian reads:
$H^{0}_{\mathbf{k}}=[M-2t(\cos(k_{x})+\cos(k_{y}))]\Gamma_{03}+\lambda\sin(k_{x})\Gamma_{31}+\lambda\sin(k_{y})\Gamma_{02}$,
where $M\geq 0$ is the energy separation between the two orbitals which plays
the role of the mass term, $t$ and $\lambda$ are the intra- and inter-orbital
hopping amplitudes. The model is invariant under time-reversal ${\cal T}$ and
inversion ${\cal P}$ symmetries, $U(1)$ spin rotation around the $z$ axis [31,
32, 33, 20]. In the following we set our energy unit so that $2t=1$ and focus
on the regime of two electrons per site, i.e. half-filling. The non-
interacting model has a continuous topological transition between a QSHI for
$M<2$ and a trivial BI for $M>2$ through the formation of a gapless Dirac
state at $M=2$ [1].

We assume a generic local interaction which preserves inversion ${\cal P}$ and
$U(1)$ spin symmetry around the $z$ axis [34, 18]:

$\mathcal{H}^{int}_{\bf i}=-\frac{g_{\scriptscriptstyle N}}{2}\hat{N}_{\bf
i}^{2}-\frac{g_{\scriptscriptstyle T}}{2}\hat{T}_{z{\bf
i}}^{2}-\frac{g_{\scriptscriptstyle S}}{2}\hat{S}_{z{\bf
i}}^{2}-\frac{g_{\scriptscriptstyle R}}{2}\hat{R}_{z{\bf i}}^{2}$ (2)

where $\hat{N}_{\bf i}\\!=\\!\tfrac{1}{2}\psi^{+}_{\bf i}\Gamma_{00}\psi_{\bf
i}$ is half of the total occupation per site, $\hat{T}_{z{\bf
i}}\\!=\\!\tfrac{1}{2}\psi^{+}_{\bf i}\Gamma_{03}\psi_{\bf i}$ and
$\hat{S}_{z{\bf i}}\\!=\\!\tfrac{1}{2}\psi^{+}_{\bf i}\Gamma_{30}\psi_{\bf i}$
are, respectively, the $z$ component of the orbital polarization and the spin
operators and $\hat{R}_{z{\bf i}}=\tfrac{1}{2}\psi^{+}_{\bf
i}\Gamma_{33}\psi_{\bf i}$; $\psi_{\bf i}$ is the Fourier transform of
$\psi_{\mathbf{k}}$. In the numerical calculations we will consider
$g_{\scriptscriptstyle N}=-(3U-5J)$, $g_{\scriptscriptstyle T}=U-5J$,
$g_{\scriptscriptstyle S}=U+J$ and $g_{\scriptscriptstyle R}=U-J$ in order to
recover the density-density version of the popular Kanamori-Hubbard [35] model
used in a variety of works to study the interplay between the Hubbard $U$ and
the Hund’s exchange $J$ and their effect on TQPTs [36, 37, 24]. We consider
non-magnetic solutions in order to study $\mathcal{T}$ symmetry preserving
transitions [38, 39]. To simplify the notation in the following we define
$\Gamma_{a=\scriptscriptstyle N,\scriptscriptstyle T,\scriptscriptstyle
S,\scriptscriptstyle R}$ as the set $\frac{1}{2}\Gamma_{\mu\nu=00,03,30,33}$
in order to highlight the different channels.

The starting point of our analysis is to rewrite the partition function of the
interacting model (1) in terms of an effective problem coupled to space- and
time-dependent (real) bosonic fields
$\mathbf{\Delta}_{q}=\Delta_{q}^{a=\scriptscriptstyle N,\scriptscriptstyle
T,\scriptscriptstyle S,\scriptscriptstyle R}$, by performing an
Hubbard–Stratonovich transformation. The partition function reads:
$\mathcal{Z}\equiv e^{-\beta\mathcal{F}}=\int\mathcal{D}\mathbf{\Delta}\
e^{-\beta\mathcal{N}F[\mathbf{\Delta}]}$ in terms of the free energy
functional

$F[\mathbf{\Delta}]=\sum_{aq}\frac{|\Delta^{a}_{q}|^{2}}{2g_{a}}-\frac{1}{\beta\mathcal{N}}\text{Tr}\ln(-{\cal
G}^{-1}_{kq})$ (3)

where $q=(\mathbf{q},i\nu_{m})$, $k=(\mathbf{k},i\omega_{n})$, with
$\mathbf{q},\ \mathbf{k}$ wave vectors in the first Brillouin zone and
$\nu_{m},\ \omega_{n}$ the bosonic and fermionic Matsubara frequencies,
respectively, $\beta$ is the inverse temperature and ${\cal N}$ the total
number of sites. Tr indicates the trace over momentum, frequency, orbital and
spin. Finally ${\cal
G}_{kq}=\left(i\omega_{n}+\mu-H^{0}_{\mathbf{k}}\delta_{\mathbf{k},\mathbf{k}-\mathbf{q}}-V_{q}\right)^{-1}$
is the interacting one-body Green’s function, where
$V_{q}=-\sum_{a}\Delta^{a}_{q}\Gamma_{a}$ is an effective time-dependent
potential depending on $\mathbf{\Delta}_{q}$.

The natural lowest-order approximation of Eq. (3) is a static mean-field (MF)
solution, where the bosonic fields $\Delta^{a}_{q}$ are approximated with
time-independent and spatially uniform quantities. The presence of the mass
term breaks explicitly the symmetry between the orbitals, leading to a finite
values of the orbital polarization $T_{z}$ already in the non-interacting
model. The non-magnetic solution ($\Delta^{\scriptscriptstyle
S}_{\scriptscriptstyle MF}=\Delta^{\scriptscriptstyle R}_{\scriptscriptstyle
MF}=0$) at half-filling ($\Delta^{\scriptscriptstyle N}_{\scriptscriptstyle
MF}=1$), reduces to the single self-consistency equation
$\Delta^{\scriptscriptstyle T}=\frac{g_{\scriptscriptstyle
T}}{\beta\mathcal{N}}\text{Tr}\left(G^{\scriptscriptstyle
MF}_{k}\Gamma_{\scriptscriptstyle T}\right)$ where $G^{\scriptscriptstyle
MF}_{k}=(i\omega_{n}-H^{0}_{\mathbf{k}}+\Delta^{\scriptscriptstyle
T}\Gamma_{\scriptscriptstyle T})^{-1}$ is the MF Green’s function. Thus, the
MF solution simply corrects the mass term $M$ so that the model describes a
continuous TQPT occurring at the critical line
$M-\frac{1}{2}\Delta^{\scriptscriptstyle T}_{\scriptscriptstyle MF}=2$ 111In
this work we consider values of the interaction which do not lead to a change
of sign in the renormalized mass term as reported in the phase diagram of Fig.
1(b). All along this line the energy gap closes through the formation of a
gapless Dirac node at the $\Gamma$ point as in the non-interacting model. As
we show in Fig. 1(a) (dotted line), the orbital polarization smoothly evolves
across the topological transition. Coherently with the above scenario, a
direct inspection of the MF free energy shows only one minimum for every value
of $g_{\scriptscriptstyle T}$ and $M$.

Figure 1: (Color online) (a) Orbital polarization
$\bar{\Delta}^{\scriptscriptstyle T}/g_{\scriptscriptstyle T}$ as a function
of $g_{\scriptscriptstyle T}$ measured with respect to the TQPT point
$g^{*}_{\scriptscriptstyle T}$. The open symbols correspond to the
fluctuation-corrected results. The dotted line is the MF solution for
$M=1.50$. (b) Phase diagram in the $M$-$g_{\scriptscriptstyle T}$ plane
comparing the topological transition line in the two approximations. The solid
lines (black and red) denote a continuous TQPT, while the dashed line (red)
marks a discontinuous one. Data from DMFT are indicated with filled symbols
(gray). The insets A, B and C show the free energy $F$ as a function of the
orbital polarization for the three points marked on the curves. (c) Static and
homogeneous components $\left\llangle{|\Delta^{a}_{q=0}|^{2}}\right\rrangle$
for $a=T,N,S,R$ of the potential fluctuations as a function of
$g_{\scriptscriptstyle T}$ across the topological transition for $M=1.70$.

In this work we overcome the limitations of the MF by approximating the exact
free energy functional with a second-order expansion in the fluctuating
fields, whose coefficients are variationally determined [41]. We underline
that this approach includes effects of higher order beyond the conventional
strategy of analyzing Gaussian fluctuations around the MF solution [42]. Most
importantly, the expectation values of the fields can change with respect to
MF, i.e. $\bar{\mathbf{\Delta}}\neq\bar{\mathbf{\Delta}}_{\scriptscriptstyle
MF}$. Assuming
$\mathbf{\Delta}_{q}\rightarrow\bar{\mathbf{\Delta}}+\mathbf{\Delta}_{q}$ we
can write

$F[\mathbf{\Delta}]\simeq
F^{(2)}[\mathbf{\Delta}]=F[\bar{\mathbf{\Delta}}]+\frac{1}{2}\sum_{abq}\Delta^{a}_{q}A^{ab}_{q}\Delta^{b}_{-q}$
(4)

where the variational principle ${\cal F}\leq{\cal
F}^{(2)}+\left\llangle{F[\mathbf{\Delta}]-F^{(2)}[\mathbf{\Delta}]}\right\rrangle$,
with ${\cal F}^{(2)}\\!=\\!-\frac{1}{\beta}\ln{\int{\cal D}\mathbf{\Delta}\
e^{-\beta\mathcal{N}F^{(2)}[\mathbf{\Delta}]}}$ leads to the new stationary
condition
$\\!\left\llangle{\partial_{\Delta^{a}_{q}}F[\mathbf{\Delta}]}\right\rrangle=0$
and
$A^{ab}_{q}\\!=\\!\left\llangle{\partial_{\Delta^{a}_{q}}\partial_{\Delta^{b}_{-q}}F[\mathbf{\Delta}]}\right\rrangle$.
The symbol $\left\llangle{\cdot}\right\rrangle$ indicates that the averages
over the possible configurations are calculated using the second-order
probability density of the fluctuating field.

The coefficients $A^{ab}$ in the free energy expansion (4) depend on the
averaged dressed Green’s function $\left\llangle{{\cal
G}_{kq}[\mathbf{\Delta}]}\right\rrangle$, which can not be calculated exactly.
In order to circumvent this problem we can introduce an auxiliary potential
$\Sigma_{k}$ implicitly determined by the condition [43, 44]:
$\frac{\partial{\cal F}^{(2)}}{\partial\Sigma_{k}}=0$, which in turn implies
that $\left\llangle{{\cal G}_{kq}[\mathbf{\Delta}]}\right\rrangle$ coincides
with an interacting Green’s function in which $\Sigma_{k}$ plays the role of a
self-energy $G_{k}=\left[i\omega_{n}+\mu-
H_{0}(\mathbf{k})-\Sigma_{k}\right]^{-1}$.

The stationary condition and the expression for $A^{ab}_{q}$ become

$\frac{\bar{\Delta}^{a}}{g_{a}}=\frac{1}{\beta\mathcal{N}}\text{Tr}\left(G_{k}\Gamma_{a}\right);\,A^{ab}_{q}=\frac{\delta_{ab}}{g_{a}}-\chi_{ab}(q).$
(5)

The first expression contains the fluctuation-corrected Green’s function,
hence it leads to corrected values of the $\bar{\mathbf{\Delta}}$, while the
second can be seen as an optimized version of the Random-phase Approximation
(RPA), as we shall discuss in the following [42]. Indeed
$\chi_{ab}(q)=-\tfrac{1}{\beta\mathcal{N}}\sum_{k}\text{Tr}\left[G_{k}\Gamma_{a}G_{k+q}\Gamma_{b}\right]$
is the susceptibility matrix in the space of the different channels which, in
presence of odd hybridization between the orbitals, has diagonal structure
$\chi_{aa}\delta_{ab}$. Moreover, the symmetries of the interaction ensure
that $\chi_{\scriptscriptstyle S\scriptscriptstyle S}=\chi_{\scriptscriptstyle
N\scriptscriptstyle N}$ and $\chi_{\scriptscriptstyle R\scriptscriptstyle
R}=\chi_{\scriptscriptstyle T\scriptscriptstyle T}$. $\Sigma_{k}$ can be
written explicitly up to second order as

$\Sigma_{k}=\bar{V}+\sum_{qa}G_{k-q}\left\llangle{|\Delta^{a}_{q}|^{2}}\right\rrangle$
(6)

where

$\left\llangle{|\Delta^{a}_{q}|^{2}}\right\rrangle=\tfrac{1}{\beta\mathcal{N}}\left[\frac{1}{g^{-1}_{a}-\chi_{aa}(q)}-g_{a}\right].$
(7)

The diagonal form in the channel index allows us to analyze the contribution
of each fluctuating term of the interaction to the potential $\Sigma_{k}$
[45]. We emphasize that, formally, the second term in (6) plays the same role
of the one-loop quantum correction of the effective Coleman-Weinberg potential
[27].

The Eqs.(5) and (6) provide a closed system of non-linear equations for
$\Sigma_{k}$ and the bosonic fields $\bar{\mathbf{\Delta}}$. Previous studies
suggest that the interaction effects on the TQPT are mainly local and that
non-local fluctuation play a minor role [46]. Thus, in order to further
simplify the treatment in the following we will assume a local
$\Sigma_{k}\simeq\Sigma(i\omega_{n})$. We solve this system iteratively using
a linear mixing algorithm which typically converges in $10-20$ steps. The BZ
is discretized with a linear grid of $20\times 20$ points and the Matsubara
axis with $L=8192$ frequencies using an effective inverse temperature
$\beta=500$. The convolution in Eq. (6) is evaluated using a FFT algorithm. We
discuss the results obtained for $J/U=1/8$ and $\lambda=0.3$.

In Fig.1(a) we show the evolution of the orbital polarization, obtained from
the self-consistent value of the bosonic field
$\bar{T}_{z}=\bar{\Delta}^{\scriptscriptstyle T}/g_{\scriptscriptstyle T}$.
The behavior at the transition point $g^{*}_{\scriptscriptstyle T}(M)$ changes
qualitatively according to the value of the bare mass $M$. For a value close
to $M=2$, i.e. the non-interacting transition point, the orbital polarization
is continuous with respect to the increasing interaction
$g_{\scriptscriptstyle T}$. This corresponds to a smooth modification of the
BI into a non-trivial insulator through the formation of a gapless state at
the TQPT. Starting from a farther point, the orbital polarization displays a
critical behavior at the transition characterized by a divergent
susceptibility $\partial_{M}\bar{T}_{z}$. Beyond this point, for any value of
$M$, the orbital polarization is characterized by a discontinuous evolution
across the topological transition. This is in stark contrast with the
continuous behavior obtained in MF for the same value of the mass, see Fig.
1(a), and it agrees with previous results obtained via Dynamical Mean-Field
Theory (DMFT) [24]. The agreement with DMFT is indeed even quantitative, as
shown by the data reported in Fig. 1(b).

These results can be summarised in a phase diagram in the plane
$g_{\scriptscriptstyle T}$-$M$ (see Fig. 1(b)) where we compare the MF and
fluctuation-corrected results for the TQPT. The two transition lines remain
close for small values of the interaction $g_{\scriptscriptstyle T}$.
Accordingly, the free energy functional $F$ displays a single minimum as a
function of $\Delta^{\scriptscriptstyle T}$ (inset A in Fig. 1(b)). However,
upon increasing the interaction strength the two curves start deviating
significantly signalling a crucial impact of the fluctuations.

A direct information about the contribution of the fluctuations in the
different channels is reported in Fig. 1(c). The fast increasing behavior with
the interaction $g_{\scriptscriptstyle T}$ in all the channels stops at the
topological transition towards the QSHI, where these quantities display a
discontinuous drop and a successive slow increase. The terms
$\left\llangle{|\Delta^{a}_{q}|^{2}}\right\rrangle$ enter, through (7), in
$\Sigma(i\omega_{n})$ giving it a dynamical nature which significantly
deviates from its static MF form $\bar{V}=-\bar{\Delta}^{\scriptscriptstyle
T}\Gamma_{\scriptscriptstyle T}$. This results in a crucial shift of the self-
consistent saddle point value of the bosonic fields. Moreover, as discussed
above, while the MF always describes a continuous transition, in the corrected
theory the boundary line is continuous up to a critical value
$g^{c}_{\scriptscriptstyle T}$ of the interaction beyond which it becomes of
first order. This reflects in the behavior of the free energy near the TQPT
point in the intermediate to strong coupling regime, i.e.
$g_{\scriptscriptstyle T}>g_{\scriptscriptstyle T}^{c}$. In the insets B and C
of Fig. 1(b) we compare the free energies of a QSHI state near the topological
transition for, respectively, the MF and the fluctuation-corrected
approximation, where two minima are found: A stable QSHI and a metastable BI.

We thus find that a QCP separates the continuous from the discontinuous regime
on the topological transition line, where we also found a divergent orbital
susceptibility. Indeed the uniform orbital susceptibility
$\mathcal{X}_{\scriptscriptstyle
T}=\partial_{M}\bar{\Delta}^{\scriptscriptstyle T}_{q=0}$ reads, using Eq.
(5),

$\mathcal{X}_{\scriptscriptstyle T}=\frac{-2g_{\scriptscriptstyle
T}\chi_{{\scriptscriptstyle T}{\scriptscriptstyle
T}}(0)}{1-g_{\scriptscriptstyle T}\chi_{{\scriptscriptstyle
T}{\scriptscriptstyle T}}(0)+\lambda},$ (8)

which reminds of the RPA result with a correction $\lambda$ that stems from
the implicit dependence of $\Sigma(i\omega_{n})$ on the orbital polarization.
This quantity also accounts for the contributions of all the other channels of
the interaction through the expression of $\Sigma(i\omega_{n})$, see Eq. (6)
[47]. At the MF level we find
$\Sigma(i\omega_{n})=-\bar{\Delta}^{\scriptscriptstyle
T}\Gamma_{\scriptscriptstyle T}$, so that $\lambda=0$ and
$\mathcal{X}_{\scriptscriptstyle T}$ reduces to the RPA form. For a Hubbard-
Kanamori interaction the RPA $\mathcal{X}_{\scriptscriptstyle T}$ diverges
only either for negative $U$ or negative $M$ (which lead to different
physics), in agreement with the continuous TQPT we always find. Note that this
result holds true also when all the coupling constants but
$g_{\scriptscriptstyle T}$ vanish and the interaction reduces to
$\frac{g_{\scriptscriptstyle T}}{2}\hat{T}_{z}^{2}$ [48].

Figure 2: (Color online) (a) The denominator of Eq. (8) as a function of the
orbital polarization $\bar{\Delta}^{\scriptscriptstyle
T}/g_{\scriptscriptstyle T}$ across the TQPT in the fluctuation-corrected
approximation and for different values of the bare mass $M$. The narrow grey
stripe indicates the minimum with its numerical uncertainty. (b) The gap
$\Lambda$ along the TQPT line as a function of the interaction
$g_{\scriptscriptstyle T}$. MF is the dotted grey line, while fluctuation-
corrected results are indicated by open symbols and solid line. The (red)
dashed line is a linear fit $A(g_{\scriptscriptstyle
T}-g^{c}_{\scriptscriptstyle T})^{\beta=1}$ ($A\simeq 0.042$) of the critical
behavior.

In order to compute $\lambda$ we consider the zero-frequency limit of
$\Sigma(i\omega_{n})$ where we obtain $\lambda\simeq\chi_{{\scriptscriptstyle
T}{\scriptscriptstyle T}}\partial_{\bar{T}_{z}}\Lambda$. Since
$\partial_{\bar{T}_{z}}\Lambda<0$, this correction is negative and tends to
enhance $\mathcal{X}_{\scriptscriptstyle T}$. For a fixed value of $M$ the
TQPT corresponds to the maximum of the response function which connects the
continuous transition to a Widom line [49, 50, 51]. The minimum of the
denominator in Eq. (8) approaches zero when we reach the QCP, as shown in Fig.
2(a). As we discussed above, in a non-interacting TQPT the spectral gap closes
at the transition. We now compute the gap in our scheme from the zero-
frequency limit of the self-energy $\Lambda={\rm
Re}\Sigma(i\omega_{n}\\!\rightarrow\\!0)-\bar{V}$. In Fig. 2(b) we report the
behavior of $\Lambda$ at the TQPT as a function of the interaction strength
$g_{\scriptscriptstyle T}$. While in MF the gap is always zero, including the
fluctuations we find a finite gap above the QCP ($g_{\scriptscriptstyle
T}>g^{c}_{\scriptscriptstyle T}$). A finite value of the gap $\Lambda$
corresponds to give a mass to the Dirac fermions at the boundary line. This is
consistent with a spontaneous symmetric mass generation process [22, 23, 21].
The presence of such finite gap (or mass) makes it impossible to continuously
connect the trivial with the non-trivial phase and leads to a first-order
TQPT. In addition, we find numerically that the critical behavior near the QCP
falls in the Gross-Neveu universality class [22, 26], with an estimated
critical exponent $\beta\simeq 1$.

In this work, using a non-perturbative analytical approach to include
interactions in the BHZ model, we have demonstrated the crucial role of
fluctuations in the different local particle-hole channels to qualitatively
change the nature of the TQPT. Within a mean-field the interactions only lead
to a renormalization of the bare mass of the model (coupled to the orbital
polarization) so that the TQPT has the same character of the non-interacting
model. Within our approach the fluctuation contributions change the MF
parameters and lead to a discontinuous TQPT for large interactions with a QCP
separating the continuous and discontinuous branches. This effect is
intrinsically related to a spontaneous gap opening (mass formation) for the
otherwise gapless Dirac nodes at TQPT point without any symmetry breaking. The
gap follows a Gross-Neveu critical behavior. This process of spontaneous mass
generation takes place through a condensed matter analog of the CW mechanism
in which one-loop quantum fluctuations lead to a mass without symmetry
breaking [27]. We expect that our mechanism can be applied also to other
models for topological phase transitions, but also to a wider class of
phenomena. A natural example is that of Lifshitz transitions in interacting
electronic systems, where the continuous deformation of the Fermi surface
topology characteristic of non-interacting systems is expected to share the
same destiny of the TQPT, i.e., to become discontinuous for large
interactions.

## Acknowledgments –

F.P. is grateful to P.Coleman and G.Perosa for illuminating discussions
specially concerning Ref. [27]. We thank M. Fabrizio, G. Mazza, G.
Sangiovanni, J.C. Budich, B. Trauzettel for useful discussions. We acknowledge
financial support of MUR via PRIN 2017 (Prot. 20172H2SC4 005), PRIN 2020
(Prot. 2020JLZ52N 002) programs, PRIN 2022 (Prot. 20228YCYY7), National
Recovery and Resilience Plan (NRRP) MUR Project No. PE0000023-NQSTI and
ICSC–Centro Nazionale di Ricerca in High Performance Computing, Big Data and
Quantum Computing, funded by European Union – NextGenerationEU (Grant number
CN00000013) - Mission 4 Component 2 Investments 1.3 and 1.4.

## References

  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Quantum Spin Hall Effect and Topological Phase Transition in HgTe Quantum Wells, Science 314, 1757 (2006).
  * Qi and Zhang [2010] X.-L. Qi and S.-C. Zhang, The quantum spin Hall effect and topological insulators, Physics Today 63, 33 (2010).
  * Moore [2010] J. E. Moore, The birth of topological insulators, Nature 464, 194 (2010), 10.1038/nature08916.
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane,  _Colloquium_ : Topological insulators. , Rev. Mod. Phys 82, 3045 (2010).
  * Bernevig and Hughes [2013] B. Bernevig and T. Hughes,  _Topological Insulators and Topological Superconductors_ (Princeton University Press, 2013).
  * Wen [2017] X.-G. Wen, Colloquium: Zoo of quantum-topological phases of matter, Rev. Mod. Phys. 89, 041004 (2017).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Topological quantum chemistry, Nature 547, 298 (2017).
  * Rachel [2018] S. Rachel, Interacting topological insulators: a review, Reports on Progress in Physics 81, 116501 (2018).
  * Nambu [2009] Y. Nambu, Nobel lecture: Spontaneous symmetry breaking in particle physics: A case of cross fertilization, Rev. Mod. Phys. 81, 1015 (2009).
  * Zhang  _et al._ [2009] Y. Zhang, Y. Ran, and A. Vishwanath, Topological insulators in three dimensions from spontaneous symmetry breaking, Phys. Rev. B 79, 245331 (2009).
  * Zhou  _et al._ [2006] S. Y. Zhou, G. H. Gweon, J. Graf, A. V. Fedorov, C. D. Spataru, R. D. Diehl, Y. Kopelevich, D. H. Lee, S. G. Louie, and A. Lanzara, First direct observation of Dirac fermions in graphite, Nat. Phys. 2, 595 (2006), arXiv:0608069 [cond-mat] .
  * Hsieh D.  _et al._ [2008] Hsieh D., Qian D., Wray L., Xia Y., Hor Y. S., Cava R. J., and Hasan M. Z., A topological Dirac insulator in a quantum spin Hall phase, Nature 452, 970 (2008), 10.1038/nature06843.
  * Wehling  _et al._ [2014] T. O. Wehling, A. M. Black-Schaffer, and A. V. Balatsky, Dirac materials, Adv. Phys. 63, 1 (2014), arXiv:1405.5774 .
  * Young and Kane [2015] S. M. Young and C. L. Kane, Dirac Semimetals in Two Dimensions, Phys. Rev. Lett. 155, 126803 (2015).
  * Zahid Hasan  _et al._ [2015] M. Zahid Hasan, S.-Y. Xu, and M. Neupane, Topological Insulators, Topological Dirac semimetals, Topological Crystalline Insulators, and Topological Kondo Insulators, in  _Topol. Insul. Fundam. Perspect._ (Wiley-VCH Verlag, 2015) wiley ed., Chap. 4, pp. 55–100, arXiv:1406.1040 .
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Ezawa  _et al._ [2013] M. Ezawa, Y. Tanaka, and N. Nagaosa, Topological Phase Transition without Gap Closing., Sci. Rep. 3, 2790 (2013).
  * Rachel [2016] S. Rachel, Quantum phase transitions of topological insulators without gap closing, Journal of Physics: Condensed Matter 28, 405502 (2016).
  * Jia  _et al._ [2022] Y. Jia, P. Wang, C.-L. Chiu, Z. Song, G. Yu, B. Jäck, S. Lei, S. Klemenz, F. A. Cevallos, M. Onyszczak, N. Fishchenko, X. Liu, G. Farahi, F. Xie, Y. Xu, K. Watanabe, T. Taniguchi, B. A. Bernevig, R. J. Cava, L. M. Schoop, A. Yazdani, and S. Wu, Evidence for a monolayer excitonic insulator, Nature Physics 18, 87 (2022).
  * Amaricci  _et al._ [2023] A. Amaricci, G. Mazza, M. Capone, and M. Fabrizio, Exciton condensation in strongly correlated quantum spin hall insulators, Phys. Rev. B 107, 115117 (2023).
  * You  _et al._ [2018] Y.-Z. You, Y.-C. He, C. Xu, and A. Vishwanath, Symmetric fermion mass generation as deconfined quantum criticality, Phys. Rev. X 8, 011026 (2018).
  * Slagle  _et al._ [2015] K. Slagle, Y.-Z. You, and C. Xu, Exotic quantum phase transitions of strongly interacting topological insulators, Phys. Rev. B 91, 115121 (2015).
  * Wang and You [2022] J. Wang and Y.-Z. You, Symmetric mass generation, Symmetry 14, 10.3390/sym14071475 (2022).
  * Amaricci  _et al._ [2015] A. Amaricci, J. C. Budich, M. Capone, B. Trauzettel, and G. Sangiovanni, First-order character and observable signatures of topological quantum phase transitions, Phys. Rev. Lett. 114, 185701 (2015).
  * Yamaji  _et al._ [2007] Y. Yamaji, T. Misawa, and M. Imada, Quantum criticalities induced by lifshitz transitions, Journal of Magnetism and Magnetic Materials 310, 838 (2007), proceedings of the 17th International Conference on Magnetism.
  * Liu  _et al._ [2021] Y. Liu, Z. Wang, T. Sato, W. Guo, and F. F. Assaad, Gross-neveu heisenberg criticality: Dynamical generation of quantum spin hall masses, Physical Review B 104, 10.1103/physrevb.104.035107 (2021).
  * Coleman and Weinberg [1973] S. Coleman and E. Weinberg, Radiative corrections as the origin of spontaneous symmetry breaking, Phys. Rev. D 7, 1888 (1973).
  * Rau  _et al._ [2018] J. G. Rau, P. A. McClarty, and R. Moessner, Pseudo-goldstone gaps and order-by-quantum disorder in frustrated magnets, Phys. Rev. Lett. 121, 237201 (2018).
  * Wu  _et al._ [2006] C. Wu, B. A. Bernevig, and S.-C. Zhang, Helical liquid and the edge of quantum spin hall systems, Phys. Rev. Lett. 96, 106401 (2006).
  * Qian  _et al._ [2014] X. Qian, J. Liu, L. Fu, and J. Li, Quantum spin hall effect in two-dimensional transition metal dichalcogenides, Science 346, 1344 (2014).
  * Blason and Fabrizio [2020] A. Blason and M. Fabrizio, Exciton topology and condensation in a model quantum spin hall insulator, Phys. Rev. B 102, 035146 (2020).
  * Budich  _et al._ [2014] J. C. Budich, B. Trauzettel, and P. Michetti, Time reversal symmetric topological exciton condensate in bilayer hgte quantum wells, Phys. Rev. Lett. 112, 146405 (2014).
  * Mazza  _et al._ [2020] G. Mazza, M. Rösner, L. Windgätter, S. Latini, H. Hübener, A. J. Millis, A. Rubio, and A. Georges, Nature of symmetry breaking at the excitonic insulator transition: ${\mathrm{ta}}_{2}{\mathrm{nise}}_{5}$, Phys. Rev. Lett. 124, 197601 (2020).
  * Hohenadler and Assaad [2013] M. Hohenadler and F. F. Assaad, Correlation effects in two-dimensional topological insulators, Journal of Physics: Condensed Matter 25, 143201 (2013).
  * Georges  _et al._ [2013] A. Georges, L. de’ Medici, and J. Mravlje, Strong Correlations from Hund’s Coupling., Annu. Rev. Condens. Matter Phys. 45, 137 (2013).
  * Werner and Millis [2007] P. Werner and A. J. Millis, High-Spin to Low-Spin and Orbital Polarization Transitions in Multiorbital Mott Systems, Phys. Rev. Lett. 99, 126405 (2007).
  * Budich  _et al._ [2012] J. C. Budich, R. Thomale, G. Li, M. Laubach, and S.-C. Zhang, Fluctuation-induced topological quantum phase transitions in quantum spin-Hall and anomalous-Hall insulators, Phys. Rev. B 86, 201407 (2012).
  * Amaricci  _et al._ [2016] A. Amaricci, J. C. Budich, M. Capone, B. Trauzettel, and G. Sangiovanni, Strong correlation effects on topological quantum phase transitions in three dimensions, Phys. Rev. B 93, 235112 (2016).
  * Amaricci  _et al._ [2017] A. Amaricci, L. Privitera, F. Petocchi, M. Capone, G. Sangiovanni, and B. Trauzettel, Edge state reconstruction from strong correlations in quantum spin hall insulators, Phys. Rev. B 95, 205120 (2017).
  * Note [1] In this work we consider values of the interaction which do not lead to a change of sign in the renormalized mass term.
  * Hertz and Klenin [1974] J. A. Hertz and M. A. Klenin, Fluctuations in itinerant-electron paramagnets, Phys. Rev. B 10, 1084 (1974).
  * Klein  _et al._ [2018] Y. Klein, M. Casula, D. Santos-Cottin, A. Audouard, D. Vignolles, G. Fève, V. Freulon, B. Plaçais, M. Verseils, H. Yang, L. Paulatto, and A. Gauzzi, Importance of nonlocal electron correlation in the BaNiS2 semimetal from quantum oscillations studies, Phys. Rev. B 97, 075140 (2018).
  * Kakehashi [2002] Y. Kakehashi, Dynamical coherent-potential approximation to the magnetism in a correlated electron system, Phys. Rev. B 65, 184420 (2002).
  * Melnikov  _et al._ [2011] N. B. Melnikov, B. I. Reser, and V. I. Grebennikov, Extended dynamic spin-fluctuation theory of metallic magnetism, Journal of Physics: Condensed Matter 23, 276003 (2011).
  * Gunnarsson  _et al._ [2015] O. Gunnarsson, T. Schäfer, J. LeBlanc, E. Gull, J. Merino, G. Sangiovanni, G. Rohringer, and A. Toschi, Fluctuation diagnostics of the electron self-energy: Origin of the pseudogap physics, Physical Review Letters 114, 10.1103/physrevlett.114.236402 (2015).
  * Crippa  _et al._ [2021] L. Crippa, A. Amaricci, S. Adler, G. Sangiovanni, and M. Capone, Local versus nonlocal correlation effects in interacting quantum spin hall insulators, Phys. Rev. B 104, 235117 (2021).
  * Fanfarillo  _et al._ [2012] L. Fanfarillo, L. Benfatto, and C. Castellani, Current-current fermi-liquid corrections to the superconducting fluctuations on conductivity and diamagnetism, Phys. Rev. B 85, 024507 (2012).
  * Roy  _et al._ [2016] B. Roy, P. Goswami, and J. D. Sau, Continuous and discontinuous topological quantum phase transitions, Phys. Rev. B 94, 041101 (2016).
  * Xu  _et al._ [2005] L. Xu, P. Kumar, S. V. Buldyrev, S.-H. Chen, P. H. Poole, F. Sciortino, and H. E. Stanley, Relation between the widom line and the dynamic crossover in systems with a liquid–liquid phase transition, Proceedings of the National Academy of Sciences 102, 16558 (2005), https://www.pnas.org/doi/pdf/10.1073/pnas.0507870102 .
  * Simeoni  _et al._ [2010] G. G. Simeoni, T. Bryk, F. A. Gorelli, M. Krisch, G. Ruocco, M. Santoro, and T. Scopigno, The widom line as the crossover between liquid-like and gas-like behaviour in supercritical fluids, Nature Physics 6, 503 (2010).
  * Sordi  _et al._ [2012] G. Sordi, P. Sémon, K. Haule, and A. M. S. Tremblay, Pseudogap temperature as a widom line in doped mott insulators, Scientific Reports 2, 547 (2012).


#  Performance Evaluation and Acceleration of the QTensor Quantum Circuit
Simulator on GPUs

Danylo Lykov15, Angela Chen25, Huaxuan Chen35, Kristopher Keipert4, Zheng
Zhang2, Tom Gibbs4, and Yuri Alexeev1  1Computational Science Division,
Argonne National Laboratory, Lemont, IL 60439  2Department of Electrical and
Computer Engineering, University of California, Santa Barbara, CA 93106
3Northwestern University, Evanston, IL 60208  4NVIDIA Corporation, Santa
Clara, CA 95050  5 These authors contributed equally to the paper

(Started 6/3/21)

###### Abstract

This work studies the porting and optimization of the tensor network simulator
QTensor on GPUs, with the ultimate goal of simulating quantum circuits
efficiently at scale on large GPU supercomputers. We implement NumPy, PyTorch,
and CuPy backends and benchmark the codes to find the optimal allocation of
tensor simulations to either a CPU or a GPU. We also present a dynamic mixed
backend to achieve optimal performance. To demonstrate the performance, we
simulate QAOA circuits for computing the MaxCut energy expectation. Our method
achieves $176\times$ speedup on a GPU over the NumPy baseline on a CPU for the
benchmarked QAOA circuits to solve MaxCut problem on a 3-regular graph of size
30 with depth $p=4$.

##  I Introduction

Quantum information science (QIS) has a great potential to speed up certain
computing problems such as combinatorial optimization and quantum simulations
[1]. The development of fast and resource-efficient quantum simulators to
classically simulate quantum circuits is the key to the advancement of the QIS
field. For example, simulators allow researchers to evaluate the complexity of
new quantum algorithms and to develop and validate the design of new quantum
circuits. Another important application is to validate quantum supremacy and
advantage claims.

One can simulate quantum circuits on classical computers in many ways. The
major types of simulation approaches are full amplitude-vector evolution [2,
3, 4, 5], the Feynman paths approach [6], linear algebra open system
simulation [7], and tensor network contractions [8, 9, 10]. These techniques
have advantages and disadvantages. Some are better suited for small numbers of
qubits and high-depth quantum circuits, while others are better for circuits
with a large number of qubits but small depth. Some are also tailored toward
the accuracy of simulation of noise in quantum computers.

For shallow quantum circuits the state-of-the-art technique to simulate
quantum circuits is currently arguably the tensor network contraction method
because of the memory efficiency for the method relative to state vector
methods that scale by $2^{N}$, where $N$ is the number of qubits. This
effectively limits the state vector methods to quantum circuits with less than
50 qubits. The challenge with the tensor network methods is determining the
optimal contraction order, which is known to be an NP-complete problem [8]. We
choose to focus on the simulation of the Quantum Approximate Optimization
Algorithm (QAOA) [11] given its importance to machine learning and its
suitability for the current state of the art with noisy intermediate state
quantum computers that generally work with circuits of short depth.

In this work we ported and optimized the tensor network quantum simulator
QTensor to run efficiently on GPUs, with the eventual goal to simulate large
quantum circuits on the modern and upcoming supercomputers. In particular, we
benchmarked QTensor on a NVIDIA DGX-2 server with a V100 accelerator using the
CUDA version 11.0. The performance is shown for the full expectation value
simulation of the QAOA MaxCut problem on a 3-regular graph of size 30 with
depth $p=4$.

##  II Methodology

###  II-A QAOA Overview

The Quantum Approximate Optimization Algorithm is a variational quantum
algorithm that combines a parameterized ansatz state preparation with a
classical outer-loop algorithm that optimizes the ansatz parameters. QAOA is
used for approximate solution of binary optimization problems [12]. A solution
to the optimization problem is obtained by measuring the ansatz state on a
quantum device. The quality of the QAOA solution depends on the depth of the
quantum circuit that generated the ansatz and the quality of parameters for
the ansatz state.

A binary combinatorial optimization problem is defined on a space of binary
strings of length $N$ and has $m$ clauses. Each clause is a constraint
satisfied by some assignment of the bit string. QAOA maps the combinatorial
optimization problem onto a $2^{N}$-dimensional Hilbert space with
computational basis vectors $\ket{z}$ and encodes $C(z)$ as a quantum operator
$\hat{C}$ diagonal in the computational basis. One of the most widely used
benchmark combinatorial optimization problems is MaxCut, which is defined on
an undirected unweighted graph. The goal of the MaxCut problem is to find a
partition of the graph’s vertices into two complementary sets such that the
number of edges between the sets is maximized. It has been shown in [12] that
on a 3-regular graph, QAOA with $p=1$ produces a solution with an
approximation ratio of at least 0.6924.

A graph $G=(V,E)$ of $N=|V|$ vertices and $m=|E|$ edges can be encoded into a
MaxCut cost operator over $N$ qubits by using $m$ two-qubit gates.

$\hat{C}=\frac{1}{2}\sum_{(ij)\in
E}1-\hat{\sigma^{z}}_{i}\hat{\sigma^{z}}_{j}$ (1)

The QAOA ansatz state $\ket{{\vec{\gamma},\vec{\beta}}}$ is prepared by
applying $p$ layers of evolution unitaries that correspond to the cost
operator $\hat{C}$ and a mixing operator $\hat{B}=\sum_{i\in
V}\hat{\sigma}^{x}_{i}$. The initial state is the equally weighted
superposition state and maximal eigenstate of $\hat{B}$.

$\ket{{\vec{\gamma},\vec{\beta}}}_{p}=\prod_{k=1}^{p}e^{-i\beta_{k}\hat{B}}e^{-i\gamma_{k}\hat{C}}\ket{+}$
(2)

The parameterized quantum circuit (2) is called the _QAOA ansatz_. We refer to
the number of alternating operator pairs $p$ as the _QAOA depth_.

The solution to the combinatorial optimization problem is obtained by
measuring the QAOA ansatz. The expected quality of this solution is an
expectation value of the cost operator in this state.

$\langle
C\rangle_{p}=\bra{{\vec{\gamma},\vec{\beta}}}_{p}C\ket{{\vec{\gamma},\vec{\beta}}}_{p}$

The expectation value can be minimized with respect to parameters
${\vec{\gamma},\vec{\beta}}$. The optimization of ${\vec{\gamma},\vec{\beta}}$
can be performed by using classical computation or by varying the parameters
and sampling many bitstrings from a quantum computer to estimate the
expectation value. Acceleration of the optimal parameters search for a given
QAOA depth $p$ is the focus of many approaches aimed at demonstrating the
quantum advantage. Examples include such methods as warm- and multistart
optimization [13, 14], problem decomposition [15], instance structure analysis
[16], and parameter learning [17].

In this paper we focus on application of a classical quantum circuit simulator
QTensor to the problem of finding the expectation value $\langle
C\rangle_{p}$.

###  II-B Tensor Network Contractions

Calculation of an expectation value of some observable in a given state
generated by some quantum circuit can be done efficiently by using a tensor
network approach. In contrast to state vector simulators, which store the full
state vector of size $2^{N}$, QTensor maps a quantum circuit to a tensor
network. Each quantum gate of the circuit is converted to a tensor. An
expectation value
$\braket{\phi}{\hat{C}}{\phi}=\braket{\psi}{\hat{U}^{\dagger}\hat{C}\hat{U}}{\psi}$
is then simulated by contracting the corresponding tensor network. For more
details on how a quantum circuit is converted to a tensor network, see [18,
19].

A tensor network is a collection of tensors, which in turn have a collection
of indices, where tensors share some indices with each other. To contract a
tensor network, we create an ordered list of tensor buckets. Each bucket (a
collection of tensors) corresponds to a tensor index, which is called _bucket
index_. Buckets are then contracted one by one. The contraction of a bucket is
performed by summing over the bucket index, and the resulting tensor is then
appended to the appropriate bucket. The number of unique indices in aggregate
indices of all bucket tensors is called a _bucket width_. The memory and
computational resources of a bucket contraction scale exponentially with the
associated bucket width. For more information on tensor network contraction,
see [20, 21, 22]. If some observable $\hat{\Sigma}$ acts on a small subset of
qubits, most of the gates in the quantum circuit $\hat{U}$ cancel out when
evaluating the expectation value. The cost QAOA operator $\hat{C}$ is a sum of
$m$ such terms, each of which could be viewed as a separate observable. Each
term generates a _lightcone_ —a subset of the problem that generates a tensor
network representing the contribution to the cost expectation value.

The expectation value of the cost for the graph $G$ and MaxCut QAOA depth $p$
is then

$\displaystyle\braket{C}_{p}({\vec{\gamma},\vec{\beta}})$
$\displaystyle=\braket{{\vec{\gamma},\vec{\beta}}}{\hat{C}}{{\vec{\gamma},\vec{\beta}}}$
$\displaystyle=\braket{{\vec{\gamma},\vec{\beta}}}{\sum_{jk\in
E}\frac{1}{2}(1-\hat{\sigma^{z}}_{j}\hat{\sigma^{z}}_{k})}{{\vec{\gamma},\vec{\beta}}}$
$\displaystyle=\frac{|E|}{2}-\frac{1}{2}\sum_{jk\in
E}\braket{{\vec{\gamma},\vec{\beta}}}{\hat{\sigma^{z}}_{j}\hat{\sigma^{z}}_{k}}{{\vec{\gamma},\vec{\beta}}}$
$\displaystyle\equiv\frac{|E|}{2}-\frac{1}{2}\sum_{jk\in
E}e_{jk}({\vec{\gamma},\vec{\beta}}),$

where $e_{jk}$ is an individual edge contribution to the total cost function.
Note that the observable in the definition of $e_{jk}$ is local to only two
qubits; therefore most of the gates in the circuit that generates the state
$\ket{{\vec{\gamma},\vec{\beta}}}$ cancel out. The circuit after the
cancellation is equivalent to calculating
$\hat{\sigma^{z}}_{j}\hat{\sigma^{z}}_{k}$ on a subgraph $S$ of the original
graph $G$. These subgraphs can be obtained by taking only the edges that are
incident from vertices at a distance $p-1$ from the vertices $j$ and $k$. The
full calculation of $E_{G}({\vec{\gamma},\vec{\beta}})$ requires evaluation of
$|E|$ tensor networks, each representing the value
$e_{jk}({\vec{\gamma},\vec{\beta}})$ for $jk\in E$.

###  II-C Merged Indices Contraction

Since the contraction in the bucket elimination algorithm is executed one
index at a time, the ratio of computational operations to memory read/write
operations is small. This ratio is also called the operational intensity or
arithmetic intensity. Having small arithmetic intensity hurts the performance
in terms of FLOPs: for each floating-point operation calculated there are
relatively many I/O operations, which are usually slower. For example, to
calculate one element of the resulting matrix in a matrix multiplication
problem, one needs to read $2N$ elements and perform $4N$ operations. The size
of the resulting matrix is similar to the input matrices. In contrast, when
calculating an outer product of two vectors, the size of the resulting matrix
is much larger than the combined size of the input vectors; each element
requires two reads and only one floating-point operation.

To mitigate this limitation, we develop an approach for increasing the
arithmetic intensity, which we call merged indices. The essence of the
approach is to combine several buckets and contract their corresponding
indices at once, thus having smaller output size and larger arithmetic
intensity. We have a group of circuit contraction backends that all use this
approach.

For the merged backend group, we order the buckets first and then find the
mergeable indices before performing the contraction. We list the set of
indices of tensors in each bucket and then merge the buckets if the set of
indices of one bucket is a subset of the other. We benchmark the sum of the
total time needed for the merged indices contraction and compare it with the
unmerged baseline results. We call this group the “merged” group and the
baseline the “unmerged” group.

###  II-D CPU-GPU Hybrid Backend

The initial tensor network contains only very small tensors of at most 16
elements (4 dimensions of size 2). We observe that the contraction sequence
obtained by our ordering algorithm results in buckets of small width for first
80% of contraction steps. Only after all small buckets are contracted,
sequence we start to contract large buckets. The GPUs usually perform much
better when processing large amount of data. We observe this behaviour in our
benchmarks on Figure 1. We therefore implement a mix backend which uses both
CPU and GPU. It combines a CPU backend and a GPU backend by dispatching the
contraction procedure to appropriate backend.

The mix or the hybrid backend uses the bucket width, which is determined by
the number of unique indices in a bucket, to allocate the correct device for
such a bucket to be computed. The threshold between the CPU backend and the
GPU backend is determined by a trial program. This program runs a small
circuit, which is used for all backends for testing, separately on a GPU
backend and a CPU backend. After the testing is complete, it iterates through
all bucket widths and checks whether at this bucket width the GPU takes less
time or not. If it finds the bucket width at which the GPU is faster, it will
output that bucket width, and the user can use this width when creating the
hybrid backend in the actual simulation. In the actual simulation, if the
bucket width is smaller than the threshold, the hybrid backend will allocate
this bucket to the CPU and will allocate it to the GPU if the width is
greater.

Since we don’t contract buckets of large width on CPU, the resulting tensors
are rather small, on the order of 1,000s of bytes. The time for data transfer
in this case is considered negligible and is not measured in our code. The
large tensors start to appear from contractions that combine these small
tensors after all the data is moved to GPU.

###  II-E Datasets for Synthetic Benchmarks

Tensor network contraction is a complex procedure that involves many
inhomogeneous operations. Since we are interested in achieving the maximum
performance of the simulations, it is beneficial to compare the FLOPs
performance to several more relevant benchmarking problems. We select several
problems for this task:

  1. 1.

Square matrix multiplication, the simplest benchmark problem which serves as
an upper bound for our FLOP performance;

  2. 2.

Pairwise tensor contractions with a small number of large dimensions and fixed
contraction structure;

  3. 3.

Pairwise tensor contractions with a large number of dimensions of size 2 and
permuted indices;

  4. 4.

Bucket contraction of buckets that are produced by actual expectation value
calculation;

  5. 5.

Full circuit contraction which takes into account buckets of large and small
width.

By gradually adding complexity levels to the benchmark problems and evaluating
the performance on each level, we look for the largest reduction in FLOPs. The
corresponding level of complexity will be at the focus of our future efforts
for optimisation of performance. The results for these benchmarks are shown in
Section III-E and Figures 6 and 7.

####  II-E1 Matrix Multiplication

We perform the matrix multiplications for the square matrices of the same size
and record the time for the operation for the CPU backend Numpy and the GPU
backends PyTorch and CuPy. We use the built-in random() function of each
backend to randomly generate two square matrices of equal size as our input,
and we use the built-in matmul() function to produce the output matrix. The
size of the input matrices ranges from $10\times 10$ to $8192\times 8192$, and
the test is done repeatedly on four different data types: float, double,
complex64, and complex128. For the multiplication of two $n\times n$ matrices,
we define the number of complex operations to be $n^{3}$, and we calculate the
number of FLOPs for complex numbers as
$8\times\frac{\text{number\\_of\\_operations}}{\text{operation\\_time}}$.

####  II-E2 Tensor Network Contraction

We have two experiment groups in benchmarking the tensor contraction
performance: tensor contractions with a fixed contraction expression and
tensor contractions with many indices where each index has a small size. We
call the former group “tncontract fixed” because we fix the contraction
expression as “abcd,bcdf$\xrightarrow{}$acf,” and we call the latter one
“tncontract random” because we randomly generate the contraction expression.
In a general contraction expression, we sum over the indices not contained in
the result indices. In this fixed contraction expression, we sum over the
common index “b” and “d” and keep the rest in our result indices. We generate
two square input tensors of shape $n\times n\times n\times n$ and output a
tensor of shape $n\times n\times n$, where $n$ is a size ranging from 10 to
100. For the “tncontract random” group, we randomly generate the number of
contracted indices and the number of indices in the results first and then
fill in the shape array with size 2. For example, a contraction formula
“dacb,ad$\xrightarrow{}$bcd” (index “a” is contracted) needs two input
tensors: the first one with shape $2\times 2\times 2\times 2$ and the second
one with shape $2\times 2$. We use the formula
$2^{\text{number\\_of\\_different\\_indices}}$ to calculate the number of
operations, and we record the contraction time and compute the FLOPs value
based on the formula used in matrix multiplication. Following the same
procedure in matrix multiplication, we use the backends’ built-in functions to
randomly generate the input tensors based on the required size and the four
data types.

####  II-E3 Circuit Simulation

For numerical evaluations, we benchmark the full expectation value simulation
of the QAOA MaxCut problem for a 3-regular graph of size 30 and QAOA depth
$p=4$. We have two properties for evaluating the circuit simulation
performance: unmerged vs. merged backend and single vs. mixed backend.

##  III Results

Figure 1:  Breakdown of mean time to contract a single bucket by bucket width.
The test is performed for expectation value as described in III-A. CPU
backends are faster for buckets of width $\leq 13-16$, and GPU faster are
better for larger buckets. This picture also demonstrates that every
contraction operation spends some time on overhead which doesn’t depend on
bucket width, and actual calculation that scales exponentially with bucket
width.  Figure 2: Distribution of bucket width in the contraction of QAOA full
circuit simulation. The y-axis is log scale; 82% of buckets have width $\leq
6$, which have relatively large overhead time. Figure 3:  Breakdown of total
time spent on bucket of each size in full QAOA expectation value simulation.
The y-value on this plot is effectively one in Figure 1 multiplied by one in
Figure 2. This figure is very useful for analyzing the bottlenecks of the
simulation. It shows that most of the time for CPU backend is spent on large
buckets, but for GPU backends the large number of small buckets results in a
slowdown.

The experiment is performed on an NVIDIA DGX-2 server (provided by NVIDIA
corporation) with a V100 accelerator using the CUDA version 11.0. The baseline
NumPy backend is executed only on a CPU and labeled “einsum” in our experiment
since we use numpy.einsum() for the tensor computation. We also benchmark the
GPU library CuPy (on the GPU only) and PyTorch (on both the CPU and GPU).

###  III-A Single CPU-GPU Backends

We benchmark the performance of the full expectation value simulation of the
QAOA MaxCut problem on a 3-regular graph of size 30 with depth $p=4$, as shown
in in Figures 1, 2, and 3. This corresponds to contraction of 20 tensor
networks, one network per each lightcone. Our GPU implementation of the
simulator using PyTorch (labeled “torch_gpu”) achieves 70.3$\times$ speedup
over the CPU baseline and 1.92$\times$ speedup over CuPy.

Figure 1 shows the mean contraction time of various bucket widths in different
backends. In comparison with ”cupy” backend, the ”einsum” backend spends less
total time for bucket width less than 16, and the threshold value changes to
around 13 when being compared to ”torch_gpu” backend. Both GPU backends have
similar and better performance for larger bucket widths. However, this
threshold value can fluctuate when comparing the same pair of CPU and GPU
backends. This is likely due to the fact that the benchmarking platform are
under different usage loads.

Figure 3 provides a breakdown of the contraction times of buckets by bucket
width. This distribution is multimodal: A large portion of time is spent on
buckets of width 4. For CPU backends the bulk of the simulation time is spent
on contracting large buckets. Figure 2 shows the distribution of bucket
widths, where 82% of buckets have width less than 7. This signifies that
simulation has an overhead from contracting a large number of small buckets.

This situation is particularly noticeable when looking at the total
contraction time of different bucket widths. Figure 3 shows that the
distribution of time vs bucket width has two modes: for large buckets that
dominate the contraction time for CPU backends and for small buckets where
most of the time is spent on I/O and other code overhead.

###  III-B Merged Backend Results

Figure 4:  Breakdown of total contraction time by bucket width in full
expectation value simulation of problem size 30. Lines with the same color use
the same type of backends. The solid lines represent the merged version of
backends, and the dashed lines denote the baseline backends. The merged GPU
backends are better for buckets of width $\geq 20$.

The “merged” groups merge the indices before performing contractions. In Fig.
4, the three unmerged (baseline) backends are denoted by dashed lines, while
the merged backends are shown by solid lines. For the GPU backends CuPy and
PyTorch, the merged group performs significantly better for buckets of width
$\geq 20$. The CuPy merged backend always has a similar or better performance
compared with the CuPy unmerged group and has much better performance for
buckets of larger width. For buckets of width 28, the total operation time of
the unmerged GPU backends is about 0.28 seconds, compared with 32 ms
(8.75$\times$ speedup) for the CuPy merged group and 8.8 ms (31.82$\times$
speedup) for the PyTorch backend. But we do not observe much improvement for
the merged CPU backend.

###  III-C Mix CPU-GPU Backend Results

Backend Name | Device | Time (second) | Speedup
---|---|---|---
Torch_CPU | CPU | 347 | 0.71$\times$
NumPy (baseline) | CPU | 246 | 1.00$\times$
CuPy | GPU | 6.7 | 36.7$\times$
Torch_GPU | GPU | 3.5 | 70.3$\times$
Torch_CPU + Torch_GPU | Mixed | 2.6 | 94.8$\times$
NumPy + CuPy | Mixed | 2.1 | 117$\times$
TABLE I: Time for full QAOA expectation value simulation using backend that
utilize GPUs or CPUs. The expectation value is MaxCut on a 3-regular graph of
size 30 and QAOA depth $p=4$. Speedup shows the overall runtime improvement
compared with the baseline CPU backend “NumPy”. “Mixed” device means the
backend uses both CPU and GPU devices.

From Figure 1 one can see that GPU backends perform much better for buckets of
large width, while CPU backends are better for smaller buckets. We thus
implemented a mixed backend approach, which dynamically selects a device (CPU
or GPU) on which the bucket should be contracted. We select a threshold value
of 15 for the bucket width; any bucket that has a width larger than 15 will be
contracted on the GPU. Figure 3 shows that for GPU backends small buckets
occupy approximately 90% of the total simulation time. The results for this
approach are shown in Table I under backend names “Torch_CPU + Torch_GPU” and
“NumPy + CuPy.” Using a CPU backend in combination with Torch_GPU improves the
performance by 1.2$\times$, and for CuPy the improvement is 3$\times$. These
results suggest that using a combination of NumPy + Torch_GPU has the
potential to give the best results.

We have evaluated the GPU performance of tensor network contraction for the
energy calculation of QAOA. The problem is largely inhomogeneous with a lot of
small buckets and a few very large buckets. Most of the improvement comes from
using GPUs on large buckets, with up to 300$\times$ speed improvement. On the
other hand, the contraction of smaller tensors is faster on CPUs. In general,
if the maximum bucket width of a lightcone is less than $\sim 17$, the
improvement from using GPUs is marginal. In addition, large buckets require a
lot of memory. For example, a bucket of width 27 produces a tensor with 27
dimensions of size 2, and the memory requirement for complex128 data type is 2
GB. In practice, these calculations are feasible up to width $\sim 29$.

###  III-D Mixed Merged Backend Results

Since the performance of the NumPy-CuPy hybrid backend is the best among all
implemented hybrid backends, cross-testing between merged backends and hybrid
backends focuses on the combination of the NumPy backend and CuPy backend.
Because of the API constraint, the hybrid of a regular NumPy backend and a
merged CuPy backend was not implemented.

Backend Name | Device | Time (seconds) | Speedup
---|---|---|---
NumPy_Merged | CPU | 383 | 0.64$\times$
NumPy (baseline) | CPU | 246 | 1.00$\times$
CuPy | GPU | 6.7 | 36.7$\times$
CuPy_Merged | GPU | 5.6 | 43.9$\times$
NumPy + CuPy | Mixed | 2.1 | 117$\times$
NumPy_Merged + CuPy_Merged | Mixed | 1.4 | 176$\times$
TABLE II: Time for full QAOA expectation value simulation using different
Merged backends, as described in Section II-C. The expectation value is MaxCut
on a 3-regular graph of size 30 and QAOA depth $p=4$. Speedup shows the
overall runtime improvement compared with the baseline CPU backend “NumPy”.

In Table II, merging buckets provide a performance boost for the CuPy backend
and Numpy + CuPy hybrid backend but not the NumPy backend. CuPy_Merged is 20%
faster than CuPy, and NumPy_Merged + CuPy_Merged is 50% faster than its
regular counterpart. However, NumPy_Merged has an significant slowdown
compared with the baseline NumPy, suggesting that combining the regular NumPy
backend with the merged CuPy backend can provide more speedup for the future.

Figure 5:  Breakdown of sum contraction time by bucket width for merged
backends. CPU backends are better for buckets of width $\leq 15$, and GPU
backends are better for larger buckets. The hybrid backend’s GPU backend
spends outperforms the regular GPU backend for buckets of width $\geq 15$.

In Fig. 5, CPU performance is better than GPU performance when the bucket
width is approximately less than 15. After 15, GPU performance scales with
width much better than that of CPU performance, providing a significant speed
boost over the CPU in the end. GPU performance of the hybrid backend is better
than that of pure GPU backend for buckets of width $\geq 15$. This speedup of
the hybrid backend is likely caused by less garbage handling for the GPU since
most buckets aren’t stored on GPU memory.

###  III-E Synthetic Benchmarks

We also benchmark the time required for the basic operations: matrix
multiplication, tensor network contraction with fixed contraction indices, and
tensor network contraction with random indices, as well as circuit
contractions.

The summary of the results is shown in Table III, which compares FLOPs count
for similar-sized problems of different types. Figures 6 and 7 show dependence
of FLOPs vs problem size for different problems. We observe 80% of theoretical
peak performance on GPU for matrix multiplication. Switching to pairwise
tensor network contraction shows similar FLOPs for GPU, while for CPU, it
results in 10$\times$ FLOPs decrease. A significant reduction in performance
comes from switching from pairwise tensor network contractions of a tensor
with few dimensions of large size to tensors with many permuted dimensions and
small size. This reduction in performance is about 10$\times$ for both CPU and
GPU. This observation suggests that further improvement can be achieved by
reformulating the tensor network operations in a smaller tensor by transposing
and merging the dimensions of participating tensors. It is partially addressed
in using the merged indices approach, where the contraction dimension is
increased. The “Bucket Contraction Merged” task shows 45% of theoretical peak
performance, which significantly improves compared to the unmerged
counterpart.

The significant reduction of performance comes when we compare bucket
contraction and full circuit contraction. It was explained in detail in
Section III-A and is caused by overhead from small buckets. It is evident from
Figure 3 that most of the time in GPU simulation is spent on overhead from
small bucket contraction. This issue is addressed by implementing the mixed
backend approach.

It is also notable that the merged approach does not improve the performance
for CPU backends which is probably due to an inefficient implementation of
original numpy.einsum().

Task | CPU FLOPs | GPU FLOPs
---|---|---
Matrix Multiplication | 50.1G | 2.38T
Tensor Network Fixed Contraction | 5.53G | 1.36T
Tensor Network Random Contraction | 640M | 97.5G
Bucket Contraction Unmerged | 241M | 61.9G
Bucket Contraction Merged | 542M | 1.14T
Lightcone Contraction Unmerged | 326M | 4.92G
Lightcone Contraction Merged | 177M | 3.1G
Circuit Contraction Mixed | 30.7G
TABLE III: Summary of GPU and CPU FLOPs for different tasks at around 100
million operations. Matrix Multiplication and Tensor Contraction tasks are
described in Section III-E. “Bucket Contraction” groups record the maximum
number of FLOPs for a single bucket. “Lightcone Contraction” groups contain
the FLOPs data on a single lightcone where the sum of operations is
approximately 100 millions, small and large buckets combined.

####  III-E1 Matrix Multiplication

The multiplication of square matrices of size 465 needs approximately 100
million complex operations according to our calculation of operations value.
The average operation time for the multiplication of two randomly generated
complex128 square matrices of size 465 is 0.3 ms on the GPU, which achieves
50$\times$ speedup compared with the operation time of 16 ms on the CPU; NumPy
produces 50G FLOPs on CPU, and the GPU backend CuPy reaches 2.38T FLOPs for
this operation. We observe that the CPU backend has an advantage in performing
small operations: for matrices of size $10\times 10$, the CPU backend NumPy
spends only 5.8 µs for the multiplication, while the best GPU backend PyTorch
spends 27 µs on the operation. When the matrix size is less than $2000\times
2000$ for the GPU backends, PyTorch outperforms CuPy, and CuPy is slightly
better for much larger operations. Moreover, the operation time for both CPU
and GPU backends decreases slightly when the size of matrices increases from
1000 to 1024 and from 4090 to 4096.

####  III-E2 Fixed Tensor Network Contraction

We use the fixed contraction formula “abcd,bcdf$\xrightarrow{}$acf” and
control the size of the tensor indices from 10 to 100. Even for the smallest
case when the number of operations is 100,000 with indices of size 10, the
slowest GPU backend is faster than the CPU backend Numpy, which spends 0.3 ms
on the contraction. For the GPU backends, we achieve 1.36T FLOPs for this
fixed contraction, which is 57% of the recorded peak performance. In
accordance with the matrix multiplication results, the CuPy backend performs
better than the PyTorch backend in the fixed tensor network contractions only
when the number of operations is greater than 1G.

####  III-E3 Random Tensor Network Contraction

We let the number of indices be any number between 4 and 25, and we set the
size of each mode to be 2. For example, we have 5 indices in total, and we
randomly generate a contraction sequence “caedb,eab$\xrightarrow{}$cde,” so
the sizes of the input tensors are $2^{5}$ and $2^{3}$, resulting in an output
tensor of size $2^{3}$. We reach 97.5 G FLOPs for the GPU backends and 640 M
for the CPU backend only when performing this random contraction. As shown in
Fig. 6, the mean FLOPs drop significantly when we use random contraction (in
green) instead of fixed contraction (in red) on the CuPy backend. On the CPU,
the gap increases with the increasing number of operations according to Fig.
7. Therefore, contractions on tensors with small numbers of indices of large
size have better performance than contractions on tensors with many indices of
small size. The ”tncontract random” group is designed to break down the
circuit simulation to tensor contraction operations, so it overlaps with the
results from the ”bucket unmerged” group in Fig. 6. From the difference in
performance of the random and the fixed tensor contraction group, we design
the merged bucket group to improve the performance of contractions. Our goal
is to make the bucket simulation curve close to the tensor contraction fixed
group (the red curve).

Figure 6:  FLOPs vs. the number of operations for all tasks on the CuPy
backend. “circuit unmerged” and “circuit merged” are results of expectation
value of the full circuit simulation of QAOA MaxCut problem on a 3-regular
graph of size 30 with depth $p=4$. “tncontract random” tests on tensors of
many indices where each index has a small size. “tncontract fixed” uses the
contraction sequence “abcd,bcdf$\xrightarrow{}$acf” for all contractions.
“matmul” performs matrix multiplication on square matrices. All groups use
complex128 tensors in the operation. We use the triangles to denote the data
at $\sim 100$ million operations, which is shown in Table III.  Figure 7:
FLOPs vs. the number of operations for all tasks on NumPy backend. Same
problem setting as Fig. 6. “tncontract random” outperforms “tncontract fixed”
as the ops value increases. Merged backend does not have an advantage on CPU
compared to the unmerged backend. We use the triangles to denote the data at
$\sim 100$ million operations, which is shown in Table III.

##  IV Conclusions

This work has demonstrated that GPUs can significantly speed up quantum
circuit simulations using tensor network contractions. We demonstrate that
GPUs are best for contracting large tensors, while CPUs are slightly better
for small tensors. Moving the computation onto GPUs can dramatically speed up
the computation. We propose to use a contraction backend that dynamically
assigns the CPU or GPU device to tensors based on their size. This mixed
backend approach demonstrated a 176$\times$ improvement in time to solution.

We observe up to 300$\times$ speedup on GPU compared to CPU for individual
large buckets. In general, if the maximum bucketwidth of a lightcone is less
than $\sim 17$, the improvement from using GPUs is marginal. It underlines the
importance of using a mixed CPU/GPU backend for tensor contraction and using
device selection for the tensor at runtime to achieve the maximum performance.
On NVIDIA DGX-2 server we found out that the threshold is $\sim 15$, but it
may change for other computing systems.

We also demonstrated the performance of the merged indices approach, which
improves the arithmetic intensity and provides a significant FLOP improvement.
Our synthetic benchmarks for various tensor contraction tasks suggest that
additional improvement can be obtained by transposing and reshaping tensors in
pairwise contractions.

The main conclusion of this paper is that we found that GPUs can dramatically
increase the speed of tensor contractions for large tensors. The smaller
tensors need to be computed on a CPU only because of overhead to move on and
off data to a GPU. We show that the approach of merged indices allows to speed
up large tensors contraction, but it does not solve the problem completely.
Where to compute tensors leads to the problem of optimal load balancing
between CPU and GPU. This potential issue will be the subject of our future
work, as well as testing of the performance of the code on new NVidia DGX
systems and GPU supercomputers using cuTensor and cuQuantum software packages
developed by NVidia.

## Acknowledgments

Danylo Lykov and Yuri Alexeev are supported by the Defense Advanced Research
Projects Agency (DARPA) grant. Yuri Alexeev and Angela Chen are also supported
in a part by the Exascale Computing Project (17-SC-20-SC), a joint project of
the U.S. Department of Energy’s Office of Science and National Nuclear
Security Administration, responsible for delivering a capable exascale
ecosystem, including software, applications, and hardware technology, to
support the nation’s exascale computing imperative. Huaxuan Chen is supported
in part by the U.S. Department of Energy, Office of Science, Office of
Workforce Development for Teachers and Scientists (WDTS) under the Science
Undergraduate Laboratory Internships Program (SULI). This work used the
resources of the Argonne Leadership Computing Facility, which is DOE Office of
Science User Facility supported under Contract DE-AC02-06CH11357.

## References

  * [1] Y. Alexeev, D. Bacon, K. R. Brown, R. Calderbank, L. D. Carr, F. T. Chong, B. DeMarco, D. Englund, E. Farhi, B. Fefferman, A. Gorshkov, A. Houck, J. Kim, S. Kimmel, M. Lange, S. Lloyd, M. Lukin, D. Maslov, P. Maunz, C. Monroe, J. Preskill, M. Roetteler, M. Savage, and J. Thompson, “Quantum computer systems for scientific discovery,” _PRX Quantum_ , vol. 2, no. 1, p. 017001, 2021.
  * [2] K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, T. Lippert, H. Watanabe, and N. Ito, “Massively parallel quantum computer simulator,” _Computer Physics Communications_ , vol. 176, no. 2, pp. 121–136, 2007.
  * [3] M. Smelyanskiy, N. P. Sawaya, and A. Aspuru-Guzik, “qHiPSTER: the quantum high performance software testing environment,” _arXiv preprint arXiv:1601.07195_ , 2016.
  * [4] T. Häner and D. S. Steiger, “0.5 petabyte simulation of a 45-qubit quantum circuit,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_. ACM, 2017, p. 33.
  * [5] X.-C. Wu, S. Di, E. M. Dasgupta, F. Cappello, H. Finkel, Y. Alexeev, and F. T. Chong, “Full-state quantum circuit simulation by using data compression,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE Computer Society, 2019, pp. 1–24.
  * [6] E. Bernstein and U. Vazirani, “Quantum complexity theory,” _SIAM Journal on Computing_ , vol. 26, no. 5, pp. 1411–1473, 1997.
  * [7] (2020) QuaC (quantum in c) is a parallel time dependent open quantum systems solver. [Online]. Available: https://github.com/0tt3r/QuaC
  * [8] I. L. Markov and Y. Shi, “Simulating quantum computation by contracting tensor networks,” _SIAM Journal on Computing_ , vol. 38, no. 3, pp. 963–981, 2008\.
  * [9] E. Pednault, J. A. Gunnels, G. Nannicini, L. Horesh, T. Magerlein, E. Solomonik, and R. Wisnieff, “Breaking the 49-qubit barrier in the simulation of quantum circuits,” _arXiv preprint arXiv:1710.05867_ , 2017\.
  * [10] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, and H. Neven, “Simulation of low-depth quantum circuits as complex undirected graphical models,” _arXiv preprint arXiv:1712.05384_ , 2017.
  * [11] E. Farhi and A. W. Harrow, “Quantum supremacy through the quantum approximate optimization algorithm,” _arXiv preprint arXiv:1602.07674_ , 2016.
  * [12] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” _arXiv preprint arXiv:1411.4028_ , 2014.
  * [13] D. J. Egger, J. Mareček, and S. Woerner, “Warm-starting quantum optimization,” _Quantum_ , vol. 5, p. 479, 2021.
  * [14] R. Shaydulin, I. Safro, and J. Larson, “Multistart methods for quantum approximate optimization,” in _2019 IEEE High Performance Extreme Computing Conference (HPEC)_. IEEE, 2019, pp. 1–8.
  * [15] R. Shaydulin, H. Ushijima-Mwesigwa, C. F. Negre, I. Safro, S. M. Mniszewski, and Y. Alexeev, “A hybrid approach for solving optimization problems on small quantum computers,” _Computer_ , vol. 52, no. 6, pp. 18–26, 2019.
  * [16] R. Shaydulin, S. Hadfield, T. Hogg, and I. Safro, “Classical symmetries and QAOA,” _arXiv preprint arXiv:2012.04713_ , 2020.
  * [17] S. Khairy, R. Shaydulin, L. Cincio, Y. Alexeev, and P. Balaprakash, “Learning to optimize variational quantum circuits to solve combinatorial problems,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 03, 2020, pp. 2367–2375.
  * [18] R. Schutski, D. Lykov, and I. Oseledets, “Adaptive algorithm for quantum circuit simulation,” _Physical Review A_ , vol. 101, no. 4, p. 042335, 2020\.
  * [19] D. Lykov, R. Schutski, A. Galda, V. Vinokur, and Y. Alexeev, “Tensor network quantum simulator with step-dependent parallelization,” _arXiv preprint arXiv:2012.02430_ , 2020.
  * [20] D. Lykov and Y. Alexeev, “Importance of diagonal gates in tensor network simulations,” _arXiv preprint 2106.15740_ , 2021.
  * [21] D. Lykov, R. Schutski, A. Galda, V. Vinokur, and Y. Alexeev, “Tensor network quantum simulator with step-dependent parallelization,” _arXiv preprint arXiv:2012.02430_ , 2020.
  * [22] R. Schutski, D. Lykov, and I. Oseledets, “Adaptive algorithm for quantum circuit simulation,” _Phys. Rev. A_ , vol. 101, p. 042335, Apr 2020.