Many-body localization in the interpolating Aubry-Andr\'e-Fibonacci model

We investigate the localization properties of a spin chain with an antiferromagnetic nearest-neighbour coupling, subject to an external quasiperiodic on-site magnetic field. The quasiperiodic modulation interpolates between two paradigmatic models, namely the Aubry-Andr\'e and the Fibonacci models. We find that stronger many-body interactions extend the ergodic phase in the former, whereas they shrink it in the latter. Furthermore, the many-body localization transition points at the two limits of the interpolation appear to be continuously connected along the deformation. As a result, the position of the many-body localization transition depends on the interaction strength for an intermediate degree of deformation of the quasiperiodic modulation. Moreover, in the region of parameter space where the single-particle spectrum contains both localized and extended states, many-body interactions induce an anomalous effect: weak interactions localize the system, whereas stronger interactions enhance ergodicity. We map the model's localization phase diagram using the decay of the quenched spin imbalance in relatively long chains. This is accomplished employing a time-dependent variational approach applied to a matrix product state decomposition of the many-body state. Our model serves as a rich playground for testing many-body localization under tunable potentials.


I. INTRODUCTION
The study of material properties most commonly begins by assuming a periodic crystalline structure within which extended Bloch waves manifest and form dispersive bands [1]. Breaking the periodicity by disorder [2] or by incommensurate potentials [3] can localize the singleparticle states of the system, leading to insulating behavior. Keeping the periodicity instead but introducing many-body interactions also breaks the Bloch picture and can lead to localization, e.g., through the formation of a charge density wave or a Mott insulator [4]. The fate of the system's conduction properties in the presence of both many-body effects and (quasi-)disorder is very rich and depends on numerous details. Specifically, the physics of quantum many-body interacting systems can be strongly altered by the presence of disorder, which can drive them from an ergodic to a many-body localized (MBL) phase [5][6][7][8][9][10][11][12]. The latter is interesting due to the fact that a closed MBL system does not thermalize at any time scale and remains robust to small perturbations, such as changing the interaction strength, the strength of disorder, and/or the temperature. Correspondingly, the dynamics of the system is frozen deep within the MBL phase, while close to the transition point between the ergodic and the MBL phase, power laws appear in transport properties. Systems exhibiting MBL are therefore interesting both from (i) the perspective of fundamental science, since they provide a test-bed for general mechanisms by which ergodicity is broken in quantum systems, * as3157@cam.ac.uk and (ii) their technological applications, since they can host types of order that cannot be present in equilibrium, which can be used for information storage, and isolation of quantum information processing devices.
Due to its complexity, the majority of works in the field of many-body localization consider one-dimensional systems, which are numerically more tractable. In one-dimensional systems, the single-particle spectrum is localized by randomly distributed disorder [2], and the many-body interactions can delocalize the system through hybridization between the localized orbits or fail to do so such that an MBL phase forms [6,7,9,10]. Noninteracting quasiperiodic systems instead show richer localization phenomena in one dimension compared to randomly disordered systems. Depending on the specific model, they are known to host (i) a pure point localization transition at a finite threshold for the Aubry-André (AA) model [13,14], (ii) critical states for the Fibonacci chain [15,16], as well as (iii) mobility edges [17][18][19]. Interestingly, also for such quasiperiodic systems it was shown theoretically [20][21][22][23][24][25][26] and experimentally [18,27,28] MBL phases appear, but the mechanism that leads to its formation depends on the fine details of the model.
For example, the quasiperiodicity in the AA model arises from an on-site cosine modulation that is incommensurate with the underlying periodic lattice spacing. The Fibonacci model instead involves two discrete onsite values that interchange throughout the system according to the Fibonacci sequence. The noninteracting Aubry-André model is known to have a metal-toinsulator transition at finite strength of quasiperiodic modulation simultaneously for all eigenstates [13,14], while the Fibonacci model always has critical eigenstates that are fractal [15,16]. Adding many-body interactions to the AA model results in the shift of the localization transition in favour of an ergodic phase [23], where transport is diffusive [29,30]. In the Fibonacci model, instead, interactions destroy the fractal critical states and introduce an MBL phase transition [24,26]. Crucially, the full details behind the MBL transition in interacting quasiperiodic models remain unknown. Hence, investigation of different interacting quasiperiodic models can provide us with more insight into the interplay between interaction and the models' exotic single-particle localization phenomena.
Interestingly, the AA and Fibonacci models can be viewed as two limits of the interpolating Aubry-André-Fibonacci model [31] (IAAF). The experimental realization of this model was reported recently in Ref. [19], wherein some of us have explored the localization phase diagram of the noninteracting IAAF model, and revealed how the critical states of the Fibonacci model form alongside a cascade of localization-delocalization transitions as the model is tuned from the AA to the Fibonacci limit. This cascade of transitions occurs nonuniformly throughout the spectrum developing many non-trivial mobility edges in the Hermitian [19] and non-Hermitian [32] version of the model. Although the physics of the noninteracting model is now well understood, the model's many-body localization properties have thus far not been explored. Specifically, how does the rich localization physics of the single-particle IAAF model with its uniform and nonuniform metal-toinsulator transitions, critical states, and mobility edges, interplay with many-body interactions? More precisely, a single-particle spectrum that is not uniformly localized (mobility edges) is predicted to exhibit a peculiar interplay with many-body interactions tending to delocalization through the ergodic parts of the single-particle spectrum [33,34].
In this work, we set out to investigate the localization phase diagram of the interacting IAAF (i-IAAF) model and answer two main questions: (i) How are the MBL transition points of the AA and Fibonacci models connected? (ii) What happens to the cascade of localizationdelocalization transitions of the noninteracting model in the presence of interactions?. In the latter, we have regions in parameter space where a nontrivial mobility edge manifests in the noninteracting limit, with coexistence of extended and localized states in the system. In this regime, we observe anomalous behavior with increasing interaction strength, where weak interactions localize the system, while stronger interactions again drive it to the ergodic phase. Our results provide a first glimpse into the complex many-body physics found in the i-IAAF model, opening the path for additional studies in this tunable setting.
The paper is outlined as follows: in Section II, we introduce the Hamiltonian of the IAAF model and discuss the main localization properties of the single-particle model. In the same section, we also review the known results related to the interacting AA and Fibonacci models. In Section III, the main results of this paper are presented; namely we (i) describe our results for the Fibonacci model at two different interactions strengths, (ii) obtain the many-body phase diagram of the i-IAAF model, and (iii) discover the anomalous impact of manybody interactions on the regime where the model exhibits mobility edges. In Section IV, we discuss in detail our results from Section III, in particular, the physical mechanisms behind the many-body phase diagram of the i-IAAF model. Our findings are summarized in Section V. Technical details are relegated to Appendices.

II. MODEL
We consider a finite spin-1/2 chain with open boundary conditions containing L sites in a quasiperiodic magnetic field. The Hamiltonian of the system is given by where S x,y,z j denote standard Pauli matrices, the in-plane coupling amplitude J = 1 sets the energy scale of the problem, and we work in units with = 1. In this parametrization, ∆ = 1 corresponds to the isotropic Heisenberg chain. The chain is subjected to a spatiallymodulated magnetic field h j = λV j (β) of strength λ ≥ 0, where the modulation function V j (β) is defined as The function V j (β) is quasiperiodic [see Fig. 1(a)] with an irrational spatial modulation frequency taken as the inverse of the golden mean, b = 2/(1 + √ 5). The parameter φ ∈ [0, 2π] acts as a global translation of the spatially-modulated potential. The tunable parameter β provides a knob by which we can interpolate between the two known limiting cases: (i) β → 0 yields the Aubry-André cosine modulation [13,14], up to a constant shift in energy, and (ii) β → ∞ corresponds to a step function switching between ±1 values according to the Fibonacci sequence [15,16]. The interpolating function (2) was introduced in Ref. [31] for proving the topological equivalence of The AA and Fibonacci models. The advantage of such an interpolation is its simplicity, more precisely, it contains a single tuning parameter which then transforms a cosine function into a steplike function, as shown in Fig. 1(a). One could think, in principle, of other interpolating functions, but the main thing to keep in mind is reaching the correct limits of the AA and Fibonacci model. For other smooth and continuous interpolations, we do not expect the main features of the model to change dramatically.
Applying the Jordan-Wigner transformation, our model (1) is mapped to an interacting one-dimensional model of spinless fermions where c † j and c j are fermionic creation and annihilation operators, n j = c † j c j is the local density operator, and t = −J/2 is the hopping amplitude. The spatiallymodulated magnetic field is transformed to a modulated on-site potential, and the term proportional to ∆ ≥ 0 acts as a nearest-neighbor repulsive interaction. We dub this model the interacting interpolating Aubry-André-Fibonacci (i-IAAF) if ∆ = 0, and IAAF in the noninteracting case.
Recently, the IAAF model was shown to exhibit interesting localization properties [19], which we briefly outline in Section II A. In this work, we concentrate on the many-body localization (MBL) properties of the i-IAAF (3) as it interpolates, with β, between the Aubry-André and Fibonacci limits, see Fig. 1

(b).
A. Noninteracting IAAF model Before considering the full interacting model (3), we summarize the salient localization properties of its noninteracting version (∆ = 0) [19], see Fig. 2. At the β = 0 limit, the noninteracting AA model is known to exhibit a pure spectrum localization transition at [13] λ C /t = 2, where all eigenstates of the model change from extended to localized. In the opposite limit, β → ∞, the noninteracting Fibonacci model exhibits critical behavior [15,48,49], namely the eigenstates decay in space as an inverse power law (critically localized eigenstates) for any value of λ.
Along the deformation (2) as a function of β, the IAAF potential becomes steeper, thus producing a stronger forcing on particles moving in such a potential. As a result, the region, where the eigenstates are purely extended, shrinks. Furthermore, starting from a localized AA phase, where every state is localized on a single site, and upon increasing β, the potential (2) pushes the states of the lowest energy band to resonance. For a chosen λ/t = 8, the first resonance occurs at β ∼ 2, where the states hybridize and become extended, leading to a mobility edge scenario, which is peculiar for MBL physics [33,50,51]. Increasing β further to β > 2, the states in the lowest band localize once again. This resonance cascade repeats for higher β until all states hybridize to form the critical states of the Fibonacci chain. In Fig. 2(b), we show the cascade for the ground state. Note that the delocalization due to hybridization followed by localization involves a doubling of the localization length, and it does not occur simultaneously for all states in the spectrum.
To visualize the aforementioned localization properties of the IAAF model, we numerically diagonalize a finite chain and evaluate the inverse participation ratio (IPR) for each eigenstate of the model (3). The IPR of an eigenstate ψ(E n ) with eigenenergy E n is given by In the regime where the n-th eigenstate ψ(E n ) is extended, the IPR is proportional to the inverse of the system length, i.e., IPR(E n ) = 1/L which approaches zero in the thermodynamic limit. On the other hand, if the n-th eigenstate is exponentially localized on N sites, its IPR is equal to 1/N and remains finite even for an infinite system size L → ∞.
In Fig. 2(a), we plot the IPR of every eigenstate as a function of β for a constant λ/t, and mark the appearance of delocalized states in the lowest band-bundle of the spectrum. In Fig. 2(c), we show the model's localization phase diagram for the lowest band-bundle of the energy spectrum, marked in Fig. 2(a). We use the averaged inverse participation ratio, IPR = (1/m) m n=1 IPR(E n ), whose value is equal to 0 when all eigenstates are extended, and 1, when all of them are localized on one site. The cascade of localization-delocalization transitions is visible, although not as clearly as for each eigenstate separately [19]. The strongest delocalization transition happens around β ∼ 2 for constant λ/t = 8. Furthermore, the localized region for 2 β 8 has a large fraction of the states localized on two neighboring sites, which follows from the fact that IPR ≈ 1/2. Before presenting the main results of this work, we discuss relevant MBL results related to the many-body Aubry-André and Fibonacci models. We also use this opportunity to summarize the computational methods that are often employed to study MBL in disordered or in quasiperiodic models. We first discuss numerical methods and emphasize their advantages and disadvantages.
We start with exact diagonalization (ED) [20]: its main advantage is that it is exact, providing solutions for eigenenergies and eigenstates of interacting Hamiltonians. The exact nature of ED comes at a price: only very small systems (∼ 24 spin-1/2 sites) are accessible and the outcome is riddled with finite-size effects, rendering the estimation of the system's behavior in the thermodynamic limit difficult to extract. Using ED, previous works investigated signatures of MBL in the behavior of the one-body density matrix [21], statistics of the many-body spectrum [8], entanglement entropy [12,[52][53][54], as well as the dynamics of entropy and imbalance after quenching the system [12,55].
Another class of numerical methods involves a description using matrix product states (MPS), which are variational Ansätze for the many-body wave function. MPS approximate many-body states by truncating the entanglement spectrum, which works especially well in one dimension due to the zero-dimensional boundary for entanglement's area law scaling. There are several MPS-based algorithms for the dynamics of quantum systems [56], such as the time-dependent density matrix renormalization group (tDMRG) [57,58]. Here, we will use the time-dependent variational principle (TDVP), as recently generalized to matrix product states [59,60]. Using this method, one can study large systems with O(100) sites [61], which is crucial for exploring MBL phases. An alternative MPS-based approach [62] employs a boundary-driven Lindblad equation. This approach is especially efficient for studying dynamical states close to pure thermal (and hence highly-entangled in the S z basis) states. Unfortunately, it is therefore less suited for strongly localized and weakly entangled systems. Hence, it was used to study the interacting Fibonacci model only in the limit of relatively weak fields [25].
We now outline relevant MBL results for the manybody AA and Fibonacci models, starting with the AA model. In Refs. [20][21][22], ED was used to analyse relatively small chains with ∼ 20 sites. Initially, an MBL transition was identified for ∆ = 1 at an estimated critical AA potential strength [20] of 2 λ C /t 5. In later works, by studying the spectrum of eigenvalues of the one-body density matrix, an MBL transition at λ C /t ≈ 4 was obtained [21] for ∆ = 1, while from the statistics of the spectrum of eigenfunctions of the Hamiltonian the transition at λ C /t 3 was inferred [22]. In Ref. [23], one of us obtained a critical AA potential strength of λ C /t ≈ 4.8 from the dynamics of a quenched imbalance using TDVP. The analysis was performed on long chains with L = 50 sites. The interacting Fibonacci model was first explored in Refs. [63,64]. Combining analytical methods, such as bosonization and renormalization group, the lowtemperature properties of the model were studied. It was found that the model shows anomalous transport properties with a scaling exponent depending on the interaction strength and the position of the Fermi level. The hightemperature regime, which we consider in this paper, was recently studied in Refs. [24][25][26]. It was shown [24] that the interacting Fibonacci model with ∆ = 1 undergoes an MBL transition at 4 λ C /t 7. Specifically, ED with exact Krylov-space approach was applied to chains of length up to 24 sites, and both static probes, e.g., spectral statistics and scaling of the entanglement entropy, as well as dynamic ones, e.g., entanglement growth and evolution of local observables after a quench, were investigated. Further results [25] using tDMRG, involved boundary-driven Lindblad dynamics, as well as unitary dynamics, suggesting that for ∆ 0.5, the model is diffusive for all potential strengths, while for ∆ 0.5 there exists an interaction-induced subdiffusive regime that persists at least up to λ/t = 3. Recently, by studying the Rényi participation entropy and the occupation number at half chain using ED, it was proposed [26] that even at a weak interaction (∆ ≈ 0.5) an MBL phase should occur within the interval 4 < λ C /t < 8.
Many-body interactions bear a somewhat opposite impact on the AA and Fibonacci models. On one hand, in the AA model, the interactions tend to delocalize the system and shift the localization transition towards higher values of λ/t compared with the noninteracting case [22]. This trend can be qualitatively understood as follows: (i) recall that the localization transition in the noninteracting limit is homogeneous throughout the whole spectrum, and occurs at λ/t = 2, where the model is self-dual; (ii) the self-duality condition is fragile and is broken by the interactions, even at the mean-field level [50]; (iii) a mobility edge thus appears, i.e., around λ/t = 2, the spectrum of the weakly interacting AA model contains both extended and localized states; and (iv) the extended states can act as an effective bath to the localized states and delocalize them [33,34].
On the other hand, in the Fibonacci model, the interactions destroy critical states and localize the system at large enough λ/t. One possible explanation for such an effect [24] involves the Fourier components of the Fibonacci potential [β → ∞ limit of Eq. (2)], namely, V (k) = L j=1 V j e −ikjb ∼ 1/k. The slow decay of the potential is responsible for the model's noninteracting critical eigenstates, which decay as a power law in real space. The interactions on mean-field level introduce new terms in the potential and therefore change the Fourier components to [24] V (k) ∼ 1/k α with α > 1. Quasiperiodic models with such a Fourier space behavior exhibit localization-delocalization transitions [65].

III. RESULTS
To explore the localization properties of the i-IAAF model (3), we examine quench dynamics. As an initial state, we take a Néel state: |ψ(t = 0) = |↑, ↓, . . . , ↑, ↓ , which is in the j S z j ≡ j ψ|S z j |ψ = 0 spin sector [Eq. (1)], and corresponds to half-filling in the particle picture [Eq. (3)]. The Néel state resembles a mid-band Bloch state that thermalizes quickly in the absence of disorder [66]. Given the initial state, we employ the timedependent variational principle (TDVP) [59,60] to simulate the time evolution of relatively-long chains (L = 50 throughout this section). Our approach is similar to that of recent works [23,55,61,67,68].
The observable that we concentrate on is the spin imbalance defined as where τ is the simulation time in units of J −1 . It quantifies the local memory of the Néel initial state; the value of the imbalance at the initial time is I(τ = 0) = 1. For a system in a localized phase, the imbalance saturates to a constant value at long times. If the system is in the ergodic phase, the imbalance vanishes in the long-time and thermodynamic limit. In the case of random on-site disorder (corresponding to Anderson localization in the noninteracting limit), the imbalance averaged over many disorder realizations decays as an inverse power law close to the MBL transition point [69] where the bar denotes disorder averaging and the exponent γ depends on the strength of the disorder. Deep inside the MBL phase, memory of the initial state persists at all times and therefore γ → 0.
In the quasiperiodic case, due to the absence of Griffiths effects, which depend on the appearance of rare regions, the decay does not exactly obey a power-law form with subdiffusive exponents at late times. Nonetheless, dynamics in the localized regime will still saturate to a finite imbalance, resulting in a vanishing γ. Note that the counterpart of disorder averaging in the quasiperiodic case is averaging over the phases φ.
This section is organized as follows: in Section III A, we report our results for the interacting Fibonacci model [the β → ∞ limit of Eq. (3)] and compare them with contemporary literature on this limit [24][25][26]. In Section III B, we obtain the phase diagram of the i-IAAF model with ∆ = 1 and ∆ = 0.3, and in Section III C, we concentrate on a specific region in the phase diagram, where anomalous localization-delocalization behavior occurs by tuning the interaction strength from ∆ = 0 to ∆ = 1.  [70].

A. Fibonacci model
We start by revisiting the interacting Fibonacci model, i.e., the β → ∞ limit of Eq. (3). We calculate the time evolution of the imbalance (4) for several values of the potential strength λ/t, and show the result in Fig. 3.
Strong interactions -The ∆ = 1 case is shown in Figs. 3(a) and (c). In Fig. 3(a), we observe that the imbalance exhibits strong oscillations as it evolves in time, showing revivals associated with short-scale oscillations of the particles within their close environment. At the same time, a clear distinction manifests in the long-time behavior, where, for λ/t ≈ 2.2, the imbalance decays roughly as an inverse power law, while for λ/t ≥ 3 it saturates to a constant value. We fit the long-time trend with an inverse power law behavior [cf. Eq. (5)], and obtain the exponent γ, see dashed lines in Fig. 3(a). Repeating this procedure for several values of λ/t, we plot the fitted γ in Fig. 3(c). Initially, the exponent decreases with increasing λ/t. This decrease saturates at λ/t ≈ 3, while at λ/t ≈ 4 the imbalance does not appear to decay on the computationally-available evolution timescales. Hence, we estimate that an MBL transition occurs within the interval 3 ≤ λ C /t < 4. As a further confirmation of our finding, in Appendix A, we show similar results for a numerically exact quench of a shorter chain with L = 16 sites and bond dimension χ = 256.
Weak interactions -We now turn to weaker interactions, in order to compare with the results presented in Ref. [25]. It was suggested that the Fibonacci model is diffusive for λ/t ≤ 3 and at low interaction strengths, ∆ < 0.5, while at ∆ = 0.5 a subdiffusive phase appears, but no MBL phase was predicted. In Figs. 3(b) and (c), we report the outcome of our numerical quench for a low-interaction strength ∆ = 0.3. Here, the imbalance decays faster and exhibits slower oscillations relative to the former ∆ = 1 case. Moreover, the appearance of an MBL phase occurs only for much higher λ/t values, where even for λ/t = 4.6, we observe a slow decay. In Fig. 3(c), we show the fitted power-law exponent as a function of the potential strength. A clear memory of the initial state remains for λ C /t 4.75. Comparison with exact numerics for a chain of L = 16 sites and χ = 256 (see Appendix B) yields a similar value of λ C /t 4.5. For even weaker interactions, i.e., ∆ = 0.1, we do not observe signatures of an MBL phase transition for potential strengths λ/t ≤ 10.
Source of oscillations -We now comment on the oscillations present in Figs. 3(a) and (b). Note that, unlike true disorder or the AA case, the amplitude of the oscillations does not reduce with averaging over a large number of different quasiperiodic potential realizations. The reason is that the Fibonacci potential has only L/2 possible configurations [24], where L is the number of sites in the system. A possible way to reduce the amplitude of the oscillations is to use random product states as the initial state instead of the Néel state used in this work, which goes beyond the scope of this work. Instead, we explore here the behavior of these oscillations stemming from the initial Néel state. Note that the persistent oscillations and the memory of the initial state can lead to negative values of γ, similar to what is observed in the AA model [23].
The oscillations are irregular for a weak potential, i.e., in the ergodic phase, while for strong quasiperiodic potentials they become more regular with an approximately constant period. Such pronounced oscillations are due to the shape of the potential, which contains many occurrences of nearest-neighbor pairs with degenerate on-site energies. In the initial Neél state, there is only one particle populating any such pair of states. That particle will dominantly hop between the nearest neighbor degenerate states and produce oscillations in the imbalance. This simplified picture is valid for strong potentials and weak interaction, and breaks in the opposite limit of weak potentials and strong interaction.
The periods of oscillations extracted from the averaged imbalance in the localized phase [see Fig. 3(a) and (b)] are T ∆=1 ≈ 2π/ √ 2J −1 and T ∆=0.3 ≈ 2πJ −1 for ∆ = 1 and ∆ = 0.3 interaction strengths, respectively. From the periods it is possible to extract an effective hopping t ∝ 1/T , see Ref. [71] for a full derivation of the imbalance in the noninteracting and clean case. The ratio of effective hopping from Figs. 3(a) and ( To explain the difference between two effective hoppings, in Appendix C we present a simple model that describes how interactions change the particles dynamics within the aforementioned degenerate nearest-neighbor pairs. This simple analytical model yieldst = t 1 + (∆/J) 2 , which agrees well with the numerical results from Figs. 3(a) and (b). For the case of weak interactions [ Fig. 3(b)], the effective hopping is given byt = t, i.e., the oscillations appear to stem from the single-particle dynamics, while in the case of strong interactions [ Fig. 3(a)], the effective hopping is renormalized tot ≈ √ 2t, suggesting that interactions increase the effective kinetic energy of the particles. To conclude this discussion, the oscillations in the MBL phase indeed stem from an ensemble of effective particles hoppings between two degenerate neighboring sites.

B. Interacting IAAF model
We now address the i-IAAF model with strong and weak interaction, ∆ = 1 and ∆ = 0.3, respectively. For the strong interaction, the MBL transition in AA and Fibonacci models occurs at different potential strengths λ C /t ≈ 5 and λ C /t ≈ 4, respectively. Note that in this case, the Fibonacci MBL transition occurs for a lower λ C /t than that of the AA. Interestingly, for the weak interaction case, the order inverts and we find in the AA case λ C /t ≈ 4 while for the Fibonacci λ C /t ≈ 4.75. It is, therefore, our goal to explore how the two transition points are connected once β is tuned from AA (β = 0) to the Fibonacci model (β → ∞), for both weak and strong interactions, cf. Fig. 1(b). To answer these questions, we again analyze the dynamics of the averaged imbalance, see Eq. (5) and Fig. 4(a) for the fitted exponent γ as a function of λ/t and β.
For ∆ = 1, starting from the MBL phase of the AA model λ/t ≥ 5 and increasing β, γ always remains around 0 and no transition to an ergodic phase is observed. On the other hand, starting from the ergodic phase of the AA, an MBL transition takes place for higher values of β. The value of λ C /t, where the transition occurs, decreases as a function of β. Such a behavior is reminis- cent of the noninteracting case, where the impact of the potential becomes more pronounced with β, and the region of extended states shrinks, see Figs. 2(b) and (c). Note, however, a crucial difference between the noninteracting and interacting cases is the absence of localizationdelocalization transitions with changing β above the critical localization point of the AA model. For ∆ = 0.3, at β = 0 and λ/t 5, the situation is similar to the strong-interaction case: the exponent always remains around zero and no transition to an ergodic phase is observed, see Fig. 4(a). On the other hand, starting from the localized phase of the AA model at 4 ≤ λ/t < 5, a transition to the ergodic phase is observed at the value of λ/t that depends on β; the transition point moves to higher values of λ/t with increasing β. Note that we explore only up to β = 6 in the phase diagram for weak interaction. For β ≥ 6, the decay of the imbalance with time does not change, but the longperiod oscillations make it difficult to properly determine the exponent, cf. Fig. 3(b), where similar fitting challenges appear in this limit. For example, the imbalance for both λ/t = 4.6 and λ/t = 7.8 in Fig. 3(b) saturates with additional short-period oscillations, but the value of the fitted exponent is not the same, see Fig. 3

(c).
Despite these technical challenges, we clearly observe that the interaction shifts the critical point toward higher values of λ/t in the AA limit, and to lower values in the Fibonacci limit. Furthermore, interpolating between the two models (i.e., for finite β), the critical points at the two limits (AA and Fibonacci) are connected with a continuous line separating the ergodic from the localized phase. As a result, we conclude that MBL transitions will occur also with changing ∆ in the intermediate regions of the interpolation.
Similarly to the Fibonacci case in the previous section, we comment briefly on the oscillatory behavior of the imbalance. While the oscillations are pronounced even at large times in the Fibonacci limit [see Figs. 3(a) and (b)], for finite β, their amplitude decays with time, see Fig. 5(b). The reason lies in the shape of the potential. As previously discussed, in the Fibonacci limit, the potential (2) allows for the appearance of pairs of neighboring sites that have the same on-site value [cf. Fig. 1(a)]. Specifically, the Rabi oscillations of particle hopping between two such sites, which appear many times throughout the chain, give rise to coherent oscillations in the imbalance that persist even at long times. In comparison, for finite β, the degeneracies of the potential characteristic of the Fibonacci limit are lifted, and, hence, particles hopping between the pairs of sites will exhibit different Rabi oscillations, which results in a reduced amplitude of oscillations in the overall imbalance. Nonetheless, there are still oscillations that are longer-lived than in the case of purely random disorder, since certain values of the potential differences between neighboring sites are more likely than others [23].

C. Anomalous localization by interaction
After analyzing the overall phase diagrams of the weakly and strongly interacting IAAF model, we now concentrate on the noninteracting delocalized-phase sliver at λ/t > 2 and β ∼ 2, and investigate its evolution when the interaction strength is increased from ∆ = 0 to ∆ = 1. We pick a point from the aforementioned region with (λ/t, β) = (4, 1.4), see the red dot in Fig. 5(a). This point lies in the ergodic phase of the strongly interacting (∆ = 1) model, cf. Fig. 4(a). Fixing the parameters λ/t and β, we calculate the time evolution of the averaged imbalance for three different interaction strengths, ∆ = 0, ∆ = 0.3, and ∆ = 1, see Fig. 5(b).
In the noninteracting and strongly interacting regimes, the imbalance decays quickly. Interestingly, for weak interactions, the decay is noticeably slower. In Fig. 5(c), we plot the exponent (5)) as a function of interaction strength, and find that it indeed does not change monotonously with ∆. Specifically, by increasing the interaction strength, the system first localizes for a weak interaction and then delocalizes at a stronger interaction,  (5), as a function of interaction strength ∆ for several points in the phase diagram [also marked in (a)]. We pick points that lie within the delocalized-phase sliver of the noninteracting model, i.e., situated at finite β ∼ 2 and at λ/t > 2. System size in all plots is L = 50, the number of different realizations of the quasiperiodic potential is 36, and the inverse power law is fitted in a time window [10,200]. The bond dimension used for all plots is χ = 64 and the time step is δτ = 0.1. Error bars are 1σ intervals based on a bootstrapping procedure.
which is indicated by a noticeably larger γ for ∆ = 0 and ∆ = 1, compared to intermediate values of ∆.
We repeat this procedure for two additional points from the noninteracting delocalized sliver, see the blue and green points in Fig. 5(a) and the result in Fig. 5(c). As expected, for ∆ = 0, the exponent γ is large, reflecting the delocalized nature of the noninteracting system, see Fig. 2(a). At finite interaction strengths, the exponent reduces in value and after a critical interaction strength ∆ C remains constant. This indicates that the system's dynamics slows down with increasing interaction strengths, until the system becomes localized.

IV. DISCUSSION
Let us now collect all the observations from Section III and discuss them from a broader perspective, focusing on qualitative mechanisms that are responsible for the localization properties of the i-IAAF model.
In the two limiting cases of AA and Fibonacci models, the localization properties have been discussed in the literature in great detail for both single-particle and manybody cases. The crucial difference between the AA and Fibonacci models is the shape of the potential modulations: in the former, the onsite energies are distributed in the interval [−λ, +λ], whereas in the latter, they are highly degenerate with only two values ±λ. This leads to a localization transition in the noninteracting AA model and a critical spectrum in the Fibonacci model. Since the single-particle properties of the two models are very dissimilar, many-body interactions introduce distinct effects in them.
As already mentioned in Section II B, in the AA model, interactions shift the localization phase transition towards higher values of λ/t due to dephasing, similar to the interacting Anderson model with random disorder. On the other hand, in the Fibonacci limit, the interactions destroy the fragile criticality caused by the numerous degenerate on-site terms, and introduce a localization transition. This can be understood qualitatively on a mean-field level, where interactions introduce new onsite terms, i.e., shift the onsite energies, and break the degeneracy of the potential. Once the degeneracy is lifted, the hopping of particles is suppressed, and the system becomes localized for large enough λ/t. It is, therefore, expected that by increasing the interaction strength ∆, the system will localize more easily, i.e., for lower values of λ/t. Such behavior is observed in Fig. 3. Now, we turn to the i-IAAF model. In its noninteracting limit, the model shows rich localization properties with a cascade of localization-delocalization transitions, as shown in Fig. 2. Let us concentrate on the first delocalization transition at finite β. Starting from a localized phase in the AA limit, the deformation of the potential from a cosine to a step-like function with increasing β brings the low-energy states into resonance. Once in resonance, the states hybridize for any finite hopping and delocalize. As shown in Fig. 2(a), all states in the lowest band-bundle are delocalized, and the singe-particle spectrum has a nontrivial mobility edge. The nontrivial mobility edges repeat at different positions in the spectrum for higher β. It is important to note that the width of the delocalized band-bundle shown in Fig. 2(a) reduces with increasing λ/t, since the hopping between the sites is effectively reduced. As a consequence, the delocalized regions (slivers) in the phase diagram shrink at higher λ/t.
Adding interactions alters the localization properties of the model as shown in Figs. 4 and 5. There are two main features of the presented many-body phase diagrams: (i) the slope of the MBL phase boundary changes sign with the strength of many-body interactions, and (ii) the delocalized slivers of the noninteracting model disappear once interactions are introduced. Point (i) is expected since interactions have opposite effects in AA and Fibonacci limits, as discussed above. On the other hand, point (ii) is surprising since one would naively expect that the existence of delocalized states in the spectrum of the noninteracting limit will act as a bath to localized states once interactions are turned on, and will eventually delocalize the whole system [33,34]. However, as already mentioned, such delocalized states at finite β are the consequence of carefully tuned degeneracy between onsite energies imposed by the deformation of the potential (2), see Figs. 1(a) and 2(a). Similar to the Fibonacci case, interactions can easily lift the degeneracy and localize these slivers.
Indeed, we recall that the width of the delocalized band-bundle reduces with increasing λ/t, and, therefore, at higher λ/t, it should feel stronger effects of interactions, i.e., the delocalized states should localize for a smaller interaction strength ∆. Such an effect is observed in Figs. 5(b) and (c), namely, at the points with higher λ/t the decaying coefficient is smaller for all values of interaction strength, indicating stronger localizing effect of interactions for higher λ/t. Furthermore, at some points of the phase diagram that lie inside the delocalized sliver, we find an anomalous behavior with tuning of the interaction strength, where weak interactions tend to localize the system, while stronger interactions delocalize it again, see (λ/t, β) = (4, 1.4) point in Fig. 5. The localization at weak interactions can be explained with the degeneracy-lifting argument above, while the delocalization could follow from the same dephasing mechanism that is present in the AA limit. Further studies are needed to better explain and describe this peculiar phenomenon.

V. CONCLUSION
To conclude, we have investigated the many-body version of the IAAF model, focusing on the half-filling sector of the model and numerically studying the time evolution of the spin imbalance using the TDVP method [59,60]. Our main criterion for localization is the saturation of the averaged imbalance at long times. In the ergodic phase, the imbalance decays over time. We have observed that for quasiperiodic models with weak interaction, the high degeneracy of the on-site potentials leads to strong oscillations with time in the averaged imbalance. It is, therefore, more difficult to locate a precise MBL transition point in such models. Nevertheless, we have constrained the transition regime, and thus explored the localization-delocalization phase diagram of the interacting IAAF model.
Our work contains three main results. First, we have shown that for strong interaction, ∆ = 1, the MBL phase transition in the Fibonacci model (β → ∞ limit of the IAAF model) is located at 3 ≤ λ/t < 4, which is lower than previously reported values [24]. Furthermore, we have demonstrated that the MBL transition can survive even for small interactions of ∆ = 0.3, and thus complemented the results of Ref. [25]. Secondly, we have presented the many-body phase diagram for the IAAF model, which exhibits a monotonous evolution of the MBL phase when going from the AA (β = 0) to the Fibonacci (β = ∞) limit. In other words, the cascade of delocalization transitions in a single-particle case [19] is destroyed for both weak and strong interactions. Lastly, we have observed an anomalous behavior when increasing the interaction strength from ∆ = 0 to ∆ = 1 in the parameter sliver where part of the spectrum is delocalized in the noninteracting limit. In other words, in the sliver both extended and localized states exist in the system, leading to a situation where weak interaction tends to localize the system, while stronger interaction drives it to the ergodic phase.
Future work will focus on energy densities away from half-filling. By studying the model only at the halffilling, we were not able to conclude if the localization properties are homogeneous throughout the spectrum, or whether an interesting cascade structure appears, as in the single-particle case. Moreover, the present study was performed at an infinite temperature, and the behavior with finite temperatures remains unknown. Therefore, one promising direction would be to analytically investigate the i-IAAF model at finite temperatures and for different chemical potential using, e.g., bosonization and renormalization-group approaches [63,64].
To further confirm the results from Fig. 3 from the main text, in Fig. 6, we show the averaged imbalance and the fitted exponent γ for a short chain of L = 16 sites and large bond dimension χ = 256. In this case, the calculation is exact, since, for a system of L sites, non-truncated MPS have a maximum bond dimension of χ = 2 L/2 = 256. The behavior of the exponent in Fig. 6(b) is essentially identical to the one presented in Fig. 3, where we used larger chains with truncation of the bond dimension.

Appendix B: Convergence of the numerical results
In this Appendix, we verify the convergence of our numerical results, namely the convergence of the obtained imbalance decay exponent, which we use to map out the phase diagrams in Fig. 4. The dynamics of the imbalance is calculated using the TDVP method, whose long-time precision depends on the chosen bond dimension χ. Since the truncation error grows with time, the bond dimension determines the maximal time below which the calculated imbalance does not differ much from its true value. Larger bond dimensions are required in the ergodic phase compared to the MBL phase to reach the same maximum time, see Refs. [59,60,55,23,67] for more details.
In our case, we use the same bond dimension for the whole phase diagram [see Fig. 4(a)], implying that we need to change the upper boundary of the fitting time interval as a function of the parameters in the phase diagram to avoid large numerical errors. To do so, we choose a few different λ/t cuts in phase diagrams shown in Fig. 4(a), and calculate the exponent for larger bond dimension of χ = 128. We assert that the result has converged when for a given fitting time-window, the exponents γ in the χ = 64 and χ = 128 cases do not differ more than the error bars from one another.
The comparison of our results for the two different bond dimensions is presented in Fig. 7. For the strong interaction, results with different χ converge for a fitting time-window [50,120] used in Fig. 7(a), and [50,200] used in Fig. 7(b). In the case of weak interaction, the convergence is reached for a fitting time-window [50,210] in Fig. 7(c), and [50,230] in Fig. 7(d). We assume that if the results converge at λ /t for a certain fitting time-window, then for all λ/t > λ /t the same fitting time-window will lead to a similar convergence. Furthermore, we assume that the convergence does not strongly depend on β, as seen in Fig. 7.
Appendix C: Imbalance oscillations in the Fibonacci limit Here, we describe the mechanism behind the strong oscillations of the imbalance mentioned in Section III A. In the Fibonacci limit, the potential modulation will contain only two values, namely ±λ, cf. Fig. 1(a). We refer to a site with on-site potential −λ (+λ) as A-site (B-site). Note that the Fibonacci sequence allows for two neighboring sites of type A, but not for two B-sites. Therefore, the combination BAAB is possible (see Fig. 8(a)), and it will occur many times throughout the chain. This configuration we dub as "a well". In the initial Néel state, there will always be two particles inside the BAAB subsystem, one at a B-site, and the other inside the well, see Fig. 8(a).
Let us now consider a system that is deeply within the MBL phase, namely when λ/t 1. The particles at isolated A-sites and at all B-sites remain localized for long times, while particles residing in their respective wells can hop from one A-site of the well to the other. To see this effect, it is sufficient to concentrate on the three sites marked with the dashed box in Fig. 8(a). Such hopping Sketch of the mechanism responsible for strong oscillations in the Fibonacci limit. A particle, initially placed in a well formed by two "A" sites (with the lower potential −λ) and bounded by "B" sites (with the higher potential +λ), hops inside the well, giving rise to the oscillations of the imbalance, see Figs changes the value of the imbalance, and is not influenced by the value of λ. Furthermore, because the particle at site B remains localized, we can write the Hamiltonian for a single particle inside the well, and incorporate the interaction with the particle at the B-site through a renormalized on-site energy, see Eq. (C1) below. The basis we choose is { |10 , |01 }, where the first (second) state represents the particle sitting on the left (right) Asite of the well. The Hamiltonian in such a basis reads with the eigenenergies (C2) Initially, the particle is placed in the left A site of the well, and the initial state is ψ 0 = |10 . We can now evolve the system state with time τ according to |ψ(τ ) = exp(−iH AA τ ) |10 , and obtain |ψ(τ ) =e −i ∆ 2 τ cos 1 + ξ 2 J 2 τ with ξ ≡ ∆/J. The expectation values for the densities on the left and right A-site of the well are readily obtained as ψ(τ )| n left |ψ(τ ) = cos 2 1 + ξ 2 J 2 τ , ψ(τ )| n right |ψ(τ ) = 1 1 + ξ 2 sin 2 1 + ξ 2 J 2 τ . (C5) The particle imbalance is then calculated as I(t) = n left − n right = 1 (1 + ξ 2 ) ξ 2 + cos 1 + ξ 2 Jτ .
This simple exercise leads to two important conclusions. First, it follows that the period of oscillations is given by namely, it reduces with increasing interaction strength. In Fig. 8(b), we show the comparison of Eq. (C7) with the period of oscillations from numerically calculated imbalance. The period obtained from numerics agrees well with the analytical prediction (C7). Second, the amplitude of oscillations decreases with increasing ξ until they completely vanish for infinite interaction strength, i.e., lim ξ→∞ I(t) = 1. Such an effect can be seen in Figs. 3(a) and (b) and Fig. 6(a).

Appendix D: Mean gap ratio
To complement our conclusions in Section III B regarding the shrinking of the ergodic phase with increasing β in the strongly interacting case, we use ED on shorter chains, and study the spectral gap statistics of the system. Specifically, the gap ratios measure the energy level repulsion [8], and provide us with a measure for the localization in the system. The gap ratio r n between two consecutive gaps is defined as 0 ≤ r n ≡ min[δ n , δ n−1 ]/ max[δ n , δ n−1 ] ≤ 1, where δ n = E n+1 −E n ≥ 0 and n labels the eigenenergies from low to high.
Taking the average over all gap ratios in the spectrum and over 36 different realizations of the quasiperiodic potential, we obtain the mean gap ratio r, see Fig. 9. Inside the ergodic phase, energy levels follow the Wigner-Dyson statistics in the thermodynamic limit, and the mean gap ratio approaches the value [73] r WD ≈ 0.53. In the MBL phase, they follow Poisson statistics and saturate at the lower value of r P ≈ 0.39. Although it is difficult to precisely locate the value of the critical disorder λ C /t from r in such short chains, the overall behavior clearly shows that for higher β, all points are shifted towards smaller values of λ/t values. This is in agreement with the TDVP results in Figs. 4(a), namely, that the ergodic phase shrinks with the increasing of β.