Displacement Transformations as a Tool to Study Many-Body Localization

We obtain eigenstates of interacting disorder Hamiltonians using unitary displacement transformations that rotate the state of the system. The method generates excited states if the strength of these transformations is chosen to optimize the energy, while decreasing the energy variance. We apply the method to analyse the localization properties of one-dimensional spinless fermions with short range interactions, reaching relatively large system sizes. We quantify the degree of localization through the size and disorder dependence of the inverse participation ratio.


INTRODUCTION
In strongly disordered interacting systems, particles are localized, a phenomenon known as manybody localization (MBL) [1][2][3][4][5][6]. MBL can be understood through the existence of an extensive number of exponentially localized integrals of motion (IOM) [7][8][9][10]. A practical computational approach to obtain the IOM and diagonalize interacting fermionic Hamiltonians was developed with the help of displacement transformations [11,12], a type of unitary transformation (see Eq. 2). A similar method to diagonalize these Hamiltonians [13] had been previously introduced in a different context.
In a previous publication [14], two of us together with Louk Rademaker improved the numerical efficiency of the displacement transformations technique by focusing on a specific reference state and considering only those transformations that affect the energy of this state. In the present work, the spotlight is again on a reference state, but instead of diagonalising the Hamiltonian, we rotate the state of the system. As we will see, this further increases the efficiency of the method. One of the reasons for this increase is that storing a state is less expensive than storing a transformed Hamiltonian. We focus on the effect of the most influential displacement transformations to gain information on the localization properties of the state. The method proposed in this work is quite general, fairly easy to implement and it may be of interest in different fields, such as quantum chemistry, nuclear physics, Rydberg atoms and molecules, quantum spins, etc.
The flow equation method [15][16][17][18][19][20][21] is a continuous versions of our model, in which a set of differential equations is solved to implement the displacement transformations. Continuous transformations have the advantage of commuting one with each other for infinitesimal strengths, while discrete implementations take into account the term to be transformed at all orders in a single step. Besides, discrete transformations can handle better the storage of the states, which can be done in a sparse representation. This opens the possibility of reaching larger system sizes.
The remainder of the paper is set up as follows. In section 2 we present the Hamiltonian considered and briefly review the definition and main property of displacement transformations. In section 3 we deduce the equations to minimize the energy with displacement transformations for a general state. We then describe how to obtain excited states in section 4 and analyse their localization properties in 5. In section 6 we discuss possible applications and extract some conclusions.

MODEL HAMILTONIAN AND DISPLACEMENT TRANSFORMATIONS
Although our method applies to any fermionic quantum Hermitian Hamiltonian, we concentrate on the onedimensional tight-binding Anderson interacting Hamiltonian where c † i (c i ) is the fermionic creation (annihilation) operator of state i and n i c † i c i is the number operator. ϕ i are the site energies, distributed uniformly in the interval (−W/2, W/2), t is the hopping energy, 〈i, j〉 refers to nearest neighbours, and U i , j the interaction. The method works for any type of interaction, and here we will only consider a nearest neighbour interaction For an operator X containing the same number of creation and annihilation operators we can define the following unitary transformation This is called a displacement transformation of strength λ for the operator X. And the total number of creation and annihilation operators in X is called the order of the transformation. The Hamiltonian given by Eq. 1 can be diagonalized by products of displacement transformations [11], ending up of the form n is the transform of the density operator under the same transformation performed on the Hamiltonian. It is trivial to show that displacement transformations can be summed up, given that X 3 X, in the following compact way where 1 is the unit operator.

ENERGY MINIMIZATION
As we have mentioned, one can in principle diagonalize the Hamiltonian by the application of successive displacement transformations [11,12], but in this process high order terms in the Hamiltonian are generated and to store them for large sizes is a real challenge. An alternative approach is to focus on a given state |Φ〉 and perform only those transformations that affect the expectation value or the variance of the Hamiltonian with respect to this state. The expectation value of the energy after a transformation D X (λ) is [14].
where A ≡ X † − X and expectation values are computed with respect to state |Φ〉. In a previous work [14], we optimized the expectation value of the energy by applying the displacement transformations to the Hamiltonian, while keeping the state fixed.
In the present study we keep the Hamiltonian fixed and apply the transformations to the state. The final state will be When the state is kept fixed, it is normal to consider a single component state (possibly in the complex basis defined by the integrals of motion), and the only transformations that affect the energy expectation value satisfy 〈A 2 〉 − 1. In the general case, we can no longer assume this condition and Eq. 5 can be written as: where To obtain these quantities for a given X, we compute the auxiliary states |Φ 1 〉 A|Φ〉 and |Φ 2 〉 A 2 |Φ〉 and from them we evaluate the action of the Hamiltonian operator on them, H|Φ 1 〉 and H|Φ 2 〉. All quantities appearing in Eqs 8-11 are scalar products of the four states just mentioned. The task that requires most of the CPU time is the evaluation of H|Φ 1 〉 and H|Φ 2 〉. The function given by Eq. 7 has period π, and there is always an extreme of the energy between − π/2 and π/2. To get excited states, we choose the extreme corresponding to the minimum value of |λ|, since this is the most effective strategy to get as many excited states as possible, as we will see. Once we have the value of λ, the calculation of the transformed state is straight forward The energy is an extreme for all eigenvalues, but not necessarily a minimum, except for the ground state. So, to avoid stability problems for excited states, we reject transformations that increase the variance of the energy. This quantity is calculated from its definition.
The two expectation values on the right hand side of this equation are evaluated according to Eq. 7.
Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 767001 The previous procedure can be implemented for any oneelectron basis. The two most obvious choices are the site basis and the solution of the one-electron problem. In the site basis, the Hamiltonian takes a very simple form and it is very efficient to evaluate its action on a wavefunction. On the other hand, wavefunctions always have a very large number of components. The second possibility, the solution of the oneelectron problem, involves a more complicated form of the Hamiltonian and wavefunctions whose effective number of components varies drastically between the localized and the extended regime. We will call this basis the non-interacting basis and we will work with it for the reasons given later on.
To avoid confusion, we refer to states in the non interacting basis with greek letters. In this basis the Hamiltonian is

EXCITED STATES
Our main aim is to obtain excite states with as little bias as possible and study their localization properties. To get these states is not a trivial task and, for example, starting from a canonical basis in the site occupation a large proportion of the final states coalesce in a few states. We have tested several different strategies to obtain as many excited states as possible. The best we have found is the following: 1. The initial state is a single component state on the noninteracting basis. Up to system size L 16 we run over all possible initial states. For larger systems, we choose initial states at random. 2. Successive displacement transformations are applied to the state, choosing the minimum absolute value of the strength λ that minimizes or maximizes the energy. Transformations that increase the variance of the energy are rejected.
As the calculation of λ, for a given X, involves a CPU time similar to the transformation of the state, it is convenient to choose first the operators X that we think are likely to be more important, i.e., that correspond to larger values of |λ|. We order the operators X according to |V α,β,c,η /ΔE|, where ΔE is the internal energy change, i.e., the energy of the states corresponding to creation operators in X minus the energy of the states corresponding to annihilation operators: We have only considered fourth order X operators. after all displacement transformations for fourth order operators has been performed, we consider again the same set of operators, which now produce displacement transformations with much smaller values of the strength λ. The whole procedure is repeated three times and we have checked that after this, convergence is basically achieved.
The displacement transformations that increase the energy variance and so are rejected depend on the state considered. They appear at late stages, when the rotations of the corresponding operators have already been performed at least once.
We can express a general state |Ψ〉 as a linear combination of configurations (Slater determinants) |Φ Θ 〉 in a given basis: Each configuration is represented by an integer, whose 0's and 1's in its binary form corresponds to the empty and occupied states of the configuration. A state is represented by a set of configurations and the corresponding weights a Θ .
In the rest of the paper, we consider a nearest neighbour interaction with U 1, which sets our unit of energy, a transfer energy t 1/2. We fix the total number of particles to half the number of sites. For this model, the ground state is always localized, while we expect to have a delocalization transition at infinite temperature for a critical disorder around W c ≈ 4 [22][23][24][25] in our units, although it could be a larger value [26][27][28].
For small system sizes, we can generate all the final states starting from a complete initial standard basis and compare the results with exact diagonalization. We start from a single component state in the non-interacting basis, i.e., a Slater determinant for a given set of one-particle states. We then apply displacement transformations following the strategy described above and store the final energy of this evolved state. We arrange these energies E in ascending order and compare them with the exact energies, E exact . The results are shown in Figure 1A for two values of the disorder W 5 (blue) and W 3 (orange) and system size L 10. The larger disorder corresponds to the localized regime, while the smaller disorder to the extended regime. The black dashed line represents E E exact and both blue and orange curves are pretty close to it. We do not appreciate large horizontal steps, which would indicate that many initial states have collapsed into the same final state. The steps seen in Figure 1A involve a few states, while with other approaches there are large proportions of states coalescing in the same final one. In the lower part of the figure, we plot E − E exact as a function of E exact to appreciate better the precision of our method. Similar results to those in Figure 1A are obtained when working in the site basis. But starting from wavefunctions corresponding to Slater determinants of non interacting wavefunctions. The problem with this approach is the very large size of the wavefunction even at the very beginning of the calculation.
In order to check the quality of the wavefunctions obtained, rather than of the energy, we have calculated the inverse participation ratio (IPR) in the site basis, I s . The IPR is defined as where a Θ are the components of the state, see Eq. 16. In this case, |Φ Θ 〉 correspond to a real space configuration, i.e., the set of occupied and empty sites. In Figure 1B, we compare the IPR Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 767001 of the states obtained in panel (a) for the exact procedure and our algorithm. The agreement for most states is good, but a few ones, mainly the least localized, present appreciable differences.
We have studied the level spacing distribution for system size L 10. We have normalized the data to the average spacing, but not to the energy dependent average spacing, since we are only interested in its behaviour at small values. The results are presented in Figure 1C for W 5 (blue) and W 3 (orange) on a semilogarithmic scale. As high order displacements have not been considered, we expect the distribution to be Poisson like, independently of whether we are in the localized or in the extended regime. This is not in contradiction with Figure 1A, since very small deviations in the energies can cause an appreciable change in the level spacing. If several states were to collapse in one, we would see a peak in the level spacing distribution at small spacings [14]. No such peak can be observed in Figure 1C. To analyse this question more quantitatively, we have fitted from the second to the fifth points by a straight line, shown as a dashed line in the figure, and we see that the first point does no deviate upwards from the fitted line.  We have also analysed the behaviour of our algorithm for a larger system, L 20. For a realization of the disorder, we consider many different initial states, i.e., Slater determinants constructed from different one-electron states, and monitor the evolution of their energy E and their energy standard deviation σ in terms of the number of displacement transformations performed. In Figure 2A we plot the evolution of the energy E for a few states as a function of the number of displacement transformations for a disorder W 5. The shaded area around each line corresponds to E ± σ/2, and the final energy should be within this interval. The energy uncertainty in the extended regime (W 3) is much larger than in the localized regime (W 5), but in neither case we observe a tendency for states to merge together, at least, at a significant rate. To quantify this question, we analyse the energy level spacing as we did for L 10. The results are shown in Figure 2C, and as before we do not appreciate a peak at very low energy spacings above the standard tendency for a Poisson distribution.

LOCALIZATION
To study localization, we start from a single configuration state in the non-interacting basis and transform it with displacement transformations as described before. At each step, we monitor the energy variance and the IPR in the non-interacting basis, I ni . The IPR is defined as in Eq. 17, but now configurations refer to occupied and empty non-interacting states. I ni quantifies the spread of the wavefunction due to interactions. To compare our results with the IPR in the site basis, we will have to take into account the increase of the single-particle wavefunction with decreasing disorder, which is a smooth dependence.
In Figure 3A, we plot in a log-log scale the evolution of I ni for several states as displacement transformations are performed on them. In the horizontal axis we represent the energy variance of the system normalized to its initial value. Each line corresponds to the evolution of a state, and each dot to the coordinates of the state after a number of displacement transformations has been performed. Curves in orange are for W 5 and in blue to W 3. System size is L 20. As displacement transformations are performed, curves evolve from right to left. The program stops when 500 transformations have been attempted (at the beginning not all transformations succeed in changing the state). In order to speed up the process, only terms larger (in absolute value) than its maximum component multiplied by 10 −4 have been kept for both the Hamiltonian and the wavefunction. We note that while orange curves have basically reached an asymptotic value, most of the blue curves are still far from saturation.
In Figure 3B, we plot the median of I ni over 1,000 states for each of the disorders shown in the legend as a function of the normalized energy variance for size L 20. The shaded region ranges between the 25% and the 75% of the IPR. The tendency is clear, as disorder decreases states become more extended and their IPR's increase.
To quantify the tendency of a state to become extended, we have calculated the slope of the curves of I ni versus the normalized variance in double logarithmic scale, as those shown in Figure 3A, at different positions and have kept the maximum (absolute) value, which we denote as δ. In Figure 4A we show the histogram of this maximum slope δ for several values of W and size L 20. Each histogram is calculated considering 600 disorder realizations and one random initial state for each realization. These histograms present a peak that increases in height and decreases in position as disorder increases. The evolution of the histograms with system size may be an adequate tool to characterize localization.
In Figure 4B we plot the average value of δ as a function of disorder for several system sizes. For disorders 2 and 3, δ drastically increases with system size, as can be seen in Figure 4C, where we represent the system size dependence of δ for the four values of the disorder considered. This is a promising tool, but we have to improve it in order to get quantitative results for the MBL transition. We are currently working on the size dependence of δ and on the relationship between the action of the displacement transformations and the dynamics.
We have also explored the effect of a cutoff in the nondiagonal terms included in the Hamiltonian on the IPR. We order these terms according to |V α,β,c,η /ΔE| and only keep terms with this quantity larger than a cutoff V c . The results are presented in Figure 5A, where we plot the median of the distribution of I ni as a function of the cutoff for several disorders and for L 20. We use the median of the distribution to avoid the bias in the most delocalized states. We note that for cutoffs smaller than 0.1 the IPR remains constant. So, as long as the IPR is concerned, a cutoff of V c 0.1 does not affect the results, while it drastically reduces computing time.
We have calculated the IPR for several system sizes and values of the disorder using the previous cutoff V c 0.1. In Figure 5B we represent the results as a function of the Hilbert space, as previously done in Ref. [23], on a double logarithmic scale. Our results are consistent with those in Ref. [23] (A precise comparison is not possible, since boundary conditions are different.) We note that, within our approximations, we can reach larger system sizes.

DISCUSSION AND OUTLOOK
Displacement transformations are useful to diagonalize interacting Hamiltonians. We have explored their application in two different procedures. When displacement transformations  rotate the state of the system, the basis employed remain fixed along the whole process, while when one rotates the Hamiltonian, the basis is updated self-consistently at every step [14]. Rotating the Hamiltonian also has the advantage of being specially suited for renormalization techniques, but these two advantages have to be confronted with the more efficient storage of the states when these are rotated, instead of the Hamiltonian. The present implementation of the method also has the advantage of degradating its performance less with increasing system size. We believe this is due to the fact that the necessary cutoffs in the present implementation are more under control than the cutoffs in the previous one.
With the current procedure, we can reach system sizes much larger than with exact diagonalization. In this work, we have calculated up to L 32 and we expect to reach larger sizes in the near future, at least for strong disorders, W > 4. DMRG and other approximate methods can reach large system sizes, but our approach produces large fractions of relatively unbias excited states as shown in a previous work [14]. We have checked that in the present implementation (rotating the wave-function instead of the Hamiltonian) we still reach most of the excited states.
Storing the displacement transformations utilised to get one state, D X (λ), one can construct an approximation to the IOM and rewrite the final state in terms of them. This is clearly an approximation, since we are only including transformations for a single state. It can be improved by considering other states.
In the site basis, the Hamiltonian takes a very simple form, but the size of the wavefunction, even the non interacting one, grows drastically with system size and renders the problem in practical for large sizes. For the calculation of the ground state, this approach is similar to that based in the ARPACK library.
Our procedure can be directly extended to any number of dimensions or to random geometries. It can also be applied to long-range interactions, and we intend to study MBL in this case. The method should be suited for molecular orbital calculations, nuclear structure determination [29], Rydberg atoms [30] and excited states in molecules [31,32]. In the future, we would like to study the possibility of determining the localization properties from the evolution with system size of the distribution of δ.
Our method is equivalent to the Hartree-Fock approximation when displacement transformations involving only a single pair of creation and annihilation operators are performed. Including more complex transformations represents the natural extension of Hartree-Fock for the construction of non-gaussian states.
Displacement transformations constitute an excellent tool to study quantum interacting systems. In this paper, we have shown that they can be used to obtain excited states and to study their localisation properties in an efficient way.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.