SciPost Phys. 6, 050 (2019)

Many-body localization in a quasiperiodic Fibonacci chain
Nicolas Macé? , Nicolas Laflorencie and Fabien Alet
Laboratoire de Physique Théorique, IRSAMC, Université de Toulouse,
CNRS, UPS, 31062 Toulouse, France? mace@irsamc.ups-tlse.fr

Abstract
We study the many-body localization (MBL) properties of a chain of interacting fermionssubject to a quasiperiodic potential such that the non-interacting chain is always delocalized and displays multifractality. Contrary to naive expectations, adding interactions inthis systems does not enhance delocalization, and a MBL transition is observed. Due tothe local properties of the quasiperiodic potential, the MBL phase presents specific features, such as additional peaks in the density distribution. We furthermore investigate
the fate of multifractality in the ergodic phase for low potential values. Our analysis isbased on exact numerical studies of eigenstates and dynamical properties after a quench.
Copyright N. Macé et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 21-12-2018
Accepted 12-04-2019
Check forupdatesPublished 29-04-2019
doi:10.21468/SciPostPhys.6.4.050

Contents
1 Introduction

2

2 The Fibonacci chain

3

3 The non-interacting Fibonacci chain

4

4 Many-body localization transition
4.1 Entanglement of eigenstates
4.2 Spectral statistics
4.3 Qualitative argument for the existence of the MBL phase

5
6
7
7

5 The Fibonacci thermal and localized phases
5.1 Local fermionic density
5.1.1 ETH phase
5.1.2 MBL phase: apparition of secondary structures
5.2 Entanglement entropy
5.3 One-particle density matrix

8
8
8
8
11
12

6 Dynamical probes of many-body localization
6.1 Imbalance
6.2 Entanglement growth

14
14
15

7 Conclusion

17

References

18
1

SciPost Phys. 6, 050 (2019)1

Introduction

The question of many-body localization (MBL) aims at extending the non-interacting Anderson localization problem [1] towards more realistic systems where interparticle interactions
cannot be neglected. Precursor works highlighted the possibility of a dynamical transition
between ergodic and localized regimes [2–5]. This scenario is now well established for one
dimensional quantum interacting systems on the lattice in the presence of disorder. A high energy transition between a delocalized phase satisfying the eigenstate thermalization hypothesis (ETH) [6, 7], and a localized MBL regime with emerging integrals of motion [8, 9] hasbeen intensively studied. In such a flourishing field of research (see recent reviews [10–12]),
numerical simulations of lattice models presenting a competition between interaction and disorder play a pivotal role in our understanding of the MBL problem [13]. Most of the initialstudies concentrated on models with a random local potential drawn from a uniform (most of
the time box-shaped) distribution as a generic representative of disorder. However, the precise form of potential chosen can have some impact as was realized in subsequent studies. For
instance, a quasiperiodic potential leads to different finite-size effects around (and possibly toa different universality class for) the ETH-MBL transition [14]. Several studies suggested thatslow dynamics observed in the ETH region are due to Griffiths-type regions [15,16], which areallowed (e.g. random box distribution) or not (quasiperiodic potential) by the type of potential
imposed. The fact that the potential couples to spin and/or charge degrees of freedom in adisordered Hubbard model leads to the existence of a full or partial MBL phase [17]. Finally, a
discrete disorder distribution can also induce MBL [18, 19], despite stronger finite size effects
observed in the case of binary distributions [19, 25].
Quite remarkably, until now the MBL phase was found to be either induced by random interactions [26–29], or directly connected to an Anderson insulator [4,5]. In this work, we want
to go further by studying a model that verifies neither of the two previous properties. Specifically, we consider the case of a potential that follows the Fibonacci quasiperiodic sequence,
which displays critical delocalized eigenfunctions in the non-interacting limit, not found forother potentials studied so far. Concretely, we study the “standard model” of MBL, a spin 1/2
XXZ chain:XX
1X + −z+(1)
H=
hi Siz ,Siz Si+1−+∆
Si Si+1 + Si− Si+1
2 iiiwhich can be recasted by a Jordan-Wigner transformation into a model of interacting fermions
on a chain:
X
X

1X †
†
H=
hi ni + hBC ,(2)ni ni+1 −ci ci+1 + ci+1
ci + ∆
2 iiiwhere hBC encodes a boundary term that depends on the choice of boundary conditions1 .
The on-site potential terms hi are usually taken to be uncorrelated random variables drawn
from a box distribution hi ∈ [−h, h], a binary distribution hi = ±h (+ with a probabilityp), or to be correlated variables, quasiperiodically varying according to the Aubry-André rulehAA
= h cos(2πωi + φ), where ω is an irrational frequency (φ a random phase). In allithese cases, the model is known to undergo a many-body localization transition as h is increased [19, 32–34, 36, 37]. Here, we focus on yet another choice for the potential, sometimes referred to as the Fibonacci potential. To the best of our knowledge this choice for the
potential has never been studied in the MBL context (note however initial [20–23] and recent [24] studies on the low-energy properties of Eq. 2). In the non-interacting case ∆ = 0,
both disordered and Aubry-André models can be Anderson localized. Free fermions subject to1

In the remainder of the paper we always choose open boundary conditions (for reasons discussed in Sec. 2),
in which case hBC = −(n1 + n L )/2 up to a constant term.

2

SciPost Phys. 6, 050 (2019)a Fibonacci sequence never exhibit Anderson localization however, but rather display criticality and multifractality (in a sense defined later) for any on-site potential strength h [38–41].
Owing to the critical, multifractal, nature of its eigenstates and to the power-law behavior of
its transport observables (discussed in Sec. 3), we can think of this model as following a lineof metal-insulator transition points, as h is varied.
Can this model host an MBL phase when interactions are added? One would naively be
tempted to answer by the negative, given that interactions generally tend to favor delocalization [42,43], and that the non-interacting eigenstates are not localized to begin with. However,as we shall show below, the model under study presents an MBL phase, which hosts specificlocal features. The plan of the paper is as follows: in Sec. 2 we recall the Fibonacci sequenceand present its basic properties, while Sec. 3 presents the well-known properties of the noninteracting model and in particular its multifractal properties. Moving on to the interacting
case in Sec. 4, we study this model in light of various standard MBL probes, and conclude
positively for the existence of an MBL transition. We then characterize more precisely thethermal and localized phases of the model, studying the eigenstates properties in Sec. 5, anddynamical properties in Sec. 6. Finally, we gather our findings and conclude in Sec. 7.

2

The Fibonacci chain
A
A

B

A
A

B
B

A

A

A

B

A

B

A

A

B

A

B

A

+h

−h

+h

+h

−h

+h

−h

+h

Figure 1: The first five Fibonacci words, depicted as horizontal chains of rectangulartiles. Applying the substitution rule σ to a word yields the next, represented justbelow it. Below the largest word is printed the corresponding sequence of on-site
potentials ±h.
The Fibonacci chain under study is a chain of interacting spinless fermions [Eq. (2)] subject
to an on-site binary potential (hi = ±h) varying according to the Fibonacci rule: hi = +h (resp.
hi = −h) if the i th letter the Fibonacci sequence is an A (resp. a B).
Fibonacci sequence The Fibonacci sequence is a quasiperiodic sequence of A and B letters.
To construct it, start from the letter A, and apply repetitively the substitution rule¨
A → ABσ:(3)
B→A
to generate words of increasing length: A → AB → ABA → ABAAB → . . . . Notice that thelength of the i th word is the i th Fibonacci number, hence the name of the model. The Fibonaccisequence is obtained by repeating the substitution procedure an infinite number of times. Fig.
3

SciPost Phys. 6, 050 (2019)

1 displays the first Fibonacci words. On the infinite sequence, it is easy to show that Ap(resp.
B) letters occur with frequency 1/τ (resp. 1/τ2 ), where τ is the golden ratio, τ = (1 + 5)/2.
Because these frequencies are incommensurate, the Fibonacci sequence cannot be periodic. It
is however close to periodicity (quasiperiodic), as we explain now.
Quasiperiodicity, samples of finite size and averaging over realizations When performing experiments or numerical simulations, we do not deal with the ideal, infinite Fibonaccisequence but with a chunk of length L. A relevant question is then: how many different wordsof length L does one encounter in the Fibonacci sequence or, in other words: how many different samples do we have at our disposal for an experiment or a simulation? This questionis crucial in practice as other potentials (e.g. the uniform box, the binary or the Aubry-André
potential) all have a random component, from which a large number of different realizations
can be drawn.
It turns out that there are L + 1 distinct samples [44], each occurring with about the same
frequency on the infinite chain. This makes the statistical sampling much less numericallydemanding, as compared to random distributions. One can also prove [45] that a binary
sequence is periodic if and only if it has strictly less than L + 1 words of length L (for L larger
than the period). Thus, the Fibonacci model is in some sense the closest it can be to periodicity.
This also means that the Fibonacci sequence does not possess Griffiths regions, anomalous rareregions where the potential is markedly different from the rest of the sample. The absence of
these regions is in common with the Aubry-André potential.
Among the L + 1 samples of length L, one is reflection symmetric around the center ofthe chain. For this sample, the Hamiltonian Eq. (2) commutes with the reflection operator,
and is therefore block-diagonal. As was already noted in the context of binary disorder [19],this block-diagonal structure makes such samples difficult to compare with the others. For
that reason we chose to disregard them. The L remaining samples (for L even) all have areflection symmetric partner. Since eigenstates and eigenspectra are identical for a sample
and its partner, we are actually left with only L/2 truly different samples. This small numberof disorder realizations naturally leads to large statistical errors2 .
For open boundary conditions, the reflection about the center of the sample is the only relevant symmetry operation. For periodic boundary conditions however, the conjugate action oftranslation and reflection can actually lead to even more equivalent samples, so that several of
the L/2 remaining samples must be abandoned. Wishing to maximize the number of available
samples, we always consider in the rest of this work systems with open boundary conditions.

3

The non-interacting Fibonacci chain

When the fermion interaction term isP
absent, ∆ = 0 in Eq. (2), the model becomes easily
tractable as it can be written as H = i, j ci† Hi, j c j , where the L × L matrix H is the singleparticle Hamiltonian. Its eigenvectors, the single-particle excitations φα and the associatedsingle-particles energies εα form the basic building blocks of the many-body physics: the manybody eigenstate of energy εα1 + εα2 + . . . is the Slater determinant of the corresponding singleparticle excitations φα1 , φα2 , . . . .
The single particle properties of the Fibonacci Hamiltonian have already been extensively
studied. It was discovered that the single-particle energy spectrum and the single-particle
2

In order to alleviate the problem, we also average over energy levels within a small window for each sample.
This is not a solution to all problems, however, since physical quantities can be correlated from one level to thenext. The gap ratios discussed in Sec. 4.2 are especially correlated since they are functions of three consecutive
energy levels, and we correspondingly observe that their value is subject to particularly large fluctuations.

4

SciPost Phys. 6, 050 (2019)

0

100

5

10

15

20

0

10−2

|φα (i)|2

510−1|hc†i cj i|

j

10−4

10

10−6

15

10−8
0

100

200

300

20

i

10−2

i

Figure 2: Two manifestations of the Fibonacci multifractality in the absence ofinteractions. Left: amplitude of a single-particle excitation (h = 3, L = 377),
Right: average one-particle density matrix of high energy many-body eigenstates(h = 3, L = 21).
wavefunctions are both multifractal [38–41, 46]. We recall that a wavefunction is multifractal
if it is locally scale invariant. As an illustration, the left panel of Fig. 2 shows a Fibonacci wavefunction, whose peak-like structure repeats at all visible scales, hinting at its multifractality.
The right panel of Fig. 2 shows the averaged one-particle density matrix of many-body eigenstates. It also exhibits a scale-invariant structure, indicating that the many-body eigenstates
inherit the multifractal properties of the single-particle excitations. The density matrix of Fig. 2(right panel) was extracted from high energy eigenstates (i.e. in the middle of the many-body
spectrum), which implies that Fibonacci multifractality is very robust, surviving even at highenergy. Multifractality is of course also visible in the low-energy many-body properties. In
particular, the ground state entanglement entropy grows logarithmically with subsystem size,with quasiperiodic oscillations [47]. We conclude this brief review of the free-fermions Fibonacci chain noting that the “structural” multifractality of the spectrum and single-particle
states result in a power-law behavior for the thermodynamic [22,48,49] and transport [50,52]
observables. For example, the conductivity at the Fermi level (T = 0) decays as σ ∼ L −α [53]3 .
This behavior contrasts with the ones present for a disordered chain (exponentially decaying
conductivity) and a periodic chain (no decay).
To summarize, the quasiperiodicity of the Fibonacci modulation results in the multifractality of the single-particle spectrum and wavefunctions, concomitant with power-law thermodynamic and transport observables, thus making the non-interacting Fibonacci chain anintermediate case between a periodic and a disordered system.

4

Many-body localization transition

Unless otherwise specified, the interaction strength is now fixed to ∆ = 1. We argue thatthe interacting Fibonacci chain exhibits a transition from a thermal to a many-body localized
phase, based on the numerical study of standard observables: the level-spacing distribution[32, 54], the half-cut von Neumann entanglement entropy [55] and its variance [56]. We use
in this section eigenstates in the so-called "infinite-temperature" limit, i.e. located in the middle
E−Eminof the spectrum at energy density ε = 0.5 where ε = Emax −E
and Emax/min denote the energyminspectrum extrema. These eigenstates are obtained via the exact shift-invert method [13] on
3Note that the T → ∞ limit is more relevant to our study. In this limit, the finite-size dependence of the
conductivity has not been computed, to our knowledge.

5

SciPost Phys. 6, 050 (2019)

large chains of length L up to L = 24. We typically consider 5000 eigenstates per sample.
Dynamical probes of localization will also be discussed in Sec. 6.

Entanglement of eigenstates

S/L

10−1

state-to-state
sample-to-sample

0.05

L = 12
L = 14
L = 16
L = 18
L = 20L = 22
L = 240.04

std(S)/L

4.1

0.03
0.02
0.01
0.00

1

2

3

4

1

5

2

3

4

5

h

h

Figure 3: Left: Average half-chain entanglement entropy divided by system size,
as a function of potential strength, for various system sizes. Right: Eigenstate-toeigenstate (solid lines) and sample-to-sample (dashed lines) standard deviation of
the entanglement entropy, divided by system size. Color coding for system size is thesame in both panels.
The half-chain von Neumann entanglement entropy of an eigenstate |ψ〉 is defined as
S = −Tr(ρA log ρA) where ρA is the reduced density matrix obtained by tracing out the Hilbert
space of the left/right half chain. In a thermal phase, the half-chain entanglement entropy isextensive and coincides with the thermodynamic entropy [57]: S ∼ 2L log 2 for ε = 0.5, while
in a localized phase it displays a sub-extensive area law [55]. Fig. 3 (left) shows the behaviorof S/L (where the overline denotes average over samples and eigenstates) as a function of h,when the system size is increased from L = 12 to L = 24. For h ® 1.5, S/L visibly tends to
its maximum value while for h ¦ 3.5 it tends to zero. This is a first hint that a many-body
localization transition occurs in this system. From the scaling we infer the critical potential
strength to be 1.5 ≤ h∗ ≤ 3.5. The important fluctuations of the average entropy prevent us
from obtaining a more precise bound on the location of the transition with this estimator. Let
us emphasize here that this noise is not due to under-sampling, but intrinsic to the model: thenumber of available samples indeed only scales modestly as L/2.
The variance of the entropy divided by system size is expected to develop a peak at the MBL
to thermal transition, and to go to zero away from it [56]. Following Ref. [14] we distinguishthe sample-to-sample variance from the eigenstate-to-eigenstate variance within a sample.
The right panel of Fig. 3 shows that both contributions indeed develop a peak around thetransition. Nevertheless, the position of the variance peak strongly drifts to larger potentials
as system size is increased, a feature previously observed in other MBL systems [33, 56, 58].
This lets us estimate a more precise lower bound of the transition location: h∗ > 1.8.
In uncorrelated disordered models, the sample-to-sample variance was found to be thelargest [14], while here the eigenstate-to-eigenstate contribution is dominating. This inversionof the hierarchy was also previously observed in the Aubry-André model in Ref [58] where it
has been attributed to the deterministic nature of the potential, preventing the appearanceof Griffiths regions. The inversion observed here is consistent with this hypothesis, since theFibonacci potential also presents this deterministic property.

6

SciPost Phys. 6, 050 (2019)4.2

Spectral statistics

0.525
Wigner-Dyson
Poisson
L = 16
L = 18
L = 20L = 22
L = 24

0.500

r

0.475
0.450
0.425
0.400
0.375
1

2

3

4

5

h
Figure 4: Gap ratio (averaged over states around energy density ε = 0.5) as a function of on-site potential strength, for various system sizes. Gap ratio value for the
Wigner-Dyson and Poisson statistics are indicated by dotted and dashed lines respectively.
One can also study the spectral statistics using the well-known gap ratios which measureenergy level repulsion [59]. Defined as r = min(g n /g n+1 , g n+1 /g n ), where g n = En − En−1
is the nth spectral gap, they are shown in Fig. 4 for system sizes 16 ≤ L ≤ 24. In the thermal phase, where energy levels follow the Wigner-Dyson statistics, they approach the valuer Wigner-Dyson ' 0.53 [60], while in the MBL phase they reach the lower value r Poisson ' 0.39 predicted by the Poisson level statistics. In Fig. 4 we show r after averaging over 10% of the statesaround energy density ε = 0.5 (except for sizes 20, 22 where 5000 states were selected, while
300 states were used for size 24). The gap ratios tend to their thermal value for h ® 1.5, andto their MBL value for h ¦ 3.5. Large sample-to-sample variation prevents us from obtaining
a more precise estimate of the transition point. Nevertheless, these results appear consistent
with the behavior of the entanglement entropy. In particular they unambiguously show thatan MBL phase is expected for h > 3.5.

4.3

Qualitative argument for the existence of the MBL phaseAt first glance, it may seem odd that an MBL phase can exist in the Fibonacci chain, giventhat the free chain is never Anderson localized, contrary (to the best of our knowledge) to allother cases of systems exhibiting MBL when adding uniform interactions. Note however a few
examples with random interactions [26–29] as well as of interaction-induced localization in
the metallic phase of two-interacting particles models in a quasiperiodic potential [30, 31]. In
this section, we argue that the criticality of the non-interacting Fibonacci chain is very fragile,
with infinitesimal perturbations enough to drive the system into the more conventional singleparticle localized or extended phases, akin to the ones of the Aubry-André model.
First, let us call W j the function that returns 0 (resp. 1) if the j th Fibonacci potential is +hP
fk e ik j/τ , with(resp. −h). We can Fourier transform this function [61]: W j = −∞<k<∞ Wk e i 2τ
k
fk =
W
sinc
,
τ
2τ

7

(4)

SciPost Phys. 6, 050 (2019)p
where τ = (1 + 5)/2 is the golden ratio. Remark that because the Fibonacci potential isfk ∼ 1/k. On the opposite end of the
discontinuous, its harmonics decay as slow as is possible: W
spectrum, the Aubry-André potential enjoys infinite differentiability, and as a consequence its
harmonics decays as fast as is possible (only the first is non-zero). In Ref. [62] Monthus studied
fk ∼ 1/k b ,the case of a free-fermion chain with a potential whose harmonics behave as W
1 ≤ b < ∞, allowing to interpolate between the free Fibonacci (bFibo = 1) and Aubry-André(bAA = ∞) chains. Ref. [62] showed that as soon as b > 1, there exist an extended and alocalized phase. Only the free Fibonacci chain b = 1 is special, in the sense that it is neverextended nor localized, but critical for any potential h.
Now, we start from the free fermions Fibonacci chain and add interactions to the picture. At
the mean field level, this amounts to modifying the on-site potentials hi → hi − ∆〈ni−1 + ni+1 〉.
This modification can only smoothen the potential, hence effectively changing the exponentof the decay of the harmonics from bFibo = 1 to b > 1. According to [62], the effective freefermion model is no longer critical but instead exhibits a localization-delocalization transition,
resembling that of the Aubry-André model, which is known to exhibit an MBL phase [34–37].
We thus argue that adding interactions to the picture may have the effect of smoothening thepotential, effectively pushing it away from the critical b = 1 line and enabling the existence ofan MBL phase at large enough potential. Note also that such a scenario is consistent with recent
work showing Anderson localization of free fermions with Fibonacci modulation subject toan additional weak random potential (albeit proceeding in a state-dependent non-monotonic
fashion) [63].

5

The Fibonacci thermal and localized phases

Having established the existence of a many-body localization transition in the Fibonacci system, we now investigate in more details the content of the thermal and localized phases.

Local fermionic density

100

L = 14
L = 16
L = 18
L = 20L = 22

L = 14
L = 18
L = 22101

P (hni)

10

1

h=5

h = 0.5

102

P (hni)5.1

10−1

δ = π/8

100

10−2
10−3
0.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

0.6

0.8

1.0

hni

hni

Figure 5: Left: density distribution P(〈ni 〉) in the thermal phase (h = 0.5). Right:
density distribution in the MBL phase (h = 5). Distribution are drawn taking into account all sites and states around energy density ε = 0.5. Dashed lines: 〈n〉 = cos2 (δ),〈n〉 = sin2 (δ), δ = π/8.
We start with the local density 〈ni 〉, where the average is taken over eigenstates at ε = 0.5.

8

SciPost Phys. 6, 050 (2019)

5.1.1

ETH phase

Fig. 5 (left) shows the local density in the ETH phase. As system size is increased, the fermionicdensity gets more sharply peaked at 〈ni 〉 = 1/2, converging towards a Gaussian with vanishing variance for large enough system size, the expected behavior for a thermal system [64].
Physically, this simply means that the fermions are increasingly getting more homogeneously
spread out over the whole chain.
5.1.2

MBL phase: apparition of secondary structures

The right panel of Fig. 5 shows the local density relatively deep in the MBL phase (h = 5)
and where one can clearly see several peaks. In the large potential limit, the density operators ni have increasingly large overlap with the local integrals of motion (commuting with the
Hamiltonian) which are known to characterize the MBL phase [65]. This naturally produces
symmetric peaks at 〈ni 〉 = 0, 1 in the distribution of density, as observed in many simulationsof MBL in spin chains (see e.g. [66]). Quite remarkably, we observe several secondary peaks
which are absent in the standard MBL phase. Note also that they are not due to finite size
effects.

Numerical characterization of the secondary density peaks The left panel of Fig. 6 shows
the site-to-site variations of the fermionic density. We observe that local density distributions
are strongly correlated with the local potential landscape. In particular, on pairs of neighbouring sites subject to the same +h potential, the density distribution is peaked around densitiesof the form cos2 (δ) and sin2 (δ), where the magic angle δ takes the values δ = π/4 , π/8(marked by vertical dashed lines on the right panel of Fig. 5 and on the left panel of Fig. 6).
We also observe a weaker peak at δ = 3π/16.
Deep in the MBL phase, to a good approximation, fermions are at most pairwise entangled.
Then, a simple ansatz for pairs of neighboring entangled sites accounting for these magicalangle densities takes the following form4 :|ψ〉 = cos δ |01〉 ± sin δ |10〉 .

(5)

We find that this ansatz reproduces better and better the features of the numerically computed
one-particle density matrix ρi, j = 〈ci† c j 〉 as h and L are increased. To understand why suchmagic angle states appear in the first place, we now turn to a toy model.
Toy model for the magic angle states Let us consider a simple model of 4 fermions on a 4
sites chain. Since the magic angle densities are observed on +h, +h couples (see left panel ofFig. 6), which are always surrounded by −h, we chose for our toy model the following sequenceof on-site
Ppotentials: {−h, +h, +h, −h}. We write the XXZ Hamiltonian as H = H0 + H1 where
H0 = − i hi ni is the dominant part deep in the MBL phase. Treating H1 perturbatively, we
obtain that 4 of the 16 eigenstates are magic angle states of the form (5), withp
δ = arctan
1 + ∆2 − ∆ .
(6)When ∆ = 1 wep obtain δ = π/8 as observed numerically. Among the remaining states, 4 are(|01〉 ± |10〉)/ 2 states (Eq. (5) with δ = π/4), and the other are product states. As shown on
the right panel of Fig. 6, the δ = π/4 states give rise to the 〈n〉 = 1/2 density peak, while themagic angle states give rise to the magic angle peaks. All of the three kind of states contributeto the peaks at 〈n〉 = 0 or 1.
4

By symmetry, the state with the role of |01〉 and |10〉 exchanged is also a valid ansatz.

9

SciPost Phys. 6, 050 (2019)

magic angle states|01i±|10i√
states2

P (hni)

0.3
0.2
0.1
0.0

0

sin2 δ

1/2

cos2 δ

1hni

Figure 6: Left: Spatial distribution of the local fermionic density (darker means larger
probability), for 5000 eigenstates at the middle of the spectrum, in the MBL phaseat h = 5, for one of the 11 Fibonacci configurations of L = 22 sites. The chosenconfiguration is printed below the horizontal axis (+: +h potential, −: −h potential).
Right: Density distribution for the toy model considered in the text.
Thus, our simple 4 sites model qualitatively reproduces the observed density distribution,and predicts the correct value for the “magic angle” δ. Moreover, adding more sites to themodel and following the same perturbative approach, we are able to account for the othersecondary peaks observed at δ ' 3π/16 (visible on the left panel of Fig. 6).
Taking a step back, we observe that the existence of magic angle states stems from thebinary nature of the Fibonacci potential, which causes the energy levels in the strong on-site
potential limit to be highly degenerate. Hopping and interaction terms then couple states
within a degenerate level, giving rise to the observed “magic angle states”. In the next paragraph, we argue that these states are also favored by the correlated nature of the Fibonacci
modulation.

Fibonacci
binary disorder

P (hni)

101

100

0.0

0.2

0.4

0.6

0.8

1.0

hni
Figure 7: Comparison of Fibonacci and random binary disorder (shuffle) density
distributions, for the system of L = 22 sites in the MBL regime at h = 5.
Secondary density peaks in a random chain It is natural to ask whether the secondary
peaks are uniquely due to the discrete (binary) nature of the Fibonacci potential or if otherspecificities are involved. To answer this, we should consider an on-site potential which isstill binary, but no longer follows the Fibonacci sequence. A possible choice is to shuffle theFibonacci potentials, keeping the same fraction of ±h but destroying the quasiperiodic long10

SciPost Phys. 6, 050 (2019)

range order. This is equivalent to drawing the potential from a random distribution, with
P(+h) = 1/τ, P(−h) = 1 − 1/τ (with τ the golden ratio). In that case we also expect amany-body localization transition [19].
The right panel of Fig. 7 compares the density distribution in the Fibonacci and randomcases, for the same potential strength. The largest secondary peaks are still visible in the random case, but are lower. Moreover, the smaller secondary peaks are not visible at all. We thus
conclude that the secondary peak structure is not a specificity of the Fibonacci modulation,
but that the Fibonacci modulation favors the appearance of these peaks compared to the random modulation. This can be easily understood. Indeed, going back to the 4 sites toy model,
one can show that the only on-site potential configurations that give rise to the magic angle
states are the one allowed by the Fibonacci rule: {−h, +h, +h, −h}, and its symmetric partner
obtained by exchanging +h and −h. On the random chain considered here, these two configurations occur with a probability 2/τ6 ' 11%, while on the Fibonacci chain they occur at thehigher frequency (1/τ)3 ' 24%, due to the highly correlated nature of the Fibonacci modulation. This explains why the observed magic angle peak is less probable, but still present in the
random case. To account for the other peaks observed in the Fibonacci case, we can take into
account more than 4 sites (say L) in the toy model. While in the quasiperiodic case, there areabout L allowed configurations of L potentials, there are exponentially more in the randomcase. Thus, assuming distinct configurations mainly contribute to different peaks, the probability that a given random configuration contributes to one of the Fibonacci secondary density
peaks is exponentially suppressed as L is increased, implying in turn that in the random caseonly the largest secondary peaks remain visible.

5.2

Entanglement entropyS(δ = π/8)
S = log 2
S = 2 log 2
Fibonacci
binary disorder

P (S)

100

10−1

10−2

10−3
0

1

2

3

4

5

S

Figure 8: Distribution of the half-chain entanglement entropy, in the MBL phase(h = 5), with L = 22 sites and around energy density ε = 0.5. For comparison,
we superimpose to the Fibonacci case the distribution obtained when shuffling the
Fibonacci potentials. Solid line: entropy of the “magic angle” states Eq. (5).
Fig. 8 displays the distribution of half-chain entanglement entropy, relatively deep in the
MBL phase, both in the case of the Fibonacci potentials and in the case of the shuffled Fibonacci
potentials (see section above). We find as expected that most eigenstates are weakly entangled,
and the probability of finding large entanglement entropies decays exponentially. Besides themain peak at zero entanglement, several other peaks can again be distinguished on the topof the background, most notably for the Fibonacci sequence. These peaks can be taken into
11

SciPost Phys. 6, 050 (2019)account by again considering simple ansatz states where fermions are entangled two-by-two
across the entanglement cut. Such states produce an entanglement entropy
S (〈n〉) = −〈n〉 log〈n〉 − (1 − 〈n〉) log(1 − 〈n〉),(7)which depends on the fermion density 〈n〉 near the cut. In this approximation, the density
peak at 〈n〉 = 1/2 [right panel of Fig. 5], also observed in the standard MBL phase [66],contributes to the S = log 2 peak of the entropy distribution (indicated by a dashed line on
Fig. 8). Similarly, the “magic angle” states contribute to a peak in the entropy distribution at
S ' 0.42 (indicated by a solid line on Fig. 8). In the binary disorder case, the magic angle
density peak is much lower than in the Fibonacci case, and accordingly, the corresponding
entropy magic angle peak is also much lower, and in fact barely visible in Fig. 8. The tailobserved for S ≥ log 2 in Fig. 8 is due to entanglement of more than two fermions. The tailis fatter in the disordered case than in the Fibonacci case, hinting that, for the same potential
strength, the Fibonacci system is more localized. Moreover, the Fibonacci entropy distribution
is peaked at S = 2 log 2, the maximal entropy for a system of 4 fermions. Although we spottedsome regularities (this peak is observed only when the on-site potential sequence is −h, h on
both sides of the cut), we were not able to use the same toy model to account for this peak.

5.3

One-particle density matrix

Sq

The one-particle density matrix is defined for a given eigenvector |ψ〉 as ρi, j = 〈ψ| ci† c j |ψ〉. It
was used to characterize the MBL transition [67], and was subsequently studied as a probe of
the MBL physics [43,68,69]. The study of the one-particle density matrix is especially relevantin the free Fibonacci case since it reveals signatures of its multifractality, as we are going todiscuss now.

3.5

S q=1
fit, Dq=1 = 0.62

3.0

S q=2
fit, Dq=2 = 0.54

2.5
2.0
1.5
1.0
101

102

L

Figure 9: Scaling of the orbitals participation entropy S q with system size, for thefree (left panel) and interacting (right panel) Fibonacci chains, at h = 1.

Free fermions In the absence of interactions, as discussed in Sec. 3, the Fibonacci model iscritical in the sense that single-particle excitations are multifractal, and that transport shows
power-law behavior (with a continuously varying exponent). As the right panel of Fig. 2
shows, the one-particle density matrix also bears signs of this multifractality. This is easily
explained by the fact that the natural orbitals defined as the eigenvectors of the density matrix,
coincide with the single-particle excitations. More quantitatively, calling φα the natural orbital
associated with the occupation number (eigenvalue of the density matrix) nα , we introduce the12

SciPost Phys. 6, 050 (2019)

participation entropy Sq of φα 5 :
 L

X
1
2q
Sq (φα ) =log|φα (i)|
.
1−q
i=1

(8)

The participation entropy is a measure of how much of the realspace volume the single-particle
orbitals occupy [70]. In particular, the q = 2 entropy is related to the inverse participation
ratio (IPR): S2 = − log IPR. In general, the entropy scales as Sq ∼ Dq × log L + bq , with0 ≤ Dq (φα ) ≤ 1 the multifractal dimension of the orbital φα . Dq = 1 signals that the orbital occupies the whole volume of the system. It is said to be extended. Dq = 0 signals acompletely localized orbital, that is almost zero everywhere except on a finite number of sites.
In the intermediate case 0 < Dq < 1, the orbital occupies an infinite number of sites but avanishing fraction of the volume. It is said to be multifractal. In the following, we study thescaling of the participation entropy averaged over orbitals:PLSq =α=1 nα Sq (φα )
.
PL
n
α=1 α

(9)

The left panel of Fig. 9 shows the scaling of the participation entropy, averaged over samples
and high energy many-body eigenstates, in the free Fibonacci case, with potential strengthh = 1. As expected, we observe non-trivial fractal dimensions 0 < Dq < 1. The multifractal dimensions vary continuously with h, starting from Dq (h = 0) = 1 and decreasing
to Dq (h → ∞) = 0.
Interacting fermions In the MBL phase, the Hamiltonian writes as a sum of commuting
localized degrees of freedom (the l-bits), and thus the natural orbitals must also be localized,
as discussed in [43, 69]. In the ETH phase however, we rather expect the natural orbitals to
be extended or multifractal, even though we are not aware of a precise generic prediction
on their nature. In the case of the disordered XXZ chain, Ref. [67] argues that the naturalorbitals are extended in the ETH phase, based on numerical computations. The right panel ofFig. 9 shows the scaling of the orbitals’ entropy with system size. The scaling is compatible
with a multifractal nature of the natural orbitals. We find multifractal dimensions substantiallylower than in the free case at the same potential strength. A possible explanation is as follows:
while in the free case Dq vanishes only in the infinite potential strength limit h → ∞, inthe interacting case Dq must vanish in the MBL phase, i.e. for some finite potential strengthh∗ . We indeed find (data not shown) that the participation entropy of the natural orbitals isessentially constant with system size (with small fluctuations at strong disorder strength) in
the MBL phase.
In conclusion, the natural orbitals are multifractal in the free Fibonacci case, as a directconsequence of the multifractality of the single-particle excitations. In the interacting case,
we find the natural orbitals to remain multifractal in the ETH phase (h = 1), at least up to
the numerically accessible length scales. This hints that some signatures of the free chain
multifractal character remain present in the ETH phase. Of course, we cannot exclude thatthe observed multifractal behavior is a finite-size effect and that it disappears on longer lengthscales (this scenario is difficult to test as we are already studing the largest possible chains).
One possibility to confirm the persistence of multifractality in the interacting Fibonacci chain
would be to study transport properties, which bear signatures of the Fibonacci multifractality
in the non-interacting case [52].
5Notice that the entropies introduced in [43, 67] differ slightly from the ones we use here. However, we expect
the scaling to be the same.

13

SciPost Phys. 6, 050 (2019)6

Dynamical probes of many-body localization

In the MBL phase, part of the information contained in an initial state remains measurable atarbitrary large times, in contrast with the ETH phase where any local detail about the initialstate is quickly lost. Hence, quench protocols are a natural and experimentally accessible [71]way of distinguishing between MBL and ETH phases.
In this part we study the quench dynamics through the properties of the time-evolved state|ψ(t)〉 = exp(−iH t) |ψ(0)〉. We use initial product states of the form |ψ(0)〉 = ci† ci† . . . |vac〉,
1 2with L/2 excitations i j , randomly chosen with the constraint that the average energy of thestate is close to the energy of a thermal state at infinite temperature6 . The time propagationis performed using Krylov-space techniques, on chains of up to L = 24 sites.

6.1

Imbalance
3.5

100

3.0

2.0

ζ

0.5

1.5

0.2
0.0

10−1

1

2

1.0

3

h
10

0

10

1

10

h

I(t)

2.5

2

103

0.5

t
Figure 10: Imbalance as a function of time for various potential strengths, for the
chains of L = 22 sites. Averaging is performed over initial high energy product statesand field configurations. Inset: exponent of the power-law decay of the imbalance
as a function of potential strength (see text).
We first examine the density imbalance (or density autocorrelation), defined as



·
L ­
4X
1
1
ψ(0) .
I(t) =
ψ(0) ni (t) −ni (0) −
L i=1
2
2

(10)

While in the ETH phase the imbalance decays from its initial value I(0) = 1 to 0, in the MBL
phase it saturates at a finite value, a signature that some information on the initial localizationof the fermions is retained even at infinite time. Fig. 10 shows the time evolution of theimbalance, for a system of L = 22 sites, and for different potential strengths h. At low h, theimbalance is indeed seen to decay to 0, while at strong enough disorder (h ¦ 2.5) it saturates
at a finite value. More quantitatively, at large times, we fit our data to a power-law decay of
the form I(t) ∼ t −ζ . The inset of Fig. 10 shows the exponent (extracted using the functional
form given by Eq. (6) of [73]) as a function of h. The exponent approaches its diffusive valueζ = 1/2 when h → 0, and monotonously decreases as a function of h. It is expected to vanish(within error bars associated to the data and fitting procedure) in the vicinity of the MBL
transition.
6

In practice we randomly generate product states and take the N closest to the infinite temperature energy,
with N chosen large enough for the imbalance variance to be small.

14

SciPost Phys. 6, 050 (2019)

6.2

Entanglement growthStarting again from a product state, the half-chain von Neumann entanglement entropy is initially 0 and increases as fermions becomes more and more entangled during the time evolution.
The manner in which entanglement propagates can also be used as a probe of localization. In
the ETH phase, we expect the entanglement entropy to grow quickly as S(t) ∝ t 1/z with z = 1
for “diffusive” [74], or z > 1 for subdiffusive systems [73], whereas the hallmark of the MBL
phase is a slow logarithmic entanglement spreading S(t) ∝ log t [75–77]. Conversely, in a(non-interacting) Anderson insulator, the entanglement entropy saturates at long time to a
non thermal size-independant value.
102

15.0
12.5
1.0

h

10

10.0

0

1/z

S(t)

101

7.5
0.5

5.0

10−1
0

5

10

2.5

15

h

10−2
10

2

10

5

10

8

10

11

1014

1017

0.0

t
Figure 11: Free fermions entanglement as a function of time for various potential
strengths, starting from random high energy initial product states on the chain of
L = 610 sites. Inset: exponent of the power-law fit S(t) ∝ t 1/z .

Free Fibonacci chain In the absence of interaction, we observe (Fig. 11) that the entanglement entropy grows as a power-law: S(t) ∝ t 1/z , with a monotonously increasing exponent,
ranging from z = 1 in the clean limit (h  1) to z → ∞ in the strong quasiperiodicity
limit (h  1). This further confirms that the free Fibonacci chain is intermediate between
an Anderson localized (S(t) = cst) and a delocalized (S(t) ∝ t) system from the point ofview of its transport properties. Remark the presence of log-periodic oscillations on top of
the dominant power-law behavior. They are observed as soon as h > 0, and their amplitude
increases monotonically with h. Such oscillations are the hallmark of an underlying discrete
scale-invariance [78], and were observed in other dynamical quantities, such as the meansquare displacement of a wave packet [79, 80]. Finally, we observe that the saturation valueof the entanglement entropy slowly decreases with h. A scaling analysis (not shown) reveals
that – as expected for a system at high energy – the entropy saturates to an extensive value,
S(t → ∞, L) = A(h)L, with a non-thermal prefactor A(h) decreasing continuously with h.
Interacting chain: Power-law growth in the ETH phase At very small disorder in the ETH
regime, the entanglement entropy is expected to grow linearly: S(t) ∝ t [74]. However, whenapproaching the localization transition, renormalization group theories [81, 82] and effectivemodels [15] predict that the dynamics can become subdiffusive in the presence of quencheddisorder: S(t) ∝ t 1/z . This subdiffusive behavior was numerically observed in the disorderedcase [10, 73, 83] not only close to the critical point but also deep in the ETH regime.
For the Fibonacci sequence, the upper left panel of Fig. 12 shows the entanglement as a

15

SciPost Phys. 6, 050 (2019)

function of time in log-log scale, for various potential strengths. At low disorder (h ® 2.5), a
clear power-law behavior is observed. The upper right panel of Fig. 12 shows the dynamical
exponent 1/z obtained by fitting the entanglement entropy to the subdiffusive form above, fordifferent system sizes. In an effort to free ourselves from finite-size effects, we checked (datanot shown) that (i) the entanglement saturation value converges to the Page prediction fora thermalized system [57] as system size is increased, and (ii) the exponent’s value does not
significantly vary with system size, for large enough systems. For h ≥ 3 we were not able toidentify time scales where the entropy has a clear power-law behavior. At very low potential
the exponent z is close to unity, but continuously increases (1/z decreases) as a function of h.
In the vicinity of the transition point, z(h) is predicted to diverge as a power-law in the caseof quench disorded systems with Griffiths regions [81, 82]. Although the Fibonacci potential
does not host Griffiths regions, owing to its deterministic, scale invariant nature, we observea non-trivial z > 1 dynamical exponent in the vicinity of the transition point. Ref. [37, 72, 84]
found a similar subdiffusive behavior using the Aubry-André Griffiths regions-free potential,
even though this has been argued to be a finite-size effect [86]. Alternatively, Griffiths effects
in the choice of the initial state have been proposed [72] to account for the slow dynamics.
Note also the recent work [85], which reports the existence of a “slow” phase with respect tothe spreading of information in the neighborhood of the Aubry-André ETH/MBL transition.
h = 0.5
h = 1.0

h = 1.5
h = 2.0h = 2.5
h = 3.0

h = 3.5
h = 4.0

L = 16
L = 18
L = 20L = 22
L = 24

1/z

S(t)

0.8
100

0.6

1.0

100

A

0.4

S(t)

10−1

0.5

10−1

0.0
100

101

102

103

1

2

3

4

h

log t

Figure 12: Left: Entanglement as a function of time, in log-log scale (upper panel)and in log-linear scale (lower panel), for the chain of L = 22 sites. Right: Best fits ofthe form S(t) ∝ t 1/z (upper panel) and S(t) = Alog t (lower panel).

Logarithmic growth and MBL regime In the MBL phase, one expects the interactions toslowly entangle the localized particles, resulting in a logarithmic growth of the entanglement
entropy [77]. The lower left panel of Fig. 12 displays the entropy as a function of time in
log-linear scale, for h ≥ 2.5. The entropy is seen to grow logarithmically as a function of timeat intermediate time scales, after an initial regime, and before reaching its saturation value.
The logarithmic growth of the entanglement entropy is not easily seen, and in order to probe
it more quantitatively, we display on the lower right panel of Fig. 12 the best fit of the entropy
to a logarithmic form (S(t) = Alog(t) + cst) as a function of potential strength, for various
system sizes. For h ® 2.5, the A prefactor increases linearly with system size, indicating thatthe entropy does not grow logarithmically with time. For h ¦ 2.5, we observe a good fit toa logarithmic form, with a size-independent A prefactor7 , in full agreement with the onset ofMBL.
7

Excluding the outlying L = 18 data, whose deviation from the trend we have no explanation for.

16

SciPost Phys. 6, 050 (2019)

In conclusion, the dynamics of high energy product states is consistent with the standard
picture of MBL: in the ETH phase, we observe a power-law decay of the density imbalance,
together with a power-law growth of the entanglement entropy, while in the MBL phase theimbalance is seen to saturate to a finite value, and the entanglement entropy to grow logarithmically with time. Moreover, the analysis of the dynamics provides us with an estimate of thetransition point (2.5 ® h∗ ® 3) which is consistent with the one obtained from the analysis of
eigenstates properties.

7

Conclusion

We have studied a model of interacting spinless fermions subject to the binary Fibonacci potential. In the non-interacting limit, the Fibonacci chain has been extensively studied as aparadigmatic model of quasiperiodicity: exhibiting multifractality, a somehow intermediate
behavior between an Anderson localized and a delocalized system.
When interactions are added, we find that the system undergoes a many-body localizationtransition at a disordered strength 2 ≤ h∗ ≤ 3.5. Our claim is supported by the exact diagonalization study of static probes (spectral statistics, scaling of the entanglement entropy, statistics
of local observables) as well as dynamic ones (entanglement growth and evolution of local
observables following a quench).
In the ETH phase (h ≤ h∗ ), we observe a sublinear scaling of the natural orbitals, compatible with a survival of the free fermions multifractality. Moreover, we find sub-diffusive
dynamics for both imbalance and growth of entanglement entropy, a result which is a bit surprising – but in line with some results on the quasiperiodic Aubry-André model [37,72,84,85]– as the Fibonacci potential is free of rare Griffiths regions, which are usually thought [15, 16]
to cause such power-law behaviors.
At large potential strength (h ≥ h∗ ), the Fibonacci MBL exhibits some specific features inthe entanglement entropy and local observables, which we relate to the binary and long-rangeordered nature of the Fibonacci modulation.
Future directions In this article we have studied the Fibonacci chain in the non-interacting
limit, and at fixed interaction strength ∆ = 1. A possible follow up would be to draw the phase
diagram in the (∆, h) plane. It would in particular be interesting – albeit technically difficult
due to the proximity to integrability – to understand the fate of the transition in the vicinityof the critical line ∆ = 0, especially given the instability of the localized phase observed in the
Aubry-André model in this limit [87]. We have argued using qualitative mean-field argumentsthat the non-interacting chain is controlled by a fixed-point unstable to the addition of interactions. To explore this picture, it would be interesting to extend the existing renormalizationgroup approaches to ground-states of aperiodic sequences [22, 23] to high energy, in the linesof [35, 81, 82].
Quasiperiodicity is usually encoded in the hopping terms [22, 48] or in the on-site potentials (this work). Another choice, which may result in new physical properties, is to encode
quasiperiodicity in the interaction term.
The Fibonacci sequence studied here belongs to a family of aperiodic sequences [88, 89],
which can be classified according to the amplitude of their geometrical fluctuations [90]. Accordingly, the low energy physics of such systems falls into the universality classes of clean ordisordered systems, or can become non-universal [23, 24]. The high-energy physics studied
here in the context of the Fibonacci chain remains to be explored for the other members of
the aperiodic sequences family. Moreover, sequences with stronger geometrical fluctuations
such as the one studied in [24] yield more finite-size samples, potentially making their scaling
17

SciPost Phys. 6, 050 (2019)analysis easier.
Many experiments conducted with phonons [91], photons [92], microwaves [93], polaritons [94] and cold atoms [95] explore the properties of non-interacting quantum particles in aquasiperiodic landscape. A cold-atomic experimental setup has been proposed [96] to realize
the Fibonacci sequence considered here, which could be used to probe the many-body physics
discussed in our work.

Acknowledgements
We thank Fabian Heidrich-Meissner, Anuradha Jagannathan, Marco Tarzia, Jean-Marc Luck,
Cécile Monthus and Marko Žnidarič for fruitful discussions. We thank the anonymous referees for their suggestions and remarks. The shift-invert [13] and time-evolution computations
were performed using the PETSc [97], SLEPc [98] and STRUMPACK [99, 100] libraries. Data
analysis and plotting was performed using the NumPy, Pandas [101] and Matplotlib [102]libraries, which are part of the SciPy [103] ecosystem.
Funding information This work benefited from the support of the project THERMOLOC
ANR-16-CE30-0023-02 of the French National Research Agency (ANR) and by the French Programme Investissements d’Avenir under the program ANR-11-IDEX-0002-02, reference ANR10-LABX-0037-NEXT. We acknowledge PRACE for awarding access to HLRS’s Hazel Hen computer based in Stuttgart, Germany under grant number 2016153659, as well as the use of
HPC resources from CALMIP (grants 2017-P0677 and 2018-P0677) and GENCI (Grant 2018A0030500225).

References[1] F. Evers and A. D. Mirlin, Anderson transitions, Rev. Mod. Phys. 80, 1355 (2008),
doi:10.1103/RevModPhys.80.1355.
[2] B. L. Altshuler, Y. Gefen, A. Kamenev and L. S. Levitov, Quasiparticle lifetime ina finite system: A nonperturbative approach, Phys. Rev. Lett. 78, 2803 (1997),
doi:10.1103/PhysRevLett.78.2803.
[3] Ph. Jacquod and D. L. Shepelyansky, Emergence of quantum chaos in finite interacting
Fermi systems, Phys. Rev. Lett. 79, 1837 (1997), doi:10.1103/PhysRevLett.79.1837.
[4] I. V. Gornyi, A. D. Mirlin and D. G. Polyakov, Interacting electrons in disorderedwires: Anderson localization and low-T transport, Phys. Rev. Lett. 95, 206603 (2005),
doi:10.1103/PhysRevLett.95.206603.
[5] D. M. Basko, I. L. Aleiner and B. L. Altshuler, Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states, Ann. Phys. 321, 1126
(2006), doi:10.1016/j.aop.2005.11.014.
[6] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046
(1991), doi:10.1103/PhysRevA.43.2046.
[7] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994),
doi:10.1103/PhysRevE.50.888.

18

SciPost Phys. 6, 050 (2019)[8] J. Z. Imbrie, Diagonalization and many-body localization for a disordered quantum spin
chain, Phys. Rev. Lett. 117, 027201 (2016), doi:10.1103/PhysRevLett.117.027201.
[9] V. Ros, M. Müller and A. Scardicchio, Integrals of motion in the many-body localized
phase, Nucl. Phys. B 891, 420 (2015), doi:10.1016/j.nuclphysb.2014.12.014.
[10] D. J. Luitz and Y. Bar Lev, The ergodic side of the many-body localization transition, Ann.
phys. 529, 1600350 (2017), doi:10.1002/andp.201600350.
[11] D. A. Abanin and Z. Papić, Recent progress in many-body localization, Ann. phys. 529,
1700169 (2017), doi:10.1002/andp.201700169.
[12] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics,
Comptes Rendus Phys. 19, 498 (2018), doi:10.1016/j.crhy.2018.03.003.
[13] F. Pietracaprina, N. Macé, D. J. Luitz and F. Alet, Shift-invert diagonalization of large many-body localizing spin chains, SciPost Phys. 5, 045 (2018),
doi:10.21468/SciPostPhys.5.5.045.
[14] V. Khemani, S. P. Lim, D. N. Sheng and D. A. Huse, Critical properties of the many-body localization transition, Phys. Rev. X 7, 021013 (2017), doi:10.1103/PhysRevX.7.021013.
[15] A. Nahum, J. Ruhman and D. A. Huse, Dynamics of entanglement and transport inone-dimensional systems with quenched randomness, Phys. Rev. B 98, 035118 (2018),
doi:10.1103/PhysRevB.98.035118.
[16] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan, D. A. Huse and M. Knap, Rareregion effects and dynamics near the many-body localization transition, Ann. phys 529,
1600326 (2017), doi:10.1002/andp.201600326.
[17] P. Prelovšek, O. S. Barišić and M. Žnidarič, Absence of full many-body localization in the disordered Hubbard chain, Phys. Rev. B 94, 241104 (2016),
doi:10.1103/PhysRevB.94.241104.
[18] B. Tang, D. Iyer and M. Rigol, Quantum quenches and many-body localization in the thermodynamic limit, Phys. Rev. B 91, 161109 (2015), doi:10.1103/PhysRevB.91.161109.
[19] J. Janarek, D. Delande and J. Zakrzewski, Discrete disorder models for many-body localization, Phys. Rev. B 97, 155133 (2018), doi:10.1103/PhysRevB.97.155133.
[20] J. Vidal, D. Mouhanna and T. Giamarchi, Correlated fermions in a onedimensional quasiperiodic potential, Phys. Rev. Lett. 83, 3908 (1999),
doi:10.1103/PhysRevLett.83.3908.
[21] J. Vidal, D. Mouhanna and T. Giamarchi, Interacting fermions in self-similar potentials,
Phys. Rev. B 65, 014201 (2001), doi:10.1103/PhysRevB.65.014201.
[22] A. P. Vieira, Low-energy properties of aperiodic Quantum spin chains, Phys. Rev. Lett. 94,
077201 (2005), doi:10.1103/PhysRevLett.94.077201.
[23] A. P. Vieira, Aperiodic quantum XXZ chains: Renormalization-group results, Phys. Rev. B
71, 134408 (2005), doi:10.1103/PhysRevB.71.134408.
[24] A. P. Vieira and J. A. Hoyos, Emergent dimerization and localization in disordered quantum chains, Phys. Rev. B 98, 104203 (2018), doi:10.1103/PhysRevB.98.104203.

19

SciPost Phys. 6, 050 (2019)[25] T. Enss, F. Andraschko and J. Sirker, Many-body localization in infinite chains, Phys. Rev.
B 95, 045121 (2017), doi:10.1103/PhysRevB.95.045121.
[26] Y. Bar Lev, D. R. Reichman and Y. Sagi, Many-body localization in system witha completely delocalized single-particle spectrum, Phys. Rev. B 94, 201116 (2016),
doi:10.1103/PhysRevB.94.201116.
[27] X. Li, D.-L. Deng, Y.-L. Wu and S. Das Sarma, Statistical bubble localization with random
interactions, Phys. Rev. B 95, 020201 (2017), doi:10.1103/PhysRevB.95.020201.
[28] P. Sierant, D. Delande and J. Zakrzewski, Many-body localization due to random interactions, Phys. Rev. A 95, 021601 (2017), doi:10.1103/PhysRevA.95.021601.
[29] P. Sierant, D. Delande and J. Zakrzewski, Many-body localization for randomly interacting bosons, Acta Phys. Pol. A 132, 1707 (2017), doi:10.12693/APhysPolA.132.1707.
[30] D. L. Shepelyansky, Two interacting particles in the Harper model, Phys. Rev. B 54, 14896(1996), doi:10.1103/PhysRevB.54.14896.
[31] A. Barelli, J. Bellissard, P. Jacquod and D. L. Shepelyansky, Double butterfly spectrum
for two interacting particles in the Harper model, Phys. Rev. Lett. 77, 4752 (1996),
doi:10.1103/PhysRevLett.77.4752.
[32] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411(2010), doi:10.1103/PhysRevB.82.174411.
[33] D. J. Luitz, N. Laflorencie and F. Alet, Many-body localization edge in the random-field
Heisenberg chain, Phys. Rev. B 91, 081103 (2015), doi:10.1103/PhysRevB.91.081103.
[34] S. Iyer, V. Oganesyan, G. Refael and D. A. Huse, Many-body localization in a quasiperiodic
system, Phys. Rev. B 87, 134202 (2013), doi:10.1103/PhysRevB.87.134202.
[35] S.-X. Zhang and H. Yao, Universal properties of many-body localizationtransitions in quasiperiodic systems, Phys. Rev. Lett. 121, 206601 (2018),
doi:10.1103/PhysRevLett.121.206601.
[36] P. Naldesi, E. Ercolessi and T. Roscilde, Detecting a many-body mobility edge with quantum quenches, SciPost Phys. 1, 010 (2016), doi:10.21468/SciPostPhys.1.1.010.
[37] M. Lee, T. R. Look, S. P. Lim and D. N. Sheng, Many-body localization in
spin chain systems with quasiperiodic fields, Phys. Rev. B 96, 075146 (2017),
doi:10.1103/PhysRevB.96.075146.
[38] S. Ostlund, R. Pandit, D. Rand, H. Joachim Schellnhuber and E. D. Siggia, Onedimensional Schrödinger equation with an almost periodic potential, Phys. Rev. Lett. 50,
1873 (1983), doi:10.1103/PhysRevLett.50.1873.
[39] P. Kalugin, A. Y. Kitaev and L. Levitov, Electron spectrum of a one-dimensional quasicrystal,Zh. Eksp. Teor. Fiz. 91, 692 (1986).
[40] J. M. Luck and D. Petritis, Phonon spectra in one-dimensional quasicrystals, J. Stat. Phys.
42, 289 (1986), doi:10.1007/BF01127714.
[41] M. Kohmoto, B. Sutherland and C. Tang, Critical wave functions and a Cantor-set
spectrum of a one-dimensional quasicrystal model, Phys. Rev. B 35, 1020 (1987),
doi:10.1103/PhysRevB.35.1020.
20

SciPost Phys. 6, 050 (2019)[42] D. L. Shepelyansky, Coherent propagation of two interacting particles in a random potential, Phys. Rev. Lett. 73, 2607 (1994), doi:10.1103/PhysRevLett.73.2607.
[43] S. Bera, T. Martynec, H. Schomerus, F. Heidrich-Meisner and J. H. Bardarson, Oneparticle density matrix characterization of many-body localization, Ann. phys 529,
1600356 (2017), doi:10.1002/andp.201600356.
[44] M. Lothaire, Algebraic combinatorics on words, Cambridge Univ. Press, Cambridge(2002), doi:10.1017/CBO9781107326019.
[45] N. P. Fogg, V. Berthé, S. Ferenczi, C. Mauduit and A. Siegel, Substitutions
in dynamics, arithmetics and combinatorics, Springer Berlin Heidelberg (2002),
doi:10.1007/b13861.
[46] Q. Niu and F. Nori, Spectral splitting and wave-function scaling in quasicrystalline and hierarchical structures, Phys. Rev. B 42, 10329 (1990), doi:10.1103/PhysRevB.42.10329.
[47] F. Iglói, R. Juhász and Z. Zimborás, Entanglement entropy of aperiodic quantum spin
chains, Europhys. Lett. 79, 37001 (2007), doi:10.1209/0295-5075/79/37001.
[48] J. M. Luck and Th. M. Nieuwenhuizen, A soluble quasi-crystalline magnetic model: The XY
quantum spin chain, Europhys. Lett. 2, 257 (1986), doi:10.1209/0295-5075/2/4/001.
[49] J. Hermisson, Aperiodic and correlated disorder in XY chains: exact results, J. Phys. A:
Math. Gen. 33, 57 (1999), doi:10.1088/0305-4470/33/1/304.
[50] S. Roche, G. Trambly de Laissardière and D. Mayou, Electronic transport properties ofquasicrystals, J. Math. Phys. 38, 1794 (1997), doi:10.1063/1.531914.
[51] J. Bellissard, Coherent and dissipative transport in aperiodic solids: An overview, In
Dynamics of Dissipation 413, Springer Berlin Heidelberg (2002), doi:10.1007/3-54046122-1_18.
[52] V. Kerala Varma, C. de Mulatier and M. Žnidarič, Fractality in nonequilibrium steady states of quasiperiodic systems, Phys. Rev. E 96, 032130 (2017),
doi:10.1103/PhysRevE.96.032130.
[53] N. Macé, A. Jagannathan, P. Kalugin, R. Mosseri and F. Piéchon, Critical eigenstates
and their properties in one- and two-dimensional quasicrystals, Phys. Rev. B 96, 045138(2017), doi:10.1103/PhysRevB.96.045138.
[54] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature,
Phys. Rev. B 75, 155111 (2007), doi:10.1103/PhysRevB.75.155111.
[55] B. Bauer and C. Nayak, Area laws in a many-body localized state and its implications for topological order, J. Stat. Mech. P09005 (2013), doi:10.1088/17425468/2013/09/P09005.
[56] J. A. Kjäll, J. H. Bardarson and F. Pollmann, Many-body localization ina disordered quantum ising chain, Phys. Rev. Lett. 113, 107204 (2014),
doi:10.1103/PhysRevLett.113.107204.
[57] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993),
doi:10.1103/PhysRevLett.71.1291.

21

SciPost Phys. 6, 050 (2019)[58] V. Khemani, D. N. Sheng and D. A. Huse, Two universality classes for the
many-body localization transition, Phys. Rev. Lett. 119, 075702 (2017),
doi:10.1103/PhysRevLett.119.075702.
[59] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature,
Phys. Rev. B 75, 155111 (2007), doi:10.1103/PhysRevB.75.155111.
[60] Y. Y. Atas, E. Bogomolny, O. Giraud and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013),
doi:10.1103/PhysRevLett.110.084101.
[61] A. Sütő, Singular continuous spectrum on a cantor set of zero Lebesgue measure for the
Fibonacci Hamiltonian, J. Stat. Phys. 56, 525 (1989), doi:10.1007/BF01044450.
[62] C. Monthus, Multifractality in the generalized Aubry-Andre quasiperiodic localization
model with power-law hoppings or power-law Fourier coefficients, Fractals, 1950007(2018), doi:10.1142/S0218348X19500075.
[63] A. Jagannathan, P. Jeena and M. Tarzia, Nonmonotonic crossover and scaling behavior in a disordered one-dimensional quasicrystal, Phys. Rev. B 99, 054203 (2019),
doi:10.1103/PhysRevB.99.054203.
[64] W. Beugeling, R. Moessner and M. Haque, Finite-size scaling of eigenstate thermalization,
Phys. Rev. E 89, 042112 (2014), doi:10.1103/PhysRevE.89.042112.
[65] D. A. Huse, R. Nandkishore and V. Oganesyan, Phenomenology of fully many-bodylocalized systems, Phys. Rev. B 90, 174202 (2014), doi:10.1103/PhysRevB.90.174202.
[66] S. P. Lim and D. N. Sheng, Many-body localization and transition by density matrix renormalization group and exact diagonalization studies, Phys. Rev. B 94, 045111 (2016),
doi:10.1103/PhysRevB.94.045111.
[67] S. Bera, H. Schomerus, F. Heidrich-Meisner and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015),
doi:10.1103/PhysRevLett.115.046603.
[68] T. L. M. Lezama, S. Bera, H. Schomerus, F. Heidrich-Meisner and J. H. Bardarson, Oneparticle density matrix occupation spectrum of many-body localized states after a global
quench, Phys. Rev. B 96, 060202 (2017), doi:10.1103/PhysRevB.96.060202.
[69] B. Villalonga, X. Yu, D. J. Luitz and B. K. Clark, Exploring one-particle orbitals in large many-body localized systems, Phys. Rev. B 97, 104406 (2018),
doi:10.1103/PhysRevB.97.104406.
[70] T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia and B. I. Shraiman, Fractal measures and their singularities: The characterization of strange sets, Phys. Rev. A 33, 1141
(1986), doi:10.1103/PhysRevA.33.1141.
[71] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk,
E. Altman, U. Schneider and I. Bloch, Observation of many-body localization ofinteracting fermions in a quasirandom optical lattice, Science 349, 842 (2015),
doi:10.1126/science.aaa7432.
[72] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider and I.
Bloch, Observation of slow dynamics near the many-body localization transition
in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017),
doi:10.1103/PhysRevLett.119.260401.
22

SciPost Phys. 6, 050 (2019)[73] D. J. Luitz, N. Laflorencie and F. Alet, Extended slow dynamical regime close
to the many-body localization transition, Phys. Rev. B 93, 060201 (2016),
doi:10.1103/PhysRevB.93.060201.
[74] H. Kim and D. A. Huse, Ballistic spreading of entanglement in a diffusive nonintegrable
system, Phys. Rev. Lett. 111, 127205 (2013), doi:10.1103/PhysRevLett.111.127205.
[75] M. Žnidarič, T. Prosen and P. Prelovšek, Many-body localization in the Heisenberg XXZ magnet in a random field, Phys. Rev. B 77, 064426 (2008),
doi:10.1103/PhysRevB.77.064426.
[76] J. H. Bardarson, F. Pollmann and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012),
doi:10.1103/PhysRevLett.109.017202.
[77] M. Serbyn, Z. Papić and D. A. Abanin, Universal slow growth of entanglement
in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013),
doi:10.1103/PhysRevLett.110.260601.
[78] D. Sornette, Discrete-scale invariance and complex dimensions, Phys. Rep. 297, 239(1998), doi:10.1016/S0370-1573(97)00076-8.
[79] I. Guarneri and G. Mantica, Multifractal energy spectra and their dynamical implications,Phys. Rev. Lett. 73, 3379 (1994), doi:10.1103/PhysRevLett.73.3379.
[80] S. Thiem, M. Schreiber and U. Grimm, Wave packet dynamics, ergodicity,and localization in quasiperiodic chains, Phys. Rev. B 80, 214203 (2009),
doi:10.1103/PhysRevB.80.214203.
[81] R. Vosk, D. A. Huse and E. Altman, Theory of the many-body localization transition in onedimensional systems, Phys. Rev. X 5, 031032 (2015), doi:10.1103/PhysRevX.5.031032.
[82] A. C. Potter, R. Vasseur and S. A. Parameswaran, Universal properties of many-body delocalization transitions, Phys. Rev. X 5, 031033 (2015), doi:10.1103/PhysRevX.5.031033.
[83] Y. Bar Lev, G. Cohen and D. R. Reichman, Absence of diffusion in an interacting system of
spinless fermions on a one-dimensional disordered lattice, Phys. Rev. Lett. 114, 100601(2015), doi:10.1103/PhysRevLett.114.100601.
[84] Y. Bar Lev, D. M. Kennes, C. Klöckner, D. R. Reichman and C. Karrasch, Transport inquasiperiodic interacting systems: From superdiffusion to subdiffusion, Europhys. Lett.
119, 37003 (2017), doi:10.1209/0295-5075/119/37003.
[85] S. Xu, X. Li, Y.-T. Hsu, B. Swingle and S. Das Sarma Butterfly effect in interacting AubryAndre model: thermalization, slow scrambling, and many-body localization (2019),
arXiv:1902.07199[86] F. Setiawan, D.-L. Deng and J. H. Pixley, Transport properties across the many-body localization transition in quasiperiodic and random systems, Phys. Rev. B 96, 104205 (2017),
doi:10.1103/PhysRevB.96.104205.
[87] M. Žnidarič and M. Ljubotina, Interaction instability of localization in quasiperiodic systems, Proc. Natl. Acad. Sci. USA 115, 4595 (2018), doi:10.1073/pnas.1800589115.
[88] S. Aubry, C. Godrèche and J. M. Luck, Scaling properties of a structure intermediate between quasiperiodic and random, J. Stat. Phys. 51, 1033 (1988),
doi:10.1007/BF01014899.
23

SciPost Phys. 6, 050 (2019)[89] C. Godrèche and J. M. Luck, Indexing the diffraction spectrum of a non-Pisot self-similar
structure, Phys. Rev. B 45, 176 (1992), doi:10.1103/PhysRevB.45.176.
[90] J. M. Luck, A Classification of critical phenomena on quasi-crystals and other aperiodic
structures, Europhys. Lett. 24, 359 (1993), doi:10.1209/0295-5075/24/5/007.
[91] W. Steurer and D. Sutter-Widmer, Photonic and phononic quasicrystals, J. Phys. D: Appl.
Phys. 40, R229 (2007), doi:10.1088/0022-3727/40/13/R01.
[92] Y. E. Kraus, Y. Lahini, Z. Ringel, M. Verbin and O. Zilberberg, Topological states
and adiabatic pumping in quasicrystals, Phys. Rev. Lett. 109, 106402 (2012),
doi:10.1103/PhysRevLett.109.106402.
[93] P. Vignolo, M. Bellec, J. Böhm, A. Camara, J.-M. Gambaudo, U. Kuhl and F.
Mortessagne, Energy landscape in a Penrose tiling, Phys. Rev. B 93, 075141 (2016),
doi:10.1103/PhysRevB.93.075141.
[94] D. Tanese et al., Fractal energy spectrum of a polariton gas in a Fibonacci quasiperiodic
potential, Phys. Rev. Lett. 112, 146404 (2014), doi:10.1103/PhysRevLett.112.146404.
[95] L. Guidoni, B. Dépret, A. di Stefano and P. Verkerk, Atomic diffusion in an
optical quasicrystal with five-fold symmetry, Phys. Rev. A 60, R4233 (1999),
doi:10.1103/PhysRevA.60.R4233.
[96] K. Singh, K. Saha, S. A. Parameswaran and D. M. Weld, Fibonacci opticallattices for tunable quantum quasicrystals, Phys. Rev. A 92, 063426 (2015),
doi:10.1103/PhysRevA.92.063426.
[97] S. Balay et al., PETSc/Tao, (2018), http://www.mcs.anl.gov/petsc.
[98] V. Hernández, J. E. Román and V. Vidal, SLEPc: A scalable and flexible toolkit forthe solution of eigenvalue problems, ACM Trans. Math. Software 31, 351 (2005),
doi:10.1145/1089014.1089019.
[99] P. Ghysels, X. S. Li, F.-H. Rouet, S. Williams and A. Napov, An efficient multicore implementation of a novel HSS-structured multifrontal solver using randomized sampling,
SIAM J. Sci. Comput. 38, S358 (2016), doi:10.1137/15M1010117.
[100] P. Ghysels, X. Sherry Li, C. Gorman and F.-H. Rouet, A robust parallel preconditioner for
indefinite systems using hierarchical matrices and randomized sampling, IEEE Int. Parallel
Distrib. Process. Symp. (IPDPS), 897 (2017), doi:10.1109/IPDPS.2017.21.
[101] W. McKinney, Data structures for statistical computing in Python, Proc. 9th Python Sci.
Conf., 51 (2010).
[102] J. D. Hunter, Matplotlib: A 2D graphics environment, Comput. Sci. Eng. 9, 90 (2007),
doi:10.1109/MCSE.2007.55.
[103] E. Jones, T. Oliphant and P. Peterson, SciPy: Open source scientific tools for Python,(2001), www.scipy.org.

24