Multifractality of correlated two-particle bound states in quasiperiodic chains
MMultifractality of correlated two-particle bound states in quasiperiodic chains
Diana Thongjaomayum, Sergej Flach,
1, 2 and Alexei Andreanov
1, 2 Center for Theoretical Physics of Complex Systems, Institute for Basic Science (IBS), Daejeon, Korea, 34126 Basic Science Program, Korea University of Science and Technology (UST), Daejeon, Korea, 34113 (Dated: March 3, 2020)We consider the quasiperiodic Aubry-Andr´e chain in the insulating regime with localised single-particle states. Adding local interaction leads to the emergence of extended correlated two-particlebound states. We analyse the nature of these states including their multifractality properties. Weuse a projected Green function method to compute numerically participation numbers of eigenstatesand analyse their dependence on the energy and the system size. We then perform a scaling analy-sis. We observe multifractality of correlated extended two-particle bound states, which we confirmindependently through exact diagonalisation.
I. INTRODUCTION
Understanding the transport properties of quantumdisordered or inhomogeneous systems has been an activetopic of research since the discovery of Anderson localisa-tion (AL). AL describes the arrest of transport in a singleparticle system due to disorder or inhomogeneous poten-tial which renders all the eigenstates in one and two spacedimensions exponentially localised. The original work ofAnderson triggered a sequence of theoretical studies andby now the single particle case is well understood. Theimportant and much harder question is the stability ofmodification of AL in the presence of many-body inter-actions. Decades of research attempts culminated in theopening of the field of Many Body Localisation.
Inter-estingly one of the strongly debated issues is the possibleexistence of ’bad’ metallic states which are non-ergodicor simply multi-fractal.
A notorious issue with MBL related studies is the com-putational complexity due to the exponential prolifera-tion of the Hilbert space dimension with increasing num-bers of particles and system size. A legitimate and com-plementary approach is therefore to consider only fewinteracting particles, which allows to increase the systemsize beyond the limits set by typical MBL models. Threemain directions with single particle localisation in one di-mension have been explored: genuine AL due to uncor-related disorder, Wannier-Stark (WSL) localisation dueto an external dc field, and Aubry-Andr´e (AAL) local-isation due to a quasiperiodic external potential. Gen-uine AL yields a nontrivial increase of the localisationlength for two interacting particles with still unsettledscaling details.
Two interacting particles yield no lo-calisation change for WSL with interaction, only affectingthe Bloch oscillation periods. At variance, AAL withquasiperiodic potentials showed an unexpected transi-tion from localisation (zero interaction) to delocalisation(non-zero interaction). These findings were later con-firmed in Ref. 14 which provided additional indicationsfor the fractal nature of the delocalised eigenstates.Are these the seeds of a bad metal and the MBL tran-sition from above? A hint might be obtained from thestriking similarity of the phase diagram of correlated metallic two-particle bound states in Fig. 4 of Ref. 13and the phase diagram of an MBL phase which was ex-perimentally assessed for interacting fermions in opticalquasiperiodic potentials in Fig. 4 of Schreiber et al inRef. 15. In the present study we attempt to add moreconclusive arguments which aim at a positive answer forthe above question for quasiperiodic potentials. We con-firm the fractal character of the two-particle spectrumand the fractalilty of some of the two-particle states. Werely on the projected Green function method, originallydeveloped to analyse the localisation length of two inter-acting particles in the AL case. The paper is organisedas follows: we introduce the tools and other necessarymeans in Sec. II. Section III benchmarks these tools inthe single particle case against the exact results and exactdiagonalisation. In Sec. IV we analyse the two interactingparticles case. This is followed by conclusions. II. SETTING THE STAGE
The starting point is a single particle placed in aquasiperiodic potential with the Aubry-Andr´e Hamilto-nian H = (cid:88) n ( | n (cid:105)(cid:104) n + 1 + h.c) + (cid:88) m h m , (1) h n = λ cos(2 παn + β ) , where λ is the strength of the potential, α is an irra-tional number ensuring quasiperiodicity of the potential.We choose α = ( √ − /
2, the golden ratio and we fix thehopping strength t = 1. Depending on the strength of thepotential λ the eigenstates are all delocalised ( λ <
2) orlocalised ( λ >
2) with localisation length ξ = 1 / ln( λ/ Finally β is aphase which can be varied to generate different realisa-tions of the quasiperiodic potential. In numerical studieswith finite system size the choice of β will affect localisedand sparse, fractal or multi-fractal extended states. Inthe present study involving critical states we use averag-ing over different values of β , that we denote as · · · , toimprove statistics. a r X i v : . [ c ond - m a t . d i s - nn ] M a r We now add the interactions and consider two inter-acting bosons. We choose the onsite Hubbard interactionof strength u . The total Hamiltonian is given by H = H ⊗ H + uP = (cid:88) n,m ( | n, m (cid:105)(cid:104) n + 1 , m | + | n, m (cid:105)(cid:104) n, m + 1 | + h.c.)+ (cid:88) n,m | n, m (cid:105) ( h n + h m ) (cid:104) n, m | + uP, (2)where | n, m (cid:105) is a basis state with two particles at site n, m , and h n is the onsite Aubry-Andr´e potential at site n given by Eq. (1). P is the projection operator definedas P | n, m (cid:105) = δ nm | n, m (cid:105) that enforces the onsite Hubbardinteraction.The authors of the work Ref. 13 used exact diagonali-sation and unitary evolution of wavepackets to study thetwo-particle properties of the model (2). The exact di-agonalisation limited the largest system sizes achievableto N ≈ to handle largesizes, up to N = 10946, of the Hamiltonian (2). Wefollow the original approach of Ref. 16. We extract therelevant two-particles properties from the projected two-particle Green function, which is obtained as a projectionof the full Green function G = ( E − H ) − onto doublyoccupied states (relying crucially on the fact that theHubbard interaction is proportional to the projector P ):˜ G = ˜ G − u ˜ G . (3)Here ˜ G = P GP and ˜ G = P G P ; G is the non-interacting two particle GF which can be obtained bystraightforward diagonalisation of the single particleHamiltonian (1). Knowing the single particle eigenen-ergies { E µ } and eigenfunctions { φ µ ( n ) } we compute G as follows: (cid:104) n, n | G ( E ) | m, m (cid:105) = (cid:88) µ,ν φ µ ( n ) φ ν ( n ) φ ∗ µ ( m ) φ ∗ ν ( m ) E − E µ − E ν = (cid:88) µ φ µ ( n ) g ( E − E µ ) φ µ ( m ) , (4) g ( E ) = 1 E − H = (cid:88) ν φ ∗ ν ( n ) φ ν ( m ) E − E ν . The reordering of the terms in the second line is doneto reduce the complexity of the computation from theoriginal O ( N ) to O ( N ), since the single particle Greenfunction g can be efficiently evaluated using tridiagonalmatrix inversion of the single particle Hamiltonian (1).This approach allows us achieve system sizes as large as N = 7000.In the insulating regime the exponential decay of theprojected Green’s function ˜ G was used to extract the two interacting particles (TIP) localisation length. Here we are aiming to investigate TIP eigenstates whichwe expect to be extended in a predominantly insulatingregion, therefore ˜ G might not decay or the decay mightnot be exponential. Consequently we adopt a differentmeasure: Interpreting the projected Green function ˜ G as a probability density function we define the participa-tion number I q =2 and its higher moments I q> as I q = ( (cid:88) k | ˜ g ( k ) | ) q / (cid:88) k | ˜ g ( k ) | q , (5)where ˜ g ( k ) = < n, n | ˜ G | n + k, n + k > . We shall use I andhigher moments that are always well defined to analysethe TIP states. To distinguish I q from the conventionalparticipation number we will refer to it as the Green func-tion participation number (GPN). However before we canproceed to the two particle case, we need to confirm that I is a valid measure of localisation of an eigenstate Ψ,similar to the conventional participation number:PN q = (cid:88) nm | Ψ nm | q , (6)e.g. that I can distinguish between extended,(multi)fractal and localised states. III. SINGLE PARTICLE: BENCHMARKING
To confirm that the above defined participation num-ber I q is a valid probe of localisation properties of eigen-states we first consider the single particle case. To achievethis we benchmark two single-particle quantities: locali-sation length – analytical and numerical ξ = 1ln (cid:0) λ (cid:1) (7)1 ξ = − lim | n − m |→∞ ln |(cid:104) n | g | m (cid:105)|| n − m | , (8)and participation number I , which is defined similarlyto its two particle version Eq. (5): I q = ( (cid:88) k | g ( k ) | ) q / (cid:88) k | g ( k ) | q , (9)where g ( k ) = < n | g | n + k > . Our aim is to confirm that I is a valid substitute for ξ for localised states and itbehaves like the conventional participation number forlocalised and extended states.To prove that we consider 3 different values of the po-tential strength λ = 1 , , .
5, which correspond to de-localised, critical and localised regimes. For each λ wecompute single particle eigenstates ψ and the Green func-tion g for energies E ∈ [ − ,
3] in steps of ∆ E = 0 . .
05 is chosen to be slightly bigger than the levelspacing δ ( N = 250) = 0 .
004 and δ ( N = 500) = 0 .
002 for I , P N = 1.0, all states extended N=250, I N=500, I N=250, PN N=500, PN Eigenvalues (a) , I = 2.5, all states localised N=250, I N=500, I N=250, N=500, e )1 =4.48Eigenvalues (b) FIG. 1. (Colour online) Benchmarking of a single particle inthe AA model: (a) λ = 1: I and PN for N = 250 , λ = 2 .
5: Localisation length ξ and I for N = 250 , ξ = 4 .
48. The par-ticipation number PN behaves similarly to I (not shown).The bottom blue/circular points in both (a) and (b) representthe spectrum of H for N = 500 and show the locations ofthe eigenstates. the data presented in Fig. 1. From the Green function g we evaluate ξ and I and from the eigenstates ψ wecompute PN . Figure 1(a) shows the results for λ = 1for which all the single particle eigenstates are extended.The plot of Fig. 1(a) shows I and PN vs E for systemsizes N = 250 , H for N = 500. We see thatboth participation numbers PN and I drop to zero inin the gaps of the spectrum of H , and increase with thesystem size for energies where eigenstates are present.We observe I > PN in general. Figure 1(b) comparesthe same quantities for λ = 2 . I and ξ against E for N = 250 , ξ ( e )1 = 1 / ln(1 . ≈ . ξ evaluated from the Green func-tion (8) is close to the exact value ξ ( e )1 for energies closeto the eigenenergies of the system, while I is systemat-ically larger than ξ , but is roughly of the same order, and does not scale with the system size N . In the gapsof the exact spectrum, I drops to zero which is expectedsince there are no eigenstates corresponding to these en-ergies and contributions from the eigenstates are negli-gible. However ξ , defined by Eq. (8), gives completelywrong value in the gaps of the single particle spectrumas seen in Fig. 1(b). This is clearly an artefact of the ex-ponential fitting of g that does not decay exponentiallyinside the gaps of the spectrum of H . The behaviour ofthe participation number PN is very similar to I (notshown). For the critical case λ = 2, the behaviour of ξ and I is similar to that of the delocalised λ = 1 case. N10 I q ( N ) f(x)=a x b q=2q=3q=4q=5q=6 FIG. 2. (Colour online) Participation number I q ( N ) (sym-bols) vs system size N = 250 to 3000 for q = 2 , , , , I q ( N ) = aN b (dashed lines) for λ = 1. Thepower-law fit works well also for λ = 2 , . This rough comparison lends support to the validityof I as a substitute for the participation number PN .To strengthen this support we look into the scaling ofthe participation numbers PN q with the power q , whichalso distinguishes extended, localised and (multi)fractalstates: PN q = aN D q ( q − where D q is the fractal di-mension of the state, and D q = 0 corresponds to lo-calised state, D q = 1 corresponds to delocalised states,0 < D q < I q and try the fit I q = aN D q ( q − for all the three regimes: λ = 1 , , .
5. We pick the en-ergy E max corresponding to the maximum of I for thelargest system size considered, N = 3000, since we wantto probe the most delocalised states in an otherwise lo-calised regime (this choice is only relevant for λ = 2 . E max to evaluate D q for smaller system sizes. For every λ we compute I q for q = 2 , , , , N = 250 to 3000. The results are shownin Fig. 2: we see a clear power law scaling of I q with N for every individual value of q . Next we fit these datafor several system sizes to extract D q for the values of λ = 1 , , .
5. Similarly we evaluate the D q from the scal-ing of PN q with system size. The PN q are computedfrom exact diagonalisation of a single particle Hamilto-nian (1). The results are summarised in Fig. 3: both D ( q ) = b ( q ) / ( q - ) GF, = 1.0(delocalised)GF, = 2.0\(critical)GF, = 2.5(localised) ED, = 1.0ED, = 2.0ED, = 2.5
FIG. 3. (Colour online) Fractal dimension D q vs q obtainedfrom Green’s function (GF) and exact diagonalisation (ED)in extended (red/yellow), critical (blue/cyan) and localised(green/magenta) regimes for the single particle case. Thedimension D q is q -independent and equal to zero(one) in theextended(localised) regime and has a non-trivial dependenceon q at the criticality, λ = 2. This a similar behaviour to thedimension D q computed from the PN q . methods agree - D q ≈ D q ≈ λ = 1 . D q ≈ D q ≈ . λ = 2 .
5, and q -dependent D q , D q for the criticalvalue λ = 2 . I can be used as a substi-tute for the localisation length ξ and the participationnumber PN q in the single particle case. We assume thatthis is also the case for two interacting particles and ver-ify this assumption self-consistently. Therefore in whatfollows we will study the behaviour of I and higher mo-ments I q> . IV. TWO INTERACTING PARTICLES:SELF-SIMILARITY OF THE SPECTRUM ANDFRACTALITY OF THE EIGENSTATES
We now turn to the case of two particles with the on-site Hubbard interaction. Earlier work has reportedthe emergence of metallic states in the single-particle in-sulating regime ( λ > N = 250 (up to N = 1000 with sparse diagonalisation) sites and analysisof the spreading of time-evolved wave packets scanned inthe entire range of interactions 0 < u <
12 for severalvalues of potential strength λ ∈ [1 . , , who performed diagonalisa-tion of systems up to N = 10946 sites and confirmed thepresence of delocalised states.We start our analysis with a cross check of Ref. 13 andevaluate the participation number I from ˜ G nm (5) for I =2.5,u=7.9, N=250 E = 0.015 FIG. 4. (Colour online) Participation number I vs energy E for two interacting particles. The peaks signal the emergenceof delocalised states in the otherwise localised spectrum. λ = 2 . u = 7 . E ∈ [ − , β , see Eq. (1). In Fig. 4 we see a minibandsstructure with the few energies where the value of I is relatively large, similarly to the findings of Ref. 13thereby lending further support to the use of I as aprobe of the extent of the eigenstates. We identified twovalues of energy, E ≈ . , E ≈ − . I achievesits local maximum (Fig. 4), suggesting the emergence ofdelocalised states at these energies.To get a better insight into the nature of these emerg-ing states we study the fine structure in the vicinity ofthe I peaks. To extract this fine structure we start witha small system size and identify the peaks of I by dis-cretising the energy range. Next we zoom into the en-ergy range around one of the peaks by using a finer en-ergy discretisation. This procedure is repeated severaltimes for increasing system sizes N . Such analysis offine details of the structure of I is possible thanks tothe usage of the projected Green functions. To be spe-cific for the peak of I at E ≈ .
8, we started with arange or energies [1 .
821 : 1 . N = 250. We observe the emergence of newpeaks which become prominent as the size is increasedto N = 500 and N = 1000 (Fig. 5(a)). Zooming inthe energy range around one peak ( E ∈ [1 . , . N = 3000, Fig. 5(b). Repeating this procedure two moretimes for the peaks marked by the black boxes, we obtainFig. 5(c-d) for N max = 7000. The largest I is observedat E = 1 . N = 7000. Upon every iterationwe observe the emergence of finer structure in I as weare zooming in energy. This strongly suggests the fractalnature of participation number I as a function of en-ergy E and consequently the spectrum of the delocalisedstates at these energies.In the original work, Ref. 13, these states were assumeddelocalised based on the analysis of wave packet spread-ing. Subsequent work in Ref. 14 performed a more de- I N=250N=500N=1000 (a) +1.821 I N=1000N=2000N=3000 (b) +1.821 I N=3000N=4000N=5000 (c) +1.8214 I N=7000N=6000N=5000N=1000 (d)
FIG. 5. (Colour online) Average Green’s function participation number I around energy E ≈ . u = 7 .
9. The energyrange is zoomed in from left to right, with the maximum system size increasing from N = 1000 (left) to N = 7000 (right) andthe resolution in energy reaching ∆ E = 3 ∗ − for the rightmost plot. The errorbars correspond to the disorder average. Thepeaks of I resolve into fine structure with subpeaks upon every iteration of zooming in. I N=1000N=500N=250 (a) I N=3000N=2000N=1000 (b) I N=5000N=4000N=3000N=1000 (c) I N=6000N=5000N=4000 (d)
FIG. 6. (Colour online) Average Green’s function participation number I around energy E ≈ − . u = 7 .
9. The energyrange is zoomed in from left to right, with the maximum system size increasing from N = 1000 (left) to N = 6000 (right) andthe resolution in energy reaching ∆ E = 10 − for the rightmost plot. The errorbars correspond to the disorder average. Thepeaks of I resolve into fine structure with subpeaks upon every iteration of zooming in. tailed analysis and confirmed this conclusion and alsoprovided some indications of fractality of these statesbased on the fitting i) inverse participation ratio in po-sition representation denoted as ξ x , ii) inverse participa-tion ratio in energy representations, ξ E (for details seeRef. 14). To clarify the fractal nature of these stateswe consider the largest I at energy E = E and com-pute I q ( N ) for q = 2 , , , , N at this energy. Assuming the multifractal ansatz forthe participation number I q ( N ) ∼ aN D q ( q − we extractthe fractal dimension D q from numerical values I q ( N ),similarly to how it was done in the single particle case,see Fig. 2. The extracted values of D q are shown as redpoints (circles) in Fig. 7 with the error bars of the fit. Weobserve that D q < q -dependent suggesting that thecorresponding eigenstates at this energy are multifractal.In Ref. 14 a power-law fit of ξ x and ξ E with system size N with N max ≈ a x,E < E = − . , .
817 and interaction u = 7 . E ≈ − . N max = 6000. The results wereaveraged over 10 disorder samples, e.g. values of β (seeEq. (1)). The results are shown in Fig. 6. We observelarger fluctuations in participation number I as com-pared to E ≈ . I on system size N is less promi-nent as compared to the global maximum of I locatedat E ≈ . D q extracted from I q shows an almost flat de-pendence on q (green circles in Fig. 7), suggesting onlyfractal but not multifractal character of the state at thisenergy.The Green function participation number results areindirect, since they do not probe the eigenstates directly.Their advantage is the much lower computational costfor larger system sizes as compared to the exact diag-onalisation. Therefore to check our predictions on thefractality of the eigenstates independently we performedsparse diagonalisation around energies E = 1 . E = − . I for N max = 7000 and N max = 6000 respectively.Among the eigenstates extracted around these two ener-gies, we systematically picked the ones with the largest P N for all system sizes N since we aimed at the mostdelocalised eigenstates embedded into the predominantlylocalised ones. The power law fits of the participationnumber moments, PN q ( N ) ∝ N D q ( q − were calculated.The resulting values of D q are shown in Fig. 7 as blue( E = 1 .
8) and green ( E = − .
8) solid lines with triangu-lar points. The dashed lines with points show D q evalu-ated from the Green function participation numbers, red D ( q ) = b ( q ) / ( q - ) ED, E=1.8ED, E=-2.8 GF, E=1.8GF, E=-2.8
FIG. 7. (Colour online) Fractal dimensions D q (extractedfrom the Green function participation number I q ) and D q (ex-tracted from the participation number P N q ) vs q at energies E ≈ . E ≈ − .
8. For E both methods predict multi-fractality, while for E the projected Green function methodunderestimates the fractality of the eigenstate. for E = 1 . E = − .
8. We see that al-though the values of D q and D q do not always agree per-fectly, nevertheless D q and D q imply at least fractalityof the eigenstates that were previously considered delo-calised. This also provides yet another evidence for thevalidity of I q as a measure of localisation of eigenstates.We elaborate further on the character of these fractalstates appearing around E , . Since the appearance ofthese states relies crucially on the interaction, we expectthem to have a peculiar spatial pattern of the wavefunc-tion amplitudes. Indeed we can construct many approxi-mate localised eigenstates with two particles separated byone or more localisation lengths ξ . Therefore the frac-tal states should have the two particles separated by atmost the single particles localisation length ξ . If we vi-sualise the amplitudes of the two particles eigenfuncion | Ψ( x , x ) | on a square lattice with coordinates x , x ,that correspond to the positions of the two particles, weexpect the fractal states to be localised along the maindiagonal x = x , the fractal structure translating intosome complicated pattern along the main diagonal. Toverify this hypothesis we plotted two exact eigenstateswith the largest P N for N = 3000 in Fig. 8 ( u = 7 . E ≈ . E ≈ − . | Ψ( x , x ) | < − on the plots. These plotsfully confirm our hypothesis outlined above, with mostweight concentrated along the main diagonal, i.e. bothparticles being close to each other. E=1.82138, u=7.9
E=-2.78277, u=7.9
FIG. 8. (Colour online) The amplitudes | Ψ( x , x ) | of eigen-states computed for N = 3000 at u = 7 . . The X -axis and Y -axis denotethe positions of the two particles - x and x - respectively.The larger amplitudes correspond to brighter colour. Valuessmaller than 10 − were discarded. The energies are E ≈ . E = − . V. CONCLUSIONS
To conclude, we have shown that previously discov-ered metallic states of two interacting particles in anAA chain in the insulating single-particle region have afractal structure. Furthermore unlike previous claims wefind that these states are multifractal. This is verified bycomputing participation numbers from projected GF aswell as from exact diagonalisation. An interesting openproblem is the fate of these multifractal states at finitedensity where many-body localisation was reported athalf-filling. As a side effect, we demonstrated that the projectedGreen functions can be used as a first probe to check thenature of eigenstates in an interacting Hamiltonian sys-tem having the advantage that larger system sizes can betargeted as compared to the computationally challengingexact diagonalisation.
ACKNOWLEDGMENTS
This work was supported by the Institute for BasicScience in Korea (IBS-R024-D1). P. W. Anderson, “Absence of diffusion in certain randomlattices,” Phys. Rev. , 1492–1505 (1958). B. Kramer and A. MacKinnon, “Localization: theory andexperiment,” Rep. Prog. Phys. , 1469–1564 (1993). D.M. Basko, I.L. Aleiner, and B.L. Altshuler, “Metalin-sulator transition in a weakly interacting many-electronsystem with localized single-particle states,” Ann. Phys. , 1126 – 1205 (2006). Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, andMaksym Serbyn, “Colloquium: Many-body localization,thermalization, and entanglement,” Rev. Mod. Phys. ,021001 (2019). Fabien Alet and Nicolas Laflorencie, “Many-body localiza-tion: An introduction and selected topics,” Comptes Rend.Phys. , 498 – 525 (2018). Hidetoshi Fukuyama, Robert A. Bari, and Hans C.Fogedby, “Tightly bound electrons in a uniform electricfield,” Phys. Rev. B , 5579–5586 (1973). Serge Aubry and Gilles Andr´e, “Analyticity breaking andanderson localization in incommensurate lattices,” Ann.Israel Phys. Soc , 18 (1980). D. L. Shepelyansky, “Coherent propagation of two inter-acting particles in a random potential,” Phys. Rev. Lett. , 2607–2610 (1994). K.M. Frahm, “Interaction induced delocalization of twoparticles: large system size calculations and dependence oninteraction strength,” Eur. Phys. J. B , 371–378 (1999). D. O. Krimer, R. Khomeriki, and S. Flach, “Two inter-acting particles in a random potential,” JETP Lett. ,406–412 (2011). Diana Thongjaomayum, Alexei Andreanov, Thomas Engl, and Sergej Flach, “Taming two interacting particles withdisorder,” Phys. Rev. B , 224203 (2019). Ramaz Khomeriki, Dmitry O. Krimer, Masudul Haque,and Sergej Flach, “Interaction-induced fractional blochand tunneling oscillations,” Phys. Rev. A , 065601(2010). Sergej Flach, Mikhail Ivanchenko, and Ramaz Khome-riki, “Correlated metallic two-particle bound states inquasiperiodic chains,” EPL , 66002 (2012). Klaus M. Frahm and Dima L. Shepelyansky, “Freed byinteraction kinetic states in the harper model,” Eur. Phys.J. B , 337 (2015). Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, Hen-rik P. L¨uschen, Mark H. Fischer, Ronen Vosk, Ehud Alt-man, Ulrich Schneider, and Immanuel Bloch, “Observa-tion of many-body localization of interacting fermions in aquasirandom optical lattice,” Science , 842–845 (2015). Felix von Oppen, Tilo Wettig, and Jochen M¨uller,“Interaction-induced delocalization of two particles in arandom potential: Scaling properties,” Phys. Rev. Lett. , 491–494 (1996). W. E. Arnoldi, “The principle of minimized iteration in thesolution of the matrix eigenvalue problem,” Quart. Appl.Math. , 17–29 (1951). Klaus M. Frahm, “Eigenfunction structure and scaling oftwo interacting particles in the one-dimensional andersonmodel,” Eur. Phys. J. B , 115 (2016). Shankar Iyer, Vadim Oganesyan, Gil Refael, and David A.Huse, “Many-body localization in a quasiperiodic system,”Phys. Rev. B87