Efficient Representation of Fully Many-Body Localized Systems Using Tensor Networks

We propose a tensor network encoding the set of all eigenstates of a fully many-body localized system in one dimension. Our construction, conceptually based on the ansatz introduced in Phys. Rev. B 94, 041116(R) (2016), is built from two layers of unitary matrices which act on blocks of $\ell$ contiguous sites. We argue that this yields an exponential reduction in computational time and memory requirement as compared to all previous approaches for finding a representation of the complete eigenspectrum of large many-body localized systems with a given accuracy. Concretely, we optimize the unitaries by minimizing the magnitude of the commutator of the approximate integrals of motion and the Hamiltonian, which can be done in a local fashion. This further reduces the computational complexity of the tensor networks arising in the minimization process compared to previous work. We test the accuracy of our method by comparing the approximate energy spectrum to exact diagonalization results for the random-field Heisenberg model on 16 sites. We find that the technique is highly accurate deep in the localized regime and maintains a surprising degree of accuracy in predicting certain local quantities even in the vicinity of the predicted dynamical phase transition. To demonstrate the power of our technique, we study a system of 72 sites, and we are able to see clear signatures of the phase transition. Our work opens a new avenue to study properties of the many-body localization transition in large systems.

Popular Summary

The thermal behavior of some quantum systems can be quite different from everyday experience. In the presence of a source of heat, most objects eventually thermalize: Their temperature begins to match that of their surroundings. In 1958, physicist Philip Anderson conjectured that a quantum system isolated from its environment might fail to thermalize in the presence of disorder, a phenomenon now known as many-body localization. Such systems can retain their initial conditions for long times. Experiments have only recently begun to probe the dynamics of many-body localization in solid-state devices and ensembles of ultracold atoms. These investigations open the door to realizing some quantum phenomena on macroscopic scales or even the development of robust memories for quantum computing. In our work, we have developed a mathematical framework that allows for the efficient simulation of these intriguing states of matter for large systems.

The total energy of a quantum-mechanical system can be described by a mathematical tool called a Hamiltonian; special features of the Hamiltonian—known as eigenstates and eigenenergies—fully encode the quantum dynamics. Evaluating these eigenstates for a system with many particles becomes exponentially harder as the size of the system is increased. By using tensor networks and the area law of entanglement, we evaluate the entire energy spectrum of a large, one-dimensional, many-body-localized system. This approach is exponentially faster than any existing method, obtains very high accuracy, and is able to capture clear signatures of the transition from the localized phase to the thermal phase that exists at lower disorder.

These findings pave the way for studying the critical properties of the novel quantum phase transition between thermal and localized phases, and may also allow us to search for many-body localization in two dimensions.

Figure 1

Tensor network $\tilde{U}$ as proposed in Ref. [48] with $n$ layers of $4\times 4$ unitaries $\{u_{x,y}\}$. Since the unitary $\tilde{U}$ is supposed to approximately diagonalize the Hamiltonian, the corresponding approximate eigenstates $|{\tilde{\psi }}_{i_1,\dots ,i_N}⟩$ are the states obtained by fixing the lower open indices in the figure to be $i_1,\dots ,i_N$.

Figure 2

Construction of the unitary $\tilde{U}$ in terms of unitaries $u_{x,1}$, $u_{x,2}$ acting on $\ell$ sites (in this example, $\ell =4$). Again, the approximate eigenstates $|{\tilde{\psi }}_{i_1,\dots ,i_N}⟩$ are the states obtained by fixing the lower open indices in the figure to be $i_1,\dots ,i_N$.

Figure 3

Sum of tensor network contractions, which yields the second term, $\sum _{i=1}^N\mathrm{tr}\mathbf{(}(H{\tilde{\tau }}_i^z{)}^2\mathbf{)}$, in Eq. (5). The multiplications from left to right in Eq. (6) correspond to top to bottom in the figure. The indices of the lower wiggly lines are to be contracted with those of the corresponding upper wiggly lines. For a given position $i$ of the ${\sigma }^z$ operator and arbitrary positions $j$, $k$ of the two-body Hamiltonian terms, all unitaries of the lower layer (i.e., $u_{n,1}$), apart from the ones directly connected to the ${\sigma }^z$ operators, cancel with their adjoints and can be replaced by identities (i.e., straight vertical lines). In the example in the figure, this corresponds to all unitaries $u_{n,1}$ for $n\ne x$. Furthermore, all unitaries of the second layer ($u_{n,2}$) which are not directly connected to the remaining ones of the first layer cancel and can be substituted by identities. This implies that the $x$th summand in Eq. (6) depends only on the unitaries $u_{x+1,1}$, $u_{x,2}$, $u_{x+1,2}$. The contraction corresponding to this term is shown in Fig. 4.

Figure 4

Decomposition of the figure of merit (5) into local terms resulting in Eq. (6). Again, the indices of the lower wiggly lines are to be contracted with those of the corresponding upper wiggly lines. The tensor network shown is obtained after replacing mutually canceling unitaries in Fig. 3 by identities (vertical lines). Those forming closed loops yield a factor of 2 each, which results in the prefactor $2^{-N+2\ell +2}$ of $f_x$ shown in the figure. Terms $j$, $k$, where $h_j$ or $h_k$ are not connected to $u_{x-1,2}$ or $u_{x,2}$, yield contributions that are independent of all unitaries $u_{x,y}$ and can thus be neglected in the local definition of our figure of merit. Note that the precise positions of ${\sigma }_i^z$, $h_j$, and $h_k$ depend on the indices $i$, $j$, and $k$ that are being summed over, and thus the graphic depicts one example configuration. For details of our contraction scheme, see Appendix pp1.

Figure 5

Comparison of the optimized tensor network $\tilde{U}$ for $\ell =2$ (a), $\ell =4$ (b), and $\ell =8$ (c) with exact diagonalization for $N=16$ and ten disorder realizations at $W=6$. The energy differences $\mathrm{\Delta }E$ were obtained by ordering the diagonal elements of ${\tilde{U}}^{†}H\tilde{U}$ in each spin-$z$ sector and subtracting them from the ordered exact eigenvalues of the Hamiltonian $H$ in the corresponding spin-$z$ sector. The plots show the concatenation of the data of all spin-$z$ sectors and disorder realizations. The optimized SCN figure of merit was, on average, $\overline{f_{\ell =2}}/2^N=1.011$, $\overline{f_{\ell =4}}/2^N=0.194$, $\overline{f_{\ell =8}}/2^N=0.0203$.

Figure 6

Comparison of the optimized tensor network $\tilde{U}$ for $\ell =8$ with exact diagonalization for $N=16$, showing the data for ten disorder realizations at $W=6$. The approximate eigenstates are given by the columns of $\tilde{U}$. (a) Distribution of the overlap of the exact and the matched eigenstates. (b) Distribution of the expectation value of ${\sigma }_i^z$ over the sites and the eigenstates obtained from exact diagonalization (light brown) and the TN (blue).

Figure 7

Comparison of the optimized tensor network $\tilde{U}$ using four layers in the ansatz of Fig. 1 with exact diagonalization for $N=16$ and the same ten disorder realizations as in Fig. 5 (i.e., $W=6$), the unitaries of the TN being $U(1)$ symmetric and real. The energy differences shown have been calculated in each spin-$z$ sector separately and concatenated, also over disorder realizations. The improvement over Fig. 5 ($\ell =2$) is minuscule. The average optimized SCN figure of merit was $\overline{f}/2^N=0.747$. (For the full parametrization of the unitaries, we were not able to produce the plot because of the limitations of virtual memory.)

Figure 8

Scaling of the SCN figure of merit as a function of system size $N$ for $W=6$ and 100 different disorder realizations, optimized using unitaries of block sizes $\ell =2$, 4, 6 and ten disorder realizations for $\ell =8$. For given $N$, the same hundred (ten) disorder realizations were taken for all values of $\ell$. As discussed in the main text, we expect $f/2^N\propto N$, which is consistent with the numerical results. The error bars denote the error of the mean calculated from the distribution over disorder realizations. The inset shows an enlargement of the data for $\ell =6$ and $\ell =8$. There, the symbol size is at least the size of the error bars.

Figure 9

Figure of merit, SCN, for $N=72$ as a function of disorder strength $W$ for $\ell =2$, 4, 6, 8. The same hundred (ten) disorder configurations were taken for $\ell =2$, 4, 6 ($\ell =8$) and all choices of $W$, adjusting only the overall prefactor of the random magnetic fields. The error bars denote the error of the mean calculated from the distribution over disorder realizations.

Figure 10

Optimized SCN averaged over disorder realizations as in Fig. 9 shown as a function of $\ell$ for various disorder strengths $W$. The error bars denote the error of the mean calculated from the distribution over disorder realizations. The inset shows the same data on a log-log plot (with symbol size at least as big as the error bars). For $W=8$, 16, the decay of the SCN is approximately exponential for $\ell =2$, 4, 6 but deviates for $\ell =8$, probably because the algorithm tends to get stuck in local minima. For $W=0.5$, 2, the decay with $\ell$ is, to a good approximation, an inverse power law (exponents $-1.3$ and $-1.9$, respectively).

Figure 11

Comparison of the eigenstates from the optimized TN $\tilde{U}$ for $\ell =6$ and exact eigenstates from exact diagonalization for $N=12$ and 100 disorder realizations at disorder strengths (a) $W=2$, (b) $W=3$, (c) $W=4$, and (d) $W=6$. We present the distribution of the expectation value of ${\sigma }_i^z$ over the sites, eigenstates, and disorder realizations from exact diagonalization (blue) and the TN (light brown).

Figure 12

Standard deviation of the half-cut entanglement entropy (${\sigma }_S$) for 100 disorder realizations averaged over the (approximate) eigenstates as a function of disorder strength for the TN with $\ell =2$, 4, 6, and exact diagonalization. The system size is $N=12$. The inset shows the average entanglement entropy ($\overline{S}$) for the four cases.

Figure 13

Standard deviation of the eigenstate-averaged half-cut entanglement entropy, followed by an average over entanglement cuts, $⟨{\sigma }_S{⟩}_{\mathrm{cuts}}$ (defined in the text), as a function of disorder strength. The system size is $N=72$, and the calculations were performed on 100 disorder realizations using the TN with $\ell =2$, 4, 6, and 10 disorder realizations using $\ell =8$. For $\ell =6$ and $\ell =8$, the entanglement entropy was sampled over at least 1% of all eigenstates, leading to an estimated relative error of about 5% (marked by error bars) of the mean and standard deviation of the entropies. The inset shows the corresponding eigenstate-averaged entanglement entropy, followed by an average over entanglement cuts ($⟨\overline{S}{⟩}_{\mathrm{cuts}}$). The errors in the inset are, at most, as big as the size of the symbols.

Figure 14

(a) The left-hand side of the equation is the same as Fig. 4. Here, $u_{x-1,2}$ and $u_{x,2}$ have been moved across the upper boundary coming back from the bottom due to the trace. The dashed orange boxes indicates how the tensors are blocked together in order to speed up the contraction. Note that the precise positions of ${\sigma }_i^z$, $h_j$, and $h_k$ depend on the indices $i$, $j$, $k$ that are being summed over; i.e., the graphic on the left-hand side shows only one example configuration. The right-hand side is obtained after using the substitutions shown in (b–f). Here, ${\sum }_{i,j,k}^{\prime }$ is the same sum as on the left-hand side, with $j=k=(x-\frac{1}{2})\ell$ and $j=k=(x+\frac{3}{2})\ell$ excluded. These terms correspond to the second and third terms on the right-hand side, respectively. The most efficient way to contract the tensor networks on the right-hand side is by contracting $A_j^{(n)}$ with $A_k^{(m)}$ ($\tilde{A}$ with $\tilde{A}$) and $B_j^{(n)}$ with $B_k^{(m)}$ ($\tilde{B}$ with $\tilde{B}$) and afterwards the resulting tensors with the $Q_i$’s. The biggest matrices to be multiplied come from the second and third terms and are of size $2^{\ell }\times 2^{\ell +1}$, corresponding to a cost of $2^{3\ell +1}$. There are $\ell$ such contractions, respectively, but to leading order $4{\ell }^3$ contractions coming from the first term (multiplication of $2^{\ell }\times 2^{\ell }$ matrices; i.e., the overall computational cost is of order ${\ell }^3{2}^{3\ell }$). (b) Definition of $Q_i$. (c) Definition of $A_{j=(x-1/2)\ell }^{(n=1)}$. The upper index ($n$) is trivial in that case. Note that $B_{j=(x+3/2)\ell }^{(n=1)}$ is defined analogously (for the Hamiltonian term $h_{j=(x+3/2)\ell }$). (d) Definition of $\tilde{A}$ for $j=(x-\frac{1}{2})\ell$. Here, $\tilde{B}$ is defined analogously. (e) Definition of $A_j^{(n=1)}$, where the index ($n$) is trivial and $(x-\frac{1}{2})\ell <j<(x+\frac{1}{2})\ell$. Note that $B_j^{(n=1)}$ is defined analogously [for $(x+\frac{1}{2})\ell <j<(x+\frac{3}{2})\ell$]. For the $j$ indices we just specified for (c–e), we have $B_j^{(n=1)}=\mathbb{1}$. (f) For $j=(x+\frac{1}{2})\ell$, $A_j^{(n)}$ and $B_j^{(n)}$ are obtained by using a singular-value decomposition of the term $h_j$ as shown on the left and blocking $u_{x-1,2}^{†}$, $h_L^{(n)}$, $u_{x-1,2}$, and $u_{x,2}^{†}$, $h_R^{(n)}$, $u_{x,2}$, respectively. This case is the only one where the index ($n$) is nontrivial.

Figure 15

The tensor network contraction shown results in $2^{-N+2l+2}M$, with the $2^{\ell }\times 2^{\ell }$ matrix $M$ described in the main text if the unitary of which the derivative is being taken is cut out on the very top or very bottom (marked in red). The tensor network can be contracted with the same computational scaling as before (${\ell }^3{2}^{3\ell }$) by using the same blocking as in Fig. 14. The resulting blocks $A_j^{(n)}$, $B_k^{(m)}$, etc. are contracted in such an order that the one which contains the unitary to be varied is contracted with last. In other words, if $u_{x,1}$ is being varied (which is contained in the upper $Q_i$), one first contracts $A_j^{(n)}$ with $A_k^{(m)}$ and the resulting tensor with the lower $Q_i$ and subsequently with the contraction of $B_j^{(n)}$ and $B_k^{(m)}$.

Figure 16

(a) The left-hand side of the equation is the reduced density matrix $\rho$ obtained by tracing out sites $(N/2)+1,\dots ,N$ according to Eq. (c1). By canceling unitaries and their adjoints, the right-hand side is obtained, cf. Eq. (c2). As explained in the text, the nontrivial spectrum of $\rho$ is the same as the one of the reduced density matrix $\sigma$ defined in (b).