Stable iPEPO Tensor-Network Algorithm for Dynamics of Two-Dimensional Open Quantum Lattice Models

Being able to accurately describe the dynamics steady states of driven and/or dissipative but quantum correlated lattice models is of fundamental importance in many areas of science: from quantum information to biology. An efficient numerical simulation of large open systems in two spatial dimensions is a challenge. In this work, we develop a tensor network method, based on an infinite projected entangled pair operator ansatz, applicable directly in the thermodynamic limit. We incorporate techniques of finding optimal truncations of enlarged network bonds by optimizing an objective function appropriate for open systems. Comparisons with numerically exact calculations, both for the dynamics and the steady state, demonstrate the power of the method. In particular, we consider dissipative transverse quantum Ising, driven-dissipative hard-core boson, and dissipative anisotropic $XY$ models in non-mean-field limits, proving able to capture substantial entanglement in the presence of dissipation. Our method enables us to study regimes that are accessible to current experiments but lie well beyond the applicability of existing techniques.

Popular Summary

In the microscopic world governed by quantum mechanics, ensembles of many interacting particles tend to be very fragile—a small perturbation can have a large impact. When numerically simulating these open quantum systems, it is therefore important to include the effects of the environment. While methods to do this are well established for small systems composed of only a few parts and for larger 1D systems, the modeling of large 2D systems is particularly challenging. Here, we present a new algorithm to predict the behavior of large interacting quantum systems of spin-1/2 particles that interact with their environment and live on a 2D square grid.

The difficulty of simulating many-body quantum systems lies in the exponential size of the space of possible states that the quantum state can explore, particularly if long-range correlations and entanglement are present. “Tensor-network” methods address this problem by restricting how much of this space the system can visit. If a system interacts strongly with its environment, then long-range correlations and entanglement are typically curtailed, making the tensor network description accurate and efficient. In this work, we present a way for such methods to properly account for the effect of any distant spatial correlations on local dynamics, which we find is crucial for producing accurate results.

Changing the number of spatial dimensions that a system can explore is often accompanied by new and unexpected phenomena. Our work provides access to the physics of open 2D quantum lattice models, with exciting opportunities to discover new physics.

Figure 1

Open quantum lattice model of interacting spins. Nearest-neighbor spins are coupled via a hopping $J$ and interact with an external bath via a coherent drive $\mathrm{\Omega }$ and/or a dissipative process $\gamma$. The open system can be modeled by describing the unit cell and its environment using a tensor network.

Figure 2

Main steps in the time evolution algorithm. (a) Vectorized form of the iPEPO. (b) Contraction of $A_{j^{\prime }}$ and $A_{l^{\prime }}$ with the dynamical map $\mathcal{L}$. (c) SVD of $e^{\tau \mathcal{L}}{A}_j{A}_l$. (d) $D^{\prime }$ singular values of relative tolerance greater than ${\epsilon }_{D^{\prime }}$ retained in the diagonal bond matrix ${\sigma }^{\prime }$. (e) Isometries $\tilde{u}$ and $\tilde{v}$, which truncate the enlarged bond from $D^{\prime }$ to $D$, giving the new bond matrix $D$. (f) Updated tensors ${\tilde{A}}_j$ and ${\tilde{A}}_l$. (g) Full environment truncation algorithm to find the isometries $\tilde{u}$ and $\tilde{v}$, which maximize the fidelity between truncated and untruncated bonds.

Figure 3

Environment of the unit cell. Panel (a) shows the tensor $a_j^{\mathrm{tr}}={\mathrm{tr}}_d(A_j)$, which is found by tracing over the physical dimension, while in panel (c), $a_j^{\mathrm{HS}}=\mathrm{tr}{A}_j^{†}{A}_j$ is given by the inner product; in each case, the bond matrices (b) are split in two. (d) Trace effective environment ${\mathcal{E}}_j^{\mathrm{tr}}$ of the iPEPO tensor $A_j$ used to calculate the $d\times d$ reduced density matrix ${\rho }_j$. (e) Hilbert-Schmidt effective environment ${\mathcal{E}}_j^{\mathrm{HS}}$ of the tensor $a_j^{\mathrm{HS}}$ used in constructing the bond environment.

Figure 4

Dynamics of the dissipative Ising model for $h_x/\gamma =0$ in (a) strong ($V/\gamma =0.2$), (b) moderate ($V/\gamma =1.2$), and (c) weak ($V/\gamma =4.0$) spin-damping regimes calculated using $\mathrm{FET}+\mathrm{WTG}$ for a range of bond dimensions $D$ and superimposed with results calculated using the EXACT method. (d) Results for a regime not applicable to the EXACT method ($V/\gamma =0.5$ and $h_x/\gamma =1.0$) but that can be treated with $\mathrm{FET}+\mathrm{WTG}$. In each case, we plot (i) the magnetization $m^x(t)$, (ii) the average purity ${\mathrm{\Pi }}_1$ of the single-site reduced density matrices, (iii) the nearest-neighbor $S_{12}^{xx}$, (iv) the next-nearest-neighbor $S_{13}^{xx}$ spin-spin correlations, and (v) the average infidelity of truncation $\mathcal{I}(t)$ at each time step.

Figure 5

Comparison between $\mathrm{FET}+\mathrm{WTG}$ (solid lines), SU (dotted lines), and EXACT (dashed lines) in the moderately damped regime of the dissipative Ising model $V/\gamma =1.2$, $h^x/\gamma =0.0$. (a) Trace distance $T_2(t)$ as a function of time and at $t\gamma =10$ (inset) for a range of bond dimensions. (b) Magnetization $m^x(t)$ and (c) nearest-neighbor $S_{12}^{xx}$, showing that $\mathrm{FET}+\mathrm{WTG}$ outperforms SU.

Figure 6

(a) Accumulation of internal correlations in time quantified by the cycle entropy $S_{\text{cycle}}(t)$, which is more effectively curtailed by $\mathrm{FET}+\mathrm{WTG}$. (b) Infidelity of truncation $\mathcal{I}(t)$, which is an order of magnitude smaller than SU at each truncation step. The results form a moderately damped regime of the dissipative Ising model $V/\gamma =1.2$, $h^x/\gamma =0.0$ with $D=4$.

Figure 7

Fate of staggered-$XY$ phase at $J/\mathrm{\Gamma }=0.3$. (a,b) Local magnetizations $m^{x,y,z}(t)$ on the $A$ and $B$ sublattices as the state evolves from the $D=1$ steady-state solution for bond dimensions $D\in [3,4,5,6]$. (c,d) Representation of a $2\times 2$ plaquette of the lattice in (c) the staggered-$XY$ phase ($D=1$ steady state) and (d) the uniform phase ($D=6$ steady state). (e) Mean-field phase diagram with a transition at $J/\mathrm{\Gamma }=\frac{1}{4}$. Inset: radii $r=|\overrightarrow{k}-\overrightarrow{j}|$ of an odd (red) and even (blue) number of steps on the lattice. (f) Correlation function $S_{j,k}^{xx}=⟨{\sigma }_j^x{\sigma }_k^x⟩$ versus distance $r$, which has a staggered form that is a remnant of the staggered-$XY$ phase. Correlations at odd step radii (red squares) are zero, and those at even step radii (blue circles) are finite and decay with $r$. Inset: exponential fit to even step correlations, giving a correlation length $\xi \approx 0.93$.

Figure 8

Tensor diagrams representing some of the steps involved in performing the left-move component of the CTMRG algorithm used to calculate the effective environment ${\mathcal{E}}^{\mathrm{HS}}$.

Figure 9

(a) Effective environment ${\mathcal{E}}_{j,l}^{\mathrm{HS}}$ of the tensors at neighboring sites $j$ and $l$. (b) Bond environment ${\mathrm{\Upsilon }}_{j,l}$, which is the contraction of ${\mathcal{E}}_{j,l}^{\mathrm{HS}}$ and the updated tensors $A_j^{\prime }$ and $A_l^{\prime }$ with enlarged bonds $\{D_j^{\prime }\}\ge D$. In panels (c)–(e), using ${\mathrm{\Upsilon }}_{j,l}$, the terms in the fidelity between the truncated ($\varphi$) and untruncated ($\rho$) density matrices are calculated by contracting with the isometries $u$, $v$ and the bond matrix $\sigma$.

Figure 10

Tensor diagrams representing some of the steps involved in finding the isometries $\tilde{u}$ and $\tilde{v}$ and the bond matrix $\tilde{\sigma }$, which maximize the fidelity between the truncated and untruncated bonds. In panel (a), the Rayleigh quotient in $R$ is proportional to ${\mathcal{F}}^2$. Panel (b) shows that $P$ is the contraction of the bond environment ${\epsilon }_{j,l}$ and the isometry $v$. In panel (c), we show that $R$ is the contraction of the bond matrix $\sigma$ and the isometry $u$. Panel (d) shows that $B$ is the contraction of ${\mathrm{\Upsilon }}_{jl}$ with the isometry $v$. In panel (e), the new (primed) isometries are found by singular value decomposition of the contraction of the maximal eigenvector $R_m$ and $v$.