Matrix product states approaches to operator spreading in ergodic quantum systems

TL;DR

This paper reviews tensor network methods, especially matrix product states, for studying operator spreading in ergodic quantum systems, comparing their accuracy and efficiency against exact simulations.

Contribution

It systematically compares TEBD and TDVP approaches in Schrödinger and Heisenberg pictures for operator spreading, highlighting their relative performance and limitations.

Findings

TDVP outperforms TEBD in Schrödinger picture.

Both methods accurately capture OTOC tails at short times.

Growth and saturation regimes are not well captured by the methods.

Abstract

We review different tensor network approaches to study the spreading of operators in generic nonintegrable quantum systems. As a common ground to all methods, we quantify this spreading by means of the Frobenius norm of the commutator of a spreading operator with a local operator, which is usually referred to as the out of time order correlation (OTOC) function. We compare two approaches based on matrix-product states in the Schr\"odinger picture: the time dependent block decimation (TEBD) and the time dependent variational principle (TDVP), as well as TEBD based on matrix-product operators directly in the Heisenberg picture. The results of all methods are compared to numerically exact results using Krylov space exact time evolution. We find that for the Schr\"odinger picture the TDVP algorithm performs better than the TEBD algorithm. Moreover the tails of the OTOC are accurately obtained both by TDVP MPS and TEBD MPO. They are in very good agreement with exact results at short times, and appear to be converged in bond dimension even at longer times. However the growth and saturation regimes are not well captured by both methods.

Figures

Equations

$$|\psi ⟩=\sum _{s_1,...,s_N}{A}^{[1]s_1}{A}^{[2]s_2}...A^{[N]s_N}|s_1,s_2,...,s_N⟩$$

$$\begin{array}{rl}U(\mathrm{d}t)& ={\mathrm{e}}^{-\mathrm{i}(H_{\mathrm{odd}}+H_{\mathrm{even}})\mathrm{d}t}\\ & \approx {\mathrm{e}}^{-\mathrm{i}{H}_{\mathrm{odd}}\mathrm{d}t}{\mathrm{e}}^{-\mathrm{i}{H}_{\mathrm{even}}\mathrm{d}t}+O(\mathrm{d}{t}^2)\\ & =\prod _j{\mathrm{e}}^{-\mathrm{i}{h}_{2j,2j+1}\mathrm{d}t}\prod _j{\mathrm{e}}^{-\mathrm{i}{h}_{2j+1,2j+2}\mathrm{d}t}+O(\mathrm{d}{t}^2)\end{array}$$

$$\sum _i{\dot{A}}^{[i]s_i}(t){\partial }_{A^{[i]s_i}}|\psi (\{A^{[k]s_k}(t)\}⟩=-i{P}_{T}H|\psi (\{A^{[k]s_k}(t)\})⟩.$$

$${\hat{V}}_i(t)=\sum _{{\alpha }_1\dots {\alpha }_L\in \{0,1,2,3\}}{v}_{{\alpha }_1,\dots ,{\alpha }_L}{\hat{\sigma }}_{{\alpha }_1}\otimes {\hat{\sigma }}_{{\alpha }_2}\otimes \cdots \otimes {\hat{\sigma }}_{{\alpha }_L}$$

$$\left[{\hat{V}}_i(t),{\hat{\sigma }}_{\alpha }^j\right]=\sum _{{\alpha }_1\dots {\alpha }_L}{v}_{\mathbf{\alpha }}\left[{\hat{\sigma }}_{{\alpha }_1}\otimes {\hat{\sigma }}_{{\alpha }_2}\otimes \cdots \otimes {\hat{\sigma }}_{{\alpha }_L},{\hat{\sigma }}_{\alpha }^j\right].$$

$$C_{ij}(t)=\frac{1}{2\cdot Z}{\Vert [{\hat{W}}_j,{\hat{V}}_i(t)]\Vert }_F^2,$$

$$\underset{t\to \infty ,|i-j|>vt}{\lim }{C}_{ij}\exp \left[\mu (v)t\right]=0,$$

$$C_{ij}=1-\frac{1}{2\cdot Z}\mathrm{Tr}\left({\hat{V}}_i(t){\hat{W}}_j{\hat{V}}_i(t){\hat{W}}_j\right),$$

$$\frac{1}{Z}\mathrm{Tr}\left({\hat{V}}_i(t){\hat{W}}_j{\hat{V}}_i(t){\hat{W}}_j\right),$$

$$\begin{array}{rl}& \frac{1}{Z}\mathrm{Tr}\left({\hat{V}}_i(t){\hat{\sigma }}_z{\hat{V}}_i(t){\hat{\sigma }}_z\right)\approx \\ & \approx \frac{1}{Z\cdot n_{\text{states}}}\sum _{\mathbf{\sigma }}^{n_{\text{states}}}⟨\mathbf{\sigma }|{\hat{V}}_i(t){\hat{\sigma }}_j^z{\hat{V}}_i(t)|\mathbf{\sigma }⟩{\sigma }_j.\end{array}$$

$$H=\sum _{i=1}^{L-1}{J}_z{\hat{\sigma }}_i^z{\hat{\sigma }}_{i+1}^z+\sum _{i=1}^L\left(h_x{\hat{\sigma }}_i^x+h_z{\hat{\sigma }}_i^z\right)$$

Full text

I Introduction

The question of quantum thermalization in closed systems receives currently a considerable amount of attention. This interest is partly due to the importance of this problem for the understanding of the foundations of statistical physics Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016); Borgonovi et al. (2016) and to the experimental progress leading to increasingly well isolated experimental realizations of quantum many-body systems in ultracold atomic gases in optical lattices Bloch et al. (2008). In general, the unitary dynamics of isolated quantum systems precludes reaching a maximally mixed state if the system is initialized in a pure state. Nevertheless, generic isolated quantum many-body systems typically are found to reach a thermal state, i.e. local observables assume values consistent with the micro-canonical ensemble at long timesKaufman et al. (2016). Considering only a small subsystem of the total system, the usual notion of thermodynamic equilibrium is recovered, as the reduced density matrix becomes equal to the corresponding thermodynamic density matrix Garrison and Grover (2018); Luitz (2016), as the rest of the system serves as a heat bath. This behavior is rooted in the local structure of eigenstates of the many body Hamiltonian which manifests itself in the eigenstate thermalization hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016); Borgonovi et al. (2016), motivated by random matrix theory considerations. It provides a precise prediction for the matrix form of local operators in the eigenbasis of the Hamiltonian, and implies thermalizationRigol et al. (2008); Luitz and Bar Lev (2016); Roy et al. (2018). The mechanism for this thermalization process is the loss of local quantum information over time, which implies that the full wave function of the initial state cannot be reconstructed from local measurements at long times Nandkishore and Huse (2015). In consequence of this loss of local information, the system becomes increasingly entangled, until the state of a subsystem reaches a maximally mixed state consistent with global constraints Kim and Huse (2013); Luitz et al. (2016); Kaufman et al. (2016). A direct local probe of the loss of local quantum information can be constructed by studying the spreading of initially local Heisenberg operators $\hat{O}_{i}(t)$ which become increasingly nonlocal over the course of time. The locality can be quantified by probing the real space support of $\hat{O}_{i}(t)$ using the norm of the commutator with another local operator $\hat{V}_{j}$. This quantity is now best known as the out-of-time-order correlator (OTOC) which was introduced to study quantum chaos Larkin and Ovchinnikov (1969); Maldacena et al. (2016), and to bound the spreading of information in systems with short ranged interactions Lieb and Robinson (1972).

Certain universal properties of the OTOC can be well understood in random unitary circuits von Keyserlingk et al. (2018); Rakovszky et al. (2018); Khemani et al. (2018); Nahum et al. (2018) where it is governed by hydrodynamic equations of motion. In these systems, a light cone structure was identified with a broadening front arising from the diffusive nature of the hydrodynamic equations. This diffusive behaviour has been found in numerically exact calculations in a noisy spin system Knap (2018). However, no exponential regime with a fixed Lyapunov exponent was found so far in such systems, nor in Hamiltonian systems with a small local Hilbert space and continuous timeLuitz and Bar Lev (2017).

While the OTOC is a powerful and universal theoretical tool, it is very difficult to calculate in practice for generic quantum many-body systems, due to its operator nature (see paragraph III.1). In the last two years several numerical methods for calculating the OTOC emerged: exact operator evolution in the Heisenberg picture Chen et al. (2016), matrix product operators (MPO) evolution in the Heisenberg picture Bohrdt et al. (2017); Xu and Swingle (2018) and an exact wave function technique in the Schrödinger picture Luitz and Bar Lev (2017). In this article, we will carefully compare these techniques and add two more MPS methods based on a stochastic sampling of the OTOC in the Schrödinger picture using both time evolving block decimation (TEBD) Vidal (2003) and the time dependent variational principle (TDVP) using matrix product states (MPS) Haegeman et al. (2011, 2016), which is currently discussed as a candidate method to extract late time hydrodynamic properties of quantum systems Leviatan et al. (2017); Kloss et al. (2018). While exact Schrödinger evolution using quantum typicality is currently the best choice to obtain the exact OTOC for Hilbert space dimensions of up to $10^{9}$ even at late times Luitz and Bar Lev (2017), MPS techniques have been recently presented as complementary approaches. In particular, MPO time evolution using TEBD can be used to extract the tails of the OTOC at very long distances (see paragraph III.1). Here, we investigate, how MPS based time evolution techniques such as TEBD and TDVP in the Schrödinger picture compare to the method of MPO evolution.

This article is structured as follows: in section II we review the different MPS methods that we used to simulate quantum dynamics. In section III, we present different ways of quantifying the spreading of operators in non-integrable systems, and the numerical approaches we choose to simulate them. Next, in section IV we compare the results obtained by the different methods, both for small and larger systems, and assess to which extent our results can be trusted despite the low amount of entanglement included in our MPS approximations.

II Matrix product states methods for the simulation of quantum time evolution

While the MPS formalism originally arose in the context of ground state physics, it has also been very successful in the description of quantum dynamics of one dimensional quantum systems Schollwöck (2011). Every quantum many body state can be brought into a MPS form, that is:

[TABLE]

where $A^{[i]s_{i}}\in\mathbb{C}^{\chi\times\chi}$ is the matrix corresponding to site $i$ and to the local state $|s_{i}\rangle$. Note that $A^{[1]s_{1}}$ ($A^{[N]s_{N}}$) are row (column) vectors of size $\chi$ to make the wavefunction coefficients scalar. The required bond dimension $\chi$ depends on the amount of bipartite entanglement contained in $|\psi\rangle$. In this section, we review different MPS based method to perform quantum time evolution.

II.1 Time evolving block decimation

The TEBD algorithm Vidal (2003); Schollwöck (2011) is a powerful and simple tool to simulate the short time dynamics of one-dimensional systems with short range interactions. It is applicable to both the Schrödinger and Heisenberg picture.

The key idea of the algorithm is to employ a Trotter decomposition of the time evolution operator in such a way that only local time evolution operators occur. For example a nearest-neighbour Hamiltonian is written as a sum over terms which involve only two neighboring sites: $\hat{H}=\sum_{i}h_{i,i+1}$ where $h_{i,i+1}$ is the part of $\hat{H}$ containing the interaction between sites $i$ and $i+1$. Since only nearest neighbor terms do not commute, we decompose the Hamiltonian into a part acting on the even bonds and a part acting on the odd bonds as $\hat{H}=H_{\mathrm{even}}+H_{\mathrm{odd}}$ where $H_{\mathrm{even}}=\prod_{j}h_{2j,2j+1}$ and $H_{odd}=\prod_{j}h_{2j+1,2j+2}$. In this way all the terms contained in $H_{\mathrm{even}}$ (resp. $H_{\mathrm{odd}}$) commute with each other. We can now apply a Trotter decomposition at first order for simplicity, although the implementation of any order is possible within this scheme. We obtain for the time evolution operator:

[TABLE]

Each term of the products in the third line of equation (2) can be written as a unitary gate linking two adjacent sites. In order to evolve an MPS with a time step $\mathrm{d}t$ we apply a layer of gates as shown in Fig. 1. It is then possible to re-express the time evolved state as an MPS of higher bond-dimension.

After each application of a unitary gate, the MPS is optimized using the Schmidt decomposition. This is done efficiently by writing the state in a Schmidt basis, and truncating the smallest Schmidt values, whose contributions to the wave-function are least important. This approximation is only valid in a regime where the state is lowly entangled, limiting the use of the method to short times for ergodic systems such as the ones we are studying in this work. Moreover the truncation process renders the time-evolution non unitary and does not preserve the norm of the state, the energy and other conserved quantities.

The TEBD algorithm is straightforwardly extended to the MPO-time evolution as illustrated on Fig. 1b, by noting that in contrast with the MPS case the (adjoint) unitary time evolution must be applied on the lower and upper legs of the MPO.

All the results obtained with the TEBD algorithm presented in this article were performed with a second order Trotter decomposition scheme.

II.2 The time dependent variational principle using matrix-product states

The time dependent variational principle (TDVP) was first introduced by Dirac Dirac (1930) for a general variational manifold. The general idea is to project the Schrödinger equation on the variational manifold of interest in such a way that the wave function after an infinitesimal time step does not leave the manifold, yielding an approximate tractable time evolution.

It has been recently formulated using MPS with a fixed bond dimension as the variational manifold Haegeman et al. (2016, 2011) and is based on the concept of the tangent space of the MPS manifold Haegeman et al. (2013). The algorithm is very similar to the density matrix renormalization group (DMRG) Östlund and Rommer (1995); Dukelsky et al. (1998) and offers several advantages with respect to the TEBD algorithm since it does not rely on truncation to keep the wave function in the manifold $\mathcal{M_{MPS,\chi}}$ of MPS with a given dimension $\chi$. Moreover it is suitable to simulate Hamiltonians with long range interactions. To derive the TDVP algorithm, we again start from an MPS wave function Eq. (1) as a variational ansatz. However, we make every tensor explicitly time dependent, i.e. $A^{[i]s_{i}}\rightarrow A^{[i]s_{i}}(t)$.

After inserting this ansatz into the Schrödinger equation, we find:

[TABLE]

where $P_{T}$ is the projector on the tangent plane of $|\psi(\{A^{[k]s_{k}}(t)\})\rangle$. This projection is necessary to ensure that we obtain closed equations of motion in the tangent space so that the time evolved MPS is confined to $\mathcal{M_{MPS,\chi}}$. The new equations of motion obtained this way can be elegantly integrated by means of a splitting method. For more details, we refer to ref. Haegeman et al. (2016).

Unlike the TEBD algorithm, the TDVP algorithm conserves energy even when the exact time evolution cannot be captured by the MPS. Moreover all the conserved quantities, which do not cause an increase of bond dimension once applied to a MPS, are respected. In other words, a quantity $\hat{O}$ commuting with the Hamiltonian is conserved by TDVP if for all states $|\phi\rangle$ expressible as a MPS of bond dimension $\chi$ we can still express $\hat{O}|\phi\rangle$ as a MPS of bond dimension $\chi$. In the case where all conserved quantities leave the manifold $\mathcal{M_{MPS,\chi}}$ invariant, TDVP appears well suited for the simulation of thermalization, since it is believed that at long times the dynamics of the system is driven by hydrodynamical equations of motion governed by conserved quantities Glorioso and Liu (2018); Blake et al. (2018); Lucas (2017). Following this line of thought, TDVP has recently been applied successfully in the context of thermalization Leviatan et al. (2017). However the accuracy of results at long times is currently under debate Kloss et al. (2018). It has also been used advantageously in disordered systems Doggen et al. (2018). In this work, we focus on short to intermediate times and will not consider the question of the relevance of TDVP in the context of hydrodymamics. However the TDVP algorithm provides a unitary time evolution. This feature will be of crucial importance when calculating OTOCs, as explained in section III.1.

III Measures of operators spreading in closed quantum systems: the out-of-time-order correlator

III.1 Definition

In ergodic isolated quantum systems, local operators in the Heisenberg picture typically spread over the course of time in a sense that their supports in real space grows. Here we study this spreading in a non-integrable spin-$\frac{1}{2}$ chain of length $L$. In this case, the growth of the support means that the expansion of an operator $\hat{V}_{i}(t)$ in the operator basis of strings of local Pauli operators $\hat{\sigma}_{i}$, $i\in\{0,1,2,3\}$ ($\hat{\sigma}_{0}=\hat{1}$):

[TABLE]

acquires increasingly long strings of non-identity Pauli operators von Keyserlingk et al. (2018); Khemani et al. (2018); Nahum et al. (2018). The growing complexity stems from the increasing nonlocality of the time evolution operator $\hat{U}(t)$ and is also related to the growth of the operator entanglement entropy Pižorn and Prosen (2009); Zhou and Luitz (2017); Dubail (2017). In order to define the latter quantity, it is useful to consider the operator as a state living in a larger Hilbert space. More precisely, if we decompose an operator $\hat{O}$ in terms of the elements of an orthonormal basis $\left\{|\phi_{i}\rangle\right\}$ as $\hat{O}=\sum_{ij}O_{ij}|\phi_{i}\rangle\langle\,\phi_{j}\,|$, we can then associate to this operator the state $|\phi_{\hat{O}}\rangle=\sum_{ij}O_{ij}|\phi_{i}\rangle\otimes|\phi_{j}\rangle$. The bipartite entanglement entropy of $|\phi_{\hat{O}}\rangle$ is called the operator entanglement entropy (assuming proper normalization). Note that the calculation of the entanglement of MPOs is identical to the MPS case: one has to find the singular values of the bond $i$ of interest and as usual $S=-\sum_{j}s_{j}^{2}\ln(s_{j}^{2})$ where $s_{j}$ are the singular values and $S$ the entanglement.

The expansion in terms of Pauli string operators in Eq. (4) is a useful measure of operator spreading but is computationally impractical due to the arising of an exponentially large number of terms in the length of the chain. However, we note that in order to study the loss of locality of quantum information, the most interesting information is contained in how the length of Pauli strings grows over time. Therefore, we instead consider commutators of the operator $\hat{V}_{i}(t)$ with nontrivial ($\alpha\in\{x,y,z\}$) local Pauli operators (i.e. $\hat{\sigma}_{\alpha}^{i}=\hat{1}\otimes\dots\otimes\hat{1}\otimes\hat{\sigma}_{\alpha}\otimes\hat{1}\otimes\dots\otimes\hat{1}$):

[TABLE]

This commutator is zero when the local Pauli operators $\hat{\sigma}_{\alpha_{j}}$ on site j are identities for all strings in the expansion $\hat{V}_{i}(t)$. Generically, when the operator support of $\hat{V}_{i}(t)$ reaches site $j$, its expansion in terms of Pauli strings will contain terms not commuting with $\hat{\sigma}_{\alpha}^{j}$. This commutator is therefore quantifying the operator spreading. However it is also an operator and is therefore usually reduced to its norm $\left\|\left[\hat{V}_{i}(t),\hat{\sigma}_{\alpha}^{j}\right]\right\|$. For computational simplicity, a standard choice for the norm is the normalized Frobenius norm given by $\left\|\hat{A}\right\|_{F}^{2}=\frac{1}{\mathcal{N}}\mathrm{Tr}\hat{A}^{\dagger}\hat{A}$, leading to the definition

[TABLE]

where $Z$ is the dimension of the Hilbert space. We have chosen the normalization in order to ensure that $C_{ij}(t)\to 1$ at long times in the case where $\hat{W}_{j}$ and $\hat{V}_{i}$ are hermitian operators, which square to identity, such as Pauli operators.

This quantity was originally proposed by Larkin and Ovchinniokv Larkin and Ovchinnikov (1969) in the context of quantum chaos. They showed that in chaotic systems with a semiclassical limit this norm of the commutator is connected to a Lyapunov exponent of the system and therefore effectively quantifies its chaoticity. Using this quantity they discussed a quantum analogue of classical chaos since in the semi-classical limit it quantifies the sensibility of classical trajectories to their initial conditions for the choice $\hat{W}_{j}=\hat{p}$ and $\hat{V}_{i}=\hat{x}$. This can be understood more intuitively by observing that the OTOC measures the effect of an initial perturbation on the value at later times of an operator located at some distance Maldacena et al. (2016). On the other hand, recent numerical studies of quantum systems with a small local Hilbert space and for $31$ sites showed that there is no regime of exponential growth Luitz and Bar Lev (2017), a discrepancy to the semiclassical case Rammensee et al. (2018), which has yet to be fully understood.

The link between the OTOC and locality of the Hamiltonian was made a few years later by Lieb and RobinsonLieb and Robinson (1972). They realized that information in systems with short range interactions can only spread within a light-cone with only exponentially suppressed leaking. This is most effectively quantified by considering the spreading of initially local operators $\hat{V}_{i}(t)$ in the Heisenberg picture. More precisely:

[TABLE]

for velocities $v>v_{LR}$, where $v_{LR}$ is called the Lieb-Robinson velocity. The function $\mu(v)$ is now referred to as velocity dependent Lyapunov exponent Khemani et al. (2018).

The OTOC has been the subject of a renewed interest in the past few years due the establishment of a duality between some strongly correlated systems and black-holes and the proposal of exactly solvable models to illustrate it Kitaev (2015). Moreover the spreading of operators is directly connected to the scrambling of local quantum information, since in chaotic systems at long times, initially local operators lose their locality and become completely scrambled Maldacena et al. (2016).

III.2 Numerical considerations

If we restrict ourselves to hermitean unitary operators, which square to identity, such as Pauli operators, the OTOC $C_{ij}$ can be expressed as:

[TABLE]

where $Z=\mathrm{dim}(\mathcal{H})$ is the dimension of the Hilbert space (see paragraph III.1). The nontrivial part of the calculation of this quantity consists of determining the correlation function

[TABLE]

which is exponentially expensive: for a spin-$\frac{1}{2}$ system of size $L$, $\hat{V}_{i}(t)$ is represented by a matrix in $\mathbb{C}^{2^{L}\times 2^{L}}$. Direct exact time evolution of the operator will therefore be very limited in system size Chen et al. (2016).

Alternatively the trace can be stochastically evaluated with typical, randomly chosen wave functions due to quantum typicality. Although the time evolution will still be exponentially expensive, larger system sizes can be achieved since a state has only $2^{L}$ components. This has been achieved using exact Krylov space time evolution Luitz and Bar Lev (2017, 2018); Park and Light (1986); Manmana et al. (2005); Moler and Van Loan (2003). Here, instead of operators, only wave functions are evolved in time by moving to the Schrödinger picture at the price of performing the time evolution forward and backwards in time, yielding an overall scaling of the method proportional to $t_{\text{max}}^{2}$. This is so far the most powerful numerically exact method to simulate the ergodic dynamics of small to intermediate system sizes up to arbitrary times and used here as a benchmark. For details of the method, see Refs. Luitz and Bar Lev, 2017, 2018.

Another approach is to use Heisenberg propagation of a matrix product operator (MPO) representation Bohrdt et al. (2017); Xu and Swingle (2018) of $\hat{V}_{j}(t)$ (see Fig. 1). To calculate $C_{i,j}$, one can use equation (8). We first evolve $\hat{V}_{j}$ using the setup of Fig.1b) with $O^{j}=\hat{V}$ and all the other operators $O^{k}$, $k\neq j$ set to identity. This way we obtain an MPO which tensors we denote as $V_{j}^{l}(t)$, the index $l$ corresponding to the site of the tensor. We also write $W_{i}$ as a MPO (which means that we place the operator $\hat{W}$ on site $i$ and identities operators on every other site). The calculation of the trace is then performed according to Fig.2.

This approach can seem to be limited to short times due to the linear growth of the operator entanglement entropy Pižorn and Prosen (2009) implying the necessity of an exponentially large bond dimension for an exact representation of the Heisenberg operator.

It has been argued recently in Ref. Xu and Swingle (2018) that this method is able to capture the early growth of the OTOC even with low bond dimension based on the following observation: first the spreading of quantum information is bounded by a light cone, implying that the operator entanglement of bipartitions with a cut outside of the light cone is small, thus leading to a small required bond dimension. Therefore in the Heisenberg picture, only tensors inside the region where the entanglement is high will be truncated within the TEBD scheme. Finally, it is assumed that the effect of the truncation propagates as a light cone, meaning that the sites with low entanglement should not be affected immediately by the effect of a truncation far away from them. In numerical simulations, the convergence of the results with bond-dimension presented in Ref. Xu and Swingle (2018) seems to support this reasoning. Our benchmarks for small systems (Fig. 3), our comparison of the contour lines of the OTOC obtained using MPO evolution to other methods (Fig. 6), as well as our analysis of the convergence with bond dimension (Fig. 10) also provide further support that MPO time evolution does indeed accurately capture the tail of the OTOC. We would like to note that this MPO technique has been applied in the past to calculate operator spreading in the one-dimensional Bose-Hubbard model in Ref. Bohrdt et al. (2017), where a discrepancy to the ballistic spreading at early times for small bond dimensions was pointed out, which was attributed to the truncation of the bond dimension.

Here, we propose a scheme based on the Schrödinger picture to MPS time evolution methods (see section II). The trace in Eq. (9) is sampled stochastically over random product states $|\sigma_{1},\dots\sigma_{L}\rangle$, which is reminiscent of minimally entangled typical thermal states (METTS) at infinite temperatureWhite (2009) ($\beta=0$). Additionally, we have the freedom to chose the basis such that the basis states are eigenstates of the operator $\hat{W}_{j}$, which we take for convenience to be $\hat{W}_{j}=\hat{\sigma}_{i}^{z}$:

[TABLE]

This way, we only have to propagate one wave function (i.e. $|\boldsymbol{\sigma}\rangle$) forward in time, apply $\hat{V}_{i}$, and propagate back to $t=0$ for each initial state $|\boldsymbol{\sigma}\rangle$. The average is performed over $n_{\text{states}}$ initial states, which are sampled uniformly from the local $\sigma_{z}$ product state basis (subject to sector constraints if required). From now on, we we will restrict ourselves to the case $\hat{V}_{i}(t)=\hat{\sigma}_{i}^{z}(t)$.

In order to evaluate equation (10), we use both the TEBD and the single site TDVP algorithms. In ref. Haegeman et al. (2016), a two-site implementation of the TDVP algorithm was proposed. However, in our case, this version of the algorithm would not be suitable for our purposes since it also relies on truncation to keep the bound dimension of the MPS fixed, hence yielding a non unitary time evolution and violating conservation of energy. However the single site algorithm does not allow to increase dynamically the bond dimension. In order to address this problem, we initialize our MPS as follows. First, we fill up the MPS with zeros in such a way that the product state, initially of bond dimension one, acquires the desired bond dimension $\chi$. Second, we bring our inflated product state in canonical (isometric) form following the usual sweeping procedure Schollwöck (2011). This state is then a proper state and the single-site TDVP reproduces the correct time evolution.

IV Results

We study a one dimensional quantum spin chain with short range interactions which has both integrable and non-integrable points as a function of the field angle: the tilted field Ising model. The Hamiltonian of the system is given by:

[TABLE]

We consider this model at a strongly nonintegrable point, with no other conservation laws besides the global conservation of energy, which is exactly respected by our TDVP approach. Following Ref. Kim and Huse (2013), we use the following parameters throughout this article: $J_{z}=1$, $h_{x}=(\sqrt{5}+1)/4=0.8090$,$h_{z}=(\sqrt{5}+5)/8=0.9045$.

IV.1 Comparison of the methods with exact results

In order to explore the domain of validity of the different methods (MPS TDVP, MPS TEBD and MPO TEBD), we compare the results for the OTOC with exact results obtained by Krylov space based exact time evolution (ETE) for chains of $L=21$ spins.

In Fig. 3, we present a detailed comparison for a system of size $L=21$ of the OTOC $C_{4,j}=\frac{1}{2}\|\left[\hat{\sigma}_{4}^{z}(t),\sigma_{j}^{z}\right]\|_{F}^{2}$ obtained from the four methods compared in this article. All panels show the numerically exact result obtained from exact time evolution (ETE) as solid lines, panel a) shows the TDVP result for the OTOC obtained from stochastic sampling of the trace in Eq. (8) using 98 random product states, panel b) shows the same calculation but using TEBD time evolution instead. In panel c), we show TEBD MPO evolution results, using a direct evaluation of the trace. Here, all calculations where performed using a maximal bond dimension of $\chi=64$. It is clear that at short times all three methods reproduce the exact result since there is no significant truncation occurring. Interestingly, the TDVP results stay close to the exact result for longer times than the OTOC obtained by TEBD. Similarly to TDVP, the MPO evolution using TEBD captures very well the regime of low values of the OTOC. Nevertheless the growth and saturation regime is not correctly reproduced by any of the methods. With the Schrödinger approach the OTOC is systematically overestimated while it saturates to an unphysical value in the case of the Heisenberg approach.

In order to make these statements more precise, we investigate the error of the different methods by considering the deviation from the exact result for $L=21$. From Fig. 3, we see that the discrepancy from the exact result is the largest at long times and long distances, independently of the choice of the method. Therefore we illustrate the errors resulting from each of the three methods at the longest spatial distance in our system from the origin at $i=4$ by the distance to the exact result $|C^{\text{exact}}_{4,21}(t)-C^{\text{MPS/MPO}}_{4,21}(t)|$ in Fig. 4. Here $C^{\text{MPS/MPO}}_{4,21}(t)$ stands for the OTOC calculated using either the MPS based methods with TEBD or TDVP time evolution or by the direct MPO based approach. We have checked that similar results are obtained for other distances.

These results are explained as follows. While TEBD in the Schrödinger picture suffers from the propagation of truncation errors since the wave function has to be propagated back to $t=0$ after a measurement, TDVP profits from the preservation of unitarity by the method, leading to a significantly smaller error compared to TEBD MPS. As for the Heisenberg picture, the low value of operator entanglement at the front of the OTOC light cone allows the MPO TEBD approach to capture the low values of the OTOC as explained in paragraph III.2. The systematic underestimation of the OTOC saturation values is due to finite bond dimensions limiting the captured operator entanglement. We note that similar results have been obtained in disordered spin chains where it has been demonstrated that TDVP performs better than TEBD Paeckel et al. (2019).

IV.2 Large systems and range of validity of the approximation

So far, we have presented results for systems small enough such that we could still compare to numerically exact results obtained by ETE. In what follows, we investigate the performance of these MPS and MPO methods for larger systems. We present in Fig. 5 the results for the OTOC $C_{4,j}(t)=\|\left[\hat{\sigma}_{4}^{z}(t),\sigma_{j}^{z}\right]\|_{F}^{2}$ as a function of time $t$ and distance $j-4$ for a system of size $L=50$ sites with bond dimension $\chi=64$ for both MPO TEBD and MPS TDVP. Again, we chose the position of the spreading operator $\hat{\sigma}_{4}^{z}(t)$ on the left of our chain with open boundaries instead of on the center, since this allows for a better resolution of the tails of the right part of the OTOC as discussed in Ref. Luitz and Bar Lev (2018). For both methods, the time step $\mathrm{d}t$ has been decreased until convergence, and we found that a significantly smaller time step for TDVP of $\mathrm{d}t=0.005$ was required compared to the MPO TEBD time step of $\mathrm{d}t=0.01$, since the splitting methods of TDVP and TEBD differ. While in TEBD the exponential of the Hamiltonian is decomposed into two-site gates using a Trotter decomposition, in TDVP the update of every tensor requires the integration of coupled differential which are solved separately for every site $i$. In TEBD, only neighboring terms do not commute, while in TDVP the differential equations involving tensors $A^{[i]s_{i}}_{kl}(t)$ at site $i$ depend on the value of the tensors $A^{[j]s_{j}}_{kl}$ on all sites $j$. Therefore, a larger time step error can be expected in TDVP due to a more severe approximation in the splitting method.

For this reason, the dependence of the error on the time step is more important and must be checked carefully. Additionally, the stochastic sampling of the trace in Eq. (8) does not admit importance sampling and is therefore costly, practically limiting the bond dimensions considered here to $\chi=64$. We evaluate the convergence in bond dimension of our results for both methods in appendix A. We find that we achieve convergence for low values of the OTOC, which is consistent with the benchmarks shown in Fig. 3.

Here, we do not consider MPS TEBD results because of the inferior accuracy of this method already identified for smaller systems as discussed in the previous section.

The representation of the results in Fig. 5 from the two methods on the same colorscale illustrates the problem observed for smaller systems in Fig. 3 that MPO TEBD (right panel) underestimates the saturation value of the OTOC, in agreement with recent results of Ref. Paeckel et al. (2019), while TDVP (left panel) reproduces the correct long time saturation value close to $1$. Next, we consider contour lines (solid lines in Fig. 5) $t(\theta,j)$ of the OTOC obtained from numerical solutions of the equation $C_{4,j}(t)=\theta$ for various thresholds. For very low thresholds, these contours capture the behavior of the tail of the OTOC, where both methods yield consistent results even at long times. At larger thresholds, the obtained contours are strikingly different: while MPS TDVP yields approximately linear contour lines, close to a linear light cone, the results obtained with MPO TEBD deviate strongly and yield a significantly slower information spreading. Due to the problems identified for MPO TEBD closer to the saturation regime of the OTOC as discussed in IV.1, we attribute this behavior to the error caused by the insufficient amount of operator entanglement included in our MPO approximation. The approximately ballistically spreading information front obtained with our MPS TDVP approach appears to be a qualitative improvement in comparison with MPO time evolution where the speed of information propagation seems to be underestimated. However, although qualitatively interesting, these results should not be trusted quantitatively at high thresholds since they are not converged in bond dimension in this region of space-time.

The results displayed in Fig. 5 can only be compared qualitatively, therefore we proceed by extracting the contours $t(\theta,j)$ of the OTOC for various values of the threshold and plot the results from MPS TDVP and MPO TEBD in the same figure panel for a direct quantitative comparison. In addition to the MPS results for $L=50$ also the exact results for $L=21$ are shown in Fig. 6. This is an important comparison, since results in other systems demonstrate that for short enough times, the OTOC does essentially not show any finite size effects Luitz and Bar Lev (2017, 2018). The contours obtained with ETE and MPS/MPO methods for different system sizes are therefore expected to coincide for short times and the contours at low thresholds should not depend on system size.

For very small thresholds ($\theta=10^{-5}$ and $\theta=10^{-3}$), see Figs. 6a and 6b, the contours obtained with MPO TEBD and MPS TDVP indeed match the exact results, in accordance with results of Fig. 3, confirming our expectations. However, some differences start to appear at higher thresholds ($\theta=0.1$ and $\theta=0.2$), see Figs. 6a and 6b, which can be expected from our study of convergence in bond dimension (see appendix A). We note that our MPS TDVP seems to yield a contour slightly closer to the exact result.

For small thresholds, it was previously observed in generic spin systems that the contours of the OTOC assume a power law shape with exponents close to unity Luitz and Bar Lev (2017). Therefore, we attempt power law fits to our numerically exact contours from ETE, yielding excellent fits. The fits are shown as gray dashed lines in Fig. 6, and should be understood as an extrapolation of the shape of the light cone from the $L=21$ results. For small thresholds ($\theta=10^{-5}$), the MPS/MPO approaches reproduce the extrapolated contours with very high accuracy, confirming the power law fit from the smaller system size and consistency with the exact result. For $\theta=10^{-3}$, the two approximated approaches are still in quite good agreement with the fit of the ETE, although some differences arise at later times. At higher thresholds, the difference is even more significant, since already short times results do not agree. This confirms our overall observation that the MPS approaches considered here reproduce the tail of the OTOC with good accuracy, while the growth and saturation regimes are not well captured.

V Conclusion

We have compared different MPS approaches to study information scrambling in a generic spin chain based on both matrix-product states (MPS) and matrix-product operators (MPO) and compared the results to an unbiased and numerically exact technique (ETE). For the calculation of the out-of-time-order correlators (OTOCs) in the Schrödinger picture based on MPS, we have shown that the use of a unitary time evolution method (TDVP) yields a significant improvement over the non-unitary truncation used in the time evolving block decimation (TEBD) algorithm. Furthermore we found that both MPO TEBD and MPS TDVP reproduce the tail of the OTOCs even at long times for low enough thresholds, while the growth and saturation regime suffers from truncation errors. The obtained shape of the light cone in large systems at low thresholds is in quantitative agreement with the ETE results at short times for smaller system size. Moreover they also match the extrapolated exact result even at late time. For larger thresholds, closer to the information front, a discrepancy from exact results is observed, which we attribute to insufficient convergence of both MPO TEBD and MPS TDVP results with bond dimension. However, our TDVP MPS approach still yields a qualitatively correct ballistic propagation of information in contrast with the results obtained with MPO TEBD where after significant truncation the spreading of information appears to halt, making the result unphysical. We also note that the asymptotic saturation value of the OTOC is correctly reproduced in our TDVP MPS approach, while strong truncation effects in the MPO TEBD approach lead to a severe underestimation of the saturation value.

We conclude that both MPS techniques in the Heisenberg and Schrödinger picture yield consistent results for the tails of the OTOC and their performance is comparable. However, the MPS TDVP approach comes at the price of introducing a stochastic sampling of the OTOC using random product states making it computationally much more expensive. An interesting future direction would be to apply our wave function approach to calculate the OTOCs in many body localized system to evaluate whether the logarithmic growth of entanglement allow us to gather reliable results in a broader region of space-time and calculate the contours OTOCs at larger thresholds.

Acknowledgements.

We would like the thank Michael Knap, Tibor Rakovszky, Johannes Hauschild and Yevgeny Bar Lev for the stimulating discussions. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 747914 (QMBDyn). We acknowledge PRACE for awarding access to HLRS’s Hazel Hen computer based in Stuttgart, Germany under grant number 2016153659. FP acknowledges the support of the DFG Research Unit FOR 1807 through grants no. PO 1370/2- 1, TRR80, the Nanosystems Initiative Munich (NIM) by the German Excellence Initiative, the DFG under Germany’s Excellence Strategy– EXC-2111-390814868, and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 771537).

Appendix A Convergence with bond dimension

A.1 Convergence with bond dimension and comparison to exact results for TDVP MPS and TEBD MPS

In any MPS calculation, the convergence of the results with the bond dimension is important. Here, we present an analysis of the convergence of the OTOC results calculated by our stochastic TDVP and TEBD methods based on MPS. We also analyze the convergence with the bond dimension of our TEBD MPO results.

In Fig. 7, for a system of size $L=21$, we show the convergence of the OTOC $C_{4,j}(t)$ with bond dimension for the MPS time evolution methods that we compare in the main text (TDVP and TEBD). Since both methods rely on a stochastic sampling of the trace in Eq. (8), we eliminate the error induced by the stochastic sampling by selecting 5 random product states $|\psi_{k}\rangle$ and then calculating the approximate OTOC $C_{4,21}^{(5)}=1-\frac{1}{2\cdot 5}\sum_{k=1}^{5}\langle\,\psi_{k}\,|\hat{\sigma}_{4}^{z}(t)\hat{\sigma}_{j}^{z}\hat{\sigma}_{4}^{z}(t)\hat{\sigma}_{j}^{z}|\psi_{k}\rangle$ with MPS TDVP, MPS TEBD and ETE for different bond dimensions (always using the same 5 product states). We plot the error given by the difference to the exact result for these $5$ states ($|C_{4,21}^{(5),\mathrm{exact}}-C_{4,21}^{(5),\mathrm{MPS}}|$) for different bond dimensions between $\chi=16$ to $\chi=128$ in Fig. 7. For clarity, we indicate the number of random product states included in this comparison in parentheses (here by $(5)$). We observe a clear convergence of the results from MPS TDVP. For MPS TEBD the error also decreases with the bond dimension, but stays always much larger than the one of TDVP. This underlines the advantage of the conservation of unitarity by TDVP time evolution.

A.2 Convergence with bond dimension for larger systems in TDVP MPS

For larger system sizes, it is difficult to obtain exact results for the OTOC as a benchmark. Therefore, a careful analysis of the dependence of the results on the bond dimension is crucial to identify the domain of validity of the methods.

In the case of the TDVP MPS approach, we compare cuts $C_{4,j}^{(15)}(t)$ of the approximate OTOC using 15 random initial product states in the $\sigma_{z}$ basis for several fixed distances. Note the same initial states are chosen for both bond dimensions in order to eliminate the importance of statistical errors in this comparison as explained above. For converged results within this approach, the mean and the error (calculated using bootstrap sampling) obtained for different bond dimension results should perfectly agree. In Fig. 8, we show the approximate OTOC $C_{4,j}^{(15)}(t)$ for $\chi=32$ and $\chi=64$ together with the errorbars of the OTOC (shaded region for $\chi=64$ and errorbars for $\chi=32$), yielding very good agreement of the results for low thresholds. At larger thresholds, the discrepancy between the two bond dimension results becomes significant as expected.

We find that the convergence in bond dimension depends significantly on the initial state and therefore we repeat this analysis for approximate OTOCs $C_{4,j}^{(1)}(t)$ using only single random product states and different bond dimensions in Fig.9. From the 15 product states included in Fig. 8, we select the states with the best and worst convergence in bond dimensions to illustrate these state to state differences. Overall convergence is only achieved only for low values of the OTOCs, confirming the observation that the tail of the OTOC is reproduced accurately, while values at larger thresholds are not converged.

A.3 Convergence with bond dimension for larger systems in TEBD MPO

In the case of the MPO TEBD approach, the study of the convergence in bond dimension is facilitated by the absence of stochastic sampling. The spacio-temporal dependency of the effect of bond dimension can be analyzed by directly looking at extracted contour lines $t(\theta,j)$ of the OTOC obtained from numerical solutions of the equation $C_{4,j}(t)=\theta$ for various thresholds and bond dimensions. We present the result of this approach in Fig. 10. The contours obtained with different bond dimension coincide very well for small thresholds (up to $\theta=10^{-3}$). However a difference between $\chi=64$ and $\chi=128$ appears already for $\theta=10^{-2}$ at times $t=30$, which we attribute to the insufficient representation of the operator entanglement of the operator in an MPO with $\chi=64$. The difference is even more striking for $\theta=0.1$, where the results are not converged and do not show the expected asymptotic behavior of a linear light cone. The breakdown appears earlier for the lowest bond-dimension as expected.