A simple tensor network algorithm for two-dimensional steady states

Understanding dissipation in 2D quantum many-body systems is an open challenge which has proven remarkably difficult. Here we show how numerical simulations for this problem are possible by means of a tensor network algorithm that approximates steady states of 2D quantum lattice dissipative systems in the thermodynamic limit. Our method is based on the intuition that strong dissipation kills quantum entanglement before it gets too large to handle. We test its validity by simulating a dissipative quantum Ising model, relevant for dissipative systems of interacting Rydberg atoms, and benchmark our simulations with a variational algorithm based on product and correlated states. Our results support the existence of a first order transition in this model, with no bistable region. We also simulate a dissipative spin 1/2 XYZ model, showing that there is no re-entrance of the ferromagnetic phase. Our method enables the computation of steady states in 2D quantum lattice systems.

Projected Entangled-Pair Operators (PEPO) are simply the operator version of Projected Entangled-Pair States (PEPS), in the same way that Matrix Product Operator (MPO) are the operator version of Matrix Product States (MPS) for the 1d case [1][2][3][4][5]. More specifically, a 2d PEPO is an operator that acts on a 2d PEPS and produces a new PEPS, and admits a tensor network description as in Supplementary Fig.(1). In principle there is no restriction on the coefficients of the tensors, so that PEPOs can represent, at least a priori, operators of any kind: generic, unitary, positive, and so on. In our case we use PEPOs to describe reduced density matrices, which are positive by construction. However, a PEPO does not need to be necessarily positive, and therefore the negative eigenvalues need to be under control in order to produce an accurate representation of a physical mixed state, as explained in the main text. Moreover, one can "vectorize" the PEPO, so that the resulting object can be treated as a 2d PEPS with double physical indices. Importantly, in this supplementary material we add diagonal weight matrices λ at the links. This is convenient in order to implement the so-called simple update, which we comment in the next section. In practice we also used a 2-site unit cell with tensors A and B at every site, and diagonal positive matrices λ 1 , λ 2 , λ 3 and λ 4 at every link, as shown in Supplementary Fig.(1). FIG. 1: From operator to vector. PEPO on an infinite 2d lattice, with a 2-site unit cell, tensors A and B at sites, and λ1, λ2, λ3 and λ4 at the links, as well as its vectorization.

Supplementary Note 2: Simple update
The time evolution generated by the Liouvillian superoperator is broken via a Trotter decomposition into small two-body gates, namely where we implemented for concreteness the first-order Trotter approximation, and we defined the 2-body gates g [i,j] acting on the different links. The action of one of these gates in a given link is accounted for by defining new approximated PEPO tensors following the scheme in Supplementary Fig.(2), called simple update [1][2][3][4][5]. This is the direct generalization of the update rule of the TEBD algorithm for 1d MPS [6,7]. The update is locally optimal in 1d, whereas only approximate in 2d because it does not take into account the effect of the environment of the link in the approximation of the tensors. Still, it is remarkably efficient, and produces good results for gapped phases with small correlation length.

Supplementary Note 3: Local observables
The calculation of local observables follows from an approximate contraction of the 2d tensor network using corner transfer matrices (CTM) [8][9][10][11][12][13][14][15]. For instance, the calculation of the (unnormalized) 1-site density matrix is done as shown in Supplementary Fig.(3). First, square-roots of the λ tensors are contracted with the tensors at every site as in Supplementary Fig.(3(a)). The partial trace is taken as in Supplementary Fig.(3(b)), which produces a tensor network as in Supplementary Fig.(3(c)). This tensor network is approximated using four CTMs C 1 , C 2 , C 3 and C 4 , as well as four half-row/column transfer matrices T a u , T a r , T a d and T a l as in Supplementary Fig.(3(d)) -which would be eight for a 4 site unit cell -. These approximating tensors are the effective environment of the site where we compute the reduced density matrix.  In Supplementary Fig.(4) we provide the basic details of how the tensors for the effective environment are computed [8][9][10][11][12][13][14][15]. In particular, for a "left move", two rows are inserted in the network and absorbed towards the left. The growth of the bond index is renormalized by an isommetry W (see Supplementary Fig.(4(c)), which can be computed according to several prescriptions [8][9][10][11][12][13][14][15]. The procedure follows by iterating directional moves along the left, right, up and down directions until convergence.

Supplementary Note 4: Operator-entanglement entropy
We define the operator-entanglement entropy as In the above equations, |ρ is the vectorized reduced density matrix, and tr E is the partial trace over the sites for which we wish to compute the entropy. In a nutshell: this is the entanglement entropy of |ρ , the vectorized density matrix, understood as a pure state. As such, this is not a measure of entanglement of the mixed state ρ. However, this is the relevant measure of correlations for our purposes, since it is upper-bounded directly by the bond dimension of the PEPO. Namely, if the PEPO has bond dimension D, then for a block of L × L sites one has which means that we can use it to quantify how large needs to be our bond dimension D for the PEPO, being this directly connected to the computational cost and the accuracy of the method [16].
In what follows we explain to procedures to compute S op (ρ): one fully taking into account the environment of the block, and one approximate taking into account some of the properties of the simple update.

Full calculation
The calculation taking into account the full environment follows the tensor contraction from Supplementary Fig.(5). In this case, we compute σ by tracing out the degrees of freedom outside the block, as shown in the figure. The corresponding contraction can be approximated in the thermodynamic limit using the CTM method explained above.

Simple calculation
For the case of a phase with small correlation length, and using the information obtained from the simple update (namely, tensors at the sites and at the links), it is indeed possible to approximate, up to a good accuracy, the tensor network in Supplementary Fig.(5) by the one in Supplementary Fig.(6). In this approximation, one does not take the surrounding environment of the block fully into account. Instead, the effect of the environment is replaced by the effect of the λ tensors surrounding the block, which amounts to a mean-field approximation of the effective environment. Moreover, one can see that in such a case the eigenvalues of σ can be approximated with good accuracy by the product of the squares of the surrounding λ tensors, i.e., and therefore the operator entanglement entropy reads with This approximation works very well in gapped phases computed via the simple update, and it is the one that we used in the main text.