# Inserting Information Bottlenecks for Attribution in Transformers

Zhiying Jiang, Raphael Tang, Ji Xin  Jimmy Lin
David R. Cheriton School of Computer Science
University of Waterloo
{zhiying.jiang, r33tang, ji.xin<EMAIL_ADDRESS>
###### Abstract

Pretrained transformers achieve the state of the art across tasks in natural
language processing, motivating researchers to investigate their inner
mechanisms. One common direction is to understand what features are important
for prediction. In this paper, we apply information bottlenecks to analyze the
attribution of each feature for prediction on a black-box model. We use BERT
as the example and evaluate our approach both quantitatively and
qualitatively. We show the effectiveness of our method in terms of attribution
and the ability to provide insight into how information flows through layers.
We demonstrate that our technique outperforms two competitive methods in
degradation tests on four datasets. Code is available at
https://github.com/bazingagin/IBA.

##  1 Introduction

Increasingly prominent is the urge to interpret deep neural networks, with the
success of these black-box models remaining vastly inexplicable both
theoretically and empirically. Within natural language processing (NLP), this
desire is particularly true for the pretrained transformer, which has
witnessed an influx of literature on interpretability analysis. Such papers
include visualizing transformer attention mechanisms Kovaleva et al. (2019),
probing the geometry of transformer representations Hewitt and Manning (2019),
and explaining the span predictions of question answering models van Aken et
al. (2019).

In this paper, we focus on prediction attribution methods. That is, we ask,
“Which hidden features contribute the most toward a prediction?” To resolve
this question, a number of methods Selvaraju et al. (2017); Smilkov et al.
(2017) generate attribution scores for features, which provide a human-
understandable “explanation” of how a particular prediction is made at the
instance level. Specifically, given an instance, these methods assign a
numerical score for each hidden feature denoting its relevance toward the
prediction.

Previous papers have demonstrated that gradient-based methods fail to capture
all the information associated with the correct prediction Li et al. (2016).
To address this weakness, Schulz et al. (2020) insert information bottlenecks
Tishby et al. (2000) for attribution, attaining both stronger empirical
performance and a theoretical upper bound on the information used.
Additionally, mutual information is unconstrained by model and task Guan et
al. (2019). Thus, we adopt information bottlenecks for attribution (IBA) to
interpret transformer models at the instance level. We apply IBA to BERT
Devlin et al. (2019) across five datasets in sentiment analysis, textual
entailment, and document classification. We show both qualitatively and
quantitatively that the method capably captures information in the model’s
token-level features, as well as insight into cross-layer behavior.

Our contributions are as follows: First, we are the first to apply information
bottlenecks (IB) for attribution to explain transformers. Second, we conduct
quantitative analysis to investigate the accuracy of our method compared to
other interpretability techniques. Finally, we examine the consistency of our
method across layers in a case study. Across four datasets, our technique
outperforms integrated gradients (IG) and local interpretable model-agnostic
explanations (LIME), two widely adopted prediction attribution approaches.

##  2 Related Work

In terms of scope, interpretability methods can be categorized as model
specific or model agnostic. Model-specific methods interpret only one family
of models, whereas

model-agnostic techniques aim for wide applicability across many families of
parametric models. We can roughly separate model-agnostic methods into three
categories: (1) gradient-based ones Li et al. (2016); Fong and Vedaldi (2017);
Sundararajan et al. (2017); (2) probing Ribeiro et al. (2016); Lundberg and
Lee (2017); Tenney et al. (2019); Clark et al. (2019); Liu et al. (2019); (3)
information-theoretical methods Bang et al. (2019); Guan et al. (2019); Schulz
et al. (2020); Pimentel et al. (2020) .

Gradient-based methods are, however, limited to models with differentiable
neural activations. They also fail to capture all the information associated
with the correct prediction Li et al. (2016). Although probing methods provide
detailed insight into specific models, they fail to capture inner mechanisms
like how information flows through the network Guan et al. (2019).
Information-theoretic methods, in contrast, provide consistent and flexible
explanations, as we show in this paper.

Guan et al. (2019) use mutual information to interpret NLP models across
different tokens, layers, and neurons, but they lack a quantitative
evaluation. Bang et al. (2019) also propose a model-agnostic interpretable
model using IB; however, they limit the information through the network by
sampling a given number of words at the beginning, which restricts the
explanation to neurons only. Our method is inspired by Schulz et al. (2020),
who use IBA in image classification.

##  3 Method

The idea of IBA is to restrict the information flowing through the network for
every single instance, such that only the most useful information is kept.

Concretely, given an input $\mathbf{X}\in\mathbb{R}^{N}$ and output
$\mathbf{Y}\in\mathbb{R}^{M}$, an information bottleneck is an intermediate
representation $\mathbf{T}$ that maximizes the following function:

$\operatorname{I}(\mathbf{Y};\mathbf{T})-\beta\cdot\operatorname{I}(\mathbf{X};\mathbf{T}),$
(1)

where $\operatorname{I}$ denotes mutual information and $\beta$ controls the
trade-off between reconstruction $\operatorname{I}(\mathbf{Y};\mathbf{T})$ and
information restriction $\operatorname{I}(\mathbf{X};\mathbf{T})$. The larger
the $\beta$, the narrower the bottleneck, i.e., less information is allowed to
flow through the network.

We insert the IB after a given layer $l$ in a pretrained deep neural network.
In this case, $\mathbf{X}=f_{l}(\mathbf{H})$ represents the chosen layer’s
output, where $\mathbf{H}$ is the input of the layer. We restrict information
flow by injecting noise into the original input:

$\mathbf{T}=\bm{\mu}\odot\mathbf{X}+(\mathbf{1}-\bm{\mu})\odot\epsilon,$ (2)

where $\odot$ denotes element-wise multiplication, $\epsilon$ the injected
noise, $\mathbf{X}$ the latent representation of the chosen layer,
$\mathbf{1}$ the all-one vector, and $\bm{\mu}\in\mathbb{R}^{N}$ the weight
balancing signal and noise. For every dimension $i$, $\bm{\mu}_{i}\in[0,1]$,
meaning that when $\bm{\mu}_{i}=1$, there is no noise injected into the
original representation. To simplify the training process, we set
$\bm{\mu}_{i}=\sigma(\bm{\alpha}_{i})$,

where $\sigma$ is the sigmoid function and $\bm{\alpha}$ is a learnable
parameter vector. In the extreme case, where all the information in
$\mathbf{T}$ is replaced with noise ($\mathbf{T}=\epsilon$), it’s desirable to
keep $\epsilon$ the same mean and variance as $\mathbf{X}$ in order to
preserve the magnitude of the input to the following layer. Thus, we have
$\epsilon\sim\mathcal{N}(\mu_{\mathbf{X}},\sigma_{\mathbf{X}}^{2})$.

After obtaining $\mathbf{T}$, we evaluate how much information $\mathbf{T}$
still contains about $\mathbf{X}$, which is defined as their mutual
information:

$\operatorname{I}(\mathbf{X};\mathbf{T})=\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T}|\mathbf{X})\|P(\mathbf{T})]],$
(3)

where $D_{KL}$ means Kullback–Leibler (KL) divergence,
$P(\mathbf{T}|\mathbf{X})$ and $P(\mathbf{T})$ represent their probability
distributions. While $P(\mathbf{T}|\mathbf{X})$ can be sampled empirically,
$P(\mathbf{T})$ has no analytical solution since it requires integrating over
the feature map $P(\mathbf{T})=\int
P(\mathbf{T}|\mathbf{X})P(\mathbf{X})\mathrm{d}\mathbf{X}$. As is standard, we
use the variational approximation
$Q(\mathbf{T})=\mathcal{N}(\mu_{\mathbf{X}},\sigma_{\mathbf{X}}^{2})$ to
substitute $P(\mathbf{T})$, assuming every dimension of $\mathbf{T}$ is
independent and normally distributed. Even though the independence assumption
does not hold in general, it only overestimates the mutual information, giving
a nice upper bound of mutual information between $\mathbf{X}$ and
$\mathbf{T}$:

$\displaystyle\operatorname{I}(\mathbf{X};\mathbf{T})$
$\displaystyle=\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T}|\mathbf{X})\|Q(\mathbf{T})]]$
$\displaystyle-D_{KL}[Q(\mathbf{T})\|P(\mathbf{T})]$ (4a)
$\displaystyle\operatorname{I}(\mathbf{X};\mathbf{T})$
$\displaystyle\leq\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T}|\mathbf{X})\|Q(\mathbf{T})]].$
(4b)

The complete derivation of Equation 4b is in Appendix A. Since we expect
$\operatorname{I}(\mathbf{X},\mathbf{T})$ to be small and mutual information
to be always nonnegative, the upper bound is a desired property.

Intuitively, the purpose of maximizing
$\operatorname{I}(\mathbf{Y};\mathbf{T})$ is to make accurate predictions.
Therefore, instead of directly maximizing
$\operatorname{I}(\mathbf{Y};\mathbf{T})$, we minimize the loss function for
the original task, e.g., the cross entropy $\mathcal{L}_{\text{CE}}$ for
classification problems after inserting the information bottleneck.

| IMDB | MNLI Matched | MNLI Mismatched | AG News | RTE
---|---|---|---|---|---
Random | 0.011 | 0.106 | 0.106 | 0.008 | 0.012
LIME | 0.038 | 0.244 | 0.260 | 0.033 | 0.014
IG | 0.090 | 0.226 | 0.233 | 0.036 | 0.043
IBA | 0.229 | 0.374 | 0.367 | 0.029 | 0.059
Table 1: Absolute probability drop for the target class after the top 11% most
important tokens removed. The larger the score, the more effective the method.

Combining the above two parts, our final loss function $\mathcal{L}$ is

$\mathcal{L}=\mathcal{L}_{\text{CE}}+\beta\cdot\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T|\mathbf{X}})\|Q(\mathbf{T})]].$
(5)

Note that we negate the sign for minimization. The $\beta$ hyperparameter
controls the relative importance between the two loss components. After the
optimization process, we obtain for every instance a compressed representation
$\mathbf{T}$.

We then calculate $D_{KL}[P(\mathbf{T}|\mathbf{X})\|Q(\mathbf{T})]$,
indicating how much information is still kept in $\mathbf{T}$ about
$\mathbf{X}$, which suggests the contribution of each token and feature. To
generate the attribution map, we sum over the feature–token axis, obtaining
the attribution score of each token.

Overall, we try to learn a compressed hidden representation $\mathbf{T}$ that
has just enough information about the input $\mathbf{X}$ to predict the output
$\mathbf{Y}$. This compression is done by adding noise, which removes the
least relevant feature-level information, with $\bm{\mu}$ controlling how much
to remove.

##  4 Experiments

Through experimentation, we analyze IBA both quantitatively and qualitatively
to understand how it interprets deep neural network across layers.

(a) IMDB

(b) MNLI Matched

(c) MNLI Mismatched

(d) AG News

(e) RTE

Figure 1: Degradation test results comparing IBA, IG, LIME, and random.

###  4.1 Experimental Setting

We compare our method on BERT with two other representative model-agnostic
instance-level methods—LIME Ribeiro et al. (2016), which explores
interpretable models for approximation and explanation, and integrated
gradients (IG) Sundararajan et al. (2017), a variation on computing the
gradients of the predicted output with respect to input features. For a simple
baseline, we also compare with “random,” whose attribution scores are assigned
randomly to tokens. On each dataset, we fine-tune BERT and apply these
interpretability techniques to the model. We note the test accuracy and
generate an attribution score for each token. Details of all parameters are
attached in Appendix D.

There is no consensus on how to evaluate interpretability methods
quantitatively Molnar (2019). LIME’s simulated evaluation leverages the ground
truth of already interpretable models like decision trees, but the ground
truth is unavailable for black-box models like neural networks. Therefore, we
follow Ancona et al. (2018) and Hooker et al. (2018) and carry out a
degradation test on IMDB Maas et al. (2011), AG News Gulli (2004), MNLI
Williams et al. (2018), and RTE Wang et al. (2018), covering sentiment
analysis, natural language inference, and text classification.

The degradation test has the following steps:

  1. 1.

Generate attribution scores $s$ for each interpretability method $f$:
$s=f(\mathcal{M},x,y)$, where $x$ is the test instance, $y$ is the target
label, and $\mathcal{M}$ is the model.

  2. 2.

Sort tokens by their attribution score in descending order.

  3. 3.

Remove top $k$ tokens to obtain $x^{\prime}$, the degraded instance; $k$ can
be preset.

  4. 4.

Test the target class probability $p(y|x^{\prime})$ with the original model on
the degraded instance.

  5. 5.

Repeat steps 3 and 4 until all tokens removed.

For the final visualization, we average all test instances at each degradation
step to compute $\bar{p}(y|x^{\prime})$. Then, we normalize the degradation
test result $\bar{p}(y|x^{\prime})$ to $[0,1]$ using the normalized
probability drop $\bar{d}=\frac{\bar{p}(y|x^{\prime})-m}{o-m}$, where $o$
means the original probability on the nondegraded instance, and $m$ means the
minimum of the fully degraded instance’s probability across all
interpretability models. In this way, the normalized probability drop
$\bar{d}$ will be independent of the original model quality and easily
comparable across models. Note that, for IBA, we perform the degradation test
on the original model, not the one with the inserted bottleneck. Thus, a large
$\beta$ does not directly cause the probability to drop. An effective
attribution map can find the most important tokens, which means
$\bar{p}(y|x^{\prime})$ after the degradation step will drop substantially.

###  4.2 Results and Analysis

Overall, the results show that our method better identifies the most important
tokens compared to other model-agnostic interpretability methods.

(a) IB after different layers.

(b) IB with different $\beta$.

Figure 2: Analysis of different layers and different $\beta$.

Figure 3: Illustrations from left to right are as follows: The before and
after comparison of inserting an information bottleneck after layer 6;
attribution for an IMDB example with the positive label; attribution for an
MNLI example with the contradiction label.

#### Quantitative Analysis.

Table 1 shows the absolute probability drop $\|\bar{p}(y|x)-o\|$ with the
first 11% of the important tokens removed. We further plot the normalized
probability drop after each percentage of the important tokens is removed, as
shown in Figure 1, indicating how much important information is lost for
prediction: the steeper the slope, the better the ability to capture important
tokens. For this experiment, we insert the information bottleneck after layer
9, and we see that removing important tokens that are identified by our method
deteriorates the probability the most on IMDB and MNLI Matched/Mismatched.

Of course, choosing the right layer to insert the information bottleneck is
crucial to the result. It also indicates which layer encodes the most
meaningful information for prediction. To investigate differences in inserting
information bottlenecks after different layers, we carry the degradation test
on 1000 random test samples across layers on IMDB, as shown in Figure 2(a)—see
Appendix B for all 12 layers. Insertion after layers 1, 8, and 9 generates
more meaningful attribution scores. At layer 1, the tokens remain distinct
(i.e., representations have not been aggregated), and it is likely that the
latent representation $\mathbf{T}$ is essentially capturing per-token
sentiment values. The big drop of $\bar{d}$ after layers 8 and 9, on the other
hand, is interesting. Recently, Xin et al. (2020) examined early exit
mechanisms in BERT and found that halting inference at layers 8 or 9 produces
results not much worse than full inference, which suggests that an abundance
of information is encoded in those layers.

Another important parameter is $\beta$, which controls the trade-off between
restricting the information flow and achieving greater accuracy. A smaller
$\beta$ allows more information through, and an extremely small $\beta$ has
the same effect of using $\mathbf{X}$ as the attribution map. As Figure 2(b)
shows, when $\beta\leq 1e-6$, the degradation curve is similar to the one
using $\mathbf{X}$ only. Appendix C shows the effects of different $\beta$ on
a specific example.

#### Qualitative Analysis.

The first plot in Figure 3 shows the before and after comparison of IB
insertion, with positive tokens highlighted. The second and third plots
visualize attribution maps for instances across layers. Consistent with our
quantitative analysis in Figure 2(a), these plots demonstrate that, for a
fully fine-tuned BERT, layers 8 and 9 seem to encode the most important
information for the prediction. For example, in the IMDB instance, liked and
intrigued have the highest attribution scores for the prediction of positive
sentiment across most layers—see layer 9 in particular. In the MNLI example,
never is mostly highlighted starting from layer 7 to predict “contradiction.”

##  5 Conclusion

In this paper, we adopt an information-bottleneck-based approach to analyze
attribution for transformers. Our method outperforms two widely used
attribution methods across four datasets in sentiment analysis, document
classification, and textual entailment. We also analyze the information across
layers both quantitatively and qualitatively.

## Acknowledgments

This research was supported in part by the Canada First Research Excellence
Fund and the Natural Sciences and Engineering Research Council (NSERC) of
Canada.

## References

  * van Aken et al. (2019) Betty van Aken, Benjamin Winter, Alexander Löser, and Felix A. Gers. 2019.  How does BERT answer questions? A layer-wise analysis of transformer representations.  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_.
  * Ancona et al. (2018) Marco Ancona, Enea Ceolini, Cengiz Öztireli, and Markus Gross. 2018.  Towards better understanding of gradient-based attribution methods for deep neural networks.  In _International Conference on Learning Representations_.
  * Bang et al. (2019) Seojin Bang, Pengtao Xie, Heewook Lee, Wei Wu, and Eric Xing. 2019.  Explaining a black-box using deep variational information bottleneck approach.  _arXiv:1902.06918_.
  * Clark et al. (2019) Kevin Clark, Urvashi Khandelwal, Omer Levy, and Christopher D. Manning. 2019.  What does BERT look at? An analysis of BERT’s attention.  In _Proceedings of the 2019 ACL Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP_.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_.
  * Fong and Vedaldi (2017) Ruth C. Fong and Andrea Vedaldi. 2017.  Interpretable explanations of black boxes by meaningful perturbation.  In _Proceedings of the IEEE International Conference on Computer Vision_.
  * Guan et al. (2019) Chaoyu Guan, Xiting Wang, Quanshi Zhang, Runjin Chen, Di He, and Xing Xie. 2019\.  Towards a deep and unified understanding of deep neural models in NLP.  In _International Conference on Machine Learning_.
  * Gulli (2004) Antonio Gulli. 2004.  _AGNews_.  http://groups.di.unipi.it/~gulli/AG_corpus_of_news_articles.html.
  * Hewitt and Manning (2019) John Hewitt and Christopher D. Manning. 2019.  A structural probe for finding syntax in word representations.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_.
  * Hooker et al. (2018) Sara Hooker, Dumitru Erhan, Pieter-Jan Kindermans, and Been Kim. 2018.  Evaluating feature importance estimates.  _arXiv:1806.10758_.
  * Kovaleva et al. (2019) Olga Kovaleva, Alexey Romanov, Anna Rogers, and Anna Rumshisky. 2019.  Revealing the dark secrets of BERT.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_.
  * Li et al. (2016) Jiwei Li, Xinlei Chen, Eduard Hovy, and Dan Jurafsky. 2016.  Visualizing and understanding neural models in NLP.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_.
  * Liu et al. (2019) Nelson F. Liu, Matt Gardner, Yonatan Belinkov, Matthew E Peters, and Noah A. Smith. 2019.  Linguistic knowledge and transferability of contextual representations.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_.
  * Lundberg and Lee (2017) Scott M. Lundberg and Su-In Lee. 2017.  A unified approach to interpreting model predictions.  In _Advances in Neural Information Processing systems_.
  * Maas et al. (2011) Andrew L. Maas, Raymond E. Daly, Peter T. Pham, Dan Huang, Andrew Y. Ng, and Christopher Potts. 2011.  Learning word vectors for sentiment analysis.  In _Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies_.
  * Molnar (2019) Christoph Molnar. 2019.  _Interpretable Machine Learning_.  https://christophm.github.io/interpretable-ml-book/.
  * Pimentel et al. (2020) Tiago Pimentel, Josef Valvoda, Rowan Hall Maudslay, Ran Zmigrod, Adina Williams, and Ryan Cotterell. 2020.  Information-theoretic probing for linguistic structure.  _arXiv:2004.03061_.
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016.  “Why should i trust you?” explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_.
  * Schulz et al. (2020) Karl Schulz, Leon Sixt, Federico Tombari, and Tim Landgraf. 2020.  Restricting the flow: Information bottlenecks for attribution.  In _International Conference on Learning Representations_.
  * Selvaraju et al. (2017) Ramprasaath R. Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra. 2017.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE international conference on computer vision_.
  * Smilkov et al. (2017) Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg. 2017.  SmoothGrad: removing noise by adding noise.  _arXiv:1706.03825_.
  * Sundararajan et al. (2017) Mukund Sundararajan, Ankur Taly, and Qiqi Yan. 2017.  Axiomatic attribution for deep networks.  In _International Conference on Machine Learning_.
  * Tenney et al. (2019) Ian Tenney, Patrick Xia, Berlin Chen, Alex Wang, Adam Poliak, R. Thomas McCoy, Najoung Kim, Benjamin Van Durme, Samuel Bowman, Dipanjan Das, and Ellie Pavlick. 2019.  What do you learn from context? Probing for sentence structure in contextualized word representations.  In _International Conference on Learning Representations_.
  * Tishby et al. (2000) Naftali Tishby, Fernando C. Pereira, and William Bialek. 2000.  The information bottleneck method.  _arXiv:physics/0004057_.
  * Wang et al. (2018) Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel Bowman. 2018.  GLUE: A multi-task benchmark and analysis platform for natural language understanding.  In _Proceedings of the 2018 EMNLP Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP_.
  * Williams et al. (2018) Adina Williams, Nikita Nangia, and Samuel Bowman. 2018.  A broad-coverage challenge corpus for sentence understanding through inference.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_.
  * Xin et al. (2020) Ji Xin, Raphael Tang, Jaejun Lee, Yaoliang Yu, and Jimmy Lin. 2020.  DeeBERT: Dynamic early exiting for accelerating BERT inference.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_.

##  Appendix A Proof of Variational Upper Bound

$\displaystyle\operatorname{I}(\mathbf{X};\mathbf{T})$
$\displaystyle=\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T}|\mathbf{X})\|P(\mathbf{T})]]$
$\displaystyle=\int_{\mathbf{X}}p(x)(\int_{\mathbf{T}}p(t|x)\log\frac{p(t|x)}{p(t)}dt)dx$
$\displaystyle=\int_{\mathbf{X}}\int_{\mathbf{T}}p(x,t)\log\frac{p(t|x)}{p(t)}\frac{q(t)}{q(t)}dtdx$
$\displaystyle=\int_{\mathbf{X}}\int_{\mathbf{T}}p(x,t)\log\frac{p(t|x)}{q(t)}dtdx$
$\displaystyle+\int_{\mathbf{X}}\int_{\mathbf{T}}p(x,t)\log\frac{q(t)}{p(t)}dtdx$
$\displaystyle=\int_{\mathbf{X}}\int_{\mathbf{T}}p(x,t)\log\frac{p(t|x)}{q(t)}dtdx$
$\displaystyle+\int_{\mathbf{T}}p(t)(\int_{\mathbf{X}}p(x|t)dx)\log\frac{q(t)}{p(t)}dt$
$\displaystyle=\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T}|\mathbf{X})\|Q(\mathbf{T})]]$
$\displaystyle-D_{KL}[Q(\mathbf{T})\|P(\mathbf{T})]$
$\displaystyle\leq\mathbb{E}_{\mathbf{X}}[D_{KL}[P(\mathbf{T}|\mathbf{X})\|Q(\mathbf{T})]]$

##  Appendix B Degradation Test across 12 Layers

Figure 4 shows the complete version of the degradation test across all 12
layers. In general, the earlier we insert the bottleneck, the larger the
probability drop is, except for layers 8 and 9, which are the only two layers
with steeper slopes than layer 1.

Figure 4: Degradation test results across all layers.

##  Appendix C Visualization of the Effects of $\beta$

Figure 5 shows the effects of different $\beta$ on a specific example. As we
can see, when $\beta$ is as small as $10^{-7}$, most information is allowed to
flow through the network and thus most parts are highlighted. In contrast,
when $\beta$ is larger, the representation is more restricted.

(a) $\beta=10^{-3}$

(b) $\beta=10^{-7}$

Figure 5: Comparison of BERT attribution maps with different values of
$\beta$.

##  Appendix D Detailed Parameters and Dataset Information

To keep as much information as possible at the beginning, $\bm{\mu}_{i}$
should be set close to 1,$\forall i$, in which case
$\mathbf{T}\approx\mathbf{X}$. So we initialize with $\alpha_{i}=5,\forall i$
and therefore $\bm{\mu}_{i}\approx 0.993$. In order to stabilize the result,
the input of the bottleneck ($\mathbf{X}$) is duplicated 10 times with
different noise added. We set the learning rate to 1 and the number of
training steps to 10. We use empirical estimation for $\beta\approx
10\times\frac{\mathcal{L}_{CE}}{\mathcal{L}_{IB}}$. For IMDB, MNLI
Matched/Mismatched, and AGNews, we insert the IB after layer 9 and $\beta$ is
set to $10^{-5}$. For RTE, we insert the IB after layer 10 and $\beta$ is set
to $10^{-4}$.

We carry out experiments on NVIDIA RTX 2080 Ti GPUs with 11GB VRAM running
PyTorch 1.4.0 and CUDA 10.0. A full technical description of our computing
environment is released alongside our codebase. For LIME, we set $N$, the
number of permuted samples drawn from the original dataset, to 100 as this
reaches the limitation of GPU memory. Similarly, the number of steps of
integrated gradients is set to 10 because it is more memory intensive. The
average time of running 25000 instances on the described GPU is about 10 hours
for IBA, 13 hours for LIME, and 2 hours for IG.

Dataset | Number of Dev/Test
---|---
IMDB | 25000
MNLI Matched | 9815
MNLI Mismatched | 9832
AG News | 7600
RTE | 277
Table 2: Dataset Details.

We use the test sets when the label is provided and use the dev sets
otherwise. See Table 2 for details. Note that “IMDB” refers to the sentiment
analysis dataset provided by Maas et al. (2011). “MNLI Matched” means that the
training set and the test set have the same set of genres while “MNLI
Mismatched” means that genres that appear in the test set don’t appear in the
training set. Detailed information of the MNLI dataset can be found in
Williams et al. (2018).


# Spectral Properties of Two Coupled Fibonacci Chains

Anouar Moustaj*1, Malte Röntgen*2,3, Christian V. Morfonios2, Peter
Schmelcher2,4, Cristiane Morais Smith1 1Institute of Theoretical Physics,
Utrecht University, Princetonplein 5, 3584CC Utrecht, The Netherlands
2Zentrum für Optische Quantentechnologien, Fachbereich Physik, Universität
Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany  3Laboratoire
d’Acoustique de l’Université du Mans, Unite Mixte de Recherche 6613, Centre
National de la Recherche Scientifique, Avenue O. Messiaen, F-72085 Le Mans
Cedex 9, France  4The Hamburg Centre for Ultrafast Imaging, Universität
Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany

###### Abstract

The Fibonacci chain, i.e., a tight-binding model where couplings and/or on-
site potentials can take only two different values distributed according to
the Fibonacci word, is a classical example of a one-dimensional quasicrystal.
With its many intriguing properties, such as a fractal eigenvalue spectrum,
the Fibonacci chain offers a rich platform to investigate many of the effects
that occur in three-dimensional quasicrystals. In this work, we study the
eigenvalues and eigenstates of two identical Fibonacci chains coupled to each
other in different ways. We find that this setup allows for a rich variety of
effects. Depending on the coupling scheme used, the resulting system (i)
possesses an eigenvalue spectrum featuring a richer hierarchical structure
compared to the spectrum of a single Fibonacci chain, (ii) shows a coexistence
of Bloch and critical eigenstates, or (iii) possesses a large number of
degenerate eigenstates, each of which is perfectly localized on only four
sites of the system. If additionally, the system is infinitely extended, the
macroscopic number of perfectly localized eigenstates induces a perfectly flat
quasi band. Especially the second case is interesting from an application
perspective, since eigenstates that are of Bloch or of critical character
feature largely different transport properties. At the same time, the proposed
setup allows for an experimental realization, e.g., with evanescently coupled
waveguides, electric circuits, or by patterning an anti-lattice with adatoms
on a metallic substrate.

  * •

Keywords: Quasicrystals, Flat Bands, Critical Eigenstates, Extented
Eigenstates.

##  1 Introduction

Aperiodic structures, in particular quasicrystals [1], have attracted the
attention of researchers for many decades [2, 3, 4, 5, 6]. Even the simplest
quasiperiodic, one-dimensional models exhibit a rich variety of behaviour,
ranging from critical eigenstate localization properties, to the appearance of
energy spectra that sometimes form a fractal set [7]. The tools used to
understand these behaviors include, but are not limited to, renormalization
procedures, multifractal analysis, or symmetry considerations [8, 9, 10]. One
of the peculiar features that these systems show is a hierarchical structure
of energy gaps [9]. This gives rise to spectral measures that are singular
continuous, such as the fractal set observed in the Fibonacci quasicrystal
[6]. This is a feature that was linked to all the critical localization
properties observed in quasicrystals as well as other systems that possess
some type of aperiodic order [7].

Although the behavior of aperiodic chains has been investigated extensively
and in great detail, comparatively little work has been dedicated to the case
where two or more chains are coupled to each other, forming an aperiodic
ladder [11, 12, 13, 14, 15]. In this work, we take a step into this realm by
analyzing a range of different coupling schemes between two identical one-
dimensional Fibonacci chains. Specifically, we study cases where the two
chains are directly coupled in a uniform, non-uniform, or quasiperiodic
manner, and evaluate the resulting spectral properties, namely the energy
eigenvalues and eigenstates. Additionally, we study a special case of an
indirect coupling, that is, two chains coupled to each other through some
intermediate sites. These cases are easily tractable, since they posses a
reflection symmetry which allows to block-diagonalize the Hamiltonian.

We find different spectra depending on the setup. In the case of uniform
coupling, the eigenvalue spectrum is identical to that of two uncoupled
Fibonacci chains, but with shifted energy eigenvalues. On the other hand, if
the two chains are coupled only through a single site, the spectrum consists
of two Fibonacci chains with an on-site defect. The structure of the
eigenvalue spectrum becomes more complex for the case of quasiperiodic
coupling, for which a richer hierarchical structure reveals itself through a
perturbative renormalization approach. If only the sites of one specific type
(A or B) are coupled to each other, we show that, for a specific value of the
interchain coupling, half of the eigenstates are critical, while the other
half are extended. Interestingly, these two different classes of eigenstates
possess different parity with respect to a corresponding reflection operation
and can thus be selectively excited by incoming waves of negative or positive
parity. This could be used to control the transport properties of this system.
Finally, we also realise couplings between the chains through some
intermediate sites. This leads to the appearance of flat bands in a
quasiperiodic lattice.

This paper is organized as follows. To be self-contained, we start by briefly
reviewing the properties of the Fibonacci chain in Section 2, followed by an
overview of the methods used to generate our results in Section 3. In Section
4, we analyze the simplest way of coupling two Fibonacci chains, namely a
uniform one. In Section 5, we investigate different cases of non-uniform
coupling. We start by connecting only sites of a specific type to each other
in Section 5.1. Then, in Section 5.2, we analyze the case where only two sites
are coupled. In Section 6, we couple the two chains in a quasiperiodic fashion
and analyze the resulting eigenvalue spectrum in terms of a renormalization
scheme. Finally, in Section 7, we consider the scenario where the two chains
are not directly coupled to each other, but through intermediate sites. Our
conclusions are presented in Section 8.

##  2 A single Fibonacci chain

We consider a general Fibonacci chain model, namely a nearest-neighbour tight-
binding chain with periodic boundary conditions, so that the chain effectively
becomes a ring. The on-site potential and hopping amplitudes are both
modulated by the Fibonacci sequence, and the corresponding Hamiltonian is
given by

$H_{F}=\sum_{i=1}^{N}v_{i}\ket{i}\bra{i}+\sum_{\langle
i,j\rangle}h_{i}\ket{i}\bra{j}$ (1)

where $\ket{i}$ denotes a basis state fully localized on the $i$-th site and
$\langle i,j\rangle$ denotes nearest-neighbors. The on-site potential $v_{i}$
and the hopping amplitude $h_{i}$ are both binary and follow the sequence of a
Fibonacci word $S_{n}$ of generation $n$. The latter is defined by the
recursion relation

$S_{n}=S_{n-1}S_{n-2},\ \ n\geq 2,$

which is a binary representation of the Fibonacci sequence. This inductive
recursion formula is expressed as a string concatenation instead of number
addition. Another typical representation is through the substitution rule
$A\to AB$ and $B\to A$. With the initial two words being $S_{0}=B$, $S_{1}=A$,
the first few words are thus

$\displaystyle S_{0}$ $\displaystyle=B$ $\displaystyle S_{1}$
$\displaystyle=A$ $\displaystyle S_{2}$ $\displaystyle=AB$ $\displaystyle
S_{3}$ $\displaystyle=ABA$ $\displaystyle S_{4}$ $\displaystyle=ABAAB$
$\displaystyle S_{5}$ $\displaystyle=ABAABABA$ $\displaystyle S_{6}$
$\displaystyle=ABAABABAABAAB\,.$

We note that the Fibonacci sequence is most commonly represented in numeric
form through the recursion formula

$F_{n}=F_{n-1}+F_{n-2},\ \ n\geq 2,$

where $F_{0}=F_{1}=1$ and $F_{n}$ is the $n^{\text{th}}$ Fibonacci number,
which is equal to the word length $|S_{n}|=F_{n}$. Another important property
of the sequence is that, in the thermodynamic limit, the ratio between the
amount of letters A and B is equal to the golden ratio $\phi$. This is more
appropriately expressed as [16]

$\lim_{n\to\infty}\frac{F_{n}}{F_{n-1}}=\frac{1+\sqrt{5}}{2}\equiv\phi\,.$

The Hamiltonian (1) can be studied in its general form or reduced to a purely
on-site or hopping Fibonacci model, where the parameters $h_{i}=h$ or
$v_{i}=v$ are uniform. These have been termed diagonal and off-diagonal models
in the literature, while Equation 1 is referred to as the mixed model. All
these models have been studied previously in their various forms [17, 18, 19,
20], but one important aspect is the equivalence between the on-site and
hopping models under a perturbative renormalization scheme [18]. This means
that we can uncover all essential features by studying the unmixed models. For
this reason, we will focus on the on-site model for our analysis. In this
case, the on-site potential $v_{i}$ of the $i$-th site is equal to either
$v_{A}$ or to $v_{B}$, depending on whether the $i$-th character of $S_{n}$ is
equal to $A$ or to $B$, and the hopping parameter is constant and equal to
$h$. We note that the number of sites in the chain is then equal to
$|S_{n}|=F_{n}$.

Figure 1: (a) Graphical representation of a Fibonacci chain with periodic
boundary conditions (see text for details). (b) Eigenstate map for such a
Fibonacci chain of $N=89$ sites. Each of the $N$ rows of the gray scale image
shows the value of $|\Psi^{(\nu)}(j)|^{1/2}$ (with each eigenstate
normalized), where $j$ is a site index, $\nu$ is an energy index, and with
each square pixel corresponding to one of the $N$ sites. The orange dots
denote the eigenenergy $\epsilon_{\nu}$ of the eigenstates. The eigenvalues
and eigenstates were obtained for $h=-1$, $v_{A}=0$ and $v_{B}=3$.

Figure 1 (a) shows a graphical representation of such a Fibonacci chain. Due
to the periodic boundary conditions, the depicted system constitutes the unit
cell of a periodic approximant of the quasicrystal (if $N\to\infty$, it
becomes a proper quasicrystal). This is a system which has long-range order,
without being periodic. Its eigenstates form a set of critical states, which
have atypical localization properties. This means that they are neither
extended nor localized [7]. The eigenvalues have a spectral measure that is
singular continuous and form a fractal set (a Cantor set of measure zero).
This feature is also observed in the eigenstates, with wavefunctions that show
multifractal properties [10]. In Figure 1 (b), the eigenstates of a Fibonacci
chain with $N=89$ are graphically depicted.

The properties of quasicrystals have been studied in a multitude of ways,
ranging from perturbative methods based on a renormalization formalism [18, 9]
to exact results using a transfer matrix approach [21] or a symmetry
perspective, which offers insights on the fragmentation of states in terms of
local spatial structures of the chain [22].

Before we continue, let us briefly comment on possible realizations of the
Hamiltonian eq. 1. We note that the following statements also hold for the
more complex models that we will present later in this manuscript. Since we
treat the matrix eigenvalue problem $H\ket{\Psi}=E\ket{\Psi}$, the setups
proposed in this work are well suited for a realization in various systems
modelled by such a problem, as long as the relevant matrix elements that
correspond to couplings and on-site potentials of our Hamiltonian can be
controlled. An immediate candidate, for instance, are systems of evanescently
coupled optical waveguides [23]. Here, each site corresponds to a monomodal
waveguide. By suitably tuning quantities such as waveguide shape, spacing, the
on-site potentials $v_{i}$ and couplings $h_{i,j}$ are tunable in a wide
range. We remark that the propagation of light in such arrays of evanescently
coupled photonic waveguides can be described by a discrete, time-dependent
Schrödinger equation, which would also allow to probe the dynamics of the
setups proposed in this work. Another realization that is directly feasible is
in terms of electric circuits, where the sites are nodes of the circuit, and
with on-site potentials and couplings being determined by how the nodes are
grounded and interconnected to each other [24, 25].

##  3 Methods

In this section, we will provide a comprehensive description of the methods
employed to calculate and illustrate the relevant quantities before delving
into the obtained results. We will begin by explaining the numerical methods
utilized in this study, which encompass the creation of all eigenstate maps.
Subsequently, we will outline the analytic methods employed to describe the
section concerning quasiperiodic coupling.

###  3.1 Numerical methods

Unless mentioned otherwise, all eigenvalues and eigenvectors were obtained by
numerically diagonalizing the corresponding Hamiltonian. In the case of two
coupled Fibonacci chains, we individually diagonalized the two Hamiltonians
$H_{\pm}$ (see section 4) and then constructed the total eigenstates by
symmetrizing/anti-symmetrizing these states. We note that this procedure
automatically provides an assignment of the eigenvalues of the total
Hamiltonian to negative/positive parity. The eigenstate maps (as in fig. 1)
and energy plots (as in fig. 3) were produced with Mathematica.

###  3.2 Analytical methods

#### Hierarchical splitting and renormalization

The main analytic tool used in this work is a perturbative renormalization
procedure based on the degenerate Brillouin-Wigner perturbation theory,
described in the works of Niu and Nori [9, 18].

Figure 2: Graphical representation of the renormalization procedure. The
topmost chain has a hierarchical distribution of hoppings, with $T_{3}\ll
T_{2}\ll T_{1}$. The blue sites in the middle chain represent the
(anti)bonding states of the unperturbed Hamiltonian of the topmost chain,
which have effective couplings $T_{2}/2$ and $T_{3}$/2. In its turn, the
lowermost chain represents the (anti)bonding states of the middle chain, with
its effective coupling given by $T_{3}/2^{2}$. Figure inspired from [9].

Consider a one-dimensional chain described by a tight-binding Hamiltonian that
incorporates a set of hierarchical hoppings, where the hopping strengths
follow the relation $T_{j}\ll T_{j-1}\ll\cdots T_{2}\ll T_{1}$ (as depicted in
the topmost chain of Figure 2). In this scenario, we can employ a perturbative
renormalization approach to calculate the energy levels of the topmost chain.

This perturbative renormalization approach involves treating the weaker
hoppings ($T_{j}$ with higher indices) as perturbations to the dominant
hopping term ($T_{1}$). This allows us to systematically incorporate the
effects of the weaker hoppings and compute the resulting energy levels of the
topmost chain. This is done by using the subset of unperturbed degenerate
eigenstates as the basis set for the perturbed Hamiltonian. We then apply
Brillouin-Wigner perturbation theory to calculate the effective couplings
between these states. In Figure 2, this process is visually represented by the
middle chain, where each blue site represents a (anti)bonding degenerate
eigenstate of the $T_{1}$ molecule. By considering these (anti)bonding states
as the basis, we can determine the effective couplings induced by the weaker
hoppings ($T_{j}$) and calculate the resulting energy levels of the topmost
chain. In principle, this process is repeated indefinitely, but in practice,
one takes a chain of finite size and imposes periodic boundary conditions to
perform the calculations. The energy levels form a cluster of $2^{j}$ values
around the base energy $E_{0}$

$E=E_{0}\pm T_{1}\pm\frac{T_{2}}{2}\pm\cdots\pm\frac{T_{j}}{2^{j-1}}\,.$ (2)

#### Brillouin Wigner perturbation theory

We shall now give a brief overview of Brillouin-Wigner perturbation theory.
Consider a Hamiltonian $H=H_{0}+H_{1}$, with $H_{1}$ acting as a perturbation
to $H_{0}$. Let $E_{0}$ be a degenerate energy eigenvalue, and $Q$ denote the
projection operator onto the corresponding degenerate subspace. We also denote
the complementary projection operator $P=\mathbbm{1}-Q$. The eigenvalue
equation $H\ket{\psi}=E\ket{\psi}$ can be rewritten such that

$P\ket{\psi}=P\frac{1}{E-H_{0}}H_{1}\ket{\psi}=P\frac{1}{E-H_{0}}PH_{1}\ket{\psi}$
(3)

where the last equality follows from $P^{2}=P$ and $PH_{0}=H_{0}P$. We now
write

$\displaystyle\ket{\psi}$ $\displaystyle=(Q+P)\ket{\psi}$
$\displaystyle=Q\ket{\psi}+P\frac{1}{E-H_{0}}PH_{1}\ket{\psi}$
$\displaystyle=Q\ket{\psi}+P\frac{1}{E-H_{0}}PH_{1}(Q+P)\ket{\psi}$
$\displaystyle\mathmakebox[width("{}={}")][c]{\vdots}$
$\displaystyle=\sum_{n=0}^{\infty}\left(P\frac{1}{E-H_{0}}H_{1}\right)^{n}Q\ket{\psi},$

where we have just consistently used Equation 3. With the above equation, we
can now easily obtain an effective Hamiltonian for the degenerate subspace
with eigenvalue $E_{0}$ through a left multiplication by $QH$,

$\begin{split}QH\ket{\psi}=EQ\ket{\psi}&=\left[QH_{0}Q+QH_{1}\sum_{n=0}^{\infty}\left(P\frac{1}{E-H_{0}}H_{1}\right)^{n}Q\right]Q\ket{\psi}\\\
&\equiv H_{\text{eff}}Q\ket{\psi}\end{split}$ (4)

Using this effective Hamiltonian, it will be possible to understand the
structure of the spectrum in the case of quasiperiodic coupling, which is
treated in Section 6.

##  4 Uniform coupling

Figure 3: (a) Two uniformly coupled Fibonacci chains. (b) The result of
writing the original system in terms of the symmetry-adapted basis Equations 7
and 8. In this basis, the setup decomposes into the two disconnected chains
$H_{+}$ and $H_{-}$. (c) A sketch of the shifting spectra for different
coupling strengths $h^{\prime}$ (see text for details). The eigenstates of
$H_{+}$ and $H_{-}$ correspond to eigenstates of the total Hamiltonian $H$
with positive and negative parity, respectively represented in (d). Here,
specific eigenstates of $H_{+}$ and $H_{-}$ are shown. At each site, the sign
of the eigenstate is represented by red/blue color, while the amplitude is
depicted by the radius of the circle.

Let us investigate the setups where two such chains—each consisting of $N$
sites and with periodic boundary conditions—are coupled to each other in
different ways. We will focus on the impact of these different coupling
schemes on the corresponding spectral properties.

The first and simplest scenario occurs when the two chains are uniformly
coupled to each other, as shown in Figure 3 (a). The setup is then described
by

$H=H_{I}+h^{\prime}\sum_{i=1}^{N}\Big{(}\ket{i_{u}}\bra{i_{l}}+\ket{i_{l}}\bra{i_{u}}\Big{)}$
(5)

with

$H_{I}=\sum_{x=u,l}\sum_{i=1}^{N}v_{i}\ket{i_{x}}\bra{i_{x}}+h\sum_{x=u,l}\sum_{\langle
i,j\rangle}\ket{i_{x}}\bra{j_{x}}$ (6)

where $\ket{i_{u}}$, $\ket{i_{l}}$ denote basis states fully localized on the
$i$-th site of the upper or lower chain, respectively. We note that such a
setup has been investigated in Ref.[26], though with a focus on the density of
states and not on the eigenvalues and eigenstates.

To understand the impact of such a uniform coupling, we employ the up/down
mirror symmetry of the setup. Due to this symmetry, the eigenstates have
definite parity under an exchange of the lower and upper chain. This fact can
be used to construct a symmetry-adapted basis $\mathcal{S}$, consisting of $N$
states of positive parity

$\ket{1_{u}}+\ket{1_{l}},\ldots{},\ket{N_{u}}+\ket{N_{l}}$ (7)

and $N$ states with negative parity

$\ket{1_{u}}-\ket{1_{l}},\ldots{},\ket{N_{u}}-\ket{N_{l}}\,.$ (8)

Written in this basis, $H^{\prime}=\mathcal{S}^{-1}H\mathcal{S}$ consists of
two isolated subsystems, $H_{+}$ and $H_{-}$. These two subsystems are shown
in Figure 3 (b). It can be seen that each of them is equal to an isolated
Fibonacci chain of $N$ sites, though with on-site potentials uniformly shifted
by an amount of plus or minus $h^{\prime}$. That is, $H_{\pm}=H_{I}\pm
h^{\prime}I$, with $I$ being the identity matrix.

Let us now explore the implications of the above, starting with the
eigenvalues of $H$. Since $H$ and $H^{\prime}$ are related by a similarity
transformation, the two Hamiltonians share the same eigenvalue spectrum.
Moreover, since $H^{\prime}$ consists of the two disconnected chains
$H_{\pm}$, the eigenvalue spectrum of $H^{\prime}$,
$\sigma(H^{\prime})=\sigma(H)$, is given by the combination of the eigenvalue
spectra of these chains 111To be precise, the spectrum of $H$ is the multiset
sum of the spectra of $H_{+}$ and $H_{-}$.. Now, because $H_{\pm}=H_{I}\pm
h^{\prime}I$, we see that the eigenvalue spectrum of $H_{+}$ ($H_{-}$) is that
of $H_{I}$ shifted upwards (downwards) by $h^{\prime}$, respectively. In other
words, the inter-chain coupling strength $h^{\prime}$ plays the role of an
“energy shift parameter” [see Figure 3 (c)]. Before we continue, we remark
that the eigenstates of $H$ can be simply constructed from those of $H_{+}$
and $H_{-}$ by symmetrizing or anti-symmetrizing these states; this is
demonstrated in Figure 3 (d).

##  5 Non-uniform coupling

Figure 4:  (a) Two Fibonacci chains that are coupled by connecting their
$A$-sites to each other through couplings with strength $h^{\prime}$. (b) The
decomposition into $H_{+}$ and $H_{-}$. (c) When choosing
$h^{\prime}=v_{B}-v_{A}$, $H_{+}$ becomes periodic, while $H_{-}$ remains a
Fibonacci chain. Thus, the positive parity eigenstates of $H$ (which
correspond to symmetrized eigenstates of $H_{+}$) are of Bloch-character,
while the negative parity eigenstates of $H$ (which correspond to anti-
symmetrized eigenstates of $H_{-}$) are critical states. (d) The eigenstate
map of the setup for this particular choice of $h^{\prime}$. Since the
eigenstates $\Psi^{(\nu)}$ of $H$ are either symmetric (orange dots; Bloch-
states) or anti-symmetric (blue dots; critical states), each of the $2N$ rows
in this eigenstate map only shows the amplitudes on one half of the sites of
the total system; the total eigenstate could be obtained through
symmetrization or anti-symmetrization. The eigenvalues and eigenstates were
obtained for $h=-1$, $v_{a}=0$ and $v_{B}=3$. The color map is the same as in
fig. 1, with black pixels corresponding to
$\max\limits_{\nu,j}\sqrt{|\Psi_{j}^{(\nu)}|}.$

Having understood the impact of coupling the two chains uniformly, we now
proceed to more complex scenarios. In all cases, we will maintain the
reflection symmetry between the two chains. Thus, we can still decompose the
total Hamiltonian into two smaller chains $H_{+}$ and $H_{-}$.

###  5.1 Coupling only $A$ or only $B$-sites

In the first case, we couple only the $A$ sites to each other, as depicted in
Figure 4(a). Repeating the same steps as above, we obtain
$H^{\prime}=\mathcal{S}^{-1}H\mathcal{S}=H_{+}\oplus H_{-}$, though now with

$H_{\pm}=\sum_{i=1}^{N}v^{\pm}_{i}\ket{i}\bra{i}+h\sum_{\langle
i,j\rangle}\ket{i}\bra{j}$ (9)

where $v^{\pm}_{A}=v_{A}\pm h^{\prime}$, while $v_{B}^{\pm}=v_{B}$ is
unchanged [see Figure 4(b)]. In a completely analogous manner, coupling only
the $B$-sites to each other will result in an energy shift of the on-site
potentials of the $B$-sites only.

A particularly interesting case occurs when coupling only the $A$-sites and
setting $h^{\prime}=v_{B}-v_{A}$. For this special choice of $h^{\prime}$,
$H_{+}$ becomes a uniform chain with zero on-site potential. However, $H_{-}$
is still a Fibonacci chain. Now, since the eigenstates of $H_{\pm}$ correspond
to positive/negative parity eigenstates of the full chain, the system features
an interesting combination of traits: while the positive parity eigenstates
are extended, the negative parity eigenstates are critical. For a quantum
system, this means that the phase diagram of such a Hamiltonian features a
special point $h^{\prime}=v_{B}-v_{A}$ , at which the system’s behavior is
highly dependent on the energies of one-particle excitations. As depicted in
Figure 4(a), depending on the energy of these parity eigenstates, they will
either form Bloch waves (orange energy levels), or critical states (blue
energy levels). It is interesting to see that a specific point in parameter
space shows a mixture of singular continuous and absolutely continuous
spectra. This provides a platform where both properties of extended and
critical states can be exploited by tuning the Fermi level. In the critical
regime, for example, (thermal) conductivities are in general very low (in some
cases, they are even lower than for conventional insulators) [27, 28]. On the
other hand, the fully extended regime provides the possibility to have phases
with high (electrical) conductivities.

Finally, let us note that the possibility of Bloch-states in coupled aperiodic
setups has also been observed in Refs. [12, 13], in which more complicated
coupling schemes have been used.

Figure 5: (a) Two Fibonacci chains coupled by connecting a single site of each
chain to each other. (b) The decomposition into $H_{+}$ and $H_{-}$. (c)
Eigenstate map for a setup of two Fibonacci chains that are coupled together
at the $10$-th site. As a result of this _defect-coupling_ (see text for
details), several eigenstates lie in gaps between quasi bands. Orange dots
correspond to eigenvalues of $H_{+}$ and blue dots to those of $H_{-}$. The
eigenvalues and eigenstates were obtained with $h=-1$, $v_{A}=0$ and
$v_{B}=3$. The color map is the same as in fig. 1, with black pixels
corresponding to $\max\limits_{\nu,j}\sqrt{|\Psi_{j}^{(\nu)}|}.$

###  5.2 Defect-coupling

Figure 6: (a) An excerpt of two Fibonacci chains coupled to each other in a
quasiperiodic manner, that is, by coupling the $B$-sites occuring in the
pattern $AABAA$ together. (b) The result of the symmetry-induced decomposition
into $H_{+}$ and $H_{-}$. (c) Since we know that the symmetry-induced
decomposition results in a different on-site potential for $H_{\pm}$, we
already assign a different $v_{B}^{\prime}$ to the coupled sites, which for
ease of analysis is set to $v_{B}^{\prime}=v_{B}+h$, resulting in the case we
study in more detail in Section 6.1. In this simplified version, $H_{-}$ is
the standard Fibonacci chain, so we only need to analyze $H_{+}$.

Another possibility is the selective coupling of only a small subset of sites.
In the extreme case, this subset consists of one site in each chain [see
Figure 5 (a)]. The result of such a coupling can be easily deduced. Using the
symmetry adapted basis, we see that $H_{\pm}$ will both be simple Fibonacci
chains, with added impurities at the sites that are coupled [see Figure 5
(b)]. Such systems have been analyzed previously, and it was found that a
single weak impurity is sufficient to render the spectrum unstable and reduce
its fractal dimension, leading to a loss of criticality in all states [29]. On
top of that, it was also recently found that this does not affect all states
equally. Using Niu’s renormalization procedure [9, 30], it was shown that the
degree of criticality loss is dependent on the renormalization path of the
site at which the impurity is placed [31]. This means that the location of the
impurity impacts which states of the unperturbed Fibonacci chain are the most
affected.

Alternatively, the impact of a single-site defect can be analyzed in a
framework of local resonators [22]. These local resonators form building
blocks of the whole chain. In an unperturbed (without impurity) chain, each
eigenstate is approximately symmetric with respect to a local parity operator.
The sites with the highest amplitudes will be the ones corresponding to the
resonator structures, which are symmetric under local parity. This approximate
symmetry depends on the contrast $c=|h|/|v_{A}-v_{B}|$, and is exact in the
limit $c\to 0$. By placing an impurity on a particular site, one creates a new
local resonator structure in the chain, and as such, new localization
properties arise, yielding states with amplitude distributions that are
radically different from the rest of the eigenstates. This can be seen in
Figure 5 (c), where the in-gap states have very strong localization, marked by
the darker patches around a particular region of the chain (see the topmost
level, for example, where there is a very dark patch around site $j=10$)
around this new local resonator block [22].

##  6 Quasiperiodic coupling

Yet another alternative way of coupling the two chains is in a quasiperiodic
manner. Out of the many possibilities, here we illustrate an immediate and
interesting one: We couple only a subset of $B$-sites to each other; namely,
those appearing in the pattern $AABAA$, as shown in Figure 6(a), where the
darker $B$ site sits in between two $A$ sites on each sides. We further set
the coupling between the chains to $h$ and choose the coupled $B$ site energy
to $v_{B}^{\prime}=v_{B}+h$, such that $H_{-}$ becomes a regular Fibonacci
chain, with on-site energies $v_{A}$ and $v_{B}$. On the other hand, $H_{+}$
now features new on-site energies: $v_{B}+2h$, distributed in a quasiperiodic
manner. For this choice of coupling, there are $F_{N-5}$ such new sites (which
is the number of $AABAA$ blocks in a chain of length $F_{N}$).

In analyzing the effective chain using a renormalization approach, we make the
assumptions that $h<0$ and $v_{A}<v_{B}$. These assumptions are made to
enforce the trifurcating structure necessary for the Fibonacci chain’s
renormalization group flow.

###  6.1 Decimation Procedure.

Figure 7: Decimation procedure for the effective chain described by $H_{+}$.
For each of the three degenerate levels, we define a decimated chain with
effective new couplings. This allows us to track how the different levels
structure themselves. In this example, we focus on the chain resulting from
the $v_{A}$ cluster. The couplings are first renormalized to three values
$|t_{s}|>|t_{b}|>|t_{w}|$. The amount of lines connecting the unperturberd
eigenstates are representative of the coupling strengths. At the next order of
perturbation theory, the unperturbed degenerate eigenstates are slightly more
complicated, and correspond to four-atom molecules, dimers, and isolated
sites. These then correspond to three chains with cyan, blue and green sites,
respectively. The next corrections correspond to the Fibonacci case, where the
levels trifurcate at each step. The hopping parameters $T_{A}$ and $T_{B}$ are
calculated in A.

We can apply the usual decimation procedure known for hierarchical chains
[18]. In this case, the novelty lies in the chain with three different on-site
energies, as shown for example in Figure 7. There are three renormalized
chains that result from the first decimation step. We can analyze them
separately to see how each energy level splits into different branches. A
sketch of the branching structure is shown in Figure 8, where one sees that
there are three main clusters. The $v_{C}$ cluster follows the Fibonacci
trifurcarting structure from the start. The $v_{B}$ cluster splits into six
levels, each of which starts trifurcating according to the Fibonacci
structure. Finally, the $v_{A}$ cluster splits into seven levels that also
trifurcate afterwards. In the thermodynamic limit, these levels keep on
trifurcating indefinitely, leading to a spectrum that is a Cantor set and
known to be singular continuous [32]. In A, we provide a comprehensive
explanation of the step-by-step renormalization process involved in handling
the binary hoppings $T_{A}$ and $T_{B}$ for each cluster. In the next section,
we show that these analytic tools lead to a good approximation and
understanding of the structure of the spectrum.

Figure 8: Energy splitting structure after each renormalization step. Energy
scales are greatly distorted and do not represent actual gap sizes. Note that
$v_{C}=v_{B}+2h$.

###  6.2 Effective Couplings and Energy Corrections

In order to calculate the effective couplings and energy corrections at each
splitting, we use the Brillouin-Wigner perturbation theory. For each cluster,
we use the effective Hamiltonian Equation 4

$H_{\rm{eff}}=QH_{0}Q+QH_{1}\sum_{n=0}^{n}\left(P\frac{1}{E-H_{0}}H_{1}\right)^{n}Q,$

where $H_{1}$ is the perturbation corresponding to the weakest coupling, $Q$
is the projector onto the eigenspace of $H_{0}$ corresponding to the cluster
of interest, while $P$ projects out of it. At first, $H_{0}$ is just composed
of on-site energies and $H_{1}$ of the coupling $h$ between the isolated
sites. This results in the three degenerate levels that form the main
clusters. In the next order of perturbation theory, $H_{0}$ denotes the
Hamiltonian of the corresponding renormalized chain, with the weakest coupling
turned off. $H_{1}$ is then the perturbation with either $t_{w}$ or $t_{b}$
turned on, depending on which subchain we are dealing with (i.e. depending on
whether $t_{b}>t_{w}$ or vice-versa). With this effective Hamiltonian for each
of the clusters, the renormalized couplings can be found by calculating
$\bra{E_{j}}H_{\rm{eff}}\ket{E_{j+1}}$, where $\ket{E_{j}}$ are the zeroth
order local eigenstates, which are either one atomic site, a two or a four-
atom molecule, depending on the situation. We now first state the results
obtained from the first-order corrections, which we shall subsequently prove.
We find the first three renormalized couplings to be given by

$\displaystyle\left(\frac{c}{1-2c},\ 1,\ c\right)h,\ \ \ \text{(Cluster A)}$
$\displaystyle\left(-\frac{c^{4}}{2},\ -c,\ c^{2}\right)h,\ \ \ \text{(Cluster
B)}$ (10) $\displaystyle\left(0,\
-\frac{1}{4}\left[\frac{c}{1-2c}\right]^{5},\
\frac{1}{16}\left[\frac{c}{1-2c}\right]^{8}\right)h,\ \ \ \text{(Cluster C)}$

where we have defined the contrast parameter $c\equiv|h/(v_{A}-v_{B})|$, which
controls how well the perturbation theory behaves. We point out that the
physics described by the theory is consistent for

$c<\frac{1}{4^{\frac{1}{3}}(1+4^{\frac{1}{6}})}\approx 0.27875,$ (11)

after which it fails to deliver reasonable results. This result will also be
proven in the following paragraphs.

#### First order renormalized hoppings

Figure 9: Decimation procedure. The topmost chain results from the
equipartition theorem applied to the coupled Fibonacci chains we started with.
It corresponds to the Hamiltonian $H_{+}$ with red sites of energies
$v_{C}=v_{B}+2h$. The white and grey sites have energies $v_{A}$ and $v_{B}$
respectively. (a) The renormalized chain that was shown as an example in
Section 6, with $|t_{s}|>|t_{b}|>|t_{w}|$. (b) The chain corresponding to the
renormalization of the grey sites. In this case, we have
$|t_{s}|>|t_{w}|>|t_{b}|$. (c) The chain corresponding to the renormalization
of the red sites, with $|t_{s}|>|t_{w}|$. This chain is already a proper
hopping Fibonacci chain.

Figure 9 shows that one can form three effective chains (a), (b) and (c),
corresponding to the clusters $v_{A}$, $v_{B}$ and $v_{C}=v_{B}+2h$
respectively. Before deriving our results, we recall that we are working in
the regime of parameters where $v_{A}<v_{B}$ and $h<0$. This is done to obtain
a hierarchy of hopping strengths that will result in the infinitely
trifurcating structure after (at most) the second renormalization step. The
$v_{C}$ cluster is a regular hopping Fibonacci chain. In fact, if one starts
from a chain of generation $N$, with $F_{N}$ sites, one ends up with one of
generation $N-5$, with $F_{N-5}$ sites. By calculating the matrix element
$\bra{E_{j}}H_{\rm{eff}}\ket{E_{j+1}}$, where $H_{\rm{eff}}$ is given by
Equation 4, between the nearest neighbouring eigenstates of the unperturbed
Hamiltonian $H_{0}$ in the $E=v_{C}$ subspace, we find that this chain has
effective renormalized hoppings, to nearest order in the contrast
$c=|h|/|v_{A}-v_{B}|$, given by

$\begin{split}t_{s}&=\frac{h^{6}}{4(v_{B}-v_{A}+2h)^{5}}=-\frac{1}{4}\left[\frac{c}{1-2c}\right]^{5}h,\\\
t_{w}&=\frac{1}{16}\frac{h^{9}}{(v_{B}-v_{A}+2h)^{8}}=\frac{1}{16}\left[\frac{c}{1-2c}\right]^{8}h.\end{split}$
(12)

Next, we consider the $v_{B}$ chain, which contains six clusters. Two clusters
correspond to the two-atom molecule bonded by $t_{s}$, while the other four
correspond to the four-atom molecule (composed of two two-atom molecules
bonded by $|t_{w}|>|t_{b}|$). These clusters will further start to trifurcate
because they end up having the structure of a hopping Fibonacci chain after
decimation. The renormalized hoppings of this chain are given by

$\begin{split}t_{b}&=\frac{-h^{5}}{2(v_{B}-v_{A})^{4}}=-\frac{c^{4}}{2}h,\\\
t_{s}&=\frac{h^{2}}{v_{B}-v_{A}}=-ch,\\\
t_{w}&=\frac{h^{3}}{(v_{B}-v_{A})^{2}}=c^{2}h.\end{split}$ (13)

Finally, we consider the $v_{A}$ cluster. It has seven clusters. The six
clusters come from the same structures as the $v_{B}$ chain, i.e two- and
four-atoms molecules, and the additional cluster comes from isolated sites.
Once again, each one of these clusters will start trifurcating at the next
steps of decimation, and one retrieves the hopping Fibonacci chain structure.
The renormalized hoppings of this chain are

$\begin{split}t_{b}&=\frac{h^{2}}{v_{A}-v_{B}-2h}=\frac{c}{1-2c}h,\\\
t_{s}&=h,\\\ t_{w}&=\frac{h^{2}}{v_{A}-v_{B}}=ch.\end{split}$ (14)

  * $E_{0}$ | $E_{1}$
---|---
$v_{A}$ | $\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}+\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}-\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$t_{s}$
0
$-t_{s}$
$\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}+\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$-\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}-\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$v_{B}$ | $\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}+\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
-$\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}-\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$t_{s}$
$-t_{s}$
$-\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}+\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$-\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}-\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}}$
$v_{B}+2h$ | $t_{s}$
0
-$t_{s}$

Table 1: The first order correction to the three clusters $v_{A}$, $v_{B}$ and
$v_{B}+2h$. Any higher order correction will further split each of these
energies into three, with corrections given by $\pm T_{A}$ and 0 (see A for
the values of $T_{A}$ in each case). Each of the three levels again will split
into three and so on, with the well known spectrum structure of the hopping
Fibonacci chain.

#### Regime of validity of the perturbation theory

Let us now study the regime of validity which results in the hopping hierarchy
discussed previously. Starting with the $v_{B}$ subspace, which is the
simplest to deal with, we see from Equation 13 that the hierarchy
$|t_{s}|>|t_{w}|>|t_{b}|$ will always hold for $c\in(0,1)$. This immediately
results in the trifurcating structure after one iteration of the
renormalization procedure. The $v_{A}$ subspace, on the other hand, restricts
the range of $c$ further. In order to have $|t_{s}|>|t_{b}|>|t_{w}|$, we must
impose

$1>\left|\frac{c}{1-2c}\right|>c$

which leads to $0<c<1/3$. Finally, the $v_{C}$ subspace will give us the final
restriction to impose the hierarchy leading to a trifurcating structure. This
means we want $|t_{s}|>|t_{w}|$, leading to

$\frac{1}{4}\left|\frac{c}{1-2c}\right|^{5}>\left|\frac{c}{1-2c}\right|^{8}$

Solving for $c>0$ yields the final restriction

$c<\frac{1}{4^{\frac{1}{3}}(1+4^{\frac{1}{6}})}.$ (15)

#### Energy corrections

We also determined the first-order energy corrections using the same
perturbation theory, and found the spectrum shown in Figure 10. The energy
levels calculated using perturbation theory (in blue) and those from numerical
direct diagonalization (in red) have the same structure, with some
discrepancies that disappear as the contrast $c\to 0$. This is illustrated in
Figure 10 (a) and (b) with $c=1/4$ and $c=1/8$, respectively. The structure of
the spectrum is well approximated by the theory, even for a small chain size
of 21 sites. The exact expressions for the energy levels can be found in Table
1.

Figure 10: Comparison between the perturbative calculations (blue) and the
direct numerical calculation (red) of the chain containing two “impurities” of
strength $V_{B}+2h$, at the center of the resonator structure $AABAA$. The
chain’s length is 21 sites. The structure is the same and the discrepancy gets
smaller as the contrast decreases. (a) Contrast $c=1/4$, which is close to the
limits of applicability of the perturbation theory. (b) Contrast $c=1/8$. The
energy is shown in arbitrary units. Figure 11: Eigenstate map for a setup
consisting of two quasiperiodically coupled Fibonacci chains, each of length
$N=55$ sites. Orange dots correspond to eigenvalues of $H_{+}$ and blue dots
correspond to eigenvalues of $H_{-}$. The color map is the same as in fig. 1,
with black pixels corresponding to
$\max\limits_{\nu,j}\sqrt{|\Psi_{j}^{(\nu)}|}.$

Now, one only needs to add the spectrum of the regular chain to this
“modified” chain to obtain the full spectrum of the two connected chains, as
shown in Figure 11.

##  7 Coupling through intermediate sites

In the following, we consider a case where the two Fibonacci chains are not
coupled directly, but instead through intermediate sites. Before we start, let
us first investigate a setup where two _periodic_ chains are coupled
indirectly, as shown in Figure 12 (a). In Figure 12 (b), we depict the band
structure of this so-called “one-dimensional Lieb lattice” [33]. What makes
this lattice interesting is that, among four dispersive bands, it also
features one completely flat band fulfilling $E(k)=\text{constant}$ for all
$k$. This defining feature of flat bands renders them dispersionless; they
suppress wave transport [34]. On the other hand, the density of states in a
flat band diverges, so that any disorder or non-linear effects may
qualitatively change the transport properties [34, 35]. Moreover, in lattices
where the single-particle Hamiltonian features a flat band, even arbitrarily
weak interactions may act non-perturbatively. This can lead to boson pair
formation [36, 37, 38] or other interesting phases, such as the Haldane
insulator [39] and Wigner crystals [40].

Figure 12:  (a) The one-dimensional Lieb-lattice and its spectrum (b). (c) the
quasiperiodic Lieb lattice and its spectrum (d). For the quasiperiodic Lieb
lattice, the spectrum was computed for a finite setup with a total number of
$88$ sites ($34$ in the upper and lower chain and $20$ $B$-sites in the
center). For all figures, $v_{A}=0$ and $v_{B}=h=1$.

In two dimensions, a classical example for a system with flat bands is the
Lieb lattice, whose structure is very similar to the system depicted in Figure
12 (a). The two-dimensional Lieb lattice has been realized in a number of
different experimental platforms, such as tailored atomic structures on
substrates [41, 42], evanescently coupled waveguide arrays [43, 44], terahertz
spoof plasmons [45], or cold atom setups [46, 47]. Interestingly, the CuO2
planes in high-temperature cuprate superconductors possess a Lieb lattice
structure, and it has been conjectured that flat bands might play a role in
their high critical temperature [34, 48, 49, 50, 51, 52].

Flat bands are also interesting from another perspective, since they are
tightly connected to the emergence of a special kind of eigenstates, the so-
called “compact localized state” (CLS) [53, 54, 55]. For the one-dimensional
Lieb lattice, a CLS is shown in Figure 12 (a), consisting of an excitation of
only four $B$ sites. The CLS thus “lives” only on a single plaquette (marked
by a dotted rectangle) and strictly vanishes outside it. In other words, it is
perfectly localized on a very small part of the setup.

The defining feature of CLSs — namely, their perfect localization — renders
these states very robust against perturbations: Since they vanish _exactly_
outside their localization domain $\mathcal{D}$, they are not affected by any
changes to the system outside $\mathcal{D}$. Due to this property, CLSs are
ideal candidates for storing information [56, 57]. The perfect localization of
CLSs might further be interesting in the context of photonic waveguide arrays,
where it allows for diffraction-free transmission of information in the form
of CLSs [58, 43].

After the above considerations, let us now couple two aperiodic Fibonacci
chains in an indirect manner. Out of the many possibilities, here we choose
one of the simplest, yet quite interesting setup, depicted in Figure 12 (c).
Each $A$-site of the upper chain is coupled by an intermediate $B$-site to its
counterpart $A$-site on the lower chain. This setup also features a
macroscopic number of CLSs, one of which is depicted in Figure 12 (c).

When comparing the one-dimensional version of the Lieb lattice [59, 60]—shown
in Figure 12 (a)—to the coupled Fibonacci chains depicted in Figure 12 (c),
one sees that they are rather similar. Due to this high similarity, we call
the coupled Fibonacci setup of Figure 12 (c) the “quasiperiodic Lieb lattice”.
The similarity in the structure of the two lattices is also visible in their
eigenvalue spectra, since both lattices feature a flat (quasi) band. The
emergence of this flat band can easily be understood, since (for an infinite
setup) there is an infinite number of plaquettes, and thus an infinite number
of degenerate CLSs. In both cases, the flat bands emerge due to destructive
interference [34].

##  8 Conclusions

In this work, we analyzed various ways of coupling two identical Fibonacci
chains to each other. The main tool that helped us to identify the features of
these systems is a symmetry adapted block-diagonalization of the Hamiltonian
$H$ into $H_{+}\oplus H_{-}$. Once we find the eigenstates of $H_{\pm}$, we
can symmetrize/anti-symmetrize them, respectively, to obtain the eigenstates
of the original Hamiltonian $H$. In addition, the eigenvalue spectrum of the
system is simply given by a multiset sum of the eigenvalue spectra of the
blocks $H_{\pm}$.

After briefly introducing the individual Fibonacci chain, we started by
exploring the effects of uniformly coupling two identical chains. We then
found that the resulting eigenvalue spectrum is just a sum of two shifted
Fibonacci spectra, which renders the behaviour of a particle in such a system
identical to that in a conventional Fibonacci chain. We have subsequently
explored the case of a nonuniform coupling, where we coupled only $A$ (or only
$B$) sites. An interesting scenario occurs when the interchain coupling is
$h^{\prime}=v_{B}-v_{A}$, since $H_{+}$ then becomes a periodic chain, for
which the eigenstates are Bloch waves. On the other hand, $H_{-}$ is still a
Fibonacci chain, such that the complete spectrum offers a mixture of critical
and fully extended eigenstates, identifiable by the parity of the
corresponding wavefunctions. The next type of coupling that we have analyzed
is between a small subset of sites, i.e. a so called defect coupling, leading
to block Hamiltonians $H_{\pm}$ which are Fibonacci chains with on-site
defects. This has been followed by another interesting and more complicated
case of quasiperiodically coupled chains. The resulting block Hamiltonians
could be thought to be like several coupled defects, but we have shown,
through a perturbative renormalization analysis, that all states in this chain
belong to the same family of critical states. Finally, we explored two
Fibonacci chains coupled to each other in the same manner as the nonuniform
coupling of Section 5, but with an intermediate site in between. This offers
the possibility of having a set of CLSs, leading to a flat bland in the energy
spectrum. Overall, the demonstrated emergence of CLSs in the quasiperiodic
Lieb lattice represents an interesting addition to the existing literature on
these phenomena in quasiperiodic setups (see, for instance, [61, 62, 63]). Two
immediate tasks for the near future would be to analyze the quasiperiodic Lieb
lattice in the context of interacting electrons, or to investigate the impact
of (correlated) disorder on the transport properties.

##  9 Acknowledgments

This publication is part of the project TOPCORE with project number
OCENW.GROOT.2019.048, which is financed by the Dutch Research Council (NWO).
This work (P.S.) is supported by the Cluster of Excellence ’Advanced Imaging
of Matter’ of the Deutsche Forschungsgemeinschaft (DFG) - EXC 2056 - project
ID 390715994.

* A.M. and M.R. contributed equally to this work.

## References

  * [1] Shechtman D, Blech I, Gratias D and Cahn J W 1984 Phys. Rev. Lett. 53(20) 1951–1953 URL https://link.aps.org/doi/10.1103/PhysRevLett.53.1951
  * [2] Janssen T 1986 Act. Crystall. Sect. A 42 261–271
  * [3] Berger C, Grenet T, Lindqvist P, Lanco P, Grieco J, Fourcaudot G and Cyrot-Lackmann F 1993 Solid State Communications 87 977 – 979 ISSN 0038-1098 URL http://www.sciencedirect.com/science/article/pii/003810989390543V
  * [4] Vieira A P 2005 Physical Review Letters 94 077201 ISSN 0031-9007 URL https://link.aps.org/doi/10.1103/PhysRevLett.94.077201
  * [5] Tanese D, Gurevich E, Baboux F, Jacqmin T, Lemaître A, Galopin E, Sagnes I, Amo A, Bloch J and Akkermans E 2014 Physical Review Letters 112 146404 ISSN 0031-9007 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.146404
  * [6] Jagannathan A 2021 Rev. Mod. Phys. 93 045001
  * [7] Maciá E 2005 Reports on Progress in Physics 69 397–441 URL https://doi.org/10.1088%2F0034-4885%2F69%2F2%2Fr03
  * [8] de Boissieu M 2019 Act. Crystall. Sect. A 75(Pt 2) 273–280.
  * [9] Niu Q and Nori F 1986 Phys. Rev. Lett. 57 2057–2060 ISSN 0031-9007
  * [10] Macé N, Jagannathan A and Piéchon F 2016 Phys. Rev. B 93(20) 205153 URL https://link.aps.org/doi/10.1103/PhysRevB.93.205153
  * [11] Moreira D A, Albuquerque E L and Bezerra C G 2006 The European Physical Journal B 54 393–398 ISSN 1434-6028
  * [12] Pal B and Chakrabarti A 2014 Physica E: Low-dimensional Systems and Nanostructures 60 188–195 ISSN 1386-9477
  * [13] Mukherjee A, Nandy A and Chakrabarti A 2017 Eur. Phys. J. B 90 52 ISSN 1434-6036
  * [14] Saha M and Maiti S K 2019 Journal of Physics D: Applied Physics 52 465304 ISSN 0022-3727 URL https://iopscience.iop.org/article/10.1088/1361-6463/ab3a0e
  * [15] Roy S, Maiti S K, Pérez L M, Silva J H O and Laroze D 2022 Materials 15 597 ISSN 1996-1944
  * [16] Kilic E 2008 European Journal of Combinatorics 29 701–711 ISSN 0195-6698 URL https://www.sciencedirect.com/science/article/pii/S0195669807000595
  * [17] Sire C and Mosseri R 1990 Journal de Physique 51 1569 – 1583
  * [18] Niu Q and Nori F 1991 Phys. Rev. B 42 10329–10341
  * [19] Kohmoto M and Banavar J R 1986 Physical Review B 34 563–566 ISSN 0163-1829
  * [20] Kohmoto M, Sutherland B and Tang C 1987 Phys. Rev. B 35 1020–1033
  * [21] Kohmoto M, Kadanoff L P and Tang C 1983 Phys. Rev. Lett. 50 1870–1872
  * [22] Röntgen M, Morfonios C V, Wang R, Dal Negro L and Schmelcher P 2019 Phys. Rev. B 99 214201
  * [23] Szameit A, Dreisow F and Nolte S 2012 Discrete optics in femtosecond laser written waveguide arrays Femtosecond Laser Micromachining: Photonic and Microfluidic Devices in Transparent Materials (Topics in Applied Physics no 123) (Springer, Berlin, Heidelberg) pp 351–388 ISBN 978-3-642-23365-4
  * [24] Lee C H, Imhof S, Berger C, Bayer F, Brehm J, Molenkamp L W, Kiessling T and Thomale R 2018 Commun Phys 1 1–9 ISSN 2399-3650
  * [25] Dong J, Juričić V and Roy B 2021 Phys. Rev. Research 3 023056
  * [26] Lazo E 2000 Multifractal Behavior of a Fibonacci Crystal Built over p Coupled Chains Instabilities and Nonequilibrium Structures VI Nonlinear Phenomena and Complex Systems ed Tirapegui E, Martínez J and Tiemann R (Dordrecht: Springer Netherlands) pp 387–392 ISBN 978-94-011-4247-2
  * [27] Archambault P and Janot C 1997 MRS Bulletin 22 48–53
  * [28] Janot, C 1996 Europhys. News 27 60–64 URL https://doi.org/10.1051/epn/19962702060
  * [29] Naumis G G and Aragón J L 1996 Phys. Rev. B 54(21) 15079–15085 URL https://link.aps.org/doi/10.1103/PhysRevB.54.15079
  * [30] Niu Q and Nori F 1990 Physical Review B 42 10329–10341 ISSN 0163-1829, 1095-3795
  * [31] Moustaj A, Kempkes S and Smith C M 2021 Phys. Rev. B 104 144201
  * [32] Sütő A 1989 Journal of Statistical Physics 56 525–531 ISSN 0022-4715
  * [33] Ramachandran A, Danieli C and Flach S 2018 Fano resonances in flat band networks Fano Resonances in Optics and Microwaves vol 219 ed Kamenetskii E, Sadreev A and Miroshnichenko A (Cham: Springer International Publishing) pp 311–329 ISBN 978-3-319-99730-8 978-3-319-99731-5
  * [34] Leykam D, Andreanov A and Flach S 2018 Adv. Phys. 3 1473052 ISSN null
  * [35] Leykam D, Bodyfelt J D, Desyatnikov A S and Flach S 2017 The European Physical Journal B 90 1 ISSN 1434-6028, 1434-6036
  * [36] Mielke A 2018 J Stat Phys 171 679–695 ISSN 1572-9613
  * [37] Pudleiner P and Mielke A 2015 Eur. Phys. J. B 88 207 ISSN 1434-6036
  * [38] Takayoshi S, Katsura H, Watanabe N and Aoki H 2013 Phys. Rev. A 88 063613
  * [39] Grémaud B and Batrouni G G 2017 Phys. Rev. B 95 165131
  * [40] Tovmasyan M 2018 Strongly Correlated Phases in Flatband Lattices Ph.D. thesis ETH Zurich
  * [41] Slot M R, Gardenier T S, Jacobse P H, van Miert G C P, Kempkes S N, Zevenhuizen S J M, Smith C M, Vanmaekelbergh D and Swart I 2017 Nat. Phys. 13 672–676 ISSN 1745-2473, 1745-2481
  * [42] Drost R, Ojanen T, Harju A and Liljeroth P 2017 Nat. Phys. 13 668–671 ISSN 1745-2481
  * [43] Vicencio R A, Cantillano C, Morales-Inostroza L, Real B, Mejía-Cortés C, Weimann S, Szameit A and Molina M I 2015 Phys. Rev. Lett. 114 245503
  * [44] Mukherjee S, Spracklen A, Choudhury D, Goldman N, Öhberg P, Andersson E and Thomson R R 2015 Phys. Rev. Lett. 114 245504
  * [45] Kajiwara S, Urade Y, Nakata Y, Nakanishi T and Kitano M 2016 Phys. Rev. B 93 075126
  * [46] Taie S, Ichinose T, Ozawa H and Takahashi Y 2020 Nat. Commun. 11 257 ISSN 2041-1723
  * [47] Taie S, Ozawa H, Ichinose T, Nishio T, Nakajima S and Takahashi Y 2015 Sci. Adv. 1 e1500854 ISSN 2375-2548
  * [48] Peotta S and Törmä P 2015 Nat. Commun. 6 8944 ISSN 2041-1723
  * [49] Julku A, Peotta S, Vanhala T I, Kim D H and Törmä P 2016 Phys. Rev. Lett. 117 045303
  * [50] Kobayashi K, Okumura M, Yamada S, Machida M and Aoki H 2016 Phys. Rev. B 94 214501
  * [51] Tovmasyan M, Peotta S, Törmä P and Huber S D 2016 Phys. Rev. B 94 245149
  * [52] Liang L, Vanhala T I, Peotta S, Siro T, Harju A and Törmä P 2017 Phys. Rev. B 95 024515
  * [53] Maimaiti W, Andreanov A, Park H C, Gendelman O and Flach S 2017 Phys. Rev. B 95 115135
  * [54] Rhim J W and Yang B J 2021 Advances in Physics: X 6 1901606 ISSN null
  * [55] Rhim J W and Yang B J 2019 Phys. Rev. B 99 045107 ISSN 2469-9950, 2469-9969
  * [56] Röntgen M, Morfonios C V, Brouzos I, Diakonos F K and Schmelcher P 2019 Phys. Rev. Lett. 123 080504
  * [57] Kempkes S N, Capiod P, Ismaili S, Mulkens J, Eek L, Swart I and Morais Smith C 2023 Quantum Front 2 1 ISSN 2731-6106
  * [58] Vicencio R A and Mejía-Cortés C 2013 J. Opt. 16 015706 ISSN 2040-8986
  * [59] Lieb E H 1989 Phys. Rev. Lett. 62 1201–1204
  * [60] Flach S, Leykam D, Bodyfelt J D, Matthies P and Desyatnikov A S 2014 Europhys. Lett. 105 30001 ISSN 0295-5075
  * [61] Sutherland B 1986 Phys. Rev. B 34 5208–5211
  * [62] Ha H and Yang B J 2021 Phys. Rev. B 104 165112 ISSN 2469-9950, 2469-9969
  * [63] Kohmoto M and Sutherland B 1986 Phys. Rev. B 34 3849–3853

##  Appendix A Determination of effective hopping amplitudes and energy
corrections of the quasiperiodic coupling case

In the next subsections, we will determine the energy corrections and the
renormalized Fibonacci hoppings for each cluster.

###  A.1 $v_{A}$ cluster

Figure 13: Decimation of the $v_{A}$ cluster’s effective chain. The cyan
rectangles represent the states corresponding to the four-atom molecules. They
correspond to the renormalized sites of the chain, and have an on-site energy
given by $E_{4M}$. The effective hoppings $T_{A}$ and $T_{B}$ are Fibonacci
distributed and further renormalization leads to the Fibonacci trifurcating
structure. Below, we have the two-atom molecular chain, where each blue site
represents a dimer. In this figure, only one of them appears, but for a larger
chain, there are Fibonacci distributed couplings, whence the trifurcating
structure starts after this point in the renormalization procedure. Finally,
the bottom chain corresponds to isolated (green) sites, forming the last
degenerate subspace of $H_{0}^{A}$. The effective chain for this subspace also
has couplings that are Fibonacci distributed.

Starting with the $v_{A}$ cluster, we can split it into the three categories
depicted in Figure 13. The 4-atom molecule leads to four energy eigenvalues,

$E^{A}_{4M}=\pm\frac{1}{\sqrt{2}}\sqrt{t_{b}^{2}+2t_{s}^{2}\pm\sqrt{(t_{b}^{2}+2t_{s}^{2})^{2}-4t_{s}^{4}}},$
(16)

where the subscript $4M$ reminds us that it is for the 4-atom molecule, while
the superscript $A$ indicates that it applies to the initial level $v_{A}$.
The 2-atom molecule and the isolated atom have eigenvalues, respectively,
equal to

$\displaystyle E^{A}_{2M}$ $\displaystyle=\pm t_{s},$ (17) $\displaystyle
E^{A}_{A}$ $\displaystyle=0.$ (18)

In order to calculate the effective hoppings using the Brillouin-Wigner
degenerate perturbation theory, we need to know the eigenstates of
$H^{A}_{0}$, which for the $E_{A}^{A}=0$ level are just the corresponding
isolated sites. For the 2-atom molecules, these are given by
$\ket{E^{A}_{2M}}=2^{-1/2}(\ket{i}\pm\ket{i+1})$ for some localized one-
particle state at site number $i$, corresponding to the sites coupled by the
strongest bond, and for the 4-atom molecules

$\ket{E^{A}_{4M}}=\frac{a_{0}\ket{i}+a_{1}\ket{i+1}+a_{2}\ket{i+2}+a_{3}\ket{i+3}}{\sqrt{a^{2}_{0}+a_{1}^{2}+a_{2}^{2}+a_{3}^{2}}}$

for the relevant sites $i$ to $i+3$, as shown in Figure 13. The coefficients
can be computed by solving the eigenvalue problem

$\begin{pmatrix}0&t_{s}&0&0\\\ t_{s}&0&t_{b}&0\\\ 0&t_{b}&0&t_{s}\\\
0&0&t_{s}&0\end{pmatrix}\begin{pmatrix}a_{0}\\\ a_{1}\\\ a_{2}\\\
a_{3}\end{pmatrix}=E_{4M}\begin{pmatrix}a_{0}\\\ a_{1}\\\ a_{2}\\\
a_{3}\end{pmatrix}.$

A solution to this problem is

$\displaystyle a_{0}$ $\displaystyle=1,$ $\displaystyle a_{1}$
$\displaystyle=\frac{E_{4M}}{t_{s}},$ (19) $\displaystyle a_{2}$
$\displaystyle=\frac{E^{2}_{4M}-t_{s}^{2}}{t_{b}t_{s}},$ $\displaystyle a_{3}$
$\displaystyle=\frac{E_{4M}}{t_{s}}\frac{E^{2}_{4M}-t_{s}^{2}-t_{b}^{2}}{t_{b}t_{s}},$

where we ommited the superscript $A$ for brevity. Each of the degenerate
subspaces can be associated to its own chain, with its own couplings in the
renormalization picture (see Figure 13).

The effective Hamiltonian for the $v_{A}$ cluster is found by setting the
perturbation Hamiltonian $H_{1}$ to be the one with all matrix elements
containing $t_{w}$. We remind the reader that the general expression was given
by Section 6.2.

####  $\bm{E_{4M}}$ Chain (cyan)

The couplings in the 4-atom molecular chain are given by the matrix element
$\bra{E^{A}_{4M,j}}H^{A}_{eff}\ket{E^{A}_{4M,j+1}}$, where $j$ is the
renormalized site index. To nearest order in $t_{w}$, they read

$\displaystyle T_{A}$ $\displaystyle=\frac{a_{3}}{N}\frac{t_{w}^{2}}{E_{4M}},$
(20) $\displaystyle T_{B}$
$\displaystyle=\frac{a_{3}}{N}\frac{t_{w}^{4}t_{s}}{E^{2}_{4M}(E^{2}_{4M}-t_{s}^{2})},$

where we defined $N\equiv a_{0}^{2}+a_{1}^{2}+a_{2}^{2}+a_{3}^{2}$. Note that
there are four different $E_{4M}$’s and hence four sets of coefficients $a$.

####  $\bm{E_{2M}}$ Chain (blue)

The 2-atom molecular chain’s couplings are obtained from the matrix elements
$\bra{E^{A}_{2M,j}}H^{A}_{eff}\ket{E^{A}_{2M,j+1}}$. The calculation yields

$\displaystyle T_{A}$
$\displaystyle=\pm\frac{t_{w}^{4}}{2t_{s}^{2}}\sum_{j=0}^{3}\frac{a_{3}^{(j)}}{N_{j}\left(\pm
t_{s}-E^{(j)}_{4M}\right)},$ (21) $\displaystyle T_{B}$
$\displaystyle=\frac{t_{w}^{6}}{2t_{s}^{3}}\left(\sum_{j=0}^{3}\frac{a_{3}^{(j)}}{N_{j}\left(\pm
t_{s}-E^{(j)}_{4M}\right)}\right)^{2},$

where we have now explicitely given a label to each of the four sets of
$a_{j}$’s, and by extention, $N_{j}$’s as well. The $\pm$ refers to the
bonding and anti-bonding energy levels.

####  $\bm{E_{A}}$ Chain (green)

Finally, the isolated atom chain has the following couplings

$\displaystyle T_{A}$ $\displaystyle=-\frac{t_{w}^{2}}{2t_{s}},$ (22)
$\displaystyle T_{B}$
$\displaystyle=-t_{w}^{2}\sum_{j=0}^{3}\frac{a_{3}^{(j)}}{N_{j}E^{(j)}_{4M}}.$

###  A.2 $v_{B}$ cluster

Figure 14: Decimation of the $v_{B}$ cluster’s effective chain. The cyan
rectangles represent the states corresponding to the four-atom molecules, like
in Figure 13. Below, we have the blue sites representing dimers, also similar
to Figure 13. The couplings in both cases follow a Fibonacci sequence and will
also trifurcate after this step in the renormalization procedure.

For the $v_{B}$ cluster’s effective chains, we only have two categories, as
shown in Figure 14. The first (cyan) corresponds to a 4-atom molecular chain,
with the same energy eigenvalues and eigenstates formulae as the previous
ones, but with a different set of $\\{t_{b},t_{s},t_{w}\\}$ [see Equation 10].
The second (blue) subchain corresponds to the 2-atom molecular chain. We
shall, once again, refer to them as the $E_{4M}$ and $E_{2M}$ chains,
respectively.

####  $\bm{E_{4M}}$ Chain

The couplings in this case are given by

$\displaystyle T_{A}$ $\displaystyle=\frac{t_{w}}{N}a_{3},$ (23)
$\displaystyle T_{B}$
$\displaystyle=\frac{a_{3}}{N}\frac{t_{w}^{2}t_{s}}{E_{4M}^{2}-t_{s}^{2}}.$

####  $\bm{E_{2M}}$ Chain

The couplings for this chain are

$\displaystyle T_{A}$
$\displaystyle=\pm\frac{t_{w}^{2}}{2}\sum_{j=0}^{3}\frac{a_{3}^{(j)}}{N_{j}\left(\pm
t_{s}-E_{4M}^{(j)}\right)},$ (24) $\displaystyle T_{B}$
$\displaystyle=\pm\frac{t_{w}^{3}}{2}\left(\sum_{j=0}^{3}\frac{a_{3}^{(j)}}{N_{j}\left(\pm
t_{s}-E_{4M}^{(j)}\right)}\right)^{2}.$