Obtaining highly excited eigenstates of the localized XX chain via DMRG-X

TL;DR

This paper benchmarks a DMRG-X algorithm variant for the localized XX chain, demonstrating high accuracy in eigenstate approximation and discussing its limitations related to entanglement and disorder.

Contribution

It introduces a small improvement to DMRG-X, assesses its performance on large systems using free fermion mapping, and analyzes entanglement effects on the algorithm.

Findings

DMRG-X accurately captures localized eigenstates for large systems

The algorithm exhibits a bias towards low-entanglement eigenstates at high entanglement

Increasing bond dimension improves convergence and reduces bias

Abstract

We benchmark a variant of the recently introduced DMRG-X algorithm against exact results for the localized random field XX chain. We find that the eigenstates obtained via DMRG-X exhibit a highly accurate l-bit description for system sizes much bigger than the direct, many body, exact diagonalization in the spin variables is able to access. We take advantage of the underlying free fermion description of the XX model to accurately test the strengths and limitations of this algorithm for large system sizes. We discuss the theoretical constraints on the performance of the algorithm from the entanglement properties of the eigenstates, and its actual performance at different values of disorder. A small but significant improvement to the algorithm is also presented, which helps significantly with convergence. We find that at high entanglement, DMRG-X shows a bias towards eigenstates with low entanglement, but can be improved with increased bond dimension. This result suggests that one must be careful when applying the algorithm for interacting many body localized spin models near a transition.

Equations

$$|\Psi ⟩=\sum _{j_1,\dots ,j_L}\sum _{0<{\gamma }_n<={\chi }_n}{B}_{{\gamma }_1}^{[1]j_1}{B}_{{\gamma }_1{\gamma }_2}^{[2]j_2}\dots B_{{\gamma }_{L-1}}^{[L]j_L}|j_1,\dots ,j_L⟩.$$

$$|\Psi ⟩$$

$$\equiv \sum _{\alpha =1}^{{\chi }_{n+1}}{\Lambda }_{\alpha \alpha }^{[n]}|{\alpha }_n{⟩}_L|{\alpha }_n{⟩}_R,$$

$$\mathcal{H}=-\sum _{i=1}^{L-1}\left(S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y\right)-\sum _{i=1}^L{h}_i{S}_i^z$$

$$S_i^{+}={(-1)}^{\sum _{j=1}^{i-1}{c}_j^{†}{c}_j}{c}_i^{†},S_i^{-}={(-1)}^{\sum _{j=1}^{i-1}{c}_j^{†}{c}_j}{c}_i,$$

$$\mathcal{H}$$

$$=$$

$$\mathcal{H}$$

$${\sigma }_E^2$$

$${\lambda }_n\sim \exp \left[-a{(\ln n)}^2\right]$$

$$P(S_E)$$

$$1-\sum _{i=1}^n{\left|⟨{\varphi }_i|\psi ⟩\right|}^2<ϵ$$

$$\mathcal{H}=\sum _{\alpha }{ϵ}_{\alpha }{a}_{\alpha }^{†}{a}_{\alpha }$$

$${\tau }_{\alpha }^{+}$$

$$=$$

$${\tau }_{\alpha }^{+}$$

$$\sum _{\beta =1}^{\alpha -1}{a}_{\beta }^{†}{a}_{\beta }$$

$${\tau }_{\alpha }^{+}$$

$$⟨\{{\sigma }_i^{\prime }\}|a_{\alpha }^{†}{a}_{\alpha }|\{{\sigma }_i\}⟩$$

$$=\sum _{\{s_i\}}{\delta }_{s_i,1}{\delta }_{s_{L+1},4}$$

$$M^{[i]}=\left[\begin{array}{cccc}\mathbb{I}& u_{\alpha }(i)S^{+}& u_{\alpha }^{\ast }(i)S^{-}& |u_{\alpha }(i){|}^2(S^z+1/2)\\ 0& -2S^z& 0& u_{\alpha }^{\ast }(i)S^{-}\\ 0& 0& -2S^z& u_{\alpha }(i)S^{+}\\ 0& 0& 0& \mathbb{I}\end{array}\right]$$

$${\mathcal{H}}_{\text{AA}}=-\sum _i\left(c_i^{†}{c}_{i+1}+\text{h.c.}\right)+W\sum _{i}V(i)c_i^{†}{c}_i$$

Full text

I Introduction

The density matrix renormalization group (DMRG) has become one of the most successful numerical tools in the study of ground state properties of one-dimensional systems White (1992); Schollwock (2011). DMRG is able to easily go to far larger sizes than accessible by methods such as exact diagonalization, and is even capable of working directly in the thermodynamic limit for translationally invariant systems McCulloch (2008). The success of DMRG has been greatly clarified by the advent of quantum information theoretic ideas such as the entanglement structure of quantum states and the language of matrix product states (MPSs) Schollwock (2011); Schuch et al. (2008). The basic premise of DMRG is the relatively small entanglement present in the ground states of most gapped local Hamiltonians. Specifically, such states obey an area law, that is, the bipartite entanglement entropy $S_{E}$ between two subsystems scales with the area of the boundary separating the regions Hastings (2007); I. Arad and Vazirani (2013). As a consequence the states can be efficiently represented with controllable error by MPSs which allow efficient computation of desired operator expectation values with controlled errors. This representability is coupled in practice with findability—which is the existence of efficient algorithms, such as the original DMRG algorithm, which find the desired representation in a time that scales polynomially —indeed linearly—with system size. While for the DMRG algorithm this scaling is an empirical fact about some unknown dynamical system111Although with possible exceptions, as discussed in Ref Eisert, there now also exists an algorithm that provably solves the problem in a time that scales polynomially in system size Landau Zeph et al. (2015). DMRG also works very well for critical points, where $S_{E}$ grows logarithmically with subsystem size which is modest enough to allow sufficient accuracy.

For thermalizing systems this cluster of ideas fails when addressed to highly excited eigenstates — eigenstates with a finite energy density corresponding to a nonzero temperature—due to the volume law entanglement that is generically present in these states. However in many body localized (MBL) systemsAnderson (1958); Basko et al. (2006); Pal and Huse (2010); Oganesyan and Huse (2007a); Nandkishore and Huse (2015); Altman and Vosk (2015) even highly excited eigenstates exhibit area law entanglementPal and Huse (2010); Bauer and Nayak (2013) which leads to the possibilityPekker and Clark (2017); Chandran et al. (2015) that excited states can be constructed efficiently via DMRG-like algorithms in a manner analogous to ground states. Being able to evade the size restrictions on many body exact diagonalization (MBED) in this fashion would be particularly useful in studying the eigenstate phase transitionHuse et al. (2013); Pekker et al. (2014); Parameswaran et al. (2016) from a volume-law obeying thermalizing phase at low disorder to an area-law obeying localized phase at strong disorder. While this MBL phase transition has been the subject of intense study numericallyOganesyan and Huse (2007b); Kjäll et al. (2014); Vosk et al. (2015); Potter et al. (2015); Devakul and Singh (2015a); Luitz et al. (2015a); Zhang et al. (2016a); Serbyn and Moore (2016); Zhang et al. (2016b); V. Khemani and Huse (2017), there is still much that is not well understood. Even for the “standard” model of MBLPal and Huse (2010), the random-field Heisenberg model, estimates of the critical disorder from numerical linked cluster expansion Devakul and Singh (2015b) predict a different value than finite size scaling estimates from exact diagonalization studies Pal and Huse (2010); Luitz et al. (2015b). Moreover, a recent studyKhemani et al. (2016) has shown that the entanglement properties of eigenstates jump discontinuously at the MBL-to-thermal transition—even while other properties look continuous—and that the critical point has far less entanglement than previously assumedGrover (2014), thereby lending encouragement for the potential use of DMRG-like methods all the way to the transition to the thermal phase.

Despite the good news on representability stemming from the area law, we note that the obstacles to constructing highly excited eigenstates of MBL systems are still formidable. One of these obstacles arises from the sheer smallness of gaps in the bulk of the many body spectrum which decrease exponentially with system size. Standard double precision arithmetic only allows a binary representation accurate to 52 bits which translates into an average spectral gap in a system of essentially the same number, say 50, spins. For systems larger than 50 spins we should expect to find exact spectral degeneracies which are artifacts and hence eigenvectors which are misleading superpositions of the underlying eigenstates. Indeed, this will start to happen even before size 50 for selected pairs of levels which are closer than average. The second obstacle is the process of optimization in a random system. Almost all work on ground state DMRG has been done on undisordered systems and based on general experience with disordered systems there is reason to worry that optimization will not work as well in the more complicated landscapes the latter exhibit. Potential obstacles here are the presence of Griffths regions and many-body resonances that may cause difficulty in finding an eigenstate.

With this caution injected, we note that the past year has nevertheless seen several studies that have tackled the problem of finding MPS representations for highly excited states. The first of these introduced an approach to finding approximations for all eigenstates of a MBL Hamiltonian at one goPollmann et al. (2016) and has been the object of more recent elaborationWahl et al. (2016). The second, which is the basis of present work, is a variant of DMRG aimed to target highly excited eigenstates, called DMRG-X Khemani et al. (2016). Rather than optimizing an MPS towards the ground state of the Hamiltonian, the DMRG-X algorithm iteratively optimizes the MPS towards an eigenstate of the Hamiltonian based on overlap with the state in the previous step. In a highly localized system, where eigenstates look very close to product states, DMRG-X was shown to converge very well starting from a product state as an initial state and reproduce the results of ED for accessible system sizes. In other works Yu et al. (2017); Lim and Sheng (2016); Serbyn et al. (2016); Kennes and Karrasch (2016), DMRG based techniques have been used to obtain excited eigenstates via energy (rather than overlap) targeting.

In this paper we return to the DMRG-X algorithm with a view to benchmarking its ability to accurately solve for eigenstates of large systems. We now clarify what we mean by accurately and what we mean by large. One can usefully think of a single DMRG step for a many body system as a many body exact diagonalization (MBED) with soft boundaries so it amplifies the size one can study with a diagonalization subroutine of given power; this is what enables ground states of systems with 100 spins to be studied on a laptop where only 14 (say) could be studied by MBED. Correspondingly we wish to understand the amplification achieved by DMRG-X in the middle of the spectrum over MBED. To gauge this meaningfully and even otherwise, we need to be able to assess the accuracy of the states obtained via DMRG-X. For any system we can quantify the accuracy by computing the variance of the energy in the states that we obtain and comparing it to the level spacing; this criterion was used in Ref. Khemani et al., 2016 at the limits of machine precision for the variance. However it is not known how to translate from a given accuracy of the variance to the accuracy of various quantities of more physical interest, which limits the utility of this metric.

To get around this limitation, in the present paper we study the the random field XX spin-1/2 chain which has an underlying free fermion representation. The free fermion character is not exploited by DMRG-X, which treats it as it does any interacting spin system. However the ability to exactly solve this model numerically for large system sizes in the fermionic version allows us to check the results of DMRG-X to larger system sizes than are accessible via MBED. Specifically, we make use a defining characteristic of localized phases, namely the existence of an emergent set of an extensive number of commuting $\mathbb{Z}_{2}$- valued local integrals of motionHuse et al. (2014); Serbyn et al. (2013) (often called “l-bits”), and use the the free fermion representation to construct the exact l-bits for very large system sizes. We then compute the expectation values of these l-bits in the DMRG-X obtained states, which allows us to directly assess how well these states approximate the exact eigenstates using a more physical (and economical) measure. A further advantage of studying the XX chain is that we are also able to study the representability of the eigenstates to large system sizes via computation of their entanglement entropies. Overall we find that DMRG-X achieves a roughly three-fold amplification over MBED when we require that all l-bits are accurate to at least 90 percent; most are vastly more accurate, but, however, DMRG-X is biased towards less entangled states when entanglement becomes too high to represent. While this particular demonstration is specific to the XX chain, we believe this analysis offers much encouragement that we can trust DMRG-X results even for interacting systems (away from the MBL phase transition). We note that similar optimism comes from the work of Serbyn et alSerbyn et al. (2016) who examine the structure of the entanglement spectrum for numerically obtained eigenstates.

In the balance of the paper, we start with a quick recapitulation of the DMRG method in the language of MPSs (Sec. II). Following this, we review the basics of the random XX chain and the accuracy measure (Sec. III), discuss representability where we also address the issue of rare Griffiths effects (Sec. IV), turn to findability where we find it important to modify the DMRG-X to use a hybrid metric for the updates (Sec. V), show that the resulting algorithm is accurate enough to go beyond ED (Sec. VI) and end with some concluding remarks (Sec. VII). Appendices describe some technical details on the locality (App. A) and representation (App. B) of the l-bit operators for the XX chain, and some unexpected wrinkles that crop up in studying the closely related problem of localization in quasiperiodic potentials (App. C).

II Brief Review of MPS/MPO/DMRG formalism

In this section, we briefly recapitulate the standard DMRG algorithmWhite (1992) implemented in the language of matrix-product states (MPSs) Schollwock (2011). This section closely follows the supplementary material in Ref. Khemani et al., 2016. A general quantum state $|\Psi\rangle$ for a one-dimensional system of $L$ sites can be written in the following matrix-product state (MPS) form:

[TABLE]

where $|j_{n}\rangle$ with $j_{n}=1,\dots,d$ is a basis of local states at site $n$ (for a spin 1/2 system, $d=2$ and $|j_{n}\rangle=|\uparrow\rangle,|\downarrow\rangle$), and the $B^{[n]}$ are rank three tensors (except on the first and last sites where they are rank two tensors). Figure 1(a) shows a pictorial representation of an MPS. The enternal legs $j_{n}$ are the “physical” spin indices whereas the internal legs $\gamma_{n}$ are the virtual indices that are contracted. Each $B^{[n]j_{n}}$ is a $\chi_{n}\times\chi_{n+1}$ matrix (at the boundaries $\chi_{1}=\chi_{L+1}=1$) and each matrix product $\prod_{i}B^{[i]j_{i}}$ in Eq. (1) produces a complex number which is the amplitude of $|\Psi\rangle$ on the basis state $|j_{1}\cdots j_{L}\rangle$.

The maximum dimension $\chi$ of the $\{B^{[n]}\}$ matrices is called the bond-dimension of the MPS and low entanglement states can be efficiently represented my MPSs of bond dimension $\chi\ll d^{L/2}$. The relationship between $\chi$ and the entanglement can be made more precise by considering the Schmidt decomposition of the state $|\Psi\rangle$. For a given bipartition of the system into left and right halves, a singular value decomposition can be used to rewrite $|\Psi\rangle=\sum_{\alpha}\Lambda_{\alpha}|\alpha\rangle_{L}|\alpha\rangle_{R}$ where the $|\alpha\rangle_{L/R}$ form orthonormal bases for the left and right halves respectively, and the entanglement entropy of the bipartition is defined through the Schmidt values $\Lambda_{\alpha}$ as $S_{E}=-\sum_{\alpha}|\Lambda_{\alpha}|^{2}\ln|\Lambda_{\alpha}|^{2}$. Following a prescription by VidalVidal (2003), it is possible to define a canonical form (Fig.1(b)) for the MPS by rewriting each matrix $B^{[n]j_{n}}$ as a product of a $\chi_{n}\times\chi_{n+1}$ dimensional complex matrix $\Gamma^{[n]j_{n}}$ and a square diagonal matrix $\Lambda^{[n]}$ such that matrices $\Lambda^{[n]}$ matrices contain the non-zero Schmidt values for a bipartition between sites $n$ and $n+1$

[TABLE]

and the states $|\alpha_{n}\rangle_{L},|\alpha_{n}\rangle_{R}$ define the orthonormal Schmidt states for the left and right halves of the bipartition respectively. This canonical form clearly relates the bond dimension $\chi$ to the number of Schmidt values contributing significantly to the entanglement entropy.

Starting from an initial random MPS, the DMRG algorithm iteratively finds the ground state $|\psi_{0}\rangle$ by sweeping through the system and variationally optimizing the MPS matrices $B^{[n]j_{n}}$ on neighboring sites to locally minimize the energy $\langle\psi_{0}|\mathcal{H}|\psi_{0}\rangle$ (keeping the rest of the chain fixed). In the commonly used two-site update which simultaneously updates two sets of matrices $B^{[n]j_{n}}$ and $B^{[{n+1}]j_{n+1}}$, an effective Hamiltonian $\mathcal{H}_{\rm eff}$ is constructed by projecting $\mathcal{H}$ to a mixed $\chi_{n}\chi_{n+2}d^{2}$ dimensional basis. Here, the local basis states $|j_{n}\rangle|j_{n+1}\rangle$ represent the two updated sites, and the eigenstates of the reduced density matrix $|\chi_{n}\rangle_{L}|\chi_{n+2}\rangle_{R}$ compactly represent the environment to the left and right of the updated sites. The ground state of $\mathcal{H}_{\rm eff}$ is then found — which is the optimal state for minimizing $\langle\psi_{0}|H|\psi_{0}\rangle$ in this subspace —and the matrices on sites $n,n+1$ are updated. The procedure is repeated for all sites until convergence is achieved. The matrix-product operator (MPO) representation of $\mathcal{H}$, defined exactly analogously to Eq. 1 but now using 4-index tensors $M$, is shown in Fig. 1(c), and the effective Hamiltonian is depicted pictorially in Fig. 1(d).

The only difference between the ground-state DMRG algorithm outlined above and the DMRG-X algorithm of Ref. Khemani et al., 2016 is in the update step. In the DMRG-X algorithm, all the $d^{2}\chi^{2}$ eigenstates of $\mathcal{H}_{\text{eff}}$ are obtained instead of just its ground state, and the matrices $B^{[n]j_{n}}$ are updated using the eigenstate of $\mathcal{H}_{\rm eff}$ with the maximum overlap with the previously found state in the iterative scheme. As will be discussed in Sec. IV, we will present a modification to this algorithm at this step, where instead of simply choosing the eigenstate with maximal overlap, we also minimize the energy variance within a small subspace of states with high overlap. The algorithm is initialized with an appropriate initial state which is perturbatively “close” to the true eigenstates of the MBL Hamiltonian, such as a product state.

III The XX Chain

The model we study using DMRG-X is the spin-1/2 XX chain with random fields,

[TABLE]

where $S^{x,y,z}_{i}$ are spin-1/2 operators on site $i$ and $h_{i}\in[-W,W]$ is chosen from a uniform distribution with width $2W$. This model is Anderson localized for any nonzero $W$, and the localization length scales as $\xi\approx 25/W^{2}$ at low $W$ Kappus M. and Wegner F. (1981), while at large $W$ the locator expansion yields $\xi\sim 1/\log W$. This model can be mapped via a Jordan-Wigner transformation onto a system of non-interacting fermions and solved exactly to very large system sizes.

Using the standard Jordan-Wigner substitutions for $S^{\pm}_{i}=S^{x}_{i}\pm iS^{y}_{i}$,

[TABLE]

we arrive at a quadratic Hamiltonian in terms of the fermionic $c^{\dagger}_{i},c_{i}$ operators

[TABLE]

up to a constant. This single particle Hamiltonian $H_{ij}$ is only of dimension $L$, and can be diagonalized numerically, $H_{ij}=\sum_{\alpha}U_{i\alpha}\varepsilon_{\alpha}U^{\dagger}_{\alpha j}$, giving

[TABLE]

up to additive constants, where $a_{\alpha}^{\dagger}=\sum_{i}c^{\dagger}_{i}U_{i\alpha}$, and similarly $a_{\alpha}=\sum_{i}U^{\dagger}_{\alpha j}c_{j}$.

In spin-language, we refer to the $\sigma^{z}_{i}=2c^{\dagger}_{i}c_{i}-1$ as physical-bit (p-bit) operators and $\tau^{z}_{\alpha}=2a^{\dagger}_{\alpha}a_{\alpha}-1$ as localized-bit (l-bit) operators, in accordance with standard MBL nomenclature Huse et al. (2014). This non-interacting Hamiltonian then takes the simple form $\mathcal{H}=\sum_{\alpha}h_{\alpha}\tau^{z}_{\alpha}$, where $h_{\alpha}=\epsilon_{\alpha}/2$. The $2^{L}$ eigenstates of $\mathcal{H}$ are then obtained by picking each l-bit to be $+1$ or $-1$, corresponding to a filled ($n_{\alpha}=1$) or empty ($n_{\alpha}=0$) fermionic state respectively, and the many-body fermionic eigenstates are constructed as $\left|\left\{n_{\alpha}\right\}\right\rangle=\prod_{\alpha}a^{\dagger}_{\alpha}n_{\alpha}\left|0\right\rangle$.

It is interesting to note that the fermionic raising and lowering operators $a^{\dagger}_{\alpha}$ and $a_{\alpha}$ are inherently nonlocal in spin language due to the ${(-1)}^{\sum_{i}c^{\dagger}_{i}c_{i}}$ chain from the Jordan-Wigner transformation. However, in MBL systems one tends to speak of local bosonic raising and lowering operators $\tau^{\pm}_{\alpha}$. Bosonic raising and lowering operators can be constructed via a reverse Jordan-Wigner transormation from the $a^{\dagger}_{\alpha}$,$a_{\alpha}$ operators, but it is not a priory clear that these operators are indeed local.222They are obviously local if one takes the definition of an MBL system to be the existence of a finite depth unitary transformation that diagonalizes the Hamiltonian. In App. A we show that, indeed, these bosonic raising and lowering operators can be explicitly shown to be local given that the original l-bit operators are local.

Our primary interest in the l-bit operators is that they can be used to gauge the accuracy of a DMRG state. An exact eigenstate should have $|\langle\tau_{\alpha}^{z}\rangle|=1$ (in the absence of degeneracy). Deviations from this tell us exactly where DMRG has failed to capture the state, and how badly it has done so. The l-bit operators themselves can be efficiently expressed as matrix product operators (MPOs) using an internal bond dimension of only 4 (App. B).

The other measure of accuracy is the total energy variance of the state, $\sigma_{E}^{2}=|\langle H^{2}\rangle-{\langle H\rangle}^{2}|$, which is also efficient to calculate for an MPS. This should be closer to zero the closer a state is to an eigenstate of the Hamiltonian. Expressed in terms of the l-bit operators, $\sigma_{E}^{2}$ contains expectation values and correlators of these operators,

[TABLE]

and is thus a related but complementary measure of accuracy.

IV Representability: $S_{E}$ and Griffiths Effects

We now turn to the question of how efficiently the MB eigenstates of the XX model can be represented, which entails a study of the entanglement properties of these states. We find that rare Griffiths-like effects cause the typical and worst-case entanglement entropy across the bonds of the system to scale differently with $L$, which in turn affects the maximum bond-dimension needed to represent these states and the scaling of the computational time of the DMRG-X algorithm. While the algorithm always scales polynomially with $L$, in practice, the power of the polynomial can get quite large as disorder is lowered.

The entanglement spectrum (ES) of a state across a cut determines how well that state can be represented as an MPS with a finite bond dimension $\chi$ on that cut. The ES in the MBL phase has been recently shown to decay as a power law Serbyn et al. (2016). However, the arguments of Ref Serbyn et al., 2016 only apply to interacting MBL systems.333The argument in Ref Serbyn et al. (2016) relies on the l-bit operators flipping many spins within a radius $r$, but in a noninteracting model the l-bit operators only flip two. In the XX model, the spectrum of the single particle entanglement Hamiltonian $\mathcal{H}_{E}$ is found to have a constant density of states on average, similar to the clean system Eisler (2009), which will lead to an entanglement spectrum decaying (on average) as

[TABLE]

for some constant $a$. This is one factor that differentiates our non-interacting model from an interacting model, making it easier to treat with DMRG-X. For now, we ignore the details of the estanglement spectrum and assume that a bond dimension of $\chi\sim e^{S_{E}}$ is needed to represent a state with entanglement entropy $S_{E}$ (which is equivalent to assuming a flat ES).

The entanglement entropy $S_{E}$ can be obtained exactly for a given many-body free-fermion eigenstate $|\psi\rangle=\left|\left\{n_{\alpha}\right\}\right\rangle$ using the correlation matrix method Eisler (2009). The idea is that the reduced density matrix for a partition $A$ is represented as a thermal density matrix for an “entanglement Hamiltonian” $\mathcal{H}_{E}$, $\rho_{\text{A}}=\text{Tr}_{\bar{A}}\rho=e^{-\mathcal{H}_{E}}$, where $\mathcal{H}_{E}$ itself is a non-interacting Hamiltonian whose single particle energy eigenvalues can be found from the matrix of two-point correlation functions of fermionic operators evaluated in the state $|\psi\rangle$. The entanglement entropy $S_{E}=-\text{Tr}\rho_{A}\ln\rho_{A}$ can then be calculated knowing the eigenvalues of $\rho_{A}$. This knowledge can also be used to obtain the (many-body) eigenvalues of $\rho_{\text{A}}$, $\lambda_{n}$, which constitute the entanglement spectrum.444Finding the largest $k$ eigenvalues of $\rho_{\text{A}}$ involves finding the $k$ lowest total eigenenergies given all the single particle energies, which is a standard computational problem.

Figure 2(left) shows the mean values of the mid-chain entanglement entropy and the maximum entanglement entropy within the whole chain. The maximum entanglement entropy $S_{E}^{\text{max}}$ for a given chain length is important as it will determine the maximum bond dimension $\chi$ needed to accurately capture that state as an MPS. While the mid-chain entanglement entropy saturates to an area law, the maximum entanglement entropy continues to grow with $L$.

The growth of $S_{E}^{\text{max}}$ with $L$ can be understood as a type of Griffiths effect. While there is no true delocalization in this model, a finite patch of length $\ell$ appears delocalized if the localization length within that patch exceeds $\ell$: $\xi\geq\ell$. Thus, the finite patch appears to have a critical disorder $W_{c}(\ell)=5/\sqrt{\ell}$ (this follows from the expression for the localization length $\xi=25/W^{2}$). The probability to get such a patch of length $\ell$ with every onsite field lying within $W_{c}(\ell)$ is $P(\ell)={(W_{c}(\ell)/W)}^{\ell}$. Typically, a system of size $L$ will have $N(\ell)=LP(\ell)$ such patches, and equating $N(\ell)\approx 1$ leads to $\ell\sim\log L/\log\log L$ at leading order. This apparently delocalized patch will have volume law entanglement $S_{E}^{\text{max}}\sim\ell$, and therefore $S_{E}^{\text{max}}\sim\log L/\log\log L$, which explains the mostly logarithmic scaling with the slight downwards curvature in the $S_{E}^{\rm max}$ plot in Figure 2(left). In particular, this analysis shows that the maximum bond dimension needed to represent a particular state which scales as $\log\chi^{\text{max}}\sim S_{E}^{\text{max}}$ will scale with $L$ (as has been previously observed Bauer and Nayak (2013); Khemani et al. (2016)), as opposed to approaching a constant as in the case of ground states of gapped systems. In App. C, we discuss how a similar effect may arise even in some quasiperiodic systems, where there are no Griffiths effects.

This analysis also helps us deduce the scaling of the time needed to carry out one sweep of the DMRG-X. If we run the algorithm with variable bond dimension $\chi$ on each bond, the locus of the maximum computational time is determined by the tail of the distribution of entanglement entropies, with the shape of the tail determined by the same Griffiths-like effects that were just discussed. If we take $S_{E}=\kappa\ell$ to be the entanglement of a thermal region of length $\ell$, assuming that the entanglement on each bond is independent, we have that

[TABLE]

with $W_{c}(\ell)$ as defined earlier. We therefore see that the tail of $P(S_{E})$ will decay slightly faster than exponential, and decay faster for higher $W$. However, the computational time required to deal with a cut of entanglement $S_{E}$ scales exponentially with $S_{E}$: taking $\chi=e^{S_{E}}$, we know that the computational cost of diagonalizing the effective Hamiltonian in a DMRG step (a $d^{2}\chi^{2}\times d^{2}\chi^{2}$ matrix) scales as $(\chi^{2})^{3}\sim e^{6S_{E}}$. Thus, the mean computational time will be dominated by bonds with entanglement where $e^{6S_{E}}P(S_{E})$ is maximized. For low $W$, this implies that the computational time will be dominated by the largest entanglement across any bond in the sample $S_{E}^{\text{max}}$. Thus, the growth of $S^{\text{max}}_{E}$ with $L$ leads to a complexity growing (almost) polynomially as $L^{b/\log\log L}$ with some constant $b$. For higher $W$, as the tail decays more quickly, the dominating time shifts to an entanglement that appears more typically in the sample. In this regime, the computational time appears to scale linearly with $L$. In the case of a purely exponentially decaying tail in $P(S_{E})$, the computation time experiences a Hagedorn-like transition at some $W$ where the dominating computational cost switches from being in the tail of the distribution to some finite value near the peak of the distribution.

Thus far we have assumed a flat ES in order to make the approximation $\chi=e^{S_{E}}$. In reality, our ES decays on average quite quickly according to Eq 10, and the relationship between $\chi$ and $S_{E}$ depends on the exact structure of the ES. However, $\chi$ should still scale exponentially with $S_{E}$ on average. While the exact values of the coefficients of the exponential will depend on the decay of the ES and the error threshold, the maximum bond dimension $\chi^{\text{max}}$ will still scale as predicted. We confirm the dependence of $\chi$ on $S_{E}$ by sampling the entanglement spectra obtained from the correlation matrix Eisler (2009). The largest values in the many-body entanglement spectra are obtained from the eigenvalues of the single particle entanglement Hamiltonian $\mathcal{H}_{E}$ via an iterative algorithm. We find the largest $\chi$ eigenvalues such that the total sum of the discarded spectrum is less than $10^{-5}$. Figure 2(right) shows a histogram of $\chi$ vs. $S_{E}$ at $W=1.0$, clearly demonstrating a rough exponential relation between the two as expected. Notice that very large $\chi$, much larger than reachable, is needed to attain this level of accuracy for low $W$. Figure 2(center) shows the max and mean $\chi$ as a function of $L$, demonstrating the expected scaling behaviour. At high $W$, $\chi$ becomes very small and the relationship between $S_{E}$ and $\chi$ is less apparent.

V Findability: Hybrid DMRG-X Algorithm

Having seen that the localized eigenstates of the XX model can be represented reasonably efficiently by MPS, we turn to the issue of finding such representations using the DMRG-X method. In performing DMRG-X for this relatively simple model, we found that the algorithm was having great difficulty converging for certain disorder realizations within a particular range of parameters. The problem is not in the ability to represent such a state as an MPS, as better MPS representations can be constructed from the exact eigenstate, using an algorithm by S.R. White Fishman and White (2015). Note that this algorithm is not the only way and may not be the optimal method to construct such MPSs, and in fact there has been much work on MPS representations for known states Silvi et al. (2013); Katsura and Maruyama (2010); Murg et al. (2012). White’s algorithm iteratively applies pairwise rotations to an initial product state MPS, truncating the bond dimension as necessary, to construct a variational approximation—not necessarily the best one—to the exact eigenstate. Given such an MPS as an initial state, performing a DMRG-X sweep will actually worsen the state when in the problematic parameter range ($W\approx 1.5$). The typical energy variance $\sigma_{E}^{2}$ of the DMRG-X obtained eigenstates is shown for $L=32$ as a function of $W$ and $\chi$ in Figure 3(left). There is a particular region $1<W<2$ in which variance remains high and the algorithm fails to converge onto an eigenstate A significant portion of samples within this region did not show a monotone decrease in variance with increasing $\chi$, meaning that the optimal state was never being reached.

A problem occurs when there is a near degeneracy between two many-body eigenstates that are not too many spin flips away. If the dimension of our restricted subspace is not sufficient to accurately capture the Hamiltonian acting on these two states, they may appear to be much closer in energy than they actually are, resulting in a false resonance in $\mathcal{H}_{\text{eff}}$. This happens as the subspace is optimized to accurately represent only the current state, not necessarily any others and therefore not the full Hamiltonian. The result of this is that some of the eigenstates of $\mathcal{H}_{\text{eff}}$ are superpositions of the true eigenstate with other nearby energy states. Thus, picking one of them based on maximal overlap does not move the overall state towards the correct answer. This is a problem in the findability.

Here, we present a small modification that can correct for this type of error. Rather than simply picking the eigenstate with highest overlap, we look at the $n$ largest overlap eigenstates. The variance can then be optimized within this $n$-dimensional subspace, via a gradient method starting with the current state as a seed. Thus, variance and overlap are both taken into account in choosing the next state, which is closer to the true eigenstate than any eigenstate of $\mathcal{H}_{\text{eff}}$.

The variance calculation can be done efficiently in MPS language by keeping track of $\mathcal{H}^{2}_{\text{eff}}$ in addition to just $\mathcal{H}_{\text{eff}}$. Note that $\mathcal{H}^{2}_{\text{eff}}$ here is $\mathcal{H}^{2}$ in the effective subspace, not simply $\mathcal{H}_{\text{eff}}$ squared. Then, the variance $\sigma_{E}^{2}=\left|\langle\mathcal{H}_{\text{eff}}^{2}\rangle-{\langle\mathcal{H}_{\text{eff}}\rangle}^{2}\right|$ can be optimized using a gradient based optimization algorithm such as the BFGS algorithm. These matrices can be further reduced to $n$ by $n$ matrices by only looking in the subspace of the $n$ states with highest overlap with the previous state. An important implementation detail is that we optimize $\log\sigma^{2}_{E}$ rather than simply $\sigma^{2}_{E}$, as we are concerned with very small values of the variance.

This hybrid algorithm works best when the MPS is already close to an eigenstate, as the overlap then provides a useful measure of which states to optimize over. It is therefore often helpful to run a few sweeps with $n=1$ (the original algorithm) before increasing $n$. Rather than keeping $n$ fixed, it is also possible to vary $n$. If we let $\left|\phi_{i}\right\rangle$ be the eigenstates of the effective Hamiltonian ordered from highest to lowest overlap with the previous state $\left|\psi\right\rangle$, then we can choose $n$ to be the smallest integer satisfying

[TABLE]

where $\epsilon$ is chosen to be a small number ($\epsilon=10^{-10}$ is used in our calculations).

Figure 3(center) shows that the application of this hybrid algorithm greatly improves the variance near the troublesome region. The peak in the variance of the original algorithm arose due to convergence issues, which have been clearly resolved with the hybrid algorithm. Also shown are the variances obtained from the MPS constructed using White’s algorithm Fishman and White (2015).555Note that it is not that the algorithm of Ref Fishman and White (2015) is failing at low $W$. Similar curves were obtained for smaller $L$ accessible by ED, and a standard MPS compression of the ED state yielded similar results. Those variances do what one might expect from a direct truncation in the entanglement spectrum, in that they perform poorly at low $W$ where there is high entanglement and monotonically do worse.

Surprisingly DMRG-X appears to be performing well even at low $W$, where it is expected to fail. Exploring further, Figure 3 shows the mean entanglement entropy of the DMRG-X states in comparison the mean of the exact eigenstates. We see that the state to which DMRG-X converges towards is significantly lower in entanglement than the true exact eigenstates. The reason for this is that DMRG-X is biased towards states which it can represent well with limited $\chi$. That is, when the available $\chi$ is small, the algorithm is biased towards states with low entanglement entropy, which lie away from the center of the spectrum in energy. At very low $W$, the only representable states for small $\chi$ are near the edges of the spectrum, which DMRG-X finds.

The primary reason this occurs is that the algorithm allows the state to “drift” with each iteraction. This is in contrast with, say, the naive energy-targeting algorithm in which the eigenstate of the effective Hamiltonian closest to a specified energy is picked at each step (which fails to converge due to small many-body level spacings Khemani et al. (2016)). The closely related algorithm ES-DMRG Yu et al. (2017), which works by selecting the eigenstate with the energy closest to the energy of the previous state, is also subject to this “drift” towards low entanglement states at low $W$.

In Fig 4 we show hybrid DMRG-X calculations, where an approximated eigenstate MPS (constructed as described earlier) has been given as seed. In contrast to previous calculations where the bond dimension $\chi$ is ramped up slowly, in these calculations each choice of $\chi$ represents an independent calculation. As seen in Fig 4(left), at high $W$, DMRG-X (both hybrid and original) is able to converge upon the chosen eigenstate very quickly. However, as $W$ is lowered and the entanglement in these eigenstates increases, neither are able to converge to the desired target. Instead, a compromise is made and the algorithm converges upon a different eigenstate with lower entanglement away from the center of the spectrum. Increasing the bond dimension $\chi$ helps extend this region of convergence. As seen in Fig 4(right), the hybrid algorithm suppresses this drift towards extremal energies somewhat, but is still subject to it.

VI Accuracy of DMRG-X

We now turn to our goal, which is to evaluate the performance of DMRG-X on system sizes larger than accessible by many-body ED. To this end we use a modest bond dimension of $\chi=32$, corresponding to a $2^{2}\chi^{2}=2^{12}$ dimensional effective Hamiltonian in the worst case. Thus the computational effort at the diagonalization step of DMRG-X is equivalent to the exact diagonalization of the full Hamiltonian of an $L=12$ system. With this restriction we examine the accuracy of the algorithm on system sizes much bigger than $L=12$ via the energy variance and, more revealingly, l-bit accuracy.

We begin by studying $L=32$. To better understand the quality of the state, Fig 5(left) shows the typical deviation of $\langle\tau^{z}_{\alpha}\rangle$ from $\pm 1$, defined by $\delta=1-\left|\langle\tau^{z}_{\alpha}\rangle\right|$. At $\chi=32$, the typical l-bit expectation value is accurate up to roughly $10^{-6}$ at its worst. Figure 5(right) shows the average number of badly captured l-bits in a sample, which we define to be those with $\delta>0.1$. Indeed, with high enough bond dimension DMRG-X is able to obtain most eigenstates with no bad l-bits, with the mean number of bad l-bits falling below 1. This is our central demonstration in this paper—that DMRG-X can obtain all of the l-bits correctly for system sizes bigger than those treatable by the ED routine utilized by the code. Evidently, if we were to use the actual maximum size treatable by ED ($L>12$) we would expect to get to even bigger actual system sizes ($L>32$) and still expect to get all of the l-bits correctly by our metric.

We now turn to the scaling with $L$ while keeping $\chi=32$ fixed. If the MPS exhibited a constant error density throughout, one would expect the energy variance to scale linearly with $L$. Figure 6(left) instead shows the variance scaling with $L$, which appears to be faster than linear. The typical l-bit errors, shown in Figure 6(right), are also very small and slowly increasing with $L$. The reason behind the super-linear growth in variance is that our states are in fact not of constant error density, as we are simply truncating each bond at $\chi$ singular values. To construct a state of constant error density, one would have to consider a variable $\chi$ on each bond such that the total discarded Schmidt weight is below an error threshold. In such a scheme one would still need to define a maximum allowable $\chi$, $\chi_{\text{max}}$, such that the computation completes within a reasonable timeframe. One can then either discard samples that would exceed this limit (which avoids the issue of potential Griffiths regions) or truncate the bonds at $\chi_{\text{max}}$ (which no longer produces constant error density states). Our fixed-$\chi$ choice corresponds to the latter with a zero error threshold.

VII Concluding Remarks

We have tested the performance of the DMRG-X algorithm for finding highly excited eigenstates of the disordered XX model. With some small improvements to the algorithm, we were able to go significantly beyond the range of system sizes accessible to exact diagonalization in the many-body Hilbert space. Most of the eigenstates are able to be obtained to high accuracy with only $\chi=32$, using the l-bit expectation value as a metric. Thus, we have successfully shown that DMRG-X is able to go well beyond exact diagonalization in the well localized regime, converging to eigenstates quickly with increasing $\chi$. For more delocalized states, DMRG-X makes a compromise to improve accuracy and convergence by biasing itself towards lower entanglement states.

There are still many ways of improving the capabilities of the algorithm. Rather than exact diagonalization of the effective Hamiltonian, shift-invert Lanczos can be used to obtain all the eigenstates within an energy window. The states with high overlap with the original state typically have very similar energy, so Lanczos is also able to find a relevant subspace of $n$ states over which to optimize variance. This theoretically allows a great increase in the range of accessible $\chi$.

The natural next step is to consider application of DMRG-X to interacting systems which exhibit MBL. Near the transition, Griffiths regions with high entanglement play an important role, thus one must be careful when applying DMRG-X near the transition. This points out the need for larger bond dimensions which may be obtained as previously mentioned.

Acknowledgements.

This research was supported by the Harvard Society of Fellows (VK), the DFG (Deutsche Forschungsgemeinschaft) Research Unit FOR 1807 through grants numbers PO 1370/2-1 (FP), the NSF-DMR via grant number 1311781 (SLS) and the Alexander von Humboldt Foundation via a Humboldt Award (SLS).

Appendix A Locality of l-bit operators

We have a Jordan Wigner transformed spin chain, diagonalized by a set of fermionic operators

[TABLE]

where each $a^{\dagger}_{\alpha}=\sum_{i}u_{\alpha}(i)c^{\dagger}_{i}$, and in the spin language, $c^{\dagger}_{i}=\prod_{j<i}{(-1)}^{(S^{z}_{j}+\frac{1}{2})}S^{+}_{i}$, and similarly for the conjugates. We want to show that the bosonic raising and lowering operators $\tau^{+}_{\alpha}$ and $\tau^{-}_{\alpha}$ are also localized, given that the fermionic ones $a^{\dagger}_{\alpha}$ and $a_{\alpha}$ are in terms of the physical $c^{\dagger}_{i}$,$c_{i}$ operators.

The corresponding bosonic raising operator for $\tau^{z}_{\alpha}=2a^{\dagger}_{\alpha}a_{\alpha}-1$ is then defined by an inverse Jordan Wigner transformation

[TABLE]

For simplicity, we will assume that $u_{\alpha}(i)$ is strictly localized within a range $\xi$, and is therefore only nonzero if $|\alpha-i|\leq\xi$ (in reality, there will be an exponential decay with lengthscale $\xi$, and our argument can be generalized for some error threshold). We can then write $\tau^{+}_{\alpha}$ as

[TABLE]

Expressing the first exponent in Equation 16 as

[TABLE]

we see that if $j<\alpha-\xi$ or $k<\alpha-\xi$, then the sum over $\beta$ contains all non-zero elements of $u_{\beta}(i)u^{*}_{\beta}(j)$, resulting in a dirac delta $\delta_{i,j}$ by completeness. This reduces the expression to

[TABLE]

where we have restricted the summation indices $j,k\leq\alpha+\xi-1$, since $\beta\leq\alpha-1$ and we have assumed $u_{\beta}(i)$ is zero for $i>\beta+\xi$.

The two sums in Eq 18 commute, and therefore the $c^{\dagger}_{j}c_{j}$ sum can be canceled out with the chain in Eq 16, leaving

[TABLE]

which is strictly localized within $\alpha\pm\xi$.

Appendix B MPO Representation of l-bit operators

The $a^{\dagger}_{\alpha}a_{\alpha}$ operators can be constructed as an MPO by hand. The methodology is similar to that used to construct the MPO representation of the Hamiltonian Schollwock (2011). We use the internal dimensions $\{s_{i}\}$ to keep track of the different terms in the expansion of $a^{\dagger}_{\alpha}a_{\alpha}$. We can define the full operator as an MPO, acting on spin in the $S^{z}$ basis,

[TABLE]

where each $s_{i}$ goes from $1$ to $4$ and $\delta_{i,j}$ is the Dirac delta.

The different terms in $a^{\dagger}_{\alpha}a_{\alpha}=\sum_{i,j}u_{\alpha}(i)u^{*}_{\alpha}(j)c^{\dagger}_{i}c_{j}$ can be taken into account by defining $M^{[i]}_{s_{i},s_{i+1}}$ appropriately:

[TABLE]

Thus, the operator $a^{\dagger}_{\alpha}a_{\alpha}$ has been expressed as an MPO of internal dimension $4$, and the expectation value in an MPS can be efficiently calculated. One can then simply take $\tau^{z}_{\alpha}=2a_{\alpha}^{\dagger}a_{\alpha}-1$.

Appendix C Griffiths Effects and Quasiperiodicity

Many of the difficulties in studying disordered systems using DMRG arise due to the existence of locally clean-looking patches, which are what the tail of the distribution $P(S_{E})$ consist of. An exponentially decaying tail in $P(S_{E})$ will lead to the typical maximum entanglement in a system of size $L$ growing as $S_{E}^{\text{max}}\sim\log L$, resulting in the computational difficulty of DMRG grow faster than linearly. To avoid the difficulties of these Griffiths region in studies of localization, one possibility is to examine quasiperiodic systems.

One commonly studied possibility is the Aubry-André model, which is given by the Hamiltonian

[TABLE]

where $V(i)=\cos(2\pi\omega(i+\phi))$, $\omega$ is an irrational number which we choose to be 1 over the golden ratio, and $\phi$ is some offset. The exact delocalization-localization critical point of this model is known to be at $W=2$, and the critical wavefunctions can be found exactly Aubry and Andre (1980). The quasiperiodic potential does not allow for large Griffiths-like regions to exist.

Despite this, sampling chains of length $L$ (over different phase offsets), there is still a clear logarithmic growth of $S^{\text{max}}_{E}$ with $L$, corresponding to an exponentially decaying tail in $P(S_{E})$, which is reminiscent of the existence of Griffiths regions!

This effect can be understood as coming from the existence of special reflection symmetry points in the quasiperiodic potential. For example, if $\phi=0$, all the sites on the left side of $i=0$ will be in exact resonance with another site on the right side, $V(-i)=V(i)$. Each resonant pair will then form an odd-even superposition and the many-body entanglement entropy across that cut will grow approximately as $S_{E}^{\text{max}}\sim L$. Thus, there exist special values of $\phi$ for which the entanglement diverges across a cut.

In an actual sample, this exact resonances will not occur, only approximate resonances. Each sample can be thought of as sampling $\phi$ $L$ times, and so one will typically be at most $\sim 1/L$ away from the exact resonance point, and sites will be detuned by $\sim 1/L$. The effective hopping between two sites distance $x$ apart scales as $\sim e^{-x/\xi}$ for localization length $\xi$. Therefore, for such a near resonance, only pairs of sites within range $x$ such that $e^{-x/\xi}\gtrsim 1/L$ will form superpositions. This leads to an entanglement going as $S_{E}^{\text{max}}\sim x\sim\log L$. Thus, a completely different mechanism in the quasiperiodic model is allowing for logarithmic growth in maximum entanglement.

To avoid a logarithmic growth of $S_{E}^{\text{max}}$, one should then avoid a quasiperiodic potential with reflection symmetry. This is possible by a potential such as $V(i)=\cos(2\pi\omega(i+\phi))+\sin(4\pi\omega(i+\phi))$, which does not exhibit any even reflection symmetry. The $S_{E}$ distribution for this model has a very sharp cutoff and the maximum entanglement therefore does not keep growing as $\log L$.