# Two-dimensional Dirac semiconductor and its material realization

Botao Fu  College of Physics and Electronic Engineering, Center for
Computational Sciences, Sichuan Normal University, Chengdu, 610068, China
Chao He  College of Physics and Electronic Engineering, Center for
Computational Sciences, Sichuan Normal University, Chengdu, 610068, China  Da-
Shuai Ma  Key Lab of advanced optoelectronic quantum architecture and
measurement (MOE), Beijing Key Lab of Nanophotonics $\&$ Ultrafine
Optoelectronic Systems, and School of Physics, Beijing Institute of
Technology, Beijing 100081, China  Zhi-Ming Yu  Key Lab of advanced
optoelectronic quantum architecture and measurement (MOE), Beijing Key Lab of
Nanophotonics $\&$ Ultrafine Optoelectronic Systems, and School of Physics,
Beijing Institute of Technology, Beijing 100081, China  Yong-Hong Zhao
College of Physics and Electronic Engineering, Center for Computational
Sciences, Sichuan Normal University, Chengdu, 610068, China  Yugui Yao
<EMAIL_ADDRESS>Key Lab of advanced optoelectronic quantum architecture and
measurement (MOE), Beijing Key Lab of Nanophotonics $\&$ Ultrafine
Optoelectronic Systems, and School of Physics, Beijing Institute of
Technology, Beijing 100081, China

###### Abstract

We propose a new concept of two-dimensional (2D) Dirac semiconductor which is
characterized by the emergence of fourfold degenerate band crossings near the
band edge and provide a generic approach to realize this novel semiconductor
in the community of material science. Based on the first-principle
calculations and symmetry analysis, we discover recently synthesised triple-
layer (TL) BiOS2 is such Dirac semiconductor that features Dirac cone at X/Y
point, protected by nonsymmorphic symmetry. Due to sandwich-like structure,
each Dirac fermion in TL-BiOS2 can be regarded as a combination of two Weyl
fermions with opposite chiralities, degenerate in momentum-energy space but
separated in real space. Such Dirac semiconductor carries layer-dependent
helical spin textures that never been reported before. Moreover, novel
topological phase transitions are flexibly achieved in TL-BiOS2: (i) an
vertical electric field can drive it into Weyl semiconductor with switchable
spin polarization direction, (ii) an extensive strain is able to generate
ferroelectric polarization and actuate it into Weyl nodal ring around X point
and into another type of four-fold degenerate point at Y point. Our work
extends the Dirac fermion into semiconductor systems and provides a promising
avenue to integrate spintronics and optoelectronics in topological materials.

Introduction.— Dirac fermion with linearly band dispersion was first
discovered in 2D grapheneNeto  _et al._ (2009). Whereafter, it was rapidly
popularized to 3D, and abundant topological semimetal phases including Weyl,
Dirac and nodal-line semimetals are harvested in plenty of materialsYan and
Felser (2017); Fang  _et al._ (2016); Weng  _et al._ (2016); Gao  _et al._
(2019); Armitage  _et al._ (2018). The unique electronic structures of
topological semimetals lead to protected surface states and novel response to
external fields and thus attract intensive research interestsWang and Wang
(2018); Nagaosa  _et al._ (2020). Actually, these quasi-fermions can be
generalized towards superconductor and metallic systemsWang  _et al._ (2016);
Chiu and Schnyder (2014). For instance, the nodal line fermion was observed on
superconducting material PbTaSe2Bian  _et al._ (2016), and the flat surface
sates of nodal-line was reported widely exist in alkaline-earth metals and
noble metalsLi  _et al._ (2016); Yan  _et al._ (2015).

In comparison with metal and semimetal, the semiconductor materials are
particularly suitable for Dirac devices due to its high tunability and
compatibility with modern electronic industry. Therefore, the introduction of
Weyl/Dirac fermion into semiconductors to develop Weyl/Dirac semiconductors
could create a new degree of freedom for the future design of semiconductor
electronic and optoelectronic devices. Recently, the chiral Weyl fermion is
theoretically predictedHirayama  _et al._ (2015) and experimentally observed
in elemental tellurium, which is an intrinsic semiconductor with band gap
about 0.36 eV, thus was dubbed as Weyl semiconductorSakano  _et al._ (2020);
Ideue  _et al._ (2019); Zhang  _et al._ (2019). Based on these considerations,
in this paper we will generalize the four-fold degenerating Dirac fermion from
semimetal to semiconductor in two-dimension. Because 2D materials possess
prior mechanical properties and small size that are more favorable for
integration and regulation.

Figure 1: A nonsymmorphic operator $\widetilde{g}=\\{g|\bm{t}\\}$ leads to
band crossings at $\mathbf{G}/2$ ($\mathbf{G}$ is the reciprocal lattice
vector) for a system with $\mathcal{PT}$ symmetry in the present of SOC. (a)
The electronic structure of system with electron filling number is $4n+2$. (b)
The electronic structure of system with electron filling number is $4n$. The
dash line indicates the Fermi level. (c) Schematic illustration for
construction of Dirac fermion in multiple-layer structure. Two Weyl fermion
locate on individual monolayers merge together to form a Dirac fermion.

Proposal of 2D Dirac semiconductors— Firstly, we review the concept of 2D
Dirac semimetal in the presence of spin-orbital coupling (SOC), as proposed by
S. M. Young at alYoung and Kane (2015). For a system with both time reverse
$\mathcal{T}$ and space inversion $\mathcal{P}$ in addition with a
nonsmmorphic symmetry $\widetilde{g}=\\{g|\bm{t}\\}$, where $g$ is point group
operator and $\bm{t}$ is fractional translation. As shown in Fig. 1(a), along
a $g$-invariant line, e.g. $\mathbf{0}$ to $\mathbf{G}$, each double-
degenerated band can be marked by the opposite and ${\bm{k}}$-dependent
eigenvalue of $\widetilde{g}$: $\pm ie^{i\bm{kt}}$. It equals to $\pm i$ at
$\mathbf{0}$ and gradually evolves into $\pm 1$ at $\mathbf{G}/2$. At
$\mathbf{G}/2$, the $\mathcal{T}^{2}=-1$ guarantees two $+1$ states or $-1$
states to be degenerate, namely Kramers degeneracy. Consequently, two branches
of double-degenerate bands have to stick together and form a four-fold
degenerate Dirac point at $\mathbf{G}/2$. If the electron filling number
satisfies $4n+2(n\geq 0)$, the Fermi level then crosses Dirac point and it
becomes an ideal Dirac semimetal with point-like Fermi surface.

Despite concise and intriguing picture, the practical and ideal 2D Dirac
semimetal materials are very rareGuan  _et al._ (2017); Kowalczyk  _et al._
(2020). One reason may be the limitation of electron filling number. By
counting the 2D nonsmmorphic materials in material databaseMounet  _et al._
(2018), we find only 16$\%$ of them have $4n+2$ electron filling number while
the rest host $4n$ electrons. The other reason may be that the second-order
SOC effect is relatively small compared with second-order electron hopping
which results in the emergence of undesirable trivial metal states at Fermi
level. This can be explicitly verified through a simple tight-binding modelSup
. From other perspective, as in Fig. 1(b) for most materials with $4n$
electrons, the Fermi level lies inside the gap between two sets of bands
including Dirac fermions at $\mathbf{G}/2$. In general, the system is a
semiconductor. Suppose the position of band edge happens to reside at or
nearby $\mathbf{G}/2$, it will become so-called “Dirac semiconductor” with
symmetry-enforced Dirac cone close to Fermi level. In such Dirac
semiconductor, the signature of Dirac fermions could be probed by angle-
resolved photoemission spectroscopy (ARPES) and transport experiments as what
have been done in Weyl semiconductorsKowalczyk  _et al._ (2020); Zhang  _et
al._ (2019); Qiu  _et al._ (2020).

In Fig. 1(c), we introduce an practical approach to construct this Dirac
semiconductor in multiple-layer systems. Suppose a monolayer system without
$\mathcal{P}$ but with $\mathcal{T}$ symmetry, with SOC being taken into
consideration each band becomes spin splitting except for that at the time-
reversal invariant momenta (TRIM) where the Kramers degeneracy occurs. It’s
worth noting that such Kramers degeneracy is recently reported host
unanticipated topological charge in 3D chiral lattice, noted as Kramers-Weyl
fermionChang  _et al._ (2018); Li  _et al._ (2019). Bearing this in mind, we
stack up two $\mathcal{P}$-breaking monolayers which are semiconductor with
band edge at $\mathbf{G}/2$, and meanwhile impose an in-plane fractional
translation. The constructed bilayer system hosts both $\mathcal{P}$ and glide
mirror symmetries. Consequently, two Kramers-Weyl points from two monolayers
will merge one Dirac point at $\mathbf{G}/2$. This tells us that Dirac
semiconductor can exist in certain stacked bilayer or some intrinsic bilayer
or multiple layer materials. More importantly, with unique multiple-layer
structure and semiconducting nature, the Dirac semiconductor will exhibit
flexibly tunable electronic and topological properties under external fields.

Material Realization— Based on above analysis, we find out four kinds of
multiple-layer semiconductorsSup  which host nonsmmorphyic space group and
proper band edge position. Here, taking TL-BiOS2 as a prototype, we will
demonstrate the emergence of Dirac fermion around the band edge and its
manipulation by external fields. The TL-BiOS2 is predicted to be a stable
direct band-gap semiconductor with high carrier mobility and large visible
light absorptionZhang  _et al._ (2018). In experiment, the high carrier
mobility as well as power conversion effciency as solar cell has been recently
observed in 2D BiOS2Huang and Yu (2020). As shown in Fig. 2(a), the TL-BiOS2
consists of one BiO layer in the middle and two BiS2 layers on the top and
bottom, forming sandwich-like structure. The optimized lattice constance is
$|\bm{a}|$=$|\bm{b}|$=3.95 Å. Although, each BiS2 layer is inversion-
asymmetric with non-centrosymmetric polar site point group $C_{4v}$. Two BiS2
layers become each other’s inversion partners and meanwhile have inter-layer
relative translation, $(\bm{a}+\bm{b})/2$. Hence the whole system hosts
centrosymmetric and nonsymmorphic space group $P4/nmm$ (No.129).

Figure 2: (a) The atomic structure of TL-Bi2OS2. (b), (c)The real-space charge
density distribution of the maxima of HVB1 and HVB2. (d) The orbital-resolved
band structure without SOC effect. The green color represents Bi’s $p$-orbital
and the red color represents S’s $p$-orbital on BiS2 layer, while the blue
color represents all orbitals from BiO layer. (e) The band structure with SOC
effect. (f) The enlarged figures around LCB and HVB with SOC. (g) The 3D Dirac
cone around $X$ with helical spin textures indicated by arrows. The different
color represent bands from different BiS2 layers. The Dirac fermions are
marked by yellow dots.

The electronic structures of TL-Bi2OS2 under generalized gradient
approximation (GGA) are displayed in Fig. 2(d)-(e). In the absence of SOC, it
has a direct band gap ($E_{g}$=0.97 eV) with both valence band maximum (VBM)
and conduction band minimum (CBM) located at X/Y point. The quivalentence of X
and Y is guaranteed by the $C_{4z}$ symmetry. Importantly, both the highest
valence band (HVB) and lowest conduction band (LCB) are double-degenerate not
only at X/Y point but along X-M-Y path. This degeneracy is protected by screw
axis $\widetilde{C}_{2x}$/$\widetilde{C}_{2y}$Fan  _et al._ (2018), in which
the tilde refers to additional $(\bm{a}+\bm{b})/2$ translation. The orbital-
resolved band structure demonstrates that both HVB and LCB primarily derive
from BiS2 layers, but the former from S-$p$ orbital while the latter from
Bi-$p$ orbital. The charge density distributions of double-degenerate HVB
(HVB1 and HVB2) states are calculated in Fig. 2(b)-(c), from which we find
despite their degeneracy in momentum-energy space, the HVB1 and HVB2 locate
separately on top and bottom BiS2 layers, respectively. Similar phenomenon
also happens for the LCB. Thus, TL-Bi2OS2 well satisfies the conditions for
the emergence of Dirac semiconductor: (i) Two $\mathcal{P}$-breaking BiS2
layers stack up via electrostatic interaction with intermediate BiO2 layer
which simultaneously host fractional translation and inversion symmetries.
(ii) The double-degenerate VBM or CBM at X/Y in TL-Bi2OS2 originate from
opposite BiS2 layers.

When considering SOC effect in Fig. 2(e), each band becomes double degenerate
because of $(\mathcal{PT})^{2}=-1$. The band gap ($E_{g}$) at X/Y reduces from
0.97 eV to 0.40 eV. More interesting, the LCB and HVB split along XM path and
exhibits Rashba-like dispersions as demonstrated in Fig. 2(f). The Rashba
parameter can be defined as $\alpha_{R}=2E_{R}/k_{R}$, which is about 1.27 eVÅ
and 2.15 eVÅ for HCB and LVB, respectively. In reality, it’s distinctive from
traditional Rashba splitting which occurs in non-centrosymmetric
materialsIshizaka  _et al._ (2011); Maaß  _et al._ (2016); Zhang and Liu
(2019); Di Sante  _et al._ (2013). Around X/Y point, two double-degenerated
bands cross each other to form a fourfold-degenerated Dirac fermion as
expected. This could be understood from following perspective. Since the HVB1
and HVB2 locate separately on opposite BiS2 layers without $\mathcal{P}$, the
SOC effect certainly gives a Rashba splitting for them. Then considering two
BiS2 layers in TL-Bi2OS2 as a whole, the nonsmmorphic operation in combination
with space inversion will stick two sets of Rashba bands together, and in
particular enables the appearance of Dirac fermion at X/Y. Therefore, although
this Dirac semiconductor host $\mathcal{PT}$ symmetry, one can identify the
spin textures of Dirac fermion. As shown in Fig. 2(g), it indeed hosts
opposite chiral spin textures on different BiS2 layers, which may provide
potential applications in spintronics.

Figure 3:  (a) The electronic structure under external vertical electric
filed, E${}_{\text{ext}}$=0.3 V/Å. (b) The enlarged view around valence band.
(c) The relation of $\Delta E$ and external electric filed. (d) The
distribution of Berry curvature around two Weyl points at X, and the Berry
phases circled them are given. (e) The electronic structure under $2\%$
biaxial strain, with spontaneous electric polarization $P_{s}$ along $y$. (f)
The enlarged view around the valence band at X. The eigenvalues of
$\widetilde{M}_{z}$ and $\widetilde{C}_{2y}$ are given. (g) The 3D view of
bands at (f), the nodal line is labelled by black line. (h) The enlarged views
around Y along different $k-$paths. The red and green dot indicates the
different Weyl fermions at X in (b) and (f), the Yellow dot represents Dirac
fermion in (h).

To capture the physics of Dirac fermion around band edge, we are going to
build an effective $k{\cdot}p$ model. For the VBM at X point, we can choose
the degenerate HVB1 and HBV2 states as the basis in the absence of SOC. Then
considering SOC effect the Hamiltonian around X is written as,

$H=\left(\begin{array}[]{cc}h_{+}&T\\\ T&h_{-}\\\ \end{array}\right),$ (1)
$h_{\pm}=\epsilon_{k}\pm(\alpha_{y}k_{y}\sigma_{x}-\alpha_{x}k_{x}\sigma_{y}),$
(2)
$\epsilon_{k}=\frac{\hbar^{2}k_{x}^{2}}{2m_{x}^{*}}+\frac{\hbar^{2}k_{y}^{2}}{2m_{y}^{*}},$
(3)

where $h_{\pm}$ describes the electron from HVB1/HVB2 that locates on
top/bottom BiS2 layer and $T$ is the inter-layer coupling term. The
$\alpha_{x,y}$ and $m_{x,y}^{*}$ are the anisotropic Rashba parameters and
effective mass along $x/y$ directions, respectively. The $\sigma_{x,y,z}$ is
Pauli matrix that represents spin degree of freedom. In Eq. 2 we consider the
predominant Rashba SOC term that originated from the interfacial polar filed
between BiO and BiS2 layer. Since the top and bottom BiS2 layers feel opposite
polar filed, they possess opposite Rashba spin texture as indicated on
$h_{\pm}$. At X point the little group is $D_{2h}$, the time reverse and space
inversion are chosen as $\mathcal{T}=i\sigma_{y}K\tau_{z}$ and
$\mathcal{P}=\sigma_{0}\tau_{y}$, then based on specific commutation relations
other symmetric operations can be readily written as
$\widetilde{C}_{2x}=\sigma_{x}\tau_{x}$,
$\widetilde{C}_{2y}=i\sigma_{y}\tau_{y}$ and
$\widetilde{M}_{z}=\sigma_{z}\tau_{x}$, where the $\tau_{x,y,z}$ is Pauli
matrix that represents layer degree of freedom. With these constraints the
inter-layer term is given as $T=t_{i}k_{x}\sigma_{0}\tau_{x}$, $t_{i}$ is
inter-layer coupling coefficient. Hereafter, the low-energy effective
Hamiltonian at X with first-order approximation is obtained as

${H_{\text{X}}=(\alpha_{y}k_{y}\sigma_{x}-\alpha_{x}k_{x}\sigma_{y})\tau_{z}+t_{i}k_{x}\sigma_{0}\tau_{x}}.$
(4)

With a straightforward solving, we get two branches of doubly degenerate
energy spectra,

${E_{\pm}^{\text{X}}=\pm\sqrt{(\alpha_{x}^{2}+t_{i}^{2})k_{x}^{2}+\alpha_{y}^{2}k_{y}^{2}}},$
(5)

which indeed depicts a quadruple degenerated Dirac fermion at
$k_{x}$=$k_{y}$=0 with Fermi velocities
$v_{x}=\sqrt{\alpha_{x}^{2}+t_{i}^{2}}$ and $v_{y}=|a_{y}|$. By fitting the
$k\cdot p$ model with our first-principle calculations, we obtain
$\alpha_{x}$=1.40 eVÅ, $\alpha_{y}$=2.15 eVÅ, $t_{i}$=0.09 eVÅ. The $t_{i}$ is
much smaller than $\alpha_{x,y}$ indicates the interlayer interaction is weak
and can be reasonably neglected. Then the Dirac fermion in Eq. 4 is able to be
decomposed into two Weyl fermions with opposite layer index,

$H_{\pm}^{W}=\pm(\alpha_{y}k_{y}\sigma_{x}-\alpha_{x}k_{x}\sigma_{y}).$ (6)

It indicates the Dirac fermion can be approximately taken as two spatially
isolated Weyl fermions with opposite layer-index and carrying opposite helical
spin texture, which differs from other known Dirac fermion in semimetal or
metal systems nd may extend the potential application of Dirac fermion. For
instance, the linearly dispersion as well as in-plane helical spin texture in
Dirac semiconductor may suppress direct backscattering and facilitate carrier
mobility, and the layer-resolved spin polarizationZhang  _et al._ (2014); Liu
_et al._ (2015a, 2016); Guan and Luo (2020) can be flexibly controlled under
external filed. Moreover, we find for two and three layers of TL-BiOS2, the
Dirac fermion still existsSup , making it easy to probe in experiment. When
extending to 3D bulk limitation, the Dirac fermion transforms into a Dirac
nodal line along $k_{z}$ direction, where each point can be taken as a Dirac
point, thus dubbed as 3D Dirac nodal-line semiconductor. We believe these
unique features in such Dirac semiconductors deserve further experimental
explorations.

Manipulation of TPTs.—

The topological phase transition (TPT) manipulated by external electric field
is of great significance, especially in lower dimensionLi and Chang (2009); Li
_et al._ (2017); Liu  _et al._ (2015b). For example, the electric-field-tuned
TPT is recently reported in ultra-thin Na3BiCollins  _et al._ (2020), a well-
known Dirac semimetal. Since Dirac fermion is protected by $\mathcal{PT}$ in
addition with other crystal symmetries, breaking either $\mathcal{P}$ or
$\mathcal{T}$ will induce various TPTsYoung and Kane (2015). The fact that the
Dirac fermion in TL-BiOS2 contains two spatially isolated Weyl points provides
a good opportunity to manipulating by vertical electric field. Also, the
semiconducting nature of TL-BiOS2 is in favor of electric field control than
semimetal or metal systems in which the electrostatic screening predominates.
Specifically, an vertical electric field breaks inversion symmetry of TL-BiOS2
by inducing stacking potential on different BiS2 layers. This process can be
approximately depicted by adding a mass term $V_{1}=V_{z}\sigma_{0}\tau_{z}$
to $H_{\text{X}}$. Neglecting interlayer coupling term, $H_{\text{X}}+V_{1}$
resolves into ${V_{z}\sigma_{0}+H_{+}^{W}}$ and
${-V_{z}\sigma_{0}+H_{-}^{W}}$, which indicate a Dirac fermion can be
successful split into two Weyl points with energy difference $\Delta
E=2V_{z}$. In Fig. 3 (a)-(c), we perform first-principle calculations with
external vertical electric filed (E${}_{\text{ext}}$) varying from 0.0 to 0.4
V/Å. We find the four-fold degenerated Dirac fermions at HVB and LCB at X/Y
split into a pair of Karamer-Weyl fermions. The split amplitude $\Delta E$ is
linearly proportional to the size of E${}_{\text{ext}}$ with a ratio of 0.96
eÅ as shown in Fig. 3(c). We already know these two Weyl fermions have layer-
dependent chiral spin textures. Moreover, we show they also demonstrate large
Berry curvature distributions and can be identified by the quantized $\pm\pi$
Berry phase surrounding them as displayed in Fig. 3(d). Hence by switching
Electric field direction, the spin polarization direction and spin Hall
coefficient can be modified, working as spin filed effect transistor (SFET)Liu
_et al._ (2013).

Strain is another effective method for engineering the electronic, magnetic
and topological properties of 2D materialsXu  _et al._ (2017); Zhao  _et al._
(2020). As we know the Dirac point in TL-Bi2OS2 is guaranteed by unique
nonsymmorphic space group, so called symmetry-enforced that robust against any
symmetry-remained perturbations. In principle, a simple uniaxial or biaxial
strain on TL-Bi2OS2 won’t change the nature of nonsymmorphic. Luckily, we find
TL-Bi2OS2 is a novel piezoelectric materialCui  _et al._ (2018); Guan  _et
al._ (2020). A very tiny biaxial extensive strain will induce a Peierls-like
atomic distortion on BiS2 layers along $y$ direction that breaks space
inversion and gives switchable in-plane electronic polarizationSup . In Fig.
3(e), we calculate the bandstructure under $2\%$ biaxial strain. One can find
the band gap at X point (0.85 eV) becomes larger than that at Y (0.70 eV).
Because the in-plane polar distortion breaks the $C_{4z}$ rotation and lifts
the inter-valley degeneracy between X and Y.

The space group of ferroelectric distorted TL-BiOS2 is $Pmn2_{1}$, which
contains four symmetry operations: $M_{x}$, $\widetilde{M}_{z}$,
$\widetilde{C}_{2y}$ and identity. At X point, the lacking of $\mathcal{P}$
and $\widetilde{C}_{2x}$ permits the emergence of a mass term
$V_{2}=m_{0}\sigma_{z}\tau_{x}$, which couples two BiS2-layers and opens a
gap. As demonstrated in Fig. 3(f), the Dirac fermion splits into two Karamer-
Weyl fermions at X point with gap $2m_{0}$=18.1 meV. Moreover, thanks to
surviving $\widetilde{M}_{z}$, all bands could be characterized by its
eigenvalues within the whole Brillouin zone. We find along X$\Gamma$ path two
bands with opposite $\widetilde{M}_{z}$ eigenvalues have to cross each other
and exhibits unique hourglass-like band connection. Actually, such band
crossing happens for arbitrary $k$-path connecting X and $\Gamma$ and forms a
closed nodal ring rather than isolated points. To confirm this, in Fig. 3(g)
we draw the 3D band spectrum and find a Weyl nodal ring surrounding X appears
in accompany with two Weyl points at X, in consistent with our symmetry
analysis. Besides, the band crossing along XM is also protected by distinctive
eigenvalues of $\widetilde{C}_{2y}$. So far as we know such glide mirror
protected hourglass nodal rings or nodal chains are mainly studied on 3D
systemsBzdusek  _et al._ (2016), our result provides a practical and tunable
Weyl nodal ring candidate in lower dimension. At Y point, despite sharing same
little group as X point, it hosts distinctive commutation relations, the zero-
order mass term $m_{0}\sigma_{z}\tau_{x}$ is forbidden by
$\widetilde{C}_{2y}=\sigma_{y}\tau_{x}$. Therefore, the quadruple degeneracy
is maintained at Y point while the $\mathcal{P}$-breaking allows first-order
term $V_{3}=w_{y}k_{y}\sigma_{x}\tau_{y}+w_{x}k_{x}\sigma_{y}\tau_{y}$, which
gives linear band splitting apart from Y along general $k-$path e.g. Y-X-Y as
shown in Fig. 3(h). Note the double band degeneracy are still maintained along
$\Gamma$Y and YM which are protected by the anticommutation relation
$\\{\widetilde{M}_{x},\widetilde{C}_{2y}\\}$=0 and by
$\mathcal{T}\widetilde{C}_{2y}=-1$, respectively. This type of Dirac point
that even exists without $\mathcal{PT}$ is very recently proposed in monolayer
SbSSnJin  _et al._ (2020).

Conclusion.— As a summary, in this work, we put forward the concept of 2D
Dirac semiconductor that host Dirac cone around band edge in nonsymmorphic
systems and then we propose an effective approach to search for such Dirac
semiconductor in certain multiple materials. We explicitly demonstrate an
practical Dirac semiconductor materials: TL-BiOS2, in which the Dirac cone is
composed of two spatially isolated Weyl cones, and provides layer-resolved
spin polarization. Moreover, with various TPTs can be induced in TL-BiOS2 by
imposing external strains or electric fields. Besides, it is worth mentioning
that starting from Dirac semiconductor, one can conveniently achieve Dirac
semimetal state by simply shifting Fermi level through imposing gate voltage,
element doping as well as band inversion. In fact, the F-substitution of O is
used to modify the Fermi level in bulk LaOBiS2 to realize high-temperature
superconducting stateYazici  _et al._ (2015), and some proposed Dirac
semimetal materials like HfGeTe can be taken as Dirac semiconductor with band
inversionGuan  _et al._ (2017). Our finding of 2D Dirac semiconductor not only
enriches topological quantum material families but also provides a new avenue
for exploring spin polarization, ferroelectric and optoelectronic properties
in topological quantum materials.

###### Acknowledgements.

This research was funded by the National Natural Science Foundation of China
(Grant No. 11874273 and No. 11974009) and the Key Project of Sichuan Science
and Technology Program (2019YFSY0044). We also thank the Sichuan Normal
University for financial support (No. 341829001). This research project was
supported by the High Performance Computing Center of Sichuan Normal
University, China. The work at BIT is supported by the National Key R&D
Program of China (Grant No. 2016YFA0300600), the National Natural Science
Foundation of China (Grants Nos. 11734003), the Strategic Priority Research
Program of Chinese Academy of Sciences (Grant No. XDB30000000).

## References

  * Neto  _et al._ (2009) A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Yan and Felser (2017) B. Yan and C. Felser, Rev. Mod. Phys. 8, 337 (2017).
  * Fang  _et al._ (2016) C. Fang, H. Weng, X. Dai, and Z. Fang, Chin. Phys. B 25, 117106 (2016).
  * Weng  _et al._ (2016) H. Weng, X. Dai, and Z. Fang, J. Phys.: Condens. Matter 28, 303001 (2016).
  * Gao  _et al._ (2019) H. Gao, J. F. Venderbos, Y. Kim, and A. M. Rappe, Annu. Rev. Mater. Rev. 49, 153 (2019).
  * Armitage  _et al._ (2018) N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Wang and Wang (2018) H. Wang and J. Wang, Chin. Phys. B 27, 107402 (2018).
  * Nagaosa  _et al._ (2020) N. Nagaosa, T. Morimoto, and Y. Tokura, Nat. Rev. Mater. 5, 621 (2020).
  * Wang  _et al._ (2016) M. X. Wang, Y. Xu, L. P. He, J. Zhang, X. C. Hong, P. L. Cai, Z. B. Wang, J. K. Dong, and S. Y. Li, Phys. Rev. B 93, 020503 (2016).
  * Chiu and Schnyder (2014) C.-K. Chiu and A. P. Schnyder, Phys. Rev. B 90, 205136 (2014).
  * Bian  _et al._ (2016) G. Bian, T.-R. Chang, R. Sankar, S.-Y. Xu, and et al., Nat. Commu. 7, 10556 (2016).
  * Li  _et al._ (2016) R. Li, H. Ma, X. Cheng, S. Wang, D. Li, Z. Zhang, Y. Li, and X.-Q. Chen, Phys. Rev. Lett. 117, 096401 (2016).
  * Yan  _et al._ (2015) B. Yan, B. Stadtmueller, N. Haag, S. Jakobs, and et al., Nat. Commu. 6, 10167 (2015).
  * Hirayama  _et al._ (2015) M. Hirayama, R. Okugawa, S. Ishibashi, S. Murakami, and T. Miyake, Phys. Rev. Lett. 114, 206401 (2015).
  * Sakano  _et al._ (2020) M. Sakano, M. Hirayama, T. Takahashi, S. Akebi, M. Nakayama, and et al., Phys. Rev. Lett. 124, 136404 (2020).
  * Ideue  _et al._ (2019) T. Ideue, M. Hirayama, H. Taiko, Takahashi, and et al., PNAS 116, 25530 (2019).
  * Zhang  _et al._ (2019) N. Zhang, G. Zhao, L. Li, P. Wang, and et al., arXiv 1906.06071 (2019).
  * Young and Kane (2015) S. M. Young and C. L. Kane, Phys. Rev. Lett. 115, 126803 (2015).
  * Guan  _et al._ (2017) S. Guan, Y. Liu, Z.-M. Yu, S.-S. Wang, Y. Yao, and S. A. Yang, Phys. Rev. Mater. 1, 054003 (2017).
  * Kowalczyk  _et al._ (2020) P. J. Kowalczyk, S. A. Brown, T. Maerkl, Q. Lu, and a. e. a. Chiu, ACS NANO 14, 1888 (2020).
  * Mounet  _et al._ (2018) N. Mounet, M. Gibertini, P. Schwaller, and et al., Nat. Nanotechnol. 13, 246 (2018).
  * (22) See Supplemental Information for the computational methods, the symmetry operators analysis, the tight-binding model, and the supplementary figures of candidates .
  * Qiu  _et al._ (2020) G. Qiu, C. Niu, Y. Wang, and et al., Nat. Nanotechnol. 15, 585 (2020).
  * Chang  _et al._ (2018) B. Chang, G.and Wieder, F. Schindler, and et al., Nat. Mater. 17, 978 (2018).
  * Li  _et al._ (2019) H. Li, S. Xu, Z. Rao, and et al., Nat. Commun. 10, 5505 (2019).
  * Zhang  _et al._ (2018) X. Zhang, B. Wang, X. Niu, Y. Li, Y. Chen, and J. Wang, Mater. Horiz. 5, 1058 (2018).
  * Huang and Yu (2020) C. Huang and H. Yu, 2D Mater. 7, 025023 (2020).
  * Fan  _et al._ (2018) X. Fan, D. Ma, B. Fu, C.-C. Liu, and Y. Yao, Phys. Rev. B 98, 195437 (2018).
  * Ishizaka  _et al._ (2011) K. Ishizaka, M. Bahramy, H. Murakawa, and a. e. a. Sakano, Nat. Mater. 10, 521 (2011).
  * Maaß  _et al._ (2016) H. Maaß, H. Bentmann, C. Seibel, C. Tusche, and et al., Nat. Commun. 7, 11621 (2016).
  * Zhang and Liu (2019) S.-H. Zhang and B.-G. Liu, Phys. Rev. B 100, 165429 (2019).
  * Di Sante  _et al._ (2013) D. Di Sante, P. Barone, R. Bertacco, and S. Picozzi, Adv. Mater. 25, 509 (2013).
  * Zhang  _et al._ (2014) X. Zhang, Q. Liu, J.-W. Luo, A. J. Freeman, and A. Zunger, Nat. Phys. 10, 387 (2014).
  * Liu  _et al._ (2015a) Q. Liu, X. Zhang, H. Jin, K. Lam, J. Im, A. J. Freeman, and A. Zunger, Phys. Rev. B 91, 235204 (2015a).
  * Liu  _et al._ (2016) Q. Liu, X. Zhang, and A. Zunger, Phys. Rev. B 93, 174119 (2016).
  * Guan and Luo (2020) S. Guan and J.-W. Luo, Phys. Rev. B 102, 184104 (2020).
  * Li and Chang (2009) J. Li and K. Chang, Appl. Phys. Lett. 95, 222110 (2009).
  * Li  _et al._ (2017) C.-H. Li, Y.-J. Long, L.-X. Zhao, L. Shan, Z.-A. Ren, J.-Z. Zhao, H.-M. Weng, X. Dai, Z. Fang, C. Ren, and G.-F. Chen, Phys. Rev. B 95, 125417 (2017).
  * Liu  _et al._ (2015b) Q. Liu, X. Zhang, L. Abdalla, A. Fazzio, and A. Zunger, Nano Lett. 15, 1222 (2015b).
  * Collins  _et al._ (2020) J. Collins, A. Tadich, W. Wu, and et al., Nature 564, 390 (2020).
  * Liu  _et al._ (2013) Q. Liu, Y. Guo, and A. J. Freeman, Nano Lett. 13, 5264 (2013).
  * Xu  _et al._ (2017) C.-Z. Xu, Y.-H. Chan, Y. Chen, P. Chen, and et al., Phys. Rev. Lett. 118, 146402 (2017).
  * Zhao  _et al._ (2020) C. Zhao, M. Hu, J. Qin, B. Xia, C. Liu, and et al., Phys. Rev. Lett. 125, 046801 (2020).
  * Cui  _et al._ (2018) C. Cui, F. Xue, W.-J. Hu, and L.-J. Li, npj 2D Mater. Appl. 2, 1 (2018).
  * Guan  _et al._ (2020) Z. Guan, H. Hu, X. Shen, P. Xiang, N. Zhong, J. Chu, and C. Duan, Adv. Electron. Mater. 6, 1900818 (2020).
  * Bzdusek  _et al._ (2016) T. Bzdusek, Q. Wu, A. Ruegg, M. Sigrist, and A. A. Soluyanov, Nature 538, 75 (2016).
  * Jin  _et al._ (2020) Y. J. Jin, B. B. Zheng, X. L. Xiao, Z. J. Chen, Y. Xu, and H. Xu, Phys. Rev. Lett. 125, 116402 (2020).
  * Yazici  _et al._ (2015) D. Yazici, I. Jeon, B. D. White, and M. B. Maple, Physica C: Superconductivity and its Applications 514, 218 (2015).


# Towards solving the BCS Hamiltonian gap in Near-Term Quantum Computers

Nahum Sá  Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150,
22290-180 Rio de Janeiro, Brazil<EMAIL_ADDRESS>Ivan S. Oliveira  Centro
Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de
Janeiro, Brazil<EMAIL_ADDRESS>Itzhak Roditi  Centro Brasileiro de Pesquisas
Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, Brazil
<EMAIL_ADDRESS>
###### Abstract

In this work, using a NISQ framework, we obtain the gap of a BCS Hamiltonian.
This could lead to interesting implications for superconductivity research.
For such task, we choose to use the Variational Quantum Deflation and analyze
the hardware restrictions that are needed to find the energy spectra on
current quantum hardware.

We also compare two different kinds of classical optimizers, Constrained
Optimization BY Linear Approximations (COBYLA) and Simultaneous Perturbation
Stochastic Approximation (SPSA), and study the effect of decoherence caused by
the presence of noise when using simulations in real devices. We implement
this method for a system with both 2 and 5 qubits. Furthermore, we show how to
approximate the gap within one standard deviation, even with the presence of
noise.

##  1 Introduction

At the current time, we may consider that we are in the Noisy Intermediate-
Scale Quantum (NISQ) [1] era of Quantum Computing which refers to devices with
50-100 qubits, without Quantum Error Correction. In other words, we are still
far from devices with perfect qubits and quantum operations. Thus, in order to
attain an acceptable performance, we need to use algorithms that are suited
for those limitations. A very valuable class of algorithms appropriate to deal
with such limitations, and that may be very useful in physics as well as in
many other fields, e.g. machine learning, are the so-called Variational
Quantum Algorithms (VQA).

Variational Quantum Algorithms is a class of hybrid classical-quantum
algorithms that shows noise resilience due to the use of parametric quantum
circuits [2], thus being able to be implemented in NISQ devices. Variational
Quantum Algorithms are universal [3] and have been proposed for various
problems, such as Combinatorial Optimization [4], Quantum Chemistry [5],
factoring [6], Machine Learning [7] and compilation [8].

One problem that can be dealt with by the use of such techniques and has the
potential to lead to interesting outcomes is superconductivity. This is fairly
well modeled using a pairing Hamiltonian.

BCS superconductivity is an appealing framework for testing VQA.
Experimentally, the simulation of pairing Hamiltonians have been shown using a
Nuclear Magnetic Resonance (NMR) setup [9]. However, due to intrinsic NMR
experimental characteristics, there are several limitations for scalability
which can be handled using VQAs.

In order to simulate Fermionic Hamiltonians, in digital quantum computers, one
needs to map fermions to qubits. There are many mapping techniques used for
quantum chemistry, for instance Jordan-Wigner [10], and Bravyi-Kitaev Mapping
[11]. One of such qubit mappings proposed by Wu et al. [12] has been only
implemented for NMR quantum computers, although it can be used for circuit-
based quantum computers and more specifically NISQ devices using Variational
Quantum Algorithms. This paper aims to use this mapping to solve the BCS
Hamiltonian. On the other hand, one can also can simulate Fermionic
Hamiltonians using analog quantum simulations [13, 14].

In order to show that this mapping is suitable for circuit-based quantum
computers, we will use the Fermion-to-Qubit mapping made by Wu et al. [12] to
solve the BCS problem with both $N=5$ and $N=2$ qubits. This also extends the
result of [9] in a setting where it is not needed to use any simplification of
the parameters of the Hamiltonian, leading to a path into solving this problem
with an arbitrary number of cooper pairs which is limited only by quantum
computer hardware.

In order to demonstrate this, our goal is to find the energy spectra. One can
use Variational Quantum Algorithms. We choose to use the Variational Quantum
Deflation algorithm [15] which could lead to restrictions of the topology of
quantum computers. Fortunately, those restrictions are met by real devices
which could be used to solve the BCS Hamiltonian. However, this method can
work with any algorithm able to find excited states on quantum computers, such
as Quantum Subspace Expansion methods [16, 17].

The paper is organized as follows: Section 2 presents the BCS Hamiltonian, the
qubit mapping used, explains the Variational Quantum Deflation algorithm, and
the ansatz used. Section 3 shows the numerical simulation results for an ideal
quantum computer and a noisy quantum computer assuming a Thermal Relaxation
Error noise model. Lastly, in Section 4, we conclude the paper and present
possible ideas for future research. In A we explain the optimizers used in
further detail. In B we further explain the methods that can be used to
measure the overlap between two states using a quantum computer and the
hardware restrictions.

##  2 Methods

First, we need to define the Fermionic Hamiltonian that we want to solve. The
general Bardeen–Cooper–Schrieffer (BCS) [18] Hamiltonian is given by:

$H_{BCS}=\sum_{m=1}^{N}\frac{\epsilon_{m}}{2}(n_{m}+n_{-m})\
+\sum_{m,k=1}^{N}V^{+}_{ml}c^{\dagger}_{m}c^{\dagger}_{-m}c_{-l}c_{l}$ (1)

where $n_{\pm m}=c^{\dagger}_{\pm m}c_{\pm m}$ is the number operator, and the
matrix elements $V^{+}_{ml}=\langle m,-m|V|l,-l\rangle$ are real and can be
calculated efficiently for any given problem. For simplicity V can be treated
as a constant which yields a good approximation to many superconductors [19].

An important property that can be obtained from the Hamiltonian is the gap,
which is defined as $2\Delta_{n}\equiv E_{n,1}-E_{n,0}$ , $n$, being the
number of Cooper pairs (see [12] for an NMR simulation experiment where the
energy spectrum of the Hamiltonian is determined). For instance, when $n=0$,
one has for the gap between the first excited state and the ground state,
$2\Delta_{0}\equiv E_{1}-E_{0}$. And for $n=1$, one has the gap between the
first and second excited states, $2\Delta_{1}\equiv E_{2}-E_{1}$. In the
present work, we will focus both on the $n=1$ gap and the $n=0$ gap.

In order to solve this problem in a quantum computer, one needs to map the BCS
Hamiltonian into a qubit Hamiltonian with a one-to-one mapping. In this paper
we will use one of the mappings presented in Wu et al. [12] which maps
Fermions to qubits, given by the following representation:

$H_{Q}=\sum_{m=1}^{N}\frac{\epsilon_{m}}{2}\sigma^{Z}_{m}+\frac{V}{2}\sum_{l>m=1}^{N}(\sigma^{x}_{m}\sigma^{x}_{l}+\sigma^{y}_{m}\sigma^{y}_{l})$
(2)

One can find the gap by finding the energy spectra of the qubit Hamiltonian
(Eq. 2) since the mapping assures that the qubit Hamiltonian has the same
energy spectra as that of the BCS Hamiltonian. The goal of this paper is to
find the energy spectra using a variational quantum algorithm called
Variational Quantum Deflation [15].

The Variational Quantum Deflation (VQD) algorithm is an extension of the
Variational Quantum Eigensolver (VQE) [5] algorithm, allowing us to
approximate excited states of a desired Hamiltonian, thus giving us access to
their energy values.

The original VQE algorithm only finds the ground energy of a given Hamiltonian
by using the decomposition of such Hamiltonian into Pauli strings and
minimizing the expected value of the Hamiltonian:

$E(\theta)=\langle\psi(\theta)|H|\psi(\theta)\rangle=\sum_{i}c_{i}\langle\psi(\theta)|P_{i}|\psi(\theta)\rangle$
(3)

Where $P_{i}$ are Pauli strings that are of polynomial size. In order to
measure this expectation value one needs to change from the computational
basis into the X and Y basis, thus this approach has low depth and is suited
for NISQ devices. By minimizing the cost function (Eq. 3) it is possible to
approximate the ground state of a Hamiltonian with high precision.

In order to get the energy of excited states, the cost function of the VQE
needs to be modified by taking into account that all eigenstates of the
Hamiltonian are orthogonal to each other. This is done by minimizing the
overlap between all eigenstates discovered using the algorithm. Thus, to
calculate the k-th excited state of a given Hamiltonian, $H$, the parameters
$\lambda_{k}$ must be optimized for a parametrized state
$|\psi(\lambda_{k})\rangle$ for a given cost function. The cost function that
needs to be minimized for the VQD is given by:

$F(\lambda_{k})=\langle\psi(\lambda_{k})|H|\psi(\lambda_{k})\rangle+\sum_{i=0}^{k-1}\beta_{i}\big{|}\langle\psi(\lambda_{k})|\psi(\lambda_{i})\rangle\big{|}^{2}$
(4)

The first term is the expected value of the energy, which can be obtained
measuring the Pauli decomposition of the Hamiltonian $H$ just like the VQE
cost function. On the other hand, the second term is not as trivial, one needs
to calculate the overlapping term between all eigenstates until the desired
k-th eigenstate, this term ensures that all eigenstates are orthogonal to each
other.

It is possible to measure the overlap between two quantum states on quantum
computers, this will be explained in more detail in B. In addition, those
methods imposes hardware restrictions, which are also discussed in the
appendix.

One of the downsides of the VQD algorithm is that the hyperparameter
$\beta_{i}$ has to be chosen, but, according to [15], it suffices to choose
$\beta_{i}>E_{k}-E_{i}$, where $E_{k}$ is the desired energy state and $E_{i}$
all energy states below the desired state, since this paper aims to find the
energy spectra of the Hamiltonian, we choose
$\beta_{i}>E_{\text{max}}-E_{\text{min}}$, where $E_{\text{max}}$ is the
highest value of the energy and $E_{\text{min}}$ the minimum value of the
energy. Another aspect to be careful with is that if the real ground state is
not satisfactorily approximated, there may be an error propagation due to the
initial state not being well-prepared.

In the following subsection, we will describe different kinds of variational
ansatz and our choice of a Hardware-Efficient Ansatz in more detail.

###  2.1 Variational Ansatz

In order to build a variational ansatz for the Variational Quantum Algorithm,
two conflicting goals should be met: The ansatz must be as general as possible
in order to create states that would approximate all eigenstates of the
Hamiltonian, but also needs to have a polynomial number of parameters in order
to express our ansatz efficiently on a quantum computer.

In the literature there are two distinct paths for constructing the
variational ansatz: 1) You can create an ansatz based on the problem’s
knowledge, in this case we can see as an example the UCCSD ansatz [20]; 2)
Since the main goal is to run the Variational Quantum Algorithm in real
hardware one needs to take into account hardware restrictions on connectivity
and the gate set and construct an ansatz that is suited for the Hardware that
the algorithm will run on, this is commonly called Hardware Efficient Ansatz
introduced by Kandala et al. [21].

Our choice of ansatz for this paper is the hardware efficient ansatz. This
kind of ansatz is depicted in FIGURE 1 and is composed by blocks of two kind
of unitary transformations. The first block is made of one or more hardware-
native parametrized one qubit gates ($R_{j}(\theta_{i})$) and depicted on the
figure as $U(\mathbb{\theta})$.

The second block is constructed using two qubit gates that can generate
entanglement, which are generally called entangling gates, the most common
entangling gate is the CNOT, however there are other entangling gates such as
iSWAP used on the Sycamore chip made by Google [22]. Since entangling gates
are the gates which have lower coherence, it is imperative to use them
according to the hardware’s connection graph in order to avoid using swap
gates, which are costly for NISQ devices and can lead to non-local errors.

Thus, the ansatz can be written as following:

$\left|\psi(\theta)\right\rangle=U_{\text{ENT}}\ U(\theta)\ \dots\
U_{\text{ENT}}\ U(\theta)\left|\psi(0)\right\rangle$ (5)

Depth, $d$, is defined as the number of both the single qubit rotations block
and the entangling block. The dependency on the number of parameters regarding
the depth is given by $2Nd$, where $N$ is the number of qubits, yielding a
polynomial number of parameters, which ensures that we are working with an
efficient ansatz construction.

Figure 1: Hardware Efficient ansatz where the unitary $U_{\mathrm{ENT}}$ is
constructed using the connection graph of the Quantum Hardware, and the
unitary $U(\mathbb{\theta})$ is made of hardware-native parametrized
rotations, in our case it is the combination of $R_{Y}$ and $R_{Z}$ rotations.

##  3 Results

In order to demonstrate our proposed method, we choose to simulate it both on
ideal and noisy quantum computers for the case of $N=2$ and $N=5$ qubits,
because the system is solvable using direct diagonalization, and we can
compare results between the result of direct diagonalization and the result
from the quantum computer. Both simulations will take into account that we
only take a finite amount of measurements (10,000 shots) from the quantum
computer.

The Hamiltonian that we aim to solve is given by equation (2) with $N=2,5$. An
interesting parameter that we can obtain for this Hamiltonian, as mentioned
before, is the gap: $2\Delta_{n}=E_{n,1}-E_{n,0}$.

This section will be divided in two subsections, in which different parameters
of the simulation will be considered. The first section 3.1 will be the
simulation performed by a perfect quantum computer, and the second 3.2 will be
the simulation subject to a noise model.

###  3.1 Ideal Quantum Computer

We will consider a perfect Quantum computer without any noise for the first
experiment. The purpose of this experiment is two-fold: firstly, to show that
the algorithm works in an ideal setting. Secondly, it will be a fair
comparison when we analyze the case with noise, which is of main importance
for running Variational Quantum Algorithms in NISQ devices.

In order to estimate the depth needed to simulate the Hamiltonian we run the
algorithm changing the depth of the ansatz and estimate the gap, this is shown
on FIGURES 2 and 3, for 2 and 5 qubits respectively. For the case of 5 qubits,
we choose to only use the COBYLA optimizer, because we have observed that it
worked better for the case of 2 qubits.

Figure 2: A plot analyzing the depth needed to solve the problem for a two-
qubit Hamiltonian using a Hardware Efficient Ansatz, where we took the
statistics of 50 runs of the algorithm with random initialization and finite
number of shots. Comparing between COBYLA and SPSA optimizers. Analysis of the
Optimizer performance, where we used $c=0.7$ and took the average of the last
25 $\lambda_{k}$ for the final $\lambda_{k}$ value on the SPSA. As the depth
of the circuit increases, the algorithm converges to the solution that is
obtained through direct diagonalization. Figure 3: A plot analyzing the depth
needed to solve the problem for a five qubit Hamiltonian using a Hardware
Efficient Ansatz, where we took the statistics of 50 runs of the algorithm
with random initialization and finite number of shots. We choose only to use
the COBYLA optimizer in order to show that our algorithm works for a higher
number of qubits.

According to the figures, we observe that for the SPSA optimizer the increase
in depth can also lead to higher variance on the samples, this is expected
since increasing the depth of the ansatz also leads to a linear increment in
the number of parameters, this shows that for this instance the gradient-free
optimizer works better for this kind of problem.

In addition, we also see that a depth of 3 is needed to find a reasonable
solution for both optimizers to solve the BCS problem for both cases, because
it is within 1 standard deviation of the statistical sample and is the
smallest depth that satisfies this condition.

After defining the depth needed to solve the problem, we benchmark the problem
of finding the $n=1$ gap of the two-qubit Hamiltonian by changing the coupling
parameter $V$, with parameters $\epsilon=\epsilon_{1}=\epsilon_{2}=3$. This is
done by calculating the first and second excited states of the Hamiltonian.
Both optimizers are able to solve the gap problem within 1 standard deviation,
which is demonstrated on FIGURE 4.

Figure 4: Measuring the gap varying V for the two-qubit Hamiltonian, where we
used $c=0.7$ and took the average of the last 25 $\lambda_{k}$ for the final
$\lambda_{k}$ value for the SPSA optimizer, and we took the statistics of 10
runs of the algorithm with random initialization.

However, only the COBYLA optimizer is able to show the expected trend on the
mean value. This is due to the need to tune the hyperparameters for each case
of the SPSA optimizer.

For the case of the five qubit Hamiltonian, we change the coupling parameter
$V$, with $\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=\epsilon_{5}=3$ and
$\epsilon_{4}=4$, and estimate the $n=0$ gap, which is done by calculating the
gap between the ground state and the first excited state of the Hamiltonian.
This is represented in FIGURE 5.

Figure 5: Measuring the gap varying V for the five qubit Hamiltonian using the
depth 3, where we took the statistics of 20 runs of the algorithm with random
initialization.

###  3.2 Noisy Quantum Computer

In this work we choose the Thermal Relaxation Error [23] as our noise model
which represents the thermalization of the qubit towards the equilibrium state
at the temperature of the environment. This is an irreversible process that is
governed by the gate time $T_{g}$, relaxation times $T_{1}$, $T_{2}$, and
device temperature, $T$, which all can be obtained during the device
calibration.

Since the device is kept at a temperature very close to absolute zero, it is
safe to assume that $T\approx 0$. This assumption implies that there is no
excited state population, which is given by the equation:

$p_{e}=\bigg{(}1+\exp\big{(}\frac{2hf}{k_{B}T}\big{)}\bigg{)}^{-1}$ (6)

Where $T$ is the device’s Temperature, $f$ is the qubit frequency, $k_{B}$ is
Boltzmann’s constant, and $h$ is Planck’s constant. Under this assumption, we
have two sources of noise, dephasing and reset to $\left|0\right\rangle$.

The relaxation error rates are defined as $\epsilon_{T_{1}}=e^{-T_{g}/T_{1}}$
and $\epsilon_{T_{2}}=e^{-T_{g}/T_{2}}$ and the reset probability is defined
as $p_{\mathrm{reset}}=1-\epsilon_{T_{1}}$.

There are two behaviors for this kind of noise: $T_{2}\leq T_{1}$, and
$T_{2}>T_{1}$.

For the case, $T_{2}\leq T_{1}$ we can represent the Thermal Relaxation Error
as a probabilistic mixture of reset operations and unitary errors.

The probability of a dephasing happening is given by
$p_{Z}=\frac{1-p_{\mathrm{reset}}}{2}\bigg{(}1-\frac{\epsilon_{T_{2}}}{\epsilon_{T_{1}}}\bigg{)}$.
The probability of a reset to $\left|0\right\rangle$ operation occurring is
given by $p_{r_{\left|0\right\rangle}}=(1-p_{e})p_{\mathrm{reset}}$. The
probability of no noise occurring is given by
$p_{I}=1-p_{Z}-p_{\mathrm{reset}}$.

Under these conditions, the error model has a Kraus representation as:

$\begin{split}K_{I}=\sqrt{p_{I}}I\\\ K_{Z}=\sqrt{p_{Z}}Z\\\
K_{\mathrm{reset}}=\sqrt{p}_{\mathrm{reset}}\left|i\right\rangle\left\langle
0\right|\end{split}$ (7)

The relaxation noise is represented as the following channel:

$\rho\mapsto\mathcal{N}(\rho)=\sum_{j}K_{j}\rho K_{j}^{\dagger}$ (8)

For the case $T_{2}>T_{1}$, the error can be represented by a Choi-Matrix
representation. Considering $T\approx 0$, the Choi-Matrix of the Thermal
relaxation noise is given by [24]:

$\begin{pmatrix}1&0&0&\epsilon_{T_{2}}\\\ 0&0&0&0\\\
0&0&p_{\mathrm{reset}}&0\\\ \epsilon_{T_{2}}&0&0&1-p_{\mathrm{reset}}\\\
\end{pmatrix}$ (9)

Using the Choi-Matrix representation, it is possible to find the Kraus
operators and consequently find the channel for the noise.

It is important to evaluate if the Variational Quantum Algorithm that we are
using is resilient to noise in order to run in NISQ devices. Thus, we follow
the same methodology as the Ideal Quantum Computer section 3.1 under the
presence of noise.

We first vary the ansatz depth in order to find the number of layers on the
ansatz for the case with two qubits, and we see a similar trend when comparing
the analysis with and without noise, however the SPSA optimizer seems to
behave not as well as the COBYLA optimizer which shows convergence when the
depth increases. This shows that just as without noise, we have a good
convergence when the depth of the ansatz is equal to 3, just like the case
without noise.

Figure 6: A plot analyzing the depth needed to solve the problem for a two
qubit Hamiltonian using a Hardware Efficient Ansatz, where we took the
statistics of 10 runs of the algorithm with random initialization and compare
between COBYLA and SPSA optimizers. Analysis of the Optimizer performance,
where we used $c=0.7$ and took the average of the last 25 $\lambda_{k}$ for
the final $\lambda_{k}$ value on the SPSA.

For the case with five qubits, we use only the COBYLA optimizer in order to
show that our method works for more than two qubits, FIGURE 7. The noise leads
to an offset which has been removed. We see that, it has a similar behavior
than the case without noise, thus choosing the depth of the ansatz equals to 3
works for the five qubit case.

Figure 7: A plot analyzing the depth needed to solve the problem for a five
qubit Hamiltonian using a Hardware Efficient Ansatz, where we took the
statistics of 50 runs of the algorithm with random initialization for the
COBYLA optimizer. The results are adjusted by a constant factor due to the
presence of noise.

After defining the depth, we implement for the same case, varying $V$ on the
Hamiltonian and using $\epsilon=\epsilon_{1}=\epsilon_{2}=3$. We observe that
for this case, both COBYLA and SPSA optimizers have demonstrated a linear
trend, which is expected when the coupling parameter $V$ is changed.

Even though we used $\epsilon=\epsilon_{1}=\epsilon_{2}$, we could choose
$\epsilon_{1}\neq\epsilon_{2}$ without any loss to our algorithm, these values
were chosen because it would demonstrate a linear trend on the BCS gap for the
case of two qubits, other conditions would lead to non-linear relation when
varying the coupling parameter.

Figure 8:  Measuring the gap varying the coupling constant, $V$, for the two
qubit Hamiltonian. We used $c=0.7$ and took the average of the last 25
$\lambda_{k}$ for the final $\lambda_{k}$ value on the SPSA.

Now we use the depth of 3, for the case varying the coupling constant, $V$,
with $\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=\epsilon_{5}=3$ and
$\epsilon_{4}=4$, and estimate the $n=0$ gap. This is represented in FIGURE 9.

Figure 9: Measuring the gap varying the coupling constant, $V$, for the two
qubit Hamiltonian. We use the COBYLA Optimizer performance, to show that our
method works with a high qubit number. The results are adjusted by a constant
factor due to the presence of noise.

##  4 Conclusion

In this paper we have explored solving the BCS Hamiltonian through the Wu et
al. [12] mapping using Near-Term Quantum Computers through a Variational
Quantum Algorithm.

The algorithm that was used is called Variational Quantum Deflation (VQD) and
we were able to obtain the $n=1$ gap of a two qubit BCS Hamiltonian within 1
standard deviation considering random initialization, showing that the
algorithm works and is not highly dependent of the initialization of the
random parameters. We also obtain the $n=0$ gap of a five qubit BCS
Hamiltonian, which shows that our methods scale for higher qubit instances.

We also analyzed both gradient and gradient-free optimizers in order to search
for variational parameters that solve our task. For the gradient-based
optimizer we chose Simultaneous Perturbation Stochastic Approximation (SPSA)
and for the gradient-free we chose constraint functions by linear
interpolation (COBYLA). Showing that even in the presence of Thermal
Relaxation noise the COBYLA optimizer have a good convergence for solving the
BCS gap problem, even though it is assumed to be a bad optimizer when the cost
function has stochastic noise.

This work shows promising ways to solve the BCS Hamiltonian in near-term
quantum computers, which could be extended for an arbitrary number of cooper
pairs. For future works, we aim to explore an ansatz that is tailor-made for
the Hamiltonian in order to avoid the Barren Plateaus problem and have a
better optimization convergence. Another direction is to analyze in more depth
the objective function landscape of the VQD algorithm and see if the barren
plateau problem is present in this Variational Quantum Algorithm.

The authors would like to thank João Ribeiro and Filipe Melo for useful
comments on this manuscript. This study was financed in part by the
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) –
Finance Code 001, and by the Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq). We acknowledge the use of IBM Quantum services for this
work. The views expressed are those of the authors, and do not reflect the
official policy or position of IBM or the IBM Quantum team. This manuscript
used Quantikz [25] for drawing the circuits and Qiskit [26] for its
simulations. All source code is available on GitHub:
https://github.com/nahumsa/volta.

## References

  * [1] Preskill J 2018 Quantum 2 79
  * [2] Reiner J M, Wilhelm-Mauch F, Schön G and Marthaler M 2019 Quantum Science and Technology 4 035005
  * [3] Biamonte J 2021 Physical Review A 103 L030401
  * [4] Farhi E, Goldstone J and Gutmann S 2014 arXiv preprint arXiv:1411.4028
  * [5] Peruzzo A, McClean J, Shadbolt P, Yung M H, Zhou X Q, Love P J, Aspuru-Guzik A and O’brien J L 2014 Nature communications 5 4213
  * [6] Anschuetz E, Olson J, Aspuru-Guzik A and Cao Y 2019 Variational quantum factoring International Workshop on Quantum Technology and Optimization Problems (Springer) pp 74–85
  * [7] Schuld M, Bocharov A, Svore K M and Wiebe N 2020 Physical Review A 101 032308
  * [8] Khatri S, LaRose R, Poremba A, Cincio L, Sornborger A T and Coles P J 2019 Quantum 3 140
  * [9] Yang X, Wang A M, Xu F and Du J 2006 Chemical Physics Letters 422 20–24 URL https://doi.org/10.1016/j.cplett.2006.02.023
  * [10] Jordan P and Wigner E 1928 Zeitschrift für Physik 47 631–651 URL https://doi.org/10.1007/bf01331938
  * [11] Bravyi S B and Kitaev A Y 2002 Annals of Physics 298 210–226
  * [12] Wu L, Byrd M and Lidar D 2002 Physical Review Letters 89 057904
  * [13] Parsons M F, Mazurenko A, Chiu C S, Ji G, Greif D and Greiner M 2016 Science 353 1253–1256 URL https://doi.org/10.1126/science.aag1430
  * [14] Bloch I, Dalibard J and Nascimbène S 2012 Nature Physics 8 267–276 URL https://doi.org/10.1038/nphys2259
  * [15] Higgott O, Wang D and Brierley S 2019 Quantum 3 156
  * [16] Colless J I, Ramasesh V V, Dahlen D, Blok M S, Kimchi-Schwartz M E, McClean J R, Carter J, de Jong W A and Siddiqi I 2018 Phys. Rev. X 8(1) 011021 URL https://link.aps.org/doi/10.1103/PhysRevX.8.011021
  * [17] Nakanishi K M, Mitarai K and Fujii K 2019 Phys. Rev. Research 1(3) 033062 URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033062
  * [18] Bardeen J, Cooper L N and Schrieffer J R 1957 Physical review 108 1175
  * [19] Waldram J 1976 Reports on Progress in Physics 39 751
  * [20] McArdle S, Endo S, Aspuru-Guzik A, Benjamin S C and Yuan X 2020 Reviews of Modern Physics 92 015003
  * [21] Kandala A, Mezzacapo A, Temme K, Takita M, Brink M, Chow J M and Gambetta J M 2017 Nature 549 242–246
  * [22] Arute F, Arya K, Babbush R, Bacon D, Bardin J C, Barends R, Biswas R, Boixo S, Brandao F G, Buell D A et al. 2019 Nature 574 505–510
  * [23] Georgopoulos K, Emary C and Zuliani P 2021 arXiv e-prints arXiv:2101.02109 (Preprint 2101.02109)
  * [24] Blank C, Park D K, Rhee J K K and Petruccione F 2020 npj Quantum Information 6 1–7
  * [25] Kay A 2018 arXiv e-prints arXiv:1809.03842 (Preprint 1809.03842)
  * [26] 2019 Qiskit: An open-source framework for quantum computing
  * [27] McClean J R, Boixo S, Smelyanskiy V N, Babbush R and Neven H 2018 Nature communications 9 1–6
  * [28] Holmes Z, Sharma K, Cerezo M and Coles P J 2021 arXiv preprint arXiv:2101.02138
  * [29] Wang S, Fontana E, Cerezo M, Sharma K, Sone A, Cincio L and Coles P J 2020 arXiv preprint arXiv:2007.14384
  * [30] Arrasmith A, Cerezo M, Czarnik P, Cincio L and Coles P J 2020 arXiv e-prints arXiv:2011.12245 (Preprint 2011.12245)
  * [31] Spall J C et al. 1992 IEEE transactions on automatic control 37 332–341
  * [32] Maryak J L and Chin D C 2001 Global random optimization by simultaneous perturbation stochastic approximation Proceedings of the 2001 American Control Conference.(Cat. No. 01CH37148) vol 2 (IEEE) pp 756–762
  * [33] Spall J C 1998 IEEE Transactions on aerospace and electronic systems 34 817–823
  * [34] Powell J 1994 A direct search optimization method that models the objective and constraint functions by linear interpolation Advances in Optimization and Numerical Analysis, Proceedings of the 6th Workshop on Optimization and Numerical Analysis, Oaxaca, Mexico ed Gomez S and Hennart J P (Dordrecht, The Netherlands: Kluwer Academic Publishers) pp 51–67
  * [35] Powell J 1998 Acta Numerica 7 287–336
  * [36] Buhrman H, Cleve R, Watrous J and de Wolf R 2001 Physical Review Letters 87 167902
  * [37] Cincio L, Subaşı Y, Sornborger A T and Coles P J 2018 New Journal of Physics 20 113022
  * [38] 2021 Amazon Braket https://aws.amazon.com/braket/
  * [39] 2021 IBM Quantum Experience https://quantum-computing.ibm.com/
  * [40] Havlíček V, Córcoles A D, Temme K, Harrow A W, Kandala A, Chow J M and Gambetta J M 2019 Nature 567 209–212 URL https://doi.org/10.1038/s41586-019-0980-2

## Appendix

##  Appendix A Classical Optimizer

There are plenty of optimization algorithms that are used on Variational
Quantum Algorithms. In this paper, we will focus on two algorithms for
optimization: The Simultaneous Perturbation Stochastic Approximation (SPSA)
and Constrained Optimization BY Linear Approximations (COBYLA).

The problem of barren plateaus [27] is present in gradient-based optimizers,
which means that the loss function gradient becomes exponentially smaller
making it harder to minimize the cost function. The barren plateaus problem
can be induced by various factors such as high ansatz expressiveness [28] or
even presence of noise [29].

Initially, the presence of barren plateaus was found only on gradient-based
optimizers, but recently there is evidence that it affects gradient-free
optimizers as well [30].

Therefore, in order to examine the performance of classical optimizers, it is
fair to compare gradient-free and gradient based optimizers. We thus choose
COBYLA as the gradient-free optimizer and SPSA as the gradient-based
optimizer.

###  A.1 SPSA

The Simultaneous Perturbation Stochastic Approximation (SPSA) [31] is an
algorithm that is robust with respect to stochastic perturbations of the
function with which it will be minimized and has been widely used in
variational quantum algorithms that deal with NISQ devices [21].

The SPSA algorithm works with only two evaluations using a symmetrical
Bernoulli distribution $\Delta_{i}$ and two zero converging sequences, $c_{i}$
and $a_{i}$. Since we want to minimize the cost function of the VQD algorithm
(eq. 4) we must take the gradient $g(\lambda_{k})$ this will be done by
sampling the cost function $F(\lambda_{k})$, since we only take the sample of
the gradient and do not calculate it analytically we name it
$\hat{g}(\lambda_{k})$ and we obtain it using the formula:

$\hat{g}(\lambda_{k}^{(i)})=\frac{F(\lambda_{k}^{(i)}+c_{n}\Delta_{i})-F(\lambda_{k}^{(i)}-c_{i}\Delta_{i})}{2c_{i}\Delta_{i}}$
(10)

where $\lambda_{k}^{(i)}$ is the $\lambda_{k}$ parameter in the
$i^{\text{th}}$ iteration of the optimization algorithm.

After evaluating the gradient, the algorithm updates the parameters using the
following rule:

$\lambda_{k}^{(i+1)}=\lambda_{k}^{(i)}-a_{i}\hat{g}(\lambda_{k}^{(i)})$ (11)

It has been proven [31, 32] that this algorithm can converge assuming
stochastic fluctuations in the cost function evaluation. The sequences $a_{i}$
and $c_{i}$ are chosen as:

$\begin{split}c_{i}=\frac{c}{i^{\gamma}}\\\
a_{i}=\frac{a}{i^{\alpha}}\end{split}$ (12)

where the values are optimally chosen [33] as
$\\{\alpha,\gamma\\}=\\{0.602,0.101\\}$. The value $a$ controls how large is
the update of the parameters and $c$ controls the gradient update, thus if
there are large statistical fluctuations on the cost function we must choose a
large $c$ in order to the gradient evaluation to be accurate, according to
experiments explained in section 3 we chose $c=0.7$.

In addition to that we use a calibration method as in [21] in order to adjust
the $a$ parameter, which starts at $a=2\pi/10$, choosing different directions
of $\Delta_{i}$, the inverse formula used to calibrate $a$ is:

$a=\frac{2\pi}{5}\frac{c}{\big{\langle}F(\lambda_{k}^{(1)}+c_{n})-F(\lambda_{k}^{(1)}-c_{i})\big{\rangle}_{\Delta_{1}}}$
(13)

Finally, in order to get the optimized $\lambda_{k}$ values, we follow [21]
and average over the last 25 iterations in order to suppress statistical
fluctuations on parameter values.

###  A.2 COBYLA

Based on a direct search algorithm that models the objective and constraint
functions by linear interpolation [34] COBYLA is a classical gradient-free
optimizer [35]. It is based on an algorithm that is iterative, in which each
iteration produces, by interpolation at the vertices of a simplex, a linear
approximation of the objective and constraint functions. The iterative change
of values is restricted by a trust region, $\rho$, to assure that the problem
has a finite solution.

It is expected that this kind of optimizer works well when stochastic
perturbations are minimal, because it supposes that function evaluation is
exact. This could, anyhow, lead to problems when working with real quantum
computers that are subject to low decoherence times. However, there is no
proof in the literature that this algorithm will not converge with stochastic
perturbations, thus we aim to test this case to see if it converges with the
presence of noise.

For the trust region, the user needs to choose an initial
($\rho_{\text{init}}$) value that represents the exploration on the objective
function landscape, which decreases by half after each iteration in order to
find the minimum of the objective function subjected to the given constraints.
For our optimization procedure, we choose $\rho_{\text{init}}=1$.

##  Appendix B Measuring the overlap of two states on a quantum computer

In order to use the Variational Quantum Deflation algorithm, we require
measuring the overlap between two quantum states. This can be done by three
techniques: SWAP test, Destructive SWAP test, and Amplitude Transition. In
this section, we will explain each method and comment its advantages and
disadvantages.

SWAP test is an algorithm that gives as a result the state overlap between two
quantum states $\rho$ and $\sigma$. There are two proposals for a SWAP test in
the literature which are equivalent to each other. The original SWAP test [36]
that uses a controlled-SWAP operation and an auxiliary qubit to measure the
overlap, and the Destructive SWAP test [37] which only uses first neighbor
connectivity, CNOT and 1 qubit gates to measure the overlap.

The SWAP test is depicted on FIGURE 10 and consists of the application of
Hadamard and Controlled-SWAP gates, and it is measured using the auxiliary
qubit. The probability on measuring the $\left|0\right\rangle$ state encodes
the overlap between $\rho$ and $\sigma$

$|\langle\rho|\sigma\rangle|^{2}=2\ \bigg{(}P(0)-\frac{1}{2}\bigg{)}$ (14)

where $P(0)$ is the probability of the outcome on the auxiliary qubit being 0.

Figure 10: Circuit representation of the SWAP test for states $\rho$ and
$\sigma$, this is done using controlled-SWAP gate and 1 qubit gates.

This is a solid approach for learning the overlap between two quantum states,
however using a controlled-SWAP gate demands high coherence times and high
connectivity between qubits, those two requirements are not available on
current quantum hardware. For NISQ devices we must use an algorithm that is
suited for the device’s restrictions on both connectivity and gate errors, the
algorithm that solve both those restrictions in certain conditions is the
Destructive SWAP test.

The Destructive SWAP test is represented on FIGURE 11 which consists of CNOT
and Hadamard Gates and no auxiliary qubit. This algorithm has depth $O(1)$
which makes it suited for NISQ devices and has demonstrated its superiority
against the original SWAP Test [37].

In order to obtain the state overlap, one ought to do measurements in the Bell
basis for the corresponding qubit in each state. Those measurements can be
achieved with a classical post-processing step that scales linearly with qubit
size. This post-processing can be interpreted as the expectation value of a
controlled-Z observable.

Figure 11: Destructive SWAP test for states $\rho$ and $\sigma$, after the
measurement on the Bell basis, it is needed to do a post-processing
classically.

It is important to analyze the topology needed for these algorithms to be used
in real quantum hardware that are available through the cloud. In order to
implement the SWAP test [36] in real hardware, even though the cost of a swap
gate is high depending on the hardware topology, we could implement it without
making various re-routing of qubits. This topology is supported by the IonQ
device provided by Amazon Braket [38].

$p_{1}$$p_{2}$$p_{3}$$a_{1}$$q_{1}$$q_{2}$$q_{3}$$\dots$$\dots$ Figure 12:
Hardware topology needed for the SWAP Test.

In order to implement the Destructive SWAP test in real hardware, it is
necessary to have a specific hardware topology that is represented in FIGURE
13. Where the first row of the graph encodes one quantum state and the second
encodes the second quantum state, and near-neighbor connection between them is
used for the Destructive SWAP test. There are actual devices that support this
kind of topology, for instance there is the Melbourne chip from IBMQ [39] and
the IonQ device provided by Amazon Braket [38].

$p_{1}$$p_{2}$$p_{3}$$q_{1}$$q_{2}$$q_{3}$$\dots$$\dots$ Figure 13: Hardware
topology needed for the Destructive Swap Test.

It is possible to measure the overlap between two states using the transition
amplitude between the two states [40]:
$\big{|}\langle\psi(\lambda_{k})|\psi(\lambda_{i})\rangle\big{|}^{2}=\big{|}\langle
0|U^{\dagger}_{\psi(\lambda_{k})}U_{\psi(\lambda_{i})}|0\rangle\big{|}^{2}$,
where $U_{\psi(\lambda_{i})}$ is the state generated by the variational
ansatz. The overlapping is obtained by measuring the frequency of obtaining
the $\left|0\right\rangle$ state. This approach leads to no restriction for
the hardware. However, it doubles the depth of the circuit, which exceed the
depth that is suited for NISQ devices. This procedure is represented on FIGURE
14.

Figure 14: Circuit representation of the transition amplitude method [40] for
measuring the overlap between $\left|\psi(\lambda_{i})\right\rangle$ and
$\left|\psi(\lambda_{k})\right\rangle$.

To summarize, if the topology of the hardware supports the FIGURE 12 topology,
using the SWAP test can be an advantage because you would be adding only 3
gates to your circuit, thus not growing the depth to the point that will make
it impossible to run on NISQ devices. If the hardware supports FIGURE 13,
using the destructive SWAP test could be an advantage because you would only
need a constant number of gates added in the end of the circuit, since the
overlap can be easily calculated using classical post-processing. It may
appear that the transition amplitude method is a superior method than the
other methods because it doesn’t imply any restrictions of the topology,
however there is a huge downside that comes with doubling the depth of the
circuit, which can make it worse than the Destructive SWAP test, and the SWAP
test in some hardware topologies.


* [37] K. Riedl.  Leveraging memory effects and gradient information in consensus-based optimization: On global convergence in mean-field law.  arXiv:2211.12184, 2022.
  * [38] M. Sion.  On general minimax theorems.  Pacific Journal of mathematics, 8(1):171–176, 1958.
  * [39] J. Von Neumann.  Zur Theorie der Gesellschaftsspiele.  Mathematische Annalen, 100(1):295–320, 1928.
  * [40] J. Von Neumann and O. Morgenstern.  Theory of games and economic behavior.  In Theory of games and economic behavior. Princeton university press, 2007.