Center for Theoretical Physics of Complex Systems, IBS School Doctor of
Philosophy December 2023

# Quantum Phase Transitions
and Dynamics in
Perturbed Flatbands

Sanghoon Lee

Ph.D. Thesis

Quantum Phase Transitions and
Dynamics in Perturbed Flatbands

Sanghoon Lee

Basic Science

UNIVERSITY OF SCIENCE AND TECHNOLOGY

December 2023

Quantum Phase Transitions and
Dynamics in Perturbed Flatbands

Sanghoon Lee

A dissertation Submitted in Partial Fulfillment of Requirements for the Degree
of Doctor of Philosophy

December 2023

UNIVERSITY OF SCIENCE AND TECHNOLOGY

Major of: Basic Science

Supervisor: Alexei Andreanov

Co-supervisor: Sergej Flach

ACKNOWLEDGEMENTS

I want to thank my supervisors, Prof. Dr. Alexei Andreanov and Prof. Dr.
Sergej Flach, for their amazing supervision and support during the Ph.D.
program on the road to attaining the ability as a researcher. I am also very
grateful to all colleagues and administration staff at IBS–PCS for the
valuable discussion and memorable time at the center. Last but not least, I
would like to thank my family for their unconditional and endless support.

ABSTRACT

In recent years, there has been a growing interest in flatband systems which
exhibit macroscopic degeneracies. These systems offer a valuable mathematical
framework for the extreme sensitivity to perturbations and interactions. This
sensitivity unveils a wide variety of exotic and unconventional physical
phenomena. Moreover, the progress in their experimental realization
contributes to the expanding landscape of exploration in this field. This
thesis aims to summarize all the works throughout the Ph.D. program. Firstly,
an in-depth exploration was conducted on the impact of weak quasiperiodic
perturbations on one-dimensional two-band all-bands-flat lattices. These
latticestight-binding Hamiltonian are diagonalized through a sequence of local
unitary transformations. By adjusting the quasiperiodic potential parameters,
the key achievement involves finding a critical-to-insulator transition and
identifying fractality edges in the flatband systems with quasiperiodic
perturbations. Next, the investigation delved into the effects of on-site
interactions among hard-core bosons in one- and two-dimensional cross-stitch
lattices. One key finding is that groundstate energy primarily depends on
compact localized states. Moreover, the presence of barriers of compact
localized states trap bosons, leading to the emergence of non-ergodic
excitation and Hilbert space fragmentation. Lastly, by using an electric
circuit, a compact localized state of the one-dimensional diamond chain was
successfully generated. Lastly, a compact localized eigenstate of the one-
dimensional diamond chain using an electric circuit was successfully generated
via local (linear and non-linear) driving. This achievement opens up a
versatile circuit platform for the generation of flatbands and holds promise
for potential applications in the field of quantum information. I hope these
collective efforts have expanded the frontiers of the field and made a
meaningful contribution to the scientific community. See abstract_kor

###### Contents

  1. 1 Introduction
    1. 1.1 Flatband systems and compact localization
      1. 1.1.1 Flatbands
      2. 1.1.2 Compact localized states
      3. 1.1.3 Construction of flatband systems
    2. 1.2 Experimental realizations
    3. 1.3 Effect of perturbation and interaction
    4. 1.4 Contributed work and outline of the thesis
  2. 2 Generating Critical States from Quasiperturbed Flatbands
    1. 2.1 Introduction
    2. 2.2 Model
      1. 2.2.1 Local unitary transformation and ABF lattice
      2. 2.2.2 Geometric symmetry of ABF lattice
    3. 2.3 Numerical methods
      1. 2.3.1 Inverse participation ratio
      2. 2.3.2 Defining a linear model for localized states
    4. 2.4 Weak perturbation
      1. 2.4.1 Extended Harper model
      2. 2.4.2 Off-diagonal Harper model
      3. 2.4.3 Lyapunov exponent and localization length
      4. 2.4.4 Case 1 – $\beta=0$ and $\lambda_{1}\neq\lambda_{2}$ at $\theta_{1,2}=\pi/4$
      5. 2.4.5 Case 2 – $\lambda_{2}/\lambda_{1}=0$ for any phase difference $\beta$
      6. 2.4.6 Case 3 – $\lambda_{2}=\lambda_{1}$ and $\beta=\pi$
      7. 2.4.7 Remark – Near zero angles $\theta_{1,2}\to 0$
    5. 2.5 Wavepacket spreading in projected model
      1. 2.5.1 Lower bound of diffusion exponent
      2. 2.5.2 Upper bound of diffusion exponent
      3. 2.5.3 Numerical results
    6. 2.6 Finite perturbation
      1. 2.6.1 Representation in semi-engtangled basis
      2. 2.6.2 Critical-to-insulator transition in finite perturbation
      3. 2.6.3 Non-degenerate compact localized states
      4. 2.6.4 Non-simple fractality edges
      5. 2.6.5 Perturbation independent fractality edges for $\lambda_{2}=0$
    7. 2.7 Conclusion
  3. 3 Trapped hard-core bosons in flatband systems
    1. 3.1 Introduction
    2. 3.2 Models
      1. 3.2.1 1D and 2D cross-stitch lattices
      2. 3.2.2 Properties of hard-core bosons
    3. 3.3 CLS groundstates
      1. 3.3.1 $\nu\leq 1/2$ – filling of the flatband CLS
      2. 3.3.2 $\nu\leq 1$ – filling the $a_{n},b_{n}$ dimers
    4. 3.4 Non-ergodic excitation in cross-stitch lattices
      1. 3.4.1 One hard-core boson between two CLSs
      2. 3.4.2 Mapping to spin-1 chain
      3. 3.4.3 Localized excitations in 2D lattices
    5. 3.5 1D diamond lattice and non-ergodic excitation
    6. 3.6 Fluxed 1D diamond lattice and ergodicity
    7. 3.7 Conclusion
  4. 4 Flatband Electric Circuits
    1. 4.1 Introduction
    2. 4.2 Models
      1. 4.2.1 One-dimensional diamond electric lattice
      2. 4.2.2 One-dimensional stub electric lattice
    3. 4.3 Experimental and numerical results
      1. 4.3.1 Results on the one-dimensional diamond lattice
      2. 4.3.2 Nonlinear driving and robustness of diamond CLS
      3. 4.3.3 Results on the one-dimensional stub lattice
    4. 4.4 Conclusion
  5. 5 Final remarks
    1. 5.1 Scientific contribution and summary
    2. 5.2 Possible future research

###### List of Figures

  1. 1.1 Compact localized state of a one-dimensional diamond chain
  2. 2.1 Perturbed ABF network and projected model
  3. 2.2 Localization lengths of extended Harper models
  4. 2.3 Phase diagram at $\theta_{1,2}=\pi/4$
  5. 2.4 More mapping to extended Harper model
  6. 2.5 More mapping to off-diagonal Harper model
  7. 2.6 Correlation of diffusion exponents and averaged IPR exponents
  8. 2.7 Log-log plots of the spreading of a wavepacket
  9. 2.8 Tunable critical-to-insulator transition
  10. 2.10 First example of perturbation-independent fractality edges
  11. 3.1 Energy bands and groundstate energy of 1D cross-stitch lattice
  12. 3.3 Energy bands and groundstate energy of 1D diamond chain
  13. 4.1 One-dimensional diamond electric circuit
  14. 4.2 One-dimensional stub electric circuit
  15. 4.3 Linear driving on the one-dimensional diamond lattice
  16. 4.4 Nonlinear driving on the one-dimensional diamond lattice
  17. 4.5 Linear driving on the stub chain

## Chapter 1 Introduction

P. W. Anderson discovered in 1958 that when an electron is in a disordered
system, the eigenstates can be spatially localized [1]. By employing a tight-
binding approach on a lattice with random on-site potentials, he showed that
an electron ceases to diffuse away from its initial position depending on a
disorder strength, resulting in the material converting from a metal to an
insulator. Such spatial localization is solely due to quantum interference
between various electron paths, $C_{n}$s, with amplitudes, $A_{n}$s,
respectively, from point $\mathbf{r}_{0}$ back to $\mathbf{r}_{0}$. For any
two paths, average of interference terms goes to zero, because
$w_{\mathrm{int}}=2\operatorname{Re}(A_{1}^{*}A_{2})$ may be positive or
negative. Then, we remain only with $w_{\mathrm{cl}}=|A_{1}|^{2}+|A_{2}|^{2}$.
If two paths are almost equal to each other, $A_{1}\approx A_{2}$, and are
related via time-reversal symmetry, constructive interference is obtained.
This explains the reduction in conductivity lead by increasing average return
probability $\langle w\rangle$ [2],

$\displaystyle\langle w\rangle=\langle|A_{1}+A_{2}|^{2}\rangle=\langle
w_{\mathrm{cl}}\rangle+\langle w_{\mathrm{int}}\rangle=\begin{dcases}\langle
w_{\mathrm{cl}}\rangle\text{ },&\text{destructive}\\\ 2\langle
w_{\mathrm{cl}}\rangle\text{ },&\text{constructive}.\end{dcases}$ (1.1)

After Anderson initially proposed the idea of a disorder-induced localization,
many mathematicians dived into the task of rigorously developing Anderson’s
intuition using a functional analysis. One of the main results is the RAGE
(acronym of Ruelle, Amrein, Georgescu, Enss) theorem [3, 4, 5, 6]. It states
that for any eigenstates $\varphi$ corresponding to discrete eigenvalues
(subspace $\mathcal{H}_{\mathrm{pp}}$ of a given Hilbert space $\mathcal{H}$
with pure point spectrum), their time evolved states remain bounded uniformly
in time in the compact ball $\\{|x|<r\\}$ with a radius $r$. In other words,
all wave amplitudes must be zero outside of the compact ball $\\{|x|>r\\}$,

$\displaystyle\lim_{r\to\infty}\sup_{t\in\mathbb{R}}\big{\|}\chi(|x|>r)e^{-iHt}\varphi\big{\|}=0\hskip
5.0pt\text{where }\chi(|x|>r)=\begin{dcases}0\text{ },&|x|\leq r\\\ 1\text{
},&|x|>r\end{dcases}.$ (1.2)

Such bounded states are called to be a _spectral localization_ of eigenstates.
One may define Anderson localization by adding an exponential bound to the
spectral localization so that a given state decays exponentially with a
localization length $\xi$ at $\mathbf{r}$ compared to the center
$\mathbf{r}_{0}$ of the localization,

$\displaystyle|\varphi(\mathbf{r})|\leq
C\exp(-|\mathbf{r}-\mathbf{r}_{0}|/\xi).$ (1.3)

Moreover, the absence of ballistic motion was proven for pure point spectrum
for exponentially localized eigenstates in a disordered system [7],

$\displaystyle\lim_{t\to\infty}\frac{\big{\|}\mathbf{r}e^{-iHt}\varphi\big{\|}^{2}}{t^{2}}=0$
(1.4)

However, the argument above does not always hold for a non-random system with
an arbitrary growth of $\langle\mathbf{r}^{2}\rangle$, even though the
spectrum is pure point with the exponentially localized states as discussed in
Ref. [8, 9]. To have a more natural notion of localization without any
transport, an additional condition of having finite
$\langle|\mathbf{r}|^{p}\rangle$ may be added to have a stronger form of
Anderson localization, for any $p$. Such localization is called a _strong
(dynamical) localization_ [10].

As it was discussed in the above paragraph, Anderson localization arises when
wave functions encounter a random disorder in a system, leading to the
prevention of long-range transport which implies the inability to propagate
freely across the entire system. Now, what we can ask is the following:

_"Is there an alternative approach other than a disordered system
to localize wavefunctions using the quantum mechanical effect?"_

###  1.1 Flatband systems and compact localization

Anderson localization originates from the enhanced return probability of a
particle induced by a random disorder. On the other hand, an alternative way
of enhancing localization can be achieved without a random disorder in some
lattices via destructive interference, known as compact localization. Such
localization is closely connected to _flatbands_ , which refer to constant
energy bands in specific arrangements of translation-invariant tight-binding
structures. This section goes through the modern definition of a flatband
system and a compact localization.

####  1.1.1 Flatbands

The modern understanding of a flatband system in a discrete lattice can be
summarized as follows. In a flatband system, there exists an energy band (or
energy bands) that exhibits a completely flat energy dispersion relation,
$E_{\mathrm{FB}}=E(\mathbf{k})$, with short-range connectivities [11, 12, 13,
14, 15] as well as some long–range connectivities [16]. Such energy band
remains independent of the momentum, resulting in a vanishing group velocity
$\nabla_{\mathbf{k}}E=0$ with effectively infinite mass. This leads to
macroscopic degeneracy at energy $E_{\mathrm{FB}}$ and a strong suppression of
transport. Moreover, as we see soon through examples, the flatband features
the eigenstates trapped over a strictly finite number of sites [11, 17],
dubbed as the _compact localized states_ (CLS), due to destructive
interference caused by network geometry. This remarkable flatness gives rise
to intriguing phenomena and opens up opportunities for unique physical results
via the _enhancement_ of perturbation and interaction effects (see Sec. 1.3),
no matter how weak they are.

Flatband energy can also manifest in continuous systems, a concept distinct
from the flatband systems discussed in this thesis. One such example is the
Landau level [18]. It describes the movement of a non-interacting charged
particle in a two-dimensional spatial domain in the presence of a magnetic
field $\mathbf{B}$ in $z$-axis. In the Landau gauge, a vector potential
$(-By,0,0)$ is a possible option. With this specific gauge, the Hamiltonian
can be expressed as follows which as a resemblance to the quantum harmonic
oscillator, shifted in the coordinate space. Then, the energy levels in the
system are equivalent to those found in the conventional quantum harmonic
oscillator,

$\displaystyle\mathcal{H}=\frac{1}{2m}\left(\hat{p}_{x}-qB\hat{y}\right)^{2}+\frac{\hat{p}_{y}^{2}}{2m}=\hslash\omega_{c}(\hat{n}+1/2).$
(1.5)

Notably, the energy remains independent of the quantum number $p_{x}$, leading
to a macroscopic degeneracy. Furthermore, in the high-energy regime, a
holographic flatband system was also achieved [19].

####  1.1.2 Compact localized states

CLS is an eigenstate of a particular tight-binding lattice that is limited to
a finite system region and has no presence of amplitudes elsewhere. These
states exist in various systems. For instance, let us consider the one-
dimensional diamond chain (Fig. 1.2) without on-site potentials and with the
same hopping strength,

$\displaystyle\mathcal{H}=\sum_{n\in\mathbb{Z}}(\hat{b}_{n}^{\dagger}+\hat{b}_{n+1}^{\dagger})(\hat{a}_{n}+\hat{c}_{n})+\mathrm{h.c}.$
(1.6)

$\hat{a}_{n}$, $\hat{b}_{n}$, and $\hat{c}_{n}$ are annihilation operators at
site $a_{n}$, $b_{n}$, and $c_{n}$, respectively. Fourier transforms of
$\hat{a}^{\dagger}_{n}$ and $\hat{a}_{n}$ are defined as follows (same for
$\hat{b}^{\dagger}_{n},\hat{c}^{\dagger}_{n}$ and their conjugates transforms
similarly),

$\displaystyle\hat{a}^{\dagger}_{n}=\frac{1}{\sqrt{N}}\sum_{k\in\mathrm{BZ}}\exp(-ikn)\hat{A}^{\dagger}_{k}\quad\text{and}\quad\hat{a}_{n}=\frac{1}{\sqrt{N}}\sum_{k\in\mathrm{BZ}}\exp(ikn)\hat{A}_{k}$
(1.7)

$\mathrm{BZ}$ stands for the first Brillouin zone. The Hamiltonian in a
$k$-space, by means of Bloch’s theorem, is the following,

$\displaystyle\mathcal{H}=\sum_{k\in\mathrm{BZ}}\hat{\psi}_{k}^{\dagger}\mathcal{H}(k)\hat{\psi}_{k},\quad\mathcal{H}(k)=\begin{bmatrix}0&1+e^{ik}&0\\\
1+e^{-ik}&0&1+e^{ik}\\\ 0&1+e^{ik}&0\end{bmatrix}.$ (1.8)

We obtain three different energy bands: a flatband, $E_{\mathrm{FB}}=0$, and
two dispersive bands, $E_{\pm}(k)=\pm 2\sqrt{2}\cos(k/2)$. The energy bands
are drawn in Fig. 1.2. $E_{\mathrm{FB}}$ is macroscopically degenerated, and
we have the CLSs, which make a complete orthogonal set in the system both in
the $k$\- and real-space ($\ket{\varnothing}$ is a vacuum state),

$\displaystyle|\psi_{\mathrm{FB}}(k)\rangle=\frac{1}{\sqrt{2}}\left(\hat{a}^{\dagger}_{k}-c^{\dagger}_{k}\right)\ket{\varnothing}$
(1.9)
$\displaystyle|\psi_{\mathrm{FB}}(n)\rangle=\frac{1}{\sqrt{N}}\sum_{k\in\mathrm{BZ}}e^{ikn}|\psi_{\mathrm{FB}}(k)\rangle=\frac{1}{\sqrt{2}}\left(\hat{a}^{\dagger}_{n}-c^{\dagger}_{n}\right)\ket{\varnothing}$
(1.10)

In this setup, the compact localization relies on the system’s geometry, as
shown in Fig. 1.2, and each CLS is an anti-symmetric state occupying one unit
cell, with zero amplitude on all other sites. The CLSs are the eigenstates of
the Hamiltonian $\mathcal{H}$, and they remain unchanged over time, except for
a possible phase factor.

Figure 1.1:  A compact localized state of a one-dimensional diamond chain at
flatband energy $E_{\mathrm{FB}}$. Because of the lattice geometry, the wave
function is trapped only at $a_{n}$ and $c_{n}$ and no leakage to $b_{n}$ and
$b_{n+1}$. $E_{\mathrm{FB}}$ has a macroscopic degeneracy and CLS can be
placed at any unit cells. The linear combination of CLSs are again an
eigenstate corresponding to $E_{\mathrm{FB}}$.

Figure 1.2:  Energy bands of a one-dimensional diamond chain where $E$ and $q$
denoting the energy and crystal momentum, respectively. At
$E=E_{\mathrm{FB}}=0$, we see a flatband energy. Other two bands are the
dispersive bands.

This counter-intuitive behavior, where the excitation remains localized,
results from the _destructive interference_. The amplitudes for tunnelling or
leaking of the CLS to the neighboring sites cancel each other out. It is
crucial to note that compact localization differs fundamentally from disorder-
induced Anderson localization since CLSs can exist even in perfect
translation-invariant systems.

Compact localized states (CLSs) are intriguing for various reasons. For
instance, the temporal robustness of CLSs against perturbations is of
significant interest due to their exact confinement within their localization
region by destructive interference. Such characteristics find a potential
application such as an information storage in quantum networks [20].
Furthermore, CLSs have the potential to realize bound states within a
scattering continuum [21, 22]. Alongside these prospects, CLSs offer different
possible applications, such as in data transmission, where the interaction of
light modes from different fibers can lead to crosstalk issues [23, 24, 25].

####  1.1.3 Construction of flatband systems

In recent years, the existence of flatband systems has been experimentally
demonstrated in a variety of photonic and dissipative condensate networks [25,
26, 27]. From a theoretical perspective, a significant amount of research has
also been conducted in the recent decade, leading to a deeper understanding of
the mechanisms involved in generating flatband systems. Here, we summarize the
current principles of generating flatband Hamiltonians.

The pioneering work of a flatband generator is based on local network
properties and classifying flatband networks using orthogonal and non-
orthogonal CLS. The main method is to apply an inverse eigenvalue problem [28,
29, 30], studied extensively for one-dimensional networks with the nearest
unit cell interaction and a specific number of unit cells.

Flatband systems can also be classified into singular and non-singular
flatband systems [31, 32, 33, 34]. Singular flatbands exhibit immovable
discontinuities in their eigenstates, resulting from band touching with other
dispersive bands. Consequently, the compact localized states cannot form a
complete set spanning the singular flatband, and we require an additional
state. This so-called non-contractible loop state is compact localized line
that circles the entire system to have a complete set [31]. The non-
contractible loop states lead to a novel bulk-boundary correspondence, where
the presence of robust boundary modes is guaranteed by the singularity of the
Bloch wave function and the non-contractible loop states. Moreover, when the
degeneracy at the band crossing point is lifted, the singular flatband becomes
dispersive. It can acquire a finite Chern number, offering a pathway for
achieving nearly flat Chern bands. In contrast, non-singular flatbands display
no singularities in their Bloch wave functions. They can be fully isolated
from other bands while maintaining their perfect flatness. It is worth noting
that all one-dimensional flat bands belong to the non-singular class [31].
Exploiting these topological characteristics, a general scheme has been
developed to systematically construct flat-band model Hamiltonians, allowing
for the design of singular or non-singular nature [33].

Another way of constructing a flatband system is to use the concept of local
and latent symmetries of a lattice structure and their role in explaining
degeneracies in energy spectra [35, 36, 37]. Latent symmetries are local
symmetries of an effective Hamiltonian derived through subsystem partitioning
and isospectral reduction [38]. This work offers new perspectives for
analyzing accidental degeneracies in terms of latent symmetries, providing
insights into the symmetries underlying physical systems.

At the same time, a non-Hermitian system, which describes an open system with
gain and loss, has also gained attention in recent decades. Non-Hermitian
systems with dispersionless energy bands were first studied using symmetries
or specific geometries inspired by Hermitian models, both theoretically and
experimentally [39, 40, 41, 42, 43, 44, 45, 46, 47]. Later, the first
systematic approach to construct non-Hermitian flat bands using one-
dimensional two-band tight-binding networks was studied in Ref. [48]. This
method extends the techniques used for constructing Hermitian flat bands and
demonstrates that the non-Hermitian case allows for fine-tuned flatbands
unprotected by symmetry. Furthermore, it can potentially reveal a greater
variety of flatband types compared to the Hermitian case.

###  1.2 Experimental realizations

The precise arrangement of lattice geometry is crucial for fine-tuned flatband
lattices. As a result, finding actual materials featuring macroscopically
degenerate flatbands in nature is difficult. Such rarity presents a challenge
for many researchers in the experimental field. Advancements in fabrication
technologies have facilitated the creation and examination of artificial
lattices containing flatbands in various physical systems. This section
provides an overview of experimental methods employed in artificial flatband
systems within diverse frameworks. For details concerning the experimental
realizations of artificial flatband lattices, readers are encouraged to
consult the papers by Leykam et al. [26, 49].

_Flatbands in electronic systems_ : In 1998, Vidal et al. [50] discovered a
completely flat spectrum in periodic electronic networks induced by a specific
value of magnetic field, which is known as the Aharonov-Bohm caging [51]. This
phenomenon was proposed to be observable in superconducting wire networks
under appropriate conditions. Subsequently, indirect evidences of flatbands
due to Aharonov-Bohm caging was reported [52, 53]. Moreover, recent
advancements in two-dimensional material fabrication enabled the creation of
artificial flatband lattices [54, 55]. In 2017, 2D Lieb lattices were
constructed with two methods via scanning tunneling microscopy: creating
electron lattices by atom removal [56] and molecule addition to a substrate
[55, 57]. These scanning tunneling microscopy setups allows electron density
measurements, offering observations into flat and dispersive band Bloch waves.

_Flatbands in optical lattices_ : Progress in cold atom research, facilitated
by laser cooling and ion trapping, has enabled precise exploration of flatband
lattices, with the Lieb lattice standing out due to its simplicity in atom
manipulation. Conversely, lattices like the dice [58, 59] and kagomee [15,
60], though easier to construct, present challenges in generating highly
excited flatbands in cold atom systems. Early experiments sought to create
optical Lieb lattices [61], focusing on the Dirac cone in dispersive bands
rather than the flatband. These optical lattices unavoidably introduced a
width to the flatband due to next-nearest neighbor hopping. However,
optimizing laser beam intensities allowed for nearly flat bands with minimal
width [62]. In 2015, Takahashi et al. [63] achieved an optical Lieb lattice
for bosonic cold atoms by dynamically adjusting the lattice, leading to the
observation of interactions that caused condensate decay into lower dispersive
bands. Subsequent experiments in 2017 furthered the understanding of optical
Lieb lattices: the introduction of fermionic cold atoms [63] and the
exploration of interaction effects on flatband energy shifts [64].

_Flatbands in photonic systems_ : Flatbands are vital in slow light
applications in photonics [65], enhancing nonlinear effects with suppressed
group velocity and possibly useful for pulse buffering. However, achieving
ideal flatbands is challenging due to balancing low group velocity with a
useful bandwidth [66, 67]. Initial efforts with photonic crystal featured by
circular rods [68] faced limitations due to difficulties in fabrication.
Later, in 2017, innovative kagome structures [69] offered improved group
velocity reduction compared to the previous research [65, 70, 71],
revitalizing interest in flatbands in slow light engineering in photonic
crystal waveguides [72, 73, 74, 75, 76]. Moreover, the femtosecond laser
writing technique enabled the fabrication of optical waveguide networks,
providing control over coupling and facilitating the study of various lattice
types [77]. While early experiments with Lieb lattices indirectly inferred
flatbands [78] because a superposition of all bands was excited by a single
waveguide input, later studies directly excited compact localized states [79,
80].

_Flatbands in exciton-polariton condensate :_ Exciton-polaritons in
semiconductor microcavities, created through strong light-matter coupling,
offer a unique platform for exploring Bose-Einstein condensation and low-power
optical switching, facilitated by their low effective mass and exciton-
mediated nonlinearity. These features extend to the creation of structured
potentials, enabling the construction of periodic lattices and flatbands
within exciton-polariton systems [81]. Below the condensation threshold, the
linear band structure is easily measurable, as reported by Jacqmin et al. [82]
in 2014 when they successfully constructed a honeycomb lattice with a
flatband. Overcoming challenges related to precise control over condensation
states at higher pump powers, Baboux et al. [83] have introduced solutions
involving manipulating the spatial structure of the optical pump to induce
condensation into a flatband. In a 2D Lieb lattice with micropillars, exciton-
polariton condensation into a flatband was also realized [84, 85].

###  1.3 Effect of perturbation and interaction

In recent years, there has been a growing interest in physical systems that
exhibit macroscopic degeneracies. As previously discussed, one key
characteristic of these systems is the presence of flatbands. The intriguing
aspects of flatband models lie in their unique behavior by having a charge
carrier with zero group velocity, which leads to negligible kinetic energy and
extreme suppression of their motion. Then, the dominant energy scale would be
a perturbation and an interaction strength, eventually breaking the
macroscopic degeneracy no matter how weak they are. This intriguing behavior
acts as a catalyst for a wide variety of exotic and unconventional correlated
phases.

For instance, unconventional Anderson localizations [86, 87, 88], an inverse
Anderson transition [89, 90] where the re-entrant of localized phase exists in
a sequence of insulator-metal-insulator, and a non-perturbative metal-to-
insulator transitions [91, 92, 93] are among the distinct manifestations
observed due to various types of perturbations in several flatband systems. In
the realm of many-body flatband systems, the effect of different types of
interactions gives rise to several notable phenomena. These include a flatband
ferromagnetism [94, 95], a frustrated magnetism [96, 97], an ergodicity
breaking [98, 99, 100, 101, 102], a quantum caging of particles [103, 104,
105, 106], a formation of tightly bounded pairs in attractive interaction
[107], a pairing formation in repulsive interaction [103, 108, 109, 110] in
bosonic system, and flatband superconductivity [111, 112, 113, 114, 115, 116,
117, 118].

###  1.4 Contributed work and outline of the thesis

In Ref. [119, 120], I studied the effect of quasiperiodic perturbations on the
full manifold of one-dimensional two-band all-bands-flat lattice models. Such
networks can be diagonalized by a finite sequence of local unitary
transformations parameterized by angles $\theta_{i}$s. Without loss of
generality, the case of two bands with bandgap $\Delta$ is focused. The
primary focus is on a two-band ABF ladder subjected to a quasiperiodic
perturbation that is comparatively weak relative to the bandgap. Weak
perturbations lead to an effective Hamiltonian with both on- and off-diagonal
quasiperiodic terms that depend on $\theta_{i}$s. For some angle values, the
effective model coincides with the extended and the off-diagonal Harper model.
By varying the parameters of the quasiperiodic potentials, localized
insulating states and an entire parameter range hosting critical states over
the full spectrum with subdiffusive transport are observed. Furthermore, we
identify and refer to the transition between these states as the _critical-to-
insulator transition._ For finite quasiperiodic potential strength, the
critical-to-insulator transition becomes energy-dependent with what we term
_fractality edges_ separating localized from critical states. Chapter 2
provides an overview of a systematic investigation into the effect of
quasiperiodic perturbations on one-dimensional ABF networks.

I also studied the effect of on-site interaction of hard-core bosons in one-
dimensional and two-dimensional cross-stitch lattices with flatbands. One of
the key findings is that the groundstate energy is only influenced by the
CLSs. The Wigner crystal arises when CLSs are completely occupied, constrained
by the mutual repulsion from hard-core boson constraint. On the other hand,
the band-insulating phase occurs when all lattice sites are fully filled,
leading to an insulating state due to the absence of available states for
conduction. In cases where the lattice is partially filled with the CLSs and
with the fully-filled dimers, a mixture of both phases is observed. Moreover,
when a closed barrier of compact localized states is present in both the one-
and two-dimensional cross-stitch lattices, hard-core bosons inside the loop
become physically trapped in the loop. This violation of weak thermalization
gives rise to the presence of non-ergodic excited states, causing Hilbert
space fragmentation due to the emergence of artificial localization arising
from CLS barriers. Chapter 3 provides an overview of the thermalization
properties of hard-core bosons in the one- and two-dimensional cross-stitch
lattices with on-site repulsive interaction.

Ref. [121] focuses on creating CLSs in an electrical diamond lattice made
solely of capacitors and inductors. An excitation of CLS mode is achieved
through local linear driving near the flatband frequency of the lattice. This
research marks a step forward in establishing a versatile circuit platform for
creating and manipulating flatband systems. The comparison between
experimental results and numerical simulations reveals a good agreement.
Moreover, the investigation includes the examination of local nonlinear
driving by substituting capacitors with varactors. The introduction of lattice
nonlinearity facilitates the generation of a nonlinear CLS continuation. At
the same time, the study considers a one-dimensional stub lattice, which falls
in a different flatband class. An important observation is that local driving
cannot isolate a single CLS due to its non-orthogonality, unlike the diamond
lattice case. My main contribution to this project include a numerical
modeling of an appropriate nonlinear capacitance of a varactor to achieve good
agreement with the experiment and the general manuscript work. Additionally, I
investigate the one-dimensional stub chain, demonstrating the inability to
isolate CLS with local driving. Chapter 4 presents an overview of generating
one-dimensional flatband electrical lattices through capacitors and inductors.

## Chapter 2 Generating Critical States
from Quasiperturbed Flatbands

###  2.1 Introduction

In condensed matter physics, an area of significant interest and exploration
in recent years has revolved around comprehending the effects of various
perturbations on single-particle localized states. While it is widely
recognized that states can become exponentially localized through the
influence of random disorder or quasiperiodic potential [1, 122, 123],
localization can also be attained even in the absence of disorder/perturbation
as it is mentioned in the introduction, which is a compact localization in a
flatband system. The fascination with flatbands stems from their extreme
sensitivity to perturbations that break their macroscopic degeneracy.
Different types of perturbation lead to an evolution of CLSs, making a system
that generates a diverse range of intriguing and unconventional phases.

Further fine-tuning of flatband systems can result in _all-bands-flat_ (ABF)
networks [50, 105], where all dispersive bands are made flat, making them even
highly sensitive to perturbations. A finite sequence of non-commuting local
unitary transformations can transform an ABF network into a diagonal
Hamiltonian, a _parent network_ of the corresponding ABF network. Any
perturbation applied to the original ABF network leads to a nontrivial
perturbation in the parent network, forming non-zero non-diagonal elements.
Recently, the effect of random disorder on ABF networks has been
systematically explored [91]. Nonperturbative metal-to-insulator transitions
were discovered in three-dimensional ABF systems, while in the one-dimensional
case, only Anderson localization occurred across the entire spectrum [91]. We
can also question about what happens if we consider quasiperiodic perturbation
on the one-dimensional ABF networks.

In this chapter, we summarize a systematic study of the impact of
quasiperiodic perturbations on the one-dimensional ABF networks, published in
Ref. [119, 120]. Without loss of generality, we focus on a two-band ABF ladder
and introduce a quasiperiodic perturbation, which is weak compared to the
bandgap. We use the smallness of the perturbation to project the Hamiltonian
onto a single sublattice, thereby deriving a new effective projected
Hamiltonian. By varying the parameters of the perturbation, we find that the
entire spectrum of the projected Hamiltonian is either localized or
critical—the _critical-to-insulator transition_ (CIT), but never extended,
using a mapping to the extended Harper model [124, 125]. Increasing the
strength of the quasiperiodic potential, the critical eigenstates are
gradually replaced by localized ones via the appearance of an energy-dependent
CIT [126, 127, 128, 129, 130] which is called as _anomalous mobility edges_
(but, in this chapter, those edges are dubbed as _fractality edges_). The
emergence of regions of critical states in one-dimensional systems is quite
different compared to the conventional Aubry-André model, which features a
metal-to-insulator transition.

The chapter is organized as follows. We start by defining and discussing the
construction of the all-bands-flat models in Sec. 2.2 and the numerical
methods in Sec. 2.3. In Sec. 2.4, we derive an effective model valid in the
limit of weak quasiperiodic perturbation and use it to chart the phase
diagram, confirmed numerically. The transport properties at the weak
perturbation is discussed in Sec. 2.5. The properties of the full model at
finite perturbation strength are investigated in Sec. 2.6, followed by
conclusions in Sec. 2.7.

###  2.2 Model

We focus on the one-dimensional ABF ladder, which consists of two flatbands
with the nearest-neighbor unit cells hopping [29]. An illustration of the
model and the scheme is provided in Fig 2.1.

####  2.2.1 Local unitary transformation and ABF lattice

ABF Hamiltonian $\mathcal{H}_{\mathrm{ABF}}$ can be constructed from a
macroscopically degenerate diagonal matrix $\mathcal{H}_{\mathrm{FD}}$ with
onsite energies $\varepsilon_{a}$ and $\varepsilon_{b}$ on the two sublattices
having the bandgap $\Delta=\absolutevalue{\varepsilon_{a}-\varepsilon_{b}}$
[105],

$\displaystyle\mathcal{H}_{\mathrm{FD}}=\sum_{n}\varepsilon_{a}\outerproduct{a_{n}}{a_{n}}+\varepsilon_{b}\outerproduct{b_{n}}{b_{n}}.$
(2.1)

The state vector $\ket{\Psi}$ of $\mathcal{H}_{\mathrm{FD}}$ is written as
linear combinations of orthonormal unit basis
$\\{\ket{a_{n}},\ket{b_{n}}:n\in\mathbb{Z}\\}$ where $a_{n},b_{n}$ are the
amplitudes at the two sublattices,

$\displaystyle\ket{\Psi}=\sum_{n\in\mathbb{Z}}a_{n}\ket{a_{n}}+b_{n}\ket{b_{n}},\quad\text{where
}\sum_{n}|a_{n}|^{2}+|b_{n}|^{2}=1.$ (2.2)

We term $\mathcal{H}_{\mathrm{FD}}$ as the parent Hamiltonian for the manifold
of ABF systems and refer to it as _fully detangled_. [88, 105, 99].

Figure 2.1:  The ABF model is fully detangled/diagonalized into a diagonal
Hamiltonian via local unitary transformations $U_{1}$ and $U_{2}$ with
appropriate angles $\theta_{1}$ and $\theta_{2}$ (left). By adding
quasiperiodic perturbation $W$ made up of two quasiperiodic fields $W_{1,2}$
in Eq. (2.5), nontrivial hoppings are created in the fully detangled basis
(right). For weak quasiperiodic fields, the first-order degenerate
perturbation theory is used to derive an effective projected model. In the
detangled basis, this corresponds to only keeping the hopping terms coupling
the sites on the same sublattices.

The construction of ABF systems is based on a sequence of unitary
transformations $U_{i}$ [88, 105, 99, 91, 119, 120] applied to the fully
detangled Hamiltonian. It is a way of achieving a set of flatband lattices
that is distinct from the method mentioned in Sec. 1.1 in Chapter 1. Each
transformation $U_{1,2}$ is a direct sum of local transformations, i.e., it
takes a block diagonal form. For a one-dimensional system with the nearest-
neighbor unit cell hopping, only two local unitary transformations
$U=U_{2}U_{1}$ are enough to produce a connected hopping network. For a one-
dimensional tight-binding system with two sites per unit cell, the most
general $2\times 2$ block is a $\mathrm{SU}(2)$ matrix,

$\displaystyle
U_{i}\\!=\\!\\!\\!\\!\\!\sum_{n,n^{\prime}\in\mathbb{N}}\\!\\!\\!\\!z_{i}|{a^{(i)}_{n}}\rangle\langle{a^{(i-1)}_{n}}|\\!+\\!w_{i}|{a^{(i)}_{n}}\rangle\langle{b^{(i-1)}_{{n}^{\prime}}}|\\!-\\!w_{i}^{*}|{b^{(i)}_{n}}\rangle\langle{a^{(i-1)}_{{n}^{\prime}}}|\\!+\\!z_{i}^{*}|{b^{(i)}_{n}}\rangle\langle{b^{(i-1)}_{n}}|.$
(2.3)

The index $i$ denotes the $i$-th local unitary transformation, and the indices
$n,n^{\prime}$ label the unit cells. In the simplest case, the blocks are
parameterized by only two angles: $\theta_{1,2}$ for $U_{1,2}$, respectively,
producing real $U_{1,2}\in\mathrm{SO(2)}$. The $2\times 2$ blocks can act on
the $2$ sites within the same unit cell, $n^{\prime}=n$, or different unit
cells, $n^{\prime}\neq n$. This generates different, non-commuting
transformations. Then, we obtain additional tunneling between sites, giving a
connected ABF network. Throughout the chapter, ABF Hamiltonian always refers
to the connected ABF network. The energy spectra is preserved due to the
property of a unitary transformation, and the local rotation allows us to
maintain the localized/delocalized property of a given eigenstate.
Specifically, we choose the $U_{1}\in\mathrm{SO(2)}$ blocks to act within the
same unit cell, while for the $U_{2}\in\mathrm{SO(2)}$ blocks, one of the
sublattice sites is taken from the neighboring unit cell, $n^{\prime}=n-1$.
Next, we add a quasiperiodic perturbation to the ABF Hamiltonian with
$\mathcal{H}=\mathcal{H}_{\mathrm{ABF}}+W$. Here, $W$ is defined as a direct
sum of $2\times 2$ matrices $W(n)$ of all $1\leq n\leq L$, where $L$ is the
number of unit cells,

$\displaystyle W=\bigoplus_{n=1}^{L}W(n),\qquad
W(n)=\begin{bmatrix}W_{1}(n)&0\\\ 0&W_{2}(n)\end{bmatrix}.$ (2.4)

The generic form of $W(n)$ is given by two quasiperiodic fields $W_{1}$ and
$W_{2}$,

$\displaystyle W_{1}(n)=\lambda_{1}\cos(2\pi\alpha n+\phi)\quad\text{and}\quad
W_{2}(n)=\lambda_{2}\cos(2\pi\alpha n+\beta+\phi).$ (2.5)

Here, the spatial frequency $\alpha$ is an irrational number,
$\alpha\in\mathbb{R}\setminus\mathbb{Q}$, $\beta$ is the phase difference
between $W_{2}$ and $W_{1}$, $\phi$ is the phase used for averaging and
$\lambda_{1}$ and $\lambda_{2}$ are the strengths of the quasiperiodic
potentials. Taking two angles $\theta_{1,2}$ into account, we can tune five
total parameters to understand how this perturbation affects different ABF
networks.

####  2.2.2 Geometric symmetry of ABF lattice

The ABF Hamiltonian
$\mathcal{H}_{\mathrm{ABF}}=U_{2}U_{1}\mathcal{H}_{\mathrm{FD}}(U_{2}U_{1})^{\dagger}$
has a rotational symmetry and it allows us to obtain the irreducible domain of
angles $(\theta_{1},\theta_{2})$. Suppose $\theta_{1}\to\theta_{1}+\pi$ and
$\theta_{2}$ remains the same. Then $U_{1}$ becomes $-U_{1}$. The ABF network
with shifted unitary transformation
$(-U_{2}U_{1})\mathcal{H}_{\mathrm{FD}}(-U_{2}U_{1})^{\dagger}$ is equal to
$\mathcal{H}_{\mathrm{ABF}}$. The same situation occurs for $\theta_{2}$
shifting to $\theta_{2}+\pi$. If both $\theta_{1}$ and $\theta_{2}$ are
shifted, then we retrieve the original ABF system as well.

The second case is when $\theta_{1}$ is reversed and shifted by $\pi$, hence
$\pi\\!-\\!\theta_{1}$. Initially,
$I\\!=\\!U^{\dagger}_{1}(\theta_{1})U_{1}(\theta_{1})$ is applied on the right
and the left end of
$U_{1}(\pi\\!-\\!\theta_{1})\mathcal{H}_{\mathrm{FD}}U^{\dagger}_{1}(\pi\\!-\\!\theta_{1})$.
By using the relation $U_{1}(\pi\\!-\\!\theta_{1})U_{1}(\theta_{1})\\!=\\!-I$,
it is easy to verify that
$U_{1}(\pi\\!-\\!\theta_{1})\mathcal{H}_{\mathrm{FD}}U^{\dagger}_{1}(\pi\\!-\\!\theta_{1})$
and $U_{1}(\theta_{1})\mathcal{H}_{\mathrm{FD}}U^{\dagger}_{1}(\theta_{1})$
are equivalent. Hence, the overall ABF system will not change. The similar
result is obtained for $\theta_{2}$ being shifted to $\pi\\!-\\!\theta_{2}$ as
well. We can further think for the case when $\theta_{1}\to\pi/2-\theta_{1}$
and $\theta_{2}$ remains at it is. However, such a case will not lead to a
symmetry. As a conclusion, it is enough to consider
$\theta_{1},\theta_{2}\in[0,\pi/2]$.

###  2.3 Numerical methods

The relation between localization and statistical properties of the energy
spectrum and eigenstates in disordered systems is widely recognized.
Eigenstates that exhibit broad spatial distribution suggest an electronic
transport, with the localization length being inversely related to the
exponent describing the decay of eigenstate components across system sites. An
exact formula for the localization length in the one-dimensional lattice model
remains unknown, although field-theoretic methods have provided expressions
for extreme cases of strong and weak random disorder [131]. Therefore,
numerical approaches are often employed to approximate the localization
length.

####  2.3.1 Inverse participation ratio

The generalized inverse participation ratio (GIPR) is a commonly used
numerical quantity. It has been utilized to measure the
localization/delocalization of the normalized eigenstate $\psi(E)$ of energy
$E$ in the analysis of a tight-binding system,

$\displaystyle\operatorname{GIPR}_{q}(E)=\sum_{n\in\mathbb{z}}|\psi_{n}(E)|^{2q}\sim
L^{-\tau_{q}}.$ (2.6)

Especially when $q=2$, the physical interpretation is well-established, and
the quantity denoted as $\operatorname{IPR}(E)$ is commonly known as the
inverse participation ratio (IPR), serving as a measure of spatial extension
for a given state. To illustrate, consider two extreme examples. In the case
of a homogeneous state $\psi(E)$, we obtain $|\psi(E)|^{2}\sim 1/L$ where $L$
is the size of a lattice. This leads to $\operatorname{IPR}(E)\sim 1/L$. On
the other hand, if $\psi(E)$ is localized at a single site, then IPR is equal
to 1,

$\displaystyle\operatorname{IPR}(E)=\begin{dcases}1/L\quad&\tau_{2}=1,\\\
1\quad&\tau_{2}=0.\end{dcases}$ (2.7)

####  2.3.2 Defining a linear model for localized states

The quantity $\tau=\tau_{2}$ is a scaling exponent of the $\mathrm{IPR}$. In
the thermodynamic limit, $\tau$ is defined explicitly as follows at energy $E$
[132],

$\displaystyle\tau=\lim_{L\to\infty}\frac{1}{\ln(1/L)}\ln\operatorname{IPR}(E).$
(2.8)

However, we cannot use this definition for a single eigenstate as a function
of system size $L$. Changing the system sizes changes the eigenenergies and
does not allow their smooth continuation from one system size to another.
Instead, we calculate the average of $\tau$ over a small energy bin to extract
the scaling behavior of $\tau$. This is achieved by rescaling the
eigenspectrum into the interval $[0,1]$ and splitting it into equidistant bins
$e$. This allows us to calculate the average of $\tau$ in a single bin $e$,
denoted as $\langle\tau(e)\rangle$ as a function of lattice size $L$. It
follows:

$\displaystyle\lim_{L\to\infty}\langle\tau(e)\rangle=\begin{dcases}0\text{,
all states localized in $e$,}\\\ \tau_{0}\text{, some states not localized in
$e$}.\end{dcases}$ (2.9)

In the thermodynamic limit, all eigenstates in the bin $e$ are localized if
$\langle\tau(e)\rangle(L\to\infty)=0$. Otherwise,
$\langle\tau(e)\rangle(L\to\infty)$ takes a finite value $\tau_{0}$ between
$0$ and $1$. If we take the entire spectrum to get the average $\tau$, the
dependence of $e$ is dropped and it is written as $\langle\tau\rangle$. To
confirm the localization of an eigenstate, we use the fact that $\tau=0$ in
the thermodynamic limit and introduce the following linear model
$\langle\hat{\tau}(e)\rangle(L)$ to be optimized with given data of
$\langle\tau(e)\rangle$, inspired by Eq. (2.8),

$\displaystyle\langle\hat{\tau}(e)\rangle(L)=\frac{\hat{\kappa}(e)}{\ln(1/L)}+\langle\tau(e)\rangle(L\to\infty),\quad\langle\tau(e)\rangle(L\to\infty)=0.$
(2.10)

The slope $\hat{\kappa}(e)$ defines the participation number (an inverse of
inverse participation ratio, $\mathrm{IPR}^{-1}(e)$) averaged over $e$, and it
is obtained by taking $\exp(-\hat{\kappa}(e))$.

An important quantity in such a statistical approach is the $R^{2}$ obtained
from the goodness-of-fit test [133]. $R^{2}$ reflects how well the data is
described with a given statistical model; $R^{2}\approx 1$ implies that the
linear model is well optimized, and the states in bin $e$ are localized in the
thermodynamic limit. On the other hand, $R^{2}\leq 0$ indicates that the set
of numerical data does not follow the linear model in Eq. (2.10) at all and
cannot be appropriately optimized. Hence, a different fitting model should be
used for extended or critical states. The definition of $R^{2}$ at an energy
bin $e$ is given as

$\displaystyle
R^{2}(e)=1-\sum_{L}\left(\langle\tau(e)\rangle(L)-\langle\hat{\tau}(e)\rangle(L)\vphantom{\overline{\langle\tau(e)\rangle}}\right)^{2}\Big{/}\sum_{L}\left(\langle\tau(e)\rangle(L)-\overline{\langle\tau(e)\rangle}\right)^{2},$
(2.11)

where $L$ is an index of different lattice sizes.
$\langle\hat{\tau}(e)\rangle(L)$ represents an optimal value of
$\langle\tau(e)\rangle(L)$ obtained from numerically optimized
$\hat{\kappa}(e)$, and the numerator of the fraction corresponds to the
residual sum of squares. On the other hand, $\overline{\langle\tau(e)\rangle}$
is an average of $\langle\tau(e)\rangle(L)$ across all lattice sizes. The
denominator of the fraction is the total sum of difference squares.

The calculation of $R^{2}$ is essential and straightforward when identifying
localized states because we expect the IPR scaling exponent $\tau$ to be zero
in the thermodynamic limit. Alternatively, the more common approach can be
used to directly extract the coefficients $\tau$ themselves. This can be
achieved, for instance, using the following scaling ansatz,

$\displaystyle\langle\mathrm{IPR}\rangle\sim L^{-\tau_{m}}.$ (2.12)

In this way, the exponents for the mean IPR, $\tau_{m}$, can be extracted.
Localized states are identical by $\tau_{m}=0$, delocalized states by
$\tau_{m}=1$, and critical states by $0<\tau_{m}<1$. This approach also allows
us to distinguish between delocalized and critical states. However, the
$R^{2}$ based approach is more convenient for identifying the localized states
because we know that the IPR scaling exponent (the intercept) $\tau$ must be
zero in the thermodynamic limit. However, diagnosing states other than the
localized states is challenging or inapplicable unless one knows the exact
value of $\tau$ in the thermodynamic limit of some states.

###  2.4 Weak perturbation

We start by examining the limit of weak quasiperiodic perturbation, i.e.,
vanishing $\lambda_{1,2}$, as compared to the bandgap $\Delta$ by applying the
first-order degenerate perturbation theory. We apply the inverse unitary
transformation $U^{\dagger}$ to $\mathcal{H}$ so that the Hamiltonian in a
fully detangled basis is
$\tilde{\mathcal{H}}=\mathcal{H}_{\mathrm{FD}}+U^{\dagger}WU$. The second term
represents hoppings solely due to the quasiperiodic perturbation. The
strongest enhancement of these hoppings occurs for $\theta_{1,2}=\pi/4$.
Without loss of generality, we set $\lambda_{2}\leq\lambda_{1}$. Factorizing
out $\lambda_{1}$ from $\tilde{\mathcal{H}}$ gives us the perturbation
$\tilde{W}=W/\lambda_{1}$, which depends on the ratio
$\lambda_{2}/\lambda_{1}\in[0,1]$ only. Then, we get the following eigenvalue
problem via first-order degenerate perturbation theory,

$\displaystyle
P_{a}U\tilde{W}U^{\dagger}P_{a}\ket{a_{n}}=\lambda_{1}\varepsilon^{(1)}_{a,n}\ket{a_{n}}.$
(2.13)

Here $\ket{a_{n}}$ are the states in the fully detangled basis, and $P_{a}$ is
a projection operator onto flatband $\varepsilon_{a}$ that is local. The above
equation describes an effective 1D tight binding problem called the _projected
model_ [91, 119, 120]. Fig. 2.1 presents the schematics of obtaining the
projected model: without perturbation, reverting the unitary transformation
$U$ gives two sublattices of decoupled sites, while one gets onsite energies
and hoppings from $U^{\dagger}WU$ due to the perturbation. The projected
model, valid for weak interactions, neglects the couplings between the
sublattices, producing two decoupled projected models. The flatband energy we
want to focus on dictates the choice of the sublattice/projected model. On
general grounds, we expect the effective model to feature both quasiperiodic
onsite energies and finite-range hopping, given that $\tilde{W}$ is local and
quasiperiodic, and $U^{\dagger}P_{a}$ is a local operator. In the fully
detangled basis, the effective problem for a particular sublattice (such as
the lower leg in Fig. 2.1) is as follows,

$\displaystyle Ea_{n}=v_{n}a_{n}+t_{n-1}a_{n-1}+t_{n}^{*}a_{n+1},$ (2.14)

where the onsite potential $v_{n}$ and the hopping $t_{n}$ are indeed
quasiperiodic,

$\displaystyle v_{n}$ $\displaystyle=v_{s}\sin(2\pi\alpha
n-\pi\alpha)+v_{c}\cos(2\pi\alpha n-\pi\alpha),$ (2.15) $\displaystyle t_{n}$
$\displaystyle=t_{s}\sin(2\pi\alpha n)+t_{c}\cos(2\pi\alpha n).$ (2.16)

The expressions for the amplitudes $v_{s,c}$ and $t_{s,c}$ of $v_{n}$ and
$t_{n}$ are quite intricate, as they contain details from both quasiperiodic
fields and the local unitary transformation. Those amplitudes are obtained as
follows (no derivation shown),

$\displaystyle v_{s}$
$\displaystyle=\sin(\pi\alpha)\left(\sin^{2}\theta_{1}\sin^{2}\theta_{2}-\cos^{2}\theta_{1}\cos^{2}\theta_{2}\right)\vphantom{\frac{\lambda_{2}}{\lambda_{1}}}$
(2.17)
$\displaystyle+\frac{\lambda_{2}\cos\beta}{\lambda_{1}}\sin(\pi\alpha)\left(\sin^{2}\theta_{1}\cos^{2}\theta_{2}-\cos^{2}\theta_{1}\sin^{2}\theta_{2}\right)$
$\displaystyle-\frac{\lambda_{2}\sin\beta}{\lambda_{1}}\cos(\pi\alpha)\left(\cos^{2}\theta_{1}\sin^{2}\theta_{2}+\sin^{2}\theta_{1}\cos^{2}\theta_{2}\right),$
$\displaystyle v_{c}$
$\displaystyle=\cos(\pi\alpha)\left(\cos^{2}\theta_{1}\cos^{2}\theta_{2}+\sin^{2}\theta_{1}\sin^{2}\theta_{2}\right)\vphantom{\frac{\lambda_{2}}{\lambda_{1}}}$
(2.18)
$\displaystyle+\frac{\lambda_{2}\sin\beta}{\lambda_{1}}\sin(\pi\alpha)\left(\sin^{2}\theta_{1}\cos^{2}\theta_{2}-\cos^{2}\theta_{1}\sin^{2}\theta_{2}\right)$
$\displaystyle+\frac{\lambda_{2}\cos\beta}{\lambda_{1}}\cos(\pi\alpha)\left(\cos^{2}\theta_{1}\sin^{2}\theta_{2}+\sin^{2}\theta_{1}\cos^{2}\theta_{2}\right),$
$\displaystyle t_{s}$ $\displaystyle=\frac{1}{4}\sin 2\theta_{1}\sin
2\theta_{2}\frac{\lambda_{2}\sin\beta}{\lambda_{1}}$ (2.19) $\displaystyle
t_{c}$ $\displaystyle=\frac{1}{4}\sin 2\theta_{1}\sin
2\theta_{2}\left(1-\frac{\lambda_{2}\cos\beta}{\lambda_{1}}\right).$ (2.20)

####  2.4.1 Extended Harper model

The projected model in Eq. (2.14) features both a quasiperiodic onsite
potential and hopping. The most generic model of this type is defined by a
self-adjoint quasiperiodic Jacobi operator $J$ acting on a vector of
$l^{2}(\mathbb{Z})$ space [134],

$\displaystyle(Ju)_{n}=v(\alpha n+\phi)u_{n}+c_{\rho}(\alpha
n+\phi)u_{n+1}+c_{\rho}^{*}(\alpha(n-1)+\phi)u_{n-1},$ (2.21)

where $\alpha$ is an irrational spatial frequency and $\phi$ is a fixed phase,
and $\rho$ is some additional parameter controlling the hopping $t$. No
complete phase diagram or transport properties — numerical or analytical —
have been established for this most generic case. However, the projected model
can be mapped onto the already studied extended Harper model for specific
parameter values, as discussed below. The definition of the extended Harper
model is given as follows,

$\displaystyle(Ju)_{n}$ $\displaystyle=2\cos(2\pi\alpha n-\pi\alpha)u_{n}$
(2.22) $\displaystyle+2\rho\left[\cos(2\pi\alpha
n)u_{n+1}+\cos(2\pi\alpha(n-1))u_{n-1}\right].$

Its spectral properties have been characterized completely in Ref. [124]. For
$2\rho<1$, the eigenstates are all localized, as guaranteed by the RAGE
theorem [3, 4, 5] (see the introduction of Chapter 1). Otherwise, for
$2\rho\geq 1$, the energy spectrum has a fractal structure [124], similar to
the Aubry-André-Harper model at its critical point. In this chapter, such
transition is called as a critical-to-insulator transition analogous to a
metal-to-insulator transition. Moreover, the corresponding eigenstates show
multifractal behavior, as demonstrated numerically in [135].

####  2.4.2 Off-diagonal Harper model

The onsite potential can be neglected for effectively infinite $2\rho$, and we
are left with the hopping terms only. The model in this limit is called the
off-diagonal Harper model (OHM) [136, 137],

$\displaystyle(Ju)_{n}=2\rho\left[\cos(2\pi\alpha
n)u_{n+1}+\cos(2\pi\alpha(n-1))u_{n-1}\right].$ (2.23)

It is self-dual under a modified Fourier transformation similar to the duality
of the Aubry-André model. [119].

$\ket{k}=\frac{1}{\sqrt{L}}\sum_{n}e^{-i2\pi\alpha nk}\ket{a_{n}}.$ (2.24)

Substituting the above transformation to the projected model, we obtain the
Hamiltonian has a similar structure in real-space and Fourier-space.
Consequently, the entire spectrum of the model is critical, and the
eigenstates are multifractal,

$\displaystyle H$ $\displaystyle=2\rho\sum_{n}\cos(2\pi\alpha
n)\left(\vphantom{\sum}{\ket{a_{n}}}\bra{a_{n+1}}+\ket{a_{n+1}}\bra{a_{n}}\right)$
$\displaystyle=2\rho\sum_{k}\cos(\pi\alpha(2k+1))\left(\vphantom{\sum}e^{-i\pi\alpha}\ket{k}\bra{k+1}+e^{i\pi\alpha}\ket{k+1}\bra{k}\right)$
(2.25)

####  2.4.3 Lyapunov exponent and localization length

After first realizing that the distribution of eigenvalues and the electron
localization in a disordered system are related in Ref. [138], Thouless came
up with a simple argument that can be applied to any random one-dimensional
tight-binding system with next nearest neighbor hopping [139],

$\displaystyle
E_{\beta}u_{n}^{\beta}=v_{n}u_{n}^{\beta}-t(u_{n+1}^{\beta}+u_{n-1}^{\beta}).$
(2.26)

His approach was to take the retarded Green function $G^{+}(E)$ of the system
and obtain the matrix element $G^{+}_{1N}(E)$, which gives the transition
amplitude of an electron from the first to the last site using the fact that
the corresponding cofactor (upper triangular matrix) is the product of
diagonal elements,

$\displaystyle
G^{+}(E)=\left(E-H+i\varepsilon\right)^{-1}=\frac{\operatorname{adj}(E-H+i\varepsilon)}{\det(E-H+i\varepsilon)},$
(2.27) $\displaystyle G_{1N}^{+}(E)=\langle
1|G^{+}(E)|N\rangle=t^{N-1}\left[\prod_{\alpha=1}^{N}(E-E_{\alpha}+i\varepsilon)\right]^{-1}.$
(2.28)

Now, let us focus on a specific eigenvalue $E_{\beta}$. Then, the explicit
calculation of $|G_{1N}^{+}(E_{\beta})|$ gives the following equation,

$\displaystyle|G_{1N}^{+}(E_{\beta})|=|u_{1}^{\beta}u_{N}^{\beta}|=|t^{N-1}|\left[\prod_{\alpha\neq\beta}^{N}\big{|}E_{\beta}-E_{\alpha}+i\varepsilon)\big{|}\right]^{-1}.$
(2.29)

If $u_{1}^{\beta}$ and $u_{N}^{\beta}$ are exponentially decaying with its
maximum value at the site $m$, the product of $u_{1}^{\beta}$ and
$u_{N}^{\beta}$ is the following,

$\displaystyle|u_{1}^{\beta}u_{N}^{\beta}|=\exp(-\lambda_{\beta}(m-1)-\lambda_{\beta}(N-m))=\exp(-\lambda_{\beta}(N-1)).$
(2.30)

The exponent $\lambda_{\beta}$ is obtained as follows, and it is called
_Thouless formula_ ,

$\displaystyle\lambda_{\beta}=\\!\\!\lim_{N\to\infty}\frac{1}{N-1}\sum_{\alpha\neq\beta}\ln|E_{\beta}-E_{\alpha}|-\ln
t=\\!\\!\int_{\mathbb{R}}\\!\\!\ln|E_{\beta}-E^{\prime}|\rho(E^{\prime})dE^{\prime}-\ln
t.$ (2.31)

In one-dimensional case, especially, an inverse of $\lambda_{\beta}$ defines
the localization length $\xi(E)$ at the energy $E_{\beta}$. For the continuum
case, the summation becomes the integral described with a (integrated) density
of state $\rho(E)$. Later on, several mathematicians rigorously argued the
Thouless formula that $\lambda_{\beta}$ turned out to be the Lyapunov exponent
of the associated transfer matrix.

The preceding analysis of a Lyapunov exponent was focused on a random system.
However, the Thouless formula may not be applicable for a quasiperiodic
system. To address this challenge, Avila came up with a theory that extends
the formula for a Lyapunov exponent to one-dimensional quasiperiodic systems
with nearest neighbor hoppings [140], so-called the _global theory_. First, we
write a transfer matrix for a quasiperiodic Jacobi operator (E. (2.21)) in the
following form,

$\displaystyle
T(E;\varphi)=\frac{1}{c_{\rho}(\varphi)}B(E;\varphi)\quad\text{and}\quad
B(E;\varphi)=\begin{bmatrix}E-v(\varphi)&-c_{\rho}^{*}(\varphi-\alpha)\\\
c_{\rho}(\varphi)&0\end{bmatrix}.$ (2.32)

The angle $\varphi$ is $\alpha n+\phi$. $B(E;\varphi)$ is a _non-singular_ ,
ensuring that the matrix is always well-defined and avoids issues with
division by zero for each element [141]. Then, the product of the transfer
matrices at a given energy $E$ gives us the update of the corresponding wave
function elements,

$\displaystyle T^{(n)}(E;\varphi)=T(E;\varphi+(n-1)\alpha)\cdots
T(E;\varphi),$ (2.33)

and the averaged asymptotic of the quasiperiodic cocycle can also be
quantified with Lyapunov exponent by definition,

$\displaystyle\lambda(E)$
$\displaystyle=\lim_{n\to\infty}\frac{1}{n}\int_{\mathbb{T}}\ln\lVert
T^{(n)}(E;\varphi)\rVert d\varphi$ (2.34)
$\displaystyle=\gamma(E)-\lim_{n\to\infty}\frac{1}{2\pi
n}\int_{0}^{2\pi}\\!\\!\\!\ln|c(\varphi)|d\varphi$ (2.35)
$\displaystyle\gamma(E)$ $\displaystyle=\lim_{n\to\infty}\frac{1}{2\pi
n}\int_{0}^{2\pi}\\!\\!\\!\ln\lVert B^{(n)}(E;\varphi)\rVert d\varphi$ (2.36)

The matrix norm is defined as the largest Lyapunov exponent. The main strategy
of the global theory is to complexify the Lyapunov exponent by making
$\varphi\to\varphi+i\varepsilon$. By then considering the asymptotic limit
$|\varepsilon|\to 0$, we can effectively compute $\lambda(E)$ and perform
extrapolation to $\varepsilon=0$. This method enables us to evaluate the
Lyapunov exponent accurately for the one-dimensional quasiperiodic system,

$\displaystyle\lambda_{\varepsilon}(E)=\lim_{n\to\infty}\frac{1}{2\pi
n}\int_{0}^{2\pi}\\!\\!\\!\ln\lVert T^{(n)}(\varphi+i\varepsilon)\rVert
d\varphi.$ (2.37)

The integral mentioned above is split into two distinct terms. The first term
represents the Lyapunov exponent of the complexified $B(E;\theta)$, while the
second term is the contribution of quasiperiodic hopping, as evident from Eq.
(2.32),

$\displaystyle\lambda_{\varepsilon}(E)=\lambda(\alpha,B(E;\varphi+i\varepsilon))-\int_{\mathbb{T}}\\!\ln|c(\varphi)|d\varphi.$
(2.38)

For the extended Harper model defined in Eq. (2.22), the Lyapunov exponent of
the complexified $B(E;\theta)$ is given in Ref. [141] and the contribution of
quasiperiodic hopping is explicitly obtained. When $\rho$ is equal to or
larger than $1/2$, the Lyapunov is zero; otherwise, we have a positive and
finite value,

$\displaystyle\lambda(E)=\left(\ln|\frac{1+\sqrt{1-4\rho^{2}}}{2}\Big{|}+\cancel{2\pi\varepsilon}\right)-\left(\ln|\rho|+\cancel{2\pi\varepsilon}\vphantom{\sum}\right).$
(2.39)

Consequently, the Lyapunov exponent remains independent of energy. This
observation suggests that we can consider the entire spectrum as our energy
bin, and then calculate the average of $\tau$, which we will denote as
$\langle\tau\rangle$.

Through numerical simulations, we calculate the inverse participation ratio
with the energy bin encompassing the entire spectrum for lattice sizes
$L=2000,4000,\dots 12000$ in steps of $2000$. $R^{2}$ values are found to be
closely approximating $1$, as depicted in Fig. 2.2. Notably, the standard
error of $\langle\tau\rangle$ is approximately $10^{-5}$ with a magnitude one
thousand times smaller than $\langle\tau\rangle$ itself, as highlighted in
Figure 2.2 (c1). In Fig. 2.2 (c2), the variance of $\tau$s and its
fluctuations are extremely small compare to $\langle\tau\rangle$.

Figure 2.2:  (a) Two different localization length are plotted. The inverse of
the Lyapunov exponent ($1/\lambda$) [Eq. (2.39)] and the participation number
(PN) diverges as $\lambda_{2}/\lambda_{1}$ approaches to the transition point
$(\lambda_{2}/\lambda_{1})_{c}\approx 0.468$ (b) Plot of $1-R^{2}$. It shows
that data follows the linear model in Eq. (2.10) very well. (c1) and (c2) show
the detailed statistics at $\lambda_{2}/\lambda_{1}=0.625$.
$\mathbb{E}[\langle\tau\rangle]$ is an average of mean scaling exponent over
different realizations, and $\mathbb{E}[\operatorname{var}(\tau)]$ is an
average of its variance over different realizations.

Some important mathematical properties have been established in the existing
literature for the off-diagonal Harper model described in Eq. (2.23). Notably,
it has been proven that there is no absolute continuous spectrum [142],
implying the absence of extended states. Additionally, it is known that there
exists no pure point spectrum [124], signifying the nonexistence of localized
states. Based on these findings, we expect the Lyapunov exponent to be zero,
indicating critical behavior for the corresponding eigenstate. Explicit
calculations of the Lyapunov exponent for the off-diagonal Harper model are
possible. First, the complexified non-singular transfer matrix is defined as
follows with
$c_{\rho}(\varphi+i\varepsilon)=2\rho\cos(2\pi(\varphi+i\varepsilon))$. At
$\varepsilon\to\infty$, $B(E;\varphi+i\varepsilon)$ is further simplified,

$\displaystyle
B(E;\varphi+i\varepsilon)=\begin{bmatrix}E&c_{\rho}(\varphi-\alpha+i\varepsilon)\\\
c_{\rho}(\varphi+i\varepsilon)&0\end{bmatrix}\approx\rho
e^{2\pi\varepsilon}M_{\varepsilon}(\varphi,\varepsilon)$ (2.40) $\displaystyle
M_{\varepsilon}(\varphi,\varepsilon)=e^{-i2\pi(\varphi-\alpha/2)}\begin{bmatrix}\mathcal{O}\left(\exp(-\varepsilon)\right)&-e^{i\pi\alpha}\\\
e^{-i\pi\alpha}&0\end{bmatrix}.$ (2.41)

Then the corresponding complexified Lyapunov exponent is obtained as follows,

$\displaystyle\gamma_{\varepsilon}(E)=\lim_{n\to\infty}\frac{1}{2\pi
n}\int_{0}^{2\pi}\\!\\!\\!\\!d\theta\ln\lVert
B^{(n)}(E,\varphi+i\varepsilon)\rVert\approx\ln|\rho|+2\pi\varepsilon.$ (2.42)

Finally, the Lyapunov exponent of the off-diagonal Harper model at energy $E$
is zero. This result indicates a diverging localization length, signifying
critical behavior for the corresponding eigenvector,

$\displaystyle\lambda(E)=\left(\cancel{\ln|\rho|}+\cancel{2\pi\varepsilon}\vphantom{\sum}\right)-\left(\cancel{\ln|\rho|}+\cancel{2\pi\varepsilon}\vphantom{\sum}\right)=0.$
(2.43)

####  2.4.4 Case 1 – $\beta=0$ and $\lambda_{1}\neq\lambda_{2}$ at
$\theta_{1,2}=\pi/4$

Let us connect the projected model with the extended and off-diagonal Harper
models. We first start our mapping by choosing $\beta=0$ and
$\lambda_{1}\neq\lambda_{2}$. The onsite potential and hopping in Eq. (2.14)
are reduced to the following expressions, which gives us the extended Harper
model,

$\displaystyle v_{n}$
$\displaystyle=\left[\frac{1}{4}\left(\frac{\lambda_{2}}{\lambda_{1}}+1\right)\cos(\pi\alpha)\right]2\cos(2\pi\alpha
n-\pi\alpha),$ (2.44) $\displaystyle t_{n}$
$\displaystyle=\left[\frac{1}{4}\left(1-\frac{\lambda_{2}}{\lambda_{1}}\right)\right]\cos(2\pi\alpha
n).$ (2.45)

$2\rho$ in Eq. (2.22) can be either larger or smaller than $1$ depending on
the values of $\alpha$ and $\lambda_{2}/\lambda_{1}$ where the critical-to-
insulator transition is at $2\rho=1$. In this case, the formula for $2\rho$ is
easily obtained as follows,

$\displaystyle
2\rho=\left|\frac{1}{\cos(\pi\alpha)}\right|\left|\frac{1-\lambda_{2}/\lambda_{1}}{1+\lambda_{2}/\lambda_{1}}\right|.$
(2.46)

For $\alpha=(\sqrt{5}-1)/2$, the phase transition point is
$(\lambda_{2}/\lambda_{1})_{c}\approx 0.468$. As we have seen from Eq. (2.39),
the localization length remains independent of energy. Then, we can take the
entire spectrum as our energy bin and calculate the average of $\tau$, denoted
as $\langle\tau\rangle$. In Fig. 2.3, we compare the participation number and
the Lyapunov exponents of the projected model.

We point out that when the quasiperiodic amplitudes are equal
($\lambda_{1}=\lambda_{2}$) in, the hopping $t_{n}$ vanishes in the projected
model, which makes it diagonal and leaves only the non-zero onsite potential.
Hence, the eigenstates of the projected model exhibit compact localization
with eigenstates occupying a single site. Since the diagonal values are all
different, the eigenenergies have no degeneracy. In numerical perspective,
$R^{2}$ is not defined in this case.

Figure 2.3:  Phase diagram of the projected model based on the values of the
exponent $\tau$ averaged over the spectrum, $\langle\tau\rangle$, computed for
$25\times 25$ values of parameters $\lambda_{2}/\lambda_{1},\beta/\pi$. The
blue and red lines on the borders indicate the critical and localized regimes,
respectively. (a) $R^{2}$ values for parameters
$\lambda_{2}/\lambda_{1},\beta/\pi$. All points have values very close to $1$,
implying that all eigenstates are localized in the thermodynamic limit in this
entire region. (b) Slope $\hat{\kappa}$ (2.10) for the parameters away from
the border of the phase diagram. The absolute value of the slope increases
closer to the border, implying larger localization length.

####  2.4.5 Case 2 – $\lambda_{2}/\lambda_{1}=0$ for any phase difference
$\beta$

If we only apply the quasiperiodic potential (2.5) to one leg of the ladder,
e.g. for example $\lambda_{2}/\lambda_{1}=0$, then the projected model (2.14)
simplifies: $\tau_{s}$ is zero and $v_{s}$ takes the following form,

$\displaystyle
v_{s}=\sin(\pi\alpha)\cos(\theta_{1}-\theta_{2})\cos(\theta_{2}+\theta_{1}).$
(2.47)

Then, the model (2.14) maps to the extended Harper model for $v_{s}\equiv 0$
only, which produces the following constraint on the angles $\theta_{1,2}$,

$\displaystyle\theta_{2}+\theta_{1}=\frac{\pi}{2}.$ (2.48)

Using this to eliminate $\theta_{2}$, we obtain the following expressions for
$v_{c}$ and $\tau_{c}$,

$\displaystyle v_{c}$
$\displaystyle=2\cos(\pi\alpha)\sin^{2}\theta_{1}\cos^{2}\theta_{1},$ (2.49)
$\displaystyle t_{c}$ $\displaystyle=\sin\theta_{1}\cos\theta_{1}\sqrt{1-\cos
4\theta_{1}}/\sqrt{8}.$ (2.50)

Note that both quantities vanish for $\theta_{1}=0,\pi/2$, producing a set of
decoupled sites. Then for $\theta_{1}\neq 0,\pi/2$ the control parameter,
$\rho$ of the EHM (2.22) is expressed as follows in terms of $\theta_{1}$,

$\displaystyle 2\rho=\left|\frac{\sqrt{1-\cos
4\theta_{1}}/\sqrt{8}}{\cos(\pi\alpha)\sin\theta_{1}\cos\theta_{1}}\right|=\left|\frac{1}{\cos(\pi\alpha)}\right|\geq
1.$ (2.51)

The absolute value above ensures that the negative $t_{c}$ is also covered
(the sign can be trivially gauged away) and also to match with the definition
of the extended Harper model [124] where non-negative $2\rho$ is assumed.

Figure 2.4:  Phase diagram $(\theta_{1},\theta_{2})$ of the projected model
(2.14) for $\lambda_{2}/\lambda_{1}=0$, e.g. one of the quasiperiodic fields
is absent. The yellow dots indicate the fully localized spectrum. The diagonal
blue line indicates the region of critical spectrum given by Eq. (2.48). The
black dotted lines on the border correspond to the case of disconnected
dimers: all eigenstates are compactly localized. (a) $R^{2}$ values for
parameters $\theta_{1},\theta_{2}$: All eigenstates are localized in the
thermodynamic limit except the diagonal. (b) Slope $\hat{\kappa}$ (2.10) for
the parameters away from the diagonal of the phase diagram. The absolute value
of the slope $\kappa$ increases upon approaching the diagonal, implying a
larger localization length.

For all the other $\theta_{1,2}$, e.g. not satisfying the above condition
(2.48), we establish the character of the spectrum — localized or not —
numerically. We scanned the full parameter region, $0<\theta_{1,2}<\pi/2$,
discretized into the $25\times 25$ grid. The exponent $\langle\tau\rangle$
averaged over the entire spectrum is computed via Eq. 2.10 for lattice sizes
$L=2000,4000,\dots 12000$ in steps of $2000$. The results are summarized in
Fig 2.4: Apart from the diagonal line, all points have the $R^{2}$-values very
close to $1$. That is, all the eigenstates of the projected model (2.14) are
localized in the thermodynamic limit whenever it does not map onto the
extended Harper model, e.g., away from the diagonal
$\theta_{1}+\theta_{2}=\pi/2$. The increase of the absolute value of the slope
$\hat{\kappa}$ (2.10) towards the diagonal line is related to the increase of
the localization length as we approach the transition to critical states.

####  2.4.6 Case 3 – $\lambda_{2}=\lambda_{1}$ and $\beta=\pi$

For $\theta_{2}=\pi/4$, equal potential strengths $\lambda_{2}=\lambda_{1}=1$
and phase difference $\beta=\pi$, $v_{s,c}=t_{s}=0$ and we are left with only
the hopping terms $t_{c}=\sin 2\theta_{1}/2$ in the model (2.14, 2.16). The
projected model becomes precisely equivalent to the off-diagonal Harper model,
and the role of $t_{c}$ is simply to rescale the energy,

$\displaystyle Ea_{n}=t_{c}\left[\cos(2\pi\alpha
n)a_{n+1}+\cos(2\pi\alpha(n-1))a_{n-1}\right].$ (2.52)

Figure 2.5:  Phase diagram $(\theta_{1},\theta_{2})$ of the projected model
(2.14) for equal strength of the quasiperiodic potentials
$\lambda_{2}=\lambda_{1}$ and phase difference $\beta=\pi$. The yellow dots
indicate the fully localized spectrum. The horizontal blue line,
$\theta_{2}=\pi/4$, the models with critical spectrum, which map to the off-
diagonal Harper mode. The black dotted lines on the border correspond to the
system of disconnected dimers and compactly localized eigenstates. (a) $R^{2}$
values for parameters $\theta_{1},\theta_{2}$. Aside from the horizontal blue
line, all points have the values very close to $1$: all eigenstates are
localized in the thermodynamic limit. (b) Slope $\hat{\kappa}$ (2.10) for the
parameters away from the horizontal of the phase diagram. The absolute value
of the slope $\kappa$ increases upon approaching the diagonal, implying a
larger localization length.

For $\theta_{2}\neq\pi/4$, we carried out a numerical scan of the spectrum of
the projected model over the region $0<\theta_{1,2}<\pi/2$, discretized into
$25\times 25$ grid. The fitting procedure of $\langle\tau\rangle$ (over the
entire spectrum) is performed via Eq. 2.10 for lattice sizes
$L=2000,4000,\dots 12000$ in steps of $2000$. The observed results follow
closely those of the case $\lambda_{2}/\lambda_{1}=0$: all the eigenstates of
the projected model (2.14) are localized in the thermodynamic limit away from
the horizontal line of the phase diagram, e.g., when the projected model does
not map onto the off-diagonal Harper model, as shown in Fig. 2.5. The increase
of the absolute value of the slope towards the horizontal line is related to
the increase of the localization length towards the line of critical states.

####  2.4.7 Remark – Near zero angles $\theta_{1,2}\to 0$

In this section, we focus on the case of $\theta_{1,2}\to 0$ and demonstrate
that all the eigenstates of the projected model (2.14) are effectively
compactly localized. To see this we approximate $\sin\theta\approx\theta$ and
$\cos\theta\approx 1$, so that the coefficients (2.15-2.16) simplify to

$\displaystyle v_{s}$
$\displaystyle=\sin(\pi\alpha)(\theta_{1}^{2}\theta_{2}^{2}-1)$ (2.53)
$\displaystyle+\frac{\lambda_{2}}{\lambda_{1}}\cos\beta\sin(\pi\alpha)(\theta_{1}^{2}-\theta_{2}^{2})-\frac{\lambda_{2}}{\lambda_{1}}\sin\beta\cos(\pi\alpha)(\theta_{2}^{2}+\theta_{1}^{2})$
$\displaystyle v_{c}$
$\displaystyle=\cos(\pi\alpha)(1+\theta_{1}^{2}\theta_{2}^{2})$ (2.54)
$\displaystyle+\frac{\lambda_{2}}{\lambda_{1}}\sin\beta\sin(\pi\alpha)(\theta_{1}^{2}-\theta_{2}^{2})+\frac{\lambda_{2}}{\lambda_{1}}\cos\beta\cos(\pi\alpha)(\theta_{2}^{2}+\theta_{1}^{2})$
$\displaystyle t_{s}$
$\displaystyle=\theta_{1}\theta_{2}\frac{\lambda_{2}\sin\beta}{\lambda_{1}}$
(2.55) $\displaystyle t_{c}$
$\displaystyle=\theta_{1}\theta_{2}\left(1-\frac{\lambda_{2}\cos\beta}{\lambda_{1}}\right)$
(2.56)

Keeping only linear terms in $\theta_{1,2}$, so that
$\theta_{1}^{2}=\theta_{2}^{2}=\theta_{1}\theta_{2}=0$, we see that the
hoppings all vanish $t_{s}=t_{c}=0$ and only the onsite potential terms are
left. Then, all the eigenstates are effectively compactly localized, but their
energies are non-degenerate because of the onsite potential. A similar
argument implies that all states are compactly localized for
$\theta_{1,2}\approx\pi/2$.

###  2.5 Wavepacket spreading in projected model

The critical states are exhibited in both the extended and off-diagonal Harper
models. Nonetheless, during our investigation of their transport properties,
there are notable qualitative distinctions in the eigenstates of each model.
Firstly, we provide an overview of the established theoretical findings
concerning the lower and upper bounds of diffusion exponents in quasiperiodic
systems. Subsequently, we see the numerical results of the projected models at
$\theta_{1,2}=\pi/4$, which are mapped to the extended and the off-diagonal
Harper model.

####  2.5.1 Lower bound of diffusion exponent

The diffusion exponent $\gamma$ is obtained from the deviation of the
particle’s position at time $t$, $\hat{X}_{\phi}(t)$. The explicit definition
is given below [143, 144, 145], and we keep the dimension of a system to be
one since we only deal with one-dimensional systems,

$\displaystyle\lim_{t\to\infty}\langle
0||\hat{X}_{\phi}(t)-\hat{X}(0)|^{2}|0\rangle\sim t^{2\gamma},$ (2.57)
$\displaystyle\hat{X}_{\phi}(t)=\exp(iH_{\phi}t)\hat{X}\exp(-iH_{\phi}t).$
(2.58)

The cases of $\gamma=1$ signify ballistic motion, $\gamma=1/2$ corresponds to
regular diffusion, and $\gamma=0$ represents a localization of an initial
wavepacket. When $\gamma$ is between $0$ and $1/2$, it represents a
subdiffusion. If $\gamma$ is between $1/2$ and $1$, it corresponds to a
superdiffusion. In the context of a one-dimensional lattice with a fractal-
like energy spectrum [146, 147, 148], there exists a rigorous lower limit for
the diffusion exponent, denoted as $\gamma^{-}$, which falls within
$0<\gamma^{-}\leq 1/2$.

####  2.5.2 Upper bound of diffusion exponent

Concerning the Aubry-André-Harper model at its critical point, the energy
spectrum has a fractal-like structure. In this scenario, some numerical
observations claimed numerically that the corresponding diffusion exponent is
almost near to $1/2$, without any theoretical proof [143, 149]. Using the
translational covariance relation that quasiperiodic systems have [150],
$U_{n}H_{\phi}U_{n}^{-1}=H_{\tau^{n}\phi}$, it is possible to obtain the upper
bound of the corresponding diffusion exponent. The operator $\tau$ is a group
action that constructs an orbit, and $U_{n}$ is a shift operator which shifts
to the $+n$ site. $H(\theta)$ is the Hamiltonian. In the case of the projected
model, $\tau^{n}\phi$ equals $\phi+\alpha n$. Furthermore, the two-point
correlation measure strongly connects with the current-current correlation
measure $m$ [150],

$\displaystyle
m(dE,dE^{\prime})=\frac{1}{L}\sum_{n,m}\delta(E-E_{n})\delta(E^{\prime}-E_{m})|\langle
E_{n}|J|E_{m}\rangle|^{2}dEdE^{\prime},$ (2.59)

where $J=i[H,x]=\nabla H$. In Ref. [151], the Stieltjes transform of the
current-current correlation $m$ allows us to calculate the diffusion exponent.

$\displaystyle
S_{m}(z_{1},z_{2})\approx\int_{\mathbb{R}^{2}}\frac{m(dE,dE^{\prime})}{(E-z_{1})(E^{\prime}-z_{2})}.$
(2.60)

The theorem 10 in Ref. [151] states that the real value of the Stieltjes form
$S$ and the diffusion exponent are interrelated as follows,

$\displaystyle\lim_{\varepsilon\to
0^{+}}\operatorname{Re}\left(\int_{\mathbb{R}^{2}}\\!\\!\\!\\!daS_{m}(a+i\varepsilon,a-i\varepsilon)\right)\sim\varepsilon^{1-2\gamma}.$
(2.61)

By means of the above theorem, we can find the upper bound by assuming
$|\langle E|[H,x]|E^{\prime}\rangle|$ is bounded by some positive constant
$M$,

$\displaystyle\int_{\mathbb{R}}S_{m}(a+i\varepsilon,a-i\varepsilon)da<\int_{\mathbb{R}^{3}}\frac{MdEdE^{\prime}da}{(E-(a+i\varepsilon))(E^{\prime}-(a-i\varepsilon))}.$
(2.62)

Applying Taylor expansion on the integrand with respect to small $\varepsilon$
and taking the real part gives us a constant leading order term,

$\displaystyle\real\left(\int_{\mathbb{R}^{3}}\frac{MdEdE^{\prime}da}{(E-(a+i\varepsilon))(E^{\prime}-(a-i\varepsilon))}\right)\sim\epsilon^{0}.$
(2.63)

Then, we obtain an upper bound $\gamma^{+}$ of the diffusion exponent of the
projected model under the assumption. Hence, the actual diffusion exponent
$\gamma$ is smaller than or equal to $\gamma^{+}=1/2$.

####  2.5.3 Numerical results

To further quantify the transport properties, we analyze the spreading of an
initially localized wavepacket for weak quasiperiodic perturbation. We use an
initial state localized on a single site in the center of the lattice. For
convenience, we assign the lattice center to be the zero coordinate. To
quantify the spreading of the initial state, we compute the root-mean-square
of the displacement $\sigma(t)$ where $\langle n\rangle$ is the average
position,

$\displaystyle\sigma(t)=\left[\sum_{n\in\mathbb{Z}}(n-\langle
n\rangle(t))^{2}\absolutevalue{\psi_{n}(t)}^{2}\right]^{1/2}\\!\\!\\!\\!\\!\\!\\!\\!\propto
t^{\gamma},\hskip 5.0pt\text{where }\langle
n\rangle(t)=\sum_{n\in\mathbb{Z}}n\absolutevalue{\psi_{n}(t)}^{2}.$ (2.64)

We also average the results over $40$ values of the phase $\phi$ (2.5) sampled
from the range $[0,\pi]$. The displacement $\sigma(t)$ provides the deviation
of the particle’s position from its average at time $t$. The absence of
spreading indicates localization. Also, the spreading stops in finite systems
once the boundaries are reached. At an intermediate time, i.e., a period
before the boundaries are reached, $\sigma(t)$ is fitted by a power law with
the exponent $\gamma$ whose value indicates the type of transport: diffusive,
subdiffusive, or ballistic [152].

$(\lambda_{2}/\lambda_{1},\beta/\pi)$ |  $\gamma\pm\Delta\gamma$ |  $(\lambda_{2}/\lambda_{1},\beta/\pi)$ |  $\gamma\pm\Delta\gamma$
---|---|---|---
$(0.00,0.00)$ |  $0.34\pm 0.01$ |  $(0.00,1.00)$ |  $0.34\pm 0.01$
$(0.15,0.00)$ |  $0.38\pm 0.01$ |  $(0.25,1.00)$ |  $0.39\pm 0.04$
$(0.25,0.00)$ |  $0.39\pm 0.04$ |  $(0.50,1.00)$ |  $0.37\pm 0.01$
$(0.35,0.00)$ |  $0.41\pm 0.03$ |  $(0.75,1.00)$ |  $0.39\pm 0.02$
$(0.46,0.00)$ |  $0.41\pm 0.03$ |  $(1.00,1.00)$ |  $0.50\pm 0.01$
Table 2.1:  Diffusion exponent $\gamma$ in Eq. (2.64) for various points on
the border of the phase diagram in Fig. 2.3. Except for
$\lambda_{2}/\lambda_{1}=1$ and $\beta=\pi$, all results show that the
critical states support clear subdiffusive transport. For the off-diagonal
Harper model, $\beta=\pi,\lambda_{1}=\lambda_{2}$, the diffusion exponent is
either diffusive or subdiffusive but very close to diffusive, similar to the
case of the Aubry-André-Harper model at the critical point [143, 149].  Figure
2.6:  Positive correlation between the exponent $\gamma$ in Eq. (2.64) (red
points) and the exponent $\tau$ averaged over the spectrum,
$\langle\tau\rangle$, in Eq. (2.10) (blue points). The values of $\gamma$ are
taken from Table 2.1. (a) $\beta=\pi$. (b) $\beta=0$.

We sampled the boundary of the phase diagram given in Fig. 2.3, which
corresponds to the extended Harper model and features critical eigenstates
(for localized cases, the spreading stops as discussed above). The results of
the fitting are provided in Table 2.1, and the details of the wavepacket
spreading are shown in Fig. 2.7 for system size $L=12801$. In all cases, we
see a clear subdiffusion [153, 154, 155], whose exponent $\gamma$ depends on
the position on the border of the phase diagram, i.e., the values of
$\lambda_{2}/\lambda_{1},\beta/\pi$. For $\lambda_{2}/\lambda_{1}=0$, $\beta$
drops out of the Hamiltonian in Eq. (2.14), and consequently, the diffusion
exponent $\gamma$ does not depend on $\beta$. Furthermore, we observe a
positive correlation between $\gamma$ and $\langle\tau\rangle$ in Fig. 2.6,
which suggests that the transport properties are strongly affected by the
details of the profiles of the critical eigenstates.

In Sec. 2.4.5, we observe that when the condition
$\theta_{1}+\theta_{2}=\pi/2$ holds, the projected model becomes equivalent,
up to an overall scaling factor, to the specific scenario where
$\theta_{1,2}=\pi/4$. Consequently, we can infer that for values of
$\theta_{1,2}$ satisfying Eq. (2.48), the wavepacket spreading exhibits
subdiffusive behavior. Similarly, in Sec. 2.4.6, when we keep
$\theta_{2}=\pi/4$, we obtain the off-diagonal Harper model. Then, the factor
$t_{c}$ in Eq. (2.52) serves as the global timescale for the almost diffusive
spreading.

Figure 2.7:  Log-log plots of the spreading of a wavepacket initially
localized on a single site using $\log_{10}$. The order of the plots follows
the data in Table 2.1.

###  2.6 Finite perturbation

In the limit of weak interaction, we saw the emergence of entirely critical
spectra for specific values of the local unitary transformation parameters
$\theta_{1,2}$, the relative potential strength $\lambda_{2}/\lambda_{1}$, and
the phase difference $\beta$. As we increase the perturbation strength and
make it finite, we expect the system to localize for large enough values of
$\lambda_{1,2}$. The open question is what happens to the critical states for
moderate values of $\lambda_{1,2}$ and whether any extended states might
emerge. Also, we do not expect any change in the localization properties of
the states that were already localized for weak quasiperiodic perturbation.

For finite potential strengths, the projected model description is no longer
valid; thus, we have to focus on the original Hamiltonian,
$\mathcal{H}=\mathcal{H}_{\mathrm{ABF}}+W$. Then, the original Hamiltonian can
be expressed as a non-Abelian Aubry-André-Harper model [156] with a
quasiperiodic block $A_{n}$ and a fixed hopping block $B$, where they are not
commutative in general: $[A_{n},B]\neq 0$,

$\displaystyle E\psi_{n}=A_{n}\psi_{n}+B\psi_{n-1}+B^{\dagger}\psi_{n+1},$
(2.65)

where $\psi_{n}=(p_{n},f_{n})^{T}$, and $A_{n}$ is a block matrix with the
quasiperiodic fields and $\theta_{1,2}$ and $B$ is a matrix only depends on
$\theta_{1,2}$. As an example, when $\theta_{1,2}=\pi/4$, the original model
looks like as follows.

$\displaystyle A_{n}=\frac{4}{\Delta}\begin{bmatrix}W_{1}(n)&0\\\
0&W_{2}(n)\end{bmatrix},\quad B=\begin{bmatrix}1&-1\\\ 1&-1\end{bmatrix}.$
(2.66)

This section focuses on the parameter regions where the projected model [Eq.
(2.14)] for weak disorder hosts critical states.

####  2.6.1 Representation in semi-engtangled basis

Sometimes it is convenient to work with the _semi-detangled_ Hamiltonian,
[105, 91, 119, 120] which is defined by inverting only the second unitary
transformation $U_{2}$ (see Sec. 2.2). This defines the semi-detangled basis
$\\{u_{n},d_{n}\\}$ and gives the new, semi-detangled Hamiltonian,

$\displaystyle\mathcal{H}_{\mathrm{SD}}=U_{1}\mathcal{H}_{\mathrm{FD}}U_{1}^{\dagger}+U_{2}^{\dagger}WU_{2}.$
(2.67)

The semi-detangled wavefunction amplitudes $u_{n}$ and $d_{n}$ are related to
the basis states of the ABF Hamiltonian $\\{p_{n},f_{n}\\}$ as follows,

$\displaystyle
u_{n}=p_{n}\cos\theta_{2}-f_{n}\sin\theta_{2}\quad\text{and}\quad
d_{n}=p_{n}\sin\theta_{2}+f_{n}\cos\theta_{2}.$ (2.68)

The model (2.67) takes the following form in the semi-detangled basis:

$\displaystyle Eu_{n}$
$\displaystyle=\left[\varepsilon_{b}\sin^{2}\theta_{1}+\varepsilon_{a}\cos^{2}\theta_{1}\right]u_{n}+\left[W_{2}(n)\sin^{2}\theta_{2}+W_{1}(n)\cos^{2}\theta_{2}\right]u_{n}$
$\displaystyle+\left[(W_{1}(n)-W_{2}(n))\cos\theta_{2}\sin\theta_{2}\right]d_{n}+\left[\Delta\cos\theta_{1}\sin\theta_{1}\right]d_{n+1},$
(2.69) $\displaystyle Ed_{n}$
$\displaystyle=\left[\varepsilon_{a}\sin^{2}\theta_{1}+\varepsilon_{b}\cos^{2}\theta_{1}\right]d_{n}+\left[W_{1}(n)\sin^{2}\theta_{2}+W_{2}(n)\cos^{2}\theta_{2}\right]d_{n}$
$\displaystyle+\left[(W_{1}(n)-W_{2}(n))\cos\theta_{2}\sin\theta_{2}\right]u_{n}+\left[\Delta\cos\theta_{1}\sin\theta_{1}\right]u_{n-1}.$
(2.70)

Throughout the section, we consider the system in semi-detangled basis.

####  2.6.2 Critical-to-insulator transition in finite perturbation

Let us set $\theta_{2}=\pi/4$ and $\theta_{1}=\theta$. In the case when both
quasiperiodic fields have same amplitude $\lambda_{2}=\lambda_{1}=\lambda$ and
the phase difference $\beta=\pi$, we have simply $W_{2}=-W_{1}=-W$. With these
conditions, the semi-detangled Hamiltonian in Eq. (2.6.1) takes the following
form,

$\displaystyle E_{u}(E,\theta)u_{n}$
$\displaystyle=W(n)d_{n}+\Delta\cos\theta\sin\theta d_{n+1},$ (2.71)
$\displaystyle E_{d}(E,\theta)d_{n}$
$\displaystyle=W(n)u_{n}+\Delta\cos\theta\sin\theta u_{n-1},$ (2.72)

where $E_{u}(E,\theta)$ and $E_{d}(E,\theta)$ are defined as

$\displaystyle E_{u}(E,\theta)$
$\displaystyle=E-\varepsilon_{b}\sin^{2}\theta-\varepsilon_{a}\cos^{2}\theta,$
(2.73) $\displaystyle E_{d}(E,\theta)$
$\displaystyle=E-\varepsilon_{a}\sin^{2}\theta-\varepsilon_{b}\cos^{2}\theta.$
(2.74)

We can eliminate one of the amplitudes $u_{n}$ or $d_{n}$ by substituting one
of the equations into the other. If we choose to eliminate $u_{n}$, we get the
following effective equation, resulting in equations with $d_{n}$ only:

$\displaystyle E_{d}(\theta)E_{u}(\theta)d_{n}$
$\displaystyle=W^{2}(n)d_{n}+\Delta^{2}\cos^{2}\theta\sin^{2}\theta d_{n}$
$\displaystyle+\Delta\cos\theta\sin\theta\left[W(n)d_{n+1}+W(n-1)d_{n-1}\right].$
(2.75)

By rearranging the above equation, we get the following eigenequation,

$\displaystyle\tilde{E}(\theta,\lambda)d_{n}$ $\displaystyle=\cos(4\pi\alpha
n)d_{n}$
$\displaystyle+K(\theta,\lambda))\left[\cos(2\pi\alpha(n-1))d_{n-1}+\cos(2\pi\alpha
n)d_{n+1}\right],$ (2.76)

where $\tilde{E}(\theta,\lambda)$ and $K(\theta,\lambda)$ are defined as,

$\displaystyle\tilde{E}(\theta,\lambda)=\frac{2}{\lambda^{2}}E_{d}(\theta)E_{u}(\theta)-\frac{1}{2}\quad\text{and}\quad
K(\theta,\lambda)=\frac{2\Delta}{\lambda}\cos\theta\sin\theta.$ (2.77)

This eigenequation looks very similar to the extended Harper model, Eq.
(2.22), except for the spatial frequency of the onsite potential that is
double that of the hopping. Also, the potential is twice as large as in the
extended Harper model.

For instance, when both $\theta_{1}$ and $\theta_{2}$ are equal to $\pi/4$,
Eq. (2.6.2) is reduced into the following equation:

$\displaystyle\tilde{E}d_{n}=\cos(4\pi\alpha
n)d_{n}\vphantom{\frac{2t}{\lambda}}+\frac{2t}{\lambda}\left[\cos(2\pi\alpha(n-1))d_{n-1}+\cos(2\pi\alpha
n)d_{n+1}\right],$ (2.78)

where the eigenvalue $\tilde{E}$ is given as

$\displaystyle\tilde{E}=\frac{2}{\lambda^{2}}\left[(E-\varepsilon)^{2}-t^{2}-\frac{\lambda^{2}}{2}\right].$
(2.79)

The non-Abelian Aubry-André-Harper model for this case has been studied
numerically [157, 156]. The critical eigenstates over the entire spectrum show
up until $\lambda$ reaches the value of the flatband bandgap
$\Delta=|\varepsilon_{b}-\varepsilon_{a}|$. Then, the CIT occurs once
$\lambda$ is equal to $\Delta$ for the entire spectrum. We remark that this
result coincides with the extended Harper model estimate of the transition,
$2t/\lambda=1$ if we neglect the differences between our model, Eq. (2.6.2),
and the true extended Harper model, Eq. (2.22). Then, we get the CIT for the
generalized case in Eq. (2.6.2) when $K(\theta,\lambda)=1$ as well.

An example for $\theta=0.1\pi/2$ is shown and summarized in Fig. 2.9. The
flatband energies are taken to be $\varepsilon_{a}=-1$ and
$\varepsilon_{b}=2$. The entire spectrum is rescaled for each $\lambda_{1}$.
Then we apply the linear model introduced in Eq. (2.10) introducing $50$
energy bins $\tilde{e}$ for the spectrum rescaled to fit in $[0,1]$. Lattice
sizes (in unit cells) are $L=3283,4181,4832,5473,6765,10946$. The exact model
has two sites for each unit cell, and the size of the Hamiltonian matrix is
$2L\times 2L$. The transition point $\lambda_{c}\approx 0.927$ based on Eq.
(2.77) is indicated with the black dashed line. It matches perfectly with the
numerical results. All states are critical to the left of $\lambda_{c}$ and
localized to the right of it. The string of eigenstates between the two
broadened flatband energies is due to the open boundary conditions.

####  2.6.3 Non-degenerate compact localized states

Suppose both $\theta_{1}$ and $\theta_{2}$ are $\pi/4$. When $\beta$ is zero
and $\lambda=\lambda_{1}=\lambda_{2}$, we obtain the following Hamiltonian in
the semi-detangled basis $u_{n},d_{n}$,

$\displaystyle\tilde{E}u_{n}$ $\displaystyle=\lambda\cos(2\pi\alpha
n)u_{n}+td_{n+1},$ (2.80) $\displaystyle\tilde{E}d_{n}$
$\displaystyle=\lambda\cos(2\pi\alpha n)d_{n}+tu_{n-1},$ (2.81)

where $\tilde{E}=E-\varepsilon$. Again, we can eliminate one of the amplitudes
through substitution. Choosing to eliminate $u_{n}$ leads to a diagonal
problem for $d_{n}$:

$\tilde{E}d_{n}=\left[\lambda\cos(2\pi\alpha
n)+\frac{t}{\tilde{E}-\lambda\cos(2\pi\alpha(n-1))}\right]d_{n},$ (2.82)

and a compact localization of the eigenstates for any $\lambda$. We recall
that we have already found compact localization for weak perturbation. This
compact localization survives for any potential strength. When $\theta_{1}$
and $\theta_{2}$ are arbitrary values, we again obtain non-degenerate compact
localization with more complicated $\tilde{E}$.

####  2.6.4 Non-simple fractality edges

Figure 2.8:  Critical-to-insulator transition at finite potential strength
$\lambda_{2}=\lambda_{1}$, $\beta=\pi$ at $\theta=0.1\pi/2$. (a) Rescaled
spectrum. The negative $R^{2}$ values are replaced with $-1$ for better visual
separation of localized and critical regions. (b) Original, non-rescaled
energy spectrum for lattice size $L=10946$. The vertical dashed black line
marks the CIT transition predicted by Eq. (2.77).

Figure 2.9:  (a) Fractality edges in the exact model at
$\beta/\pi=1,\lambda_{2}/\lambda_{1}=0.5$. Every $R^{2}$ lower than or equal
to zero is revalued to $-1$ for the purpose of clear distinction between the
localized and critical regions in the figure. The $-1$ $R^{2}$ value implies
that states are critical. For $R^{2}\approx 1$, the states are localized. (b)
Original, non-rescaled energy spectrum for lattice size $L=10946$.

We now consider all the other values of $\beta,\lambda_{1,2}$ where we
observed critical states for weak perturbation. In the semi-detangled basis
$u_{n},d_{n}$, we obtain the following Hamiltonian,

$\displaystyle\tilde{E}u_{n}$
$\displaystyle=\frac{W_{1}(n)+W_{2}(n)}{2}u_{n}+\frac{W_{1}(n)-W_{2}(n)}{2}d_{n}+td_{n+1},$
$\displaystyle\tilde{E}d_{n}$
$\displaystyle=\frac{W_{1}(n)+W_{2}(n)}{2}d_{n}+\frac{W_{1}(n)-W_{2}(n)}{2}u_{n}+tu_{n-1}.$

$\tilde{E}$ is $E-\varepsilon$ and $W_{1,2}$ are the quasiperiodic fields. The
above effective equations are a mixture of the one given in Sec. 2.6.2 and
Sec. 2.6.3, implying that an eigenstate’s behavior is at most critical.

To look into the properties of the states, we have to resort to numerical
analysis on the full model in Eq. (2.66) to analyze their localization
properties. Flatband energies are fixed to $\varepsilon_{a}\\!=\\!-1$ and
$\varepsilon_{b}\\!=\\!2$ with bandgap $\Delta=3$. We consider several values
of $\beta,\lambda_{1,2}$ to probe the different regions of the phase diagram
with critical states and repeat the computation of the IPR of the eigenstates.
Then we apply the linear model introduced in Eq. (2.10) and introduce $50$
energy bins for the eigenenergy range rescaled to fit between $0$ and $1$,
which is referred to as the rescaled energy. Then, it allows us to distinguish
localized and critical states as a function of eigenenergy. The lattice sizes
(in unit cells) are $L\\!=\\!2584,3283,4181,4832,5473$, and $6765$. Since the
exact model has two sites for each unit cell, the size of the Hamiltonian
matrix is $2L\times 2L$.

The generic observation is that the critical spectrum of the projected model
is replaced with a mixed one, comprising partially critical and partially
localized states depending on the eigenenergy. We dub the border between
critical and localized states as _fractality edges_ by analogy with mobility
edges separating localized and extended states [126, 127, 128, 129]. For
instance, let $\beta\\!=\\!\pi$ and $\lambda_{2}/\lambda_{1}\\!=\\!0.5$.
Figure 2.9 shows the emergence of fractality edges for a finite strength of
perturbation $\lambda_{1,2}$. We note that the critical states extend to
values $\lambda_{1}\\!\geq\\!\Delta$, at variance with the case in Sec. 2.6.2
where all states localize for $\lambda_{1}\\!>\\!\Delta$.

####  2.6.5 Perturbation independent fractality edges for $\lambda_{2}=0$

Figure 2.10:  Fractality edges at finite potential strength $\lambda_{1}$ and
$\lambda_{2}=0$, for $\theta=0.01\pi/2$. (a) Fractality edges in the rescaled
spectrum. (b) Fractality edges in the original, non-rescaled energy spectrum
with lattice size $L=10946$. The bottom black line is the flatband
$\varepsilon_{a}=-1$, while the top black black is $E\approx-0.997$. The
dashed black line is the upper bound $E\approx-0.998$ given by Eq. (2.88).

Figure 2.11:  Fractality edges at finite potential strength $\lambda_{1}$ and
$\lambda_{2}=0$, for $\theta=0.68\pi/2$. (a) Fractality edges in the rescaled
spectrum. (b) Fractality edges in the original, non-rescaled energy spectrum
with lattice size $L=10946$. The bottom black line is the flatband
$\varepsilon_{a}=-1$, while the top black line is $E\approx 1.05$. The dashed
black line is the upper bound $E\approx 0.393$ given by Eq. (2.88).

Now, let us consider a case where we remove one of the quasiperiodic
potentials, $\lambda_{2}=0$ and set
$\theta_{2}=\pi/2-\theta,\theta_{1}=\theta$. The model in the semi-detangled
basis takes the following form,

$\displaystyle\frac{Eu_{n}}{\cos\theta\sin\theta}$
$\displaystyle=\left[\frac{\varepsilon_{a}}{\tan\theta}+\frac{\left(\varepsilon_{b}+W_{1}(n)\right)}{1/\tan\theta}\right]u_{n}+W_{1}(n)d_{n}+\Delta
d_{n+1},$ (2.83) $\displaystyle\frac{Ed_{n}}{\cos\theta\sin\theta}$
$\displaystyle=\left[\frac{\left(\varepsilon_{b}+W_{1}(n)\right)}{\tan\theta}+\frac{\varepsilon_{a}}{1/\tan\theta}\right]d_{n}+W_{1}(n)u_{n}+\Delta
u_{n-1}.$ (2.84)

After multiplying Eq. (2.84) by $\tan\theta$ and taking the difference of the
above equations, we get a single equation without quasiperiodic potentials:

$\displaystyle E_{u}(E,\theta)u_{n}+\tan\theta E_{g}u_{n-1}=\tan\theta
E_{d}(E,\theta)d_{n}+E_{g}d_{n+1},$

where $E_{u}(E,\theta)$ and $E_{d}(E,\theta)$ are

$\displaystyle E_{u}(E,\theta)$
$\displaystyle=\left[\frac{E}{\cos\theta\sin\theta}-\left(\frac{E_{a}}{\tan\theta}+E_{b}\tan\theta\right)\right],$
(2.85) $\displaystyle E_{d}(E,\theta)$
$\displaystyle=\left[\frac{E}{\cos\theta\sin\theta}-\left(\frac{E_{b}}{\tan\theta}+E_{a}\tan\theta\right)\right].$
(2.86)

For $\theta_{1,2}=\pi/4$, we obtain a simpler model. There are fractality
edges at the flatband energies $\varepsilon_{a,b}$ for any perturbation
strength $\lambda_{1}$ and all the eigenstates in between are critical:

$\displaystyle(2E-\sigma)u_{n}+E_{g}u_{n-1}=(2E-\sigma)d_{n}+E_{g}d_{n+1},$
(2.87)

where $\sigma$ is defined as $\varepsilon_{a}+\varepsilon_{b}$. This result
can be rationalized with a simple assumption that the critical states appear
when the hopping is larger than the onsite potential, $|2E-\sigma|\leq E_{g}$,
implying $\varepsilon_{a}\leq E\leq\varepsilon_{b}$ for the critical states.
We can apply the same logic for $\theta\neq\pi/4$. Then, the above inequality
for energy $E$ is modified as follows,

$\displaystyle|E_{u}|<E_{g}\tan\theta\quad\text{and}\quad|E_{d}|\tan\theta<E_{g}.$
(2.88)

However, the simple constraint does not explain the entire range of fractality
edges but only covers the narrower range of the critical states. Examples are
shown and summarized in Fig. 2.11 and 2.11, where $\theta=0.01\pi/2$ and
$\theta=0.68\pi/2$, respectively. We fixed the flatband energies to
$\varepsilon_{a}=-1$ and $\varepsilon_{b}=2$. The narrowness of the critical
region in Fig. 2.11 is due to the smallness of the angle $\theta_{1}$;
remember that for $\theta_{1}=0$ we have trivial flatbands due to the model
decoupling into disconnected dimers. Since the fractality edges are potential
strength independent, we only considered part of the spectrum around the edges
and used it for the rescaled spectrum. We apply the linear model introduced in
Eq. (2.10), and introduce $50$ energy bins $e$ for the spectrum rescaled to
fit in $[0,1]$. The lattice sizes (in unit cells) are
$L=3283,4181,4832,5473,6765,10946$. The exact model has two sites for each
unit cell, and the size of the Hamiltonian matrix is $2L\times 2L$. In both
examples, two fractality edges are independent of the potential strength
$\lambda_{1}$, with critical states in between and localized states outside.
The lower fractality edge coincides with the flatband energy
$\varepsilon_{a}$. By symmetry, setting $\lambda_{1}=0$ and varying
$\lambda_{2}$, one would get the similar spectral behavior with the upper
fractality edge would be at $\varepsilon_{b}$.

###  2.7 Conclusion

In this chapter, we summarized our study on the effect of quasiperiodic
perturbation, composed of two potentials, on a general model of the ABF
manifold with two bands. First, we considered the case of weak quasiperiodic
perturbation. We identified the ABF submanifolds of models with spectra
containing critical states and exhibiting subdiffusive and almost diffusive
transport for weak quasiperiodic perturbation. The critical states with
subdiffusive transport were found for the ABF models satisfying the condition
$\theta_{1}+\theta_{2}=\pi/2$ for the angles of the unitary transformation
generating the model. The other constraint is the absence of one of the
quasiperiodic fields. The other submanifold of models featuring critical
states with almost diffusive transport is given by $\theta_{2}=\pi/4$ and
arbitrary $\theta_{1}$. The quasiperiodic potentials must have equal strengths
and the relative phase $\beta=\pi$, e.g., be of opposite signs. All such
models can be mapped onto the off-diagonal Harper model. Outside of these
manifolds, all the states are localized in the thermodynamic limit, as
suggested by numerical analysis. Specifically, for small angles
$\theta_{1,2}$, we demonstrated compact localization of the spectrum.

For a finite quasiperiodic potential, the perturbed ABF models, which map to
the extended Harper model, display fractality edges in the spectrum separating
critical from localized states. These edges are independent of the potential
strength. On the other hand, the models that map onto the off-diagonal Harper
model show no fractality edges. Instead, they exhibit a critical-to-insulator
transition. The transition point is derived analytically and confirmed
numerically and depends on the angle $\theta_{1}=\theta$ with fixed
$\theta_{2}=\pi/4$.

The presented ABF manifold, perturbed by quasiperiodic potentials, could be
potentially implemented using ultracold atomic gases loaded onto optical
lattice potentials where ABF networks have been realized, e.g. in the Creutz
ladder [158, 159] or the diamond chain [93, 160]. A second promising platform
is light propagation in photonic lattices where the diamond chain with its
Aharonov-Bohm cages was experimentally obtained [161]. Another promising
future platform is electric circuits. The main issue for the experimental
realization of ABF systems is achieving both positive and negative hopping. A
possible candidate is a stack of LC electric circuits [162] in which the
topology of wiring can generate positive and negative hopping, even complex
hopping. We also expect by adding extra LC resonators to each site might
create the onsite quasiperiodic modulation of the ABF system.

## Chapter 3 Trapped Hard-core Bosons
in Flatband Systems

###  3.1 Introduction

Recently, there has been considerable interest in physical systems with
macroscopic degeneracies. One particular example of such a system is a
flatband system [163, 49, 164] with compact localized states, where the
eigenstates are confined to a finite number of sites [11, 17], which is a
result of destructive interferences caused by the network geometry as we have
mentioned in Chapter 1. The extreme sensitivity of macroscopically degenerate
flatbands to perturbations and interactions gives rise to a diverse range of
interesting and exotic phases. One example is ergodicity breaking by
considering a fine-tuned interaction [98, 99, 105, 106, 100]. In this
contenxt, flatband systems offer a platform for investigating intriguing
phenomena related to quantum ergodicity breaking.

Thermalization has fascinated physicists as it describes the evolution of
quantum many-body systems from reversible microscopic dynamics toward
equilibrium. One intriguing aspect is the tendency of all pure states within a
specific energy shell to exhibit thermal-like behavior [165]. In search of an
explanation, an eigenstate thermalization hypothesis (ETH) was proposed [166],
suggesting that the thermalization in isolated quantum systems can be
attributed to the assumption that every energy eigenstate possesses thermal
properties. ETH has been widely discussed, but weak ETH emerged where _nearly_
all energy eigenstates exhibit thermal properties to some extent [167]. It
should be noted that weak ETH alone cannot determine the presence or absence
of thermalization for physically realistic initial states [167]. Nevertheless,
weak ETH does assure that an initial state lacking significant overlap with
rare states will undergo thermalization. It is known that weak ETH holds even
for a wide range of translation-invariant systems, regardless of their
integrability [168, 169]. However, many-body localization induced by the
disorder can violate weak ETH [170].

This chapter summarizes the work on the thermalization behavior of hard-core
bosons in the one- and two-dimensional cross-stitch lattice with on-site
repulsive interaction having the strongest collision effect. It is known that
the presence of macroscopic degeneracy in the flatband energy allows for the
amplification of the effects of interactions and perturbations. We
specifically focus on the interplay between compact localized states (CLSs)
and the extreme limit of strong repulsion imposed by hard-core constraints.
Our findings reveal the band insulating phase, the Wigner crystal phases, and
the violation of weak ETH. The non-ergodic behaviors observed in the absence
of disorder are of particular significance, highlighting the presence of
strictly non-ergodic excited states. This non-ergodic behavior is closely
linked to the concept of Hilbert space fragmentation. The related work is in
preparation.

The chapter is organized as follows. We start by defining one- and two-
dimensional cross-stitch lattices and briefly mention the properties of hard-
core bosons in Sec. 3.2. In Sec. 3.3, we explore the process of site filling
to obtain the band insulating phase and Wigner crystal phases. Moving on to
Sec. 3.4, we investigate the non-ergodic excited states in a closed CLS
barrier. We briefly point out that, in the diamond lattice, similar phenomena
occur as the cross-stitch lattice in Sec. 3.5. Then, we make a conclusion of
the chapter in the last section.

###  3.2 Models

This section defines the one- and two-dimensional cross-stitch lattices, which
this chapter focuses on. Additionally, the properties of hard-core bosons are
going to be discussed.

####  3.2.1 1D and 2D cross-stitch lattices

The Hamiltonian for the one-dimensional cross-stitch chain, which includes
vertical hopping denoted as $t$, is expressed as follows, with the intra-unit
cell hopping term $\hat{h}_{n}$ and the inter-unit cell hopping term
$\hat{v}_{n}$ defined accordingly,

$\displaystyle\mathcal{H}=-\sum_{n}\hat{h}_{n}+\hat{h}_{n}^{\dagger}+t(\hat{v}_{n}+\hat{v}_{n}^{\dagger}),$
(3.1)
$\displaystyle\hat{h}_{n}=(\hat{a}_{n}^{\dagger}+\hat{b}_{n}^{\dagger})(\hat{a}_{n+1}+\hat{b}_{n+1})\quad\text{and}\quad\hat{v}_{n}=\hat{a}_{n}^{\dagger}\hat{b}_{n}.$
(3.2)

Here, $\hat{a}_{n}^{\dagger}$ and $\hat{b}_{n}^{\dagger}$ represent the
creation operators at site $A_{n}$ and $B_{n}$, respectively. When $t$ is
sufficiently small, it results in a ground state characterized by a flatband.
The corresponding energy bands are depicted in Fig. 3.2. Likewise, the two-
dimensional cross-stitch lattice can be defined as

$\displaystyle\mathcal{H}=-\sum_{n,m}\hat{h}_{n,m}+\hat{h}_{n,m}^{\dagger}+t(\hat{v}_{n,m}+\hat{v}_{n,m}^{\dagger}),$
(3.3)
$\displaystyle\hat{h}_{n,m}=(\hat{a}_{n,m}^{\dagger}\\!+\hat{b}_{n,m}^{\dagger})(\hat{a}_{n+1,m}\\!+\hat{b}_{n+1,m}\\!+\hat{a}_{n,m+1}\\!+\hat{b}_{n,m+1}),$
(3.4) $\displaystyle\hat{v}_{n,m}=\hat{a}_{n,m}^{\dagger}\hat{b}_{n,m}.$ (3.5)

####  3.2.2 Properties of hard-core bosons

A hard-core boson is a particle that obeys a set of mixed commutation
relations, such that
$[\hat{c}_{i},\hat{c}_{j}]=[\hat{c}_{i}^{\dagger},\hat{c}_{j}^{\dagger}]=[\hat{c}_{i},\hat{c}_{j}^{\dagger}]=0,$
for all $i\neq j$, and exhibits fermion-like behavior at the same site, as
described by the anti-commutation relation
$\\{\hat{c}_{i},\hat{c}_{i}\\}=\\{\hat{c}_{i}^{\dagger},\hat{c}_{i}^{\dagger}\\}=0$
and $\\{\hat{c}_{i},\hat{c}_{i}^{\dagger}\\}=1$. $\hat{c}_{n}^{\dagger}$ is a
hard-core boson creation operator at site $n$. If we merge all the commutation
properties, we get the following single expression,

$\displaystyle[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=\delta_{ij}(1-2\hat{a}_{i}^{\dagger}\hat{a}_{i}).$
(3.6)

These commutation relations define a bosonic system characterized by strong
repulsion at short distances, making it theoretically unfavorable for states
with more than two particles to occupy a single site. Such particle can be
achieved at low densities and low temperature in certain regimes of large
scattering length in quasi-one-dimensional systems [171, 172, 173]. It is
known that, in the one-dimensional hard-core bosonic gas system, the true
condensation at zero temperature cannot be achieved [174, 175, 176]. Instead,
due to the presence of a quasi-long-range order in the system, quasi-
condensates at finite momentum can emerge [177, 178]. In three-dimensional
systems, Mott insulating phase and Bose-Einstein condensation can exist [179].
Our study focuses specifically on the on-site interaction, disregarding the
interaction term.

###  3.3 CLS groundstates

This section demonstrates two findings on the one-dimensional cross-stitch
lattice. Firstly, we establish that occupying CLSs results in the groundstate
configuration. Additionally, we illustrate that fully-filled dimers do not
contribute to the overall groundstate eigenenergy. The subsequent discussion
will provide a comprehensive understanding of these observations.

####  3.3.1 $\nu\leq 1/2$ – filling of the flatband CLS

First, consider a CLS in the one-dimensional cross-stitch lattice at the
$m$-th unit cell. To describe the ground state of the flatband, we propose an
ansatz, where $P_{\mathrm{FB}}$ represents the set of all CLSs occupying
certain unit cells,

$\displaystyle\ket{\mathrm{GS}(P_{\mathrm{FB}})}=\\!\\!\\!\prod_{n\in
P_{\mathrm{FB}}}\\!\\!\\!\ket{\mathrm{CLS}}_{n}\quad\text{and}\quad\ket{\mathrm{CLS}}_{m}=\frac{\hat{a}_{m}^{\dagger}-\hat{b}_{m}^{\dagger}}{\sqrt{2}}\ket{\varnothing}.$
(3.7)

Applying the ansatz to the Hamiltonian gives the following filling procedure,

$\displaystyle\hat{h}_{n}\ket{\mathrm{CLS}}_{m}=\hat{h}_{n}^{\dagger}\ket{\mathrm{CLS}}_{m}=0,$
(3.8)
$\displaystyle(\hat{v}_{n}+\hat{v}_{n}^{\dagger})\ket{\mathrm{CLS}}_{m}=-\delta_{nm}\ket{\mathrm{CLS}}_{m}.$
(3.9)

Then, the proposed ansatz in Eq. (3.7) represents the ground state. The energy
of the flatband groundstate is derived as follows,

$\displaystyle E_{\mathrm{GS}}=t|P_{\mathrm{FB}}|=2t\nu N,$ (3.10)

where $\nu$ is a filling fraction which is $\nu\leq 1/2$, and $N$ is the total
number of unit cells. Moreover, $E_{\mathrm{GS}}$ is macroscopically
degenerated. At $\nu=1/2$, CLSs are completely occupied, and the mutual
repulsion from hard-core boson constraint prevents hopping to other sites.
Hence, a Wigner crystal configuration is attained.

####  3.3.2 $\nu\leq 1$ – filling the $a_{n},b_{n}$ dimers

For hard-core bosons, when the number of particles in some flatband systems
exceeds the critical filling fraction of the flatband-induced Wigner crystal,
an extra particle forms a pair with an existing CLS, as discussed in Ref.
[180, 110]. In the context of the cross-stitch lattices, introducing an extra
particle above the critical filling fraction has the same effect as replacing
a CLS with a fully-filled dimer, as expressed mathematically below,

$\displaystyle\hat{a}^{\dagger}(\hat{a}^{\dagger}-\hat{b}^{\dagger})\propto\hat{a}^{\dagger}\hat{b}^{\dagger}\propto\hat{b}^{\dagger}(\hat{a}^{\dagger}-\hat{b}^{\dagger}).$
(3.11)

Now, let us replace some of these CLSs with fully-filled dimers. Let $P_{d}$
denote the set of unit cells occupied by these fully-filled dimers. Then the
flatband groundstate $\ket{\mathrm{GS}(P_{\mathrm{FB}})}$ is composed of
states that are not elements of $P_{d}$. From that, our new groundstate ansatz
is proposed as follows,

$\displaystyle\ket{\mathrm{GS}}=\prod_{k\in
P_{d}}\hat{a}_{k}^{\dagger}\hat{b}_{k}^{\dagger}\ket{\mathrm{GS}(P_{\mathrm{FB}})}.$
(3.12)

The size of $P_{d}$ is expressed with the filling fraction $\nu$ where
$1/2<\nu\leq 1$ and $N$ is total number of unit cells, $|P_{d}|=(2\nu-1)N$.

Next, let us consider the filling procedure for the $n$-th unit cell belongs
to $P_{d}$ which include cases either $n+1\in P_{d}$ or $n+1\notin P_{d}$. In
either scenario, the ansatz in Eq. (3.12) represents an eigenstate since the
following relations hold,

$\displaystyle\hat{h}_{n}\ket{\mathrm{GS}}=\hat{h}_{n}^{\dagger}\ket{\mathrm{GS}}=0\quad\text{and}\quad\hat{v}_{n}\ket{\mathrm{GS}}=\hat{v}_{n}^{\dagger}\ket{\mathrm{GS}}=0.$
(3.13)

The energy $E_{\mathrm{GS}}$ is determined solely by the contribution of CLSs
and exhibits macroscopic degeneracy,

$\displaystyle E_{\mathrm{GS}}=t(N-|P_{d}|)=2tN(1-\nu).$ (3.14)

At $\nu=1$, every site is entirely occupied by hard-core bosons, resulting in
an insulating phase due to the lack of available states for conduction. Hence,
we see a band insulator phase. Otherwise, we observe a combination of the band
insulator and the Wigner crystal phases. The diagram of $E_{\mathrm{GS}}$
depending on $\nu$ is drawn in Fig. 3.2.

Figure 3.1:  (Left) Energy bands of 1D cross-stitch lattice. Here we set
$t=-5$. (Right) Changing of the groundstate energy as a function of a filling
fraction $\nu$. Up to $\nu=1/2$, we see the decreasing $E_{\mathrm{GS}}$. Then
we have increasing $E_{\mathrm{GS}}$.

Figure 3.2:  The time evolution of the initial wavefunction with four hard-
core bosons on the one-dimensional cross-stitch lattice is plotted. Here, we
set $t=-5$. At time zero, two CLSs are located at the second and the 11th unit
cells (yellow) with two hard-core bosons positioned at the seventh and the
15th sites (red). The CLSs are fixed as time evolves and two hard-core bosons
are strictly confined.

###  3.4 Non-ergodic excitation in cross-stitch lattices

This section presents our finding of non-ergodic excitations on the cross-
stitch lattices. Initially, we examine the scenario where one hard-core boson
is situated between two CLSs. Next, we map our model onto a spin-1 chain. From
that, we aim to explain the non-ergodicity of bosonic excitations.

####  3.4.1 One hard-core boson between two CLSs

Let us examine the scenario where two CLSs are positioned at the $\pm N$-th
unit cells. Then, we introduce a single hard-core boson between these two
CLSs. For instance, we consider the case where the boson is created at the
$a_{q}$-th site,

$\displaystyle\ket{\psi}_{1}=\hat{a}_{q}^{\dagger}\ket{\mathrm{FB}},\quad\ket{\mathrm{FB}}=\ket{\mathrm{CLS}}_{N}\otimes\ket{\mathrm{CLS}}_{-N},\quad\text{where
}|q|<N.$ (3.15)

Applying the Hamiltonian to $\ket{\psi}_{1}$, we encounter four distinct cases
that require evaluation. In the first case, when $n\neq\pm N$ and $n\neq q$,
all hoppings become zero as no particles are present at those sites,

$\displaystyle\hat{h}_{n}\ket{\psi}_{1}=\hat{h}_{n-1}\ket{\psi}_{1}=\hat{h}_{n}^{\dagger}\ket{\psi}_{1}=\hat{h}_{n-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.16)
$\displaystyle\hat{v}_{n}\ket{\psi}_{1}=\hat{v}_{n}^{\dagger}\ket{\psi}_{1}=0.$
(3.17)

In the second case, when $n=q$ and $q\pm 1\neq\pm N$, the term
$\hat{v}_{q}\ket{\psi}_{1}$ is zero as no particle is created at the
$b_{q}$-th site and no multiple counting is allowed at the $a_{q}$-th site.
Similarly, acting with $\hat{h}_{q}$ and $\hat{h}_{q-1}^{\dagger}$ on
$\ket{\psi}_{1}$ yields zero values because no particles are created at the
$a_{q\pm 1}$-th and $b_{q\pm 1}$-th sites. For $\hat{h}_{q-1}$,
$\hat{h}_{q}^{\dagger}$ and $\hat{v}_{q}^{\dagger}$, we do observe particle
hops to nearby vacant sites,

$\displaystyle\hat{h}_{q}\ket{\psi}_{1}=\hat{h}_{q-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.18)
$\displaystyle\hat{h}_{q-1}\ket{\psi}_{1}=(\hat{a}_{q-1}^{\dagger}+\hat{b}_{q-1}^{\dagger})\ket{\mathrm{FB}},$
(3.19)
$\displaystyle\hat{h}_{q}^{\dagger}\ket{\psi}_{1}=(\hat{a}_{q+1}^{\dagger}+\hat{b}_{q+1}^{\dagger})\ket{\mathrm{FB}},$
(3.20)
$\displaystyle\hat{v}_{q}\ket{\psi}_{1}=0\quad\text{and}\quad\hat{v}_{q}^{\dagger}\ket{\psi}_{1}=\hat{b}_{q}^{\dagger}\ket{\mathrm{FB}}.$
(3.21)

Moving to the third case, $n=q$ and $q+1=N$, applying $\hat{h}_{q}$ and
$\hat{h}_{q}^{\dagger}$ on $\ket{\psi}_{1}$ results in zero due to the
properties of a CLS. Similarly, acting with $\hat{h}_{q-1}^{\dagger}$ on
$\ket{\psi}_{1}$ results in zero because no particles are present at the
$a_{q-1}$-th and $b_{q-1}$-th sites. Again, acting with $\hat{v}_{q}$ on
$\ket{\psi}_{1}$ also gives zero, since no particle exists at the $b_{q}$-th
site and multiple counting is prohibited. This same logic applies when
$q-1=-N$,

$\displaystyle\hat{h}_{q}\ket{\psi}_{1}=\hat{h}_{q}^{\dagger}\ket{\psi}_{1}=\hat{h}_{q-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.22)
$\displaystyle\hat{h}_{q-1}\ket{\psi}_{1}=(\hat{a}_{q-1}^{\dagger}+\hat{b}_{q-1}^{\dagger})\ket{\mathrm{FB}},$
(3.23)
$\displaystyle\hat{v}_{q}\ket{\psi}_{1}=0\quad\text{and}\quad\hat{v}_{q}^{\dagger}\ket{\psi}_{1}=\hat{b}_{q}^{\dagger}\ket{\mathrm{FB}}.$
(3.24)

In the fourth case, $n=N$, acting with
$\hat{h}_{N}$,$\hat{h}_{N-1}$,$\hat{h}_{N}^{\dagger}$, and
$\hat{h}_{N-1}^{\dagger}$ on $\ket{\psi}_{1}$ yield zero values, which are
attributed to the nature of a CLS. Additionally, $\ket{\psi}_{1}$ is the
eigenstate of $\hat{v}_{N}+\hat{v}_{N}^{\dagger}$. Similar results are
observed in the case of $n=-N$,

$\displaystyle\hat{h}_{N}\ket{\psi}_{1}=\hat{h}_{N-1}\ket{\psi}_{1}=\hat{h}_{N}^{\dagger}\ket{\psi}_{1}=\hat{h}_{N-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.25)
$\displaystyle(\hat{v}_{N}+\hat{v}_{N}^{\dagger})\ket{\psi}_{1}=-\ket{\psi}_{1}.$
(3.26)

In summary, the presence of CLSs at both ends leads to eliminating the lattice
outside of the boundary, resulting $\ket{\mathrm{FB}}$ as an effective vacuum
state. Hence, we are left with an effective finite Hamiltonian
$\tilde{\mathcal{H}}$, which describes the system confined within two CLSs.

$\displaystyle\tilde{\mathcal{H}}\ket{\psi}_{1}$
$\displaystyle=-(t\hat{b}_{q}^{\dagger}\\!+\hat{a}_{q+1}^{\dagger}\\!+\hat{b}_{q+1}^{\dagger}\\!+\hat{a}_{q-1}^{\dagger}\\!+\hat{b}_{q-1}^{\dagger})\ket{\mathrm{FB}},$
(3.27) $\displaystyle\tilde{\mathcal{H}}\ket{\psi}_{1}$
$\displaystyle=(t\hat{b}_{q}^{\dagger}+\hat{a}_{q-1}^{\dagger}+\hat{b}_{q-1}^{\dagger})\ket{\mathrm{FB}},\quad
q+1=N$ (3.28) $\displaystyle\tilde{\mathcal{H}}\ket{\psi}_{1}$
$\displaystyle=(t\hat{b}_{q}^{\dagger}+\hat{a}_{q+1}^{\dagger}+\hat{b}_{q+1}^{\dagger})\ket{\mathrm{FB}},\quad
q-1=-N$ (3.29)

To illustrate this, let us consider a scenario where two CLSs are created with
three vacant unit cells between them. Then, we add one particle in one of the
vacant unit cells. In this configuration, the resulting effective Hamiltonian
becomes the truncated version of the original Hamiltonian, representing the
confined system within this boundary. The effective vacuum state
$|\tilde{\varnothing}\rangle$ is $\ket{\mathrm{FB}}$,

$\displaystyle\tilde{H}\cong\begin{bmatrix}0&t&1&1&0&0\\\ t&0&1&1&0&0\\\
1&1&0&t&1&1\\\ 1&1&t&0&1&1\\\ 0&0&1&1&0&t\\\
0&0&1&1&t&0\end{bmatrix}\quad\text{and}\quad|\tilde{\varnothing}\rangle=\ket{\mathrm{FB}}.$
(3.30)

Since a single hard-core boson is unable to escape the confinement imposed by
the CLSs, we can expect that when the number of bosons is less than the number
of sites within the unoccupied unit cells, these bosons located in the vacant
unit cells will also be trapped, and unable to leak out. Fig. 3.2 presents an
example of this confinement, depicting two hard-core bosons residing within
eight vacant unit cells of the one-dimensional cross-stitch lattice. The
numerical simulation in Fig. 3.2 clearly shows that the hard-core bosons
become trapped in the specified region, unable to move beyond the boundaries
set by the CLSs as time evolves.

####  3.4.2 Mapping to spin-1 chain

To provide a more rigorous explanation for the non-ergodic excited states
involving multiple particles in the one-dimensional many-body cross-stitch
system, we establish a mapping of the hard-core bosons to a spin-1 chain. To
accomplish this mapping, we introduce the following ladder operators:
$S_{n}^{-}=\hat{a}_{n}^{\dagger}+\hat{b}_{n}^{\dagger}$ and
$S_{n}^{+}=\hat{a}_{n}+\hat{b}_{n}$. After performing the necessary
calculations, we derive the expression for $S^{z}$ in the following form,

$\displaystyle S_{z}=1-(a_{n}^{\dagger}a_{n}+b_{n}^{\dagger}b_{n}).$ (3.31)

Then, we get the following triplet states at the $n$-th unit cell,

$\displaystyle\ket{+}_{n}:=|S_{n}=1,m_{n}$
$\displaystyle=+1\rangle=\ket{\varnothing},$ (3.32)
$\displaystyle\ket{0_{t}}_{n}:=|S_{n}=1,m_{n}$
$\displaystyle=\hphantom{+}0\rangle=\frac{\hat{a}_{n}^{\dagger}+\hat{b}_{n}^{\dagger}}{\sqrt{2}}\ket{\varnothing},$
(3.33) $\displaystyle\ket{-}_{n}:=|S_{n}=1,m_{n}$
$\displaystyle=-1\rangle=\sqrt{2}\hat{a}^{\dagger}_{n}\hat{b}^{\dagger}_{n}\ket{\varnothing}.$
(3.34)

To clarify, the state $\ket{+}$ corresponds to an empty particle state
$\ket{\varnothing}$ in the language of a hard-core boson at the $n$-th unit
cell in the one-dimensional cross-stitch lattice. Furthermore, CLSs are
represented as singlet states $\ket{0_{s}}$,

$\displaystyle
S^{-}_{n}\ket{\mathrm{CLS}}_{n}=0=S^{+}_{n}\ket{\mathrm{CLS}}_{n}.$ (3.35)

Then we can rewrite the one-dimensional cross-stitch Hamiltonian (Eq. (3.1),
but in the many-body context) in spin notation as follows,

$\displaystyle\mathcal{H}$
$\displaystyle=-\\!\\!\sum_{n}S_{n}^{-}S_{n+1}^{+}\\!+S_{n+1}^{-}S_{n}^{+}\\!+t(S_{n}^{-}S_{n}^{+}\\!\\!+S_{n}^{z}-1)$
(3.36) $\displaystyle=-\\!\\!\sum_{\langle
i,j\rangle}\\!\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}\right)\\!-\\!t\\!\sum_{n}\\!\left((S_{n}^{z})^{2}-S_{n}^{z}-1\right).$
(3.37)

First, let us examine the (ground) state $\ket{0_{s}\cdots 0_{s}}$, which
occurs when CLSs are present in every unit cell. This results in the outcome
described in Eq. (3.10),

$\displaystyle\mathcal{H}\ket{0_{s}\cdots 0_{s}}=tL\ket{0_{s}\cdots 0_{s}}.$
(3.38)

Let us focus on the following state $\ket{+\cdots+}$ full of triplet states.
This state corresponds to every site being empty, leading to
$\mathcal{H}\ket{+\cdots+}=0\ket{+\cdots+}$. Next, suppose we have a sequence
$\Omega$ of triplet states and a linear combination of triplet and singlet
states. Then, for the given initial state,
$\ket{+\cdots+0_{s},\Omega,0_{s}+\cdots+}$, $\ket{+}$ and $\ket{0_{s}}$ remain
unchanged and $\Omega$ hops to the another possible sequence
$\Omega^{\prime}$. The system described by the Hamiltonian in Eq. (3.37)
corresponds to a spin-1 XY model [181] (and a spin-1/2 XYZ creutz ladder
[182]). What is particularly intriguing about this model is the existence of
quantum many-body scars [183, 181, 184] and Hilbert space fragmentation [182,
181, 184]. The pivotal role played by CLSs involves the decomposition of the
Hilbert space into disjoint Krylov subspaces, leading to the emergence of a
truncated system, as shown in Fig. 3.2. However, it is important to note that
the trapped bosons within the CLS barriers in Fig. 3.2 do not constitute a
quantum many-body scar because the state is a linear combination of triplet
and singlet states. Instead, their presence gives rise to localized-like
states induced by the CLS singlet state.

####  3.4.3 Localized excitations in 2D lattices

Same phenomena as in the one-dimensional cross-stitch lattice are also
observed in the two-dimensional cross-stitch model. The idea of non-ergodic
excited states in the one-dimensional cross-stitch lattice expands to a 2D
scenario, resulting in the confinement of particles in a closed loop formed by
CLSs. The calculations for this case are straightforward and equivalent to
those performed for the one-dimensional cross-stitch lattice.

Consider a configuration where multiple CLSs are arranged in a closed loop
denoted as $\mathcal{S}$. Now, let us create a single boson at the
$a_{q,k}$-th site and examine whether this boson can escape the confines of
the loop,

$\displaystyle\ket{\psi}_{1}=\hat{a}_{q,k}^{\dagger}\ket{\mathrm{FB}},\quad\ket{\mathrm{FB}}=\\!\\!\\!\\!\prod_{(i,j)\in\mathcal{S}}\\!\\!\\!\frac{\hat{a}_{i,j}^{\dagger}-\hat{b}_{i,j}^{\dagger}}{\sqrt{2}}\ket{\varnothing}.$
(3.39)

Applying the Hamiltonian on $\ket{\psi}_{1}$, we get four distinct cases that
need to be considered. In the first case, when $(n,m)\notin\mathcal{S}$ and
$(n,m)\neq(q,k)$, all hoppings become zero since there are no particles
present at those sites,

$\displaystyle\hat{h}_{n,m}\ket{\psi}_{1}=\hat{h}_{n-1,m}\ket{\psi}_{1}=\hat{h}_{n,m-1}\ket{\psi}_{1}=0,$
(3.40)
$\displaystyle\hat{h}_{n,m}^{\dagger}\ket{\psi}_{1}=\hat{h}_{n-1,m}^{\dagger}\ket{\psi}_{1}=\hat{h}_{n,m-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.41)
$\displaystyle\hat{v}_{n,m}\ket{\psi}_{1}=\hat{v}_{n,m}^{\dagger}\ket{\psi}_{1}=0$
(3.42)

In the second case, when $(n,m)=(q,k)$ and none of the neighboring sites
$(q+1,k),(q,k+1),(q-1,k),(q,k-1)$ are elements of $\mathcal{S}$, the term
$\hat{v}_{q,k}\ket{\psi}_{1}$ is zero. This results from the absence of a
particle at site $b_{q,k}$-th site and the principle that no multiple counting
is allowed. Similarly, acting with $\hat{h}_{q,k}$ and
$\hat{h}_{q-1,k}^{\dagger}$ on $\ket{\psi}_{1}$ results zero because there are
no particles created at sites $a_{q\pm 1,k}$-th and $b_{q\pm 1,k}$-th sites.
For $\hat{h}_{q,k}^{\dagger}$ and $\hat{v}_{q,k}^{\dagger}$, we obtain the
hoppings to the nearby empty sites,

$\displaystyle\hat{h}_{q,k}\ket{\psi}_{1}=\hat{h}_{q-1,k}^{\dagger}\ket{\psi}_{1}=\hat{h}_{q,k-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.43)
$\displaystyle\hat{h}_{q-1,k}\ket{\psi}_{1}=(\hat{a}_{q-1,k}^{\dagger}+\hat{b}_{q-1,k}^{\dagger})\ket{\mathrm{FB}},$
(3.44)
$\displaystyle\hat{h}_{q,k-1}\ket{\psi}_{1}=(\hat{a}_{q,k-1}^{\dagger}+\hat{b}_{q,k-1}^{\dagger})\ket{\mathrm{FB}},$
(3.45)
$\displaystyle\hat{h}_{q,k}^{\dagger}\ket{\psi}_{1}\\!=(\hat{a}_{q+1,k}^{\dagger}\\!\\!+\hat{b}_{q+1,k}^{\dagger}\\!\\!+\hat{a}_{q,k+1}^{\dagger}\\!\\!+\hat{b}_{q,k+1}^{\dagger})\ket{\mathrm{FB}},$
(3.46)
$\displaystyle\hat{v}_{q,k}\ket{\psi}_{1}=0\quad\text{and}\quad\hat{v}_{q,k}^{\dagger}\ket{\psi}_{1}=\hat{b}_{q,k}^{\dagger}\ket{\mathrm{FB}}$
(3.47)

Now, let us move to the third case. It corresponds to the case when
$(n,m)=(q,k)$, and we designate some or all of the unit cells at $(q+1,k)$,
$(q,k+1)$, $(q-1,k)$, and $(q,k-1)$ as elements of $\mathcal{S}$. Then, the
movement of particles to other sites depends on the configuration of
$\mathcal{S}$. For particles initially placed within $\mathcal{S}$, they
remain confined within it. Likewise, $\hat{v}_{q,k}\ket{\psi}_{1}$ equals zero
because there are no particles at site $b_{q,k}$-th site and no multiple
counting is allowed,

$\displaystyle\left\\{\begin{aligned}
&\hat{h}_{q,k}\ket{\psi}_{1},\quad\hat{h}_{q-1,k}\ket{\psi}_{1},\quad\hat{h}_{q,k-1}\ket{\psi}_{1}\\\
&\hat{h}_{q,k}^{\dagger}\ket{\psi}_{1},\quad\hat{h}_{q-1,k}^{\dagger}\ket{\psi}_{1},\quad\hat{h}_{q,k-1}^{\dagger}\ket{\psi}_{1}\end{aligned}\right\\}=\text{Depends
on the shape }\mathcal{S},$ (3.48)
$\displaystyle\hat{v}_{q,k}\ket{\psi}_{1}=0\quad\text{and}\quad\hat{v}_{q,k}^{\dagger}\ket{\psi}_{1}=\hat{b}_{q,k}^{\dagger}\ket{\mathrm{FB}}.$
(3.49)

In the fourth case, where $(n,m)\in\mathcal{S}$, acting with
$\hat{h}_{n,m}$,$\hat{h}_{n-1,m}$,$\hat{h}_{n,m-1}$, and their conjugates on
$\ket{\psi}_{1}$ yields zero values. This is attributed to the nature of a
CLS. Additionally, $\ket{\psi}_{1}$ serves as an eigenstate of
$\hat{v}_{n,m}+\hat{v}_{n,m}^{\dagger}$,

$\displaystyle\hat{h}_{n,m}\ket{\psi}_{1}=\hat{h}_{n-1,m}\ket{\psi}_{1}=\hat{h}_{n,m-1}\ket{\psi}_{1}=0,$
(3.50)
$\displaystyle\hat{h}_{n,m}^{\dagger}\ket{\psi}_{1}=\hat{h}_{n-1,m}^{\dagger}\ket{\psi}_{1}=\hat{h}_{n,m-1}^{\dagger}\ket{\psi}_{1}=0,$
(3.51)
$\displaystyle(\hat{v}_{n,m}+\hat{v}_{n,m}^{\dagger})\ket{\psi}_{1}=-\ket{\psi}_{1}$
(3.52)

In brief, the presence of the CLS loop $(n,m)\in\mathcal{S}$ effectively
erases again the region outside of $\mathcal{S}$. Hence, $\ket{\mathrm{FB}}$
behaves as vacuum state, yielding an effective Hamiltonian
$\tilde{\mathcal{H}}$ that characterizes the system confined within
$\mathcal{S}$.

The non-ergodic excitation can again be explained using the spin-1
representation, where $\mathbf{n}$ and $\mathbf{m}$ are unit-cell indices, and
the angular bracket means that the indices correspond to the nearest neighbor,

$\displaystyle\mathcal{H}=\\!\\!\sum_{\langle\mathbf{n},\mathbf{m}\rangle}S_{\mathbf{n}}^{-}S_{\mathbf{m}}^{+}+S_{\mathbf{n}}^{-}S_{\mathbf{m}}^{+}-\sum_{\mathbf{n}}t(S_{\mathbf{n}}^{-}S_{\mathbf{n}}^{+}+S_{\mathbf{n}}^{z}-1),$
(3.53)