S J Thomson
email: steven.thomson@fu-berlin.de

M Schiró

Local Integrals of Motion in Quasiperiodic Many-Body Localized Systems

Local Integrals of Motion in Quasiperiodic Many-Body

Localized Systems

Introduction

The study of quenched random disorder in non-interacting quantum systems dates back to the seminal work of Anderson [START_REF] Anderson | Absence of diffusion in certain random lattices[END_REF] who first understood how diffusion could eventually break down in the presence of strong enough randomness, leading to localization of the single-particle wavefunctions. It was later realized [START_REF] Abrahams | Scaling theory of localization: Absence of quantum diffusion in two dimensions[END_REF] that this localization transition occurs in three dimensional systems, while lower dimensional non-interacting quantum systems are always localized for any small disorder. The complex interplay of quantum many-body interactions and random disorder was studied thereafter [START_REF] Fleishman | Interactions and the anderson transition[END_REF], but has been the subject of renewed attention over the last decade since the theoretical prediction of localization at finite temperatures in interacting systems [START_REF] Basko | Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states[END_REF], followed by rapid experimental verification [START_REF] Smith | Many-body localization in a quantum simula-tor with programmable random disorder[END_REF][START_REF] Choi | Exploring the many-body localization transition in two dimensions[END_REF][START_REF] Hess | Non-thermalization in trapped atomic ion spin chains[END_REF] and further theoretical progress, summarized in Refs. [START_REF] Altman | Universal dynamics and renormalization in many-bodylocalized systems[END_REF][START_REF] Nandkishore | Many-body localization and thermalization in quantum statistical mechanics[END_REF][START_REF] Alet | Many-body localization: An introduction and selected topics[END_REF][START_REF] Abanin | Colloquium: Many-body localization, thermalization, and entanglement[END_REF]. Now known as many-body localization (MBL), this is primarily a dynamical effect whereby highly excited eigenstates may spontaneously fail to thermalize due to disorder, and is typically understood in terms of the existence of an extensive number of locally conserved quantities known as integrals of motion [START_REF] Serbyn | Local conservation laws and the structure of the many-body localized states[END_REF][START_REF] Huse | Phenomenology of fully many-bodylocalized systems[END_REF][START_REF] Ros | Integrals of motion in the many-body localized phase[END_REF][START_REF] Imbrie | Local integrals of motion in many-body localized systems[END_REF] (LIOMs, or l-bits) which prevent the system from thermalizing, in contrast to ergodic systems which conserve only global quantities such as the total energy or number of particles. The breakdown of MBL is understood in terms of resonances, i.e. spatially separated sites with on-site energies which differ by less than the hopping amplitude, facilitating transport in the system. Proliferation of these resonances can lead to so-called avalanche effects, which can lead to transport throughout the entire system and consequently destroy localization [START_REF] De Roeck | Stability and instability towards delocalization in many-body localization systems[END_REF][START_REF] Thiery | Many-body delocalization as a quantum avalanche[END_REF]. Understanding the role of these resonances and their distribution within a given system is crucial in order to understand the stability of disorder-induced localization.

While many-body localization in systems subject to a random external potential is by now largely understood in terms of LIOMs, in some cases rigorously [START_REF] Imbrie | On many-body localization for quantum spin chains[END_REF], there is no corresponding rigorous understanding of MBL phenomena in systems with quasiperiodic potentials of the sort routinely realised in current experiments, despite the fact that both experimental [START_REF] Schreiber | Observation of many-body localization of interacting fermions in a quasirandom optical lattice[END_REF][START_REF] Bordia | Coupling identical one-dimensional many-body localized systems[END_REF][START_REF] Lüschen | Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems[END_REF][START_REF] Bordia | Probing slow relaxation and many-body localization in two-dimensional quasiperiodic systems[END_REF] and numerical evidence exists for MBL phenomena in such systems [START_REF] Iyer | Many-body localization in a quasiperiodic system[END_REF][START_REF] Lee | Many-body localization in spin chain systems with quasiperiodic fields[END_REF][START_REF] Chandran | Localization and symmetry breaking in the quantum quasiperiodic ising glass[END_REF][START_REF] Nag | Many-body mobility edges in a one-dimensional system of interacting fermions[END_REF][START_REF] Khemani | Two universality classes for the many-body localization transition[END_REF]. Quasiperiodic potentials may be loosely understood as intermediate between periodic and random potentials, and are characterised by a set of discrete but non-periodic Fourier components. First studied contemporaneously with Anderson localization [START_REF] Harper | Single band motion of conduction electrons in a uniform magnetic field[END_REF] and further investigated in the decades since [START_REF] Aubry | Analyticity breaking and anderson localization in incommensurate lattices[END_REF], they exhibit many crucial differences with respect to random disorder as the potential is strongly correlated in real space. This dramatically changes the spatial structure of resonances already in single-particle models [START_REF] Szabó | Non-power-law universality in one-dimensional quasicrystals[END_REF], enhancing the probability of encountering resonances between spatially separated sites, and could have important implications for the stability of many-body localized phases of matter in quasicrystalline materials [START_REF] Žnidarič | Interaction instability of localization in quasiperiodic systems[END_REF]. Furthermore rare-region Griffiths effects [START_REF] Vojta | Quantum griffiths effects and smeared phase transitions in metals: Theory and experiment[END_REF], which are responsible for much of the exotic phenomenology of MBL systems on both sides of the transition and possibly also for the critical behavior [START_REF] Vosk | Theory of the many-body localization transition in one-dimensional systems[END_REF][START_REF] Agarwal | Anomalous diffusion and griffiths effects near the many-body localization transition[END_REF][START_REF] Bar Lev | Absence of diffusion in an interacting system of spinless fermions on a one-dimensional disordered lattice[END_REF][START_REF] Schiró | Toy model for anomalous transport and griffiths effects near the many-body localization transition[END_REF][START_REF] Schulz | Phenomenology of anomalous transport in disordered one-dimensional systems[END_REF], are absent in quasiperiodic systems. This points towards a qualitatively different physics at play in the quasiperiodic MBL problem, as well as towards a different mechanism responsible for its eventual instability towards full delocalisation, an issue which is still far from being fully understood [START_REF] Turkeshi | Destruction of localization by thermal inclusions: Anomalous transport and Griffiths effects in the Anderson and André-Aubry-Harper models[END_REF]. Recently there have been several studies presenting contrasting results for the nature of the MBL transition in quasiperiodic systems, with the level statistics obtained from exact diagonalization [START_REF] Khemani | Two universality classes for the many-body localization transition[END_REF] giving a rather different result for the transition point and critical exponent than renormalization group calculations [START_REF] Zhang | Universal properties of many-body localization transitions in quasiperiodic systems[END_REF] and studies based around the construction of LIOMs using time-averaged local observables [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF], as first proposed in Ref. [START_REF] Chandran | Constructing local integrals of motion in the many-body localized phase[END_REF].

In this work we explicitly construct the LIOMs for a model of interacting one dimensional fermions in a quasiperiodic potential and use them to gain insight on the nature of the MBL phase and its possible transition towards delocalization. To this extent we develop a new computational implementation of the established flow equation method, leveraging the high degree of parallelization possible with modern computer hardware and tensor algebra operations. We show that this tensor flow equation leads to a significant improvement upon previous implementations [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] in both accuracy, thanks to the inclusion of additional running couplings in the flow parametrization, and transparency. Using this new method we compute explicitly the LIOM interactions and real-space support, showing how they contain clear fingerprints of the quasiperiodic structure of the potential in the form of local dips or peaks at distances given by the associated asymptotic numbers. Furthermore we show that the richer structure of the tensor flow equation and of its LIOM description allows us to build up a diagnostic for the MBL phase transition, so far out of reach for previous implementations of the method, that we estimate here for the quasiperiodic problem using finite-size scaling. We end by discussing our results in the contexts of other recent related works [START_REF] Zhang | Universal properties of many-body localization transitions in quasiperiodic systems[END_REF][START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF][START_REF] Aramthottil | Finite-size scaling analysis of many-body localization transition in quasi-periodic spin chains[END_REF].

We wish to emphasise here at the outset that by comparison with previous works making use of truncated flow equation methods [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Kelly | Exploring many-body localization in quantum systems coupled to an environment via wegner-wilson flows[END_REF], this new work offers the following major advantages:

• Inclusion of high-order off-diagonal terms in the ansatz for the running Hamiltonian, which leads to a significant improvement in the accuracy and stability of the method.

• Ability to compute the many-body contributions to the LIOM operator expansion, which in turn allows us to construct a diagnostic for the MBL phase transition.

• Formulation in terms of tensors and tensor contractions allows us to make use of highly efficient numerical linear algebra routines in order to implement the technique in an entirely new way that is not only faster to compute, but is more transparent and offers a clear path towards systematic extension of the method to higher orders and other systems of interest.

Model

We start from a one-dimensional fermionic Hamiltonian of the form:

H = i h i : n i : +J i (: c † i c i+1 : +H.c.) + ∆ 0 i : n i n i+1 : (1)
where the on-site terms h i are drawn from a quasiperiodic potential given by h i = W cos(2πi/φ+ θ) where φ is some irrational number, θ is a random phase and W is the amplitude of the quasiperiodic potential (i.e. plays the role of 'disorder strength'), with nearest-neighbour hopping J and interactions ∆ 0 . For convenience, we shall refer to the choice of phase θ as a 'disorder realization' and to averages over different values of θ as 'disorder averaging', though it should be noted that quasiperiodic potentials are not true disorder as they are entirely deterministic and there is no randomness. The notation : ... : signifies vacuum normal-ordering, which is necessary to enforce a consistent ordering of operators -see Appendix A for details.

In the absence of interactions, i.e. ∆ 0 = 0, Eq. 1 reduces to the well known Aubry-André (or Aubry-André-Harper) model [START_REF] Harper | Single band motion of conduction electrons in a uniform magnetic field[END_REF][START_REF] Aubry | Analyticity breaking and anderson localization in incommensurate lattices[END_REF][START_REF] Grempel | Localization in an incommensurate potential: An exactly solvable model[END_REF], which despite being a single-particle system exhibits a great deal of non-trivial physics. This model undergoes a localization transition at a quasiperiodic potential strength of W c /J = 2, which can be shown analytically through the self-duality of the Hamiltonian. This is unlike non-interacting one-dimensional systems subject to a random external potential, which are always Anderson localized and do not exhibit a delocalization transition. Above this critical value of W c /J, the single-particle wavefunctions decay exponentially in space, but with quasiperiodic revivals due to the structure of the underlying potential. Precisely at the critical point, the single particle spectrum is fractal in nature, and takes the form of a Cantor set [START_REF] Yao | Critical Behavior and Fractality in Shallow One-Dimensional Quasiperiodic Potentials[END_REF]. The Aubry-André model does not have a mobility edge, and although variations exist which do exhibit single-particle mobility edges [START_REF] Ganeshan | Nearest neighbor tight binding models with an exact mobility edge in one dimension[END_REF][START_REF] Duthie | Localization in quasiperiodic chains: A theory based on convergence of local propagators[END_REF], we will not consider those modifications in the present work.

The interacting Aubry-André model is comparatively less understood, although it has received some attention in recent years motivated by the MBL problem [23-27, 31, 39]. In the limit of weak interactions, the bosonic Aubry-André model in one dimension has been shown to exhibit an inverted mobility edge, with the high-energy excitations appearing exponentially localized while the low-energy excitations remain extended [START_REF] Lellouch | Localization transition in weakly interacting bose superfluids in one-dimensional quasiperdiodic lattices[END_REF]. The low-energy properties of bosonic versions of the interacting Aubry-André model have been studied recently in both one [START_REF] Yao | Lieb-liniger bosons in a shallow quasiperiodic potential: Bose glass phase and fractal mott lobes[END_REF] and two [START_REF] Sbroscia | Observing Localization in a 2D Quasicrystalline Optical Lattice[END_REF][START_REF] Gautier | Strongly Interacting Bosons in a Two-Dimensional Quasicrystal Lattice[END_REF] spatial dimensions and shown to exhibit unusual fractal properties. Interacting quasiperiodic spin chains in one dimension have been studied using tensor network techniques [START_REF] Doggen | Many-body delocalization dynamics in long aubry-andré quasiperiodic chains[END_REF] to simulate their non-equilibrium dynamics, and recently with exact diagonalization [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF] in order to construct time-averaged local integrals of motion. Fermionic Aubry-André models have also been studied in the context of many-body localization using quantum quench spectroscopy to identify mobility edges [START_REF] Naldesi | Detecting a many-body mobility edge with quantum quenches[END_REF]. Here we will focus on the one-dimensional fermionic model (which maps onto the spin model via a Jordan-Wigner transform), however the method we outline in the following section is not unique to fermions and may also be applied to interacting bosons with only minor modifications.

Method

Here we will first review the use of continuous unitary transforms, also known as flow equation methods, before detailing the new implementation we use in the present manuscript which we find offers many practical advantages over previous implementations, while still sharing the same underlying philosophy. The new implementation makes use of the Python package PyFlow [START_REF] Thomson | PyFlow: A python package for diagonalising quantum systems using continuous unitary transforms[END_REF], developed by one of the present authors, and running on graphics processing units (GPUs).

Brief Review of Continuous Unitary Transforms

The use of continuous unitary transforms to diagonalize Hamiltonians has a long and varied history. The technique was originally proposed in condensed matter physics by Wegner [START_REF] Wegner | Flow-equations for hamiltonians[END_REF] under the name 'flow equations' (later popularized by Kehrein [START_REF] Kehrein | The flow equation approach to many-particle systems[END_REF] and coworkers [START_REF] Kehrein | Diagonalization of system plus environment hamiltonians[END_REF][START_REF] Kehrein | Scaling and decoherence in the nonequilibrium kondo model[END_REF][START_REF] Moeckel | Interaction quench in the hubbard model[END_REF][START_REF] Hackl | Real time evolution in quantum many-body systems with unitary perturbation theory[END_REF][START_REF] Hackl | A unitary perturbation theory approach to real-time evolution problems[END_REF]), independently in high-energy physics by Glazek and Wilson under the name 'similarity transform' [START_REF] Wilson | Renormalization of hamiltonians[END_REF][START_REF] Glazek | Perturbative renormalization group for hamiltonians[END_REF], and also in mathematics under the names 'isospectral flow' and 'double bracket flow' [START_REF] Brockett | Dynamical systems that sort lists, diagonalize matrices, and solve linear programming problems[END_REF][START_REF] Chu | The projected gradient method for least squares matrix approximations with spectral constraints[END_REF][START_REF] Chu | A list of matrix flows with applications[END_REF]. Since then the method has also been generalized to time-dependent systems in a variety of forms [START_REF] Kehrein | Diagonalization of system plus environment hamiltonians[END_REF][START_REF] Tomaras | Scaling approach for the time-dependent kondo model[END_REF][START_REF] Verdeny | Accurate effective hamiltonians via unitary flow in floquet space[END_REF][START_REF] Rosso | Dissipative flow equations[END_REF][START_REF] Thomson | Flow Equations for Disordered Floquet Systems[END_REF] including driven and dissipative scenarios, however here we will focus on Hamiltonians with no explicit time dependence.

The core of the method is the idea of using a series of infinitesimal unitary transformations to diagonalise a particular system of interest. This is spiritually similar to the well-known Schrieffer-Wolff [START_REF] Schrieffer | Relation between the anderson and kondo hamiltonians[END_REF][START_REF] Bravyi | Loss, Schrieffer-wolff transformation for quantum many-body systems[END_REF] transform, where a Hamiltonian may be diagonalized to leading order by careful choice of a suitable unitary transformation:

H = e S He -S = H -[H, S] + ..., (2)
and choosing the generator S such that [H, S] is equal to the off-diagonal terms of the Hamiltonian which we wish to remove. In general, this procedure will also generate new higher-order terms which must then be removed by further transforms, or treated in some other perturbative manner. Here, rather than making a single 'large' unitary transform, we instead make a series of infinitesimal transforms, each of which can be made arbitrarily accurate:

H(l + dl) = e η(l)dl H(l)e -η(l)dl = H(l) + dl [η(l), H(l)]. ( 3
)
where l is a fictitious 'flow time' which parameterises the transform, and η(l) is some scaledependent anti-Hermitian generator chosen to diagonalize the specific problem of interest in the l → ∞ limit. In the spirit of renormalization group techniques, the transformation of the Hamiltonian can be written in terms of a single so-called flow equation:

dH(l) dl = [η(l), H(l)] (4)
such that the eventual 'fixed point' of the flow is a diagonal Hamiltonian, with [η(l → ∞), H(l → ∞)] = 0. The properties of the final fixed-point Hamiltonian are controlled by the generator η(l), which we can freely choose to be any anti-Hermitian operator such that the overall transform is unitary. The choice of generator is far from unique and various options exist in the literature which result in a diagonal fixed-point Hamiltonian [START_REF] Monthus | Flow towards diagonalization for many-body-localization models: adaptation of the toda matrix differential flow to random quantum spin chains[END_REF][START_REF] Savitz | Stable unitary integrators for the numerical implementation of continuous unitary transformations[END_REF], however a common choice is the so-called 'Wegner generator':

η(l) = [H 0 (l), V (l)] (5)
where H 0 represents the diagonal terms in the Hamiltonian and V = H -H 0 contains the off-diagonal terms which we want to vanish in the l → ∞ limit. The Wegner generator is suitable for most problems as it a robust choice that can be stably numerically integrated, although it has the significant disadvantage for sparse models (i.e. those with nearest-neighbour couplings only) that it does not preserve the sparsity of the initial Hamiltonian, and in the early stages of the flow it typically generates long-range off-diagonal couplings which must be kept track of until they later decay to zero. Other choices of generator exist which do preserve the sparsity of nearest-neighbour models [START_REF] Monthus | Flow towards diagonalization for many-body-localization models: adaptation of the toda matrix differential flow to random quantum spin chains[END_REF], however they are typically much less numerically stable [START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF], and so we do not consider them here. One interesting frontier is the development of adaptive generators, following an intriguing proposal for adaptive Schrieffer-Wolff transforms [START_REF] Wurtz | Variational Schrieffer-Wolff transformations for quantum many-body dynamics[END_REF], however this is still an active area of development and in the present work we focus on the established Wegner generator. Flow equations in disordered systems have been the subject of increasing attention over the last few years with a variety of implementations including non-interacting systems with long-range couplings [START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF][START_REF] Quito | Localization transition in one dimension using wegner flow equations[END_REF], numerically exact studies of interacting systems [START_REF] Pekker | Fixed points of wegnerwilson flows and many-body localization[END_REF][START_REF] You | Bulk geometry of the many body localized phase from wilson-wegner flow[END_REF], analytical treatments proposing efficient generators [START_REF] Monthus | Flow towards diagonalization for many-body-localization models: adaptation of the toda matrix differential flow to random quantum spin chains[END_REF] and approximate methods for interacting systems [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Kelly | Exploring many-body localization in quantum systems coupled to an environment via wegner-wilson flows[END_REF] which made use of a truncation in operator space to limit the number of running couplings to a computationally tractable number even for large systems, enabling studies of MBL phenomenology not only in one-dimension, but also in coupled chains and even in two spatial dimensions.

Here we will build on the latter approach, employing a truncation in operator space together with an entirely new implementation more suited to the efficient linear algebra capabilities of modern computer hardware and extending the truncation to include new higher-order terms which, as we shall show, dramatically increases the accuracy of the method and enables us to do away with some of the complexities required in earlier works. Moreover, this new implementation allows us to compute new quantities which were out of the reach of previous implementations, including the identification of a delocalization transition in a quasiperiodic many-body system.

Tensor Flow Equations

To illustrate the new implementation, which we call the Tensor Flow Equation (TFE) technique, here we start from a very general normal-ordered (fermionic) Hamiltonian with quadratic and quartic terms:

H = ij H (2) ij + kqlm H (4) kqlm (6) = ij H (2) ij : c † i c j : + kqlm H (4) kqlm : c † k c q c † l c m : (7)
For simplicity we will assume this Hamiltonian to describe a one-dimensional chain of length L, however this procedure can easily be generalized to two-dimensional models, as with previous implementations of flow equations [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF]. In all of the following, the bracketed superscripts indicate how many fermionic operators are associated with a term. The couplings in this Hamiltonian may be short-or long-ranged, and may be homogeneous or disordered, therefore this form incorporates a large range of systems of interest1 . This Hamiltonian can be represented schematically with the following diagram:

H = i j H (2) + k q l m H (4)
where H (2) is a rank-2 tensor representing the quadratic part of the Hamiltonian, H (4) is a rank-4 tensor representing the quartic part of the Hamiltonian, the 'out' arrows signify fermionic creation operators (c † i ), and the 'in' arrows signify fermionic annihilation operators (c j ). By construction, the end result of our procedure will be a Hamiltonian which is diagonal in terms of the fermionic number operators and will take the form:

H = i H(2) ii ñi + ij H(4) iijj ñi ñj + ... (8)
H = i i H(2) + i i j j H (4)
where the tilde notation indicates that all quantities are given in the diagonal basis and the ellipsis indicates possible higher-order terms which we shall discuss later. For flow equations, the key ingredient is the calculation of commutation relations between different parts of the Hamiltonian. For example, the canonical Wegner generator is given by: η = [H 0 , V ] = H (2) , V (2) + H (4) , V (2) + H (2) , V (4) + ...

where V represents the off-diagonal part of the Hamiltonian, which we can choose freely. Let's take these three commutators individually:

i) H (2) , V (2) = ijk H (2) ik V (2) kj -V (2)
ik H

(2) kj [START_REF] Alet | Many-body localization: An introduction and selected topics[END_REF] which is just a series of matrix contractions which can be represented schematically by i k

H (2) k j V (2) - i k V (2)
k j H (2) where the line joining the two legs marked k represents a sum over this index, i.e. a matrix/tensor contraction. These tensor contractions can be efficiently computed numerically by most standard linear algebra packages, and can be effectively parallelised over multiple cores -if the arrays fit in memory, then graphical processing units (GPUs) are ideal for this, as modern GPUs typically have 10 3 -10 4 cores, compared to even a cluster computer which may have only ∼ 10 2 cores per node, and this can result in a significant speed increase -see Appendix C for benchmarks demonstrating this. We can follow the same procedure for the higher order terms, where the commutator (see Appendix A) can be written as:

ii) H (4) , V (2)

= ijkql H (4) ljkq V (2) il + H (4) ijlq V (2) kl + H (4) ilkq V (2) lj + H (4) ijkl V
(2) lq [START_REF] Abanin | Colloquium: Many-body localization, thermalization, and entanglement[END_REF] This can be represented graphically as the sum of all possible two-point tensor contractions where each index is from a different tensor:

l j k q H (4) i l V (2) + i j l q H (4) k l V (2) + i l k q H (4) l j V (2) + i j k l H (4) l q V (2)
where the contracted indices are summed over. This reproduces the same results as obtained in Appendix A using a lengthier algebraic procedure. After evaluating the contractions, this can be written graphically as:

j k q i H (4) ljkq V (2) il + i j q k H (4) ijlq V (2) kl + i k q j H (4) ilkq V (2) lj + i j k q H (4) iljk V (2) lq
where here we use the Einstein summation convention such that any repeated indices are summed over. By swapping the order of indices and keeping in mind fermionic anticommutation relations for normal-ordered operators, we can collect together the terms with the correct ± prefactors.

i j k q -H (4) ljkq V (2)
il -H

ijlq V

(2)

kl + H (4) ilkq V (2) lj + H (4) iljk V (2) lq
Here we can start to see why the tensor notation is useful -the contractions and the signs can be read off from the graphical notation. This also enables us to generalise the above method easily to bosons, without having to work out Wick's Theorem explicitly in full to compute the products of normal-ordered operators (see Appendix A), as the commutation relations are automatically taken into account when re-ordering the indices. By the same token, this graphical notation allows straightforward systematic extension to even higher order terms, as commutators involving six or more fermionic operators can be broken down into sums of elementary two-point contractions and dealt with as above.

Similarly, the third part of Eq. 3 can be written as:

iii) H (2) , V (4)

= i j k q V (4) ljkq H (2) il + V (4) ijlq H (2) kl -V (4) ilkq H (2) lj -V (4)
iljk H

(2) lq

In the end this leads to an anti-Hermitian generator of the form:

η = i j η (2) + k q l m η (4)
We can then compute the flow of the Hamiltonian order-by-order using the standard flow equation for H: 2) , H (2) + η (4) , H (2) + η (2) , H (4) + ...

dH dl = [η, H] = η (
which we can evaluate in exactly the same way as we computed η above, using tensor contractions. The ellipsis (...) represent terms containing more than four fermionic operators which will be generated by the flow. The lowest order neglected term will be obtained by contracting H (4) with η (4) to generate a new term H (6) , which should then be added back into the ansatz for the running Hamiltonian, which will then in turn lead to another new term H (8) , and so on. The source term for the lowest order neglected term, H (6) comes from [η (4) , H (4) ] = [[H (4) , V (2) ] + [H (2) , V (4) ], H (4) ], therefore it is at most quadratic in the microscopic interaction strength. For weak interactions, each of the newly generated terms will be smaller than the previous order. In the following, we shall assume that we work with sufficiently weak interactions that all newly-generated terms above fourth order are negligible; we will later on examine the validity of this assumption. At this point, we proceed to numerically solve Eq. 12 directly. In contrast to previous work using flow equations to study disordered many-body systems [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF], here we do not analytically write down explicit equations of motion for the coupling constants, but instead generate the flow of the couplings numerically using the above tensor contraction procedure. This has the significant advantage that it avoids a great deal of algebraic complexity and is far more systematic and transparent, allowing us to more easily incorporate higher-order terms and even extend this formalism further than in the present work, something which was extremely challenging in our previous formulation of a flow equation technique [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF].

Numerical Procedure

Here we outline the specific steps required to diagonalize a Hamiltonian using the above method. In this work we will always work in real space as we are primarily interested in (quasi)-disordered systems, but it should be noted that this procedure may be used for momentum-space Hamiltonians as well, which may provide advantages for systems where the momentum is a well-defined quantum number.

1. Construct the Hamiltonian H order-by-order as an L × L matrix plus an L × L × L × L tensor, i.e. storing L 2 + L 4 floating point numbers in memory. This memory storage could be lowered significantly by using symmetries, e.g. J ij = J ji but most linear algebra packages require the full tensor to be constructed before they can be efficiently multiplied, which requires a lot of CPU time to rebuild at every flow time step. If sufficient memory is available, it is faster to store the full matrices/tensors rather than a reduced representation of them. This is an area where there is scope for technical improvement which may facilitate the study of larger system sizes.

2. Compute the generator η = [H 0 , V ] order-by-order using efficient linear algebra packages to do the tensor contractions specified above, also storing η as L 2 + L 4 floating point numbers in memory.

3. Compute the right hand side (RHS) of the flow equation [η, H] by performing the tensor contractions specified above.

4. Use a suitable numerical integration routine to integrate from H(l) → H(l +dl) with the RHS as constructed above. Here we make use of the JAX package [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF], employing the ode routine with a Runge-Kutta 4th order algorithm and an adaptive step size. When storing the transform, we use a logarithmic grid in flow time, since most of the fastest changes occur at the very start of the flow and the flow-time dynamics slow significantly as l becomes large. This logarithmic grid allows us to more efficiently store the full transform in memory, if necessary, as we focus on storing the greatest number of points in the region where the couplings vary quickly.

5.

Repeat until all off-diagonal elements decay (l → ∞ limit). In practice, we numerically integrate to some large finite value of l. At each step we consider any couplings smaller than some threshold ε to be zero, and typically we use ε = 10 -3 . When all couplings -or at least the majority of them -decay below this threshold, we stop the transform at some finite (usually large) value of l. We have checked and using a smaller cutoff does not lead to noticeably different results, but do however come at greatly increased computational cost -see Appendix B. The use of very small thresholds is only advised if extremely high precision is required, for example computing level statistics where the spacings become exponentially small in the system size.

6. Note: if we want to store the full unitary transform (e.g. in case we want to reverse the transform later, which is required for computing local integrals of motion and nonequilibrium dynamics), we have to store q × (L 2 + L 4 ) floating point numbers, where q is the total number of flow time steps -for large system sizes, this leads to very large memory requirements, and becomes the primary computational limiting factor of this method. This is not a factor if we only wish to diagonalize the Hamiltonian and do not need to store the reverse transform, as in this case we need only store the Hamiltonian itself at each stage (L 2 + L 4 floating point numbers), which is typically much more manageable.

Error Estimation: Invariants of the Flow

As the flow equation method involves a truncation in operator space to a polynomialrather than exponential -number of Hamiltonian coefficients, it is necessarily an approximate method. For small systems, we can compare the results with numerically exact methods (as shown in Appendix B), however the value of this technique lies in its ability to reach larger system sizes than exact methods can access. It is therefore necessary that we develop some form of self-consistent error estimate which can give us confidence that the results can be relied upon even when comparison with exact numerics is not possible. We can achieve this goal using so-called invariants of the flow. These invariants are quantities which are unchanged under the action of a unitary transform, and therefore by comparing them at the start and end of the flow (i.e. computing the invariants of the l = 0 initial Hamiltonian and the l → ∞ final Hamiltonian) we can get a measure of whether the transform has been applied accurately or whether the neglected terms lead to a violation of the unitarity and hence the introduction of an error.

Here we follow the procedure laid out in previous works [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Monthus | Flow towards diagonalization for many-body-localization models: adaptation of the toda matrix differential flow to random quantum spin chains[END_REF] and make use of the second invariant of the flow, defined as:

I 2 (l) = 1 2 L Tr H(l) 2 (13)
for a system of size L. We are interested in whether this quantity changes throughout the flow, and so we define the following quantity of interest:

δI 2 = I 2 (l = 0) -I 2 (l → ∞) I 2 (l = 0) (14)
If δI 2 = 0, the the transform is exact and perfectly unitary. If, on the other hand, we find that δI 2 > 0, this means that terms not included in our ansatz for the running Hamiltonian have at some point during the flow become significant, and we have 'lost' some information by not including them, leading to a deviation from unitarity and thus an error in the final result. In practice, for any finite value of the microscopic interaction strength ∆ 0 > 0, we will find a non-zero value of δI 2 , and as such the reliability of our results rests upon this value being sufficiently small. We will show results for this in the next section.

Numerical Results

For all of the following, we consider the Hamiltonian given in Eq. 1 with hopping amplitude J = 1, nearest-neighbour interactions ∆ 0 = 0.1 and simulate systems of size L = 36, averaged over N s = 128 disorder realizations. Throughout this work we will use open boundary conditions. We consider three different quasiperiodic potentials h i = W cos(2πi/φ + θ) each generated by a different choice of incommensurate ratio φ corresponding to the golden

(φ = (1 + √ 5)/2), silver (φ = 1 + √ 2)
and bronze ratios (φ = (3 + √ 13)/2). These metallic means can all be generated from the recurrence relation a n = ma n-1 + a n-2 , where m is an integer, by computing the ratio a n /a n-1 in the limit of n → ∞. For m = 1 this gives the golden ratio, and the series a n is given by the Fibonacci numbers F n , while for m = 2 one obtains the silver ratio and for m = 3 the bronze ratio. This sequence of metallic means can be continued indefinitely to larger integer values of m.

Under the action of the Wegner generator, the initially short-ranged Hamiltonian will acquire long-range hopping and interaction terms over the course of the flow (the former of which will decay to zero in the l → ∞ limit), and as such we make the following ansatz for the scale-dependent running Hamiltonian:

H(l) = i h i (l) : n i : + ij J ij (l) : c † i c j : + ij ∆ ij (l) : n i n j : + ijkq Γ ijkq (l) : c † i c j c † k c q : (15)
where, as before, we assume that for sufficiently weak interactions any newly generated higherorder terms are negligible, and in the final term we have at least i = j or k = q such that this term is not diagonal in terms of fermionic number operators. We split this running Hamiltonian into diagonal and off-diagonal components, given by:

H(l) = H 0 (l) + V (l) (16)
H 0 (l) = i h i (l) : n i : + ij ∆ ij (l) : n i n j : (17)
V (l) = ij J ij (l) : c † i c j : + ijkq Γ ijkq (l) : c † i c j c † k c q : (18)
We emphasize that as compared to previous implementations of the flow-equation method [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Thomson | Flow Equations for Disordered Floquet Systems[END_REF][START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF] we have now introduced an off-diagonal interaction term, encoded in the tensor Γ ijkq (l). While this will eventually decay to zero at long flow times its presence throughout the flow will play an important role, as we are going to discuss in the following.

The end result of our diagonalization procedure will be a Hamiltonian of the form:

H = i hi ñi + i =j ∆ij ñi ñj + ... (19)
where the ... refers to neglected higher order terms, and the tilde means all quantities are given in the l → ∞ basis. Sample flows of the coefficients h i , J ij , ∆ ij and Γ ijkq are shown in Fig. 1 for a small system to illustrate how these quantities vary throughout the flow. In particular, the variation of all quantities with flow time l is smooth and well-controlled, with the dynamics occurring in two stages: the first stage is at the early flow times where the coefficients all change rapidly, before settling in to a second stage characterized by a slow asymptotic approach to the diagonal Hamiltonian, which by analogy with renormalization group we also refer to as the 'fixed-point Hamiltonian'. In particular, note that during the first stage of the flow, off-diagonal coefficients which are initially zero (e.g. long-range hopping terms) generically become non-zero, but decay back to zero during the course of the flow. The reason for these two distinct stages of the flow is that the initial Hamiltonian contains only nearest-neighbour couplings, and this represents something of an 'unstable fixed point' for the flow equation procedure. The first stage of the process is a rapid flow away from this unstable fixed point, which involves the generation of new couplings, before the second stage of the flow takes over to drive the system to the final diagonal fixed point.

One particularly important point to note before moving on is that several of the interaction terms ∆ ij in Fig. 1b) increase dramatically in the early stages of the flow, before eventually decaying back to much smaller values as the flow time l increases. This rapid increase for small values of l is what led to problems with divergent terms encountered in previous works [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] when attempting to access delocalized phases at small disorder strengths and/or large interaction strengths. In the present formalism, the inclusion of the off-diagonal quartic terms (denoted Γ ijkq in Eq. 15) facilitates this later decay of the density-density interaction coefficients, contrary to what was possible with previous implementations of this technique. In effect, retaining these terms allows them to act like a 'reservoir' for some of the information contained in the initial Hamiltonian, allowing them to become non-zero during the early stages of the flow to 'store' some of this information, before eventually decaying at later stages due

ij (l), with the on-site terms h i (l) indicated by the solid lines with circular markers and the off-diagonal terms J ij (l) indicated by dashed lines with square markers. The black dashed lines indicate the single-particle eigenvalues of the system, to which the H [START_REF] Abrahams | Scaling theory of localization: Absence of quantum diffusion in two dimensions[END_REF] ii terms converge in the l → ∞ limit. Panel (b) shows the flow of the quartic terms H (4) ijkq , again with diagonal terms (i = j, k = q) shown as solid lines and off-diagonal terms shown as dashed lines. For reference, by comparison with exact diagonalization the mean relative error in the many-body eigenvalues of this final Hamiltonian is 8.3 × 10 -3 . to the structure of the generator which guarantees that off-diagonal terms must vanish in the l → ∞ limit. Now that we have seen an example of the flow equation method in action for a small system where we can visualize the behaviour of the couplings, we will move to much larger systems and will investigate two main quantities of interest:

1. The form of the interactions ∆ in the fixed-point Hamiltonian; in random systems, these decay exponentially with distance, however their real-space form in quasiperiodic systems has not been studied.

2. The real-space support of the LIOM density operators ñi ; in the diagonal basis, these operators act on a single lattice site, however we can also express them in the initial microscopic basis to extract their real-space support and see how 'local' the LIOMs really are.

Before presenting these results, however, we will briefly introduce our error estimation measures in order to quantify the accuracy of our methods.

Flow Invariant

First, we will briefly confirm the accuracy of the method by checking the conservation of the flow invariant for the system sizes we will consider in this section. In previous works [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Monthus | Flow towards diagonalization for many-body-localization models: adaptation of the toda matrix differential flow to random quantum spin chains[END_REF], the flow invariant has been computed from an analytic formula which assumes full diagonalisation is achieved. Here we instead compute the flow invariant numerically. We do so by building Here we take the median value of δI 2 , over N s ∈ [START_REF] Wilson | Renormalization of hamiltonians[END_REF]2048] samples, and the error bars show the median absolute deviation. The deviation from unitarity is extremely small in the localized phase where δI 2 appears independent of system size. Deep in the delocalized phase (W/J 2), there is a slow upwards drift of δI 2 with increasing system size. the full Hamiltonian in the basis of all two-particle states (using QuSpin), and compute Eq. 14 directly. This works because the trace can be computed for any orthonormal basis, and the basis of all two-particle states is the smallest which takes into account inter-particle interactions, allowing us to compute this measure numerically exactly even for moderate-to-large system sizes.

The results are shown in Fig. 2, averaged over disorder realizations. Here we use the median rather than the mean, as the median is less sensitive to rare individual disorder realizations which require an abnormally long flow time to converge.

We find that the flow invariant is well conserved for all system sizes at large quasi-disorder strengths, with the largest deviations occurring when W/J is small. In this limit, very large flow times are required to fully diagonalize the problem; loosely speaking, the longer the flow runs for, the larger the fourth-order terms in the Hamiltonian become, and correspondingly the larger the neglected terms of sixth-order and higher will also be. In other words, deep in the delocalized phase where W/J ≤ 2 we may expect the approximate form of the running Hamiltonian to no longer be sufficient, and we see the consequences of this manifest as a deviation from unitary, i.e. a larger value of δI 2 . At intermediate and large values of W/J, however, we see that the flow invariant for all system sizes is conserved to a high degree of accuracy, with deviations on the order of less than 1% for W/J ≥ 2. Furthermore, the flow invariant in this regime does not display significant variation with system size, again providing a strong indication that our method does not lose accuracy as L increases.

We provide further benchmarks and comparisons with exact diagonalization in Appendix B. Having established the accuracy of the numerical technique, we will now turn to the physics of many-body localization in quasiperiodic systems.

LIOM Interactions

We first investigate the form of the LIOM interactions which appear in the fixed-point Hamiltonian, ∆ij . In particular we will examine how these interactions retain a fingerprint of the underlying quasiperiodic potential, contrast these results with the LIOM interactions obtained from a random potential, and show how the deterministic distribution of resonances in the quasiperiodic potentials lead to dramatically different results from the case of purely random disorder.

In Fig. 3 we plot the real-space decay of the LIOM interactions, where the [...] notation signifies the median (typical) magnitude across N s = 128 disorder realizations, and we averaged over distances r = |i -j|. Panels (a), (b) and (c) refer respectively to the golden (φ = (1 + √ 5)/2), silver (φ = 1 + √ 2) and bronze ratio (φ = (3 + √ 13)/2), while for comparison we plot in panel (d) the results for the random disorder case. The overall form of the fixed-point interactions is consistent with the exponential decay seen in the random case [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Rademaker | Explicit local integrals of motion for the many-body localized state[END_REF][START_REF] Rademaker | Many-body localization from the perspective of integrals of motion[END_REF], however there is some additional structure on top of this exponential decay which does not vanish when the number of samples N s is increased; moreover, the same features are present for all disorder strengths and for all system sizes we have tested (L ∈ [START_REF] Serbyn | Local conservation laws and the structure of the many-body localized states[END_REF][START_REF] Wilson | Renormalization of hamiltonians[END_REF]).

In particular we see that the LIOM interactions show dips at specific values of the distance r. For the golden ratio case in panel (a) we have checked that these dips are located at distances given by the Fibonacci numbers, r = F n , as shown by the dashed lines. It is known from earlier work [START_REF] Žnidarič | Interaction instability of localization in quasiperiodic systems[END_REF] that for quasiperiodic potentials generated using the golden ratio, single-particle resonances are enhanced at distances r = F n where F n is a Fibonacci number. Curiously, the role of these resonances in the many-body problem appears to be the suppression of the LIOM interactions at these particular distances. This observation seems to be in contrast with the suggestion of Ref. [START_REF] Žnidarič | Interaction instability of localization in quasiperiodic systems[END_REF] where single-particle resonances at distances r = F n are shown to be enhanced by the quasiperiodic potential, leading to delocalization - (d) Localization length computed from each of the different potentials: all follow a qualitatively similar pattern, with the results from the quasiperiodic potential obtained using the bronze ratio displaying a marginally more localized behaviour than the other two quasiperiodic potentials, and the random box system being more localized than all quasiperiodic potentials at small W/J but less localized at large W/J.

naively, one would expect this to manifest as an increase of the LIOM interactions at distances r = F n , rather than the decrease that we observe here. We have checked carefully (see Appendix B) and verified that this effect is not due to convergence errors, and we shall show in Sec. 4.4 that the single-particle component of the LIOMs indeed behaves as indicated by Ref. [START_REF] Žnidarič | Interaction instability of localization in quasiperiodic systems[END_REF], but that this suppression at r = F n again appears (albeit weakly) in the interaction terms of the LIOM operator expansion. Alternatively, it may be that at distances r = F n ± 1, the LIOM interactions are enhanced by these resonances, leading to the appearance of a dip at r = F n .

The structure of local dips appears also for the other metallic means we considered as shown in panels (b) and (c). Analogously to the Fibonacci series, the silver ratio can be approximated by successive ratios of the Pell numbers P n = 1, 2, 5, 12, 29..., and similarly we see that in Fig. 3b) the fixed-point interactions are suppressed at distances r = P n (gray lines). Curiously, we also see additional dips at distances equal to half the Pell-Lucas numbers, given by 2, 6, 14, 34... and so on. This arises from the fact that the Pell numbers and Pell-Lucas numbers can be combined to give the closest rational approximations to √ 2: the numerators of the approximations are given by the Pell-Lucas numbers, and the denominators given by the Pell numbers, such that the sequence begins 1/1, 3/2, 7/5, 17/2... with each successive fraction giving a closer approximation to √ 2. Due to the close connection between the silver ratio and √ 2, sites with distances r corresponding to these numbers also exhibit a high degree of spatial correlation, resulting in the additional dips seen in Fig. 3b) (blue lines). As with previous metallic means, also the bronze ratio can be approximated by a series of numbers with the recurrence relation a n = ma n-1 + a n-2 , here with m = 3, which gives rise to the sequence 1, 3, 10, 33, 109.... While the rapid increase of this series would suggest that the resonances which arise due to the bronze ratio are typically further apart than for the gold and silver ratios, and that consequently the bronze ratio may lead to more strongly localized behaviour for small systems, in Fig. 3b) we in fact see additional dips at distances unrelated to these numbers. Presumably these are linked to other good rational approximations for the bronze ratio, however we have been unable to verify which approximation they correspond to. Nonetheless, this illustrates the point that the structure of resonances in quasiperiodic potentials is extremely complex.

Finally, in Fig. 3(d) we show the fixed-point couplings ∆ij for a random potential, with on-site energies h i = W cos(2πi/φ+θ i ) where we again set φ = (1+ √ 5)/2 but crucially in this case we choose the phase to be random on each site, θ i ∈ [-π, π]. We choose this form in order to keep the overall distribution of on-site terms the same between random and quasiperiodic potentials. The dips seen in the quasiperiodic cases are entirely absent here, and instead we observe a precise exponential decay. This confirms that the 'dips' seen previously are entirely due to the deterministic structure of resonances in quasiperiodic potentials, and that for randomly-distributed resonances these features do not appear. we have also checked with a more conventional box distribution of on-site energies (with h i ∈ [-W, W ], see Fig. 4) and found the same qualitative behaviour. It is also interesting to note that the results obtained for the LIOM interactions bear some similarity with recent strong-disorder renormalization group results on quantum spin chains [START_REF] Agrawal | Universality and quantum criticality in quasiperiodic spin chains[END_REF]. In that case it was shown that for quasiperiodic potentials generated by metallic means the couplings would flow to a self-similar binary sequence given by the respective numbers.

As shown in Fig. 4, the fixed-point interactions in the random case are typically more localized than the quasiperiodic systems, which all have an almost identical real-space decay profile. This is due to the spatially correlated nature of the quasiperiodic potentials, which enhances the probability to find resonances between nearby sites, leading to less localized behaviour. In Fig. 4 we also show results for the random potential generated using a box distribution. Here we find that for W/J > 2, it seems that Griffiths effects play a significant role, leading to the random box potential exhibiting slightly less localized behaviour, however for W/J < 2 the spatial correlations of the quasiperiodic potentials win out, leading to the quasiperiodic potentials being less localized. Taken together, the results of studying these different potentials suggest that careful real-space tuning of the resonant sites could be employed to enhance the stability of many-body localization in quasiperiodic potentials without resorting to the commonly studied case of random box disorder, which brings with it rare-region effects which may lead to enhanced delocalization.

Localization Length

As the decay of the fixed-point interactions is approximately exponential in all cases, up to some modulations due to the specific quasiperiodic potential used, we can fit this exponential decay profile to extract a localization length, defined via

∆ij med ∝ exp -2 |i -j| ξ ( 20
)
where ξ is the localization length, and the 2 in the exponent is from the convention established in studies of Anderson localization, see e.g. Ref. [START_REF] Billy | Direct observation of anderson localization of matter waves in a controlled disorder[END_REF]. Although this fitting procedure is somewhat approximate, as it ignores the large deviations in ∆ij seen when r = |i -j| corresponds to distances typical of single-particle resonances, it is nonetheless instructive.

The results are shown in Fig. 4d) for all three quasiperiodic potentials considered in this work, and the case of random disorder for comparison. We show results for the random potential generated using the cosine-like distribution with a random phase, as in the previous section, and for comparison we also show the results from a random potential generated using a box distribution. We can see that the potentials generated using the golden and silver ratios are very similar, while the potential generated from the bronze ratio is slightly more localized for all values of W/J. This confirms the conjecture mentioned in Sec. 4.2 that the more widely separated resonances generated by this potential may lead to more localized behaviour than for the other quasiperiodic potentials studied here, however this effect is very weak. Interestingly, the localization length obtained from the random box disorder is smaller than the localization length computed for any of the quasiperiodic potentials at small W/J, but is larger than the localization lengths obtained for the quasiperiodic potentials at large W/J. This illustrates the very different roles of resonances in both types of potential. For small W/J, the deterministic resonances generated by quasiperiodic potentials dominate, and they lead to delocalization -this can already be seen in the non-interacting Aubry-André model, where the proliferation of these resonances leads to a delocalized phase for W/J < 2 which does not exist in the Anderson model. On the other hand, at large W/J the near-resonant sites generated by the quasiperiodic potential are in general too far apart in energy to lead to a proliferation of resonances, and so instead we see that the effects of rare resonant regions in the case of the random box potential win out and lead to less localized behaviour in this regime.

Real-space support of LIOMs

We can also directly compute the real-space support of the local integrals of motion which characterize the fixed point Hamiltonian. Again, the diagonal Hamiltonian is given by:

H = i hi : ñi : + ij ∆ij : ñi ñj : . ( 21
)
We can start from the LIOMs in the diagonal basis, given by the local operators ñi , and applying the inverse of the unitary transform used to diagonalise the Hamiltonian we can express these operators in the original microscopic basis:

ñi = j α (i) j : n j : + j =k β (i) jk : c † j c k : + j =k γ (i) jk : n j n k : + j =k∨l =m ξ (i) jklm : c † j c k c † l c m : ... ( 22
)
which we can obtain by numerically solving the equation of motion dn i (l)/dl = [η(l), n i (l)] backwards from l → ∞ to l = 0, with the initial condition α )/2 (the bronze ratio), and (d) a random phase potential, as defined in the text. The grey and blue dashed lines show the same distances r as indicated in Fig. 3, here measured with respect to the central site on which the LIOM is defined. couplings are zero. This requires that we store the full generator at every stage of the flow, η(l), which may consume a large amount of memory for large system sizes2 . In contrast to previous work, the inclusion of higher-order off-diagonal terms in the running Hamiltonian (Eq. 15) allows us to include higher-order terms in the above ansatz for the density operator.

By looking at the prefactors, we can get some idea of how 'local' these LIOMs really are. We show the results for the lowest-order coefficients α (i) j in Fig. 5 for a LIOM on the central lattice site (i = L/2) and for quasiperiodic potentials given by the golden ratio, silver ratio, bronze ratio and finally for a random potential. Interestingly, the coefficients α (i) j also exhibit non-monotonic features at distances associated with enhanced resonances, as in Fig. 3, however in this case these distances are associated with local peaks rather than dips. This is in agreement with the results of Ref. [START_REF] Žnidarič | Interaction instability of localization in quasiperiodic systems[END_REF] who showed that resonances are enhanced at these distances: these enhanced resonances manifest here as peaks in the real-space support of the LIOMs. Interestingly, this same structure can be seen in the real-space behaviour of the eigenstates of the non-interacting Aubry-André model (not shown), suggesting that the quadratic component of the LIOMs inherits this structure from the single-particle eigenstates. As before, this structure is entirely absent in the case of a random potential, shown in Fig. 5d), where the LIOM support decays exponentially with no visible non-monotonic features.

We can also look at the coefficient of the 2-body interaction term in Eq. 22, γ

jk . We show this in Fig. 6, averaged over distance r = |j -k| similarly to the plots of ∆ij in Fig. 3. The behaviour is very similar to that of the fixed-point interaction coefficients, displaying an approximately exponential decay that becomes steeper at strong (quasi)-disorder strengths. There are hints of some weak additional structure on top of the exponential decay, as seen in )/2 (the bronze ratio), and (d) a random potential with on site energies generated using a random phase on each site. Unlike the interaction terms in the effective Hamiltonian H, the 2-body terms here do not have a clear structure on top of the exponential decay, though there are a small number of very weak features at the same distances as highlighted in Fig. 3 (dashed lines). previous quantities, however any such features are far weaker than those visible in the fixed point interaction coefficients ∆ij discussed previously. Again, this structure is absent in the case of the random phase potential, suggesting that its origin is in the deterministic structure of resonances in quasiperiodic systems.

Evidence for a delocalization phase transition

Having established that we can compute LIOMs in the real-space basis, we can now use the form of these LIOMs to get some insight as to whether there exists a delocalization transition. Although the LIOMs studied in this work are very different to the ones obtained in Ref. [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF] (there obtained using exact diagonalization to compute the non-equilibrium dynamics and extract time-averaged local quantities, and here computed by reversing the unitary transform which directly diagonalizes the Hamiltonian -see Ref. [START_REF] Chandran | Constructing local integrals of motion in the many-body localized phase[END_REF] for a discussion of this point) the key elements of the analysis can nonetheless be retained.

Following Ref. [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF], we define two ratios, f 2 and f 4 , which are the relative weights of the quadratic and quartic terms in the LIOM operator expansion. with

f 2 = j |α (i) j | 2 + jk |β (i) jk | 2 ||n|| 2 (23)
f 4 = jk |∆ (i) jk | 2 + jkpq |ξ (i) jkpq | 2 ||n|| 2 (24) 1.0 1.5 2.0 2.5 3.0 0 1 f n a) b) c) d) L = 12 L = 16 L = 24 L = 36 L = 45 L = 64 -2 0 0 1 W c = 2.40 ± 0.04 ν = 2.42 ± 0.24 1.0 1.5 2.0 2.5 3.0 W/J 0 1 f n -4 -2 0 W -W c W c L 1/ν 0 1 W c = 2.49 ± 0.44 ν = 2.19 ± 0.08
||n|| 2 = j |α (i) j | 2 + jk |β (i) jk | 2 + jk |∆ (i) jk | 2 + jkpq |ξ (i)
jkpq | 2 , e.g. f 2 is the sum of the squares of the quadratic terms divided by the sum of the squares of all terms, and f 4 is the same but for the quartic terms instead. These quantities tell us about the relative weight of the quadratic and quartic terms in the operator expansion, which is related to the role of 2-body collisions and the growth of entanglement. The ability to explicitly compute these n-body terms allows us to avoid issues with other observables, such as entanglement or transport dynamics, which may strongly reflect the underlying single-particle critical point at W c /J = 2. Here instead, we are able to explicitly calculate true many-body effects. This is a significant advancement over earlier applications of the flow equation method to disordered systems [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF] where only the quadratic terms in Eq. 22 were considered, and the following analysis was not possible.

In Fig. 7 we show both f 2 and f 4 versus disorder strength W/J for a variety of different system sizes, averaged over N s ∈ [START_REF] Wilson | Renormalization of hamiltonians[END_REF]2048] disorder realizations depending on system size. Of particular note is that both f 2 and f 4 display a strong variation with system size for small W/J values that extends up to the largest system sizes we consider here, although this variation slows significantly for the largest systems we simulate, which give virtually identical results. This implies that large system sizes are required in order to obtain accurate behaviour in the delocalized phase, which is in line with previous suggestions that large system sizes and/or simulation times are required to accurately estimate the MBL transition point [START_REF] Panda | Can we study the many-body localisation transition?[END_REF], and illustrates a key advantage of our approach, namely that it allows accurate simulations of large system sizes and provides direct access to the fixed-point Hamiltonian and the local integrals of motion without requiring the dynamical evolution or time-averaging procedure used in other works.

Following the standard procedure for phase transitions, we can make use of a finite size scaling argument to collapse all of the data onto a single curve and extract a crossing point which corresponds to the transition. As in Ref. [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF] we find good collapse of data using the scaling form Φ(L, W -W c ) = Φ(L 1/ν |W -W c |/W c ) where W c is the critical disorder strength and ν is the critical exponent, which we interpret as evidence that the system undergoes a delocalization transition at a critical value of quasiperiodic potential strength W c /J. We use the Python package pyfssa to compute the data collapse and extract these parameters [START_REF] Melchert | autoscale. py-a program for automatic finite-size scaling analyses: A user's guide[END_REF][START_REF] Sorge | [END_REF].

In Fig. 7, we show the result of the scaling analysis which (approximately) collapses the data onto a single curve and locates a transition point. The scaling of both f 2 and f 4 give results consistent within their error bars, suggesting that the transition in the interacting system has moved to a larger critical disorder strength than in the non-interacting system, as expected from previous work. The uncertainty in the critical value of W c /J resulting from the data collapse [shown in Fig. 7b) and d)]is smaller than the spacing between data points, indicating that the resolution of our dataset is the dominant source of uncertainty. The values of W c /J and ν are both consistent with the results reported in Refs. [START_REF] Zhang | Universal properties of many-body localization transitions in quasiperiodic systems[END_REF][START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF], although our lower interaction strength consequently leads to a lower value of W c /J. Interestingly, our results support the claim [START_REF] Zhang | Universal properties of many-body localization transitions in quasiperiodic systems[END_REF][START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF] that the critical exponent ν is well above the Luck bound ν ≥ 1 [START_REF] Luck | A classification of critical phenomena on quasi-crystals and other aperiodic structures[END_REF] and also fulfills the Harris criterion [START_REF] Harris | Effect of random defects on the critical behaviour of ising models[END_REF], ν ≥ 2/d, thus implying that the quasiperiodic MBL transition is perturbatively stable to the addition of weak randomness.

The central limitation of our results is the quantity of available data for the finite-size scaling, due to the time taken to diagonalize large systems using this method, and in particularly the slow approach to convergence in systems with near-degenerate single-particle eigenvalues (see Appendix B for details on this point), something which is almost guaranteed in large systems with quasiperiodic potential landscapes due to the deterministic structure of resonances which distinguishes quasiperiodic systems. Future developments in increasing the efficiency of the linear algebra operations, and in particular deploying these calculations on clusters of GPUs rather than CPUs, could allow for faster, more accurate data collected on even larger system sizes, which may result in an improved estimate of the critical point which we have tentatively identified in the present work, and may enable a more in-depth study of how this critical point moves with increasing interaction strength.

In the late stages of completion of this work, we became aware of a recent preprint (Ref. [START_REF] Aramthottil | Finite-size scaling analysis of many-body localization transition in quasi-periodic spin chains[END_REF]) which conducted a careful finite-size scaling analysis of the MBL transition in a quasiperiodic system and concluded that, for system sizes accessible to ED, the best data collapse was achieved with the assumption of a Berezinskii-Kosterlitz-Thouless (BKT) scaling form, as also seen in recent numerical studies on random systems. This result is a surprise in a quasiperiodic system, as the rare ergodic seed regions which lead to avalanaches and a BKT-like transition are not currently known to exist in deterministic potentials, and is in contrast with Refs. [START_REF] Khemani | Two universality classes for the many-body localization transition[END_REF][START_REF] Zhang | Universal properties of many-body localization transitions in quasiperiodic systems[END_REF][START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF]. The slow drift of critical disorder strength with increasing system size interpreted in Ref. [START_REF] Aramthottil | Finite-size scaling analysis of many-body localization transition in quasi-periodic spin chains[END_REF] as evidence of a BKT-like transition provides further support for the idea that finite-size effects are much weaker in quasiperiodic potentials -in this respect, further studies on large systems will need to be done in order to ascertain whether the most accurate scaling form is the logarithmic drift of system size indicative of a BKT transition, or a more generic system-size-dependent critical disorder which is also shown in Ref. [START_REF] Aramthottil | Finite-size scaling analysis of many-body localization transition in quasi-periodic spin chains[END_REF] to result in a good collapse of the data. In the present work, we interpret the weak dependence on system size for increasingly large systems as a reflection the slow growth of the number of single-particle resonances at distances r = F n (where F n is the nth Fibonacci number for φ equal to the golden ratio) with system size. Here we have assumed the same scaling form as Ref. [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF] and found broadly consistent results, however it would be interesting to revisit this assumption in light of the recent results of Ref. [START_REF] Aramthottil | Finite-size scaling analysis of many-body localization transition in quasi-periodic spin chains[END_REF]. We further suggest a connection to the work of Ref. [START_REF] Doggen | Many-body delocalization dynamics in long aubry-andré quasiperiodic chains[END_REF] who found that the localization transition seen in the non-equilibrium dynamics was strongly affected by the choice of potential: based on our results which indicate the importance of the single-particle resonances, we conjecture that the choice of quasiperiodic potential may strongly affect the non-equilibrium dynamics of finite-size systems due to the different distribution of resonances in different quasiperiodic potentials.

Conclusion

In this manuscript we have studied localization phenomena in a many-body quantum system subject to a quasiperiodic potential, using a new computational implementation of a diagonalization method based around continuous unitary transforms. We have shown that this method represents a significant improvement in accuracy over previous implementations of a truncated flow equation approach, and that the many advances incorporated in this new method allows us to access physics that was out of reach of previous approaches, in particular an estimate of the MBL transition point in a many-body system subject to a quasiperiodic potential via explicit calculation of the local integrals of motion which characterize the system. We have also shown that quasiperiodic potentials result in effective Hamiltonians which contain a strong 'fingerprint' of the structure of the underlying potential, in stark contrast to systems with random disorder which display no such clear feature.

By explicit construction of the local integrals of motion which characterize the MBL phase, we have been able to study their real-space support, gain insight as to how the struction of the LIOMs reflects the underlying potential, and have also used the fraction of quadratic and quartic components of the LIOM operator expansion to tentatively identify a phase transition between localized and ergodic phases. By numerically simulating these quantities for large system sizes, we find critical exponents which are consistent with both large-scale renormalization group studies [START_REF] Zhang | Universal properties of many-body localization transitions in quasiperiodic systems[END_REF] and recent exact numerics which constructed the LIOMs in a different manner [START_REF] Singh | Local integrals of motion and the quasiperiodic many-body localization transition[END_REF]. In particular, we find a critical exponent ν > 2/d, where d is the spatial dimension, implying the perturbative stability of the quasiperiodic MBL transition to weak randomness [START_REF] Harris | Effect of random defects on the critical behaviour of ising models[END_REF]. Our results for the largest system sizes we have studied are almost identical: as finite-size effects are known to be weaker in quasiperiodic potentials than in true disordered systems, this suggests that we are able to simulate system sizes large enough that our results are an accurate reflection of the thermodynamic limit, and that finite-size effects play only a small role here.

As with all approaches based around diagonalizing a static Hamiltonian, it is important to note that this is not unequivocal evidence for a true phase transition. The MBL transition is a dynamical phenomenon which occurs when highly excited eigenstates spontaneously fail to thermalize: here we have studied the local integrals of motion which characterize the diagonal Hamiltonian, both in terms of their real-space structure and their effective interactions, however the precise localization properties of any particular initial state will depend on its decomposition in terms of these LIOMs. In other words, a key future goal for this method will be its extension to also transform arbitrary initial states in order to express them explicitly in terms of the LIOMs computed here, in a manner which is efficient and computationally tractable.

In this work we have focused on static properties of the local integrals of motion and the corresponding fixed-point Hamiltonian H, however the method presented can be readily extended to compute the non-equilibrium dynamics following a quench, as in previous works using the flow equation method [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF], enabling the study of operator spreading [START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF] via a decomposition of a local operator similar to that used here (Eq. [START_REF] Bordia | Probing slow relaxation and many-body localization in two-dimensional quasiperiodic systems[END_REF]. It is also possible to compute correlation functions using this method [START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF], which provides an alternative way to define a localization length. The applications of this method go beyond localization in one-dimensional quasiperiodic systems. We have shown in the present manuscript that our results are valid on both sides of the localization transition, in contrast to previous truncated flow equation implementations, and as such this method is also well-suited for the study of more conventional MBL transitions in random systems. As the incorporation of higherorder terms in the ansatz of Eq. 15 is now possible in a systematic manner, it may also be possible to study more strongly interacting systems with high accuracy simply by including more terms in Eq. 15. As all flow equation methods based on Wegner-type generators result in the generation of new long-range couplings throughout the flow, there is essentially no performance cost to including long-range terms in the initial Hamiltonian, meaning that this method is also extremely well-suited to the study of localization in systems with longrange couplings, e.g. building on the results of Ref. [START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF]. The lack of explicit dependence on the geometry of the system similarly means that an L = 64 chain, for example, can be diagonalized with essentially the same complexity as an L = 8 × 8 two-dimensional system, facilitating the construction of local integrals of motion in two-dimensional systems. The method is not restricted to fermions either, and can be generalized to bosonic systems simply by taking into account bosonic rather than fermionic commutation relations: this may allow an advantage over tensor network based methods where the size of the local bosonic Hilbert space features explicitly in the construction of matrix product state wavefunctions, which can quickly become prohibitive for weakly-interacting systems where large local occupations are possible. Here, by contrast, the dimension of the local Hilbert space does not enter explicitly in the diagonalization nor the construction of the LIOMs.

The biggest computational advantage of this method, however, is the significant speed increase that can be obtained from using graphical processing units (GPUs) rather than the more conventional CPUs used in the majority of current numerical works in many-body physics. As we demonstrate in Appendix C, for large systems even a modest GPU can result in a significant speed increase. In the present manuscript, we have taken only the first steps in this direction: we expect that the continued development of this method on cutting-edge GPU hardware will enable its extension into parameter regimes which are currently inaccessible to all other theoretical methods. For completeness, and to motivate others to adopt flow equation methods, we also present a small sample code in Appendix D illustrating the application of this technique to a non-interacting system. icity Breaking in Quantum Matter), as well as support from the ANR grant "NonEQuMat" (ANR-19-CE47-0001). SJT also acknowledges support from an NVIDIA Academic Hardware Grant. Initial computations were performed on the Collège de France IPH-JE computing cluster. We gratefully acknowledge useful discussions with B. Ware and feedback on the manuscript from L. Sanchez-Palencia, as well as valuable comments from C. Bertoni. Code and data are available at [START_REF] Thomson | PyFlow: A python package for diagonalising quantum systems using continuous unitary transforms[END_REF]92].

A Consistent ordering of operators

In the main text, we have assumed a consistent order of fermionic operators in the Hamiltonian, however we did not go into detail about how this is achieved. In general swapping the order of two fermionic operators can result in the generation of new terms, e.g.:

c † i c j c † k c q = c † i (c † k c j + δ jk )c q = c † i c † k c j c q + δ jk c † i c q (25)
Therefore if this term arose when evaluating dH/dl, one could arrange the operators in the form c † i c j c † k c q and conclude that this term renormalizes the quartic terms only, or one could rearrange the operators into the form c † i c † k c j c q + δ jk c † i c q and conclude that this term renormalized both quartic and quadratic parts of the Hamiltonian. While not necessarily inconsistent, this arbitrariness complicates the computation of high-order commutators. The key ingredient in avoiding this ambiguity is known as normal-ordering, which enforces a consistent ordering of the operators in each term. A detailed description of normal ordering in the context of unitary flow methods is given in Refs. [START_REF] Kehrein | The flow equation approach to many-particle systems[END_REF][START_REF] Wegner | Flow equations and normal ordering: a survey[END_REF], and we briefly summarize the main points here. While this prescription is often thought of as simply subtracting the expectation value from an operator, the fundamental idea is rather more sophisticated. For fermions, it consists of expressing a fermionic anticommutation relation as:

{c † i , c j } = G ij + Gji ( 26
)
where we define the contractions used in normal-ordering as:

G ij = c † i c j (27) Gji = c j c † i (28)
and the expectation values are computed in some reference state, which is typically the ground state of the corresponding non-interacting model, but can also be a more complex mixed state described by some density matrix. To calculate the commutation relations of normal-ordered strings of operators, we use the following theorem [START_REF] Kehrein | The flow equation approach to many-particle systems[END_REF]:

: O 1 (A) :: O 2 (A ) :=: exp kl G kl ∂ 2 ∂A l ∂A k O 1 (A)O 2 (A ) : (29)
where A and A are the set of labelled operators in the expression O 1 and O 2 which in our case are just strings of fermionic operators. We can evaluate this by expanding the exponential -for simple expressions, the expansion can be exact, as all derivatives are zero above some order greater than the number of operators in O 1 or O 2 .

The precise choice of reference state in which the contractions are computed depends on the calculation in question: in Refs. [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF], we used a uniformly half-filled product state with n i = 0.5∀i to approximate an arbitrary state of high energy density, although in the framework of that calculation the precise state we chose made little difference: incorporation of normal ordering corrections was only necessary in order to prevent divergent interacting terms at small disorder strengths, as in Refs. [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] the main role of normal-ordering corrections was to terminate the flow early if interaction terms began to diverge. Here, however, due to the more sophisticated form of running Hamiltonian which we employ, these divergences are no longer an issue, and we are free to choose a simpler form of normal ordering, which we have verified the accuracy of via a series of numerical checks B. Here we use vacuum normal ordering, commonly thought of as 'moving all dagger operators to the left' in any composite string of fermionic operators. Vacuum normal ordering implies that G ij = 0 and Gji = δ ij , such that we still have the general identity G ij + Gji = δ ij that satisfies fermionic anticommutation relations.

A.1 Example: Quadratic terms

: c † α c β :: c † γ c δ : =: 1 + G αδ ∂ 2 ∂c † α ∂c δ + Gβγ ∂ 2 ∂c β ∂c † γ + G αδ Gβγ ∂ 4 ∂c † α ∂c δ ∂c β ∂c † γ c † α c β c † γ c δ : (30) =: c † α c β c † γ c δ : +G αδ : c β c † γ : + Gβγ : c † α c δ : +G αδ Gβγ (31)
So the commutation relation is given by:

[: c † α c β :, : c † γ c δ :] = -(G αδ + Gαδ ) : c β c † γ : +( Gβγ + G γβ ) : c † α c δ : +(G αδ Gβγ -G γβ Gδα ) (32) = δ βγ : c † α c δ : -δ αδ : c † γ c β : +(G αδ Gβγ -G γβ Gδα ) (33)
which, after applying the vacuum normal-ordering identity G ij = 0, reduces back to the same result as one would obtain without normal-ordering.

A.2 Example: Higher-order terms

The higher-order terms proceed in essentially the same way, giving:

[: c † α c β c † γ c δ :, : c † µ c ν :] = -(G αν + Gνα ) : c † µ c β c † γ c δ : -(G γν + Gνγ ) : c † α c β c † µ c δ : + ( Gβµ + G µβ ) : c † α c ν c † γ c δ : +( Gδµ + G µδ ) : c † α c β c † γ c ν : + (G αν Gβµ -G µβ Gνα ) : c † γ c δ : +(G αν Gδµ -G µδ Gνα ) : c † γ c β : + (G γν Gβµ -G µβ Gνγ ) : c † α c δ : +(G γν Gδµ -G µδ Gνγ ) : c † α c β : (34)
After applying vacuum normal-ordering identities for the contractions, this results in:

[: c † α c β c † γ c δ :, : c † µ c ν :] = -δ αν : c † µ c β c † γ c δ : -δ γν : c † α c β c † µ c δ : +δ βµ : c † α c ν c † γ c δ : +δ δµ : c † α c β c † γ c ν :
which is exactly the same result as obtained from the (much simpler) graphical notation introduced in Sec. 3.2. Even higher order terms may be computed in similar ways, though we do not go into further detail here.

B Checks B.1 Comparison with ED

As a first check of our method, for small systems we benchmarked the many-body eigenvalues obtained using the flow equation (FE) approach against numerically exact diagonalization (ED) using the QuSpin library [START_REF] Weinberg | QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems part I: spin chains[END_REF][START_REF] Weinberg | QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems. Part II: bosons, fermions and higher spins[END_REF]. Here we present some of these results for small systems in order to demonstrate the accuracy of our method. In the diagonal l → ∞ basis, the manybody eigenstates are simply product states of LIOMs, e.g. |ψ = |0001, 0010, 0011... , and can be straightforwardly constructed by applying the Hamiltonian H to all 2 L possible product states. Note that this is a very demanding check of our method, as it involves all many-body eigenstates, not simply the half-filled states or those in the middle of the excitation spectrum. Once the eigenvalues E n are obtained, the relative error in the n-th eigenvalue is defined as:

n = E ED n -E F E n E ED n ( 35
)
where the superscripts ED and F E refer to the eigenvalues obtained with exact diagonalization and flow equation methods respectively. The results are shown in Fig. 8 for system size L = 12 and three different interaction strengths, each in a quasiperiodic potential generated using the golden ratio and averaged over N s = 2048 values of random phase θ. In the localised phase, the median relative error is extremely small, but diverges at low values of W/J when the system becomes delocalised. We note that the relative error computed for the same system using the older method of Refs. [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] (dashed lines in Fig. 8) is larger in the localised phase, but smaller in the delocalised phase: this is a combination of the normal-ordering procedure used in Ref. [START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] to prevent divergences, and the fact that fewer terms were kept in the older implementation, resulting in fewer divergent terms and correspondingly a smaller error. W/J = 1.00 W/J = 2.00 W/J = 3.00 )/2 (the golden ratio), a single value of the phase (θ = 0). Note that although the data for small values of W/J may appear flat on this scale, they are in fact decaying exponentially with a shallow gradient. (b) Fraction of off-diagonal couplings at each stage in the flow which are above the threshold cutoff of 10 -3 . There is an initial increase as new couplings are generated in the first few flow time steps, followed by a sharp decrease at small flow times, and then a slow decay towards zero.

B.2 Convergence checks

To evaluate the performance of diagonalization methods based on continuous unitary transforms, one key metric is the how quickly the transform converges. For the Wegner generator used in the present manuscript, the convergence is determined by how widely-separated the single-particle eigenvalues are, with more widely separated eigenvalues leading to faster convergence, which implies that in general the convergence will be faster for larger quasi-disorder amplitudes -see Refs. [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF][START_REF] Kehrein | The flow equation approach to many-particle systems[END_REF][START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF] for further discussion of this point. We first examine the convergence properties of the non-interacting Aubry-André model before looking at the convergence of the interacting system given in Eq. 1.

To give some idea of the flow times required in order to obtain the eigenvalues to a given accuracy, in Fig. 9a) we show the relative error of the on-site terms of a non-interacting system, computed with respect to the exact single-particle eigenvalues obtained by directly diagonalizing the matrix, while in Fig. 9b) we show the fraction F of off-diagonal matrix elements which are above the cutoff of 10 -3 .

Here we use a system size L = 64 and three different quasi-disorder strengths, with each quasi-disorder potential representing a single random realization obtained using the golden ratio. We integrated to a maximum flow time l max = 10 3 , but stop the flow early if all offdiagonal elements have decayed below the cutoff (as seen for W/J = 3). For all values of W/J the error exhibits a sharp decrease at the beginning of the flow, with larger values of W/J typically decaying more quickly, before exhibiting a slow exponential decay towards zero in the later states of the flow. For quasi-disorder realizations with near-degenerate eigenvalues, this flow can be extremely slow and require very large flow times in order to obtain convergence to a high accuracy. We note that the late-time convergence is slowest around the critical point of the non-interacting system W/J = 2 where the spectrum is fractal in nature [START_REF] Yao | Critical Behavior and Fractality in Shallow One-Dimensional Quasiperiodic Potentials[END_REF] and therefore the eigenvalues are very close together. As we are primarily interested in the range W/J > 2, )/2 (the golden ratio), one single value of the phase (θ = 0), and several different values of the cutoff ε. In each case, the flow is run until all off-diagonal elements have decayed to below the cutoff value. We see that for large cutoffs the ∆ij terms decay approximately exponentially with distance, but when smaller cutoffs are used, the characteristic 'dips' at distances r = F n appear, where the F n are the Fibonacci numbers. This confirm that these features are not convergence errors.

in the main text we used a maximum flow time of l max = 25, which we found to be a good compromise between speed and accuracy, with a typical error in the single-particle eigenvalues of around 10 -2 . Beyond this point we see diminishing returns by going to larger flow times, and in any case relative error in the presence of interactions is primarily limited by the truncated form of Eq. 15 and the strength of the microscopic interactions rather than the maximum flow time, so we do not find larger flow times to lead to a meaningful increase in accuracy.

We now turn to the convergence of the interaction terms. By comparison with the findings of Ref. [START_REF] Žnidarič | Interaction instability of localization in quasiperiodic systems[END_REF] who noted that single-particle resonances are enhanced in quasiperiodic systems generated using the golden mean at distances corresponding to the Fibonacci numbers F n , one might ask whether the 'dips' seen in the fixed-point couplings ∆ij discussed in Sec. 4.2 could be due to convergence errors, i.e. the associated hopping terms J ij may not have decayed sufficiently under the action of the transform, leading to an incorrect result for ∆ij at these distances. In Fig. 10 we show that this is not the case, and that on the contrary we do not see these dips until all couplings have sufficiently decayed, here for approximately ε ≤ 10 -2 of their initial value. For still smaller cutoffs there are some quantitative changes, however these features persist. In the results presented in the main manuscript, we use a flow time such that the majority of couplings have decayed below ε = 10 -3 : we do not insist that all couplings decay below this threshold, as in this case our simulations become limited by unusual realizations of the quasperiodic potential where near-degeneracies can cause a small number of couplings to flow extremely slowly. Loosely speaking, requiring all off-diagonal elements to decay by another order of magnitude requires a flow time (and therefore computation time) which is itself an order of magnitude larger, and this quickly becomes prohibitive once we go beyond ε = 10 -3 . While it may be the case that additional effects not taken into account in the present calculation, e.g. higher-order terms or an alternative form of normal-ordering, could modify our results, the features seen in the fixed-point couplings ∆ seem robust and it compared with the method presented in this manuscript. Here we reproduce the results from Fig. 3a) [solid lines], and compare with the same quantity computed using the older flow equation implementation [dashed lines]. We can clearly see that while the broad features are the same, the new method leads to significantly modified short-range behaviour, and consequently far more localized fixed-point interactions ∆ij .

is hard to anticipate any modifications that would destroy them.

B.3 Comparison with previous implementations of flow equations

An earlier implementation of continuous unitary transforms for the study of MBL phenomena was developed in Refs. [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] using a more approximate form of the running Hamiltonian which allowed us to analytically obtain flow equations for the running couplings, but which by construction did not allow for the generation of off-diagonal quartic terms during the flow -we refer the reader to Refs. [START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Quasi-many-body localization of interacting fermions with long-range couplings[END_REF] for further details regarding this earlier implementation. The numerous and varied improvements presented in this manuscript lead to a significant improvement over the older method, which we shall demonstrate here. In Fig. 11 we show the results for the fixed-point couplings ∆ij obtained using both the old and new flow equation implementations, both for the same system size of L = 36 and averaged over N s = 128 realizations of the quasiperiodic potential generated using the golden ratio, as in Fig. 3a). There is a clear difference between the results, particularly at short distances. Due to the inclusion of higher-order off-diagonal terms, the short-range behaviour is captured much more accurately by the new implementation: in the previous implementation, the short-range couplings can exhibit divergences at small values of the quasiperiodic potential amplitude (i.e. in the delocalized phase). The new implementation does not exhibit the same problem, and this change in the short-range behaviour leads to a noticeably faster decay of the couplings, corresponding to a 'more local' unitary transform. Nonetheless, at longer distances many of the main features present in the results obtained using the new method are also present in the older implementation, notably the 'dips' at distances r = F n where the F n are the Fibonacci numbers.

C CPU/GPU Speed Comparison

All comparisons in this Appendix were performed on a standard laptop, a Dell G5 5590 with a 6-core (12-thread) Intel Core i7 9750-H CPU with 16Gb RAM and a 6Gb NVIDIA GTX1660 Ti GPU. Even with modest hardware, the speed increase from making use of the high degree of parallelization possible with a GPU can be significant. These comparisons made use of the PyTorch library [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] and the torchdiffeq package [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Chen | Learning neural event functions for ordinary differential equations[END_REF] to facilitate the straightforward use of sophisticated integration routines on a GPU: to facilitate a like-for-like comparison, precisely the same code was run on the CPU by setting device='cpu' in the PyTorch settings. This is in contrast to the code used in the main part of this manuscript which used custom-coded tensor contraction routines using Numba's just-in-time (JIT) compilation feature to achieve the fastest possible performance [START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF]. Here we use a flow time l max large enough that all couplings have decayed below the threshold of ε/J = 10 -3 , and all code uses single-precision floats. We used a quasiperiodic potential generated using the golden ratio, and fixed the phase to be θ = 0 for all simulations. The diagonalization for each system size was run several times (except for L = 24, for which we conducted only one run due to the long time required on the CPU) and then we took the average value of the total time taken: error bars showing the standard deviation over different runs are smaller than the plot markers. These benchmarks show the time taken to fully diagonalize the Hamiltonian and to transform a single local operator from the initial basis to the final diagonal basis. The time shown is the full execution time of the script, measured using the dateTime module.

As shown in Fig. 12, the relative speed of the GPU over the CPU increases sharply with system size, with simulations for L = 24 which took approximately 2 hours on the CPU taking around 15 minutes on the GPU. We expect this trend to continue to larger system sizes, with GPUs exhibiting a clear performance advantage for the diagonalization of large many-body systems, and perhaps enabling further studies of MBL phenomena in two-dimensional models where diagonalization with a CPU would be unfeasible. The limiting factor in the utility of GPUs is their relatively small memory, which becomes an issue if one wishes to store the full unitary transform (typically O(10 3 × L 4 ) floating point numbers, which can rapidly reach 10 2 Gb or more in size) in order to be able to reverse the transform from the diagonal basis back into the initial microscopic basis: to achieve this on large system sizes, it may be necessary to develop more sophisticated algorithms which can run on GPU clusters and take advantage of their shared memory, or develop a more efficient way of passing the unitary transform data from the GPU to the CPU during the initial diagonalization, and quickly re-loading it back to the GPU when applying the reverse transform.

D Sample Code for Non-interacting Systems

As a simple, non-optimized example to demonstrate the applicability of this method, here we provide a sample Python code for diagonalising a non-interacting system using this technique. Here we prioritize readability over speed, and so the following is adapted from the code used in the main manuscript to be a short, self-contained example with most of the sophisticated features stripped out.

First we import the required modules [START_REF] Harris | Array programming with NumPy[END_REF][START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF][START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF] and set the Hamiltonian parameters: and GPU (orange). As the system size increases, the relative speed-up possible with a GPU increases significantly.

2 from scipy . integrate import odeint 3 import matplotlib . pyplot as plt For the purposes of this example, we chose an l max that should be sufficiently large to ensure convergence to the diagonal Hamiltonian, however in the code used to obtain our main results we simply set a large value of l max and stop the flow once all off-diagonal couplings have decayed below some threshold value which we choose to be small. Next, we initialize the Hamiltonian by defining matrices for the on-site terms H 0 and the off-diagonal terms V :

1 # Non -interacting matrices 2 H0 = np . zeros (( n , n ) , dtype = np . float64 ) 3 if dis_type == ' random ': This should produce an output similar to that shown in Fig. 13. We invite the interested reader to experiment with this sample code, testing different potentials and system sizes to gain a feel for how this procedure works.

Figure 1 :

 1 Figure1: Sample flows of the coefficients in Eq. 15 for a small system of size L = 4, with J = 1 as the unit of energy, quasi-disorder strength W/J = 1 and nearest-neighbour interaction strength ∆/J = 0.1, with incommensurate ratio φ = (1 + √ 5)/2 (the golden ratio) and phase θ = 0. Panel (a) shows the flow of the quadratic terms, H

Figure 2 :

 2 Figure2: Conservation of the flow invariant, as defined in Eq. 14. Here we take the median value of δI 2 , over N s ∈[START_REF] Wilson | Renormalization of hamiltonians[END_REF] 2048] samples, and the error bars show the median absolute deviation. The deviation from unitarity is extremely small in the localized phase where δI 2 appears independent of system size. Deep in the delocalized phase (W/J 2), there is a slow upwards drift of δI 2 with increasing system size.

  = 1.00 W/J = 1.20 W/J = 1.40 W/J = 1.60 W/J = 1.80 W/J = 2.00 W/J = 2.20 W/J = 2.40 W/J = 2.60 W/J = 2.80 W/J = 3.00

Figure 3 :

 3 Figure 3: Comparison of the fixed-point Hamiltonian interaction terms for a variety of different potentials: (a) a quasiperiodic potential generated using φ = (1 + √ 5)/2 (the golden ratio), (b) with φ = 1 + √ 2 (the silver ratio), (c) with φ = (3 + √ 13)/2 (the bronze ratio), and (d) a random potential, as defined in the text. All three quasiperiodic potentials exhibit approximately exponential decay, but with sharp dips at distances r = |i -j| corresponding to numbers from their respective series expansions (grey and blue dashed lines, detailed in the main text). These dips are absent entirely in the case of a random potential, illustrating the how the different structure of the resonances in quasiperiodic systems leads to qualitatively different behaviour from random disorder.

Figure 4 :

 4 Figure 4: (a-c) Comparison of different potentials at fixed disorder strength. In each case, all three quasiperiodic potentials behave almost indistinguishably (except for the features discussed in Sec. 4.2), the random phase potential (black dashed line) is typically more localized than the quasperiodic potentials, and the random box potential (purple dashed line) exhibits clear differences for W/J < 2 and W/J > 2. (d) Localization length computed from each of the different potentials: all follow a qualitatively similar pattern, with the results from the quasiperiodic potential obtained using the bronze ratio displaying a marginally more localized behaviour than the other two quasiperiodic potentials, and the random box system being more localized than all quasiperiodic potentials at small W/J but less localized at large W/J.

j

  (l → ∞) = δ ij and all other

Figure 5 :

 5 Figure 5: Real space support of a LIOM in the centre of a chain of length L = 36, averaged over N s = 128 disorder realizations. Similarly to Fig. 3, the panels show (a) a quasiperiodic potential generated using φ = (1 + √ 5)/2 (the golden ratio), (b) with φ = 1 + √ 2 (the silver ratio), (c) with φ = (3 + √ 13)/2 (the bronze ratio), and (d) a random phase potential, as defined in the text. The grey and blue dashed lines show the same distances r as indicated in Fig. 3, here measured with respect to the central site on which the LIOM is defined.

Figure 6 :

 6 Figure 6: Decay of the 2-body interaction term γ (i) jk in Eq. 22 for the central site i = L/2 of a system with size L = 36, averaged over distance r = |j -k|. As in Fig. 3, panel (a) shows results for φ = (1 + √ 5)/2 (the golden ratio), (b) with φ = 1 + √ 2 (the silver ratio), (c) with φ = (3 + √ 13)/2 (the bronze ratio), and (d) a random potential with on site energies generated using a random phase on each site. Unlike the interaction terms in the effective Hamiltonian H, the 2-body terms here do not have a clear structure on top of the exponential decay, though there are a small number of very weak features at the same distances as highlighted in Fig.3(dashed lines).

Figure 7 :

 7 Figure7: Relative weights of the quadratic (f 2 ) and quartic (f 4 ) terms in the LIOM operator expansion. Panels (a) and (c) show the raw data, including the error bars marking the standard error used as input to the pyfssa package. Note that the data for the largest system sizes L = 36, 48, 64 are almost on top of one another, suggesting that these system sizes may be large enough to observe the behaviour in the thermodynamic limit. Panels (b) and (d) show the rescaled data, exhibiting good collapse in the region close to W c /J ≈ 2.4.

10 Figure 8 :

 108 Figure 8: Median relative error in the many-body eigenenergies obtained from the flow equation method for L = 12 and three different microscopic interaction strengths ∆ 0 , using φ = (1 + √ 5)/2 (the golden ratio) and averaged over N s = 2048 random values of the phase θ, computed with respect to the results obtained from exact diagonalization. The error bars represent the median average deviation over different random values of the phase θ. The dashed lines represent the result of the older implementation used in Refs [42, 43].

Figure 9 :

 9 Figure9: (a) Mean relative error in the diagonal matrix elements of the noninteracting Aubry-André model throughout the flow, computed with respect to the exact single-particle eigenvalues for a system of size L = 64 for φ = (1 + √ 5)/2 (the golden ratio), a single value of the phase (θ = 0). Note that although the data for small values of W/J may appear flat on this scale, they are in fact decaying exponentially with a shallow gradient. (b) Fraction of off-diagonal couplings at each stage in the flow which are above the threshold cutoff of 10 -3 . There is an initial increase as new couplings are generated in the first few flow time steps, followed by a sharp decrease at small flow times, and then a slow decay towards zero.

- 4 Figure 10 :

 410 Figure10: Convergence of the interaction terms in the running Hamiltonian, ∆ij , shown for a system of size L = 16, φ = (1 + √ 5)/2 (the golden ratio), one single value of the phase (θ = 0), and several different values of the cutoff ε. In each case, the flow is run until all off-diagonal elements have decayed to below the cutoff value. We see that for large cutoffs the ∆ij terms decay approximately exponentially with distance, but when smaller cutoffs are used, the characteristic 'dips' at distances r = F n appear, where the F n are the Fibonacci numbers. This confirm that these features are not convergence errors.

Figure 11 :

 11 Figure 11: Comparison of the truncated flow equation method from Refs.[START_REF] Thomson | Time evolution of many-body localized systems with the flow equation approach[END_REF][START_REF] Thomson | Dynamics of disordered quantum systems using flow equations[END_REF] compared with the method presented in this manuscript. Here we reproduce the results from Fig.3a) [solid lines], and compare with the same quantity computed using the older flow equation implementation [dashed lines]. We can clearly see that while the broad features are the same, the new method leads to significantly modified short-range behaviour, and consequently far more localized fixed-point interactions ∆ij .

Figure 12 :

 12 Figure12: Comparison of time taken to diagonalize a system of size L for CPU (blue) and GPU (orange). As the system size increases, the relative speed-up possible with a GPU increases significantly.

4 5 12 #

 412 dis_type = ' QP ' # Specify potential : ' random ' or ' QP ' Generate list of flow time steps to store 13 dl_list = np . linspace (0 ,10 ,100 , endpoint = True )

4 for i in range ( n ) : 5 # 11 #Figure 13 : 6 #

 4511136 Figure 13: A sample flow of a non-interacting Hamiltonian with L = 4, W/J = 2 and J = 1 using the code described in Appendix. D.

While the Hamiltonian used here is very general, one must however take care that the generator of the unitary transform is non-zero: for translationally invariant systems, the Wegner generator which we use in this manuscript vanishes and the unitary transform we employ here reduces to the identity. In such cases, one must use a different generator or work in momentum space rather than real space.

In principle, one could recompute η(l) = [H0, V ] on the fly starting from H(l → ∞) and avoid having to store the full transform, however as the new initial conditions for all off-diagonal couplings are all zero, this equation is numerically unstable.

import numpy as np

We note that it is not so straightforward to run Cython-compiled code on a GPU, and as such we opted to focus on writing code which was largely CPU/GPU-agnostic using Numba.

We use odeint here as it results in very simple code: for the results in the main text, we used ode with the dopri5 integrator.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No. 101031489 (Ergod-

Now we can define a function to compute the commutator of two matrices, [A, B] = AB -BA. Here we use NumPy's einsum routine as its inner workings are quite transparent due to the explicit labelling of indices, however tensordot is also a good (and often faster) choice, and the best performance can be found by using Numba or Cython to write compiled code which can be efficiently parallelized 3 :

1 # Function to contract square matrices ( matrix multiplication ) 2 def comm (A , B ) :

We can also define a function to compute the generator η = [H 0 , V ] and then the right hand side of the flow equation dH/dl = [η, H]:

1 # Function to compute the RHS of the flow equation 2 def nonint_ode (H , l ) : And finally we pass this function to the integrator odeint and integrate the differential equation for the Hamiltonian up to our chosen l max value 4 . Note that odeint does not accept matrix inputs, so we first flatten the n × n matrix into a list of length n 2 .

1 sol = odeint ( nonint_ode ,( H0 + V0 ) . reshape ( n **2) , dl_list )

Now we have solved the problem and diagonalized the Hamiltonian. We can compare the results with the exact solution, which here can be computed straightforwardly using NumPy's eigenvalue function:

2 print ( ' Flow eigenvalues ' , eig ) 3 print ( ' NumPy eigenvalues ' , np . sort ( np . linalg . eigvalsh ( H0 + V0 ) ) )

This will produce an output of the form: Note that by the nature of (quasi)-random systems, the result will be different each time the code is run. Depending on the system size and the precise potential generated, a longer flow time may be required in order for the results to properly converge. We can also plot the flow of the Hamiltonian parameters to visualize how they behave during the diagonalization process:

1 sol = sol . reshape ( len ( sol ) ,n , n ) 2 for i in range ( n ) :