# Comparison of MPS based real time evolution algorithms for Anderson Impurity
Models

Daniel Bauernfeind<EMAIL_ADDRESS>Institute of Theoretical and
Computational Physics
Graz University of Technology, 8010 Graz, Austria  Markus Aichhorn  Institute
of Theoretical and Computational Physics
Graz University of Technology, 8010 Graz, Austria  Hans Gerd Evertz
<EMAIL_ADDRESS>Institute of Theoretical and Computational Physics
Graz University of Technology, 8010 Graz, Austria

###### Abstract

We perform a detailed comparison of two Matrix Product States (MPS) based time
evolution algorithms for Anderson Impurity Models. To describe the bath, we
use both the star-geometry as well as the commonly employed Wilson chain
geometry. For each bath geometry, we use either the Time Dependent Variational
Principle (TDVP) or the Time Evolving Block Decimation (TEBD) to perform the
time evolution. To apply TEBD for the star-geometry, we use a specially
adapted algorithm that can deal with the long-range coupling terms. Analyzing
the major sources of errors, one expects them to be proportional to the system
size for all algorithms. Surprisingly, we find errors independent of system
size except for TEBD in chain geometry. Additionally, we show that the right
combination of bath representation and time evolution algorithm is important.
While TDVP in chain geometry is a very precise approach, TEBD in star geometry
is much faster, such that for a given accuracy it is superior to TDVP in chain
geometry. This makes the adapted version of TEBD in star geometry the most
efficient method to solve impurity problems.

##  I Introduction

Describing strongly correlated electron materials is among the most difficult
tasks in solid state physics. One breakthrough in this field was the
development of the Dynamical Mean Field Theory (DMFT) Metzner and Vollhardt
(1989); Georges and Kotliar (1992); Georges  _et al._ (1996). DMFT accounts
for local electronic correlations by a self-consistent mapping of a lattice
problem, describing the low-energy subspace of the material, onto an Anderson
Impurity Model (AIM) Anderson (1961). The subsequent solution of this impurity
problem is the most important part of a DMFT calculation. At present,
Continuous Time Quantum Monte Carlo (CTQMC) Werner  _et al._ (2006); Gull  _et
al._ (2011) is the work-horse method for DMFT real material calculations. Many
other approaches exist, like the Numerical Renormalization Group Wilson
(1975); Bulla  _et al._ (2008), Configuration Interaction Lu  _et al._ (2014);
Zgid  _et al._ (2012); Mejuto Zaera  _et al._ (2017), and also methods based
on the Density Matrix Renormalization Group (DMRG) and the related Matrix
Product States (MPS) White (1992); Schollwöck (2011).
One of the major advantages of MPS based impurity solvers is that they can
give access to real-frequency spectral functions by employing real-time
evolution. While these methods provide excellent results for the single
orbital case, adding more and more orbitals is a very challenging task. To
overcome this issue, the Fork Tensor Product States (FTPS) Bauernfeind  _et
al._ (2017) solver was recently developed and applied to several materials
Bauernfeind  _et al._ (2017, 2018). This new approach has allowed to resolve
multiplets in the Hubbard bands of the single particle excitation spectrum for
real materials, which have also been observed experimentally Bauernfeind  _et
al._ (2018) but are inaccessible from the imaginary time results from standard
CTQMC methods. It is therefore important to identify the most efficient
methods for real-time evolution of impurity models.

A recent review Paeckel  _et al._ (2019) reported an extensive comparison of
several MPS-based time evolution algorithms for different types of models. It
did not, however, include impurity problems. They are special in the sense
that most of the degrees of freedom are non-interacting. The conclusion of
Ref. Paeckel  _et al._ , 2019 was that methods perform better or worse
depending on the model or observable studied. In all cases, the Time Dependent
Variational Principle (TDVP) Haegeman  _et al._ (2011, 2016) was among the
most reliable approaches, while the global Krylov method although accurate
tended to be too time-consuming.

Using tensor networks, impurity problems usually have been transformed to the
Wilson chain-geometry representation of the bath - essentially a nearest
neighbor tight-binding chain. The so-called star geometry on the other hand
involves hopping processes from the impurity to all bath sites that become
long ranged when mapped onto a linear chain. Although the reformulation of
DMRG as a variational principle on the space of MPS allows to deal with such
long-range terms, the chain geometry was believed to be the superior
representation of impurity models. Surprisingly, Wolf _et al._ Wolf  _et al._
(2014) demonstrated that for MPS, the star geometry is in fact a more economic
representation. Wolf _et al._ used the global Krylov method to perform the
time evolution which can be, as mentioned above, very time-consuming. As an
alternative, some of us published a modified Time Evolving Block Decimation
(TEBD) Vidal (2003, 2004); Daley  _et al._ (2004) algorithm for the star
geometry in Ref. Bauernfeind  _et al._ , 2017. There is also the possibility
to use the TDVP, which has been argued to be suited to long-range terms
Haegeman  _et al._ (2011). Therefore, open questions remain regarding the time
evolution methods used for impurity problems, for example: Does the bath
geometry affect the accuracy of the result? If so, does this depend on the
time evolution algorithm? And arguably the most important one: Which method is
best?

In the present paper, we use TDVP and TEBD as time evolution algorithms for
AIMs in the star and chain geometry and perform an in-depth comparison of the
four possibilities. We compare the quantity of interest in DMFT calculations,
i.e., the impurity Green’s function. We show that to obtain precise results,
the correct combination of bath representation and time evolution algorithm is
crucial. Specifically we demonstrate that the adapted TEBD Vidal (2003, 2004);
Daley  _et al._ (2004) is more accurate in the star geometry, while TDVP gives
better results in the chain geometry. We find that the TEBD approach in star-
geometry is much faster than TDVP in the chain geometry when a specific
precision of the Green’s function is prescribed. Additionally, we focus on the
behavior of the algorithms as a function of bath size. Surprisingly, all
algorithms except TEBD in chain-geometry have errors nearly independent of
system size. This very favorable scaling is especially important for DMFT
calculations, as one wants to use a large number of sites to represent the
bath hybridization well.

This paper is structured as follows. In Sec. II we introduce AIMs in star-
geometry and in chain-geometry representation of the bath and discuss how to
obtain bath parameters from the bath hybridization. In Sec. III, we briefly
introduce MPS and how to obtain Green’s functions using real-time evolution.
In Sec. IV we discuss the different time evolution methods, including a
discussion of the various error sources. Sec. V contains results, first for
the non-interacting case, and in Sec. V.2 for the interacting case. Finally,
Sec. VI contains the conclusions.

##  II Anderson Impurity Models

In the star geometry, an impurity described by a local interacting Hamilton
${H_{\text{loc}}}$ is coupled to a bath of free fermions via hopping terms
from the impurity to _every_ bath site. For a one-band model this results in:

$\displaystyle H_{\text{star}}$ $\displaystyle=H_{\text{loc}}+H_{\text{bath}}$
(1) $\displaystyle H_{\text{loc}}$
$\displaystyle=Un_{I,\uparrow}n_{I,\downarrow}+\epsilon_{I}\left(n_{I,\uparrow}+n_{I,\downarrow}\right)$
$\displaystyle H_{\text{bath}}$
$\displaystyle=\sum_{k,\sigma}H_{k,\sigma}=\sum_{k,\sigma}\epsilon_{k}n_{k,\sigma}+V_{k}\left(c^{\dagger}_{I,\sigma}c_{k,\sigma}+c^{\dagger}_{k,\sigma}c_{I,\sigma}\right).$

Here, $c_{k,\sigma}^{\dagger}$ $(c_{k,\sigma})$ create (annihilate) an
electron at bath site $k$ with spin $\sigma$,
$n_{k,\sigma}=c^{\dagger}_{k,\sigma}c_{k,\sigma}$ are the usual particle
number operators and we label the impurity degrees of freedom by an index $I$.
If one puts the sites of $H_{\text{star}}$ on a one-dimensional manifold (like
an MPS - see below), the hopping terms become long-range.

Using a Lanzcos-like tridiagonalization, one can map the star geometry onto
the so-called Wilson chain Wilson (1975); Bulla  _et al._ (2008) or just
chain-geometry, where the impurity couples to the first bath site only:

$\displaystyle H_{\text{chain}}$
$\displaystyle={H_{\text{loc}}}+t_{0}\left(c^{\dagger}_{I,\sigma}c_{1,\sigma}+c^{\dagger}_{1,\sigma}c_{I,\sigma}\right)+$
(2)
$\displaystyle\sum_{i,\sigma}t_{i}\left(c^{\dagger}_{i,\sigma}c_{i+1,\sigma}+c^{\dagger}_{i+1,\sigma}c_{i,\sigma}\right)+\bar{\epsilon}_{i}n_{i\sigma}.$

###  II.1 Determination of Bath parameters

The bath can also be described by the continuous (DMFT Georges  _et al._
(1996)) hybridization function $\Delta(\omega)$. By tracing out the bath in
Eq. (1), one finds:

$\displaystyle\Delta(\omega)$
$\displaystyle=\sum_{k}\frac{V_{k}^{2}}{\omega-\epsilon_{k}+i0^{+}}$ (3a)
$\displaystyle-\frac{1}{\pi}\mathcal{I}\Delta(\omega)$
$\displaystyle\approx\sum_{k=1}^{N}V_{k}^{2}\frac{\eta}{\pi\left(\left(\omega-\epsilon_{k}\right)^{2}+\eta^{2}\right)}.$
(3b)

A continuous bath corresponds to an infinite sum. In the last line, we used
only a finite number of $N$ bath sites, and approximated the delta-peaks by
Lorentzians of finite width $\eta$. Usually, the hybridization function
$\Delta(\omega)$ is given and in order to map it to an AIM of finite size, one
needs to find values $\epsilon_{k}$ and $V_{k}$, such that Eq. (3b) is as good
approximation. In the present study, we choose to split the $\omega$-axis into
$N$ equidistant intervals $I_{k}$ of size $\Delta\epsilon$. We describe each
interval using a single bath site with on-site energy $\epsilon_{k}=\min
I_{k}+\frac{\Delta\epsilon}{2}$, where $\min I_{k}$ is the minimum of interval
$I_{k}$. The hopping amplitude $V_{k}$ can then be computed from Bulla  _et
al._ (2008); Wolf  _et al._ (2014):

$V_{k}^{2}=\int_{I_{k}}-\frac{1}{\pi}\mathcal{I}\Delta(\omega).$ (4)

A subsequent basis transformation into the Lanzcos-basis of
${H_{\text{bath}}}$ yields the bath parameters in the chain geometry Bulla
_et al._ (2008); Wolf  _et al._ (2014). Note that for particle hole symmetry,
the on-site energies $\bar{\epsilon}_{i}$ of the chain geometry are exactly
zero Bulla  _et al._ (2008). For both geometries, we enforce particle hole
symmetry in the bath parameters. In the following, we will give all results in
units of the half bandwidth of the bath spectral function
$-\frac{1}{\pi}\mathcal{I}\Delta(\omega)$.

##  III Matrix Product States

Figure 1: Graphical representation of an MPS used for an AIM with $N=3$ bath
sites. _Top:_ To separate the spin species, we place the impurity in the
middle (green circles) and attach the spin-up (spin-down) bath to the left
(right) colored in blue (red). _Bottom:_ Labeling of sites in the star- as
well as in the chain-geometry. The bath sites are arranged in ascending order
with increasing distance from the impurity.

MPS are an efficient parametrization of quantum mechanical states as a product
of local matrices. Consider a system consisting of $2N+2$ sites with a local
basis $\ket{s_{i}}$ at site $i$:

$\ket{\psi}=\sum_{\\{s\\}}c_{s_{1},\cdots s_{2N+2}}\ket{s_{1}\cdots
s_{2N+2}}.$ (5)

In an MPS, the coefficient $c_{s_{1},\cdots s_{2N+2}}$ is factorized into a
product of matrices using repeated Singular Value Decompositions (SVDs)
Schollwöck (2011):

$c_{s_{1},\cdots s_{2N+2}}=A^{s_{1}}\cdot A^{s_{2}}\cdots A^{s_{2N+1}}.$ (6)

Each $A_{i}^{s_{i}}$ is a rank-$3$ tensor, except the two tensors at the edges
which are of rank-$2$. Since $s_{i}$ represents the local basis $\ket{s_{i}}$,
this index is called physical. The _matrix_ indices which are summed over are
called bond-indices. The dimension of the bond indices (bond dimension) is the
number of Schmidt values kept during the calculation, implying that some
Schmidt values are discarded. The sum of the square of all discarded Schmidt
values is called truncated weight $t_{w}$ Schollwöck (2011). MPS and other
tensor networks are often depicted using a graphical representation as in Fig.
1. Since the Hamiltonian of an AIM connects the two spin species only via
interactions on the impurity (see Eq. (1)), it has turned out to be favorable
to separate them in the MPS, and to use a local Hilbert space of dimension two
(empty and occupied) for each site Saberi  _et al._ (2008); Ganahl  _et al._
(2015); Bauernfeind  _et al._ (2017). Thus, we place the impurity in the
middle of the chain and connect the spin-up (spin-down) bath to its left
(right), as shown in Fig. 1.
Using this arrangement of sites, we calculate the ground state
$\ket{\psi_{0}}$ and its energy $E_{0}$ using DMRG. The central object of
interest in DMFT calculations is the Green’s function of the impurity. In the
following, we will hence focus our attention on the greater Green’s function
of the impurity (omitting the spin index):

$\displaystyle G^{>}(t)$
$\displaystyle=\bra{\psi_{0}}c_{I}e^{-iHt}c_{I}^{\dagger}\ket{\psi_{0}}e^{iE_{0}t}$
(7)
$\displaystyle=\left(e^{iH\frac{t}{2}}c_{I}^{\dagger}\ket{\psi_{0}}\right)^{\dagger}\left(e^{-iH\frac{t}{2}}c_{I}^{\dagger}\ket{\psi_{0}}\right)e^{iE_{0}t}.$

To calculate Eq. (7), we first apply the operator $c^{\dagger}_{I}$ onto the
ground state. Next, as indicated by brackets in the the second line of Eq.
(7), we perform two separate time evolutions up to time $\frac{t}{2}$ Ganahl
_et al._ (2015); Kennes and Karrasch (2016) and calculate the overlap. We
checked that the results are the same for the lesser Green’s function and also
checked that non-particle hole symmetric models show a behavior very similar
to the results presented below.

##  IV Time Evolution Algorithms

In the following, we discuss the different time evolution algorithms with
special attention to our adapted TEBD approach in star geometry, and we
analyze the main sources of error.

###  IV.1 TEBD

Several slightly different formulations of TEBD Vidal (2003, 2004); Daley  _et
al._ (2004); Paeckel  _et al._ (2019) exist. The common strategy is to split
the full time evolution operator of some small time step $\tau$ into
manageable parts which can be applied to evolve the state forward in time.
Specifically, one employs Suzuki-Trotter breakups Suzuki (1990) to obtain a
decomposition into operators acting on two sites only. The action of such a
_gate_ merges the MPS tensors of the two sites, which are then separated again
using a SVD combined with a truncation. Apart from the truncation, the main
error of this approach is caused by the (in our case second order) Trotter
breakup:

$\displaystyle e^{\tau(A+B)}$
$\displaystyle=e^{-iA\frac{\tau}{2}}e^{-iB\tau}e^{-iA\frac{\tau}{2}}+\frac{\tau^{3}}{12}\mathcal{C}+\mathcal{O}(\tau^{4})$
(8) $\displaystyle\mathcal{C}$
$\displaystyle=\left(\frac{1}{2}\big{[}A,[A,B]\big{]}+\big{[}B,[A,B]\big{]}\right).$

According to Eq. (8), the total error of the time evolution is not only
determined by the time step $\tau$, but also by the matrix elements of the
double commutators $\mathcal{C}$. Since the latter are different for the star
and the chain geometry, the error of TEBD in both geometries can be very
different. This will be discussed in the following two subsections.

#### Chain Geometry

To employ TEBD in the chain-geometry, we use the standard second order breakup
between even and odd terms Schollwöck (2011). We write the Hamiltonian of the
chain geometry $H_{\text{chain}}=\sum_{i}h_{i,i+1}$ as a sum of local terms
$h_{i,i+1}$ acting on nearest neighbor sites $i$ and $i+1$ only and define:

$\displaystyle H_{\text{chain}}$
$\displaystyle=H_{\text{even}}+H_{\text{odd}}$ $\displaystyle H_{\text{even}}$
$\displaystyle=\sum_{i:\text{even}}h_{i,i+1}\text{, }$ $\displaystyle
H_{\text{odd}}$ $\displaystyle=\sum_{i:\text{odd}}h_{i,i+1}$ $\displaystyle
e^{-iH_{\text{chain}}{\Delta t}}$ $\displaystyle\approx
e^{-iH_{\text{odd}}\frac{{\Delta t}}{2}}e^{-iH_{\text{even}}{\Delta
t}}e^{-iH_{\text{odd}}\frac{{\Delta t}}{2}}.$ (9)

In the particle hole symmetric case, the on-site energies $\bar{\epsilon_{i}}$
are exactly zero, removing any ambiguity of how to distribute the on-site
terms among even- and odd parts of the Hamiltonian. For particle hole symmetry
($\bar{\epsilon_{i}}=0$), we evaluated the double commutators of Eq. (8). They
turn out to be various sums over hopping terms multiplied by three amplitudes:

$\mathcal{C}_{\text{chain}}=-\frac{1}{2}\sum_{i:\text{even}}t_{i}t_{i+1}t_{i+2}\left(c^{\dagger}_{i}c_{i+3}+c^{\dagger}_{i+3}c_{i}\right)+\cdots.$
(10)

Additionally, the interaction $U$ gives two terms that couple to the neighbors
of the impurity, independent of $N$. For the scaling with bath size the latter
terms can be neglected, since they are of order $\mathcal{O}(1)$. Each of the
individual sums of Eq. 10 scales linearly with system size $N$, they contain
$N$ terms of order $1$, and the hopping amplitudes $t_{i}$ are independent of
$N$ 111For example a semi-circular bath spectral function can be represented
by $t_{i}=t$ for all sites. Note that we use the Lanzcos tri-diagonalization
instead.. Thus, one should expect the overall error of TEBD in chain geometry
to scale like the bath size $N$.

#### Star Geometry

Figure 2: Application of the time evolution operator of Eq. (12) onto the
spin-down bath. We use swap gates depicted by two crossing arcs to evolve the
long-range hopping terms. The green line visualizes the position of the
impurity as it is moved through the bath during the time evolution. All gates
except the one with $H_{N,\sigma}$ are calculated with a time step
$\frac{{\Delta t}}{2}$, reflecting the use of the second order breakup. Note
that with this approach one can combine every swap gate with an actual time
evolution gate, thus avoiding any additional SVDs.

To perform the TEBD time evolution of an AIM in star geometry (Eq (1)), we
employ an approach we recently developed (see Refs Bauernfeind  _et al._ ,
2017; Bauernfeind, 2018), which is based on swap-gates Orús (2014); Verstraete
_et al._ (2009) and iterative Suzuki-Trotter decompositions. First, we first
split off ${H_{\text{loc}}}$ from the time evolution operator for some small
time-step ${\Delta t}$ in a second order Trotter decomposition:

$e^{-iH{\Delta t}}\approx e^{-i{H_{\text{loc}}}\frac{{\Delta
t}}{2}}\left(\prod_{\sigma}e^{-i\sum_{k}H_{k,\sigma}{\Delta
t}}\right)e^{-i{H_{\text{loc}}}\frac{{\Delta t}}{2}},$ (11)

where $H_{k,\sigma}$ was defined in Eq. (1). Then we split off the term
containing $H_{1,\sigma}$, then $H_{2,\sigma}$ and so on, until we obtain:

$e^{-i\sum_{k}H_{k,\sigma}{\Delta
t}}\approx\left(\prod_{k=1}^{N}e^{-iH_{k,\sigma}\frac{{\Delta
t}}{2}}\right)\left(\prod_{k=N}^{1}e^{-iH_{k,\sigma}\frac{{\Delta
t}}{2}}\right).$ (12)

The Trotter errors of these approximations will be discussed below. Note that
Eq. (12) is a product of operators acting on two sites only, but all terms
except $H_{1,\sigma}$ involve sites that are not nearest neighbors in the MPS.
To be able to apply such non-nearest neighbor gates, we use swap-gates. Their
purpose is to swap the degrees of freedom of two sites in the MPS-tensor
network. For a fermionic Hilbert space with local dimension of $2$ (empty
$\ket{0}$ or occupied $\ket{1}$), a swap-gate $S$ is given by:

$S=\ket{00}\bra{00}+\ket{01}\bra{10}+\ket{10}\bra{01}-\ket{11}\bra{11}.$ (13)

Note the minus sign in the last term, resulting from the exchange of two
fermions Bauernfeind (2018). The process of applying all the gates in Eq. (12)
is depicted in Fig. 2 and described in the following.

Instead of a gate with $H_{1,\sigma}$ only, we apply a combined gate $S\cdot
e^{-iH_{1,\sigma}\frac{{\Delta t}}{2}}$, meaning that we first time evolve and
swap afterwards. The impurity degrees of freedom are then located at what was
previously the first bath site and is now a nearest neighbor of the second
bath site. We continue by applying a combined two-site gate $S\cdot
e^{-iH_{2,\sigma}\frac{{\Delta t}}{2}}$ after which the impurity and bath site
$3$ are nearest neighbors and so on. When the impurity arrives at site $N-1$,
we time evolve with $H_{N,\sigma}$ (without swapping). Since we used a second
order decomposition, we have to re-apply all gates with $k=N-1\cdots 1$. This
time, though, we have to swap first and time evolve afterwards, since
otherwise we would have to take care of an additional fermionic sign in the
hopping terms Bauernfeind (2018). This means that we apply a two-site gate
$e^{-iH_{N-1,\sigma}\frac{{\Delta t}}{2}}\cdot S$ followed by
$e^{-iH_{N-2,\sigma}\frac{{\Delta t}}{2}}\cdot S$, etc., until the impurity is
back in the middle of the MPS and every term in Eq (12) is dealt with.

Now let us evaluate the Trotter errors due to these various decompositions. In
the first breakup, we split off ${H_{\text{loc}}}$ from the rest of the
Hamiltonian, i.e., $A={H_{\text{loc}}}$ and $B={H_{\text{bath}}}$ in Eq. (8).
Evaluation of the double commutators yields:

$\displaystyle\mathcal{C}_{{H_{\text{loc}}}}=-2\sum_{\sigma}K_{\sigma}n_{I,\sigma}\sum_{k}V_{k}^{2}$
(14)
$\displaystyle+\sum_{k,\sigma}V_{k}\left(c^{\dagger}_{I,\sigma}c_{k,\sigma}+c^{\dagger}_{k,\sigma}c_{I,\sigma}\right)K_{\sigma}\left(\frac{1}{2}K_{\sigma}-\epsilon_{k}\right)$
$\displaystyle+\sum_{k,k^{\prime},\sigma}V_{k}V_{k^{\prime}}K_{\sigma}\left(c^{\dagger}_{k,\sigma}c_{k^{\prime},\sigma}+c^{\dagger}_{k^{\prime},\sigma}c_{k,\sigma}\right)$
$\displaystyle-\sum_{k,k^{\prime},\sigma}V_{k}V_{k^{\prime}}U\left(c^{\dagger}_{I,\sigma}c_{k^{\prime},\sigma}-c^{\dagger}_{k^{\prime},\sigma}c_{I,\sigma}\right)\left(c^{\dagger}_{I,\bar{\sigma}}c_{k,\bar{\sigma}}-c^{\dagger}_{k,\bar{\sigma}}c_{I,\bar{\sigma}}\right),$

with $K_{\sigma}=Un_{I\bar{\sigma}}+\epsilon_{0}$ and $\bar{\sigma}$ the spin
in opposite direction of spin $\sigma$.
Next, we split off $H_{1,\sigma}$ from $\sum_{k>1}H_{k,\sigma}$, i.e.,
$A=H_{1,\sigma}$ and $B=\sum_{k>1}H_{k,\sigma}$. The error of this
decomposition (with $j=1$):

$\displaystyle\mathcal{C}_{H_{j,\sigma}}$
$\displaystyle=\frac{1}{2}\sum_{k>j}V_{j}^{2}V_{k}\left(c^{\dagger}_{I,\sigma}c_{k,\sigma}+c^{\dagger}_{k,\sigma}c_{I,\sigma}\right)$
(15) $\displaystyle-
V_{j}\left(c^{\dagger}_{I,\sigma}c_{j,\sigma}+c^{\dagger}_{j,\sigma}c_{I,\sigma}\right)\sum_{k>j}V_{k}^{2}$
$\displaystyle+\sum_{k>j}V_{j}V_{k}\left(c^{\dagger}_{j,\sigma}c_{k,\sigma}+c^{\dagger}_{k,\sigma}c_{j,\sigma}\right)\left(\frac{1}{2}\epsilon_{j}-\epsilon_{k}\right).$

The next decomposition separating $H_{2,\sigma}$ from $\sum_{k>2}H_{k,\sigma}$
has the exact same form as above now with $j=2$, and hence an error of
$\mathcal{C}_{H_{j,\sigma}}$ for $j=2$. Iterating over all decompositions, the
total error of the breakup used to separate the time evolution operator in the
star geometry becomes:

$\mathcal{C}_{\text{star}}=\mathcal{C}_{H_{\text{loc}}}+\sum_{j\neq
N,\sigma}\mathcal{C}_{H_{j,\sigma}}.$ (16)

Let us analyze how this error scales with the number of bath sites $N$. From
Eq. (3a), we find that $V_{k}\sim\frac{1}{\sqrt{N}}$. Hence,
$\mathcal{C}_{H_{\text{loc}}}$ in Eq. (14) scales linearly with $N$, since
each of the last two lines contains $N^{2}$ terms of order $1$ multiplied by
$V_{k}^{2}\sim\frac{1}{N}$.
Similarly, the other terms in Eq.16 give a growth of the error not faster than
$N$ (two terms in $\mathcal{C}_{H_{\text{loc}}}$ also have $\sim N$ scaling
and the next lowest order is a $\sqrt{N}$-scaling). We hence should expect the
error in the star geometry to scale linearly with $N$.
This might seem surprising at first, since to obtain Eq. (12), we performed
$N$ individual decompositions and one might therefore expect the total error
to scale at least $\sim N^{2}$. However, since the hopping parameters are
$N$-dependent themselves ($V_{k}\sim\frac{1}{\sqrt{N}}$), the overall error
improves by a factor of $\frac{1}{N}$. Note that the actual errors in the
calculations will depend on the matrix elements of the error terms Eq.(10) or
Eq. (16).

###  IV.2 TDVP

Instead of approximating the time evolution operator, TDVP directly solves the
time dependent Schrödinger equation, albeit only in a restricted subspace - in
the space of MPS of fixed bond dimension. TDVP constructs a projection
operator $P_{T_{\ket{\psi}}}$ that projects the right hand side of the
Schrödinger equation onto the tangent space $T_{\ket{\psi}}$ of the current
MPS $\ket{\psi}$:

$\frac{\partial}{\partial t}\ket{\psi}=-iP_{T_{\ket{\psi}}}H\ket{\psi}.$ (17)

This equation then results in a set of equations that can be integrated with
an approach very similar to DMRG, replacing the ground-state search by Krylov
time propagation Haegeman  _et al._ (2016); Paeckel  _et al._ (2019). In the
present publication we employ two-site TDVP (2TDVP in Ref. Paeckel  _et al._ ,
2019) using the second order integrator by sweeping left-right-left with half
time step $\frac{\Delta t}{2}$. For additional details on TDVP we refer to the
existing literature Haegeman  _et al._ (2011, 2016); Paeckel  _et al._ (2019).
Apart from MPS-matrix truncation and finite time step $\Delta t$, TDVP has an
additional parameter, namely when to terminate the Krylov series. We stop
creating new Krylov vectors when the total contribution of two consecutive
vectors to the matrix exponential is less than $10^{-15}$.

The errors of TDVP are, apart from the truncation, a time step error similar
to TEBD of order $\Delta t^{3}$ and an error due to the projection of the
Schrödinger equation Paeckel  _et al._ (2019). We expect the latter to
strongly depend on the bath representation, because it is exactly zero for
Hamiltonians with nearest neighbor terms only Paeckel  _et al._ (2019);
Haegeman  _et al._ (2016). As a function of bath size, we expect an error
linear in $N$, since TDVP approximates $N$ coupled equations by integrating
them one after the other.

##  V Results

Figure 3: Absolute differences $|\Re G^{>}_{\text{exact}}(t)-\Re
G^{>}_{\text{DMRG}}(t)|$ for chain- (C) and star (S) geometries and various
bath sizes ($N=29,59,109$ from top to bottom) at $U=0$. The bath for this
calculation was obtained from a semi-circular spectral function
$-\frac{1}{\pi}\mathcal{I}\Delta(\omega)=\frac{1}{2\pi}\sqrt{1-\omega^{2}}$.
The Trotter time-step was ${\Delta t}=0.05$. The truncated weight was set to
$10^{-12}$ during DMRG as well as during time evolution and the bond
dimensions were not restricted to any maximal value.  Figure 4: Maximum error
for times $20<t<30$ as a function of bath size. Top: $U=0$ calculation shown
in Fig. 3. Bottom: $U=1$ calculation shown in Fig. 6. Note that the errors of
TEBD-S and TDVP-C have been multiplied by a factor of 50 to make the graphs
comparable.

To compare the different algorithms which have different sources of errors, we
use the following strategy. We first fix the parameters (time step $\Delta t$
and truncated weight $t_{w}$) and compare the error of the Green’s function
without taking the computation time into account. In the next step, we then
address the real question of computation time versus error.

###  V.1 U=0

Since interactions only affect the impurity degrees of freedom, it is
reasonable to expect that most of the errors from the approximations in the
time evolution is already present in the non-interacting case. Starting with
$U=0$ has the advantage of giving us access to the exact Green’s function by
diagonalization of the hopping matrix $T_{ij}$:

$\displaystyle
H=\sum_{ij}c^{\dagger}_{i}T_{ij}c_{j}=\sum_{ijk}\underbrace{c^{\dagger}_{i}U^{\dagger}_{ik}}_{c^{\dagger}_{k}}E_{k}\underbrace{U_{kj}c_{j}}_{c_{k}}=\sum_{k}E_{k}n_{k},$
(18)

where $H=H_{\text{star}}$ or $H=H_{\text{chain}}$, respectively. Note that in
the following, we compare each calculation to the exact Green’s function for
the system defined by its hopping matrix $T_{ij}$ for the finite size bath.
In Fig. 3 we compare the four time evolution schemes for various system sizes
and for fixed parameters $\Delta t$ and truncation. The truncated weight was
chosen very small ($10^{-12}$), such that the main error source is the
approximation of the time evolution operator, not the truncation of the tensor
network. The magnitude of the error suggests that TDVP in chain geometry
(TDVP-C) is the best algorithm followed by the TEBD in star geometry (TEBD-S).
This does not consider the necessary bond dimension $m$ though, which was very
different depending mostly on the bath geometry. While in star geometry the
bond dimensions were $48$ and $62$ for TDVP and TEBD respectively, they grew
to $145$ and $220$ in chain geometry. This difference has drastic effects on
the computation time discussed below, which scales $\sim m^{3}$.
Three observations are especially interesting in Fig. 3. First, for TDVP the
chain geometry gives more precise results, whereas for TEBD the star geometry
gives a lower error. For TDVP this behavior may be related to the projection
error being zero in the chain geometry Haegeman  _et al._ (2016).

Second, in all cases, the error seems to be about constant in time. This is
true until the particle added in the calculation of the Green’s function is
reflected at the boundaries of the finite size system (chain geometry) and
then reaches the impurity again. Somewhat surprisingly, the same time scale
also holds for the star geometry, even though the picture of a particle
traveling towards the end of the bath applies only to the chain geometry. On
the other hand the star geometry and the chain geometry are equivalent by a
unitary transformation that can explain this time scale. The increase in error
in the top plot of Fig. 3 (N=29) beginning at times $t>60$ is due to such a
reflection.
Third, as a function of bath size, only the error of TEBD-C shows the expected
linear scaling with $N$ as also shown in Fig. 4 (top). The error of all other
algorithms are almost exactly constant in the bath size. For TEBD-S and for
TDVP, the error even seems to become slightly smaller as $N$ is increased.
This is very surprising considering Eq. (16) and the $N$ successive
approximations of TDVP.

Figure 5: Matrix elements $|U_{ik}|$ (see Eq. (18)) in the star geometry
(left) and in the chain geometry (right) for $N=159$ for the semi-circular
bath as in Fig. 3. While the non-zero entries in star geometry are strongly
concentrated around the diagonal, $U_{ik}$ in the chain geometry has $\sim
N^{2}$ relevant non-zero elements. Note that for the star geometry, the terms
$U_{ik}=V_{k}$ for $i=I$ (see Eq. (22)) are located at the bottom of the plot,
but are barely visible.

We can understand the observed behavior for TEBD in both geometries by
examining the matrix elements of the Trotter error terms in Eq. (16) and Eq.
(10). We start with the many-body ground state for $U=0$, i.e., the filled
Fermi-sea (FS):

$\ket{\psi_{0}}=\prod_{k\in FS}c^{\dagger}_{k}\ket{0},$ (19)

where $k$ labels the eigenstates of the hopping matrix $T_{ij}$ (Eq. (18)). As
we have seen above, the leading errors of the Trotter decompositions
correspond to hopping terms from one site to another. The error for the
Green’s function is then given by the matrix element of Eq. (15) (star) and
Eq. (10) (chain) with the state $c_{I}^{\dagger}\ket{\psi_{0}}$. The exact
time evolution operator for $U=0$ is simply $U(t)=\prod_{k}e^{-iE_{k}t}n_{k}$,
i.e., a time dependent phase factor for each occupied $k$. To find the
magnitude and the number of terms contributing to the total error it hence
suffices to evaluate (using the exact ground state (the Fermi-sea) for
$\bra{\psi_{0}}$):

$\bra{\psi_{0}}c_{I}c^{\dagger}_{i}c_{j}c_{I}^{\dagger}\ket{\psi_{0}}.$ (20)

The indices $i$ and $j$ are determined by the various hopping terms appearing
in Eq. (15) and Eq. (10). In $k$-space, we find:

$\displaystyle\bra{\psi_{0}}c_{I}c^{\dagger}_{i}c_{j}c_{I}^{\dagger}\ket{\psi_{0}}=$
$\displaystyle\sum_{k>FS}U^{\dagger}_{kI}U_{Ik}\sum_{k^{\prime}\in
FS}U^{\dagger}_{k^{\prime}i}U_{jk^{\prime}}$ (21) $\displaystyle+$
$\displaystyle\sum_{k,k^{\prime}>FS}U^{\dagger}_{Ik^{\prime}}U_{kI}U^{\dagger}_{ki}U_{jk^{\prime}},$

where $U_{ik}$ are the matrix elements of the unitary transformation in Eq.
(18), and the index value $I$ again denotes the impurity. Note that this
expression is valid for both the star- and the chain geometry. Differences
between them are encoded in the different unitary transformations $U_{ik}$
diagonalizing the hopping matrix $T_{ij}$. In star geometry, the bath states
with energy $\epsilon_{k}$ are already very close to the eigenstates of
$T_{ij}$. This implies that most entries in $U_{ik}$ are nearly zero, except
for a few values around $i=k$. Indeed, as Fig. 5 (left plot) demonstrates,
$U_{ik}$ in star geometry has relevant non-zero entries only around the
diagonal. This means that in star geometry there are only $\sim N$ relevant
terms in the matrix products in Eq. 21. In other words, not all terms of
$\mathcal{C}_{\text{star}}$ contribute and we expect an error independent of
system size for the leading order $({\Delta t})^{3}$. To make this point more
clear, let us look at the case where only the diagonal of $U_{ik}$
contributes:

$U_{ik}^{\text{star}}\approx\left\\{\begin{array}[]{lr}\delta_{i,k}\text{, for
}i\neq I\\\ \sim V_{k}\text{, for }i=I\end{array}\right.$ (22)

We note that the approximation $U_{Ik}\sim V_{k}$ is in good agreement with
the true form of $U_{ik}$. Because of Eq. (22), $\mathcal{C}_{\text{star}}$
does not scale with system size anymore, since $\delta_{ik}$ always removes at
least one summation. This remains true, when $U_{ik}^{\text{star}}$ contains a
finite width band of relevant values instead of just the diagonal. For the
chain geometry on the other hand $U_{ik}$ is a full matrix with $\sim N^{2}$
entries of similar size, as shown in Fig. 5 (right plot). Therefore, no such
simplification occurs and the error should indeed scale linearly with $N$.
For TDVP, such arguments do not hold, since its error is independent of the
bath geometry used. The evaluation of the corresponding commutators is far
from trivial and we hence leave this point open for future studies.

###  V.2 Finite $U$

Figure 6: Absolute difference of $\Re G^{>}_{\text{ref}}(t)$ between high
precision reference calculations and calculations with lower precision for
various bath sizes and $U=1$. The bath was obtained from a semi-circular
spectral function
$-\frac{1}{\pi}\mathcal{I}\Delta(\omega)=\frac{1}{2\pi}\sqrt{1-\omega^{2}}$.
Because of the results at $U=0$, we chose TEBD for the star geometry and TDVP
for the chain geometry for the reference calculations. These used a truncated
weight of $10^{-14}$ ($10^{-13}$) and a Trotter time-step ${\Delta t}=0.005$
(${\Delta t}=0.005$) for TEBD (TDVP) respectively. Parameters for the other
calculations were ${\Delta t}=0.05$ and truncated weight $10^{-12}$, the same
as in Fig.3. The bond dimensions of the MPS were not restricted to a maximal
value.

In Fig. 6 we show a comparison similar to Fig. 3, now for $U=1$. Since in the
interacting case we do not have access to the exact solution, we compare to
reference calculations with very high precision done separately for each bath
size (see Fig. 6 for details). Overall, we find that the errors of lower
precision calculations in Fig. 6 are almost identical to the ones obtained in
the non-interacting case. In particular, we again find that only the error of
TEBD-C scales appreciably with system size (see Fig. 4).
The necessary bond dimension on the other hand changes drastically. Using the
semi-circular bath spectral function, the maximal value was
$m_{\text{star}}=364$ $(128)$ compared to $m_{\text{chain}}=260$ $(214)$ for
TEBD (TDVP) respectively.
So far, we only compared the error for a given set of parameters (i.e.,
truncation $t_{w}$ and time step $\Delta t$). The actual quantity of interest
is the computation time for given accuracy (or vice versa). Although all
algorithms scale$~{}\sim m^{3}$, their computation times are very different
and comparisons are not straightforward for several reasons.

First, in star geometry the bond dimensions is strongly peaked around the
center bath site, whereas in chain geometry it is more flat. Therefore, the
maximal bond dimension is not a good indicator of actual computation times.
Second, a single TDVP step is generally much more expensive than a TEBD step.
For example, in the calculations shown in Fig. 6, TEBD-S is faster than TDVP-S
although the maximal bond dimension is larger by a factor of three in TEBD-S.
Third, TDVP generally allows for much larger time steps $\Delta t$ for a give
accuracy.

Additionally, the advantage of the different geometries will likely depend on
the actual bath parameters. For DMFT, most of the calculations are performed
close to the self consistent point. We therefore studied the computation times
necessary to reach a prescribed precision at the self consistent point of the
Bethe lattice at $U=2$ ($U=4$ in the units of Ref. Wolf  _et al._ , 2014) for
a time evolution up to $t=15$.

Alg. | $\Delta t$ | $t_{w}$ | Error ($10^{-4}$) | Wall time (s)
---|---|---|---|---
TEBD-C | $0.005$ | $10^{-10}$ | $2.5$ | 2370
TDVP-S | $0.005$ | $10^{-12}$ | $1.0$ | 4434
TEBD-S | 0.1 | $10^{-9}$ | 1.7 | 35
TEBD-S | 0.01 | $10^{-11}$ | 0.14 | 384
TEBD-S | 0.01 | $10^{-13}$ | 0.01 | 1655
TDVP-C | 0.5 | $10^{-8}$ | 2.6 | 457
TDVP-C | 0.1 | $10^{-10}$ | 0.2 | 1856
TDVP-C | 0.05 | $10^{-12}$ | 0.02 | 7130
Table 1: Comparison of computation times for the $U=2$ Bethe lattice self
consistent bath. We chose example parameters for the time step $\Delta t$ and
the truncated weight $t_{w}$ to approximately obtain errors of different
orders of magnitude $10^{-4}$ to $10^{-6}$. We used the largest time step out
of $\\{0.005,0.01,0.05,0.1,0.5\\}$ for which this error can be achieved with
low enough truncated weights and then used the largest truncated weight at
this time step with a similar error. For TEBD-C and TDVP-S only an error of
about $10^{-4}$ was possible with the parameters studied. Reference
calculations were again performed with TDVP-C in chain-geometry ($\Delta
t=0.005$, $t_{w}=10^{-13}$) and TEBD-S in star geometry ($\Delta t=0.005$,
$t_{w}=10^{-14}$). Wall times are reported for a Intel(R) Core(TM) i7-7740X
CPU using a single thread and for a single calculation of $G^{>}(t)$ up to
$t=15$ including DMRG for $N=59$.

Results are shown in Tab. 1. The two algorithms with large errors in Figs. 3
and 6, TEBD-C and TDVP-S are not able to obtain errors smaller than
approximately $10^{-4}$ for the parameters studied and are comparable in
computation time (TDVP-S has about half the error of TEBD-C but double the
computation time). Additionally, they are slow compared to the other two
approaches confirming behavior seen in Fig. 3 and Fig. 6.
Let us now compare the two favorable algorithms, TDVP-C and TEBD-S. While from
Fig. 3 and Fig. 6 on would expect TDVP-C to be the better algorithm, Tab. 1
clearly shows that TEBD in star geometry is actually superior to TDVP-C. Its
computation times are lower by about a factor of $5$ or more for the same
error. Conversely, for similar computation times, the error of TEBD-S is about
one order of magnitude smaller than TDVP-C. This shows that TEBD in star
geometry is the best approach to calculate impurity Green’s functions and
should also be considered the time evolution algorithm of choice for general
non-equilibrium impurity problems as clever arrangement of sites allows time
evolutions to be performed up to surprisingly long times Rams and Zwolak
(2019).

##  VI Conclusions

We compared tensor network time evolution algorithms for Anderson Impurity
Models for different bath representations. Specifically, we used TEBD and TDVP
for the star- as well as the chain-geometry. TDVP is readily applicable for
the long-range hybridizations present in the star geometry. For TEBD this is
not the case and we explained in some detail the adapted TEBD approach using
swap-gates first published in Ref. Bauernfeind  _et al._ , 2017. Its major
advantage is that each actual time evolution gate can be combined with a swap
gate, involving no additional computational cost and thus preserving the
simplicity of TEBD. For TEBD, we additionally performed an analytic
calculation of the leading order of the error due to the Suzuki-Trotter
decomposition in both the star- and chain-geometry. This and the
approximations of TDVP led us to expect that the error in the Green’s function
should be proportional to the system size $N$ in all four cases. Surprisingly,
we found that only TEBD in chain-geometry shows this behavior. We were able to
find an analytical explanation of the better scaling of TEBD in star geometry
from the fact that the bath states in star geometry are already a good
approximation to the single particle eigenbasis and therefore, most error
terms do not contribute. For DMFT calculations, such a favorable scaling with
system size is especially important, since one has to make sure to use large
bath sizes $N$ to represent the bath hybridization well enough to reach the
correct self-consistent point.
Regarding the magnitude of the error, it is important to use the best
combination of time evolution algorithm and bath representation. TDVP is more
precise in the chain geometry likely due to the absence of projection error in
the chain geometry. For TEBD it turned out to be vice versa, i.e., star
geometry has a lower error. With given set of parameters (time step $\Delta t$
and truncated weight $t_{w}$) TDVP in chain geometry has the lowest error. On
the other hand, the actual quantity of interest is the computation time for a
given maximal error. With this metric, we found that TEBD in star geometry is
the most favorable algorithm, faster than TDVP in chain geometry by about a
factor of $5$ or more. Combined with its general simplicity and stability,
this makes TEBD in star geometry at present the best approach to solve
impurity problems using real-time evolution.

## References

  * Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * Georges and Kotliar (1992) A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Werner  _et al._ (2006) P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
  * Gull  _et al._ (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Wilson (1975) K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
  * Bulla  _et al._ (2008) R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
  * Lu  _et al._ (2014) Y. Lu, M. Höppner, O. Gunnarsson, and M. W. Haverkort, Phys. Rev. B 90, 085102 (2014).
  * Zgid  _et al._ (2012) D. Zgid, E. Gull, and G. K. L. Chan, Phys. Rev. B 86, 165128 (2012).
  * Mejuto Zaera  _et al._ (2017) C. Mejuto Zaera, N. M. Tubman, and K. B. Whaley, arXiv e-prints (2017), arXiv:1711.04771 [cond-mat.str-el] .
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * Bauernfeind  _et al._ (2017) D. Bauernfeind, M. Zingl, R. Triebl, M. Aichhorn, and H. G. Evertz, Phys. Rev. X 7, 031013 (2017).
  * Bauernfeind  _et al._ (2018) D. Bauernfeind, R. Triebl, M. Zingl, M. Aichhorn, and H. G. Evertz, Phys. Rev. B 97, 115156 (2018).
  * Paeckel  _et al._ (2019) S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, arXiv e-prints (2019), arXiv:1901.05824 [cond-mat.str-el] .
  * Haegeman  _et al._ (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Wolf  _et al._ (2014) F. A. Wolf, I. P. McCulloch, and U. Schollwöck, Phys. Rev. B 90, 235131 (2014).
  * Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
  * Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
  * Daley  _et al._ (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
  * Saberi  _et al._ (2008) H. Saberi, A. Weichselbaum, and J. von Delft, Phys. Rev. B 78, 035124 (2008).
  * Ganahl  _et al._ (2015) M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Held, and F. Verstraete, Phys. Rev. B 92, 155132 (2015).
  * Kennes and Karrasch (2016) D. Kennes and C. Karrasch, Computer Physics Communications 200, 37 (2016).
  * Suzuki (1990) M. Suzuki, Physics Letters A 146, 319 (1990).
  * Note (1) For example a semi-circular bath spectral function can be represented by $t_{i}=t$ for all sites. Note that we use the Lanzcos tri-diagonalization instead.
  * Bauernfeind (2018) D. Bauernfeind,  _Fork Tensor Product States: Efficient Multi-Orbital Impurity Solver for Dynamical Mean Field Theory_ , Ph.D. thesis, Graz University of Technology (2018).
  * Orús (2014) R. Orús, The European Physical Journal B 87, 280 (2014).
  * Verstraete  _et al._ (2009) F. Verstraete, J. I. Cirac, and J. I. Latorre, Phys. Rev. A 79, 032316 (2009).
  * Rams and Zwolak (2019) M. M. Rams and M. Zwolak, arXiv e-prints (2019), arXiv:1904.12793 [cond-mat.str-el] .


# Revisiting Flat bands and localization

Yasuhiro Hatsugai  Department of Physics, University of Tsukuba, 1-1-1
Tennodai, Tsukuba, Ibaraki 305-8571, JAPAN

###### Abstract

Flat bands imply lack of itinerancy due to some constraints that, in
principle, results in anomalous behaviors with randomness. By a molecular
orbital (MO) representation of the flat band systems, random MO models are
introduced where the degeneracy due to the flat bands is preserved even with
randomness. The zero modes of the chiral symmetric system with sublattice
imbalance belong to the class. After explaining the generic flat band
construction by MOs, several examples are discussed with numerical
demonstration as sawtooth lattice in one dimension and hyper-Pyrochlore
lattice in any $d$-dimensions that extends the Kagome ($d=2$) and Pyrochlore
($d=3$) lattices to general dimensions.

###### keywords:

Flat bands, molecular orbitals , Kagome , Pyrochlore , randomness

††journal: Journal of LaTeX Templates††journal: Annals of Physics

##  1 Introduction

Massless Dirac fermions with singular dispersion are sources of non-trivial
topology and get focused substantially in relation to topological phases and
Anderson localization as well. Flat bands are another singular dispersion and
have a long history of studies. Any projection operator of an eigen state in
momentum space has degenerate eigen values (0 and 1), that gives a trivial
example of the flat band. However, constructing a tight-binding Hamiltonian
with strictly short range hopping is not trivial. A tight-binding Hamiltonian
of a simplified $sp^{3}$ network of an amorphous solid is one of the oldest
examples of the flat band systems[1]. The origin of the flat bands in the
model was clarified based on the molecular orbital (MO) construction[2] and
the MO representation is applied to the simplified model of silicene as a two-
dimensional analogue[3]. The systems with flat bands also appeared
occasionally in various studies from quite different points of view such as
fermion doubling[4], ferromagnetism[5, 6, 7] and ferrimagnetism[8]. The flat
band is also observed experimentally in frustrated Kagome metal[9]. The flat
bands of local Hamiltonian are due to some local constraints, which imply
possible non trivial behaviors in Anderson localization such as
multifractality[10, 11]. One of the well known origin of the flat bands is
imbalance of the sublattices in a chiral symmetric system. Appearance of
localized boundary modes of graphene near the Zigzag edges belong also to the
class. A hopping Hamiltonian on the Penrose tiling is also in this class since
the Penrose tilling is bipartite and there exist sublattices imbalance[12].
Although the system is aperiodic, the Hamiltonian hosts macroscopic degeneracy
(ring states) due to the constraint. This is an example of ”flat band” without
translational symmetry.

Recently we are proposing a general scheme to construct a strictly local
Hamiltonian with flat bands [2, 3, 13, 14]. It includes the above chiral
symmetric sublattice imbalance case. The construction is independent of the
translational symmetry and can be applied to random systems. We are here
introducing such a class of systems as a random MO model and present some of
numerical results as a demonstration. As for the numerical evaluation of the
probability density $p_{i}$ and its correlation function $C_{ij}$, we do not
need diagonalizaton. Also the numerical studies can be done using a sparse
matrix technique since the MO is local. It enables us to treat large systems.

##  2 Molecular Orbital representation

Let us first define $M$ molecular orbitals (MOs) $|i)$, ($i=1,\cdots,M$) as

$\displaystyle|i)=\sum_{n=1}^{N}|n\rangle\psi_{n,i}=|{\cal
A}\rangle\psi_{i},\quad|{\cal A}\rangle=(|1\rangle,\cdots,|N\rangle),$

where $\psi_{i}$ is an $N$-dimensional column vector of the $i$-th molecular
orbital and $|{\cal A}\rangle$ is a set of orthonormalized atomic basis,
$\langle n|m\rangle=\delta_{nm}$, $\langle{\cal A}|{\cal A}\rangle=I_{N}$ ($N$
is the number of atoms).

A Hamiltonian written by the $M$ molecular orbitals (MOs) is [2, 3, 13, 14]

$\displaystyle{\cal H}$ $\displaystyle=\sum_{i,j}^{M}|i)h_{i,j}(j|=|{\cal
M}){h}({\cal M}|=|{\cal A}\rangle{H}\langle{\cal A}|$ $\displaystyle H$
$\displaystyle=\Psi h\Psi^{\dagger},$

where $H\in M(N,N)$, ${h}\in M(M,M)$ 111 $M(n,m)$ is a set of $n\times m$
matrices. and $h_{i,j}=\\{h\\}_{i,j}$, is a hopping between the molecular
orbitals $|i)$ and $|j)$. Here $|{\cal M})$ is a set of the MOs’ basis

$\displaystyle|{\cal M})$ $\displaystyle=(|1),\cdots,|M))=|{\cal
A}\rangle\Psi,\quad\Psi=(\psi_{1},\cdots,\psi_{M})\in M(N,M).$

We assume that $\det h\neq 0$ and the MOs $\psi_{i}$, $i=1,\cdots M$ are
linearly independent (${\rm rank}\,\Psi=M$) 222If $\det h=0$, one can reduce
the independent MOs. . It implies the overlap matrix ${\cal O}=({\cal M}|{\cal
M})=\Psi^{\dagger}\Psi$ is positive definite ($\det O\neq 0$).

If $N>M$, this Hamiltonian ${\cal H}$ has at least $N-M$ degenerate zero modes
since the zero energy condition, $({\cal M}|\varphi\rangle=0$ for
$|\varphi\rangle=|{\cal A}\rangle\varphi$, gives $M$ linear conditions for
$N$-dimensional vector $\varphi$. 333 It also obeys from a simple algebra for
the secular equation $\det_{N}(\lambda I_{N}-H)=\det_{N}(\lambda I_{N}-\Psi
h\Psi^{\dagger})=\lambda^{N-M}\det_{M}(\lambda I_{M}-h{\cal O})$ [2, 13, 14].
444 When the system is periodic, there exists flat bands. One may also
consider $H$ as a Hamiltonian in momentum representation and the label $n$ as
an index of the atomic species within the unit cell in momentum space. Then
the zero energy flat band is $N-M$ fold degenerated.  555 This MO construction
of the zero mode flat bands reminds us of the Maxwell relations for
stability/instability of mechanical systems[15, 16].  This is true even when
$\Psi$ is random, which we propose as a random molecular orbital model[2, 17,
18].

By using an orthnormalized basis of the MOs,

$\displaystyle|\bar{\cal M})$ $\displaystyle=|{\cal M}){\cal O}^{-1/2},\
(\bar{\cal M}|\bar{M})=I_{M},$

the projection to the non-zero energy states, ${\cal P}_{1}$ is written as

$\displaystyle{\cal P}_{1}$ $\displaystyle=|\bar{\cal M})(\bar{\cal M}|=|{\cal
A}\rangle\Psi{\cal O}^{-1}\Psi^{\dagger}\langle{\cal A}|.$

Now we have an orthogonal decomposition of the one-particle Hilbert space as

$\displaystyle 1$ $\displaystyle={\cal P}_{0}+{\cal P}_{1},\ \ {\cal
P}_{i}{\cal P}_{j}=\delta_{ij}{\cal P}_{i},$

where ${\cal P}_{1}$ and ${\cal P}_{0}$ are projections to the non-zero/zero
energy states respectively (${\cal H}{\cal P}_{1}={\cal H}$ and ${\cal H}{\cal
P}_{0}=0$). Then the complementary projection to the zero modes is written as

$\displaystyle{\cal P}_{0}$ $\displaystyle=1-{\cal P}_{1}=|Z\rangle\langle
Z|,\ \langle Z|Z\rangle=I_{N-M}$ $\displaystyle|Z\rangle$
$\displaystyle=(|\varphi_{1}\rangle,\cdots,|\varphi_{N-M}\rangle),\ \langle
M|\varphi_{i}\rangle=0$ $\displaystyle|\varphi_{\ell}\rangle$
$\displaystyle=|{\cal A}\rangle\varphi_{\ell},\ H\varphi_{\ell}=0,$

where $|\varphi_{\ell}\rangle$, $\ell=1,\cdots,N-M$ are orthonormalized states
of the zero mode,
$\langle\varphi_{\ell}|\varphi_{\ell^{\prime}}\rangle=\delta_{\ell\ell^{\prime}}$.

To discuss properties of localization of the wave function $\varphi_{\ell}$ of
the zero energy, the probability distribution
$p_{i}^{\ell}=C^{\ell}_{ii}=|\varphi_{i,\ell}|^{2}$ and the correlation
function $C^{\ell}_{ij}=\varphi_{i,\ell}\varphi_{j,\ell}^{*}$ can be useful.
However, the exact degeneracy of the flat band implies that each
$\varphi_{\ell}$ is not uniquely defined and one may take another set of the
zero energy states,
$|Z^{\prime}\rangle=(|\varphi_{1}^{\prime}\rangle,\cdots,|\varphi_{N-M}^{\prime}\rangle$.
Of course $|Z\rangle$ and $|Z^{\prime}\rangle$ are unitary equivalent as
$|Z^{\prime}\rangle=|Z\rangle U$, $U\in U(N-M)$. Then what is invariant
(independent of the transformation $U$) is an average or a trace over the
degenerate states defined as

$\displaystyle C_{ij}$
$\displaystyle={\color[rgb]{0,0,0}\sum_{\ell=1}^{N-M}C^{\ell}_{ij}=}\sum_{\ell=1}^{N-M}\varphi_{{\color[rgb]{0,0,0}i,\ell}}\varphi_{{\color[rgb]{0,0,0}j,\ell}}^{*}=\langle
i|{\cal P}_{0}|j\rangle=\delta_{ij}-\tilde{C}_{ij}$
$\displaystyle\tilde{C}_{ij}$ $\displaystyle=(\Psi{\cal
O}^{-1}\Psi^{\dagger})_{ij}.$

Its invariance is clear by ${\cal P}_{0}=|Z\rangle\langle Z|=|Z\rangle\langle
Z^{\prime}|$.  The local probability of the site $i$ for the zero mode space
defined as $p_{i}=C_{ii}/M$, ($\sum_{i}p_{i}=1$) is the average of the local
probability distribution of $p_{i}^{\ell}$. This is the only well-defined
quantity that is independent of the choice of the states. One may wonder
effects of the randomness might be washed out by the average. However it is
not the case as shown in the concrete examples discussed in the next section.
Since the projection into the zero energy space ${\cal P}_{0}$ is independent
of the hopping between the MOs, $h_{ij}$, we assume $h=I_{M}$ without loss of
generality as far as the degenerate zero energy space is concerned.

Note also that $C_{ij}$ is also written as a correlation function of the
filled zero modes $|z\rangle_{F}$ of the second quantized fermions as

$\displaystyle C_{ij}$
$\displaystyle={{}_{F}\langle}z|c_{j}^{\dagger}c_{i}|z\rangle_{F}$
$\displaystyle|z\rangle_{F}$
$\displaystyle=\prod_{\ell=1}^{N-M}\bm{c}^{\dagger}\varphi_{\ell}|0\rangle{\color[rgb]{0,0,0}{}_{F}},$

since
${{}_{F}\langle}z|c_{i}c_{j}^{\dagger}|z{\rangle_{F}}=\delta_{ij}-\langle
i|{\cal P}_{0}|j\rangle$ where $c_{i}$, ($i=1,\cdots N$) is a fermion
annihilation operator (${\\{}c_{i},c_{j}^{\dagger}{\\}}=\delta_{ij}$),
$\bm{c}^{\dagger}=(c_{1},\cdots,c_{N})$ and
$|0\rangle{\color[rgb]{0,0,0}{}_{F}}$ is their vacuum
($c_{i}|0\rangle{\color[rgb]{0,0,0}{}_{F}}=0$).  See appendix A.

##  3 Examples

###  3.1 Sawtooth lattice in 1D

Figure 1:  (a): Sawtooth lattice and its energy dispersion for $a=b=c=1$. (b)
Integrated density of states $N(\epsilon)$ ($\sigma=0.2$) (c) ”Charge”
$n_{i}=C_{ii}$ of the zero energy space ($N=2M=8192$) with periodic boundary
condition. Normalized probability density $p_{i}=n_{i}/M$. (d) ”Correlation
function” $\tilde{C}_{i,i+j}$ of the zero energy space. The black lines are
for 21 different $i$’s and the red one is their average.

The simplest example with a flat band is a one-dimensional sawtooth lattice (
Fig.1(a)) This model has been discussed in many papers (See.[19], for
example). The Hamiltonian is written by $M$-MO,
$|j)=|\xi_{j}\rangle\xi_{j}^{1}+|\eta_{j}\rangle\eta_{j}+|\xi_{j+1}\rangle\xi_{j}^{2}$,
($j=1,\cdots,M)$ for an $N=2M$ site system as $H=\sum_{j=1}^{M}|j)(j|$ where
$|\xi_{j}\rangle$ and $|\eta_{j}\rangle$ are orthonormalized atomic states and
$\xi_{j}^{1},\xi_{j}^{2},\eta_{j}\in\mathbb{C}$ are arbitrary. We assume
periodic boundary condition $|\xi_{M+1}\rangle=|\xi_{1}\rangle$ for
simplicity. Since the Hamiltonian is written by $M$ MOs and the number of the
atomic sites is $N=2M$, it has $N-M=M$ degenerate zero energy states. If the
system is translational invariant as $\xi_{j}^{1}=a,\xi_{j}^{2}=b,\eta_{j}=c$,
one has $H=\sum_{k_{n}}|k_{n})(k_{n}|$,
$|k_{n})=(|\xi_{k_{n}}\rangle,|\eta_{k_{n}}\rangle)\psi_{k_{n}}$
${}^{t}\psi_{k_{n}}=(a+ce^{ik_{n}},b)$ in the momentum representation
$|\xi_{j})=M^{-1/2}\sum_{k_{n}}|\xi_{k_{n}})e^{ik_{n}j}$,
$|\eta_{j})=M^{-1/2}\sum_{k_{n}}|\eta_{k_{n}})e^{ik_{n}j}$, $k_{n}=\frac{2\pi
n}{M},n=1,\cdots,M$. Now its non-zero energy band is given by
$\psi_{k_{n}}^{\dagger}\psi_{k_{n}}=|a+be^{ik_{n}}|^{2}+|c|^{2}$ and the rest
is the flat band at the zero energy.

We have calculated probability distribution $p_{i}$ and its correlation
function $C_{ij}$ for the random MO sawtooth model in Fig.1(b)-(d). The random
MOs are defined by $\xi_{i}^{1,2}=1+X_{i}^{1,2}$ and $\eta_{i}=1+X_{j}$ where
$X_{i}^{1,2}$ and $X_{i}$ are random variables of the normal distribution with
variance $\sigma^{2}$ and its mean $0$. In the Fig.1(b), the integrated
density of states is shown for $\sigma=0.2$. The charge $n_{i}$ is shown for
$\sigma=0.2$ for $N=8192$ system in Fig.1. It is very spiky and behaves
randomly in a microscopic scale. We have also calculated inverse participation
ratios (IPR), $I=\sum_{i}p_{i}^{4}$, ($p_{i}=n_{i}/M$, $\sum_{i}p_{i}=1$) for
different system sizes ($N=8$ to $N=8192$). It behaves ${\cal O}(1/N)$ that
suggests the space is extended since it is an average of the IPR of the
degenerated zero energy states. Also the correlation function
$\tilde{C}_{i,i+j}$ is plotted in Fig.1(d) averaged over several $i$’s. It
suggests the average $\tilde{C}_{ij}$ is exponentially localized, that is, the
zero energy space is short-ranged and its localization length is roughly
several lattice spacing. It implies that the zero mode space is intrinsically
localized due to the local nature of the MO within microscopic space and the
macroscopic degeneracy makes the whole zero mode space extended as is shown in
the IPR.

###  3.2 Hyper-Pyrochlore lattice in $d$-dimensions

Figure 2:  (a)Unit vectors $e_{i}$, $j=1,\cdots,d$ for a series of hyper-
Pyrochlore lattices[2]. Equivalent $(d+1)$ points $E^{(d)}_{j}$,
$j=1,\cdots,d+1$ are distributed on the Sphere $S^{d-1}$. (b) Several lattice
points and its unit vectors are shown for $d=2$ (Kagome) and $d=3$
(Pyrochlore).

We have also proposed a series of systems with flat bands in $d$-dimensions as
a hyper-Pyrochlore lattice[2]. This is dual to the $d$-dimensional graphene
that is an extension of 4-$d$ graphene[20]. Energy bands of the hyper-
Pyrochlore are given by the (shifted) bands of $d$-dimensional graphene with
extra $(d-1)$-fold degenerated flat bands at zero energy[2]. As the
$d$-dimensional graphene is a minimum model for the massive/massless Dirac
fermions in $d$-dimension, the hyper-Pyrochlore lattice is also minimum as a
system with flat bands in general dimension.

Let us start from $(d-1)$-dimensional sphere $S^{d-1}$ (See Fig.2(a)) and
$d+1$ points, $E^{(d)}_{1},\cdots,E^{(d)}_{d+1}$ which are equivalently
distributed on the unit sphere ($E^{d}_{j}\cdot E^{(d)}_{j^{\prime}}=-d^{-1}$,
$j=1,\cdots,d+1$). 666They are recursively constructed from $d=1$ by
$\displaystyle E^{(d)}_{j}$
$\displaystyle=\left\\{\begin{array}[]{cc}\alpha_{d}E_{j}^{(d-1)}+\beta_{d}E^{(d)}_{d+1}&j=1\cdots,d\\\
^{t}(0,\cdots,0,1_{d},0,\cdots)&j=d+1\end{array}\right.,$ where
$\alpha_{d}=\sqrt{1-d^{-2}}$ and $\beta_{d}=-d^{-1}$. Here we have assumed the
system is embedded in sufficiently large dimension.  Then the unit
translations, $e_{1},\cdots,e_{d}$, are given by
$e_{j}=e^{(d)}_{j}=E^{(d)}_{j}-E^{(d)}_{d+1}$. The unit cell of the hyper-
Pyrochlore lattice includes $d+1$ atoms which locate at
$-e_{1}/2,\cdots,-e_{d+1}/2$, ($e_{d+1}=0$ is supplemented).

Figure 3:  (a) Integrated density of states for Kagome lattice without
randomness ($\xi^{j}_{r}=\eta^{j}_{r}=1$). (b) ”Correlation function”
$\tilde{C}_{i,i+j}$ of the zero energy space ($N=3\times 32^{2}$). (c)
”Charge” $n_{i}=C_{ii}$ of the zero energy space ($N=3M=3\times 256^{2}$) with
periodic boundary condition. Normalized probability density $p_{i}=n_{i}/M$.
($\sigma=0.3$) (d) ”Charge” $n_{i}=C_{ii}$ of the zero energy space
($N=3M=3\times 256^{2}$) with periodic boundary condition. Normalized
probability density $p_{i}=n_{i}/M$. ($\sigma=3.0$)

With periodic boundary condition of the linear size $L$, $(d+1)L^{d}$ atomic
positions of the system with $L^{d}$ unit cells are

$\displaystyle r_{j}$ $\displaystyle=r-e_{j}/2,\ (j=1,\cdots,d+1),$

where $r=n_{1}e_{1}+\cdots+n_{d}e_{d}$ and $n_{j}=1,\cdots,L$,
($j=1,\cdots,d$). See Fig.2(b). The Hamiltonian of the hyper-Pyrochlore is
written as

$\displaystyle H$ $\displaystyle=\sum_{r}(|B_{r})(B_{r}|+|R_{r})(R_{r}|).$

Two MOs, $|B_{r})$ and $|R_{r})$ at the unit cell $r$ are defined as

$\displaystyle|B_{r})$
$\displaystyle=\sum_{j=1}^{d+1}|r_{j}\rangle\xi_{r}^{j}$
$\displaystyle|R_{r})$
$\displaystyle=\sum_{j=1}^{d+1}|r+e_{j}\rangle\eta_{r}^{j}{\color[rgb]{0,0,0},}$

where $|r_{j}\rangle$ and $|r_{j}+e_{j}\rangle$ are orthonormalized atomic
states at the sites. The weights (coefficients) within the MO,
$\xi_{r}^{j},\eta_{r}^{j}\in\mathbb{C}$, can be any. As for the numerical
demonstration, we take $\xi_{r}^{j}=1+X_{r}^{j}$ and
$\eta_{r}^{j}=1+Y_{r}^{j}$ where $X_{r}^{j}$ and $Y_{r}^{j}$ are random
variables of the normal distribution with variance $\sigma^{2}$ and its means
$0$. Since the number of the atomic sites is $N=(d+1)L^{d}$ and that of the
MOs is $M=2L^{d}$, the system has $N-M=(d-1)L^{d}$ zero energies. If the
system is translationally invariant, that is, $\xi_{r}^{j}$ and $\eta_{r}^{j}$
are $r$-independent, the system has (at least) $(d-1)$-fold degenerate flat
bands (at the zero energy).  In a translationally invariant clean system, the
Hamiltonian in momentum space is written by the MOs in momentum space. Then
the flat band is give by the momentum independent zero mode(see footnote in
page 2). In this case, the MOs are momentum dependent and may not be always
linearly independent. If the MOs are not independent, $\det{\cal O}=0$, it
causes extra degeneracies such as band touching[21, 2, 13].

Numerical results for the case $d=2$, that is, the random MO Kagome lattice is
shown in Fig.3. They are comparable to the results of the one-dimensional
sawtooth lattice. In Fig.3(a), the integrated density of states for the random
MO Kagome model is shown ($\sigma=0.2$). The Dirac cones for the clean case is
also shown for reference. The “correlation function” $\tilde{C}_{i,i+j}$ of
the zero energy space is shown in Fig.3(b) for a rather small system
($N=3\times 32^{2}$) $i$ is at the center of the system. It also shows the
correlation function is short-ranged over several lattice spacing. Charge
distributions/probability distribution $n_{i}$ are shown in Fig.3(c),(d) for
$\sigma=0.3$ and $\sigma=3.0$ respectively ($N=3\times 256^{2}$). They show
the distribution is very spiky over several to 10 lattice spacing. We have
also calculated IPR, $I=\sum_{i}p_{i}^{4}$, ($p_{i}=n_{i}/M$,
$\sum_{i}p_{i}=1$) for $N=3\times L^{2}$, $L=2^{4}$ to $2^{7}$. This is the
average of the IPR’s over the degenerate zero energy states this is unitary
invariant . They behave ${\cal O}(1/L^{2})$ that is consistent with the
extended states. It suggests the zero mode space is microscopically quite
spiky but extended globally. It is similar to the results of the one-
dimensional sawtooth lattice. 777 When one uses the method in Ref.[22] for
calculating a singularity spectrum (multifractal analysis), after averaging
over a box of the size $\ell$, taking an extrapolation $\ell\to 0$ is needed.
Using the average over the box size $\ell$ up to several lattice spacing,
numerical multifractal analysis of the wave functions of the random MO Kagome
lattice has been performed by exact diagonalizations[18]. It suggests
multifractality of the flat band wave functions in a microscopic scale (over
several lattice spacing). To obtain conclusive results, further analysis is
needed.  Similarly to the results of the sawtooth lattice, the zero mode space
of the random MO Kagome lattice is intrinsically short-ranged (this is also
true without any randomness). It is consistent with the spiky behavior of the
probability distribution within a microscopic scale. However, due to the
macroscopic degeneracy, the space as a whole is extended as shown by the IPR.

At the mobility edge of the Anderson localization as the critical point, one
may naturally expect multifractal behaviors of the wave function since the
localization length diverges and the system is scale free. Away from the
critical point, such a multifractal behavior is cut off at the localization
length. It implies that that the wave function of localized states is
multifractal-like up to the localization length. Even for the present Kagome
flat band system, one may expect standard Anderson localization except the
flat bands. As discussed, the zero mode flat band space, that is, the average
of the flat bands are complementary to the sum of them those are localized and
their localization centers are randomly chosen. The spiky structure observed
in $C_{ij}$ is reflecting this multifractal-like behavior within the
localization length.

###  3.3 Bipartite lattice

If the lattice sites are divided into two sublattices A and B and the hopping
matrix elements are only non-zero between them, the lattice is bipartite and
the corresponding Hamiltonian written in the atomic basis is chiral symmetric.
By the choice of the MO in Ref.[13], it is written by the MOs as

$\displaystyle{\cal H}_{C}$ $\displaystyle=\sum_{a\in A,b\in
B}(|a\rangle_{A}t_{ab}{{}_{B}\langle}b|+|b\rangle_{B}t_{ab}^{*}{{}_{A}\langle}a|)$
$\displaystyle=|{{\cal A}}\rangle D_{AB}\langle{\cal B}|+|{\cal B}\rangle
D_{BA}\langle{\cal A}|$ $\displaystyle=|{\cal A},{\cal B}\rangle
H_{C}\langle{\cal A},{\cal B}|$ $\displaystyle H_{C}$
$\displaystyle=\left(\begin{array}[]{cc}O_{AA}&D_{AB}\\\
D_{BA}&O_{BB}\end{array}\right)$ $\displaystyle|{\cal A},{\cal B}\rangle$
$\displaystyle=(|{\cal A}\rangle,|{\cal B}\rangle),\ \ |{\cal
A}\rangle=(|1\rangle_{A},\cdots,|N_{A}\rangle_{A}),\ \ |{\cal
B}\rangle=(|1\rangle_{B},\cdots,|N_{B}\rangle_{B}),$

where $D_{BA}=D_{AB}^{\dagger}\in M(N_{A},N_{B})$, $\\{D_{AB}\\}_{ab}=t_{ab}$
and $H_{C}$ is chiral symmetric ${\\{}H_{C},\Gamma{\\}}=0$ with $\Gamma={\rm
diag}\,(I_{A},-I_{B})$. It is also written by $M=2N_{B}$ MOs, $\Psi$ as [13]
888 $O_{\alpha\beta}\in M(N_{\alpha},N_{\beta})$ is a zero matrix
($\alpha,\beta\in\\{A,B\\}$).

$\displaystyle{\cal H}_{C}$ $\displaystyle=|{\cal
M})h({\color[rgb]{0,0,0}{\cal M}}|{\color[rgb]{0,0,0},}\ |{\cal M})=(|{\cal
A}\rangle D_{AB},\ |{\cal B}\rangle)=|{\cal A},{\cal B}\rangle\Psi$
$\displaystyle\Psi$ $\displaystyle=\left(\begin{array}[]{cc}D_{AB}&O_{AB}\\\
O_{BB}&I_{N_{B}}\end{array}\right)\in
M(N_{A}+N_{B},2N_{B}){\color[rgb]{0,0,0},}$ $\displaystyle h$
$\displaystyle=\left(\begin{array}[]{cc}O_{BB}&I_{N_{B}}\\\
I_{N_{B}}&O_{BB}\end{array}\right)\in M(2N_{B},2N_{B}).$

Since the number of total atomic sites is $N=N_{A}+N_{B}$, the Hamiltonian has
$Z=N-M=N_{A}-N_{B}$ zero modes if $N_{A}>N_{B}$. (The case, $N_{B}>N_{A}$, is
discussed similarly.) It implies that the chiral symmetric (random) system is
regarded as the (random) MO model. One of the historical examples of the flat
band due to this mechanism is the ring states of the Penrose tiling[12] and
discussed from a generic point of view[23]. Another typical example is the
zero energy state of the random Dirac fermions[24], where $N_{A}-N_{B}=1$ and
the zero mode wave function is discussed in relation to its multifractality.

##  4 Discussion

As has been discussed previously[10, 11], the states near the flat band energy
become singular by randomness and their multifractality has been implied. In
the random MO models presented in examples of the sections 2 and 3,
macroscopic degeneracy of the flat bands is preserved even with randomness.
This macroscopic degeneracy makes the situation unclear although the zero mode
space obtained numerically is very local and singular in a microscopic scale.
It is also consistent with the numerical calculation in Ref.[24] where the
zero mode is not degenerated ($Z=N-M=1$). Then introducing extra MOs
associated with the randomness can be interesting. If one introduces
macroscopic $M^{\prime}$ extra random MOs ($M+M^{\prime}\lesssim N$), the
dimension of the zero energy space can be ${\cal O}(1)$. Also the MO
representation are useful for diagonalization of random systems within the
projected flat band subspace.

##  5 Acknowledgment

We thank discussion with T. Mizoguchi and T. Kuroda. This work is supported in
part by the JSPS KAKENHI, Grant No. JP17H06138.

##  Appendix A Fermion correlation function

Let us here discuss fermion correlation functions for a filled fermi state. We
label the Fock state $|\Psi\rangle_{F}$ specified by the multiplet, $\Psi\in
M(N,M)$, of the dimension $M$ for the dimension $N$ Hilbert space as

$\displaystyle|\Psi\rangle_{F}$
$\displaystyle=\prod_{\ell=1}^{M}({c}^{\dagger}\psi_{\ell})|0\rangle_{F}=({c}^{\dagger}\psi_{1})\cdots({c}^{\dagger}\psi_{M})|0\rangle_{F}=\prod_{\ell=1}^{M}(c_{j}^{\dagger}\psi^{j}_{\ell})|0\rangle_{F}$
$\displaystyle\Psi$ $\displaystyle=(\psi_{1},\cdots,\psi_{M}),\
\psi_{\ell}=\left(\begin{array}[]{c}\psi_{1,\ell}\\\ \vdots\\\
\psi_{N,\ell}\end{array}\right),\ {c}=\left(\begin{array}[]{c}c_{1}\\\
\vdots\\\ c_{N}\end{array}\right),$

where $|0\rangle_{F}$ a vacuum, $c_{i}|0\rangle_{F}=0$, and the ordering of
the fermions are implicitly assumed as shown above. The overlap of such two
fermionic states $|\Phi\rangle_{F}$ and $|\Psi\rangle_{F}$ is directly
evaluated as

$\displaystyle{{}_{F}\langle}\Psi|\Phi{\rangle_{F}}$
$\displaystyle=\det\nolimits_{M}\Psi^{\dagger}\Phi.$

Using this, the correlation function,
$g^{U}_{ab}={{}_{F}\langle}\Psi|c_{a}c_{b}^{\dagger}|\Phi{\rangle_{F}}$, is
written as

$\displaystyle g^{U}_{ab}$
$\displaystyle={{}_{F}\langle}\Psi|c_{a}c_{b}^{\dagger}|\Phi\rangle_{F}=\det\nolimits_{M+1}\Psi_{a}^{\dagger}\Phi_{b}$

where $\Phi_{a}=(e_{a},\Phi)\in M(N,M+1)$, $\Psi_{b}=(e_{b},\Psi)\in M(N,M+1)$
are obtained from $\Phi$ and $\Psi$ by adding column vectors $e_{a}$ and
$e_{b}$ respectively where $\\{e_{i}\\}_{j}=\delta_{ij}$.

Using a formula for the determinant 999 As for $A\in M(N,N)$, $B\in M(N,M)$,
$C\in M(M,N)$, $D\in M(M,M)$ and $\det_{M}D\neq 0$, the following relation
holds,  $\displaystyle\det\nolimits_{N+M}\left(\begin{array}[]{cc}A&B\\\
C&D\end{array}\right)$
$\displaystyle=\det\nolimits_{N+M}\left(\begin{array}[]{cc}I_{N}&-BD^{-1}\\\
O_{MN}&I_{M}\end{array}\right)\left(\begin{array}[]{cc}A&B\\\
C&D\end{array}\right)=\det\nolimits_{N+M}\left(\begin{array}[]{cc}A-BD^{-1}C&O_{NM}\\\
C&D\end{array}\right)$
$\displaystyle=\det\nolimits_{N}(A-BD^{-1}C)\det\nolimits_{M}{D}$ , it is
calculated as

$\displaystyle g^{U}_{ab}$ $\displaystyle=$
$\displaystyle\det\nolimits_{M+1}\left(\begin{array}[]{c}e_{a}^{\dagger}\\\
\Psi^{\dagger}\end{array}\right)(e_{b},\Phi)=\det\nolimits_{M+1}\left(\begin{array}[]{cc}e_{a}^{\dagger}e_{b}&e_{a}^{\dagger}\Phi\\\
\Psi^{\dagger}e_{b}&\Psi^{\dagger}\Phi\end{array}\right)$ $\displaystyle=$
$\displaystyle\det\nolimits_{M+1}\left(\begin{array}[]{cc}e_{a}^{\dagger}e_{b}-e_{a}^{\dagger}\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}e_{b}&0\\\
\Psi^{\dagger}e_{b}&\Psi^{\dagger}\Phi\end{array}\right)$ $\displaystyle=$
$\displaystyle\det\nolimits_{1}[e_{a}^{\dagger}e_{b}-e_{a}^{\dagger}\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}e_{b}]\det\nolimits_{M}(\Psi^{\dagger}\Phi)$
$\displaystyle=$
$\displaystyle\det\nolimits_{1}[e_{a}^{\dagger}\big{[}I_{N}-\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}\big{]}e_{b}]\det\nolimits_{M}(\Psi^{\dagger}\Phi)$
$\displaystyle=$
$\displaystyle\big{[}I_{N}-\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}\big{]}_{ab}\det\nolimits_{M}(\Psi^{\dagger}\Phi).$

It implies

$\displaystyle g_{ab}$
$\displaystyle\equiv\frac{{{}_{F}\langle}\Psi|c_{a}c_{b}^{\dagger}|\Phi{\rangle_{F}}}{{{}_{F}\langle}\Psi|\Phi{\rangle_{F}}}=g^{U}_{ab}/\det\nolimits_{M}(\Psi^{\dagger}\Phi)=[g]_{ab}$
$\displaystyle g$
$\displaystyle=I_{N}-\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}$
$\displaystyle\tilde{g}_{ba}$
$\displaystyle=\frac{{{}_{F}\langle}\Psi|c_{b}^{\dagger}c_{a}|\Phi{\rangle_{F}}}{{{}_{F}\langle}\Psi|\Phi{\rangle_{F}}}=\big{[}\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}\big{]}_{ab}$

Similarly the 4-point correlation function is evaluated as

$\displaystyle g^{U}_{ba;cd}$
$\displaystyle={{}_{F}\langle}\Psi|c_{b}c_{a}c_{c}^{\dagger}c_{d}^{\dagger}|\Phi{\rangle_{F}}$
$\displaystyle=\det\nolimits_{M+2}\left(\begin{array}[]{c}e_{ab}^{\dagger}\\\
\Psi^{\dagger}\end{array}\right)(e_{cd},\Phi),\quad e_{ab}=(e_{a},e_{b}),\
e_{cd}=(e_{c},e_{d})$
$\displaystyle=\det\nolimits_{2}[e_{ab}^{\dagger}e_{cd}-e_{ab}^{\dagger}\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}e_{cd}]\det\nolimits_{M}(\Psi^{\dagger}\Phi)$
$\displaystyle=\det\nolimits_{2}[e_{ab}^{\dagger}\big{[}I_{N}-\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}\big{]}e_{cd}]\det\nolimits_{M}(\Psi^{\dagger}\Phi)$
$\displaystyle=\det\nolimits_{2}\bigg{[}\left(\begin{array}[]{c}e_{a}^{\dagger}\\\
e_{b}^{\dagger}\end{array}\right)\big{[}I_{N}-\Phi(\Psi^{\dagger}\Phi)^{-1}\Psi^{\dagger}\big{]}({e_{c}},{e_{d}})\bigg{]}\det\nolimits_{M}(\Psi^{\dagger}\Phi)$

It implies the Wick theorem

$\displaystyle g_{ba;cd}$ $\displaystyle\equiv$
$\displaystyle\frac{{{}_{F}\langle}\Psi|c_{b}c_{a}c_{c}^{\dagger}c_{d}^{\dagger}|\Phi{\rangle_{F}}}{{{}_{F}\langle}\Psi|\Phi{\rangle_{F}}}=g^{U}_{ab;cd}/\det\nolimits_{M}(\Psi|\Phi)$
$\displaystyle=$
$\displaystyle\det\nolimits_{2}\left(\begin{array}[]{c}e_{a}^{\dagger}\\\
e_{b}^{\dagger}\end{array}\right)g({e_{c}},{e_{d}})$ $\displaystyle=$
$\displaystyle g_{ac}g_{bd}-g_{ad}g_{bc}.$

A generic case for the $2M$-point function follows straightforwardly as

$\displaystyle g_{a_{M}\cdots a_{2}a_{1};b_{1}b_{2}\cdots b_{M}}$
$\displaystyle\equiv\frac{{{}_{F}\langle}\Psi|c_{a_{M}}\cdots
c_{a_{1}}c^{\dagger}_{b_{1}}\cdots
c^{\dagger}_{b_{M}}|\Phi{\rangle_{F}}}{{{}_{F}\langle}\Psi|\Phi{\rangle_{F}}}$
$\displaystyle=\det\nolimits_{M}e_{a_{1},\cdots,a_{M}}^{\dagger}ge_{b_{1},\cdots,b_{M}}$
$\displaystyle=\det\nolimits_{M}g^{R},$

where $e_{i_{1},\cdots,i_{M}}=(e_{i_{1}},\cdots,e_{i_{M}})\in M(N,M)$ and
$\\{g^{R}\\}_{i,j}=g_{a_{i}b_{j}},\ (i,j=1,\cdots,M)$.

## References

  * [1] D. Weaire, M. F. Thorpe, Electronic properties of an amorphous solid. i. a simple tight-binding theory, Phys. Rev. B 4 (1971) 2508–2520.  doi:10.1103/PhysRevB.4.2508.
URL https://link.aps.org/doi/10.1103/PhysRevB.4.2508

  * [2] Y. Hatsugai, I. Maruyama, ZQ topological invariants for polyacetylene, kagome and pyrochlore lattices, EPL (Europhysics Letters) 95 (2) (2011).
  * [3] Y. Hatsugai, K. Shiraishi, H. Aoki, Flat bands in the weaire–thorpe model and silicene, New Journal of Physics 17 (2) (2015).
  * [4] E. Dagotto, E. Fradkin, A. Moreo, A comment on the nielsen-ninomiya theorem, Physics Letters B 172 (3) (1986) 383 – 386.  doi:https://doi.org/10.1016/0370-2693(86)90274-1.
URL http://www.sciencedirect.com/science/article/pii/0370269386902741

  * [5] A. Mielke, Ferromagnetic ground states for the hubbard model on line graphs, Journal of Physics A: Mathematical and General 24 (2) (1991) L73–L77.  doi:10.1088/0305-4470/24/2/005.
URL https://doi.org/10.1088/0305-4470/24/2/005

  * [6] H. Tasaki, Ferromagnetism in the hubbard models with degenerate single-electron ground states, Phys. Rev. Lett. 69 (1992) 1608–1611.  doi:10.1103/PhysRevLett.69.1608.
URL https://link.aps.org/doi/10.1103/PhysRevLett.69.1608

  * [7] H. Katsura, I. Maruyama, A. Tanaka, H. Tasaki, Ferromagnetism in the hubbard model with topological/non-topological flat bands, EPL (Europhysics Letters) 91 (5) (2010) 57007.  doi:10.1209/0295-5075/91/57007.
URL https://doi.org/10.1209/0295-5075/91/57007

  * [8] E. H. Lieb, Two theorems on the hubbard model, Phys. Rev. Lett. 62 (1989) 1201–1204.  doi:10.1103/PhysRevLett.62.1201.
URL https://link.aps.org/doi/10.1103/PhysRevLett.62.1201

  * [9] M. Kang, S. Fang, L. Ye, H. C. Po, J. Denlinger, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, J. G. Checkelsky, R. Comin, Topological flat bands in frustrated kagome lattice cosn, Nature Communications 11 (1) (2020) 4004.  doi:10.1038/s41467-020-17465-1.
URL https://doi.org/10.1038/s41467-020-17465-1

  * [10] M. Goda, S. Nishino, H. Matsuda, Inverse anderson transition caused by flatbands, Phys. Rev. Lett. 96 (2006) 126401.  doi:10.1103/PhysRevLett.96.126401.
URL https://link.aps.org/doi/10.1103/PhysRevLett.96.126401

  * [11] J. T. Chalker, T. S. Pickles, P. Shukla, Anderson localization in tight-binding models with flat bands, Phys. Rev. B 82 (2010) 104209\.  doi:10.1103/PhysRevB.82.104209.
URL https://link.aps.org/doi/10.1103/PhysRevB.82.104209

  * [12] M. Kohmoto, B. Sutherland, Electronic states on a penrose lattice, Phys. Rev. Lett. 56 (1986) 2740–2743.
  * [13] T. Mizoguchi, Y. Hatsugai, Molecular-orbital representation of generic flat-band models, EPL (Europhysics Letters).
  * [14] T. Mizoguchi, Y. Hatsugai, Systematic construction of topological flat-band models by molecular-orbital representation, Physical Review B 101 (23) (2020).
  * [15] J. Maxwell, L. on the calculation of the equilibrium and stiffness of frames, Phil. Mag. 27 (1864) 294–299.
  * [16] J. Paulose, B. G.-g. Chen, V. Vitelli, Topological modes bound to dislocations in mechanical metamaterials, Nature Physics 11 (2) (2015) 153–156.  doi:10.1038/nphys3185.
URL https://doi.org/10.1038/nphys3185

  * [17] Y. Hatsugai, K. Kohda, unpublished.
  * [18] K. Kohda, Master thesis, University of Tsukuba (2017) (Japanese).
  * [19] Y. Kuno, T. Mizoguchi, Y. Hatsugai, Flat band quantum scar, Phys. Rev. B 102 (2020) 241115.  doi:10.1103/PhysRevB.102.241115.
URL https://link.aps.org/doi/10.1103/PhysRevB.102.241115

  * [20] M. Creutz, Four-dimensional graphene and chiral fermions, Journal of High Energy Physics 2008 (04) (2008) 017–017.  doi:10.1088/1126-6708/2008/04/017.
URL https://doi.org/10.1088/1126-6708/2008/04/017

  * [21] D. L. Bergman, C. Wu, L. Balents, Band touching from real-space topology in frustrated hopping models, Phys. Rev. B 78 (2008) 125104\.  doi:10.1103/PhysRevB.78.125104.
URL https://link.aps.org/doi/10.1103/PhysRevB.78.125104

  * [22] A. Chhabra, R. V. Jensen, Direct determination of the f($\alpha$) singularity spectrum, Phys. Rev. Lett. 62 (1989) 1327–1330.  doi:10.1103/PhysRevLett.62.1327.
URL https://link.aps.org/doi/10.1103/PhysRevLett.62.1327

  * [23] B. Sutherland, Localization of electronic wave functions due to local topology, Phys. Rev. B 34 (1986) 5208–5211.  doi:10.1103/PhysRevB.34.5208.
URL https://link.aps.org/doi/10.1103/PhysRevB.34.5208

  * [24] Y. Hatsugai, X.-G. Wen, M. Kohmoto, Disordered critical wave functions in random-bond models in two dimensions: Random-lattice fermions ate=0without doubling, Physical Review B 56 (3) (1997) 1061–1064.