Many-body localization in a quasiperiodic Fibonacci chain

We study the many-body localization (MBL) properties of a chain of interacting fermions subject to a quasiperiodic potential such that the non-interacting chain is always delocalized and displays multifractality. Contrary to naive expectations, adding interactions in this systems does not enhance delocalization, and a MBL transition is observed. Due to the 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 is based on exact numerical studies of eigenstates and dynamical properties after a quench.


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][3][4][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] has been intensively studied. In such a flourishing field of research (see recent reviews [10][11][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 initial studies 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 to a different universality class for) the ETH-MBL transition [14]. Several studies suggested that slow dynamics observed in the ETH region are due to Griffiths-type regions [15,16], which are allowed (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 a disordered 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][27][28][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 for other potentials studied so far. Concretely, we study the "standard model" of MBL, a spin 1/2 XXZ chain: which can be recasted by a Jordan-Wigner transformation into a model of interacting fermions on a chain: where h BC encodes a boundary term that depends on the choice of boundary conditions 1 . The on-site potential terms h i are usually taken to be uncorrelated random variables drawn from a box distribution h i ∈ [−h, h], a binary distribution h i = ±h (+ with a probability p), or to be correlated variables, quasiperiodically varying according to the Aubry-André rule h AA i = h cos(2πωi + φ), where ω is an irrational frequency (φ a random phase). In all these cases, the model is known to undergo a many-body localization transition as h is increased [19,[32][33][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][21][22][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 to 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][39][40][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 line of 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 specific local features. The plan of the paper is as follows: in Sec. 2 we recall the Fibonacci sequence and 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 the thermal and localized phases of the model, studying the eigenstates properties in Sec. 5, and dynamical properties in Sec. 6. Finally, we gather our findings and conclude in Sec. 7.

The Fibonacci chain
The Fibonacci chain under study is a chain of interacting spinless fermions [Eq. (2)] subject to an on-site binary potential (h i = ±h) varying according to the Fibonacci rule: h i = +h (resp. h i = −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 to generate words of increasing length: A → AB → ABA → ABAAB → . . . . Notice that the length of the i th word is the i th Fibonacci number, hence the name of the model. The Fibonacci sequence is obtained by repeating the substitution procedure an infinite number Figure 1: The first five Fibonacci words, depicted as horizontal chains of rectangular tiles. Applying the substitution rule σ to a word yields the next, represented just below it. Below the largest word is printed the corresponding sequence of on-site potentials ±h.
of times. Fig. 1 displays the first Fibonacci words. On the infinite sequence, it is easy to show that A (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 Fibonacci sequence but with a chunk of length L. A relevant question is then: how many different words of 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 question is 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 numerically demanding, 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 rare regions 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 of the 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 a reflection 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 number of disorder realizations naturally leads to large statistical errors 2 .
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 of translation 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. When the fermion interaction term is absent, ∆ = 0 in Eq. (2), the model becomes easily tractable as it can be written as H = i,j c † i H i,j c j , where the L × L matrix H is the singleparticle Hamiltonian. Its eigenvectors, the single-particle excitations φ α and the associated single-particles energies α form the basic building blocks of the many-body physics: the many-body eigenstate of energy α 1 + α 2 + . . . is the Slater determinant of the corresponding single-particle excitations φ α 1 , φ α 2 , . . . .

The non-interacting Fibonacci chain
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 wavefunctions are both multifractal [38][39][40][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 high energy. Multifractality is of course also visible in the low-energy many-body This is not a solution to all problems, however, since physical quantities can be correlated from one level to the next. 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.
properties. In particular, the ground state entanglement entropy grows logarithmically with subsystem size, with quasiperiodic oscillations [47]. We conclude this brief review of the freefermions 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,51] observables. For example, the conductivity at the Fermi level (T = 0) decays as σ ∼ L −α [52] 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 an intermediate case between a periodic and a disordered system.

Many-body localization transition
Unless otherwise specified, the interaction strength is now fixed to ∆ = 1. We argue that the 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,53], the half-cut von Neumann entanglement entropy [54] and its variance [55]. We use in this section eigenstates in the so-called "infinite-temperature" limit, i.e. located in the middle of the spectrum at energy density = 0.5 where = E−E min Emax−E min and E max/min denote the energy spectrum extrema. These eigenstates are obtained via the exact shift-invert method [13] on 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
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 is extensive and coincides with the thermodynamic entropy [56]: S ∼ L 2 log 2 for = 0.5, while in a localized phase it displays a sub-extensive area law [54]. Fig. 3 (left) shows the behavior of 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: the number 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 [55]. Following Ref. [14] we distinguish the 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 the transition. Nevertheless, the position of the variance peak strongly drifts to and sample-to-sample (dashed lines) standard deviation of the entanglement entropy, divided by system size. Color coding for system size is the same in both panels. larger potentials as system size is increased, a feature previously observed in other MBL systems [33,55,57]. This lets us estimate a more precise lower bound of the transition location: In uncorrelated disordered models, the sample-to-sample variance was found to be the largest [14], while here the eigenstate-to-eigenstate contribution is dominating. This inversion of the hierarchy was also previously observed in the Aubry-André model in Ref [57] where it has been attributed to the deterministic nature of the potential, preventing the appearance of Griffiths regions. The inversion observed here is consistent with this hypothesis, since the Fibonacci potential also presents this deterministic property.

Spectral statistics
One can also study the spectral statistics using the well-known gap ratios which measure energy level repulsion [58]. Defined as r = min(g n /g n+1 , g n+1 /g n ), where g n = E n − E n−1 is the n th 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 value r Wigner-Dyson 0.53 [59], 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 states around 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, and to 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 that an MBL phase is expected for h > 3.5.

Qualitative argument for the existence of the MBL phase
At first glance, it may seem odd that an MBL phase can exist in the Fibonacci chain, given that the free chain is never Anderson localized, contrary (to the best of our knowledge) to all other cases of systems exhibiting MBL when adding uniform interactions. Note however a few examples with random interactions [26][27][28][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 single-particle 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 +h (resp. −h). We can Fourier transform this function [60]: where τ = (1 + √ 5)/2 is the golden ratio. Remark that because the Fibonacci potential is discontinuous, its harmonics decay as slow as is possible: W k ∼ 1/k. On the opposite end of the 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. [61] Monthus studied the case of a free-fermion chain with a potential whose harmonics behave as W k ∼ 1/k b , 1 ≤ b < ∞, allowing to interpolate between the free Fibonacci (b Fibo = 1) and Aubry-André (b AA = ∞) chains. Ref. [61] showed that as soon as b > 1, there exist an extended and a localized phase. Only the free Fibonacci chain b = 1 is special, in the sense that it is never extended 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 h i → h i − ∆ n i−1 + n i+1 . This modification can only smoothen the potential, hence effectively changing the exponent of the decay of the harmonics from b Fibo = 1 to b > 1. According to [61], the effective free fermion 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][35][36][37]. We thus argue that adding interactions to the picture may have the effect of smoothening the potential, effectively pushing it away from the critical b = 1 line and enabling the existence of an 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 to an additional weak random potential (albeit proceeding in a state-dependent non-monotonic fashion) [62].

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. We start with the local density n i , where the average is taken over eigenstates at = 0.5. Fig. 5 (left) shows the local density in the ETH phase. As system size is increased, the fermionic density gets more sharply peaked at n i = 1/2, converging towards a Gaussian with vanishing variance for large enough system size, the expected behavior for a thermal system [63]. Physically, this simply means that the fermions are increasingly getting more homogeneously spread out over the whole chain.

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 n i have increasingly large overlap with the local integrals of motion (commuting with the Hamiltonian) which are known to characterize the MBL phase [64]. This naturally produces symmetric peaks at n i = 0, 1 in the distribution of density, as observed in many simulations of MBL in spin chains (see e.g. [65]). 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 densities of the form cos 2 (δ) and sin 2 (δ), 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 magical angle densities takes the following form 4 : We find that this ansatz reproduces better and better the features of the numerically computed one-particle density matrix ρ i,j = c † i c j as h and L are increased. To understand why such magic angle states appear in the first place, we now turn to a toy model.
When ∆ = 1 we 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 the magic angle states give rise to the magic angle peaks. All of the three kind of states contribute to the peaks at n = 0 or 1.
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 the model and following the same perturbative approach, we are able to account for the other secondary 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 the binary nature of the Fibonacci potential, which causes the energy levels in the strong onsite 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.  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 other specificities are involved. To answer this, we should consider an on-site potential which is still binary, but no longer follows the Fibonacci sequence. A possible choice is to shuffle the Fibonacci potentials, keeping the same fraction of ±h but destroying the quasiperiodic longrange 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 a many-body localization transition [19].
The right panel of Fig. 7 compares the density distribution in the Fibonacci and random cases, 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 the higher 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 are about L allowed configurations of L potentials, there are exponentially more in the random case. 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 case only the largest secondary peaks remain visible.   . We find as expected that most eigenstates are weakly entangled, and the probability of finding large entanglement entropies decays exponentially. Besides the main peak at zero entanglement, several other peaks can again be distinguished on the top of the background, most notably for the Fibonacci sequence. These peaks can be taken into account by again considering simple ansatz states where fermions are entangled two-by-two across the entanglement cut. Such states produce an entanglement entropy

Entanglement entropy
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 [65], 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 tail observed for S ≥ log 2 in Fig. 8 is due to entanglement of more than two fermions. The tail is 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 spotted some 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.

One-particle density matrix
The one-particle density matrix is defined for a given eigenvector |ψ as ρ i,j = ψ| c † i c j |ψ . It was used to characterize the MBL transition [66], and was subsequently studied as a probe of the MBL physics [43,67,68]. The study of the one-particle density matrix is especially relevant in the free Fibonacci case since it reveals signatures of its multifractality, as we are going to discuss now. Free fermions In the absence of interactions, as discussed in Sec. 3, the Fibonacci model is critical 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 the participation entropy S q of φ α 5 : The participation entropy is a measure of how much of the realspace volume the single-particle orbitals occupy [69]. In particular, the q = 2 entropy is related to the inverse participation ratio (IPR): S 2 = − log IPR. In general, the entropy scales as S q ∼ D q × log L + b q , with 0 ≤ D q (φ α ) ≤ 1 the multifractal dimension of the orbital φ α . D q = 1 signals that the orbital occupies the whole volume of the system. It is said to be extended. D q = 0 signals a completely localized orbital, that is almost zero everywhere except on a finite number of sites. In the intermediate case 0 < D q < 1, the orbital occupies an infinite number of sites but a vanishing fraction of the volume. It is said to be multifractal. In the following, we study the scaling of the participation entropy averaged over orbitals: 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 strength h = 1. As expected, we observe non-trivial fractal dimensions 0 < D q < 1. The multifractal dimensions vary continuously with h, starting from D q (h = 0) = 1 and decreasing to D q (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,68]. 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. [66] argues that the natural orbitals are extended in the ETH phase, based on numerical computations. The right panel of Fig.  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 substentially lower than in the free case at the same potential strength. A possible explanation is as follows: while in the free case D q vanishes only in the infinite potential strength limit h → ∞, in the interacting case D q must vanish in the MBL phase, i.e. for some finite potential strength h * . We indeed find (data not shown) that the participation entropy of the natural orbitals is essentially 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 direct consequence 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 that the observed multifractal behavior is a finite-size effect and that it disappears on longer length scales (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 [51].

Dynamical probes of many-body localization
In the MBL phase, part of the information contained in an initial state remains measurable at arbitrary large times, in contrast with the ETH phase where any local detail about the initial state is quickly lost. Hence, quench protocols are a natural and experimentally accessible [70] 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(−iHt) |ψ(0) . We use initial product states of the form |ψ(0) = c † i 1 c † i 2 . . . |vac , with L/2 excitations i j , randomly chosen with the constraint that the average energy of the state is close to the energy of a thermal state at infinite temperature 6 . The time propagation is performed using Krylov-space techniques, on chains of up to L = 24 sites. We first examine the density imbalance (or density autocorrelation), defined as

Imbalance
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 localization of the fermions is retained even at infinite time. Fig. 10 shows the time evolution of the imbalance, for a system of L = 22 sites, and for different potential strengths h. At low h, the imbalance 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 [72]) 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.

Entanglement growth
Starting 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 grows quickly as S(t) ∝ t 1/z with z = 1 for "diffusive" [73], or z > 1 for subdiffusive systems [72], whereas the hallmark of the MBL phase is a slow logarithmic entanglement spreading S(t) ∝ log t [74][75][76]. Conversely, in a (non-interacting) Anderson insulator, the entanglement entropy saturates at long time to a non thermal size-independant value.

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 of view 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 [77], and were observed in other dynamical quantities, such as the mean square displacement of a wave packet [78,79]. Finally, we observe that the saturation value of 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 [73]. However, when approaching the localization transition, renormalization group theories [80,81] and effective models [15] predict that the dynamics can become subdiffusive in the presence of quenched disorder: S(t) ∝ t 1/z . This subdiffusive behavior was numerically observed in the disordered case [10,72,82] 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 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, for different system sizes. In an effort to free ourselves from finite-size effects, we checked (data not shown) that (i) the entanglement saturation value converges to the Page prediction for a thermalized system [56] 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 to identify 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 case of quench disorded systems with Griffiths regions [80,81]. Although the Fibonacci potential does not host Griffiths regions, owing to its deterministic, scale invariant nature, we observe a non-trivial z > 1 dynamical exponent in the vicinity of the transition point. Ref. [37,71,83] 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 [85]. Alternatively, Griffiths effects in the choice of the initial state have been proposed [71] to account for the slow dynamics. Note also the recent work [84], which reports the existence of a "slow" phase with respect to the spreading of information in the neighborhood of the Aubry-André ETH/MBL transition.
Logarithmic growth and MBL regime. In the MBL phase, one expects the interactions to slowly entangle the localized particles, resulting in a logarithmic growth of the entanglement entropy [76]. The lower left panel of Fig. 12 displays the entropy as a function of time in loglinear scale, for h ≥ 2.5. The entropy is seen to grow logarithmically as a function of time at 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) = A log(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 that the entropy does not grow logarithmically with time. For h 2.5, we observe a good fit to a logarithmic form, with a size-independent A prefactor 7 , in full agreement with the onset of MBL.
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 the imbalance 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 the transition point (2.5 h * 3) which is consistent with the one obtained from the analysis of eigenstates properties.

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 a paradigmatic 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 localization transition 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,71,83,84] -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 in the entanglement entropy and local observables, which we relate to the binary and long-range ordered 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 vicinity of the critical line ∆ = 0, especially given the instability of the localized phase observed in the Aubry-André model in this limit [86]. We have argued using qualitative mean-field arguments that the non-interacting chain in controlled by a fixed-point unstable to the addition of interactions. To explore this picture, it would be interesting to extend the existing renormalization group approaches to ground-states of aperiodic sequences [22,23] to high energy, in the lines of [35,80,81].
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 [87,88], which can be classified according to the amplitude of their geometrical fluctuations [89]. Accordingly, the low energy physics of such systems falls into the universality classes of clean or disordered 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 analysis easier.
Many experiments conducted with phonons [90], photons [91], microwaves [92], polaritons [93] and cold atoms [94] explore the properties of non-interacting quantum particles in a quasiperiodic landscape. A cold-atomic experimental setup has been proposed [95] to realize the Fibonacci sequence considered here, which could be used to probe the many-body physics discussed in our work. libraries, which are part of the SciPy [102] 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 ANR-10-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 2018-A0030500225).