Unraveling long-time quantum dynamics using flow equations

The study of many-body quantum dynamics in strongly-correlated systems is extremely challenging. To date few numerical methods exist which are capable of simulating the non-equilibrium dynamics of two-dimensional quantum systems, in part reflecting complexity theoretic obstructions. In this work, we present a new technique able to overcome this obstacle, by combining continuous unitary flow techniques with the newly developed method of scrambling transforms. We overcome the prejudice that approximately diagonalizing the Hamiltonian cannot lead to reliable predictions for relatively long times. To the contrary, we show that the method works well in both localized and delocalized phases, and makes reliable predictions for a number of quantities including infinite-temperature autocorrelation functions. We complement our findings with rigorous incremental bounds on the truncation error. This approach shows that in practice, the exploration of intermediate-scale time evolution may be more feasible than is commonly assumed, challenging near-term quantum simulators.

Taming the exponential complexity of many-body quantum systems remains one of the biggest challenges in modern physics.Exact numerical simulations provide the gold standard in accuracy, however, the computational cost quickly becomes prohibitive above a few tens of particles, and even rapid developments in computing power cannot outpace the exponential scaling of the complexity of fully solving a many-body quantum system.While there are efficient methods able to estimate the ground states of various quantum systems captured by local Hamiltonians -including tensor network and quantum Monte Carlo techniques -the issue of complexity becomes even more of an obstacle for time evolution.Time evolution of a given quantum state under the action of a local Hamiltonian is BQP complete in worst case complexity.For this reason, one cannot hope to find universal classical methods that can accurately and efficiently simulate this evolution for all time and all local Hamiltonians [1].While the ultimate goal may be the development of flexible and reliable quantum simulators [2][3][4] able to directly realize many models of interest in the near-term we must continue to rely upon classical computers in order to simulate quantum matter.
To that end, many highly effective numerical techniques have been developed to study many-body quantum systems subject to controlled and clear approximations.Leading the charge are tensor network methods [5,6], instances of variational methods that build on tensor networks, particularly matrix product state (MPS) approaches in one dimension and projected entangled pair states (PEPS) techniques in two dimensions.These methods work well for ground states and short time evolution, but are limited in the way they can capture dynamics, a state of affairs sometimes dubbed the 'entanglement barrier'.This core limitation stems from the generation of entanglement, as highly entangled systems require large bond dimensions, giving rise to computationally intractable situations.Quantum Monte Carlo techniques [7] are also widely usedincluding for non-equilibrium dynamics [8,9] -however, they suffer from the well-known sign problem and stability issues.Dynamical mean-field theory can also capture quantum dynamics [10,11], but again stability matters arise.These obstacles all reflect the computational hardness of the task.
In this work, we develop a radically different approach to addressing the issue of time evolution in closed quantum systems.
+ l q 0 5 J 5 s 5 h T 9 w P n 8 A i X y N S A = = < / l a t e x i t > (b)  Combining the established method of continuous unitary transforms (CUTs) -also known as 'flow equations' [12][13][14][15][16][17][18][19] -with the newly developed method of scrambling transforms, we present a flexible and powerful approach to diagonalizing large Hamiltonians and computing time evolution to very long times.The key ingredient in our work is the use of scrambling transforms to improve the convergence properties of CUT-based methods, significantly improving their accuracy and validity.We demonstrate the potential of this technique by computing the dynamics of disordered quantum systems in one and two dimensions.The limitation is a very different one compared to tensor network approaches: here, it is not entanglement that provides a limitation, but the accuracy of the approximate transform used.Heuristically, we find that in practice this restriction is less severe than overcoming the entanglement barrier.
We will focus on a generic system of interacting fermions, given by the Hamiltonian  where : ... : represents normal-ordering with respect to the vacuum, and |L| =: L is the system size in the total number of modes.We make no assumptions as to the form of the couplings, nor the dimensionality of the system: The complexity of the calculation is set by the total number of lattice sites L, not by their geometry or size of the local Hilbert spaces.A twodimensional system can be unfolded onto a one-dimensional system with long-range hopping, as sketched in Fig. 2, which does not pose a problem for CUT-based techniques.A similar procedure can be followed in three dimensions.Flow equation methods diagonalise the Hamiltonian by successively applying infinitesimal unitary transforms dU (l) = exp(−η(l)dl) = 1 − η(l)dl, where η(l) is the generator and l represents a fictitious "flow time" such that l = 0 is the initial Hamiltonian, and the parametrized Hamiltonian H(l) := U † (l)HU (l) becomes diagonal in the limit l → ∞, where the full unitary transform U (l) = T l exp(− l 0 η(l ′ )dl ′ ) is a time-ordered integral over flow time l.The diagonalization procedure can be recast as solving the equation of motion dH/dl = [η(l), H(l)] [17,18].We store H (2) as a matrix with O(L 2 ) entries and H (4) as a tensor of order four with O(L 4 ) real entries, and employ a similar procedure for the generator η(l) =: η (2) (l) + η (4) (l).This allows the relevant commutators to be computed efficiently as the sum of all one-point contractions of pairs of matrices or tensors [20].The main consequence of fermionic statistics is the minus signs which arise when computing these contractions; the method can be applied to bosons with minor changes.
A common choice of generator is η(l) := [H 0 , V (l)], where H 0 (l) and V (l) are, respectively, the diagonal and off-diagonal parts of the Hamiltonian.In what follows, we will make use of the symbol V whenever referring to offdiagonal elements.This is often known as the Wegner generator [17,18].The diagonalisation can be seen from the fact that the squared ∥V (l)∥ 2  2 is non-increasing in the fictitious time l as d∥V (l)∥ 2  2 /dl = −2∥η(l)∥ 2 2 ≤ 0 (see, e.g., Ref. [21]).Convergence relies upon the model in question having a clear separation of energy scales in the initial basis.Models where this is not true -such as homogeneous systems, or disordered systems with many near-degeneracies (known as resonances) -cannot be fully diagonalised by this generator, as they act like unstable fixed points.Perturbing the Hamiltonian away from this fixed point can allow the flow to begin, however, small perturbations can result in long convergence times while large perturbations improve convergence but risk changing Dashed black lines close to the origin in panels (a) and (b) are fits with the form χ ∝ 1/L 3 , valid for large systems and strong disorder.Insets show the un-normalised complexity (i.e., the numerator of Eq. ( 4)), which tends to a constant in strongly-disordered one-dimensional chains, but grows in two dimensions even for strong disorder.the underlying physics.Here, we resolve this by introducing scrambling transforms, which are targeted unitary transforms aimed at lifting degeneracies which the Wegner procedure alone is unable to resolve.As they are unitary, they cannot change the underlying physics: they simply act to 'prepare' the Hamiltonian in a basis more amenable to being diagonalized by the conventional Wegner flow.The procedure is sketched in Fig. 1.The (infinitesimal) scrambling transform takes the form dS(l) = exp(−λ(l)dl), with a generator λ(l) given by with δh = ε|h i (l) − h j (l)|, where ε > 0 is the threshold parameter which controls how easily the scrambling transform triggers.For ε = 0, this reduces to the Toda-Mielke generator [21,22].Here, we use ε = 0.5.The full scrambling transform S(l) can be written as a time-ordered integral over dS(l).It is employed at the beginning of the flow, and during the diagonalization procedure if degeneracies are encountered.
The particular scrambling transform used here is quadratic and does not induce any new higher order terms, however, the action of the Wegner generator will typically lead to the generation of new terms containing six or more fermionic operators, similar to the way that such terms can arise in renormalization group procedures.The central approximation of the CUT  technique is that the Hamiltonian must be truncated, and terms above a certain order neglected.We shall present rigorous error bounds later; for the moment, we emphasize that in cases where the method is insufficiently exact, higher order terms can be systematically included until the desired precision is reached, at a cost polynomial in the system size.We will investigate initial local Hamiltonians of the form ], using open boundary conditions, with J = 1 and ∆ 0 = 0.1.In one dimension, this Hamiltonian maps onto the XXZ chain via a Jordan-Wigner transform.We diagonalise these Hamiltonians in both one and two dimensions, for two different choices of h i : random disorder (h i ∈ [−d, d]) and quasi-periodic (QP) potentials (in one spatial dimension h i := d cos(2πi/ϕ + θ), with ϕ := (1 + √ 5)/2 and θ a (real) randomly chosen phase that plays the role of a 'disorder realization'; in two dimensions given by h i := d(cos(2πi x /ϕ+θ)+cos(2πi y /ϕ 2 +θ 2 )), where (i x , i y ) represent the coordinates of lattice site i, ϕ 2 = 1 + √ 2 and the θ 2 is another random phase).For simplicity, we shall refer to d as the 'disorder strength' in both cases.The end point is an (approximately) diagonal Hamiltonian where R represents neglected higher-order terms, typically of order O(∆ 2 0 ) and higher.The interaction coefficients have been shown to decay exponentially with distance in (quasi)disordered systems, ∆ ij ∝ e −|i−j|/ξ [20,[23][24][25][26][27].
Once the Hamiltonian has been diagonalized, it is possible to obtain a closed-form solution (within a given truncation scheme) to the Heisenberg equation of motion for any operator O expressed in the diagonal basis.The operator must first be transformed according to the flow equation dO/dl = [η(l), O(l)], where η(l) collectively denotes both the scrambling and Wegner generators.This transformed operator also contains valuable information about the locality of the unitary transform, and can be used to extract both a localization length and as a measure of the "complexity" of the diagonalization procedure, which can be linked to the existence of Lieb-Robinson bounds in flow time [28].Specifically, the transformed creation operator takes the form k cq , and higher-order terms are neglected.A measure of the complexity of the transformed operators is given by the fraction of non-zero terms which appear in this operator expansion.In practice, we choose a cutoff value ϵ = 10 −6 below which we consider terms to be zero.The complexity in this sense is defined as where (A ∪ B) represents the set of all coefficients A i and B ijk in the operator expansion of c † i .We also define χ(ϵ) = |{x ∈ (A ∪ B)|x 2 > ϵ 2 }|.Results are shown in Fig. 3, demonstrating a qualitative difference between one and two spatial dimensions.In one dimension, we find a phase where χ(ϵ) tends to a constant, and χ(ϵ) ∝ (1/L) 3 for large system sizes, indicating a 'low complexity' situation at strong disorder, as well as a higher complexity phase at small values of d where χ(ϵ) increases rapidly with system size, suggestive of thermalisation.In two dimensions, we find that χ(ϵ) always increases, although in the case of the QP potential at large values of d, it increases sufficiently slowly that the normalised complexity χ(ϵ) still vanishes.By contrast, for small values of d in two dimensions, the complexity χ(ϵ) remains much larger than zero for all system sizes studied here.This suggests a slow crossover from a high complexity phase -consistent with the expectation of thermalisation at small values of d -to a low complexity phase with anomalous thermalisation properties.
Previous work which have used CUT methods to compute non-equilibrium dynamics [25,26,29] employed a computationally costly inversion of the unitary transform in order to obtain time-evolved operators in the original basis.Here, we circumvent this limitation and directly obtain the infinite-temperature autocorrelation function, a highly nontrivial and challenging quantity to compute.The thermal expectation value of any arbitrary operator O is given by ⟨O⟩ = Tr[exp(−βH)O]/Tr[exp(−βH)], where β = 1/T is the inverse temperature (in units of k B = 1).In the limit of T → ∞, the expectation value becomes a uniform average over eigenstates, which in the diagonal basis are trivial product states.We approximate this average for large systems by randomly sampling N s ∈ [50,256] half-filled eigenstates, with an associated standard error that scales as 1/ √ N s .Specifically, we compute the dynamical autocorrelation function To minimize boundary effects, we choose i to be in the center of the system.We benchmark the performance of this approximation by comparison with exact diagonalization, making use of dynamical quantum typicality [30][31][32][33][34][35][36][37] to compute the infinite temperature correlation function.This is a highly demanding quantity that can be extremely challenging to compute with other methods, but can be obtained very efficiently with the flow equation approach.Results for system size L = 16 (L 2 = 4 × 4 for two dimensions) are shown in Fig. 4, demonstrating excellent agreement with exact diagonalization in both one and two dimensions, even to timescales as long as tJ = 10 5 .This is beyond the naive expectation that the accuracy should break down beyond timescales tJ ∼ 1/∆ 2 0 , corresponding to the typical inverse magnitude of the terms cut off by the truncation.We also show results for larger system sizes, demonstrating that they remain smooth and under control.Fig. 5 shows results for system sizes up to L = 100 (L 2 = 10 × 10 in two dimensions), along with a linear fit indicating the L → ∞ behavior.Strikingly, very little dependence on system size is observed in one dimension, although in two dimensions there is a slow trend towards decreasing values of C(t) as the system size increases, except for strong quasiperiodic potentials.This is consistent with the expectation that many-body localization may be ultimately unstable in two dimensions, although these results suggest one must go to very large system sizes and long times to see signs of this potential instability.For weak random disorder in two dimensions, the linear fit for large values of L breaks down, suggesting our results are likely to overestimate the long-time value of C(t) as L → ∞.Additionally, at any finite order of truncation there may still exist higher-order processes which could contribute to thermalisation on very long timescales.Nonetheless, for a given truncation scheme we can make precise statements about the validity of this technique.
To do so, we develop an incremental bound on the error in the unitary transform.If at each flow time step we discard all newly-generated terms above fourth-order, we obtain where ] represents the terms of the flow which are kept and A(l) = [η (4) (l), H (4) (l)] + T represents the truncation error, where the higher-order terms T are assumed to be negligible in what follows.The norm of the truncation error A(l) at each infinitesimal time step is upper bounded by using the sub-multiplicativity of the ∥.∥ 2 norm [38].The total truncation error in flow time can be written as an integral of Eq. ( 7) such that the mean value of this squared energy difference is 2d 2 /3, and that the largest parts of the interaction tensor remain proportional to the initial interaction strength (as new terms induced by the flow should always be smaller than the initial interactions), the error can be approximated as In the case of weak (quasi)disorder, the disorder bandwidth d is replaced by the effective bandwidth d ≥ d induced by the scrambling transform [39].A numerical analog can be computed by replacing the Hilbert-Schmidt (or Frobenius) norms in Eq. ( 7) with tensor Frobenius norms; the typical truncation error at each flow time step is well below one percent [39].
The above analysis indicates that energy differences below O(J∆ 2 0 /d 2 ) cannot be reliably resolved, implying that the method will break down on timescales on the order t ∝ d 2 /(J∆ 2 0 ) when oscillations at corresponding frequencies ω ∼ O(J∆ 2 0 /d 2 ) become relevant for the dynamics.The accuracy of the method can be systematically improved by incorporating additional higher-order terms into the truncated Hamiltonian, allowing accurate simulations of quantum dynamics to even longer times (proportional to 1/∆ 3 0 at the next order of approximation) with a computational cost that remains polynomial in the system size.Future developments in massively parallel implementations of the tensor flow equation method [20] used in this work, as well as advances in computer hardware, will facilitate extension of this method to larger system sizes, longer timescales, stronger interactions, and additional physical systems (including both driven [27,40] and dissipative [41] systems, which have been previously studied with CUT-based techniques).Scrambling transforms may be of interest in a variety of other contexts, as they are essentially a way of transforming a highly entangled system into a simpler representation that is easier to simulate.
We end the discussion by briefly comparing the findings with tensor network methods.Standard tensor network methods are challenged in time evolution by the exponentially growing bond dimension that is required to accommodate the states faithfully in time, for the linear growth of entanglement.Some steps have already been taken to allow tensor networks to access longer times [42][43][44][45], e.g., by means of folding techniques [44] or adaptive mode transformations [43]: The latter overcomes the prejudice that quantum states have to be represented in real space: One can co-rotate the frame of mode transformations, so that only the entanglement between these effective modes need to be accommodated.The ideas introduced here show that one can go further than that, by overcoming the prejudice that the fermionic mode transformation needs to be linear: accepting small truncation errors in the procedure, one can even fully diagonalise the Hamiltonian to a good approximation.There are good reasons to believe that this is a favourable mindset to reach long simulation times.First connections between flow equation and tensor network methods have been made [46], anticipating the time-dependent variational principle based on a differential geometric picture [47].It is conceivable that the ideas introduced here can be further merged with tensor network techniques, as one could think of final Hamiltonians that are not treated as fully diagonal ones.There is also the intriguing possibility of combining scrambling transforms and other CUT-based techniques with tensor network approaches, such as entanglement-based continuous unitary transforms [48], which may allow tensor network methods to break through the entanglement barrier and access much longer timescales than are currently available.
In this work, we have introduced a flow-based method equipped with scrambling transforms that allow for simulating interacting fermionic quantum many-body systems to good accuracy for intermediate and long times.Adding modern machine learning techniques into the mix could also allow for the development of efficient data-driven scrambling transforms tailored for specific problems.Such a classical development can also be seen as a challenge to dynamical quantum simulators [2,4] which aim to probe non-equilibrium properties of quantum matter beyond the reach of classical computers.These are exciting avenues for future progress.

ACKNOWLEDGEMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No.101031489 (Ergodicity Breaking in Quantum Matter) and under the Quantum Flagship (PASQuanS2, Millenion), as well as support from the NVIDIA Corporation through the Academic Hardware Grant Program, the BMBF (FermiQP, MuniQC-Atoms), and the DFG (CRC 183).S. J. T. thanks J. Catton for providing the artwork used in Fig. 1.Data and code are available at Refs. [40,49].

METHODS
This work made use of the PyFlow library [40], developed and maintained by S. J. T. , based on Ref. [20].

A. Computing and integrating the flow equation
All commutators computed in this work follow the scheme of Ref. [20], where the representation of the Hamiltonian in terms of a quadratic component (stored in memory as a matrix) and quartic component (stored as a tensor) allow the commutators to be recast in terms of matrix/tensor contractions, which are highly-optimised linear algebra operations that can be performed efficiently on modern computing hardware.A complete description is contained in the Supplemental Information [39].We use vacuum normal-ordering, such that higher-order terms in the running Hamiltonian have no feedback onto lower-order terms.The incorporation of additional non-perturbative corrections due to different choices of normal-ordering has previously been done in the time-independent scenario [50], but in the non-equilibrium setting is left for future work.This would require specifying a particular reference state with respect to which the corrections are computed, and depending on the structure of this state, the resulting corrections may be hard to vectorize, complicating their implementation on graphics processing units (GPUs).Calculations for all system sizes with more than a total of 16 lattice sites were performed on GPUs (specifically, NVIDIA RTX A5000 GPUs with 24Gb RAM and NVIDIA RTX 2080Ti GPUs with 12Gb RAM) using single precision arithmetic.
The flow equation dHl/dl = [η(l), H(l)] is solved using a mixed 4th/5th order Runge-Kutta integration method as implemented in the JAX library [51], making use of an adaptive stepsize algorithm for high accuracy.The maximum integration time used was l max = 1000, and the integration is stopped before then if the Hamiltonian is diagonalized to the target accuracy, which we choose to be when max[|V (2) |] < 10 −6 and max[|V (4) |] < 10 −3 .Results for longer integration times showed no significant increase in accuracy, despite incurring a significantly higher computational cost.This is because the running Hamiltonian H(l) approaches full diagonalisation only asymptotically at large values of l, so the use of larger values of l max leads to diminishing returns.

B. Computing the dynamics
The transformed number operator is reconstructed from the transformed creation/annihilation operators for large fictitious time with the creation operator given by transformed creation/annihilation operators and the annihilation operator obtained by taking its Hermitian conjugate Multiplying these together allows us to reconstruct the number operator including terms up to sixth-order in the fermionic creation/annihilation operators for the diagonal basis, c † i and ci .The number operator can then be time-evolved in the diagonal basis according to the Heisenberg equation of motion, neglecting newly-generated higher-order terms, resulting in a closed-form solution.This step is performed on CPUs rather than GPUs due to memory limitations, and is a prime candidate for future efficiency improvements.At long times, near-degenerate single-particle eigenvalues can still lead to divergent terms in this solution (consistent with the expectation that the simulation of a BQPhard problem will eventually run into accuracy issues on a classical computer), however, these terms are strongly suppressed, arising only very rarely and at very long times.To avoid these rare scenarios dominating the averaged data, in the two dimensional dynamical data we exclude disorder realizations where the maximum value of |C(t)| > 1.1.(Alternatively, we could have used the typical rather than mean value of C(t).)See the Supplementary Information for full details of the calculation and where divergent terms arise from.In one dimension, this procedure is not necessary as the divergent terms are rare enough to have essentially no effect, as can be seen from the data in Fig. 4. It is possible to subtract the divergent terms when they are encountered, however, we do not employ this procedure here.The long-time average can be obtained directly by setting all off-diagonal terms in the transformed number operator to zero (as when time-evolved, they acquire oscillating phases which average to zero).For systems with greater than 36 lattice sites in total, we neglect the sixth order contributions and keep only the quadratic and quartic terms when computing the dynamics.For the systems considered here, the sixth-order terms have a negligible effect, which can be seen from the qualitative agreement between small and larger system sizes.

C. Rescaling the correlation function
As the norm of the number operator n i is not precisely conserved by the unitary transform, we rescale the correlation function for each disorder realization according to the ansatz C(t) →c 1 (C(t) − c 2 ), where c 1 and c 2 are determined by minimizing the error with respect to the short-time dynamics of the non-interacting system (as many-body interactions are essentially irrelevant at very short times).This is computationally efficient, as we get the exact dynamics of the non-interacting system essentially for free in this formalism by just retaining the quadratic components of the Hamiltonian and relevant observables.The rescaling employed in this work is justified a posteriori by the clear agreement between the rescaled C(t) and the exact result, computed for system sizes small enough for the comparison to be practical.For small enough systems, an alternative would be to construct the operator as a matrix in the full Hilbert space and renormalize it by hand, however, this is not practical for systems as large as those considered here.We emphasize that the norm is preserved to high accuracy for sufficiently strong disorder, and the effects of this rescaling are most important for the weakly disordered systems.This is independent of any error introduced in the eigenvalues, and reflects the difficulty in simultaneously preserving the unitary evolution of both the Hamiltonian and the number operator within the same truncation scheme.The norm of the operator could in principle be exactly preserved by constructing the unitary transforms subject to additional constraints [18], however, in practice this is challenging to implement.This underscores the need for further work in developing more flexible genera-tors for the types of continuous unitary transform developed here, perhaps in concert with machine learning approaches to design data-driven generators tailored for specific problems subject to specific hard-to-satisfy constraints.flow equation into an efficient numerical procedure.Any Hamiltonian may be written in a generic form as where the superscripts indicate how many operators appear in each term, e.g., H (2) (l) contains all quadratic terms, H (4) (l) all quartic terms and so on.This can be stored in memory as a matrix (H (2) ) and a fourth-order tensor (H (4) ), plus higher-order terms as required.A similar decomposition can be performed for the generator η(l):=η (2) (l) + η (4) (l) + ... .(14) This means that the flow equation can be computed as Each of these commutators can be rewritten as the sum of all one-point contractions between the arrays η (n) (l) and H (m) (l), as shown in Fig. 7a) and Fig. 7b), with the result being an array of order (n + m − 2) [20].Note that if any arrays are present with an order greater than 2, the above commutator cannot be written as a closed expression and some form of truncation must be employed.The source term in the flow (Eq.( 15)) which generates sixth order terms is given by the elements of which we can estimate to be of order O(J 0 ∆ 2 0 ) or smaller.A more detailed error analysis is given in the main text.For weak interactions, even the induced sixth-order terms will be extremely small, and the eighth-order and higher terms will be essentially negligible.The exception to that is if any elements of H (4) (l) become of order one during the flow, at which point the higher order terms may be required to retain accuracy.The chances of this happening are heavily suppressed by the scrambling transform.We shall see later that higher-order terms in this framework enter with increasing powers of the interaction strength, and so by working at weak interactions we can ensure that the higher-order terms are essentially negligible.
Two-point and higher contractions can also be taken into account [50], however, we shall not consider these in the present manuscript as they are only non-zero if a state other than the vacuum is chosen for the normal-ordering procedure.This can be important for the incorporation of sophisticated non-perturbative corrections [50] or for finite-temperature contributions [18], but in the present case we will take the simpler approach of using vacuum normal-ordering.
The computational cost of this technique is set by the number of lattice sites, and by the highest-order of terms included in Eq. (1).In this work, we shall consider running Hamiltonians with up to fourth order (O(L 4 )) and sixth order (O(L 6 )) terms included.Crucially, the overall computational cost does not depend on the geometry or dimensionality of the lattice sites: for example, a one-dimensional system of size L = 64 can be diagonalized with essentially the same computational cost and accuracy as a two-dimensional system of size L = 8 × 8, or a three dimensional system of size L = 4 × 4 × 4.

Error analysis
The relative error in the i-th eigenvalue E i is defined as where the superscripts F E and ED refer to the eigenvalues obtained by flow equations and exact diagonalisation respectively.We compute the median error within each disorder realization (as the error distribution has significant skewness, and a small number of outliers dominate the mean error), and then average this over disorder realizations to obtain the results shown in Fig. 8.We find that the error remains small and under control for all system sizes and disorder strengths considered here.As expected, the error is larger in two dimensions than in one dimension, however, the error shrinks with increasing system size and remains within acceptable limits.We can also define a second error metric, this time a self-consistent measure of the truncation error accrued over the course of the entire flow, based on a modified version < l a t e x i t s h a 1 _ b a s e 6 4 = " A schematic of how the commutator between an order two and four array may be computed in terms of one-point contractions.The most important thing is that the indices are put in a consistent order -using appropriate (anti)commutation relations -before the four resulting arrays are added together.  of Eq. ( 8) from the main text, where the norms ∥.∥ F are now tensor Frobenius norms (i.e., square root of the sum of the squares of the entries of each tensor), and the integral is normalised by the total flow time, giving a measure of the average truncation error per timestep over the course of the entire diagonalization procedure.We approximate this integral numerically using the trapezoidal rule.Results are shown in Fig. 9, demonstrating that the errors due to the truncation of the running Hamiltonian remain small and under control throughout the procedure.
Overall, there are two sources of errors: This is on the one hand the dominant truncation error in the non-linear transformations along the fictitious time of the flow, bounded from above in Frobenius norm above.On the other hand, there is the error arising from pinching, when replacing matrices and tensors that are to a very good approximation diagonal for the final fictitious time l max by exactly diagonal matrices and tensors as b) The same quantity in two dimensions.Note that the error bars are much larger for random disorder than for the quasi-periodic potential, reflecting the greater likelihood of rare disorder-free regions (often known as 'resonant regions' in the MBL literature) which act as effectively delocalised, and would require higher-order terms in the truncation to fully capture.These terms are typically absent in deterministic quasi-periodic potentials.

Computing the flow of operators
Once we have diagonalized the Hamiltonian, we will wish to compute various observables in the diagonal basis, and so we need a way to transform arbitrary operators into the same basis as the Hamiltonian.This can be done straightforwardly, similarly to how the Hamiltonian itself is transformed.Rather than directly transforming fermionic number operators, as in previous works, here we will work directly with fermionic creation/annihilation operators.These are cheaper to compute by a factor of L, and once we know the form of either the creation or annihilation operator, the other follows straightforwardly by taking the Hermitian conjugate.
The flow of the creation operator is given by which can be computed at the same time as the flow of the Hamiltonian, meaning that η(l) does not need to be stored and can be discarded after each time step.Under the action of the flow, the creation operator takes the form which can be stored in memory as a unit order array (vector) plus an order three array (tensor), in contrast to the number operator which must be stored as a larger order two array plus an order four array.To reconstruct the number operator at any point during the flow, we simply multiply Eq. ( 22) with its Hermitian conjugate using the definition Note that in addition to being cheaper by a factor of L as compared with previous work which directly transformed the number operator up to quartic order [20,50,52], this also allows us immediate access to the sixth-order contributions to the number operator, which were not included in previous work.This approach is therefore faster, more efficient, and more accurate.In practice, we compute the creation operator to order (O − 1), where O ∈ {4, 6} is the order of the effective Hamiltonian.For O = 6, this allows us to reconstruct the number operator up to terms including ten fermionic operators, in principle allowing us to take into account extremely weak effects that would only appear at high orders in perturbation theory.In practice, we do not go beyond sixth order, as otherwise the calculation becomes too memory-intensive for the present hardware, and for the largest system sizes considered in this work we do not go beyond fourth order.This also allows us to transform (product) states from the microscopic basis to the diagonal basis, and vice-versa.We can write a generic product state vector as |ψ⟩ = Π i c † i |0⟩ and transform each of the creation operators individually, multiplying the resulting expressions together at the end to obtain the transformed state.If the series of unitary transforms used to diagonalise the Hamiltonian has been stored, one can take any basis state (which are simply product states of zeros and ones in the diagonal basis) and reverse the transform in order to construct any desired eigenstate in the original microscopic basis.The resulting expression is likely to be quite complex, however, and it is unclear whether this procedure can be made of practical use.

Time evolution in the diagonal basis
In the basis where the Hamiltonian is diagonal, operators can be time-evolved using the Heisenberg equations of motion This can again be handled on a term-by-term basis.For an operator O containing an odd number of fermionic operators, one obtains the expression While for an operator O containing an even number of fermionic operators, one obtains Once the operator has been time-evolved, we can compute its expectation value with respect to any basis state, or mixture of basis states.If the sequence of unitary transforms used to diagonalise the Hamiltonian has been stored, the transform can also be reversed in order to obtain the time-evolved operator back in the original basis.This is useful if one wishes to study operator spreading in real time, or compute the non-equilibrium dynamics starting from an initial state that is particularly simple in the original basis, e.g., some form of unentangled state like a Néel state.
To be specific, we will be particularly interested in computing time-evolved number operators.The number operators take the form Higher order terms can be included, however, for the systems considered in this work they do not lead to any measurable difference in the results, and so we do not include them here.In all of the following, we drop the tilde notation for clarity, however, it is to be understood that from this point onwards, all operators are defined in the diagonal (tilde) basis.The number operator can now be time evolved according to the Heisenberg equation dn i (t)/dt = i[ H, n i (t)] and the contributions can be evaluated step-by-step.We could alternatively time-evolve the creation/annihilation operators, allowing us to avoid having to store the sixth-order term in computer memory, at the cost of having to perform a more complex series of summations when evaluating expectation values.The first contribution is given by the quadratic components of both the Hamiltonian and the number operator This is the only contribution to the quadratic part of the number operator, and it can be integrated to obtain with The first term is constant, as the terms that govern time evolution commute.There are two contributions to the quartic term.The first involves the (diagonal) quartic Hamiltonian [ H(4) , n The second involves the quadratic part of the Hamiltonian and the quartic part of the number operator as which together with the previous contribution implies  for all times t ≥ 0. One can proceed similarly for the quartic term in the fermionic operators.
In the very rare case of exact single-particle degeneracies, the [∆ ik (β ij (t)−β ij (0))/( hi − hj )] and [∆ ik (β kq (t)−β kq (0))/( hk − hq )] terms are replaced by ∆ ik β ij (0)t and ∆ ik β kq (0)t respectively.All higher-order terms can be evaluated similarly.As in the calculation of the commutators used in the flow equations, these can be efficiently computed using the graphical notation.In this work, we compute only the leading order phase shift applied to the sixth order terms coming from [ H(2) , n i ], which leads to Sixth-order terms are not included for system sizes L > 36, as the memory cost becomes prohibitive, but the results for larger system sizes are entirely consistent with the trends observed for smaller systems, demonstrating that these terms have only a minor effect in all observed cases.With this done, we have computed the time-evolved number operator at some arbitrary time t ≥ 0, without reference to any state, and without requiring step-by-step numerical evolution.Now we can move on and use this to compute observables.
Computing the infinite-temperature correlation function In previous works which have used flow equation techniques to compute the non-equilibrium dynamics of observables, the procedure has typically involved transforming the desired operator into the diagonal basis (where time evolution becomes easier), performing the time evolution, and then reversing the unitary transform to write the time-evolved operator back in the original basis, where appropriate expectation values could then be computed.This process was cumbersome, as it involved two unitary transforms (each with associated small but non-negligible numerical error), required storing the entire series of infinitesimal unitary transforms in memory, and the reverse transform had to be performed separately for each individual timestep.Altogether, this led to time evolution using flow equation techniques to be extremely computationally demanding, limiting their use.
Here, we demonstrate that infinite-temperature expectation values may be computed using a much more efficient approach.The thermal expectation value of any arbitrary operator O reflecting the Gibbs state is given by ⟨O⟩ = Tr[exp(−βH)O]/Tr[exp(−βH)], (39) where β = 1/T is the inverse temperature (in units of k B = 1).In the infinite temperature limit, β → 0, both exponentials reduce to the identity.In particular, the denominator simply becomes the trace of the identity, which gives the Hilbert space dimension, and the infinite-temperature expectation value becomes where the {|ψ i ⟩} are the basis state vectors.This is simply the mean value of the expectation value of the operator O taken over the entire Hilbert space.
While there are elegant ways to approximate this using dynamical typicality (where the sum over states is replaced by the expectation value in a single randomly chosen pure state), here we will take a simpler approach.We will replace the Hilbert space size D H with a smaller integer N s < D H , and approximate the mean value by For a sufficiently large value of N s , this converges to the true value.In the following, we will choose N s = min(D H , 512) and restrict ourselves to the zero magnetization (half-filled) sector.We will work in the basis in which the Hamiltonian is diagonal, and the basis vectors are simply given by all possible binary strings |0, 1, 0, . ..⟩.This allows us to avoid having to explicitly compute all eigenstates in the initial microscopic basis, which would be exponentially costly in the system size.This can be done by inserting the identity U U † = 1 like so reducing the problem to computing the transformed operator Õ = U † OU , and then computing its expectation value across N s product states.
In this work, we are interested in computing the infinite-temperature autocorrelation function, defined as where the factor of 4 is chosen for convenience such that C i (0) = 1, and where n i (t) is obtained as previously described.We shall choose site i to be in the middle of the system in all the following so that we are far from the boundaries and finite-size effects should be minimized.We will benchmark the performance of this approximation by comparison with exact diagonalization, making use of dynamical quantum typicality [30][31][32][33][34][35][36][37] to efficiently compute the infinite temperature correlation function with an error exponentially small in the Hilbert space dimension D H .In this case, the trace over basis states will be replaced by the expectation value taken with respect to a single state vector where C > 0 is a normalization constant, a k and b k are random real variables chosen from Gaussian distributions of zero mean and variance 1/2.The state is chosen randomly for each disorder realization.This is essentially a way of approximating the trace in Eq. ( 40) by uniformly randomly sampling contributions from all basis states, ⟨(: n (2) (t) : + : n < l a t e x i t s h a 1 _ b a s e 6 4 = " j b L x X J a k I 5 C A S K n e v l L w w 6 o o + i A = " > A A A B 7 X i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R a h X s q u F O 1 F K H j p s Y L 9 g H Y p 2 T T b x m a T J c k K Z e l / 8 O J B E a / + H 2 / + G 9 N 2 D 9 r 6 Y O D x 3 g w z 8 4 K Y M 2 1 c 9 9 v J b W x u b e / k d w t 7 + w e H R 8 X j k 7 a W i S K 0 R S S X q h t g TT k T t G W Y 4 b Q b K 4 q j g N N O M L m b + 5 0 n q j S T 4 s F M Y + p H e C R Y y A g 2 V m o 3 y v z W v R w U S 2 7 F X Q C t E y 8 j J c j Q H B S / + k N J k o g K Q z j W u u e 5 s f F T r A w j n M 4 K / U T T G J M J H t G e p Q J H V P v p 4 t o Z u r D K E I V S 2 R I G L d T f E y m O t J 5 G g e 2 M s B n r V W 8 u / u f 1 E h P W / J S J O D F U k O W i M O H I S D R / H Q 2 Z o s T w q S W Y K G Z v R W S MF S b G B l S w I X i r L 6 + T 9 l X F u 6 5 U 7 6 u l e i 2 L I w 9 n c A 5 l 8 O A G 6 t C A J r S A w C M 8 w y u 8 O d J 5 c d 6 d j 2 V r z s l m T u E P n M 8 f H b i O J w = = < / l a t e x i t > H(l = 0) < l a t e x i t s h a 1 _ b a s e 6 4 = " w EA t y Z M a M r q U k B 7 d W R Z D g J / I l H U = " > A A A B + H i c b V B N S 8 N A E N 3 U r 1 o / G v X o Z b E I 9 V I S K d p j w U u P F e w H N K F s t p t 2 6 W Y T d i d C L P 0 l X j w o 4 t W f 4 s 1 / 4 7 b N Q V s f D D z e m 2 F m X p A I r s F x v q 3 C 1 v b O 7 l 5 x v 3 R w e H R c t k 9 O u z p O F W U d G o t Y 9 Q O i m e C S d Y C D Y P 1 E M R I F g v W C 6 d 3 C 7 z 0 y p X k s H y B L m B + R s e Q h p w S M N L T L r a r A H s T Y 4 z K E 7 G p o V 5 y a s w T e J G 5 O K i h H e 2 h / e a O Y p h G T Q A X R e u A 6 C f g z o o B T w e Y l L 9 U s I X R K x m x g q C Q R 0 / 5 s e f g c X x p l h M N Y m Z K A l + r v i R m J t M 6 i w H R G B C Z 6 3 V u I / 3 m D F M K G P + M y S Y F J u l o U p g K b R x c p 4 B F X j I L I D C F U c X M r p h O i C A W T V c m E 4 K 6 / v E m 6 1 z X 3 p l a / r 1 e a j T y O I j p H F 6 i K X H S L m q i F 2 q i D K E r R M 3 p F b 9 a T 9 W K 9 W x + r 1 o K V z 5 y h P 7 A + f w B x 0 p J K < / l a t e x i t > H(l ! 1) < l a t e x i t s h a 1 _ b a s e 6 4 = " j b L x X J a k I 5 C A S K n e v l L w w 6 o o + i A = " > A A A B 7 X i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R a h X s q u F O 1 F K H j p s Y L 9 g H Y p 2 T T b x m a T J c k K Z e l / 8 O J B E a / + H 2 / + G 9 N 2 D 9 r 6 Y O D x 3 g w z 8 4 K Y M 2 1 c 9 9 v J b W x u b e / k d w t 7 + w e H R 8 X j k 7 a W i S K 0 R S S X q h t g T T k T t G W Y 4 b Q b K 4 q j g N N O M L m b + 5 0 n q j S T 4 s F M Y + p H e C R Y y A g 2 V m o 3 y v z W v R w U S 2 7 F X Q C t E y 8 j J c j Q H B S / + k N J k o g K Q z j W u u e 5 s f F T r A w j n M 4 K / U T T G J M J H t G e p Q J H V P v p 4 t o Z u r D K E I V S 2 R I G L d T f E y m O t J 5 G g e 2 M s B n r V W 8 u / u f 1 E h P W / J S J O D F U k O W i M O H I S D R / H Q 2 Z o s T w q S W Y K G Z v R W S MF S b G B l S w I X i r L 6 + T 9 l X F u 6 5 U 7 6 u l e i 2 L I w 9 n c A 5 l 8 O A G 6 t C A J r S A w C M 8 w y u 8 O d J 5 c d 6 d j 2 V r z s l m T u E P n M 8 f H b i O J w = = < / l a t e x i t > H(l = 0) < l a t e x i t s h a 1 _ b a s e 6 4 = " w EA t y Z M a M r q U k B 7 d W R Z D g J / I l H U = " > A A A B + H i c b V B N S 8 N A E N 3 U r 1 o / G v X o Z b E I 9 V I S K d p j w U u P F e w H N K F s t p t 2 6 W Y T d i d C L P 0 l X j w o 4 t W f 4 s 1 / 4 7 b N Q V s f D D z e m 2 F m X p A I r s F x v q 3 C 1 v b O 7 l 5 x v 3 R w e H R c t k 9 O u z p O F W U d G o t Y 9 Q O i m e C S d Y C D Y P 1 E M R I F g v W C 6 d 3 C 7 z 0 y p X k s H y B L m B + R s e Q h p w S M N L T L r a r A H s T Y 4 z K E 7 G p o V 5 y a s w T e J G 5 O K i h H e 2 h / e a O Y p h G T Q A X R e u A 6 C f g z o o B T w e Y l L 9 U s I X R K x m x g q C Q R 0 / 5 s e f g c X x p l h M N Y m Z K A l + r v i R m J t M 6 i w H R G B C Z 6 3 V u I / 3 m D F M K G P + M y S Y F J u l o U p g K b R x c p 4 B F X j I L I D C F U c X M r p h O i C A W T V c m E 4 K 6 / v E m 6 1 z X 3 p l a / r 1 e a j T y O I j p H F 6 i K X H S L m q i F 2 q i D K E r R M 3 p F b 9 a T 9 W K 9 W x + r 1 o K V z 5 y h P 7 A + f w B x 0 p J K < / l a t e x i t > H(l !1)< l a t e x i t s h a 1 _ b a s e 6 4 = " x Y P u v q A R h U m h K z R T 5 l B v 4 t K d 4 k w = " > A A A B + H i c b V B N S 8 N A E J 3 U r 1 o / G v X o Z b E I n k o i R X s s e O m x g m m F N p b N Z t M u 3 W z C 7 k a o o b / E i w d F v P p T v P l v 3 L Y 5 a O u D g c d 7 M 8 z M C 1 L O l H a c b 6 u 0 s b m 1 v V P e r e z t H x x W 7 a P j r k o y S a h H E p 7 I + w A r y p m g n m a a 0 / t U U h w H n P a C y c 3 c 7 z 1 z h P P i v D s f y 9 a c k 8 2 c w h 8 4 n z + H 9 4 1 H < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " j R e N x e + C R 4 I L A v q P e Y 9 e p K 5 Q H 1 1 E B N 1 E I E P a J n 9 I r e n C f n x X l 3 P p a t O S e b O U F / 4 H z + A P / h k 9 A = < / l a t e x i t > Final Hamiltonian < l a t e x i t s h a 1 _ b a s e 6 4 = " b B 4 +

Figure 1 .
Figure 1.A cartoon illustration of the scrambling process.a) The conventional CUT process, using a single unitary transform U to diagonalise a Hamiltonian H, smoothly transforming from the initial basis (l = 0) to the diagonal basis (l → ∞).b) The scrambling transform S first induces an 'effective disorder' even in completely clean systems, which allows established CUT techniques to then take over and efficiently diagonalise the 'scrambled' Hamiltonian S † HS with a second unitary transform U .
arXiv:2308.13005v1 [quant-ph] 24 Aug 2023 < l a t e x i t s h a 1 _ b a s e 6 4 = " y Q DJ R M v 3 V m U w 5 H E q / v L v u F W N y w s = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g a d m V o B 4 D X j z G R x 6 Q L G F 2 0 p s M m Z 1 d Z m a F E P I J X j w o 4 t U v 8 u b f O E n 2 o I k F D U V V N 9 1 d Y S q 4 N p 7 3 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R U y e Z Y t h g i U h U O 6 Q a B Z f Y M N w I b K c K a R w K b I W j m 5 n f e k K l e S I f z T j F I K Y D y S P O q L H S g + u 6 v X L F c 7 0 5 y C r x c 1 K B H P V e + a v b T 1 g W o z R M U K 0 7 v p e a Y E K V 4 U z g t N T N N K a U j e g A O 5 Z K G q M O J v N T p + T M K n 0 S J c q W N G S u / p 6 Y 0 F j r c R z a z p i a o V 7 2 Z u J / X i c z 0 X U w 4 T L N D E q 2 W B R l g p i E z P 4 m f a 6 Q G T G 2 h D L F 7 a 2 E D a m i z N h 0 S j Y E f / n l V d K 8 c P 1 L t 3 p X r d T u 8 z i K c A K n c A 4 + X E E N b q E O D W A w g G d4 h T d H O C / O u / O x a C 0 4 + c w x / I H z + Q N T o o 0 7 < / l a t e x i t > ...

Figure 2 .
Figure 2.An illustration of how a two-dimensional lattice can be mapped onto a one-dimensional chain with correlated long-range hopping, which can be easily handled with CUT-based techniques.

Figure 3 .
Figure 3.The complexity χ of the transformed creation operator c † i in the middle of the system.Error bars show the standard deviation over disorder realizations.Panels (a) and (b) show the results in one dimension (L = 8, 10, 12, 16, 24, 36, 48, 64, 100) for random and quasi-periodic potentials respectively, while panels (c) and (d) show the same in two dimensions (L 2 = 9, 16, 25, 36, 49, 64, 100).Dashed black lines close to the origin in panels (a) and (b) are fits with the form χ ∝ 1/L 3 , valid for large systems and strong disorder.Insets show the un-normalised complexity (i.e., the numerator of Eq. (4)), which tends to a constant in strongly-disordered one-dimensional chains, but grows in two dimensions even for strong disorder.

Figure 4 .
Figure 4. Infinite-temperature correlation functions shown for a variety of disorder types, strengths, dimensionalities and system sizes.The grey dot-dashed vertical line indicates the approximate timescale beyond which accuracy cannot be guaranteed, however, the results typically remain reasonable until much longer timescales.In all plots, the result from exact diagonalisation (ED) is shown as a pale dashed line.The equivalent result from the flow equation (FE) method is shown by circular markers, and the result obtained using flow equations for a much larger system is indicated by a dark solid line.Black dashed lines show the long-time average computed directly without explicit time evolution, valid at strong (quasi)disorder only.The results are averaged over Ns ∈ [64, 256] disorder realizations, depending on system size.Error bars indicate the variance over disorder realizations, for clarity shown only for the flow equation results for L = 16 (L 2 = 4 × 4 in two dimensions).In both D = 1 and D = 2, the quasi-periodic potential exhibits most much robust localization at large values of d/J, but by contrast also exhibits more complete thermalisation at low values of d/J due to the underlying single-particle phase transition at d/J = 2.

Figure 5 .
Figure 5. Finite-size scaling results for the long-time average of the infinite-temperature correlation function C(t), averaged over a time window of Jt ∈ [50, 10 3 ], and then averaged over Ns ∈ [32, 1024] disorder realizations depending on system size.(a-d) Results for various system sizes, disorder types and strengths, and dimensionalities.Error bars indicate the variance over disorder realizations.The black squares indicate the results obtained from exact diagonalisation, shown for small system sizes only.System sizes are L = 8, 10, 12, 16, 24, 36, 48, 64 and 100 in one dimension, and L 2 = 9, 16, 25, 36, 49, 64 and 100 in two dimensions.(e) The result of linearly extrapolating to L → ∞.Error bars indicate the uncertainty in the fits [shown as dashed lines in panels (a-d)] to extract the scaling behavior; lines are smoothed guides to the eye.

Figure 6
Figure 6.a) The induced disorder bandwidth following the application of the scrambling step, with hopping J = 1 and different values of the cutoff ϵ and the microscopic (random) disorder strength d.System size if L = 36, and the results are averaged over 50 disorder realisations.Error bars indicate the standard deviation over disorder realisations.b) A comparison of the median relative error ε in the eigenvalues for the method with (blue) and without (orange) the scrambling transform, for an interacting system of size L = 10 with J = 0.5, interaction strength ∆0 = 1.0, maximum flow time lmax = 50 and averaged over 192 disorder realisations.Error bars indicate the median absolute deviation.

K b 8 6 j 8 ++Figure 7
Figure 7. a) A schematic of how the commutator between two matrices may be computed in terms of one-point contractions.Creation and annihilation operators are indicated by up and down arrows respectively, and connected indices are summed over.The indices of the resulting arrays are rearranged into a consistent order (taking into account the appropriate (anti)commutation relations), and then added together to get the final result.This is equivalent toi,j,k,q AijB kq [c † i cj, c † k cq] = i,j,k A ik B kq c † k cq + AijB ki cjc † k .b)A schematic of how the commutator between an order two and four array may be computed in terms of one-point contractions.The most important thing is that the indices are put in a consistent order -using appropriate (anti)commutation relations -before the four resulting arrays are added together.

Figure 8 .
Figure 8.The relative error in the eigenvalues computed with flow equations and with exact digaonalisation, for a variety of system sizes and disorder strengths, averaged over Ns ∈ [256, 1024] disorder realizations depending on system size.a) The relative error in one dimension.b) The relative error in two dimensions.

Figure 9 .
Figure 9.The average truncation error at each flow time step, as defined in the text, averaged over Ns ∈ [32, 1024] disorder realizations.Solid lines represent random disorder, and dashed lines show the results for the quasi-periodic potential.a) The truncation error in one dimension.b)The same quantity in two dimensions.Note that the error bars are much larger for random disorder than for the quasi-periodic potential, reflecting the greater likelihood of rare disorder-free regions (often known as 'resonant regions' in the MBL literature) which act as effectively delocalised, and would require higher-order terms in the truncation to fully capture.These terms are typically absent in deterministic quasi-periodic potentials.
ijkq (t) = exp[i( hi − hj + hk − hq )t]Γ ijkq (0) + δ ij ∆ ik (β kq (t) − β kq (0)) hk − hq + δ kq ∆ ik (β ij (t) − β ij (0) hi − hj .(35)Once the Hamiltonian is in diagonal form, one can also bound the error in the time evolution due to the imperfect diagonal elements.To upper bound the error in β If the numbers {h j } are the approximations of the exact {h ′ j } that are being used in the time evolution in the above Hamiltonian, then the error made can easily be upper bounded as follows.With jk (t) = exp(i(h j − h k )t)β j |h j − h ′ j |(36)being the largest element-wise error, one obtainsexp(i(h ′ j − h ′ k )t) − exp(i(h j − h k )t) β jk (0) = exp(i(h ′ j − h ′ k − h j + h k )t) − 1 exp(i(h j − h k )t)β k − h j + h k ) 2 t 2 2| exp(i(h j − h k )t)|β