# Movable but unavoidable nodal lines through high-symmetry points in two-
dimensional materials

V Damljanović  Institute of Physics Belgrade, University of Belgrade,
Pregrevica 118, 11080 Belgrade, Serbia<EMAIL_ADDRESS>
###### Abstract

In two-dimensional materials electronic band contacts often give non-trivial
contribution to materials topological properties. Besides band contacts at
high-symmetry points (HSP) in the Brillouin zone (BZ), like those in graphene,
there are nodal lines which form various patterns in the reciprocal space. In
this paper we have found all movable nodal lines, which shape depends on the
model, that pass through HSPs in the presence of time-reversal symmetry. Cases
with and without spin-orbit coupling are included by studying all eighty layer
groups and their double extensions. Eight single and six double groups,
including three symmorphic, necessarily host Dirac and Weyl nodal lines that
extend through the whole BZ, respectively. Our research might be of interest
in designing new materials with interesting physical properties.

  * February 2023

##  1 Introduction

Physics of two-dimensional (2D) materials became one of the focal points of
solid-state science, since the electric field effect has been reported in
atomically thin graphene in 2004. Brillouin zone (BZ) of 2D materials is an
even-dimensional manifold, so topological invariants are used in explaining
properties of 2D materials - such as integer quantum Hall effect. Non-zero
topological charge originates from singularities of electron wave functions in
the reciprocal space, and these singularities appear at BZ points where two or
more electron bands touch [1].

In 2D, bands can touch at isolated points or lines. Band touchings can be
essential [2, 3] or accidental [4]. The former are induced by the crystal and
time-reversal symmetry (TRS), they are pinned to a high symmetry point (HSP)
or line (HSL) and they are intact by symmetry preserving perturbations.
Accidental band contacts change their position in the BZ by perturbations that
keep the symmetry, they cannot be removed if protected by topology, or can be
removed in the opposite case (truly accidental band contacts). So far, the
existence of such movable but unavoidable nodal points (MNP) and lines (MNL)
was almost exclusively linked to non-symmorphic symmetries of a crystal [5, 6,
7, 8, 9]. We show below that MNLs exist in some symmorphic groups too.

Nodal lines are further classified according to dispersion in the directions
perpendicular to the line. For higher order nodal lines, the splitting is of
order higher then one. In 3D non magnetic materials with spin-orbit coupling
(SOC) such nodal lines are of order two and three [10]. Weyl (WNL) and Dirac
nodal lines (DNL) have linear splitting, with two-fold and four-fold
degeneracy on the line, respectively. Existence of WNL or DNL in bulk band
structure leads to interesting physical properties, such as _drumhead_ edge
states, which are nearly dispersionless and hence exhibits strong electron
correlations [11, 12]. Existence of nodal lines in 3D is connected with
particular symmetry element [13, 14, 15, 16, 17, 18], and also analyzed by use
of point groups [19], all 230 space groups [20] and with 3D magnetic groups
[21]. These publications are augmented with exhaustive tables of all possible
quasiparticles in 3D magnetic and non-magnetic crystals [22, 23, 24, 25, 26].

In 2D WNLs are predicted in $\mathrm{GaTeI}$ family [27] and in ferromagnetic
layers [28, 29, 30]. For determination of properties of 2D materials using
symmetry, characters or irreducible (co)representations of layer groups are
necessary. Those are tabulated for layer double little groups [31], for
generators of layer single groups [32] and for maximal symmetry of hexagonal
materials [33]. Recently, full tabulation of layer groups, single and double,
gray and ordinary is published [34]. We used [34] together with tables for
(3D) space groups from Bilbao Crystallographic Server [35, 36] as a starting
point of our research.

Here we present all MNLs that pass through HSPs in layer groups with TRS, with
and without SOC. We show that these MNLs are unavoidable for crystals
belonging to eight layer groups (describing systems without SOC), and six
layer double groups (for crystals with SOC). We give dispersions and effective
Hamiltonians near mentioned HSPs, and have illustrated our theory with few
examples. An interesting curiosity are three symmorphic layer groups that
support MNLs when SOC is included. Finally, we compare our results with those
published in the literature on related topics.

##  2 Method

For $\mathbf{k}_{0}$ being a BZ point (HSP), $\mathbf{q}$ a small wave vector,
$R$ allowed irreducible (co)representation of little group
$G_{\mathbf{k}_{0}}$ of $\mathbf{k}_{0}$ and
$g=\left\\{\widehat{h}\right|\left.\mathbf{\tau}_{\hat{h}}\right\\}$ an
element (in Seitz notation) of $G_{\mathbf{k}_{0}}$, the following formula
holds [37]:

$\widehat{H}\left(\mathbf{k}_{0}+\mathbf{q}\right)=\widehat{R}^{{\dagger}}\left(g\right)\widehat{H}\left(\mathbf{k}_{0}+\widehat{h}^{\prime}\mathbf{q}\right)\widehat{R}\left(g\right),$
(1a) where $\widehat{h}^{\prime}$ is the operator reduction of $\widehat{h}$
to 2D BZ. If, in addition, $-\mathbf{k}_{0}$ belongs to the star of
$\mathbf{k}_{0}$, then TRS (denoted by $\theta$) gives:
$\widehat{H}^{\prime*}\left(\mathbf{q}\right)=\widehat{R}^{{\dagger}}\left(\theta
g_{s}\right)\widehat{H}^{\prime}\left(-\widehat{h}_{s}^{\prime}\mathbf{q}\right)\widehat{R}\left(\theta
g_{s}\right),$ (1b)

where
$\widehat{H}^{\prime}(\mathbf{q})=\widehat{H}(\mathbf{k}_{0}+\mathbf{q})-\widehat{H}(\mathbf{k}_{0})$,
$g_{s}=\left\\{\widehat{h}_{s}\right|\left.\mathbf{\tau}_{\hat{h}_{s}}\right\\}$
and $\widehat{h}_{s}\mathbf{k}_{0}$ differs from $-\mathbf{k}_{0}$ by a
reciprocal lattice vector. By differentiating these relations with respect to
$\mathbf{q}$ at $\mathbf{q}=0$ once, twice,… _etc._ we get relations between
Hamiltonian derivatives, and from there we derive dispersions in the vicinity
of HSP. We have considered Hamiltonians up to and including second (third)
order for case without (with) SOC and we have isolated groups and HSPs that
are not located at HSL and that have line-like solutions for vanishing
splitting terms. In all such cases, two-component Hamiltonians contain only
one Pauli matrix (in addition to the unit matrix) which we choose to be
$\widehat{\sigma}_{1}$. Since $n$-th order Hamiltonian derivatives transform
as $[\Gamma_{2DPV}^{\otimes n}]\otimes R^{*}\otimes R$, and since in all cases
of interest (except for double group $74^{D}$-BZ center),
$\Gamma_{2DPV}^{\otimes 3}=4\Gamma_{2DPV}$ (or $\Gamma_{2DPV}^{\otimes
3}=2\Gamma_{1}+3\Gamma_{2DPV}$, for $74^{D}$-BZ center), we could exclude the
possibility that some higher order derivatives couple to any other Pauli
matrix, than the ones already contained in the linear or quadratic (linear,
quadratic or cubic for $74^{D}$-BZ center) terms in Hamiltonian Taylor
expansion, which would cause energy splitting. $\Gamma_{2DPV}$ denotes 2D
polar-vector representation that acts in 2D reciprocal space, $\Gamma_{1}$ is
the unit representation, while $[\Gamma_{2DPV}^{\otimes n}]$ is the
symmetrised $n$-th power of $\Gamma_{2DPV}$. It turned out that four-component
Hamiltonians, that do not correspond to fortune teller (FT) dispersion [38,
39], do not support MNLs. Since gray layer single and double groups have
multidimensional little group correps of dimensions two and four, our method
exhausts all possibilities.

We now prove that, if MNL exists in the vicinity of a HSP, it must continue
further. For $\mathbf{k}_{0}$ now being touching point of two bands located
near HSP, its little group contains only horizontal reflection plane as a
nontrivial element. Both bands can be approximated with
$E_{1}=p_{1}q_{1}+p_{2}q_{2}$ and $E_{2}=c_{1}q_{1}+c_{2}q_{2}$ respectively,
so that $p_{1}\neq c_{1}$ and $p_{2}\neq c_{2}$ (at $\mathbf{q}=0$ two bands
touch by initial assumption, $\mathbf{k}_{0}$ is now not a HSP but lies close
to it). Now $E_{1}=E_{2}$ has solution
$q_{2}=q_{1}(p_{1}-c_{1})/(c_{2}-p_{2})$ which means that the line continues.
We can repeat the same reasoning for the another point close to
$\mathbf{k}_{0}$ and so on. In other words, MNL can not be interrupted.

Finally, guided by the result described in the previous paragraph, we have
included groups that host FT states, described in [38] and [39] for non-SOC
and SOC cases, respectively, and also experimentally observed [40]. Among
conclusions of previous paragraph is that FT states cannot be the only
features at the Fermi level. This holds no matter if Chern number of FT states
is integer or not.

##  3 Results

In what follows, $c$’s denote real parameters, while $q_{1}$ and $q_{2}$ are
projections of $\mathbf{q}$ to orthonormal basis vectors $\mathbf{e}_{1}$ and
$\mathbf{e}_{2}$, respectively. The orientation of basis vectors is arbitrary
since no little group of HSPs in question contains horizontal rotation axis or
vertical reflection plane, except for groups hosting FT dispersion. However,
depending on the parameters of the model, sometimes it is more convenient to
switch to other orthonormal basis $\mathbf{e}_{1}^{\prime}$,
$\mathbf{e}_{2}^{\prime}$, which will be defined below. In these cases
projections of $\mathbf{q}$ to this new basis are denoted with
$q_{1}^{\prime}$ and $q_{2}^{\prime}$. Finally $V_{l,m}^{(2)}$
($V_{l,m,n}^{(3)}$) are second (third) order derivatives of certain components
of the Hamiltonian. These are real parameters, symmetric under any permutation
of indices, but with no additional symmetry. The superscript $D$ denotes
double group, which is the symmetry of the system with SOC.

Groups and HSPs hosting MNL when SOC is neglected are shown in the Table 1,
those with SOC in the Table 2. For groups $5$, $7$, $36$, $48$, $52$(points
$(0,\pm 1/2)$, $(\pm 1/2,0)$), $4^{D}$, $5^{D}$, $35^{D}$ and $74^{D}$(points
$(0,\pm 1/2)$, $(\pm 1/2,0)$, $\pm(1/2,1/2)$), the Hamiltonian and dispersion
are:

$\widehat{H}=\left(E_{0}+\sum_{l,m=1}^{2}V_{l,m}^{(2)}q_{l}q_{m}\right)\widehat{\sigma}_{0}+\left(c_{1}q_{1}+c_{2}q_{2}+\sum_{l,m,n=1}^{2}V_{l,m,n}^{(3)}q_{l}q_{m}q_{n}\right)\widehat{\sigma}_{1},$
(1ba)
$E_{1,2}=E_{0}+\sum_{l,m=1}^{2}V_{l,m}^{(2)}q_{l}q_{m}\pm|c_{1}q_{1}+c_{2}q_{2}+\sum_{l,m,n=1}^{2}V_{l,m,n}^{(3)}q_{l}q_{m}q_{n}|.$
(1bb) For new basis given by:
$\displaystyle\mathbf{e}_{1}^{\prime}=\frac{c_{1}\mathbf{e}_{1}+c_{2}\mathbf{e}_{2}}{\sqrt{c_{1}^{2}+c_{2}^{2}}},$
(1bc)
$\displaystyle\mathbf{e}_{2}^{\prime}=\frac{-c_{2}\mathbf{e}_{1}+c_{1}\mathbf{e}_{2}}{\sqrt{c_{1}^{2}+c_{2}^{2}}},$
the Hamiltonian and dispersion are:
$\widehat{H}=\left(E_{0}+\sum_{l,m=1}^{2}V_{l,m}^{\prime(2)}q_{l}^{\prime}q_{m}^{\prime}\right)\widehat{\sigma}_{0}+\sqrt{c_{1}^{2}+c_{2}^{2}}\left(q_{1}^{\prime}+c_{0}q_{2}^{\prime
3}\right)\widehat{\sigma}_{1},$ (1bd)
$E_{1,2}=E_{0}+\sum_{l,m=1}^{2}V_{l,m}^{\prime(2)}q_{l}^{\prime}q_{m}^{\prime}\pm\sqrt{c_{1}^{2}+c_{2}^{2}}|q_{1}^{\prime}+c_{0}q_{2}^{\prime
3}|,$ (1be) with
$c_{0}=(c_{1}^{3}V_{222}^{(3)}-3c_{1}^{2}c_{2}V_{122}^{(3)}+3c_{1}c_{2}^{2}V_{112}^{(3)}-c_{2}^{3}V_{111}^{(3)})/(c_{1}^{2}+c_{2}^{2})^{2}$.
Since $q_{1}^{\prime}+c_{0}q_{2}^{\prime 3}=0$ has one solution, one MNL
passes through $\mathbf{q}=0$. In (1bd) and (1be) we have neglected energy
splitting terms that contain products of $q_{1}^{\prime}$ with
$q_{2}^{\prime}$, in comparison to $q_{1}^{\prime}$. Such HSPs were mentioned
in [41] as having linearity rank one.

Table 1:  Layer single groups hosting unavoidable accidental nodal lines
through HSP. The notation for layer and space groups as well as primitive
basis vectors are according to [42] and [43], respectively. The coordinates of
HSPs in the last column, are in the primitive basis
$\left\\{\mathbf{b}_{1},\mathbf{b}_{2}\right\\}$ of the reciprocal 2D lattice.

  * Layer single group | Corresponding space group | Diperiodic plane | $HSP$
---|---|---|---
5 | $p\ 1\ 1\ b$ | 7 | $P\ 1\ c\ 1$ | $C_{s}^{2}$ | $y=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},\pm\frac{1}{2})$
7 | $p\ 1\ 1\ 2/b$ | 13 | $P\ 1\ 2/c\ 1$ | $C_{2h}^{4}$ | $y=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},\pm\frac{1}{2})$
33 | $p\ b\ 2_{1}\ a$ | 29 | $P\ c\ a\ 2_{1}$ | $C_{2v}^{5}$ | $y=0$ | $(\pm\frac{1}{2},\pm\frac{1}{2})$
36 | $c\ m\ 2\ e$ | 39 | $A\ e\ m\ 2$ | $C_{2v}^{15}$ | $x=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
43 | $p\ 2/b\ 2_{1}/a\ 2/a$ | 54 | $P\ 2_{1}/c\ 2/c\ 2/a$ | $D_{2h}^{8}$ | $y=0$ | $(\pm\frac{1}{2},\pm\frac{1}{2})$
45 | $p\ 2_{1}/b\ 2_{1}/m\ 2/a$ | 57 | $P\ 2/b\ 2_{1}/c\ 2_{1}/m$ | $D_{2h}^{11}$ | $x=0$ | $(\pm\frac{1}{2},\pm\frac{1}{2})$
48 | $c\ 2/m\ 2/m\ 2/e$ | 67 | $C\ 2/m\ 2/m\ 2/e$ | $D_{2h}^{21}$ | $z=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
52 | $p\ 4/n$ | 85 | $P\ 4/n$ | $C_{4h}^{3}$ | $z=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
|  |  |  |  |  | $(\pm\frac{1}{2},\pm\frac{1}{2})$

Table 2:  Layer double groups hosting unavoidable accidental nodal lines
through HSP. The notation for layer and space groups as well as primitive
basis vectors are according to [42] and [43], respectively. The coordinates of
HSPs in the last column, are in the primitive basis
$\left\\{\mathbf{b}_{1},\mathbf{b}_{2}\right\\}$ of the reciprocal 2D lattice.

  * Layer double | Corresponding space | Diperiodic |
---|---|---|---
group | double group | plane | $HSP$
$4^{D}$ | $p\ 1\ 1\ m$ | $6^{D}$ | $P\ 1\ m\ 1$ | $C_{s}^{1}$ | $y=0$ | $(0,0)$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
|  |  |  |  |  | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},\pm\frac{1}{2})$
$5^{D}$ | $p\ 1\ 1\ b$ | $7^{D}$ | $P\ 1\ c\ 1$ | $C_{s}^{2}$ | $y=0$ | $(0,0)$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
$29^{D}$ | $p\ b\ 2_{1}\ m$ | $26^{D}$ | $P\ m\ c\ 2_{1}$ | $C_{2v}^{2}$ | $x=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},\pm\frac{1}{2})$
$33^{D}$ | $p\ b\ 2_{1}\ a$ | $29^{D}$ | $P\ c\ a\ 2_{1}$ | $C_{2v}^{5}$ | $y=0$ | $(0,\pm\frac{1}{2})$
$35^{D}$ | $c\ m\ 2\ m$ | $38^{D}$ | $A\ m\ m\ 2$ | $C_{2v}^{14}$ | $x=0$ | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
$74^{D}$ | $p\ \overline{6}$ | $174^{D}$ | $P\ \overline{6}$ | $C_{3h}^{1}$ | $z=0$ | $(0,0)$
|  |  |  |  |  | $(\pm\frac{1}{2},0)$
|  |  |  |  |  | $(0,\pm\frac{1}{2})$
|  |  |  |  |  | $\pm(\frac{1}{2},\frac{1}{2})$

For group $52$ near BZ corners, the following relations hold:

$\widehat{H}=\left(E_{0}+c\mathbf{q}^{2}\right)\widehat{\sigma}_{0}+\left(c_{1}q_{1}^{2}+c_{2}q_{1}q_{2}-c_{1}q_{2}^{2}\right)\widehat{\sigma}_{1},$
(1bca)
$E_{1,2}=E_{0}+c\mathbf{q}^{2}\pm|c_{1}q_{1}^{2}+c_{2}q_{1}q_{2}-c_{1}q_{2}^{2}|.$
(1bcb) In the new basis:
$\displaystyle\mathbf{e}_{1}^{\prime}=\frac{\mathbf{e}_{1}+\lambda\mathbf{e}_{2}}{\sqrt{1+\lambda^{2}}},$
(1bcc)
$\displaystyle\mathbf{e}_{2}^{\prime}=\frac{-\lambda\mathbf{e}_{1}+\mathbf{e}_{2}}{\sqrt{1+\lambda^{2}}},$
where $\lambda=(c_{2}+\sqrt{c_{2}^{2}+4c_{1}^{2}})/(2c_{1})$, the following
holds: $\widehat{H}=\left(E_{0}+c\mathbf{q}^{\prime
2}\right)\widehat{\sigma}_{0}-\sqrt{c_{2}^{2}+4c_{1}^{2}}q_{1}^{\prime}q_{2}^{\prime}\widehat{\sigma}_{1},$
(1bcd) $E_{1,2}=E_{0}+c\mathbf{q}^{\prime
2}\pm\sqrt{c_{2}^{2}+4c_{1}^{2}}|q_{1}^{\prime}q_{2}^{\prime}|.$ (1bce)

Since $q_{1}^{\prime}q_{2}^{\prime}=0$ has two solutions, two MNL intersect at
right angle for $\mathbf{q}=0$.

For group $74^{D}$ near the BZ center ($\mathbf{k}_{0}=0$):

$\widehat{H}=\left(E_{0}+c\mathbf{q}^{2}\right)\widehat{\sigma}_{0}+\left(c_{1}q_{1}(q_{1}^{2}-3q_{2}^{2})+c_{2}q_{2}(q_{2}^{2}-3q_{1}^{2})\right)\widehat{\sigma}_{1},$
(1bcda)
$E_{1,2}=E_{0}+c\mathbf{q}^{2}\pm|c_{1}q_{1}(q_{1}^{2}-3q_{2}^{2})+c_{2}q_{2}(q_{2}^{2}-3q_{1}^{2})|.$
(1bcdb) New basis is given by:
$\displaystyle\mathbf{e}_{1}^{\prime}=\frac{c_{2}\mathbf{e}_{1}+c_{1}B\mathbf{e}_{2}}{\sqrt{c_{2}^{2}+c_{1}^{2}B^{2}}},$
(1bcdc)
$\displaystyle\mathbf{e}_{2}^{\prime}=\frac{-c_{1}B\mathbf{e}_{1}+c_{2}\mathbf{e}_{2}}{\sqrt{c_{2}^{2}+c_{1}^{2}B^{2}}},$
where
$B=1+2\sqrt{1+c_{2}^{2}/c_{1}^{2}}cos[(1/3)arccos(1/\sqrt{1+c_{2}^{2}/c_{1}^{2}})]$
is a real solution of the cubic equation
$(c_{1}/c_{2})^{2}B^{3}-3(c_{1}/c_{2})^{2}B^{2}-3B+1=0$. In the new basis the
Hamiltonian and dispersion are: $\widehat{H}=\left(E_{0}+c\mathbf{q}^{\prime
2}\right)\widehat{\sigma}_{0}-\frac{c_{1}^{2}(B^{2}-2B)-c_{2}^{2}}{\sqrt{c_{2}^{2}+c_{1}^{2}B^{2}}}q_{2}^{\prime}(q_{2}^{\prime
2}-3q_{1}^{\prime 2})\widehat{\sigma}_{1},$ (1bcdd)
$E_{1,2}=E_{0}+c\mathbf{q}^{\prime
2}\pm\frac{c_{1}^{2}(B^{2}-2B)-c_{2}^{2}}{\sqrt{c_{2}^{2}+c_{1}^{2}B^{2}}}|q_{2}^{\prime}(q_{2}^{\prime
2}-3q_{1}^{\prime 2})|.$ (1bcde)

Here we note that $c_{1}^{2}(B^{2}-2B)-c_{2}^{2}>0$. Since
$q_{2}^{\prime}(q_{2}^{\prime 2}-3q_{1}^{\prime 2})=0$ has three solutions,
three MNL intersects at $\mathbf{q}=0$ at angles of $\pi/3$.

Groups 33, 43, 45, $29^{D}$ and $33^{D}$ host fortune teller (FT) dispersion
which is described in more detail in [38, 39]. Here we only give the
dispersion:

$E_{j,l}=j|c_{1}|q_{1}|+lc_{2}|q_{2}||,(j,l\in\left\\{+,-\right\\}),$ (1bcde)

where $c_{1}$ and $c_{2}$ are (only in this formula) positive quantities.
Bands are ordered in the following way: $E_{-,+}\leq E_{-,-}\leq E_{+,-}\leq
E_{+,+}$. It follows that MNLs $E_{-,-}=E_{+,-}$ are given by
$c_{1}q_{1}=c_{2}q_{2}$ and $c_{1}q_{1}=-c_{2}q_{2}$. Two MNLs intersect at
$\mathbf{q}=0$ at an angle that is dependent on $c_{1}$ and $c_{2}$. Those
MNLs arising from FT dispersion are at constant energy only in linear
approximation. When quadratic corrections to the Hamiltonian are accounted
for, the MNL remain unsplit but overall energy shift is quadratic in $q_{||}$
\- the wave vector along MNL.

##  4 Illustrative examples

We first consider the tight-binding model from the $s$-orbitals, without SOC,
on the structure shown in Figure 1a), that belongs to layer group $p112/b$
($7$). Occupied Wyckoff position is $2d$ with $z=1.3$Å. The Hamiltonian and
its eigenvalues are:

$\displaystyle\widehat{H}(\mathbf{k})=$ (1bcdfa)
$\displaystyle\left(\begin{array}[c]{cc}f_{0}+2f_{1}cos(\mathbf{k}\cdot\mathbf{a}_{1})&2e^{-\frac{i}{2}\mathbf{k}\cdot\mathbf{a}_{2}}(f_{2}cos(\frac{1}{2}\mathbf{k}\cdot\mathbf{a}_{2})+f_{3}cos(\mathbf{k}\cdot\mathbf{a}_{1}-\frac{1}{2}\mathbf{k}\cdot\mathbf{a}_{2}))\\\
2e^{\frac{i}{2}\mathbf{k}\cdot\mathbf{a}_{2}}(f_{2}cos(\frac{1}{2}\mathbf{k}\cdot\mathbf{a}_{2})+f_{3}cos(\mathbf{k}\cdot\mathbf{a}_{1}-\frac{1}{2}\mathbf{k}\cdot\mathbf{a}_{2}))&f_{0}+2f_{1}cos(\mathbf{k}\cdot\mathbf{a}_{1})\end{array}\right),$
(1bcdfd)
$E_{1,2}(\mathbf{k})=f_{0}+2f_{1}cos(\mathbf{k}\cdot\mathbf{a}_{1})\pm
2|f_{2}cos(\frac{1}{2}\mathbf{k}\cdot\mathbf{a}_{2})+f_{3}cos(\mathbf{k}\cdot\mathbf{a}_{1}-\frac{1}{2}\mathbf{k}\cdot\mathbf{a}_{2})|.$
(1bcdfe) Here $f_{0}$, $f_{1}$, $f_{2}$ and $f_{3}$ are hopping real
parameters for zeroth, first, second and third neighbors, respectively. The
two energies coincide for:
$\mathbf{k}\cdot\mathbf{a}_{2}\in\left\\{-2arctan\left(\frac{f_{2}+f_{3}cos(\mathbf{k}\cdot\mathbf{a}_{1})}{f_{3}sin(\mathbf{k}\cdot\mathbf{a}_{1})}\right)+2m\pi\right\\},$
(1bcdff)

where $m$ is any integer. Equation (1bcdff) is well defined for all possible
values of hopping parameters. The full band structure for particular values of
parameters is shown in Figure 1b), where band touching lines are clearly
visible. Position of lines in the reciprocal space are shown in Figure 1c). By
varying parameters, the shape of lines changes but the lines always pass
through points $(0,1/2)$ and $(1/2,1/2)$ and equivalent to them, as shown in
Table 1 and Figure 1c), d).

Figure 1:  Tight binding model for $p112/b$: a) - crystal structure (all
nuclei are of the same type, black and gray circles denotes nuclei above and
below the drawing plane, respectively, visualization by VESTA [44]), b) -
electronic band structure for $f_{0}$, $f_{1}$, $f_{2}$ and $f_{3}$ equal to
$4.2$, $1.3$, $0.9$ and $0.5$ in arbitrary units ($a.u.$), respectively (red
arrows denote MNLs), c) - position of MNL in the reciprocal space, d) - same
as c) for $f_{0}$, $f_{1}$, $f_{2}$ and $f_{3}$ equal to $3.6$, $3.1$, $2.7$
and $0.9$ $a.u.$, respectively. In c) and d) numbers in brackets denote
coordinates of special points in the reciprocal basis (blue and black
circles). Mutually equivalent points are denoted with circles of the same
color. Figure 2:  Tight binding model for symmorphic $p11m$: a) - crystal
structure (nuclei of both types are located at one plane, visualization by
VESTA [44]), b) - electronic band structure for $f_{0}$, $f_{1}$, $z_{1}$,
$z_{2}$ and $z_{3}$ equal to $2$, $3$, $4+2i$, $3-6i$ and $1+2i$ $a.u.$,
respectively (red arrows denote MNLs), c) - position of MNL in the reciprocal
space, d) - same as c) for $f_{0}$, $f_{1}$, $z_{1}$, $z_{2}$ and $z_{3}$
equal to $4.2$, $1.3$, $-2.4+1.2i$, $3.4-5.6i$ and $6.1+1.2i$ $a.u.$,
respectively. In c) and d) numbers in brackets denote coordinates of special
points in the reciprocal basis (blue, orange, green and black circles).
Mutually equivalent points are denoted with circles of the same color.

We next consider the tight-binding model on the structure shown in Figure 2a)
that belongs to symmorpic layer (double) group $p11m$ ($4^{D}$). We choose
basis functions
$\left|s_{1}\uparrow\right\rangle,\left|s_{1}\downarrow\right\rangle,\left|s_{2}\uparrow\right\rangle,\left|s_{2}\downarrow\right\rangle$,
where $s_{1},s_{2}$ are $s$-orbitals for atoms $1$ and $2$, respectively,
while $\uparrow$ and $\downarrow$ denote spinors for up and down spin,
respectively. Since horizontal reflection plane is symmetry of the system, the
_up_ and _down_ states are decoupled from each other, so that Hamiltonian is
block-diagonal:

$\widehat{H}=\left(\begin{array}[c]{cc}\widehat{H}_{\uparrow}&\widehat{0}\\\
\widehat{0}&\widehat{H}_{\downarrow}\end{array}\right),$ (1bcdfga) with
$\widehat{H}_{\uparrow}(\mathbf{k})=\widehat{H}_{\downarrow}^{*}(-\mathbf{k})=\left(\begin{array}[c]{cc}f_{0}&z_{1}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}}+z_{2}+z_{3}e^{i\mathbf{k}\cdot(\mathbf{a}_{2}-\mathbf{a_{1}})}\\\
z_{1}^{*}e^{i\mathbf{k}\cdot\mathbf{a}_{1}}+z_{2}^{*}+z_{3}^{*}e^{-i\mathbf{k}\cdot(\mathbf{a}_{2}-\mathbf{a_{1}})}&g_{0}\end{array}\right).$
(1bcdfgb) Here $f_{0}$ and $g_{0}$ are real, on-site hopping parameters for
sites $1$ and $2$, respectively, while $z_{1}$, $z_{2}$ and $z_{3}$ are
complex hoping parameters for first, second and third neighbors, respectively.
The full band structure is:
$E_{2}^{\uparrow,\downarrow}(\mathbf{k})=\frac{f_{0}+g_{0}}{2}+\frac{1}{2}\sqrt{(f_{0}-g_{0})^{2}+4|z_{1}e^{\mp
i\mathbf{k}\cdot\mathbf{a}_{1}}+z_{2}+z_{3}e^{\pm
i\mathbf{k}\cdot(\mathbf{a}_{2}-\mathbf{a}_{1})}|^{2}},$ (1bcdfgc)
$E_{1}^{\uparrow,\downarrow}(\mathbf{k})=\frac{f_{0}+g_{0}}{2}-\frac{1}{2}\sqrt{(f_{0}-g_{0})^{2}+4|z_{1}e^{\mp
i\mathbf{k}\cdot\mathbf{a}_{1}}+z_{2}+z_{3}e^{\pm
i\mathbf{k}\cdot(\mathbf{a}_{2}-\mathbf{a}_{1})}|^{2}}.$ (1bcdfgd) Both nodal
lines $E_{1}^{\uparrow}=E_{1}^{\downarrow}$ and
$E_{2}^{\uparrow}=E_{2}^{\downarrow}$ are given by the equation:
$h_{3}sin(\mathbf{k}\cdot\mathbf{a}_{1})+h_{2}sin(\mathbf{k}\cdot\mathbf{a}_{2})+h_{1}sin[\mathbf{k}\cdot(\mathbf{a}_{2}-\mathbf{a}_{1})]=0,$
(1bcdfge) where $h_{1}=-i(z_{2}^{*}z_{3}-z_{2}z_{3}^{*})$,
$h_{2}=-i(z_{1}^{*}z_{3}-z_{1}z_{3}^{*})$ and
$h_{3}=-i(z_{1}^{*}z_{2}-z_{1}z_{2}^{*})$ are real parameters. Equation
(1bcdfge) has the following solution:
$\mathbf{k}\cdot\mathbf{a}_{2}\in\left\\{\Phi-
arcsin\left(\frac{h_{3}}{h_{1}}sin\Phi\right)+2m\pi,\Phi+arcsin\left(\frac{h_{3}}{h_{1}}sin\Phi\right)+(2m+1)\pi\right\\},$
(1bcdfgf) with $m$ being any integer and
$sin(\Phi)=\frac{h_{1}sin(\mathbf{k}\cdot\mathbf{a}_{1})/h_{3}}{\sqrt{(h_{1}^{2}+h_{2}^{2})/h_{3}^{2}+2h_{1}h_{2}cos(\mathbf{k}\cdot\mathbf{a}_{1})/h_{3}^{2}}},$
(1bcdfgg)
$cos(\Phi)=\frac{h_{1}cos(\mathbf{k}\cdot\mathbf{a}_{1})/h_{3}+h_{2}/h_{3}}{\sqrt{(h_{1}^{2}+h_{2}^{2})/h_{3}^{2}+2h_{1}h_{2}cos(\mathbf{k}\cdot\mathbf{a}_{1})/h_{3}^{2}}}.$
(1bcdfgh)

The argument of $arcsin$ in (1bcdfgf) has modulus less or equal to one except
if $(h_{1}/h_{3})^{2}$ and $(h_{2}/h_{3})^{2}$ are both less then one. In that
case $-1\leq cos(\mathbf{k}\cdot\mathbf{a}_{1})\leq t_{2}$ or $t_{1}\leq
cos(\mathbf{k}\cdot\mathbf{a}_{1})\leq 1$ must hold
($t_{1,2}=-h_{1}h_{2}/h_{3}^{2}\pm\sqrt{(1-h_{1}^{2}/h_{3}^{2})(1-h_{2}^{2}/h_{3}^{2})}$).
Band structure shown in Figure 2b) for particular values of parameters,
confirm existence of MNL. In this case $h_{1}^{2}$ and $h_{2}^{2}$ are both
less then $h_{3}^{2}$, while in the case shown in Figure 2d) this does not
hold, with the clear difference in accidental degeneracy pattern. In any case,
MNLs pass through HSPs as shown in Table 2 and in Figure 2c), d).

##  5 Discussion and conclusions

We compare our results with the ones on related topics published in the
literature. Analysis of all 3D crystallographic point double groups, for
systems with SOC, that shows in which cases movable or unmovable nodal lines
exist is reported [45]. Their results can be used for symmorphic little
groups, since then allowed representations are just a product of point group
irreducible representations with the phase factors that represent
translations. Point groups $\underline{C}_{s}^{D}$ and
$\underline{C}_{3h}^{D}$ hosts one and three MNL, respectively [45]. This is
in accordance with our results in Table 2 (little group for $74^{D}$ in BZ
center is $\underline{C}_{3h}^{D}$, for remaining symmorphic groups it is
$\underline{C}_{s}^{D}$). On the other hand, HSPs for non-symmorphic $29^{D}$
and $33^{D}$ from Table 2 host two intersecting MNL, while [45] predicts one
nodal line along axis of order 2 for $\underline{C}_{2v}^{D}$. This difference
is because $29^{D}$ and $33^{D}$ are non-symmorphic and because band
degeneracy at HSPs hosting FT states is four, while [45] considers two-band
models.

Origin of doubly degenerate nodal lines that connect time-reversal invariant
momenta in 3D crystals with SOC is reported in [46]. List of symmorphic 3D
space groups that host such lines [46] contains all symmorphic corresponding
space groups in Table 2, which support our findings. Compound
$\mathrm{PbTaSe_{2}}$ is a symmorphic 3D metal with SOC that exhibits nodal
lines in experiments [47], as predicted in [46]. To avoid Kramers degeneracy
in the whole BZ, authors of [46] considered non-centrosymmetric 3D crystals.
They show that nodal line through HSP exists if it’s little groups have at
least one impropper rotation (symmetry element with determinant $-1$). In our
cases this is always fulfilled, since at least horizontal plane must be
present in order to avoid band repulsion due to non-crossing rule.

To illustrate their own results, authors of [48] used tight binding model with
SOC on a structure that belongs to group $74^{D}$ ($p\overline{6}$). Three
nodal lines intersects symmetrically at BZ-center and also pass through middle
of BZ edges, thus confirming our prediction. Numerical calculations [49]
suggest that boron bilayer adopts structure with the symmetry $p\overline{6}$.
Since boron is light element, reported band structure is without SOC [49].
Future research could show if it is possible to enhance SOC using bilayer
boron as a starting material.

Group theoretical conditions for semi-Dirac dispersion in non-magnetic 2D
materials with negligible SOC are reported [50]. All groups and HSPs from [50]
are here included in Table 1. Tight-binding model for group $c\ 2/m\ 2/m\ 2/e$
gives line of degeneracy which depends on numerical values of tight-binding
parameters [50], hence it is a MNL as predicted here. Moreover, the dispersion
near $(0,1/2)$ point in [50], matches the one for group $c\ 2/m\ 2/m\ 2/e$
given here in (1be). Similarly, tight-binding model for group $p\ 1\ 1\ b$
published in [51], exhibits MNL across the whole BZ with the same dispersion
near HSPs indicated in Table 1. Lines of accidental degeneracy, noticed in
[51], gained full explanation in the present work. Group $c\ m\ 2\ e$ is
omitted in [50], due to a bug in REPRES program on Bilbao Crystallographic
Server, which is fixed recently. According to [50], authors used REPRES [35]
for finding allowed representations of layer groups via their corresponding
space groups. Finally, one can point out that groups from [50] give semi-Dirac
dispersion (quadratic in one direction of the BZ that split linearly in the
perpendicular direction) while the bands do not form cones. Formation of semi-
Dirac cones requires Hamiltonian of the form:
$c_{1}q_{1}\widehat{\sigma}_{j}+c_{2}q_{2}^{2}\widehat{\sigma}_{l}$ (or, in
more general case:
$c_{1}q_{1}\widehat{\sigma}_{j}+c_{2}q_{2}^{n}\widehat{\sigma}_{l}$), with two
different Pauli matrices ($j\neq l$) and with $n$-natural number greater then
one. We have omitted term with the unit matrix $\widehat{\sigma}_{0}$, since
it does not lead to band splitting. In this case the band splitting is:
$\pm\sqrt{c_{1}^{2}q_{1}^{2}+c_{2}^{2}q_{2}^{4}}$ (in more general case:
$\pm\sqrt{c_{1}^{2}q_{1}^{2}+c_{2}^{2}q_{2}^{2n}}$), with $\mathbf{q}=0$ as
the only touching point. Detail analysis of all HSP of all gray layer single
or double groups does not give any HSP with semi-Dirac cones [41].

In summary we have determined all layer groups that host nodal lines through
high-symmetry BZ points, which shape depends on the model parameters (MNL). In
our analysis we have included both cases without and with SOC in the presence
of TRS. Although we didn’t restrict our search to Weyl or Dirac nodal lines
only, it turned out that all MNL we found are either DNL (in cases without
SOC) or WNL (in cases with SOC). Since little group representations
responsible for MNLs are the only allowed in listed HSPs, the existence of
MNLs is guaranteed by the layer groups from Tables 1 and 2. On the other hand,
it is not guaranteed that MNLs are at constant energy, nor that they lie close
to the Fermi level, although such cases are not necessarily excluded. Our
finding extend the knowledge of band degeneracies in 2D materials and might be
useful for designing 2D materials which sensitivity to external conditions
like temperature, pressure, strain _etc_ is yet to be studied.

Author acknowledge funding by the Ministry of Science, Technological
Development and Innovation of the Republic of Serbia provided by the Institute
of Physics Belgrade, University of Belgrade.

## References

## References

  * [1] Xiaolong Feng, Jiaojiao Zhu, Weikang Wu, and Shengyuan A. Yang.  Two-dimensional topological semimetals.  Chinese Physics B, 30(10):107304, nov 2021.
  * [2] L. P. Bouckaert, R. Smoluchowski, and E. Wigner.  Theory of $\mathrm{Brillouin}$ zones and symmetry properties of wave functions in crystals.  Phys. Rev., 50:58–67, Jul 1936.
  * [3] Conyers Herring.  Effect of time-reversal symmetry on energy bands of crystals.  Phys. Rev., 52:361–365, Aug 1937.
  * [4] Conyers Herring.  Accidental degeneracy in the energy bands of crystals.  Phys. Rev., 52:365–373, Aug 1937.
  * [5] Steve M. Young and Charles L. Kane.  Dirac semimetals in two dimensions.  Phys. Rev. Lett., 115:126803, Sep 2015.
  * [6] Benjamin J. Wieder and C. L. Kane.  Spin-orbit semimetals in the layer groups.  Phys. Rev. B, 94:155108, Oct 2016.
  * [7] J. Zhang, Y.-H. Chan, C.-K. Chiu, M. G. Vergniory, L. M. Schoop, and A. P. Schnyder.  Topological band crossings in hexagonal materials.  Phys. Rev. Materials, 2:074201, Jul 2018.
  * [8] Moritz M. Hirschmann, Andreas Leonhardt, Berkay Kilic, Douglas H. Fabini, and Andreas P. Schnyder.  Symmetry-enforced band crossings in tetragonal materials: Dirac and $\mathrm{Weyl}$ degeneracies on points, lines, and planes.  Phys. Rev. Materials, 5:054202, May 2021.
  * [9] Andreas Leonhardt, Moritz M. Hirschmann, Niclas Heinsdorf, Xianxin Wu, Douglas H. Fabini, and Andreas P. Schnyder.  Symmetry-enforced topological band crossings in orthorhombic crystals: Classification and materials discovery.  Phys. Rev. Materials, 5:124202, Dec 2021.
  * [10] Zhi-Ming Yu, Weikang Wu, Xian-Lei Sheng, Y. X. Zhao, and Shengyuan A. Yang.  Quadratic and cubic nodal lines stabilized by crystalline symmetry.  Phys. Rev. B, 99:121106, Mar 2019.
  * [11] Chen Fang, Hongming Weng, Xi Dai, and Zhong Fang.  Topological nodal line semimetals.  Chinese Physics B, 25(11):117106, nov 2016.
  * [12] Shuo-Ying Yang, Hao Yang, Elena Derunova, Stuart S. P. Parkin, Binghai Yan, and Mazhar N. Ali.  Symmetry demanded topological nodal-line materials.  Advances in Physics: X, 3(1):1414631, 2018.
  * [13] Chen Fang, Yige Chen, Hae-Young Kee, and Liang Fu.  Topological nodal line semimetals with and without spin-orbital coupling.  Phys. Rev. B, 92:081201, Aug 2015.
  * [14] Youngkuk Kim, Benjamin J. Wieder, C. L. Kane, and Andrew M. Rappe.  Dirac line nodes in inversion-symmetric crystals.  Phys. Rev. Lett., 115:036806, Jul 2015.
  * [15] Zihao Gao, Meng Hua, Haijun Zhang, and Xiao Zhang.  Classification of stable $\mathrm{Dirac}$ and $\mathrm{Weyl}$ semimetals with reflection and rotational symmetry.  Phys. Rev. B, 93:205109, May 2016.
  * [16] Bohm-Jung Yang, Troels Arnfred Bojesen, Takahiro Morimoto, and Akira Furusaki.  Topological semimetals protected by off-centered symmetries in nonsymmorphic crystals.  Phys. Rev. B, 95:075135, Feb 2017.
  * [17] Zhida Song, Tiantian Zhang, and Chen Fang.  Diagnosis for nonmagnetic topological semimetals in the absence of spin-orbital coupling.  Phys. Rev. X, 8:031069, Sep 2018.
  * [18] Heqiu Li, Chen Fang, and Kai Sun.  Diagnosis of topological nodal lines with nontrivial monopole charge in the presence of rotation symmetries.  Phys. Rev. B, 100:195308, Nov 2019.
  * [19] Shingo Kobayashi, Youichi Yamakawa, Ai Yamakage, Takumi Inohara, Yoshihiko Okamoto, and Yukio Tanaka.  Crossing-line-node semimetals: General theory and application to rare-earth trihydrides.  Phys. Rev. B, 95:245208, Jun 2017.
  * [20] Lin Wu, Feng Tang, and Xiangang Wan.  Symmetry-enforced band nodes in 230 space groups.  Phys. Rev. B, 104:045107, Jul 2021.
  * [21] Jian Yang, Chen Fang, and Zheng-Xin Liu.  Symmetry-protected nodal points and nodal lines in magnetic materials.  Phys. Rev. B, 103:245141, Jun 2021.
  * [22] Feng Tang and Xiangang Wan.  Exhaustive construction of effective models in 1651 magnetic space groups.  Phys. Rev. B, 104:085137, Aug 2021.
  * [23] Zhi-Ming Yu, Zeying Zhang, Gui-Bin Liu, Weikang Wu, Xiao-Ping Li, Run-Wu Zhang, Shengyuan A. Yang, and Yugui Yao.  Encyclopedia of emergent particles in three-dimensional crystals.  Science Bulletin, 67(4):375–380, 2022.
  * [24] Gui-Bin Liu, Zeying Zhang, Zhi-Ming Yu, Shengyuan A. Yang, and Yugui Yao.  Systematic investigation of emergent particles in type-$\mathrm{III}$ magnetic space groups.  Phys. Rev. B, 105:085117, Feb 2022.
  * [25] Zeying Zhang, Gui-Bin Liu, Zhi-Ming Yu, Shengyuan A. Yang, and Yugui Yao.  Encyclopedia of emergent particles in type-$\mathrm{IV}$ magnetic space groups.  Phys. Rev. B, 105:104426, Mar 2022.
  * [26] Feng Tang and Xiangang Wan.  Complete classification of band nodal structures and massless excitations.  Phys. Rev. B, 105:155156, Apr 2022.
  * [27] Weikang Wu, Yalong Jiao, Si Li, Xian-Lei Sheng, Zhi-Ming Yu, and Shengyuan A. Yang.  Hourglass $\mathrm{Weyl}$ loops in two dimensions: Theory and material realization in monolayer $\mathrm{GaTeI}$ family.  Phys. Rev. Mater., 3:054203, May 2019.
  * [28] Baojie Feng, Run-Wu Zhang, Ya Feng, Botao Fu, Shilong Wu, Koji Miyamoto, Shaolong He, Lan Chen, Kehui Wu, Kenya Shimada, Taichi Okuda, and Yugui Yao.  Discovery of $\mathrm{Weyl}$ nodal lines in a single-layer ferromagnet.  Phys. Rev. Lett., 123:116401, Sep 2019.
  * [29] Lei Jin, Xiaoming Zhang, Tingli He, Weizhen Meng, Xuefang Dai, and Guodong Liu.  Ferromagnetic two-dimensional metal-chlorides $\mathrm{MCl(M=Sc,Y,}$ and $\mathrm{La)}$: Candidates for $\mathrm{Weyl}$ nodal line semimetals with small spin-orbit coupling gaps.  Applied Surface Science, 520:146376, 2020.
  * [30] Feng Zhou, Ying Liu, Minquan Kuang, Peng Wang, Jianhua Wang, Tie Yang, Xiaotian Wang, Zhenxiang Cheng, and Gang Zhang.  Time-reversal-breaking $\mathrm{Weyl}$ nodal lines in two-dimensional $\mathrm{A_{3}C_{2}(A=Ti{,}Zr{,}}$ and $\mathrm{Hf)}$ intrinsically ferromagnetic materials with high $\mathrm{Curie}$ temperature.  Nanoscale, 13:8235–8241, 2021.
  * [31] D. B. Litvin and T. R. Wike.  Character Tables and Compatibility Relations of the Eighty Layer Groups and Seventeen Plane Groups.  Plenum Press, 1991.
  * [32] I Milosevic, B Nikolic, M Damnjanovic, and M Krcmar.  Irreducible representations of diperiodic groups.  Journal of Physics A: Mathematical and General, 31(15):3625, apr 1998.
  * [33] V Damljanović, R Kostić, and R Gajić.  Characters of graphene’s symmetry group $\mathrm{Dg80}$.  Physica Scripta, 2014(T162):014022, sep 2014.
  * [34] Božidar Nikolić, Ivanka Milošević, Tatjana Vuković, Nataša Lazić, Saša Dmitrović, Zoran Popović, and Milan Damnjanović.  Irreducible and site-symmetry-induced representations of single/double ordinary/grey layer groups.  Acta Crystallographica Section A, 78(2):107–114, Mar 2022.
  * [35] Mois I. Aroyo, Asen Kirov, Cesar Capillas, J. M. Perez-Mato, and Hans Wondratschek.  Bilbao Crystallographic Server. II. Representations of crystallographic point groups and space groups.  Acta Crystallographica Section A, 62(2):115–128, Mar 2006.
  * [36] Luis Elcoro, Barry Bradlyn, Zhijun Wang, Maia G. Vergniory, Jennifer Cano, Claudia Felser, B. Andrei Bernevig, Danel Orobengoa, Gemma de la Flor, and Mois I. Aroyo.  Double crystallographic groups and their representations on the $\mathrm{Bilbao}$ $\mathrm{Crystallographic}$ $\mathrm{Server}$.  Journal of Applied Crystallography, 50(5):1457–1477, Oct 2017.
  * [37] J. L. Mañes.  Existence of bulk chiral fermions and crystal symmetry.  Phys. Rev. B, 85:155118, Apr 2012.
  * [38] Vladimir Damljanović, Igor Popov, and Radoš Gajić.  Fortune teller fermions in two-dimensional materials.  Nanoscale, 9:19337–19345, 2017.
  * [39] V Damljanović, N Lazić, A Šolajić, J Pešić, B Nikolić, and M Damnjanović.  Peculiar symmetry-protected electronic dispersions in two-dimensional materials.  Journal of Physics: Condensed Matter, 32(48):485501, sep 2020.
  * [40] Marek Kopciuszyński, Mariusz Krawiec, Lucyna Żurawek, and Ryszard Zdyb.  Experimental evidence of a new class of massless fermions.  Nanoscale Horiz., 5:679–682, 2020.
  * [41] N Lazić, V Damljanović, and M Damnjanović.  Fully linear band crossings at high symmetry points in layers: classification and role of spin–orbit coupling and time reversal.  Journal of Physics A: Mathematical and Theoretical, 55(32):325202, aug 2022.
  * [42] V Kopsky and D. B. Litvin.  International Tables of Crystallography Volume E: Subperiodic Groups.  Kluwer Academic Publishers, 2002.
  * [43] T Hahn.  International Tables of Crystallography Volume A: Space-Group Symmetry.  Springer, 2005.
  * [44] Koichi Momma and Fujio Izumi.  VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data.  Journal of Applied Crystallography, 44(6):1272–1276, Dec 2011.
  * [45] Andy Knoll and Carsten Timm.  Classification of $\mathrm{Weyl}$ points and nodal lines based on magnetic point groups for spin-$\frac{1}{2}$ quasiparticles.  Phys. Rev. B, 105:115109, Mar 2022.
  * [46] Ying-Ming Xie, Xue-Jian Gao, Xiao Yan Xu, Cheng-Ping Zhang, Jin-Xin Hu, Jason Z. Gao, and K. T. Law.  Kramers nodal line metals.  Nature Communications, 12:3064, May 2021.
  * [47] Guang Bian, Tay-Rong Chang, Raman Sankar, Su-Yang Xu, Hao Zheng, Titus Neupert, Ching-Kai Chiu, Shin-Ming Huang, Guoqing Chang, Ilya Belopolski, Daniel S. Sanchez, Madhab Neupane, Nasser Alidoust, Chang Liu, BaoKai Wang, Chi-Cheng Lee, Horng-Tay Jeng, Chenglong Zhang, Zhujun Yuan, Shuang Jia, Arun Bansil, Fangcheng Chou, Hsin Lin, and M. Zahid Hasan.  Topological nodal-line fermions in spin-orbit metal $\mathrm{PbTaSe_{2}}$.  Nature Communications, 7:10556, Feb 2016.
  * [48] Zeying Zhang, Weikang Wu, Gui-Bin Liu, Zhi-Ming Yu, Shengyuan A. Yang, and Yugui Yao.  Encyclopedia of emergent particles in 528 magnetic layer groups and 394 magnetic rod groups.  Phys. Rev. B, 107:075405, Feb 2023.
  * [49] Shao-Gang Xu, Baobing Zheng, Hu Xu, and Xiao-Bao Yang.  Ideal nodal line semimetal in a two-dimensional boron bilayer.  The Journal of Physical Chemistry C, 123(8):4977–4983, 2019.
  * [50] V Damljanović and R Gajić.  Existence of $\mathrm{semi-Dirac}$ cones and symmetry of two-dimensional materials.  Journal of Physics: Condensed Matter, 29(18):185503, mar 2017.
  * [51] V Damljanović.  An example of diperiodic crystal structure with $\mathrm{semi-Dirac}$ electronic dispersion.  Optical and Quantum Electronics, 50(7):272, jun 2018.


# Simulating quantum backflow on a quantum computer

Arseni Goussev and Jaewoo Joo  School of Mathematics and Physics, University
of Portsmouth, Portsmouth PO1 3HF, United Kingdom

###### Abstract

Quantum backflow is a counterintuitive effect in which the probability density
of a free particle moves in the direction opposite to the particle’s momentum.
If the particle is electrically charged, then the effect can be viewed as the
contrast between the direction of electric current and that of the momentum.
To date, there has been no direct experimental observation of quantum
backflow. However, the effect has been simulated numerically (using classical
computers) and optically (using classical light). In this study, we present
the first simulation of quantum backflow using a real quantum computer.

##  I Introduction

Imagine a bead of mass $M$ and electric charge $Q$ constrained to move without
friction along a rigid ring of radius $R$ (Fig. 1).

Figure 1: A particular manifestation of the quantum backflow effect.

Suppose the bead moves freely (i.e., in the absence of any external forces
other than the constraint force) in the counterclockwise direction, so that
the time derivative of the azimuthal angle $\theta$ is non-negative,
$d\theta/dt\geq 0$. The laws of classical physics guarantee that, at every
point of the ring, momentum density and electric current density are both
nonnegative, i.e., both are either zero or point tangentially to the ring in
the direction of increasing $\theta$. Interestingly, the situation can be
drastically different if the motion of the bead is governed by the laws of
quantum mechanics: There are quantum states of the bead for which its
instantaneous momentum and electric current point in the opposite directions
(Fig. 1). This counterintuitive scenario is one particular manifestation of a
broad class of phenomena concerned with the classically forbidden flow of
probability commonly referred to as quantum backflow (QB).

QB was first mentioned in the context of the arrival time problem in quantum
mechanics Allcock (1969); Kijowski (1974). The first in-depth analysis of QB
for a free particle on a line was carried out by Bracken and Melloy Bracken
and Melloy (1994). In particular, they showed that the effect is weak: Only a
small amount of probability – less than 4$\%$ 111The numerical estimate of
this bound, commonly referred to as the Bracken-Melloy constant, has later
been improved Eveson  _et al._ (2005); Penz  _et al._ (2006) and currently
stands at 0.0384517. Recently, the first analytical bound has been reported in
Ref. Trillo  _et al._ (2023), according to which the value of the Bracken-
Melloy constant lies between 0.0315 and 0.0725. – can be transported in the
direction opposite to the particle’s momentum. Interestingly, QB can be
significantly more pronounced for rotational motion. Thus, the classically
forbidden probability transfer can reach values up to approximately 0.116816
in the case of a particle rotating on a ring (Fig. 1) Goussev (2021) and can
be arbitrarily high in two-dimensional systems Strange (2012); Paccoia  _et
al._ (2020); Barbier  _et al._ (2023). Recently, the problem of QB on a ring
has also been considered for the case of a massless Dirac fermion Di Bari  _et
al._ (2023). The literature on QB is substantial, and reviewing it goes beyond
the scope of the present paper. The reader is referred to Ref. Yearsley and
Halliwell (2013) for an elementary introduction and to Refs. Bracken (2021);
Di Bari  _et al._ (2023) for an extensive list of references to more recent
results in the area.

As of today, QB has not been observed experimentally. A promising experimental
scheme that could lead to the observation of QB in Bose-Einstein condensates
was proposed in Ref. Palmero  _et al._ (2013), but, to our knowledge, has
never been realized in practice.

While a direct experimental observation of QB remains an open challenge, there
has been exciting progress in simulating the QB effect using classical light
Eliezer  _et al._ (2020); Daniel  _et al._ (2022); Ghosh  _et al._ . The
simulations utilize the analogy between the dynamics of quantum particles and
the transverse spreading of light beams. In this analogy, the momentum and
probability current of a quantum particle are represented by the transverse
wave vector and Poynting vector of the light beam, respectively. Thus, the
optical experiments reported in Refs. Eliezer  _et al._ (2020); Daniel  _et
al._ (2022) can be regarded as simulating QB for a particle on a line, whereas
the experiment in Ref. Ghosh  _et al._ is the optics counterpart of QB for a
two-dimensional rotational motion.

In this paper, we demonstrate how QB for a particle on a ring (Fig. 1) can be
simulated using a quantum computer. The demonstration involves the following
two steps: (i) we use $N$ qubits to encode a particle state $|\psi\rangle$
comprised of $2^{N}$ eigenstates with non-negative angular momentum, and (ii)
pass $|\psi\rangle$ through a quantum circuit designed to measure the
probability (or electric) current at a given point on the ring. A negative
readout is a manifestation of QB. Theoretical results derived in this paper
are applicable to the case of arbitrary $N$, whereas the actual experimental
demonstration, utilizing the IBM-Q quantum computer IBM Quantum , is performed
for the cases of $N=1$ and $N=2$.

The paper is organized as follows. In Sec. II, we specify the system and
formulate the QB effect. In Sec. III, we introduce a finite-dimensional
probability current operator and derive its decomposition in terms of sums of
tensor products of one-qubit operators. This decomposition is what makes the
following quantum-computer simulation of QB possible. In Sec. IV, we construct
a concrete example, valid for arbitrary $N$, of a quantum state exhibiting QB.
We later use this state in our quantum simulation. Section V presents an
experimental simulation of QB performed on the IBM-Q quantum computer. In Sec.
VI we summarize our work and make concluding remarks. Throughout the paper, we
set the particle mass and electric charge, the ring radius, and the Planck
constant equal to unity, i.e. $M=Q=R=\hbar=1$.

##  II Quantum backflow for circular motion

The Hamiltonian of the particle-on-a-ring system (Fig. 1) is
$\hat{H}=\frac{1}{2}\hat{L}^{2}$, where $\hat{L}=-i\frac{d}{d\theta}$ is the
angular momentum operator. Stationary states $|m\rangle$ that simultaneously
diagonalize $\hat{H}$ and $\hat{L}$, namely

$\hat{H}|m\rangle=\frac{m^{2}}{2}|m\rangle\,,\qquad\hat{L}|m\rangle=m|m\rangle\,,$

are given by

$\langle\theta|m\rangle=\frac{e^{im\theta}}{\sqrt{2\pi}}\qquad(m\in\mathbb{Z})\,.$
(1)

The states are orthonormal,

$\langle m|n\rangle=\int_{-\pi}^{\pi}d\theta\,\langle
m|\theta\rangle\langle\theta|n\rangle=\delta_{mn}\,,$

and form a complete basis.

Now suppose that the particle is in a state $|\psi\rangle$ given by a
superposition of $2^{N}$ stationary states of the lowest possible energy and
non-negative angular momentum:

$|\psi\rangle=\sum_{m=0}^{2^{N}-1}a_{m}|m\rangle\,.$ (2)

The expansion coefficients $a_{m}$ satisfy the normalization condition,

$\sum_{m=0}^{2^{N}-1}|a_{m}|^{2}=1\,.$ (3)

By construction, any angular momentum measurement performed on $|\psi\rangle$
is guaranteed to return a non-negative result. However, the probability (or
electric) current

$J=\left.\real\left\\{\langle\psi|\theta\rangle\left(-i\frac{d}{d\theta}\right)\langle\theta|\psi\rangle\right\\}\right|_{\theta=\theta_{0}}$
(4)

at some fixed point $\theta=\theta_{0}$ on the ring can be negative. This is
the essence of the QB effect.

##  III Current operator

We now introduce a probability current operator and derive its representation
in terms of tensor products of one-qubit gates. Without any loss of
generality, and in order to simplify the calculations to follow, we set
$\theta_{0}=0$. Then, substituting Eq. (2) into Eq. (4) and making use of and
Eq. (1), we obtain

$\displaystyle J$
$\displaystyle=\frac{1}{2\pi}\real\sum_{m,n=0}^{2^{N}-1}a_{m}^{*}na_{n}$
$\displaystyle=\frac{1}{4\pi}\sum_{m,n=0}^{2^{N}-1}a_{m}^{*}(m+n)a_{n}\,.$ (5)

Alternatively, $J$ can be written as

$J=\frac{1}{4\pi}\langle\psi|\hat{\mathcal{J}}_{N}|\psi\rangle\,,$

where the operator

$\hat{\mathcal{J}}_{N}=\sum_{m,n=0}^{2^{N}-1}|m\rangle(m+n)\langle n|$ (6)

represents the (scaled) probability current at $\theta_{0}=0$ for quantum
states in the $2^{N}$-dimensional subspace of the Hilbert space spanned by
$|0\rangle$, $|1\rangle$, $\ldots$, $|2^{N}-1\rangle$.

In the rest of this section, we show how the operator $\hat{\mathcal{J}}_{N}$,
for any $N=1,2,3,\ldots$, can be decomposed into a sum of tensor products of
the following two-dimensional operators: the identity operator $\hat{I}$ and
the Pauli gates $\hat{X}$ and $\hat{Z}$. The main result of this section is
given by Eq. (17). This decomposition is essential for one’s ability to
simulate QB on a quantum computer.

###  III.1 $N=1$ case

We begin by considering the $N=1$ case. In the matrix representation defined
by

$|0\rangle=\begin{pmatrix}1\\\
0\end{pmatrix}\,,\qquad|1\rangle=\begin{pmatrix}0\\\ 1\end{pmatrix}\,,$

the current operator reads

$\hat{\mathcal{J}}_{1}=\begin{pmatrix}0&1\\\ 1&2\end{pmatrix}\,.$ (7)

Clearly, $\hat{\mathcal{J}}_{1}$ can be decomposed as

$\hat{\mathcal{J}}_{1}=\hat{I}+\hat{X}-\hat{Z}\,,$ (8)

where

$\hat{I}=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\,,\quad\hat{X}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\,,\quad\hat{Z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}\,.$

###  III.2 $N=2$ case

We now turn to the $N=2$ case. Writing

$|0\rangle=\begin{pmatrix}1\\\ 0\\\ 0\\\
0\end{pmatrix}\,,\quad|1\rangle=\begin{pmatrix}0\\\ 1\\\ 0\\\
0\end{pmatrix}\,,\quad\text{etc.}\,,$

we have

$\hat{\mathcal{J}}_{2}=\begin{pmatrix}0&1&2&3\\\ 1&2&3&4\\\ 2&3&4&5\\\
3&4&5&6\end{pmatrix}\,.$

The last matrix can be rewritten as follows:

$\hat{\mathcal{J}}_{2}=\begin{pmatrix}\hat{\mathcal{J}}_{1}&\hat{\mathcal{J}}_{1}+2\hat{C}_{1}\\\
\hat{\mathcal{J}}_{1}+2\hat{C}_{1}&\hat{\mathcal{J}}_{1}+4\hat{C}_{1}\end{pmatrix}\,,$
(9)

where $\hat{\mathcal{J}}_{1}$ is given by Eqs. (7) and (8), and

$\hat{C}_{1}=\begin{pmatrix}1&1\\\ 1&1\end{pmatrix}\,.$

Then,

$\displaystyle\hat{\mathcal{J}}_{2}$ $\displaystyle=\begin{pmatrix}1&1\\\
1&1\end{pmatrix}\otimes\hat{\mathcal{J}}_{1}+\begin{pmatrix}0&2\\\
2&4\end{pmatrix}\otimes\hat{C}_{1}$
$\displaystyle=\hat{C}_{1}\otimes\hat{\mathcal{J}}_{1}+2\hat{\mathcal{J}}_{1}\otimes\hat{C}_{1}\,.$
(10)

Noticing that

$\hat{\mathcal{J}}_{1}=\hat{C}_{1}-\hat{Z}\,,$ (11)

we find

$\hat{\mathcal{J}}_{2}=3\hat{C}_{1}\otimes\hat{C}_{1}-\hat{C}_{1}\otimes\hat{Z}-2\hat{Z}\otimes\hat{C}_{1}\,.$
(12)

Finally, using

$\hat{C}_{1}=\hat{I}+\hat{X}\,,$ (13)

we arrive at the following decomposition:

$\hat{\mathcal{J}}_{2}=3(\hat{I}+\hat{X})\otimes(\hat{I}+\hat{X})-(\hat{I}+\hat{X})\otimes\hat{Z}-2\hat{Z}\otimes(\hat{I}+\hat{X})\,.$
(14)

###  III.3 General case

We now generalize the method of Sec. III.2 to construct an explicit
decomposition of $\hat{\mathcal{J}}_{N}$, for arbitrary $N$. To this end, we
first write $\hat{J}_{N}$, defined by Eq. (6), in the matrix form. Following
the convention adopted in Secs. III.1 and III.2, we take $|m\rangle$ to be
represented by the $2^{N}$-dimensional column vector with the $n^{\text{th}}$
element equal to $\delta_{mn}$. Then,

$\hat{\mathcal{J}}_{N}=\begin{pmatrix}0&1&2&\cdots&2^{N}-1\\\
1&2&3&\cdots&2^{N}\\\ 2&3&4&\cdots&2^{N}+1\\\
\cdots&\cdots&\cdots&\cdots&\cdots\\\
2^{N}-1&2^{N}&2^{N}+1&\cdots&2^{N+1}-2\end{pmatrix}\,.$

It is easy to see that the last matrix can be written as [cf. Eq. (9)]

$\hat{\mathcal{J}}_{N}=\begin{pmatrix}\hat{\mathcal{J}}_{N-1}\;&\hat{\mathcal{J}}_{N-1}+2^{N-1}\hat{C}_{N-1}\\\\[5.69046pt]
\hat{\mathcal{J}}_{N-1}+2^{N-1}\hat{C}_{N-1}\;&\hat{\mathcal{J}}_{N-1}+2^{N}\hat{C}_{N-1}\end{pmatrix}\,,$
(15)

where $\hat{C}_{N-1}$ is the $2^{N-1}\times 2^{N-1}$ matrix of ones. It
follows from Eq. (15) that

$\hat{\mathcal{J}}_{N}=\hat{C}_{1}\otimes\hat{\mathcal{J}}_{N-1}+2^{N-1}\hat{\mathcal{J}}_{1}\otimes\hat{C}_{N-1}\,.$
(16)

Substituting

$\hat{C}_{N-1}=\underbrace{\hat{C}_{1}\otimes\hat{C}_{1}\otimes\cdots\otimes\hat{C}_{1}}_{N-1\text{
times}}=\hat{C}_{1}^{\otimes(N-1)}$

into Eq. (16), we get the following recurrence relation [cf. Eq. (10)]:

$\hat{\mathcal{J}}_{N}=\hat{C}_{1}\otimes\hat{\mathcal{J}}_{N-1}+2^{N-1}\hat{\mathcal{J}}_{1}\otimes\hat{C}_{1}^{\otimes(N-1)}\,.$

It is straightforward to verify that the solution to this recurrence relation
is [cf. Eq. (10)]

$\hat{\mathcal{J}}_{N}=\sum_{n=0}^{N-1}2^{n}\hat{C}_{1}^{\otimes(N-1-n)}\otimes\hat{\mathcal{J}}_{1}\otimes\hat{C}_{1}^{\otimes
n}\,.$

Substituting Eq. (11) into the last expression, and using
$\sum_{n=0}^{N-1}2^{n}=2^{N}-1$, we obtain [cf. Eq. (12)]

$\hat{\mathcal{J}}_{N}=(2^{N}-1)\hat{C}_{1}^{\otimes
N}-\sum_{n=0}^{N-1}2^{n}\hat{C}_{1}^{\otimes(N-1-n)}\otimes\hat{Z}\otimes\hat{C}_{1}^{\otimes
n}\,.$

Finally, using Eq. (13), we arrive at the following explicit decomposition of
$\hat{\mathcal{J}}_{N}$ [cf. Eq. (14)]:

$\begin{split}&\hat{\mathcal{J}}_{N}=(2^{N}-1)(\hat{I}+\hat{X})^{\otimes N}\\\
&-\sum_{n=0}^{N-1}2^{n}(\hat{I}+\hat{X})^{\otimes(N-1-n)}\otimes\hat{Z}\otimes(\hat{I}+\hat{X})^{\otimes
n}\,.\end{split}$ (17)

Equation (17) constitutes the main result of this section.

##  IV Example of a backflowing state

Decomposition (17) allows us to devise a quantum computing circuit for
measuring the probability current at a fixed point on the ring. What we also
need for a quantum simulation of QB is a quantum state that would give rise to
a substantially negative probability current. In this section, we present an
explicit example of such a state.

Consider the state $|\psi\rangle$ defined by Eq. (2) with

$a_{m}=\frac{3m-2(2^{N}-1)}{\sqrt{(2^{N+1}+1)(2^{N}-1)2^{N-1}}}\,.$ (18)

The state is normalized. Indeed, using identities
$\sum_{m=0}^{2^{N}-1}m=(2^{N}-1)2^{N-1}$ and
$\sum_{m=0}^{2^{N}-1}m^{2}=\frac{1}{3}(2^{N+1}-1)(2^{N}-1)2^{N-1}$, it is
straightforward to verify that the normalization condition (3) is fulfilled.

Now, substituting Eq. (18) into Eq. (5) and performing a straightforward
calculation, taking into account that
$\sum_{m=0}^{2^{N}-1}m^{3}=2^{2N-2}(2^{N}-1)^{2}$, we find

$J=-\frac{1}{4\pi}\frac{2^{N}(2^{N}-1)}{2^{N+1}+1}\,.$

Clearly $J<0$, for all $N\geq 1$, meaning that $|\psi\rangle$, defined by Eq.
(18) exhibits QB. In particular,

$J=-\frac{1}{10\pi}\simeq-0.031831\quad\text{for}\quad N=1$ (19)

and

$J=-\frac{3}{9\pi}\simeq-0.106103\quad\text{for}\quad N=2\,.$ (20)

It is also interesting to note that, for the state defined by Eq. (18),
$J\to-\infty$ as $N\to\infty$. This example shows that, just as in the
particle-on-a-line case Bracken and Melloy (1994), the instantaneous
probability current for non-negative angular momentum states in a ring is
unbounded from below.

##  V Implementation of quantum backflow on a quantum computer

Equipped with the decomposition (17) and the explicit backflowing state
example (18), we proceed to simulating QB on the IBM-Q quantum computer IBM
Quantum .

Figure 2: Quantum circuits for computing
$\langle\hat{\cal{V}}_{k}\rangle\equiv\langle\psi|\hat{\cal
V}_{k}|\psi\rangle$ for (a) $N=1$ and (b) $N=2$.

The quantum circuits used in our simulation are schematically shown in Figs.
2(a) and 2(b) for the cases of $N=1$ and $N=2$, respectively. The angular
momentum states $|m\rangle$, with $m=0,\ldots,2^{N}-1$, of the particle on the
ring are represented by system qubits that are denoted by $S$ in the $N=1$
case, and by $S_{1}$ and $S_{2}$ in the $N=2$ case. More specifically, we have

$|0\rangle=|0\rangle_{S}\quad\text{and}\quad|1\rangle=|1\rangle_{S}$

for $N=1$, and

$\displaystyle\begin{split}&|0\rangle=|0\rangle_{S_{1}}|0\rangle_{S_{2}}=|00\rangle\\\
&|1\rangle=|0\rangle_{S_{1}}|1\rangle_{S_{2}}=|01\rangle\\\
&|2\rangle=|1\rangle_{S_{1}}|0\rangle_{S_{2}}=|10\rangle\\\
&|3\rangle=|1\rangle_{S_{1}}|1\rangle_{S_{2}}=|11\rangle\end{split}$ (21)

for $N=2$.

The simulation consists of the following four stages. First, we prepare the
system qubit(s) in state $|0\rangle_{S}$ for $N=1$, and in states
$|0\rangle_{S_{1}}$ and $|0\rangle_{S_{2}}$ for $N=2$. Second, we apply a set
of gates to the system quibits to create the target input state
$|\psi\rangle$, Eq. (2), with the expansion coefficients given by Eq. (18).
Third, we act with single-quibit gates $\hat{U}_{k}$ for $N=1$, and with
$\hat{U}_{k1}$ and $\hat{U}_{k2}$ for $N=2$, on $|\psi\rangle$ to obtain a new
state $|\psi^{2}\rangle$. The gates $\hat{U}_{k}$, $\hat{U}_{k1}$, and
$\hat{U}_{k2}$ are to be defined below. Finally, the qubits of
$|\psi^{2}\rangle$ are measured in the $\hat{Z}$ basis and the measurement
outcome probabilities are used to compute the expectation value of operators
$\hat{\cal V}_{k}$. As we explain in detail below, the operators $\hat{\cal
V}_{k}$ are elementary $N$-quibit gates allowing one to represent the current
operator as

$\hat{\cal J}_{N}=\lambda_{0}\hat{I}^{\otimes
N}+\sum_{k=1}^{K_{N}}\lambda_{k}\hat{\cal V}_{k}$ (22)

with $K_{1}=2$, $K_{2}=7$, and $\lambda_{k}$’s being some real numbers.

We now provide further details of the outlined simulation procedure and
present the experimental results, separately in the $N=1$ and $N=2$ case.

###  V.1 $N=1$ case

Here we detail the experimental procedure, schematically illustrated in Fig.
2(a). Hereinafter, we use the notation
$\langle\cdot\rangle=\langle\psi|\cdot|\psi\rangle$, with $|\psi\rangle$ given
by Eqs. (2) and (18).

Let us rewrite Eq. (8) as $\hat{\cal
J}_{1}=\hat{I}+\hat{\cal{V}}_{1}-\hat{\cal{V}}_{2}$ with
$\hat{\cal{V}}_{1}=\hat{X}$ and $\hat{\cal{V}}_{2}=\hat{Z}$. This
representation corresponds to Eq. (22) with $N=1$, $K_{1}=2$,
$\lambda_{0}=\lambda_{1}=1$ and $\lambda_{2}=-1$. Consequently, we have

$\langle\hat{\cal
J}_{1}\rangle=1+\langle\hat{\cal{V}}_{1}\rangle-\langle\hat{\cal{V}}_{2}\rangle\,.$
(23)

The following procedure allows us to experimentally measure
$\langle\hat{\cal{V}}_{1}\rangle$ and $\langle\hat{\cal{V}}_{2}\rangle$.

We first initialize the system qubit, $S$, in state $|0\rangle_{S}$. We then
transform $|0\rangle_{S}$ into the desired backflowing state

$\displaystyle|\psi\rangle$
$\displaystyle=\frac{1}{\sqrt{5}}\left(-2|0\rangle_{S}+|1\rangle_{S}\right)$
$\displaystyle\simeq-0.894427|0\rangle_{S}+0.447214|1\rangle_{S}\,,$

whose expansion coefficients are obtained from Eq. (18) by setting $N=1$ and
$m=0,1$. The transformation is achieved by acting on $|0\rangle_{S}$ with the
$Y$-rotation operator
$\hat{R}_{y}(\alpha)\equiv\exp\left(-i\frac{\alpha}{2}\hat{Y}\right)$:

$|\psi\rangle=\hat{R}_{y}(\alpha)|0\rangle_{S}$

with $\alpha\simeq 5.35589$. Next, we apply gate $\hat{U}_{k}$ to the system
qubit to obtain

$|\psi^{2}\rangle=\hat{U}_{k}|\psi\rangle\,.$ (24)

Here, $\hat{U}_{1}$ is the Hadamard gate,

$\hat{U}_{1}=\hat{H}_{d}\equiv\frac{1}{\sqrt{2}}(\hat{X}+\hat{Z})\,,$

and $\hat{U}_{2}$ is the identity operator,

$\hat{U}_{2}=\hat{I}\,.$

This choice of the transformations $\hat{U}_{1}$ and $\hat{U}_{2}$ allows us
to experimentally determine the sought expectation values
$\langle\hat{\cal{V}}_{1}\rangle$ and $\langle\hat{\cal{V}}_{2}\rangle$ by
measuring $|\psi^{2}\rangle$ in the $\hat{Z}$ basis. This works as follows.

The probabilities that the measurement will collapse $|\psi^{2}\rangle$ onto
$|0\rangle_{S}$ and $|1\rangle_{S}$ are given by

$\displaystyle
P_{0}=\langle\psi^{2}|0\rangle_{S}\langle{0}|\psi^{2}\rangle\,,$
$\displaystyle
P_{1}=\langle\psi^{2}|1\rangle_{S}\langle{1}|\psi^{2}\rangle\,,$

respectively, and so

$P_{0}-P_{1}=\langle\psi^{2}|\hat{Z}|\psi^{2}\rangle\,.$

For $k=1$, we have

$|\psi^{2}\rangle=\hat{U}_{1}|\psi\rangle=\hat{H}_{d}|\psi\rangle\,,$

which means that

$P_{0}-P_{1}=\langle\psi|\hat{H}_{d}\hat{Z}\hat{H}_{d}|\psi\rangle=\langle\psi|\hat{X}|\psi\rangle=\langle\hat{\cal{V}}_{1}\rangle\,,$

where we have used the identity $\hat{X}=\hat{H}_{d}\hat{Z}\hat{H}_{d}$. For
$k=2$, we have

$|\psi^{2}\rangle=\hat{U}_{2}|\psi\rangle=|\psi\rangle\,,$

implying that

$P_{0}-P_{1}=\langle\psi|\hat{Z}|\psi\rangle=\langle\hat{\cal{V}}_{2}\rangle\,.$

This is how, by measuring the difference $P_{0}-P_{1}$ in the quantum circuit
shown in Fig. 2(a), we can evaluate the expectation values
$\langle\hat{\cal{V}}_{k}\rangle$ experimentally.

We used the IBM-Q machine named $ibmq_{-}lima$ IBM Quantum  to perform 8000
projective measurements for each $\langle\hat{\cal{V}}_{k}\rangle$. The
measurement results are $P_{0}\simeq 0.122875$ and $P_{1}\simeq 0.877125$ for
$k=1$, implying $\langle\hat{\cal V}_{1}\rangle\simeq-0.75425$, and
$P_{0}\simeq 0.789625$, $P_{1}\simeq 0.210375$ for $k=2$, implying
$\langle\hat{\cal V}_{2}\rangle\simeq 0.57925$. Thus, in view of Eq. (23), the
experimentally obtained value of probability current in the $N=1$ case is

$J={1\over 4\pi}\langle\hat{\cal J}_{1}\rangle\simeq-0.026539\,.$

The experimental value is negative, signifying QB, and is reasonably close to
the theoretical one, given by Eq. (19).

###  V.2 $N=2$ case

The quantum circuit used for simulating QB in the $N=2$ case is shown in Fig.
2(b). It utilizes two system qubits, $S_{1}$ and $S_{2}$, and seven two-qubit
operators $\hat{\cal{V}}_{k}$, with $k=1,\ldots,7$. According to Eq. (14), the
(scaled) probability current operator $\hat{\cal J}_{2}$ can be represented by
Eq. (22) with $N=2$, $K_{2}=7$, $\hat{\cal V}_{1}=\hat{I}\otimes\hat{X}$,
$\hat{\cal V}_{2}=\hat{X}\otimes\hat{I}$, $\hat{\cal
V}_{3}=\hat{Z}\otimes\hat{I}$, $\hat{\cal V}_{4}=\hat{I}\otimes\hat{Z}$,
$\hat{\cal V}_{5}=\hat{X}\otimes\hat{X}$ $\hat{\cal
V}_{6}=\hat{Z}\otimes\hat{X}$, $\hat{\cal V}_{7}=\hat{X}\otimes\hat{Z}$, and
$\lambda_{0}=\lambda_{1}=\lambda_{2}=\lambda_{5}=3$,
$\lambda_{3}=\lambda_{6}=-2$, $\lambda_{4}=\lambda_{7}=-1$. Thus, the
expectation value of $\hat{\cal J}_{2}$ reads

$\langle\hat{\cal J}_{2}\rangle=3+3\langle\hat{\cal
V}_{1}\rangle+3\langle\hat{\cal V}_{2}\rangle-2\langle\hat{\cal
V}_{3}\rangle-\langle\hat{\cal V}_{4}\rangle+3\langle\hat{\cal
V}_{5}\rangle-2\langle\hat{\cal V}_{6}\rangle-\langle\hat{\cal
V}_{7}\rangle\,.$ (25)

The expectation values $\langle\hat{\cal{V}}_{k}\rangle$, $k=1,\ldots,7$, are
measured as follows.

We begin by initializing the system qubits in the $|00\rangle$ state [see Eq.
(21)]. Then, we apply to $|00\rangle$ a sequence of one- and two-qubit
transformations aiming to generate the following (entangled) backflowing
state:

$\displaystyle|\psi\rangle$
$\displaystyle={1\over\sqrt{6}}(-2|00\rangle-|01\rangle+|11\rangle)$
$\displaystyle\simeq-0.816497|00\rangle-0.408248|01\rangle+0.408248|11\rangle\,.$

The expansion coefficients of $|\psi\rangle$ are obtained from Eq. (18) by
setting $N=2$ and $m=0,1,2,3$. In particular, the overlap probabilities
between $|\psi\rangle$ and the basis states $|00\rangle$, $|01\rangle$,
$|10\rangle$, and $|11\rangle$ are $P_{00}\simeq 0.66667$,
$P_{01}=P_{11}\simeq 0.166667$, and $P_{10}=0$, respectively.

Experimentally, the transformation from $|00\rangle$ to $|\psi\rangle$ is
performed via successive application of three operators Lubasch  _et al._
(2020); Joo and Moon : an $\hat{R}_{y}$ gate, a CNOT gate $CX_{S_{1},S_{2}}$,
and a controlled-$\hat{R}_{y}$ gate. Thus, as illustrated in Fig. 2(b),

$|\psi\rangle=CR_{S_{2},S_{1}}(\alpha_{1})\,CX_{S_{1},S_{2}}\left(\hat{R}_{y}(\alpha_{0})\otimes\hat{I}\right)|00\rangle$

with $\alpha_{0}\simeq 7.51414$ and $\alpha_{1}\simeq 4.7124$. In the actual
experiment, carried out on $ibmq_{-}lima$ IBM Quantum , the controlled-
rotation gate $CR_{S_{2},S_{1}}(\alpha_{1})$ was realized by
$\hat{R}_{y}(\pm\alpha_{1}/2)$ and two CNOT gates Nielsen and Chuang (2000).
The accuracy of the experimental state preparation can be characterized by
following probabilities of the overlap between $|\psi\rangle$ and the two-
qubit basis states: $P_{00}\simeq 0.66125$, $P_{01}\simeq 0.1375$,
$P_{10}\simeq 0.03325$, and $P_{11}\simeq 0.168$, which are to be compared
against the corresponding theoretical values stated above.

Once the backflowing state $|\psi\rangle$ has been prepared, the rest of the
simulation follows the same steps as in the $N=1$ case. The only modification
is that the transformation from $|\psi\rangle$ to $|\psi^{2}\rangle$ [see Fig.
2(b)] is performed by means of two single-qubit gates, $\hat{U}_{k1}$ and
$\hat{U}_{k2}$:

$|\psi^{2}\rangle=\hat{U}_{k1}\otimes\hat{U}_{k2}|\psi\rangle$

[cf. Eq. (24) in the $N=1$ case]. The single-qubit gate pairs
$\\{\hat{U}_{k1},\hat{U}_{k2}\\}$ are $\\{\hat{I},\hat{H}_{d}\\}$ for $k=1$
and $k=6$, $\\{\hat{H}_{d},\hat{I}\\}$ for $k=2$ and $k=7$,
$\\{\hat{I},\hat{I}\\}$ for $k=3$ and $k=4$, and
$\\{\hat{H}_{d},\hat{H}_{d}\\}$ for $k=5$.

Each expectation value $\langle\hat{\cal{V}}_{k}\rangle$, for $k=1,\ldots,7$,
was obtained by performing 4000 two-qubit independent measurements on the
$ibmq_{-}lima$ machine IBM Quantum . The expectation value results from the
experiments are $\langle\hat{\cal V}_{1}\rangle\simeq 0.6215$,
$\langle\hat{\cal V}_{2}\rangle\simeq-0.339$, $\langle\hat{\cal
V}_{3}\rangle\simeq 0.5975$, $\langle\hat{\cal V}_{4}\rangle\simeq 0.389$,
$\langle\hat{\cal V}_{5}\rangle\simeq-0.7395$, $\langle\hat{\cal
V}_{6}\rangle\simeq 0.5015$, and $\langle\hat{\cal V}_{7}\rangle\simeq 0.229$.
Thus, in view of Eq. (25), the experimentally obtained value of probability
current in the $N=2$ case is

$J={1\over 4\pi}\langle\hat{\cal J}_{2}\rangle\simeq-0.0944585\,.$

The experimental value is negative, signifying QB, and is very close to the
theoretical one, given by Eq. (20).

##  VI Conclusion

In this paper, we report a quantum simulation of QB within a circular
geometry. The system under consideration – a particle moving freely in a
circular ring – is well-suited for the simulation, as its discrete angular
momentum eigenstates can be effectively modeled using qubits. More
specifically, $N$-qubit circuits can be employed to simulate QB for quantum
states that are superpositions of $2^{N}$ angular momentum eigenstates. We
explicitly design such circuits for the cases when $N=1$ and $N=2$, and
subsequently implement them on the IBM-Q quantum computer. Our quantum
simulations demonstrate negative probability current for quantum states
comprised solely of non-negative angular momentum states, thereby confirming
the presence of QB. The simulated probability current values are reasonably
close to the corresponding theoretical values.

The quantum simulations presented in this paper have thus far been limited to
one and two qubits, i.e., to the cases of $N=1$ and $N=2$, due to us having
access only to noisy quantum devices. Nonetheless, we have established a
comprehensive theoretical framework for simulating QB with an arbitrary number
of qubits (i.e., for any value of $N$). To elaborate, a quantum simulation of
QB involves two main stages: (i) the preparation of a backflowing state, and
(ii) the subsequent measurement of probability current. Concerning (i), we
have devised a specific example of a backflowing state valid for arbitrary
$N$; this example is detailed by Eqs. (2) and (18). As for (ii), we have
derived a universally applicable decomposition of the probability current
operator, given by Eq. (17), which can be readily employed to construct a
current-measuring circuit for arbitrary $N$. In light of these findings, our
theoretical results provide a clear pathway for simulating QB using any number
of qubits.

The exploration of novel applications for quantum computing has been an active
area of research Gerritsma  _et al._ (2010); Somma ; Bauer  _et al._ . The
study reported in the present paper expands the domain of potential
applications of quantum computing and uncovers new connections between
different areas of quantum physics.

###### Acknowledgements.

JJ acknowledges the support from the Basic Science Research Program through
the National Research Foundation of Korea (NRF) funded by the Ministry of
Education, Science and Technology (NRF-2021M3H3A1038085) and the Institute for
Information & communications Technology Promotion (IITP) grant funded by the
Korea government(MSIP) (No. 2019-0-00003, Research and Development of Core
technologies for Programming, Running, Implementing and Validating of Fault-
Tolerant Quantum Computing System). We acknowledge the utilization of IBM
Quantum services in conducting 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.

## References

  * Allcock (1969) G. R. Allcock, “The time of arrival in quantum mechanics III. The measurement ensemble,” Ann. Phys. (N. Y). 53, 311 (1969).
  * Kijowski (1974) J. Kijowski, “On the time operator in quantum mechanics and the Heisenberg uncertainty relation for energy and time,” Rep. Math. Phys. 6, 361 (1974).
  * Bracken and Melloy (1994) A. J. Bracken and G. F. Melloy, “Probability backflow and a new dimensionless quantum number,” J. Phys. A: Math. Gen. 27, 2197 (1994).
  * Note (1) The numerical estimate of this bound, commonly referred to as the Bracken-Melloy constant, has later been improved Eveson  _et al._ (2005); Penz  _et al._ (2006) and currently stands at 0.0384517. Recently, the first analytical bound has been reported in Ref. Trillo  _et al._ (2023), according to which the value of the Bracken-Melloy constant lies between 0.0315 and 0.0725.
  * Goussev (2021) A. Goussev, “Quantum backflow in a ring,” Phys. Rev. A 103, 022217 (2021).
  * Strange (2012) P. Strange, “Large quantum probability backflow and the azimuthal angle–angular momentum uncertainty relation for an electron in a constant magnetic field,” Eur. J. Phys. 33, 1147 (2012).
  * Paccoia  _et al._ (2020) V. D. Paccoia, O. Panella, and P. Roy, “Angular momentum quantum backflow in the noncommutative plane,” Phys. Rev. A 102, 062218 (2020).
  * Barbier  _et al._ (2023) M. Barbier, A. Goussev, and S. C. L. Srivastava, “Unbounded quantum backflow in two dimensions,” Phys. Rev. A 107, 032204 (2023).
  * Di Bari  _et al._ (2023) L. Di Bari, V. D. Paccoia, O. Panella, and P. Roy, “Quantum backflow for a massless Dirac fermion on a ring,” Phys. Lett. A 474, 128831 (2023).
  * Yearsley and Halliwell (2013) J. M. Yearsley and J. J. Halliwell, “An introduction to the quantum backflow effect,” J. Phys. Conf. Ser. 442, 012055 (2013).
  * Bracken (2021) A. J. Bracken, “Probability flow for a free particle: new quantum effects,” Phys. Scr. 96, 045201 (2021).
  * Palmero  _et al._ (2013) M. Palmero, E. Torrontegui, J. G. Muga, and M. Modugno, “Detecting quantum backflow by the density of a Bose-Einstein condensate,” Phys. Rev. A 87, 053618 (2013).
  * Eliezer  _et al._ (2020) Y. Eliezer, T. Zacharias, and A. Bahabad, “Observation of optical backflow,” Optica 7, 72 (2020).
  * Daniel  _et al._ (2022) A. Daniel, B. Ghosh, B. Gorzkowski, and R. Lapkiewicz, “Demonstrating backflow in classical two beams’ interference,” New J. Phys. 24, 123011 (2022).
  * (15) B. Ghosh, A. Daniel, B. Gorzkowski, and R. Lapkiewicz, “Azimuthal backflow in light carrying orbital angular momentum,” arXiv:2304.13124 .
  * (16) IBM Quantum, https://quantum-computing.ibm.com.
  * Lubasch  _et al._ (2020) M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, “Variational quantum algorithms for nonlinear problems,” Phys. Rev. A 101, 010301(R) (2020).
  * (18) J. Joo and H. Moon, “Quantum variational PDE solver with machine learning,” arXiv:2109.09216 .
  * Nielsen and Chuang (2000) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, Cambridge, 2000).
  * Gerritsma  _et al._ (2010) R. Gerritsma, G. Kirchmair, F. Zähringer, E. Solano, R. Blatt, and C. F. Roos, “Quantum simulation of the Dirac equation,” Nature 463, 68 (2010).
  * (21) R. D. Somma, “Quantum simulations of one dimensional quantum systems,” arXiv:1503.06319 .
  * (22) C. W. Bauer _et al._ , “Quantum Simulation for High Energy Physics,” arXiv:2204.03381 .
  * Eveson  _et al._ (2005) S. P. Eveson, C. J. Fewster, and R. Verch, “Quantum Inequalities in Quantum Mechanics,” Ann. Henri Poincaré 6, 1 (2005).
  * Penz  _et al._ (2006) M. Penz, G. Grübl, S. Kreidl, and P. Wagner, “A new approach to quantum backflow,” J. Phys. A: Math. Gen. 39, 423 (2006).
  * Trillo  _et al._ (2023) D. Trillo, T. P. Le, and M. Navascués, “Quantum advantages for transportation tasks - projectiles, rockets and quantum backflow,” npj Quantum Inf. 9, 69 (2023).