Time evolution of many-body localized systems in two spatial dimensions

Many-body localization is a striking mechanism that prevents interacting quantum systems from thermalizing. The absence of thermalization behavior manifests itself, for example, in a remanence of local particle number configurations, a quantity that is robust over a parameter range. Local particle numbers are directly accessible in programmable quantum simulators, in systems of cold atoms, even in two spatial dimensions. Yet, the classical simulation aimed at building trust in quantum simulations is highly challenging. In this work, we present a comprehensive tensor network simulation of a many-body localized systems in two spatial dimensions using a variant of an infinite projected entangled pair states algorithm. The required translational invariance can be restored by implementing the disorder into an auxiliary spin system, providing an exact disorder average under dynamics. We can quantitatively assess signatures of many-body localization for the infinite system: Our methods are powerful enough to provide crude dynamical estimates for the transition between localized and ergodic phases. Interestingly, in this setting of finitely many disorder values, which we also compare with simulations involving noninteracting fermions and for which we discuss the emergent physics, localization emerges in the interacting regime, for which we provide an intuitive argument, while Anderson localization is absent.

Figure 1

Initial state expressed in terms of iPEPS for (a) the physical state vector $|{\psi }_0{\rangle }_p$, which is a Néel state; (b) the auxiliary state vector $|{\psi }_0{\rangle }_a$, which is a product of equal superposition states; and (c) the overall initial state vector $|{\mathrm{\Psi }}_0\rangle$, which is the tensor product of the previous two states. The red patterns correspond to the classical interaction between the physical and the auxiliary states, which is required for introducing the disorder. All three states are iPEPS with bond dimension $D=1$ and the lattice extends indefinitely in all directions. A choice of a two-site unit cell in a checkerboard pattern is enough to exactly represent this configuration.

Figure 2

Averaged local particle number for even and odd sites in a free fermionic model for system sizes $L=10,20,40$ (markers). Averages are taken over 100 realizations. Blue shades are for discrete disorder (spin-1) and red shades for continuous disorder. The black dotted lines are iPEPS results for an infinite system. They are presented only to give a qualitative prediction of how the results using the free fermion simulation will change in the thermodynamic limit. A one-to-one quantitative comparison should not be made with the iPEPS results, since the plot shown is for $D=4$. However, the agreement between the two techniques is striking for the large disorder case $h=100$ since the entanglement growth is much slower here and the bond dimension does not play much of a role within the time scale presented here.

Figure 3

Single-particle spectrum for the Anderson model with discrete and continuous disorder for $h=100$ averaged over 100 realizations.

Figure 4

Cumulative inverse participation ratio for the Anderson model with discrete and continuous disorder for different disorder strengths. Lines are guides to the eye. When the energy levels for a certain value $\epsilon$ are not populated, no line is drawn.

Figure 5

Real-time evolution of the Heisenberg Hamiltonian starting from a Néel state. Left: Expectation value of the particle number as a function of time for the two different sites for the case with no disorder $h=0$. Big inset: Renyi entropies of the reduced density matrix of one site as a function of time for $\alpha =1$ and $1/2$. The entropies start saturating to their maximum value. Small inset: Accumulated local truncation error of one site. Right: The same evolution as above, but now with a disorder strength $h=2$ and for $d_A=2$. The simulation time has been extended slightly compared to the clean case without disorder. One also notices a slow down in the growth of local Renyi entropies. The simulations are done for $\delta t=0.1$ and 0.01 for $D=4$ and $\delta t=0.01$ for $D=5$. The results are consistent with different Trotter sizes as well as different bond dimensions, building trust in our simulations.

Figure 6

Similar to Fig. 5 but for stronger disorder and more disorder levels. Top: Relatively strong disorder $h=6$ but only two levels of disorder $d_A=2$. Middle: More levels of disorder $d_A=5$ but weak disorder strength $h=2$. Bottom: Many disorder levels $d_A=5$ as well as relatively strong disorder $h=6$. In all the plots, we show simulations with bond dimensions $D=4$ and 5 up to good agreement. This is also consistent with the growth of Renyi entropies, which are shown in the insets.

Figure 7

Particle imbalance $\mathcal{I}$ for various configurations of disorder dimensions $d_A$ and strengths $h$. We show the dynamics of the longest available times for the localized case (blue circles, $h=6$ and $d_A=5$), which is up to three hoppings. Also shown are the cases in which the particles do not localize (yellow circles with $h=6$, $d_A=2$, red circles with $h=2$, $d_A=5$, and green circles with $h=4$, $d_A=5$). The dynamics can be extrapolated using different polynomials such as linear, quadratic, fourth, and fifth degree fit, and one can notice the imbalance dropping to zero in all these cases. Also shown are the residuals corresponding to each fit (dashed lines). The linear fit has the largest error, while the fifth degree polynomial fits in this sense most accurately.

Figure 8

Crude estimate of the phase diagram as being assessed from dynamical localization as a function of the disorder dimension $d_A$ and the disorder strength $h$. The criterion to assign an ergodic or localized phase is whether the achievable simulation in time or the polynomial interpolation exhibits a localization of the imbalance in time. A disorder strength of $h=6J$ with at least four levels of disorder seems necessary to give rise to many-body localization in 2D.