Construction of Many-Body Eigenstates with Displacement Transformations

TL;DR

This paper introduces a practical method using displacement transformations to compute local integrals of motion, enabling accurate construction of many-body eigenstates beyond Gaussian approximations, including excited states.

Contribution

The authors develop a systematic approach combining energy minimization and exact diagonalization to efficiently compute many-body eigenstates and local integrals of motion.

Findings

High-precision ground state calculations for large systems

Efficient construction of multiple excited states

Extension of Hartree-Fock and configuration interaction theories

Abstract

Many-body eigenstates beyond the gaussian approximation can be constructed in terms of local integrals of motion (IOM), although their actual computation has been until now a daunting task. We present a new practical computation of IOMS based on displacement transformations. It represents a general and systematic way to extend Hartree-Fock and configuration interaction theories to higher order. Our method combines minimization of energy and energy variance of a reference state with exact diagonalization. We show that our implementation is able to perform ground state calculations with high precision for relatively large systems. Since it keeps track of the IMO's forming a reference state, our method is particularly efficient dealing with excited states, both in accuracy and the number of different states that can be constructed.

Figures

Equations

$$H=\sum _X{V}_X(X^{†}+X)$$

$$X=n_{{\alpha }_1}\cdots n_{{\alpha }_k}{c}_{{\beta }_1}^{†}{c}_{{\gamma }_1}{c}_{{\beta }_2}^{†}\cdots c_{{\gamma }_l}$$

$$\tilde{H}=\sum _i{ϵ}_{\alpha }{\tau }_{\alpha }+\sum _{\alpha ,\beta }{ϵ}_{\alpha ,\beta }{\tau }_{\alpha }{\tau }_{\beta }+\sum _{\alpha ,\beta ,\gamma }{ϵ}_{\alpha ,\beta ,\gamma }{\tau }_{\alpha }{\tau }_{\beta }{\tau }_{\gamma }+\cdots$$

$$|\Phi ⟩={\tilde{c}}_{{\beta }_1}^{†}\cdots {\tilde{c}}_{{\beta }_N}^{†}|0⟩=U^{†}{c}_{{\beta }_1}^{†}\cdots c_{{\beta }_N}^{†}|0⟩.$$

$${\mathcal{D}}_X(\lambda )=\exp \{\lambda (X^{†}-X)\},$$

$${\mathcal{D}}_X(\lambda )=1+\sin (\lambda )(X^{†}-X)+(\cos (\lambda )-1)(X^{†}X+X{X}^{†}).$$

$${\mathcal{D}}_X(-\lambda )H(\lambda ){\mathcal{D}}_X(\lambda )$$

$${\mathcal{D}}_X(-\lambda )H(\lambda ){\mathcal{D}}_X(\lambda )=H+\lambda [H,A]$$

$$+\frac{{\lambda }^2}{2!}[[H,A],A]+\frac{{\lambda }^3}{3!}[[[H,A],A],A]+\cdots$$

$$\tan (2\lambda )=\frac{2V_X}{\Delta {ϵ}_X}$$

$$\Delta {ϵ}_X$$

$$V_X\tan (\lambda )(X^{†}X-X{X}^{†})$$

$$\{{\mathcal{D}}_X{c}_i^{†}{\mathcal{D}}_{-X},{\mathcal{D}}_X{c}_j{\mathcal{D}}_{-X}\}=\{c_i^{†},c_j\}={\delta }_{i,j}.$$

$$\frac{\mathrm{d}⟨{\Phi }_0|{\mathcal{D}}_X^{†}(\lambda )H{\mathcal{D}}_X(\lambda )|{\Phi }_0⟩}{\mathrm{d}\lambda }\equiv \frac{\mathrm{d}⟨H(\lambda )⟩}{\mathrm{d}\lambda }=0$$

$$⟨H(\lambda )⟩$$

$$-⟨AHA⟩-⟨H⟩=⟨{\Phi }_X|H|{\Phi }_X⟩-⟨{\Phi }_0|H|{\Phi }_0⟩\equiv \Delta E_X$$

$$V_{Y,Z}\equiv ⟨{\Phi }_0|Y^{†}HZ|{\Phi }_0⟩$$

$$⟨{\Phi }_0|[H,A]|{\Phi }_0⟩=-2V_{X,1}$$

$$⟨H(\lambda )⟩=⟨H⟩-\sin (2\lambda )V_{X,1}+{\sin }^2\lambda \Delta E_X$$

$$\tan 2{\lambda }_X=\frac{2V_{X,1}}{\Delta E_X}.$$

$${\sigma }^2\equiv ⟨{\Phi }_0|H^2|{\Phi }_0⟩-⟨{\Phi }_0|H|{\Phi }_0{⟩}^2$$

$$⟨H^2(\lambda )⟩$$

$$⟨{\Phi }_X|H^2|{\Phi }_0⟩=\sum _n⟨{\Phi }_X|H|{\Phi }_n⟩⟨{\Phi }_n|H|{\Phi }_0⟩$$

$$⟨{\Phi }_X|H^2|{\Phi }_X⟩-⟨{\Phi }_0|H^2|{\Phi }_0⟩=$$

$$E_X^2-E_0^2+\sum _{Z_n\ne 1,X}\left(V_{X,Z_n}^2-V_{Z_n,1}^2\right)$$

$$\sigma (\lambda {)}^2$$

$$H=\sum _i{ϵ}_i{n}_i+\sum _{⟨i,j⟩}t{c}_j^{†}{c}_i+\frac{1}{2}\sum _{⟨i,j⟩}U{n}_j{n}_i$$

$${\sigma }_{L/2}^2=⟨{\Phi }_0|n_{L/2}^2|{\Phi }_0⟩-⟨{\Phi }_0|n_{L/2}|{\Phi }_0{⟩}^2,$$

$$t{c}_L^{†}{c}_{L+1}+U{n}_L{n}_{L+1}+h.c.$$

$${\tilde{n}}_{\delta }$$

Full text

I Introduction

Integrals of motion (IOM) have gained a renewed interest in the context of disordered many-body systems. In the non-interacting Anderson insulator Anderson:1958fz in $d=1,2$ dimensions the single particle wavefunctions are exponentially localized for any degree of disorder. In the presence of interactions particles will remain localized, which is known as many-body localization (MBL) for strong disorderBaskoAA ; Nandkishore:2015kt ; Znidaric:2008cr ; Pal:2010gr ; Bardarson:2012gc ; DeLuca:2013ba ; Abanin2017 ; Alet . It was realized that MBL can be understood through the existence of an extensive number of exponentially localized IOMs Chandran ; Huse:2014co ; Serbyn:2013cl . Inevitably, this observation led to a rush of new methods to compute the IOMs in the MBL-phase,2015NuPhB.891..420R ; Imbrie2016 ; YouQiXu2016 ; 2016arXiv160707884P ; 2016arXiv160801328G ; 2016arXiv160803296O ; 2016arXiv160300440I ; Friesdorf:2015dd ; 2016arXiv160609509H ; ImbrieReview ; Peng2019 ; Varma2019 ; Mierzejewski . We published our own method based on displacement transformations of the Hamiltonian written in second quantization, which was one of the first practical computational approaches Rademaker:2016jf ; ROS17 .

The presence of these localized IOMs in the fully many-body localized phases prevents thermalization. The question of whether a many-body quantum system thermalizes has been cast into the Eigenstate Thermalization Hypothesis (ETH):Deutsch:1991ju ; 1996JPhA…29L..75S ; Srednicki:1994dl ; Rigol:2011bf the expectation value of any local observable in an eigenstate with a given energy density is equal to its expectation value in the Gibbs ensemble with corresponding temperature. If, however, there exist local density IOMs (that can be expressed as the sum of local operators) the corresponding thermal state will be a so-called generalized Gibbs ensemble.Rigol:2006df It has been shown that any finite-ranged translationally invariant Hamiltonian will thermalize towards their corresponding generalized Gibbs ensemble.2015arXiv151203713D Therefore, whether and how a system thermalizes is directly related to the structure of its integrals of motion.

The existence of IOMs has a deep consequence for the structure of many-body eigenstates. Written in the basis of IOMs, all many-body eigenstates are just product states defined by the occupation of each IOM. Given the basis transformation $U$ from the physical basis to the basis of IOMs, eigenstates can be written as $U|\psi_{SD}\rangle$ where $|\psi_{SD}\rangle$ is a Slater determinant of the form $c_{\beta_{1}}^{\dagger}\cdots c_{\beta_{N}}^{\dagger}|0\rangle$.

Though an accurate computation of all IOMs is in practice very tedious, we propose here that one can focus on a specific reference (product) state. By only considering those displacement transformations that affect that reference state, we increase the efficiency and thus accuracy of the numerical implementation. In this work we show that this indeed works, and we find accurate results for ground states and excited states for systems up to $L=30$ sites. In the lowest possible order, the method can work up to much larger sizes.

White White introduced a similar method to ours to diagonalize the Hamiltonian through what he calls canonical transformations (we will see that our displacement transformation are canonical), sharing many ingredients with our approach, although no explicit mention of the IOM and their general applicability was made. Both methods deal directly with the Hamiltonian, instead of the wavefunctions, unlike most other numerical approaches. They manipulate the second quantized Hamiltonian, which is a sum of abstract operator terms with the corresponding numerical coefficient. This presents several advantages. For example, one can efficiently separate the low-energy and high-energy orbitals from the strongly correlated intermediate ones. Then, fixing the occupancy of low-energy levels to 1 and of high-energy levels to 0, one can remove them from the problem. A second advantage is that basically the same approach is used for the ground state and for excited states.

The so called flow equation methodWegner ; Wegner2 and similarity renormalizationWilson are continuous versions of unitary transformations, in which a set of differential equations is solved to implement the displacement transformations.2016arXiv160707884P ; Monthus2016 ; Savitz ; Monthus3 ; Monthus4 ; Thomson Continuous transformations have the advantage of commuting one with each other, while discrete implementations take into account the term to be transformed at all orders in a single step.

The remainder of the paper is set up as follows. In Sec. II and III we briefly review the concept of IOMs and our method of displacement transformations. The increase in efficiency gained by focusing on a specific reference state is explained in Sec. IV. We proceed with comments on the practical implementation (Sec. V) and some numerical results for a one-dimensional interacting disordered system in Sec. VI. In Sec. VII we discuss possible applications and related work.

II General quantum Hamiltonian

We consider a general fermionic quantum Hermitian Hamiltonian, which can always be written in the form

[TABLE]

with the terms $X$ of the form

[TABLE]

where $c_{i}^{\dagger}$ ($c_{i}$) is the fermionic creation (annihilation) operator of state $i$ and $n_{i}=c_{i}^{\dagger}c_{i}$ is the number operator, and all $\alpha,\beta,\gamma$ are different. We will assume particle conservation and so the number of creation and annihilation operators will be the same in each term of the Hamiltonian. Terms containing creation and annihilation operators are called non-diagonal, while terms composed of density operators only are called diagonal or classical.

In principle, our goal is to find a unitarity transformation $U$ that transforms the Hamiltonian $H$ into a classical Hamiltonian $\tilde{H}=U^{\dagger}HU$, i.e., one that only contains density operators. This ”diagonalized” Hamiltonian can be written in the form

[TABLE]

where $\tau_{\alpha}=U^{\dagger}n_{\alpha}U=\tilde{c}_{\alpha}^{\dagger}\tilde{c}_{\alpha}$ are the transformed number operators and correspond to the integrals of motion (IOMs). Now $\tilde{c}_{\alpha}^{\dagger}$ and $\tilde{c}_{\alpha}$ are the transformed creation and annihilation operators. Any many-body eigenstate can be constructed by application of the transformed creation operators on the vacuum state $|0\rangle$,

[TABLE]

In the basis of IOMs, therefore, all many-body eigenstates are ‘product states’ in the transformed creation operators. This also implies that any Slater determinant of the original fermions becomes a many-body eigenstate under the unitary transformation $U$.

III Displacement transformations

Two of us showed earlierRademaker:2016jf that the Hamiltonian given by Eq. (1) can be diagonalized by a unitary transformation consisting of a product of displacement transformations. A displacement transformation is defined as

[TABLE]

where $X$ is any given quantum term. It is trivial to show that ${\cal D}_{X}(\lambda)$ is equal to

[TABLE]

The fact that the displacement transformation can be summed in a compact way has important practical consequences, as we will see.

To start with, any operator transformed under ${\cal D}_{X}(\lambda)$ can be cast in a closed expression. Defining $A\equiv X^{\dagger}-X$ and focusing in the Hamiltonian $H$, we arrive at

[TABLE]

This is in contrast with the typical infinite series expansion of the transformed operator in terms of commutators of $A$ and the operator:

[TABLE]

Under a displacement transformation ${\cal D}_{X}(\lambda)$, the prefactor of the term $X^{\dagger}+X$ in the new Hamiltonian is changed and depends on the prefactors in the old Hamiltonian of the term itself and of the density operators corresponding to the creation and annihilation operators of $X$. Both contributions can cancel each other if the strength of the transformation $\lambda$ is chosen as (see Appendix A)

[TABLE]

where $V_{X}$ is the coefficient of the operator $X^{\dagger}+X$ in the Hamiltonian and $\Delta\epsilon_{X}$ is the energy difference between the two configurations involved in the transition $X$. That is, if $X$ is given by Eq. (2), then

[TABLE]

where $\epsilon_{\alpha_{1},\beta_{1}}$ is the coefficient of the term $n_{\alpha_{1}}n_{\beta_{1}}$ in the Hamiltonian.

After transforming the Hamiltonian under ${\cal D}_{X}(\lambda)$, with $\lambda$ given by Eq. (9), the operator $X$ only remains in the classical combination

[TABLE]

It is interesting to note that the condition of canceling out the quantum term is equivalent to extremizing (maximizing/minimizing) the coefficient of the classical term, since the derivative with respect to $\lambda$ of its contributions is equal zero when Eq. (9) is satisfied (see Appendix A). This is just a consequence of the conservation of the trace of $H^{2}$ under the displacement transformations.

One can calculate any observable by constructing the corresponding operator and transforming it with the same displacement transformations that diagonalize the Hamiltonian. As the many-body states take a simple form in this new basis, Eq. (4), it is easy to compute the expectation value and the variance of the operator.

One can in principle diagonalize the interacting Hamiltonian by the application of successive displacement transformations, starting with those involving operators of the form $c_{j}^{\dagger}c_{i}$ and continuing with higher (larger number of terms) order operators. In practice this method generates and exponentially large number of terms, as was shown in our earlier workRademaker:2016jf ; ROS17 . As we will show in the following sections, it is more efficient to concentrate on certain expectation values with respect to a given state. This allows us to study with more accuracy certain parts of the spectrum.

We pause here to note that our displacement transformations are canonical transformations, as defined in (6). A transformation is canonical if when applied to the creation $c^{\dagger}_{i}$ and annihilation $c_{j}$ operators preserves their anticommutation relations. Our displacement transformations ${\cal D}_{X}$ are indeed canonical transformations since

[TABLE]

In fact our displacement transformations are the main type of canonical transformation considered in Ref. White, .

IV Reference States

IV.1 Energy minimization

Since the diagonalization of the full Hamiltonian is a huge task, it is more efficient to focus on a given state $|\Phi_{0}\rangle$ and perform only those transformations that affect the expectation value or the variance of the Hamiltonian with respect to this state. This is the core novel development of this work.

Let us begin by choosing transitions $X$ and their corresponding $\lambda$ so that the displacement transformation defined by them minimizes (or maximizes) the expectation value of the energy

[TABLE]

where our state $|\Phi_{0}\rangle$ is of the form $|\Phi_{0}\rangle=c_{\beta_{1}}^{\dagger}\cdots c_{\beta_{N}}^{\dagger}|0\rangle$, where $|0\rangle$ is the vacuum and the set of indices $\beta_{i}$ is fixed and define our state. The creation and annihilation operators are self-consistently redefined at each step (at each displacement transformation). Strictly speaking we minimize the energy if we are considering the ground state, while we find a local extremum of the energy if we are dealing with an excited state. If we are interested in the ground state, $|\Phi_{0}\rangle$ is our running estimate of it and the indices $\beta_{i}$, instead of being fixed, have to be chosen to minimize the expectation value of the energy.

Taking expectation values with respect to state $|\Phi_{0}\rangle$ in (III) we arrive at

[TABLE]

In order to have $A|\Phi_{0}\rangle\neq 0$, all density operators in $X$ must correspond to occupied states in $|\Phi_{0}\rangle$, while creation operators must have the opposite occupancy than annihilation operators in $|\Phi_{0}\rangle$. Displacement transformations of operators $X$ not satisfying these two conditions do not have to be considered as long as the state $|\Phi_{0}\rangle$ is involved. This is at the heart of the increase in efficiency of our implementation.

For the transformations affecting the expectation value of $H$, we have $\langle A^{2}\rangle=-1$ and substituting in (IV.1) we get

[TABLE]

We have the freedom to choose without loss of generality the operator $X$ (as opposed to $X^{\dagger}$) such that $X|\Phi_{0}\rangle\neq 0$, and we will denote the resulting state as $|\Phi_{X}\rangle\equiv X|\Phi_{0}\rangle$. Then

[TABLE]

This quantity is the energy difference between the states $|\Phi_{X}\rangle$ and $|\Phi_{0}\rangle$ and so only depends on the classical terms of $H$. Let us define the expectation values of the Hamiltonian multiplied by different operators $Y$ and $Z$ as

[TABLE]

From now on we will consider a real Hamiltonian, and so

[TABLE]

where ${\bf 1}$ denotes that the unity operator is acting on $|\Phi_{0}\rangle$. We note that all terms of the Hamiltonian with the same quantum part as $X$ and with density operators of states occupied in $|\Phi_{0}\rangle$ will contribute to $V_{X,{\bf 1}}$. With the previous definitions, $\langle H(\lambda)\rangle$ becomes

[TABLE]

And the value of $\lambda$ minimizing (or maximizing) $\langle H(\lambda)\rangle$ is

[TABLE]

This equation has two solutions for $\lambda_{X}$ differing by $\pi/2$. They correspond to the two possible ways to associate the bonding and antibonding states with the original states. In order to maximize the overlap with the initial state, it is convenient to choose the solution satisfying $|\lambda|<\pi/4$.

One can transform $H$ reiteratively with different displacement transformations up to a given order until all coefficients $V_{X,1}$ are smaller than a given cutoff.

Considering only second order displacement transformations, our method is equivalent to the standard Hartree-Fock approximation. Two important advantages of our method is that it can naturally be extended by considering higher order transformations and that it can be applied efficiently to excited states, as we will see.

IV.2 Energy variance minimization

The previous scheme of getting extreme values for the energy works fairly well for ground state calculations and in some cases for excited states, but not in general. Considering the MBL problem, for example, energy minimization (maximization) is efficient in the localized regime, but not so in the extended regime, where it does not converge for some states, that end up oscillating in a cyclic way between different configurations. The action of a given transformation is undone by the action of two other transformations maximally overlapped with the first one. The problem arises because we are not minimizing any function. The energy is an extreme for all eigenvalues, but not necessarily a minimum, except for the ground state. We can get around this problem by minimizing the variance of the energy

[TABLE]

This procedure is practical if we are able to calculate the value of $\lambda$ that minimizes $\sigma^{2}$ without having to obtain the operator $H^{2}$ explicitly. Under a displacement transformation $\langle H^{2}(\lambda)\rangle$ is given, in analogy with Eq. (19), by

[TABLE]

Inserting the identity operator $\sum_{n}|\Phi_{n}\rangle\langle\Phi_{n}|$, where $|\Phi_{n}\rangle$ are a basis of the many-body space, which can be obtained by applying all possible combinations of $N/2$ creation operators, we can rewrite $\langle\Phi_{X}|H^{2}|\Phi_{0}\rangle$ as

[TABLE]

where $Z_{n}$ is the operator that takes the system from $|\Phi_{0}\rangle$ to $|\Phi_{n}\rangle=Z_{n}|\Phi_{0}\rangle$ and we can assume that it does not contain density operators. With a similar procedure, we get

[TABLE]

Substituting Eqns. (22-24) in (21) we arrive at

[TABLE]

This expression can be minimized with respect to $\lambda$ numerically. We note that the terms in the third line of Eq. (IV.2) contain the same parameters that appear in the expression for the energy, Eq. (19), and produce the same contribution to $\sigma(\lambda)^{2}$ at both energy minima (differing in $\lambda_{X}$ by $\pi/2$). The terms in the first and second line of Eq. (IV.2) have half the period than the other terms and, in general, unbalance the contribution to $\sigma(\lambda)^{2}$ in the two energy minima. This implies that when minimizing the variance we cannot restrict $\lambda$ to values satisfying the condition $|\lambda|<\pi/4$.

When obtaining excited states, we can get extremes of the energy, but with the constraint that the energy variance never increases. In this way the procedure is stable. Once there are no displacement transformations that minimize/maximize the energy up to a given cutoff, we can minimize the energy variance, which results in a further diagonalization of the Hamiltonian not affecting directly the reference state.

IV.3 Perturbation treatment of the remaining Hamiltonian

In practice it is necessary to establish a cut-off for the strength of the transformations $\lambda$ (or for the coefficients of the off-diagonal terms of the Hamiltonian) to be performed. We then apply perturbation theory to the remaining small off-diagonal Hamiltonian and to any other operator we may want to calculate. However, the level spacing of many-body states is so small that direct perturbation theory becomes impractical due to the uncertainty in the energy denominators. A convenient alternative is to ‘project’ the Hamiltonian onto a subspace close to the reference state $|\Phi_{0}\rangle$ by constructing the matrix formed by the elements $V_{Y,Z}$ where $Y$ and $Z$ are (purely quantum) operators of a given order. This matrix is then diagonalized.

If we are calculating the expectation value of some other operator $\cal O$, we construct the corresponding matrix $\langle\Phi_{0}|Y^{\dagger}{\cal O}Z|\Phi_{0}\rangle$ and rotate it with the unitary transformation $U_{1}$ that diagonalizes the projected Hamiltonian matrix. The element $\langle\Phi_{0}|U_{1}^{\dagger}{\cal O}U_{1}|\Phi_{0}\rangle$ is then our estimate of the expectation value of ${\cal O}$.

This final matrix diagonalization of the ‘projected’ Hamiltonian is similar to the so-called configuration interaction approach, usually based on a Hartree-Fock basis. The advantage of our procedure is that we can go naturally beyond the Hartree-Fock method by performing displacement transformations of higher order. Furthermore, using the displacement transformation we have constructed a nearly diagonal Hamiltonian, even before the final diagonalization procedure. Note that since Hartree-Fock methods construct so-called gaussian states, our method allows for the construction of non-gaussian states.Nooijen2000

as our method controls the IOM’s explicitly, it is possible to carefully select which states are included in the final diagonalization, depending on the quantities to be calculated or the nature of the problem. For example, one could selesct states closest in energy to the reference state, instead of selecting in terms of the number of excitations. Or alternatively, one could try to keep states that are more likely to affect the operator to be calculated.

IV.4 Finite density

To overcome as much as possible the problem of the large amount of high order terms generated by the displacement transformations, it is highly convenient to consider a finite density of particles and neglect terms according to the number of excitations, instead of the number of particles. We consider as the ”vacuum” our reference state $|\Phi_{0}\rangle$, which as we have seen can be either the ground state or any excited state. Creation and annihilation operators for occupied states in $|\Phi_{0}\rangle$ are transformed into annihilation and creation operators, respectively, of the corresponding hole excitations. Operators for unoccupied states are kept the same. We keep terms containing up to a given number of these new operators.

V Practical implementation

In order to test the viability of our method, we have applied it to the important problem of many-body localization. Consider the following interacting Hamiltonian for spinless fermions on a one-dimensional chain of length $L$,

[TABLE]

where $c_{i}^{\dagger}$ is the creation operator on site $i$ and $n_{i}=c_{i}^{\dagger}c_{i}$ is the number operator. We consider a nearest neighbor interaction with $U=1$, which sets our unit of energy, a transfer energy $t=1/2$ (to use the same values as for spins models) and a disordered site energy $\epsilon_{i}\in[-W,W]$. We assume the total number of particles is equal to half the number of sites.

In our code, each term of the Hamiltonian is represented by a real coefficient corresponding to its strength and by a long integer with every two bits indicating every site operator (creation, annihilation, density or no operator). Operator multiplication is a very slow process to implement numerically. To overcome this shortcoming and taking into account that a displacement transformation $\mathcal{D}_{X}$ only modifies the operators for the sites involved in $X$, we construct a table with the outcome of the transformation for any possible combination of operators in the sites involved in $X$. Given a general product of site operators, we extract the ones at these positions, look for the outcome in the table and insert the new operators in the original product.

For this model, the ground state is always localized, while we expect to have a delocalization transition at infinite temperature (i.e., for states chosen at random) for a critical disorder around $W=3.5$ in our units.

To evaluate the degree of localization of the eigenstates of the Hamiltonian, we calculate the variance of the density operator at site $L/2$, $n_{L/2}$. We apply to this operator the same transformations as to the Hamiltonian. It is always expressed in our current basis. When the procedure is finished, it is easy to evaluate the variance of this operator with respect to our reference state $|\Phi_{0}\rangle$,

[TABLE]

since the operator is already in our final basis, where $|\Phi_{0}\rangle$ is a single product of creation operators acting on the vacuum.

We have performed transformations up to fourth order and we have kept terms in the Hamiltonian containing up to two particle and two hole excitations with respect to the reference state.

V.1 Ground state

To obtain a state as close as possible to the ground state we choose the following strategy. We first get rid of quantum terms of second order by performing all possible displacement transformation of this order (this is equivalent to be in the one-electron basis). Then we choose the set of creation operators $\{\beta_{i}\}$ such that the state $|\Phi_{0}\rangle=c_{\beta_{1}}^{\dagger}\cdots c_{\beta_{N}}^{\dagger}|0\rangle$ minimizes the energy of the present Hamiltonian. $|\Phi_{0}\rangle$ is our reference state and we calculate $V_{X,{\bf 1}}$ and $\Delta E_{X}$ for all possible $X$, up to a given order, and choose the one that reduces most the energy of this state. After performing the corresponding displacement transformation, we check whether a different set of creation operators reduces further the energy of the new rotated Hamiltonian. This step is cyclically repeated until no displacement transformation can reduce the energy any further. A third stage looks for transformations that reduce the energy variance, keeping in this case fixed the set of creation operators defining our state. We finally ‘project’ the Hamiltonian over configurations differing from the reference state by one or two electron-hole excitations and diagonalize it.

V.2 Excited states

The basis of IOMs provides an unambiguous way to specify a many-body eigenstate, since each can be labeled by a specific choice of the IOMs. A great advantage of our method is that, as we will see, we can get a large proportion of the excited states of the system, many more than with other available techniques. To construct an excited state, we select a given set of IOMs (i.e., a set of creation operators $\{\beta_{i}\}$) and keep performing displacement transformations that minimize (or maximize) the expectation value of the energy with respect to this state. Unlike for ground state calculations, the set of creation operators defining our reference state is kept fixed. By varying systematically the set of IOMs chosen, we can get all the excited states of the system for small enough sizes. For large sizes we choose a random set of IOM to construct a state.

There are many possible minimization strategies and which is the most suitable may depend on the quantity to be studied. Here we have concentrated in trying to obtain as many different excited states as possible, i.e., to avoid that states end up collapsing in the same final state. We found that a good strategy to this end is to perform first the transformations that most increase (or reduce) the energy, with the constraint that $|\lambda|<\pi/4$ and provided that they do not increase the variance of the energy $\sigma^{2}$. Once there are no transitions changing the energy, we reduce the energy variance as much as possible. Reducing the variance from the beginning tends to effectively interchange the IOM indices and may change the reference state, resulting in a bias in our choice of excited states.

VI Numerical results

VI.1 Ground state

Let us first focus on the error of our estimate of the ground state energy as compared with results from exact diagonalization. In Fig. 1 we present the geometrical mean of the absolute error in the ground state energy per site as a function of disorder for several system sizes and two implementations of the method. Empty symbols and dashed lines correspond to results for 2nd order calculations, while solid symbols and continuous curves to results for transformations up to fourth order. In both cases the residual Hamiltonian was projected over configurations differing from the ground state by one or two electron-hole excitations and then diagonalized. System sizes are 12, 14 and 16. In all cases, the cutoff strength for $\lambda$ was $10^{-3}$.

We see in Fig. 1 that the errors are very small for the (small) system sizes that we can diagonalize exactly. The error increases with system sizes, but does so relatively slowly. Results for second and forth order are almost comparable, probably due to the final diagonalization performed, which for these system sizes can deal with all one and two particle excitations.

Next, we studied the variance of the occupation of the site at position $L/2$ in order to test the accuracy of our method for other quantities, beside the energy, and to prove the applicability of the method to relatively large system sizes. The value of the variance should be 0 in the strongly localized limit and 1/4 in the completely extended case, since we are considering 1/2 occupation. In Fig. 2 we plot this variance as a function of disorder for several system sizes ranging from 12 to 30. Up to size 16, exact results are also plotted (thick lines). Results for ground state calculations, in which we focus for the moment, correspond to the lower set of curves. Our calculation has been done up to 4th order with diagonalization of the final Hamiltonian projected over one and two particle excitations. For large system sizes we only include the first 8000 configurations ordered according to the ratio $|V_{X}/\Delta E_{X}|$, and this is the most important limiting factor at the moment. As expected for the ground state, there is no sign of the many-body delocalization transition around $W\approx 3.5$.

The results from exact diagonalization indicate that for the ground state the variance of the number occupation is basically independent of system size. This allows us to estimate the error in this quantity for system sizes that cannot be solved exactly. One can appreciate that our results for small system sizes are in excellent agreement with exact results, and for large system sizes we have a small deviation towards lower values of the variance. For ground state calculations, our procedure can deal with system sizes much larger than exact diagonalization procedures and the quality of the results is quite acceptable.

The CPU time per sample for the present implementation of the algorithm up to fourth order transformations grows as $L^{10}$, implying a present practical limit at approximately $L=30$.

VI.2 Excited states

A crucial quantity when studying excited states is the ratio between the number of states generated and the total number of states, which we will call success ratio. We try to generate all states by systematically varying the occupation number of the different IOMs keeping the number of particles equal to half the number of sites. Nevertheless, one can expect that sometimes the algorithm starting from different states ends up in the same final state, since the evolving dynamic is not exact and some states may be more stable than others under this algorithm. We want to estimate the success ratio of our procedure. Once the level spacing is of the order of the energy uncertainty, it is difficult to calculate this quantity. We have designed a method to estimate the success ratio through the level statistics, as we explain in what follows.

We have studied the level spacing distribution for system sizes 10, 12 and 14, for which we can get all states by exact diagonalization and at the same time we can generate states corresponding to all possible combinations of IOM. We have calculated the level spacing distribution, normalized to the average spacing of the 20 nearest neighbors, for exact diagonalization results and for our method up to second and up to fourth order. The results are presented in Fig. 3. Panels a and b correspond to disorders $W=5$ and 2, respectively. The red curves correspond to exact results, while the blue curve correspond to our results up to second order with a cut-off $\lambda_{\rm min}=10^{-4}$ and with diagonalization of the projected Hamiltonian. The black curves in Fig. 3 are a fit of our results to the Poisson distribution. The exact results show Poisson statistics for $W=5$ and tend to Wigner-Dyson for $W=2$, in agreement with previous works Pal:2010gr and with the existence of an extended to localized transition for a critical disorder between 2 and 5. Our results show Poisson statistics for all disorders considered, as expected since we are cutting the number of density operators kept in the Hamiltonian and all terms are needed in order to reproduce the correct level repulsion. We observe a peak in our results at very small spacings which we interpret as due to different states collapsing in the same one. The area of this peak is a measure of the success rate. We take the success rate as one minus the area between our curve and the Poisson distribution. In Fig. 3 c we plot the success ratio as a function of disorder for our results up to second order for system sizes $L=10$ (red), 12 (black) and 14 (blue). Fourth order results are indistinguishable from second order results as long as the level spacing is concerned.

The succes ratio in Fig. 3 c is rather large, specially in the localized regime, where we only miss less than 4 %  of the states. It is weakly size dependent. To our knowledge, it is the larger success rate achieved in this type of calculations. Obtaining highly excited states of interacting systems is a very hard task and missing states often results in a biased selection of other states ratio .

Even for moderate size systems, it is not possible to ascertain the precision of a numerical method for excited states from the energy eigenvalues, since there the average spacing is smaller than the error in the energy. So, we estimate this precision from the temperature dependence of the variance of the number operator $\sigma_{L/2}^{2}$, whose quality depends on both the precision in the energy eigenvalues and in the expectation value of the operator. In Fig. 4 we show $\sigma_{L/2}^{2}$ as a function of temperature for four values of the disorder and a system size $L=12$. The continuous line corresponds to exact calculations and the dashed line to our procedure up to second order.

We note that our procedure is excellent for all $T$ in the localized regime, blue and green curves, while for lower values of the disorder is good at low $T$, but get worse at hight values of the temperature, where a delocalization transition is expected.

It is also interesting to know how our procedure behaves for excited states in relatively large systems. To this aim, we have computed $\sigma_{L/2}^{2}$ at infinite temperature, i.e., chosen excited states at random. In Fig. 2 we plot $\sigma_{L/2}^{2}$ at infinite $T$, upper set of curves, as a function of disorder for system sizes between $L=12$ and 30. The continuous curves are for results from exact diagonalization, for sizes $L=12$ and 16, while symbols, joint by dashed curves, correspond to our results up to second order. The exact results are almost independent of system size for large disorder, but they increases with size for small disorders, reflecting the many-body delocalization transition around $W_{\rm c}\approx 3.5$. Our results agree with exact results fairly well in the localized regime, but not so in the extended region where $\sigma_{L/2}^{2}$ does not increase with system size as it should. We think that this is due to the fact that excited states with very far apart resonances are difficult to select by our method (and probably by most methods) and there is a bias against these more extended states. If we consider transitions up to fourth order and include more excitations in our finite density approach, we should be able to reach this type of states. It is, in fact surprising that the simplest possible implementation of our procedure, that is, including transitions up to second order only, produces relatively good results in the localized regime.

VII Outlook

We have developed a new strategy to perform displacement transformations, by focusing on a specific reference state. Only transformations affecting expectation values for the reference state are included, leading to a higher efficiency. A final diagonalization of the subspace of states up to a fixed number of electron-hole excitations leads to high accuracy results.

In practice, the displacement transformation are chosen such that they minimize or maximize the energy expectation value with respect to a reference state. This applies equally well to ground states as to excited states. In the latter case, we also minimize the energy variance with respect to the reference state.

Renormalization group ideas along the lines developed for matrix product states in the so-called density matrix renormalization group can be naturally incorporated into our scheme White ; Thiery . One can start diagonalizing with displacement transformations a Hamiltonian for a small system of size $L$ and, at the same time, keep track of how the operators

[TABLE]

are transformed. These operators are included to increase the system size $L$ with an extra site to $L+1$. Once the original system is in diagonal form, we add a new site at $L+1$, that is, we add to the Hamiltonian the transformed operators of Eq. (28) plus the site energy $\epsilon_{L+1}n_{L+1}$. The procedure is repeated iteratively.

One can consider a large system right from the beginning, but perform only displacement transformations involving the first site. Once this site enters solely in diagonal terms, one can fix its occupation so that further transformations will not affect it. The same procedure is now applied to site two, and so on. Alternatively, one can start transforming the whole Hamiltonian and ‘freeze’ sites of either high or low energy along the diagonalization process. With ‘freezing’ we mean substituting the corresponding number operator $n_{i}$ by either 1 or 0.

We also need to credit an earlier approaches to construct many-body eigenstates beyond the gaussian approximation. NooijenNooijen2000 similarly developed a method to construct many-body states using transformations involving fourth order operators, though this method has not been implemented numerically.

Unlike for matrix product states, our procedure can be directly extended to any number of dimensions or even to amorphous systems with random links. In particular, it should be specially suited for molecular orbital calculations and for nuclear structure determination Szalay . The suitability of the method to deal with excited states should be relevant in the study of Rydberg atoms Rydberg ; Rydberg2 and in excited states in molecules Molecules ; Rydberg3 .

Displacement transformations can also be used to obtain even and odd normalized zero modes in random interacting Majorana models Monthus1 .

Finally, we want to remark that our method is a natural extension of the Hartree-Fock approximation for the construction of non-gaussian states. As an extension of Hartree-Fock theory specially adequate in the localized regime, this approach can be relevant in quantum chemical and nuclear matter calculations. We therefore expect that our method can lead to a new quantitive results in many fields.

Acknowledgements.

M.O. and A.M.S. were supported by Fundación Séneca grant 19907/GERM/15 and by Spanish MINECO Grant No. FIS2015-67844-R. L.R. was supported by the SNSF through an Ambizione Grant.

Appendix A Action of the displacement transformations

A general displacement transformation is defined through Eq. (5). A density operator of an state appearing in any creation or annhilation operator in $X$ transforms as

[TABLE]

the upper (lower) sign applies when $\delta$ corresponds to a creation (anhilation) operator in $X$. The operator $X$ itself transforms as

[TABLE]

We can choose $\lambda$ in order to cancel all the terms $X^{\dagger}+X$ appearing in the transformed Hamiltonian. The sum of the coefficients of all the classical terms with number operators in the quantum part of $X$ is equal to the energy difference $\Delta\epsilon_{X}$, Eq. (III) between the two configurations involved in the transition $X$. As $V_{X}$ is the coefficient of the term $X^{\dagger}+X$ in the Hamiltonian, condition Eq. (9) give us the strength of the displacement transformation that cancel all these terms.

It is straight forward to deduce from Eqs. (A) and (30) that the condition for canceling the quantum term $X^{\dagger}+X$ is equivalent to maximizing/minimizing the classical term $X^{\dagger}X-XX^{\dagger}$.

One can in principle diagonalize the interacting Hamiltonian by the application of successive displacement transformations, starting with those involving operators of the form $c_{j}^{\dagger}c_{i}$ and continuing with higher (larger number of terms) order operators. The transformation of quantum terms different from $X^{\dagger}+X$ generates other new terms, but they are always of higher order than the original term. To transform away terms of a given order, it is usually convenient to perform the displacement transformations that most reduce the energy first, although other schemes are also possible, for example, in the order of decreasing strength $V_{X}$. Although some times are generated terms with larger strengths, the general trend is a systematic decrease of the maximum remaining strength.