Exploring many body localization and thermalization using semiclassical method

TL;DR

This paper evaluates the Discrete Truncated Wigner Approximation (DTWA) for studying many-body localization and thermalization, showing it captures key dynamical features in a 1D spin chain and could be useful for complex quantum systems.

Contribution

The study demonstrates that DTWA can effectively reproduce dynamical signatures of MBL and thermal phases, especially at short times, and extends its potential use to higher-dimensional and long-range interacting systems.

Findings

DTWA accurately reproduces dynamical signatures in both phases.

DTWA captures logarithmic entanglement growth in MBL.

Good agreement with exact solutions at short times.

Abstract

The Discrete Truncated Wigner Approximation (DTWA) is a semi-classical phase space method useful for the exploration of Many-body quantum dynamics. In this work we investigate Many-Body Localization (MBL) and thermalization using DTWA and compare its performance to exact numerical solutions. By taking as a benchmark case a 1D random field Heisenberg spin chain with short range interactions, and by comparing to numerically exact techniques, we show that DTWA is able to reproduce dynamical signatures that characterize both the thermal and the MBL phases. It exhibits the best quantitative agreement at short times deep in each of the phases and larger mismatches close to the phase transition. The DTWA captures the logarithmic growth of entanglement in the MBL phase, even though a pure classical mean-field analysis would lead to no dynamics at all. Our results suggest the DTWA can become a useful method to investigate MBL and thermalization in experimentally relevant settings intractable with exact numerical techniques, such as systems with long range interactions and/or systems in higher dimensions.

Figures

Equations

$$\hat{H}=J\sum _{i=1}^{N-1}{\hat{\mathbf{s}}}_i\cdot {\hat{\mathbf{s}}}_{i+1}+\sum _{i=1}^N{h}_i{\hat{s}}_i^z,$$

$$S_{\gamma }\left(\hat{\rho }\right)=\frac{1}{1-\gamma }{\log }_2\left[\mathrm{tr}\left\{{\hat{\rho }}^{\gamma }\right\}\right],$$

$$S_2\left({\hat{\rho }}_{i,j}\right)=-{\log }_2\left\{\frac{1}{4}\left[1+g_{z,0}^2+g_{0,z}^2+g_{z,z}^2+8{\left|g_{+,-}\right|}^2\right]\right\},$$

$$⟨\hat{O}⟩=\mathrm{tr}\left\{\hat{\rho }\hat{O}\right\}=\int \int \mathrm{d}p\mathrm{d}qW(q,p)\mathcal{O}(q,p,t).$$

$$⟨\hat{O}\left(t\right)⟩=\sum _{\mu \in 4^N}{w}_{\mu }\left(0\right){\mathcal{O}}_{\mu }\left(t\right).$$

$$\frac{d{\mathbf{r}}_{i,\mu }}{dt}=2J{\mathbf{r}}_{i,\mu }\times \left[\left({\mathbf{r}}_{i+1,\mu }+{\mathbf{r}}_{i-1,\mu }\right)+h_i{\mathbf{e}}_z\right],$$

$$MSE=\frac{{\sum }_{i=1}^M(f_i-g_i{)}^2}{M}$$

Full text

I Introduction

Since the seminal discovery of Anderson localization Anderson (1958) it has been known that the transport of an isolated quantum particle can be suppressed by the presence of disorder. While the conditions of this disorder driven localization have become well established for non-interacting particles Lee and Ramakrishnan (1985), our understanding of localization of interacting systems is still in rapid development (see Nandkishore and Huse (2015) for a recent review). This latter phenomenon, now known as many-body localization (MBL) Gornyi et al. (2005); Basko et al. (2006); Pal and Huse (2010a), has gained significant interest motivated by its fundamental connection to failure of thermalization in closed systems Pal and Huse (2010a); Žnidarič et al. (2008); Canovi et al. (2011); Cohen et al. (2016), emergent integrability Serbyn et al. (2013); Huse et al. (2014); Ros et al. (2015); Imbrie (2016); Chandran et al. (2016); Geraedts et al. (2017), new types of quantum order with possible applications for improved quantum memories Huse et al. (2013); Pekker et al. (2014); Vosk and Altman (2014), due to its novel response properties Khemani et al. (2015); Gopalakrishnan et al. (2015) and entanglement structures Bardarson et al. (2012a); Geraedts et al. (2016), and due to the emergence of experimental platforms for investigation of such physics Kondov et al. (2015); Schreiber et al. (2015); Smith et al. (2016); y. Choi et al. (2016).

In spite of intensive effort, our current ability to numerically simulate MBL systems faces severe limitations. In particular, currently numerical techniques consist mostly of exact diagonalization, the strong disorder renormalization group Vosk and Altman (2013); Vosk et al. (2015); Potter et al. (2015), and variants of the density matrix renormalization group (DMRG) Pekker and Clark (2017); Khemani et al. (2016), and are largely limited to one dimensional systems with short range interactions. Given the experimental and theoretical interest in higher dimensional systems and/or systems with long range interactions 111We note that it has been argued De Roeck and Huveneers (2017) that MBL may only exist in one dimension with short range interactions. Even if this argument is correct, however, the obstruction to MBL comes from ‘rare regions’ and only manifests itself on asymptotically long length scales and timescales, and so the quantum dynamics of disordered systems in higher dimensions and/or with longer range interactions remains experimentally relevant and theoretically interesting, there is a pressing need for more widely applicable numerical tools for exploring these more general situations Altman and Vosk (2015).

A recently developed semi-classical phase space method known as Discrete Truncated Wigner Approximation (DTWA) may help to overcome these limitations, as it has been been shown to correctly capture many aspects of the non-equilibrium behavior of many-body quantum spin lattices Schachenmayer et al. (2015a, b); Pucci et al. (2016, ); Piñeiro Orioli et al. . As DTWA is based on solving the system time evolution by classical mean field equations, it has the advantage that its computational complexity increases only polynomially with system size. Also, DTWA is not slowed down by the growth of quantum entanglement, which can hinder, for instance, DMRG methods. Given these advantages, we seek to investigate the utility of DTWA for addressing questions regarding many body localization and thermalization. In this work, we benchmark DTWA by studying the 1D nearest neighbor Heisenberg model in random uniaxial field, which is one of the most thoroughly studied systems exhibiting MBL (see e.g. Luitz et al. (2015)). We focus on dynamical signatures of localization that are both suitable to be computed by DTWA, and have also proven experimentally useful in characterizing MBL Schreiber et al. (2015); Smith et al. (2016). We note that, insofar as DTWA is a ‘mean field’ method, one dimensional systems with short range interactions are expected to be the least favorable setting for this method. Nevertheless, upon benchmarking DTWA to exact diagonalization, we find that DTWA works not just in the short time limit (when it is formally exact), but even at moderately long times. At intermediate times, DTWA is able to reproduce the main dynamical features of both the MBL phase and the thermal phase. The agreement persists up to longer times as one goes deeper into either phase. We attribute this to both extremes being in some sense classical - the deep MBL regime because there is very little entanglement, and the strongly delocalized regime because there is so much entanglement that reduced density matrices look thermal Deutsch (1991); Srednicki (1994); Rigol et al. (2008). We also note that DTWA performs better for some diagnostics than others, and discuss under what circumstances we expect the method to faithfully capture the quantum dynamics.

This paper is organized as follows. In Sec. II we introduce the model studied and the dynamical quantities that characterize its MBL behavior. In Sec. III we explain the semiclassical phase-space DTWA method. In Sec. IV, we present and discuss our results, concluding with some final remarks in Sec. V.

II Spin model and MBL signatures

In the following we study the dynamical properties of the 1D Heisenberg model with nearest-neighbor interactions and a random uniaxial field described by the Hamiltonian

[TABLE]

where $\hat{s}_{i}^{\alpha}=\tfrac{1}{2}\hat{\sigma}_{i}^{\alpha}$ are spin 1/2 operators ($\hat{\sigma}_{i}^{\alpha}$ are Pauli matrices), $J$ is the interaction strength, and the values of the random potential $h_{i}$ (measured in units of $J$) are uniformly distributed $h_{i}\in\left[-h,h\right]$ with $h$ the disorder strength. Previous studies have shown that this model features an MBL transition at $h_{c}\approx 3.5$ in the thermodynamic limit and at infinite temperature Pal and Huse (2010b). At finite temperature the eigenstates exhibit a localization edge between $h\approx 2$ and $h\approx 3.5$ Luitz et al. (2015); Serbyn et al. (2015).

Here we study the dynamics when the system is initialized in the Néel state, $\left|\Psi(0)\right\rangle=\left|\uparrow\downarrow\uparrow\downarrow\ldots\right\rangle$. For this setting, the critical point of $h_{c}\approx 3.5$ should be taken only as rough reference value. For instance, in the thermodynamic limit the Néel state is expected to localize below the infinite temperature critical point $h_{c}$ due to its average energy eigenstates crossing the localization edge at lower disorder strengths. In addition, finite sizes effects may affect the disorder strengths at which signatures of the MBL transition are observed Potter et al. (2015); Kjäll et al. (2014).

We probe the dynamics by initializing the system in a Néel state and by then computing dynamical observables which have been previously used to characterize the phase Pal and Huse (2010b); Bardarson et al. (2012b); Kjäll et al. (2014); Vasseur et al. (2015); Schreiber et al. (2015); Smith et al. (2016). These include the imbalance $I(t)$, Quantum Fisher Information (QFI) denoted by $\mathcal{F}(t)$, Rényi entropy $S_{2}(t)$ over a two-spin subsystem, and entanglement (Von Neumann) entropy $S_{1}(t)$ of one half of the system. The imbalance is a measure of the Néel state’s staggered order across the spin chain. It is defined as $I\left(t\right)=\sum_{i}(-1)^{i+1}\left\langle\hat{\sigma}_{i}^{z}\left(t\right)\right\rangle$. The QFI associated to an operator $\hat{\mathcal{O}}$ for a spin chain in a pure state is given by $\mathcal{F}_{\hat{\mathcal{O}}}\equiv 4\left(\left\langle\hat{\mathcal{O}}^{2}\right\rangle-\left\langle\hat{\mathcal{O}}\right\rangle^{2}\right)=4\Delta\hat{\mathcal{O}}$. Here $\hat{\mathcal{O}}=\sum_{i=1}^{N}\hat{O}_{i}$, and each $\hat{O}_{i}$ is a spin $1/2$ operator acting on spin $i$. The QFI is a witness of multi-partite entanglement if $\mathcal{F}_{\hat{\mathcal{O}}}/N>1$ Hyllus et al. (2012); Tóth and Apellaniz (2014); Apellaniz et al. (2017). Due to our particular choice of the initial state, it is natural to choose the imbalance as the operator $\hat{\mathcal{O}}$, so that QFI is proportional to its variance $\mathcal{F}(t)=4\Delta I\left(t\right)$. Finally, both Rényi entropy $S_{2}$ and the von-Neumann entanglement entropy $S_{1}$ can be thought as belonging to a continuous family of entropy measures defined as,

[TABLE]

where $\gamma\geq 1$. The Rényi entropy corresponds to $S_{\gamma=2}$, while the von-Neumann entropy $S_{1}=\lim_{\gamma\to 1}S_{\gamma}$. Given that the total spin $\hat{S}_{z}$ is a conserved quantity, and the Néel state is an eigenstate of $\hat{S}_{z}$, the Rényi entropy over a two-spin subsystem with spins labeled as $i$ and $j$ can be particularly expressed as,

[TABLE]

where $g_{\alpha,\beta}=\left\langle\hat{\sigma}_{i}^{\alpha}\hat{\sigma}_{j}^{\beta}\right\rangle$ are two particle moments of Pauli matrices ($\hat{\sigma}^{0}$ is the identity matrix).

III Semi-classical phase-space method

The DTWA approximately describes the dynamics of quantum observables as statistical averages over classical phase-space trajectories Schachenmayer et al. (2015a, b). It starts from a phase-space formulation of quantum mechanics. The method was inspired by the well known Truncated Wigner Approximation derived for systems with continuous degrees of freedom, described e.g. by the canonical $\left(q,p\right)$ coordinates along each dimension Polkovnikov (2010). In a continuous phase-space, the state of the system (the density matrix $\hat{\rho}$) and any observable $\hat{O}$ are represented by real functions of phase-space coordinates, $W(q,p)=\left(\hat{\rho}\right)_{W}$ and $\mathcal{O}(q,p)=\left(\hat{O}\right)_{W}$ respectively, and any expectation value is computed as,

[TABLE]

The function $W(q,p)$ is known as the Wigner quasi-probability distribution, and the function $\mathcal{O}(q,p)$ is known as the Weyl symbol of operator $\hat{O}$. Either the Wigner function, the Weyl symbol, or both can vary in time depending on the dynamical representation. The previous formula is in the Heisenberg picture, where the operators are the ones that change in time.

A similar formulation can be deduced for systems with discrete degrees of freedom, such as the spin chain of $N$ qubits considered in Eq. (1) Wootters (1987). In the discrete case, the phase-space consists of $4^{N}$ points instead of the continuous $\left(q,p\right)$ plane, and a formula analogous to Eq. (4) reads,

[TABLE]

Up to this moment, no approximation has been made: if both $w_{\mu}\left(0\right)$ and $\mathcal{O}_{\mu}\left(t\right)$ are known for every phase-space point, then the previous formula gives the exact quantum solution. The total number of phase-points $\mu$ can become very large, so that a first practical approximation should be done. Any time the Wigner function can be seen as a (positive) probability distribution on phase-space, Eq. (5) can be computed via a Monte-Carlo sampling over a number of phase-space points large enough to obtain an accurate estimate. The initial state in consideration fulfills this positive probability condition.

To compute the Weyl symbols $\mathcal{O}_{\mu}\left(t\right)$ efficiently, a more important approximation is required. There is a hierarchy of equations that determine the time evolution of these symbols Pucci et al. (2016), which can be truncated up to the number of spins in the operators that are considered independent. To first order, i.e. if only single particle operators are considered independent, the equations for $\mathcal{O}_{\mu}\left(t\right)$ can be approximated to obey the classical mean field equations,

[TABLE]

where $\mathbf{r}_{i,\mu}\equiv\left\{\mathcal{S}^{\alpha}_{i}\right\}=\left\{\left(\hat{\sigma}_{i}^{\alpha}\right)_{W}\right\}$ is the Bloch vector of the $i$-th spin. The index $\mu$ accounts for the specific phase-space point, which in turn defines the initial value of $\mathbf{r}_{i,\mu}$. Equation (6) is equivalent to assume that all higher order Weyl symbols like $\mathcal{S}^{\alpha,\beta}_{i,j}=\left(\hat{\sigma}_{i}^{\alpha}\hat{\sigma}_{j}^{\beta}\right)_{W}$ are approximated as products of single particle Weyl symbols, i.e., $\mathcal{S}^{\alpha,\beta}_{i,j}\approx\mathcal{S}^{\alpha}_{i}\mathcal{S}^{\beta}_{j}$ for $i\neq j$. This approximation should be exact at least at short times for initially uncorrelated spin states as the Néel state in consideration. Thus, the DTWA can be described as a semi-classical method where the time evolution is approximated by classical trajectories whose initial conditions are sampled according to the exact initial Wigner function.

It should be noticed that there is a crucial difference between simple mean field solutions and DTWA. In a simple mean field approach, the vectors $\mathbf{r}_{i,\mu}$ in Eq. (6) represent expectation values and not Weyl symbols, and so the set of initial conditions used in Eq. (6) is very different. For the case considered here wherein we start from a Néel state, a simple mean field approach would predict no dynamics at all. Also, the mean field approach predicts zero correlations between single site observables at any time. Instead, DTWA does develop such correlations in time, further extending its applicability to capture the exact quantum evolution.

IV Results

In the literature, a commonly used dynamical witness of the MBL phase is the growth of the half chain entanglement entropy $S_{1}\left(\hat{\rho}_{N/2}(t)\right)$ Bardarson et al. (2012b); Altman and Vosk (2015). However, this quantity is hard to measure experimentally. It is also computationally difficult to estimate within the DTWA. Instead, there are other indicators of localization that are more accessible both to experiments and to the DTWA. One of them is the loss of any initial local structure of the system in the ergodic phase. In the localized phase, on the contrary, there is an extensive set of Local Integrals of Motion (LIOM) that are able to partially preserve the initial order. For instance, in the very large disorder limit, the Pauli operators $\hat{\sigma}_{i}^{z}$ become a well defined set of LIOM, and any eigenstate of these operators will remain frozen. Therefore, the initial imbalance $I(t)$ remains constant for very large disorder strengths. On the other hand, in the thermal phase, any initial imbalance is expected to decay to zero during the dynamics. At weaker disorder strengths, but still in the localized phase, the LIOM become less local, so that the initial imbalance is no longer a conserved quantity, and cannot stay at its initial value; however, it is still partially conserved and should saturate to a non-zero value in the long term Pal and Huse (2010b); Vasseur et al. (2015). Figure 1 shows this tendency in the exact solutions. We observe a change of behavior around $h\approx 3$, which is in close agreement with the critical value expected in the thermodynamic limit Luitz et al. (2015); Serbyn et al. (2015). The DTWA results are in agreement with the exact dynamics of imbalance at short times, namely $Jt\lesssim 1$. Furthermore, in the intermediate regime $Jt\lesssim 100$ we find good qualitative agreement away from the critical region. We note that the DTWA imbalance decays for all values of disorder strength at sufficiently long times, even deep in the localized phase. This is due to the non-linear nature of the mean field equations used in the DTWA which always leads to long time thermalization, similar to Ref. Oganesyan et al. (2009). Nevertheless the imbalance decay is largely suppressed in the localized regime, and the overall slow-down can be used as probe of localization. In fact, the discrepancy with the exact quantum dynamics is delayed further in time as disorder strength increases. In the thermal phase, this decay is real and dictated by the quantum evolution and the DTWA captures well the imbalance dynamics. The agreement in the ergodic phase can be understood by the fact that even though the DTWA neglects higher order correlations between spins, single particle observables like imbalance are not sensitive to those correlations. The role of those higher order correlations when tracing out all but a single spin reduces to a high temperature bath which can also be mimicked by the statistical mixture of classical trajectories. Another difference between the DTWA results and the exact ones are transitory oscillations at intermediate times which are not present in the DTWA solution. These oscillations are fine features generated either by boundary effects or sporadic few-body clusters, and they are not fundamental characteristics of the two phases.

A similar tendency of agreement between DTWA and exact results is seen for two-site Rényi entropy, as can be seen in Fig. 2. This quantity is a measure of the entanglement between the two-sites-subsystem and the rest of the system. Equation (3) shows that $S_{2}$ is determined by a handful of one and two sites operators. In the long term, both the MBL ($h\geq 3$) and thermal regimes exhibit a monotonic growth of $S_{2}$, but this growth is slow and logarithmic in the localized regime, while in the thermal regime $S_{2}$ quickly saturates to a final value. This behavior has strong similarities with known results of entanglement entropy evaluated by tracing over half of the system (see for example Bardarson et al. (2012b)). The main difference is that the upper bound $\max\left(S_{2}\right)=2$ is much lower for a two-spin subsystem. In fact, the saturation values reached by the two-site $S_{2}$ in the thermal phase are not that different than the ones reached by the same observable in the localized phase. In general terms, the DTWA is able to reproduce the main features exhibited by the different regimes, though it overestimates the logarithmic growth in the MBL phase. As with imbalance, the disparity between the DTWA and the exact results gets smaller as the disorder increases. From Eq. (3) one expects that most of the arguments that explain the ability of DTWA to reproduce the imbalance evolution still hold for $S_{2}$. However, given that the subsystem is composed of two sites, there can be an extra dependence on spin-spin correlations that are not well approximated by DTWA. This explains the discrepancy with the exact results at intermediate times in the localized phase. After thermalization, the spin-spin correlations lose their significance and the agreement with DTWA is recovered.

Contrary to imbalance and two-site Rényi Entropy, the QFI is a quantity that witnesses the build up of spin-spin correlations. Evidence of MBL has already been found using a QFI associated to the particle imbalance in a recent trapped ion experiment that simulated the 1D Ising model Smith et al. (2016). The fact that the long term behavior of entropy and QFI is logarithmic in time for the MBL phase is probably a coincidence restricted to short range models, and it may not be generally guaranteed. We observe the logarithmic growth of QFI in the DTWA results, showing a strong analogy of this semi-classical approach to the exact quantum system even to the level of spin fluctuations. However, a quantitative agreement between the exact and DTWA curves is now only possible in the large disorder limit, that is DTWA accurately captures the growth of QFI deep in the localized phase, but not deep in the thermal phase. Moreover the time scale over which DTWA accurately captures the growth of QFI grows longer as one goes deeper into the localized phase.

We present a summary of the suitability of the DTWA for reproducing several dynamical indicators of MBL in Fig. 4. The curves correspond to a Mean Squared Error (MSE) of the form,

[TABLE]

where $M$ is the number of computed points of the quantities (which lie in the interval $Jt\in\left[0,120\right]$), and $f_{i}$ ($g_{i}$) are the exact (DTWA) values. Although there are other ways to statistically analyze the discrepancy between exact and DTWA results, we have found that Eq. 7 is a very simple and common approach that express the visually qualitative closeness between the curves in the figures presented so far. The units of the MSE have been left unspecified as there are very different quantities displayed in the same figure, and what is relevant is the tendency of the agreement for each curve as the disorder strength varies. All MSE curves have maxima around the critical disorder strength and go to smaller values in both thermal and localized extremes. Other absolute error measures like Root Mean Squared Deviation or Mean Absolute Error behave in a similar fashion. We also included in the inset of Fig. 4 an example of a relative or normalized statistical measure of the error, namely the Coefficient of Variation of the Root-Mean-Squared-Deviation (CVRMSD), which is defined as $\varepsilon=\sqrt{MSE}/\bar{f}$, where $\bar{f}=\sum_{i=1}^{M}f_{i}$. The CVRMSD has the advantage of being dimensionless, so all curves can be put in the same scale. However, the CVRMSD has the disadvantage of being too sensitive to fluctuations in the exact solution with respect to the DTWA for quantities like the QFIM and Rényi Entropies in the localized phase, which remain very close to zero during the dynamics. We have already seen that DTWA is not well suited for reproducing those fluctuations, but only to show the average tendency of the relevant MBL observables. As it is the average tendency which is the most relevant aspect for MBL, we think that the MSE provides a better measure for the suitability of DTWA for exploring MBL dynamics.

V Final remarks

We have shown that the semi-classical picture provided by DTWA is able to reproduce the general features of the quantum dynamics in both the localized and thermal phases. The agreement persists up to time scales that grow longer as one goes deeper into either phase. Particular diagnostics we have considered include the memory of the initial conditions in imbalance, growth of Rényi entropy for a two spin subsystem, and quantum Fisher information. While the discrepancy between exact methods and DTWA grows larger at long times, the DTWA works reasonably well over a range of intermediate timescales (which can be made longer by going deeper into the respective phases). It is worth noting that the DTWA can capture the main features of the dynamics over the range of timescales relevant for experiments e.g. Schreiber et al. (2015). It is also worth noting that in this paper we have focused on the short-range 1D Heisenberg model, which is presumably the ‘worst case’ scenario for a method of this type Schachenmayer et al. (2015a, b); Pucci et al. (2016, ); Piñeiro Orioli et al. . The DTWA is expected to become more accurate for problems with higher coordination number (either long range interacting problems or higher dimensional problems). These long range interacting and/or high dimensional problems are formidably difficult to tackle using exact techniques. However, when treated using DTWA, the computational cost scales only polynomially with system size. As such, we expect DTWA can become a powerful method for investigating localization and thermalization in these more complex and less explored settings, which are inaccessible to existing numerical techniques, but are of relevance for ongoing experiments.

Acknowledgements.

We thank Robert Lewis-Swan, Asier Piñeiro-Orioli, Shainen Davidson, and Anatoli Polkovnikov for insightful discussions. We acknowledge support from NSF-PHY-1521080, JILA-NSF-PFC-1125844, DARPA (W911NF-16-1-0576 through ARO), MURI-AFOSR and AFOSR. R.N. is supported by the Foundational Questions Institute (fqxi.org; grant no. FQXi-RFP-1617) through their fund at the Silicon Valley Community Foundation.