On the stability of the infinite Projected Entangled Pair Operator ansatz for driven-dissipative 2D lattices
Dainius Kilda, Alberto Biella, Marco Schiró, Rosario Fazio, Jonathan Keeling
SSciPost Physics Submission
On the stability of the infinite Projected Entangled PairOperator ansatz for driven-dissipative 2D lattices
D. Kilda , A. Biella , M. Schir´o , R. Fazio , J. Keeling Division of Chemistry and Chemical Engineering, California Institute of Technology,Pasadena, CA, USA Universit´e Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France JEIP, USR 3573 CNRS, Coll`ege de France, PSL Research University, 11 Place MarcelinBerthelot, 75321 Paris Cedex 05, France Institut de Physique Th´eorique, Universit´e Paris Saclay, CNRS, CEA, 91191Gif-sur-Yvette, France Abdus Salam ICTP, Strada Costiera 11, I-34151 Trieste, Italy Dipartimento di Fisica, Universit`a di Napoli ’Federico II’, Monte S. Angelo, I-80126Napoli, Italy SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews, KY169SS, United Kingdom* [email protected] 10, 2021
Abstract
We present calculations of the time-evolution of the driven-dissipative XYZ modelusing the infinite Projected Entangled Pair Operator (iPEPO) method, intro-duced by [A. Kshetrimayum, H. Weimer and R. Or´us, Nat. Commun. 8, 1291(2017)]. We explore the conditions under which this approach reaches a steadystate. In particular, we study the conditions where apparently converged calcula-tions may become unstable with increasing bond dimension of the tensor-networkansatz. We discuss how more reliable results could be obtained.
Contents
A.1 Summary of the iPEPS algorithm 9A.1.1 Time evolution: simple update 10A.1.2 Contraction: corner transfer matrix 11A.2 Extension to iPEPO 151 a r X i v : . [ c ond - m a t . o t h e r] F e b ciPost Physics Submission A.3 Benchmarking with ground state calculations 15
References 17
Tensor network approaches have provided a route to efficient numerical simulations acrossa wide range of physical problems [1–4]. In one dimension, matrix product states (MPS)were originally introduced in the context of one-dimensional quantum ground states [5, 6].They have subsequently been extended to finite temperatures [7, 8], and to open quantumsystems and density-matrix evolution [9, 10]. Such methods have been very fruitful in explor-ing the nonequilibrium steady states (NESS) of driven-dissipative one-dimensional systems,using matrix product operators (MPO) [11–21]. These methods can also be extended beyondone dimension, either by mapping a finite two-dimensional lattice onto a one-dimensionalchain [22]—see Ref. [21] for a driven-dissipative implementation—or via the projected entan-gled pair state (PEPS) algorithm [23–26]. The PEPS approach represents the two-dimensionallattice directly as a tensor network [3, 4], and allows a direct simulation of an infinite (trans-lationally invariant) lattice (iPEPS).In a significant development, Kshetrimayum et al. [27] presented results adapting theiPEPS algorithm to simulate open quantum systems on infinite 2D lattices. By using aninfinite projected entangled pair operator (iPEPO) algorithm, they calculated the NESS ofthe dissipative XYZ and transverse field Ising models — i.e. finding the steady state of amany-body Lindblad master equation. The ability to routinely apply such methods to two-dimensional open quantum systems is potentially very powerful. While one-dimensional sys-tems have shown a rich variety of collective behaviour, symmetry-breaking phase transitionsgenerally do not occur in open one-dimensional systems, while they can in two dimensions.Several alternative approaches to approximately simulate two-dimensional open systems havebeen proposed, including cluster mean field theory [18], corner space renormalization [28, 29],and neural network states [30–34]. However, so far, these methods have generally been re-stricted to small systems (or small clusters), making it challenging to extract critical behavior.As such, the ability to routinely use iPEPO could be extremely powerful to numerically ex-plore phase transitions and critical behavior in driven dissipative systems.Here, we explore in detail the stability of the iPEPO algorithm, introduced by Kshetri-mayum et al [27]. We find that while at short times the algorithm shows reasonable timeevolution, the behaviour at long times varies. In particular, we find that the algorithm onlyreaches a steady state in some parameter regimes, and close to dissipative critical points [35]it can fail to reach a steady state. The regimes where we fail to find a steady state corre-spond closely to the regimes where Ref. [27] found a larger value of their parameter ∆, whichmeasures how close the state they find is to a steady state. Moreover, we find that for someparameters, increasing bond dimension of the iPEPO representation does not systematicallyimprove the accuracy of the results. On the contrary, it can in some cases destabilize a fixedpoint obtained at a lower bond dimension. We also suggest some possible alternatives to thesimple-update iPEPO algorithm, which could help to alleviate the problem.2 ciPost Physics Submission
Our paper is organised as follows. In Sec. 2 we apply the iPEPO algorithm to calculateNESS of the dissipative XYZ model in 2D, and analyse whether a steady state can be found.Section. 3 concludes with some comments on alternative tensor network approaches for com-puting NESS in 2D. We also provide an extended appendix, Sec. A, which summarises ourimplementation of the iPEPO algorithm; this implementation can be found at [36]. Since thecore of the iPEPO and iPEPS algorithms is similar, we also present results benchmarkingour implementation against the prototypical applications of iPEPS: the ground states of thetransverse field Ising model and the hardcore Bose-Hubbard model.
In this section, using the iPEPO implementation described and benchmarked in Appendix A,we discuss finding the NESS of the dissipative spin-1 / ∂ t ρ = − i [ H XY Z , ρ ] + κ (cid:88) j (cid:16) σ − j ρσ + j − σ + j σ − j ρ − ρ σ + j σ − j (cid:17) , (1) H XY Z = (cid:88) (cid:104) i,j (cid:105) (cid:16) J x σ xi σ xj + J y σ yi σ yj + J z σ zi σ zj (cid:17) , (2)where σ x,y,zj are Pauli matrices at lattice site j , σ ± = ( σ x ± iσ y ), and J x,y,z are spin-spincoupling constants. J y M x (a) J y M x (b) Figure 1: The magnetization order parameter M x of the dissipative XYZ model as a functionof coupling strength J y , for J x = 0 . J z = 1 and D = 4. Energies given in units of κ = 1. Thered highlighted areas indicate parameter regimes where the iPEPO algorithm fails to reach asteady state. (a) Results computed using timesteps δt = 10 − , − run until a steady stateis found. (In the red regions, the run is stopped after N = 1000 steps with δt = 10 − followedby N = 2000 steps with δt = 10 − .) (b) Results calculated using a large timestep δt = 10 − and stopping after N = 1000 steps.We begin by computing the time evolution of the dissipative XYZ model for the sameparameters considered by Ref. [27]. Figure 1(a) shows magnetization averaged over the twosites i = A, B of iPEPO unit cell, M x = ( | (cid:104) σ xi = A (cid:105) | + | (cid:104) σ xi = B (cid:105) | ), as a function of couplingstrength J y using iPEPO bond dimension D = 4. We find that iPEPO algorithm onlyreaches a steady state for some values of J y , while in the red highlighted areas no steady3 ciPost Physics Submission state is found—the results continue to change with time. Where a steady state is found,our results closely match Ref. [27]. The red regions in our figure—where no steady state isreached—correspond to points where the Kshetrimayum et al [27] report a large error in theirsteady state result. As observed in [27], these regions occur near the critical points, where onecan expect correlation lengths to diverge. Figure 1(b) shows that if we use a large timestep, δt = 10 − , and deliberately stop the simulation early — i.e. after N = 1000 steps — one canreproduce results similar to those presented by Kshetrimayum et al [27] in the red region.However, because the simulation was stopped artificially early, the results in Fig 1(b) do notcorrespond to an actual steady state, and also inevitably contain significant Trotter errorsdue to the large timestep size. As we discuss further below, the failure to reach a steadystate that we observe occurs specifically in the SU time evolution. That is, it is completelyunaffected by the corner transfer matrix (CTM) contractions needed to compute observables.As a result, none of the results in the rest of this paper depend on the CTM contraction orenvironment bond dimension.We have encountered similar issues of failing to reach a steady state in other parameterregimes of the dissipative XYZ model, as well as for other systems such as the dissipativetransverse field Ising model in 2D. This raises important question about the practicality ofthe iPEPO algorithm as a tool to find the NESS of open quantum systems. In the following,to keep our discussion concise, we will restrict our attention to the dissipative XYZ model.Our goal below will be to understand when the iPEPO algorithm does and does not reach asteady state, focusing entirely on the SU time evolution, by studying the relative change ofsingular values. To measure this, we define (cid:15) Λ = | Λ n − Λ n − | max δt | Λ n | max , (3)in terms of the set of singular values Λ n at timestep n . As such, (cid:15) Λ is a measure of the largest change of singular value, rescaled for ease of comparison between different timesteps. step (a) J y = 1.2J y = 1.5 step (b) step (c) Figure 2: The evolution of (cid:15) Λ (from Λ [ U ] as defined in Eq. (3)) at J y = 1 . J y = 1 . δt = 10 − ,(b) δt = 10 − , (c) δt = 10 − . Energies given in units of κ = 1.In Figs. 2(a-c) we plot (cid:15) Λ against simulation step, for timestep sizes δt = 10 − , − , − respectively. We show both J y = 1 .
2, in the parameter regime where iPEPO fails to reach asteady state (blue line), and J y = 1 .
5, where a steady state is found (green line). At J y = 1 . (cid:15) Λ quickly decreases, indicating that we approach a steady state.4 ciPost Physics Submission However, at J y = 1 . (cid:15) Λ undergoes noisy oscillations throughout the time evolution,for all timestep sizes, never approaching zero. step IC-1IC-2 IC-3IC-4 step (a) step (b) step (c) Figure 3: The evolution of (cid:15) Λ at J y = 1 . (cid:104) σ x (cid:105) = 1, (cid:104) σ y (cid:105) = 1, (cid:104) σ z (cid:105) = − δt = 10 − , (b) δt = 10 − , (c) δt = 10 − . All other parameters are the same as in Fig. 1.Energies given in units of κ = 1.We next explore if using different initial conditions affects whether a steady state is found.Figures 3(a-c) show that the evolution of (cid:15) Λ at J y = 1 . (cid:104) σ x (cid:105) = 1, (cid:104) σ y (cid:105) = 1, (cid:104) σ z (cid:105) = − step J y = 1.39J y = 1.35J y = 1.34J y = 1.33 J y = 1.32J y = 1.31J y = 1.25J y = 1.2 step (a) step (b) Figure 4: The evolution of (cid:15) Λ at selected values of J y during the adiabatic parameter sweepfrom J y = 1 . J y = 1 . J y = 0 .
01 and usingtimesteps (a) δt = 10 − , (b) δt = 10 − . All other parameters are the same as in Fig. 1.Energies given in units of κ = 1.Another possible way to choose initial conditions for the problematic parameter regime isan adiabatic parameter sweep. We first calculate the NESS for a value of J y where iPEPO does5 ciPost Physics Submission reach a steady state, then change J y in small steps, using the steady state at each value of J y as the initial state for the next value. This strategy can bypass highly entangled intermediatestates where a very high D may be needed. In our case, we start at J y = 1 .
4, and graduallyreduce J y in steps of ∆ J y = 0 .
01 to J y = 1 .
2. Figures 4(a,b) show the evolution of (cid:15) Λ atselected values of J y during the parameter sweep, with timesteps δt = 10 − , − respectively.We observe that for J y ≥ . (cid:15) Λ shows a decreasing trend, indicating that iPEPO finds asteady state, while for J y ≤ .
32 we again find noisy oscillations. Smaller timesteps δt andsmaller sweeping steps ∆ J y (not shown) lead to the same conclusion. step = 1.0= 3.0= 5.0= 5.1 = 5.2= 5.3= 6.0= 7.0 step (a) step (b) Figure 5: The evolution of (cid:15) Λ at J y = 1 . κ during the adiabaticparameter sweep from a strong dissipation regime with κ = 8 to a weak dissipation regimewith κ = 1 in steps of ∆ κ = 0 . δt = 10 − , (b) δt = 10 − . All otherparameters are the same as in Fig. 1. Energies given in units of κ = 1.A similar strategy is to start from a strong dissipation regime (i.e. large κ ) where weknow the steady state is approximately factorizable, and then perform an adiabatic sweep tolower κ . Figures 5(a,b) show the time evolution of (cid:15) Λ at J y = 1 . κ , ina sweep starting from κ = 8 reducing κ in steps of ∆ κ = 0 .
1, with timesteps δt = 10 − , − respectively. Similarly to the J y sweep, we find a steady state exists for κ ≥ .
2, but beyondthis we again find noisy oscillations. To summarise this section, for those points where asteady state is not found, this result appears to be robust to a variety of initial states andsimulation protocols.
Figure 6 presents the effect of changing iPEPO bond dimensions D. Panels (a-c) show thetime evolution of (cid:15) Λ at J y = 1 . δt = 10 − , − , − respectively. Each panelshows simulations performed using different bond dimensions 3 ≤ D ≤
6; no steady state isfound for any of these values of D . Panels (d-f) show the same quantities but for J y = 1 . D = 4. In this case, notably, while iPEPO reaches asteady state for D = 3 ,
4, no steady state is found for D = 5 ,
6. To explore this further, Panels(g-i) show the behavior at J y = 1 . ≤ D ≤
15. We observe thatfor D = 12 a steady state is found for timesteps δt = 10 − , − . However, increasing thebond dimension further to D = 14 ,
15 leads again to noisy oscillations. These results suggestthat while larger bond dimension may eventually yield a meaningful steady state, spurious6 ciPost Physics Submission step D = 3D = 4 D = 5D = 6 step (a) step (b) step (c) step (d) step (e) step (f) step D = 10D = 12 D = 14D = 15 step (g) step (h) step (i) Figure 6: Evolution of (cid:15) Λ with bond dimension at J y = 1 . J y = 1 . δt = 10 − , middle (b,e,h) δt = 10 − , and right (c,f,i) δt = 10 − . Panels (a-f) show D = 3 , , ,
6, while panels (g-i) show D = 10 , , ,
15. Allother parameters are the same as in Fig. 1. Energies given in units of κ = 1.7 ciPost Physics Submission steady states can arise at small bond dimension which then change as the bond dimensionincreases further. In addition, we note that while we can run the SU time evolution for D = 15, the CTM calculations for this bond dimensions would require over 128 GB of RAM,making it challenging to find observables without a support of high performance distributed-memory calculations and quantum symmetries. As noted above, the results shown in Fig. 6do not depend on any CTM contraction, so this issue does not arise within the calculationswe present. From the results of the last section, we may conclude that the SU iPEPO algorithm at lowbond dimensions is not always stable, reaching a steady state only in some parameter regimes,typically away from dissipative critical points. In other regimes, the algorithm failed to reacha steady state for all bond dimensions D that we could access. Moreover, in some cases, evenwhen a steady state is found for a given value of D , this may change as the bond dimensionis increased, switching instead to noisy time-dependent dynamics. While we believe thatthere exists a value of D allowing for a faithful representation of the steady-state densitymatrix (when spatial correlations decay exponentially), this study is unable to conclude what typical value of the bond dimension is needed for a prototypical driven-dissipative latticemodel. Overall, the results shown here suggest that a significant caution is required whenextending the SU iPEPS algorithm in 2D to Liouvillian evolution. Below we discuss a numberof alternative approaches which may be employed instead.One alternative approach is to adapt the Full Update (FU) iPEPS algorithm [3, 24, 26, 37]to the Lindblad time evolution of mixed states. For closed systems, the FU scheme [3,24, 38] achieves an optimal truncation by using a variational update scheme that computesthe full environment at every step. Since the Liouvillian evolution involves non-Hermitianoperators, an issue here is to find a reliable non-Hermitian algorithm that could substitute thealternating least-squares scheme used in the two-site variational minimization in the standardFU algorithm. There are, however, recent works on time evolution in closed systems [39–42] forwhich FU presents problems with stability, meaning SU can be more accurate. However, veryrecent work by McKeever and Szymanska [43] has shown that a variation on full update—fullenvironment truncation—can indeed improve the stability of iPEPO.A closely related idea is to consider a global variational search algorithm that targets thenull eigenstate | ρ (cid:105) of either Liouvillian L or a Hermitian positive semidefinite object L † L ;such approaches were successfully used in one dimension [15,16], and extension of these meth-ods to iPEPO has been discussed in Ref. [44]. Solving the variational problem with L † L isparticularly appealing since it allows reusing the standard and robust Hermitian optimiza-tion algorithms. However, even when L = (cid:80) l L l,l +1 contains only nearest-neighbour terms,the product L † L = (cid:80) l,r L † l,l +1 L l + r,l + r +1 will introduce highly nonlocal terms. While this ismanageable in 1D [15], for 2D the nonlocal couplings may easily lead to unfeasibly large bonddimensions. This may perhaps be adressed by truncating the range of these nonlocal terms ashas been discussed in 1D [45]. In addition, variational optimization iPEPO approaches wouldrequire computationally expensive tensor contractions involving both the iPEPS representing | ρ (cid:105) and the iPEPO representing either L or L † L .A more promising approach, also discussed in Ref. [44] may be to extend novel variational8 ciPost Physics Submission iPEPS techniques for ground state calculations in 2D introduced in Refs. [46, 47], optimizingiPEPS tensors using tangent space methods or by solving a local generalized eigenvalue prob-lem. Notably, both approaches avoid the need to construct a full PEPO for the Hamiltonian.Adapting these algorithms to either L or L † L could dramatically reduce the computationalcosts that limit the practical use of variational iPEPS methods. The global variational opti-mization could also offer a potentially much more robust way of finding the NESS of L thanone could hope to achieve with the standard iPEPS algorithm relying on two-body updates. Acknowledgements
We acknowledge helpful discussions with Marzena Szyma´nska, Conor McKeever, and AndrewDaley. We are grateful for comments from Augustine Kshetrimayum on an earlier version ofthis manuscript.
Author contributions
The calculations presented here were performed by DK, using codedeveloped by DK with contributions from AB. The project was initially conceived by JK andRF. All authors contributed to the writing of the manuscript.
Funding information
DK acknowledges support from the EPSRC Condensed Matter Cen-tre for Doctoral Training (EP/L015110/1). JK, RF and AB acknowledge the Kavli Institutefor Theoretical Physics, University of California, Santa Barbara (USA) for the hospitalityand support during the early stages of this work; this research was supported in part by theNational Science Foundation under Grant No. NSF PHY11-25915. AB acknowledges fundingby LabEx PALM (ANR-10-LABX-0039-PALM).
A Implementation of iPEPO algorithm
As noted above, iPEPO is a simple extension of the iPEPS algorithm to nonequilibrium steadystates of Lindblad superoperators. As with other extensions of tensor network methods [10]the main idea is to frame the density matrix evolution through superoperators applied tovectorized many body density matrices, which we denote | ρ (cid:105) (cid:93) . In this appendix we first brieflysummarize the iPEPS algorithm and then discuss its extension to open quantum systems.While these algorithms are standard and have been described in full by Or´us in [3], wesummarise them here to provide a self-contained description of the method that we apply tothe open quantum systems. A.1 Summary of the iPEPS algorithm
The basic idea behind PEPS is to parameterize the quantum state tensor Ψ k ,k ,...,k N by atwo-dimensional array of interconnected rank-5 tensors (see Fig. 7). Each individual tensorrepresents a single site of the quantum many-body system, with one vertical leg correspondingto the local Hilbert space of dimension d , and four in-plane legs corresponding to the bondsbetween different lattice sites. We denote the bond dimension of PEPS by D , which limitsthe amount of entanglement that can be captured by PEPS.9 ciPost Physics Submission For translationally invariant systems, one may use the infinite PEPS (iPEPS) ansatz [24]working directly in the thermodynamic limit. We can construct an iPEPS by choosing a unitcell and representing its sites with tensors. We will consider problems with a square unit cell.Since we use a Trotterised time evolution that propagates pairs of sites, we will need onlytwo on-site tensors A and B to define iPEPS. We next describe the two main ingredients ofthe iPEPS approach: the imaginary time propagation of iPEPS and the calculation of theenvironment needed to extract observables. Λ [U] Λ [D] Λ [R] Λ [U] Λ [D] Λ [L] Λ [R] Λ [L] Λ [L] Λ [D] Λ [R] Λ [U] Γ [A] Γ [B] Γ [B] Γ [A] L U RD
Figure 7: The iPEPS time evolution using Eq. (5) involves propagating four different bonds‘U’, ‘D’, ‘R’, and ‘L’, indicated by different colours. The SU algorithm uses Vidal form withΓ [ A,B ] site tensors and Λ [ U,R,D,L ] diagonal bond matrices to represent iPEPS with a two-siteunit cell. A.1.1 Time evolution: simple update
Time evolution can be performed by the Simple Update (SU) method, which follows essentiallythe same main steps as the imaginary time infinite Time Evolving Block Decimation (iTEBD)algorithm [48–51]. In two dimensions we perform Trotter decomposition by splitting ourHamiltonian into four terms H U , H D , H R and H L , describing respectively the ‘U’ (up), ‘D’(down), ‘R’ (right), and ‘L’ (left) bonds of the lattice: H = H U + H D + H R + H L (4)The first order Trotter decomposition of the time evolution operator U ( δτ ) = e − Hδτ thenreads: U ( δτ ) = e − δτH U e − δτH R e − δτH D e − δτH L + O ( δτ ) (5)where δτ is the imaginary timestep. Similarly to iTEBD in one dimension, in SU we representiPEPS using Vidal form: i.e. the iPEPS with a two-site unit cell is fully specified by two Γ [ A,B ] site tensors and four Λ [ U,R,D,L ] diagonal matrices that store the singular values of iPEPS bonds,as seen in Fig. 7. We denote the local Hilbert space dimension by d , and the bond dimensionby D . The SU then consists of the following steps:1. Absorb Λ [ R,D,L ] tensors on the external bonds into Γ [ A,B ] to obtain Q [ A,B ] .2. Decompose each of Q [ A,B ] into subtensors v A,B and X A , Y B using an exact SVD orQR/LQ decompositions. The original rank-5 tensors Γ [ A,B ] had the dimensionalityof dD giving rise to a large computational cost O ( d D ) of the update procedure.10 ciPost Physics Submission However, an update performed using the new rank-3 subtensors has a substantiallyreduced cost of O ( d D ) since the dimensions of v A,B are considerably smaller andequal to d × dD × D and d × D × dD respectively.3. Contract the two-body propagator e − δτH U with v A,B and Λ [ U ] to form θ tensor.4. Decompose θ tensor into ˜v A,B and ˜Λ [ U ] tensors using SVD. To prevent the bond di-mension of our tensors from growing indefinitely, we must truncate ˜v A,B and ˜Λ [ U ] byretaining D largest singular values and discarding the rest.5. Recover the updated rank-5 tensors ˜ Q [ A,B ] by contracting the rank-3 subtensors ˜v A,B with X A , Y B respectively.6. To restore Λ [ R,D,L ] on the external bonds, we divide each of Q [ A,B ] by Λ [ R ] , Λ D , andΛ [ L ] . This procedure brings iPEPS back to its original Vidal form, with updated tensors˜Γ [ A,B ] and ˜Λ [ U ] for each ‘U’ bond on the lattice.All other steps are the same as in the one dimensional iTEBD algorithm. To find theground state we propagate iPEPS for N imaginary timesteps δτ until a steady state is achievedwith respect to the spectrum of singular values in Λ [ U,R,D,L ] .The SU algorithm is both simple and very efficient, with computational cost O ( D d ) ofthe time evolution. However, SU is suboptimal since it employs local truncations withouttaking into account the full environment of a unit cell. For MPS, this issue can be resolvedrelatively easily by transforming the tensor network into a canonical form, which orthonor-malizes other bonds surrounding the bond being truncated. This solution is not possiblein 2D, since there is no known canonical form for PEPS. To achieve an optimal truncation,one must use a variational update scheme that computes the full environment at every step.This procedure, known as the Full Update (FU) [3, 24], is considerably more expensive andbears the computational cost of O ( N χ D + N χ D ) where N is the number of steps of theimaginary time evolution. In practice, SU has been applied extensively to various modelsand yields sufficiently accurate results for systems with large gaps and sufficiently short cor-relation lengths [52]. Due to its simplicity and efficiency, it allows shorter computation timesand significantly higher bond dimensions than FU, and thus remains popular. However, thesuboptimal truncation becomes an issue near quantum critical points when correlation lengthsbecome long, and in these cases FU should be used instead. A.1.2 Contraction: corner transfer matrix
For the iPEPS representation to be of practical use, we must be able to extract expecta-tion values from it. Unlike MPS where we could evaluate overlaps exactly at a polynomialcost, the exact contraction of two PEPS is an exponentially hard problem that scales as O ( e L ) with PEPS size L [3]. Fortunately, there exist various computational algorithms thatcan perform this contraction approximately with high precision. For infinite systems, thesemethods typically proceed by computing the approximate environment of an iPEPS unit cell.This effective environment consists of a small set of tensors that represent the infinite ten-sor network surrounding the unit cell. Possibly the most successful technique for computingiPEPS environments is the corner transfer matrix (CTM) method [25, 53, 54], which will bethe method we use. In this section we will explain the details of the CTM algorithm.11 ciPost Physics Submission Since observables for quantum states involve the overlap of two copies of the state, thestarting point for CTM is the contraction of two iPEPS, which produces an infinite 2D networkmade of reduced tensors a . Each reduced tensor a results from the contraction of (cid:2) M A (cid:3) † and M A iPEPS tensors by their physical indices, except at the sites where any operator is applied,leading to a different reduced tensor a O . Supposing the bond indices of M A had dimension D , the reduced tensors a now have bonds of dimension D .For simplicity of exposition, we will start by considering a one-site unit cell. However,methods based on Trotter decomposition into even and odd bonds modify the translationalinvariance from one-site to two-site. Therefore, in practice we always use the two-site versionof the CTM algorithm. Let us now subdivide this network into a 1 × a tensors, and its environment that contains the remaining infinite tensor network in whichthe unit cell in embedded (see Fig. 8). The key idea is to represent the environment by aset of four corner matrices { C , , , } , and four transfer tensors { T , , , } – these tensors areconnected by new virtual bond indices of size χ . Similarly to the bond dimension of MPSand PEPS, the environmental bond dimension χ is the parameter that controls the accuracyof the CTM approximation of environment. The goal of CTM algorithm is to obtain theenvironmental tensors by performing a series of coarse graining moves:1. Initialize the CTM tensors, e.g. using a random-number initialization, or a mean fieldenvironment with χ = 1.2. Perform four coarse graining moves in the left, right, up, and down directions. The leftmove involves the following steps, illustrated graphically in Fig. 8(a):3. Insertion: insert an extra column into the CTM network that contains the unit celltensor a , and the transfer tensors T , .4. Absorption: absorb the new column into the left side of the CTM environment bycontracting their respective tensors. This increases the environmental bond dimensionby χ → D χ : to prevent the bonds from growing indefinitely we must implement anappropriate truncation scheme.5. Renormalization: truncate the environmental tensors by inserting appropriate isometries ZZ + = I that reduce the bond dimensions D χ → χ by projecting onto a relevantsubspace, as shown in Fig. 8(b).6. Repeat steps (2-5) to let the CTM environment grow in all four directions until itconverges. Convergence is typically achieved when the eigenspectrum of each cornermatrix { C , , , } reaches the fixed point.The CTM algorithm for a two-site unit cell, containing two a and b tensors, follows thesame steps as the one-site algorithm outlined above. The environment is now specified by aset of four corner matrices { C , , , } , and eight transfer tensors { T a , , , , T b , , , } . The maindifference is that in the ’Insertion’ step we now insert two new columns instead of just one,as shown in Fig. 8(c). There are now two ’Absorption’ and ’Renormalization’ steps in thealgorithm: the absorption of each column is followed by renormalization to reduce the bonddimension D χ → χ . The two-site algorithm also needs two types of isometries ZZ + = I and W W + = I in the ’Renormalize’ step, to obtain renormalized transfer tensors ˜ T a,b and cornertensors ˜ C , in Fig. 8(d). 12 ciPost Physics Submission (c) (b) C ' C ' T ' Z + ZZ + Z C ∼ T ∼ C ∼ C T T T T C T a bba C T C T b a T ba C T T T T C T a bba C T C T b a T ba C T T T T a bba C T T b a T T a b T T b a T T a b T T b a a bb a (d) C ' C ' T ' Z + ZW + Z a T ' WZ + b T ∼ C ∼ C ∼ T ∼ ab (a) C T C C T C T T (3) C T C T T T C C T T a (4)(5) C ' C ' T C T C T T ' C T C C T C T T ∼∼∼ Figure 8: (a) The main steps (3-5) of the left move of the CTM algorithm for a one-siteunit cell containing tensor a . (b) The renormalization step (5) is done by inserting isometries ZZ + = I into the left edge of CTM network. (c) The insertion step of the left move of theCTM algorithm for a two-site unit cell containing tensors a , b . (d) The renormalization stepfor a two-site unit cell is done by inserting two types of isometries, ZZ + = I and W W + = I .13 ciPost Physics Submission (c) C T T T T C T a bba C T C T T ba (a) a bb a cut 1 Q B Q A R B R A SVD V A V B †† (b) R B R A R B R A -1-1 SVD R B R A V † U W + W Λ Λ C T T T T C T a bba C T C T T ba a bb a cut 2 Q B Q A T T T T a bb a ba baa b a b Figure 9: (a,b) The steps to compute
W, W + isometries for the bonds split by the ’cut-1’. (c)To compute Z, Z + isometries for the bonds split by the ’cut-2’, we perform a translationally-invariant shift by inserting two extra rows of tensors in the green box, and repeat the stepsin (a).Clearly, the crucial step of CTM algorithm is calculating the isometries. Several differentmethods exist, for instance the ones described in Refs. [25,37,55], which we have implementedand tested in the process of developing our iPEPS code. The prescription we have found towork best was the one in Ref. [37], which achieved a smoother convergence and a more efficientrepresentation of the environment than its predecessors. Figure 9 illustrates graphically thecalculation of Z, W isometries to be used in the ’Renormalize’ step of the left move. InFig. 9(a,b) we compute
W, W + isometries to be inserted into the bonds split by the ’cut-1’. Inthe first stage in Fig. 9(a), we contract the lower and upper parts of the network, producing thetensors Q A and Q B respectively. In the second stage, also in Fig. 9(a), we decompose Q A and Q B using an exact SVD to obtain the R A,B and V † A,B tensors. In the third stage in Fig. 9(b),we form the product I = R A (cid:2) R − A R − B (cid:3) R B , and decompose (cid:2) R − A R − B (cid:3) = U Λ V † using SVD,this time truncating to the χ dominant singular values. The R A,B matrices and the SVDmatrices are then combined in the symmetrized fashion I ≈ (cid:2) R A U Λ / (cid:3) (cid:2) Λ / V † R B (cid:3) = W W + to construct the isometries W = (cid:2) R A U Λ / (cid:3) and W + = (cid:2) Λ / V † R B (cid:3) . To obtain the Z, Z + isometries for the bonds split by the ’cut-2’, we perform a translationally-invariant shift byinserting two extra rows of tensors in the green box, as shown in Fig. 9(c). Contracting thelower and upper parts of the network then gives tensors Q A,B for the cut-2. We can nowcompute
Z, Z + repeating exactly the same steps as for the W, W + before. Once all isometries Z, Z + and W, W + are available, one can finally use them to carry out the renormalization inFig. 8(d).The CTM algorithm described above has a computational complexity of O ( χ D ciPost Physics Submission + χ D ). Once we have found a converged CTM environment, it can be used to computevarious observables. A.2 Extension to iPEPO
To extend the above method to an open quantum system, as noted above, one may firstrepresent the density matrix ρ as a Projected Entangled Pair Operator (PEPO), and thenreshape (vectorize) it into a PEPS | ρ (cid:105) (cid:93) by combining both physical indices at each site. Theproblem of computing NESS for an infinite 2D lattice thus becomes equivalent to the problemof finding the ground states of Hamiltonians using the SU iPEPS algorithm, by replacingimaginary time Hamiltonian propagation with the real time Liouvillian propagation. Todistinguish between the iPEPS representing wavefunctions and vectorized density operators,we will refer to the density matrix version as the iPEPO algorithm.The main steps of the iPEPO algorithm are the same as in Sec. A.1, except for twodifferences that we discuss next. The first difference is the propagator. The imaginary timetwo-body propagators U α ( δτ ) = e − δτH α with Hamiltonian H α for a bond α ∈ { U, R, D, L } are replaced by the real time two-body propagators U α ( δt ) = e − δt L α , where L α is the two-body Liouvillian for a bond α ∈ { U, R, D, L } and t is the real time. The second differenceis that observables are calculated using (cid:104) O (cid:105) = Tr [ Oρ ], instead of (cid:104) O (cid:105) = (cid:104) Ψ | O | Ψ (cid:105) . Similarly,the correct normalization Tr [ ρ ] = 1 in contrast to (cid:104) Ψ | Ψ (cid:105) = 1. As such, when extractingobservables from iPEPO, the reduced tensors that are contracted to find the environmentcome from tracing out local indices instead of computing inner products. We may then applythe CTM method from Sec. A.1.2 to observables.We have implemented our iPEPO code in Fortran [36], including the CTM algorithm,the functionality required for computing local observables and two-point correlators, the SUprocedure for both NESS and ground state calculations, as well as the FU procedure forground state calculations (for benchmarking purposes only). In our calculations we graduallydecrease the time step δt during the simulation to reduce the effects of Trotter error whilekeeping the computational cost low. To determine when the calculation has reached a steadystate we require that the spectrum of singular values contained in each diagonal bond matrixΛ ∈ (cid:8) Λ [ U,D,R,L ] (cid:9) in Fig. 7 stops changing within some accuracy (cid:15) . More specifically, we takethe largest difference between singular values in diagonal matrices Λ n and Λ n − , at timesteps n and n − | Λ n | max and by timestep size δt : (cid:15) Λ = | Λ n − Λ n − | max δt | Λ n | max . (3)For a steady state we expect (cid:15) Λ to approach zero (or more precisely, a value dependingon machine precision), as the eigenvalue spectrum should cease changing. To determinenumerically when to stop, we define a steady state as being reached when (cid:15) Λ < (cid:15) for eachΛ ∈ (cid:8) Λ [ U,D,R,L ] (cid:9) . A.3 Benchmarking with ground state calculations
Due to the similarity between the iPEPS and iPEPO algorithms, we may benchmark ourimplementation against the known numerical results from ground state calculations of twomodels: the transverse field Ising model, and the hardcore Bose-Hubbard model.15 ciPost Physics SubmissionTransverse-field Ising model
Our first test problem is the transverse field Ising modelon an infinite square lattice, H = − J (cid:88) (cid:104) i,j (cid:105) σ zi σ zj − g (cid:88) i σ xi . (6)Here σ x,zi are Pauli matrices at site i , J is the nearest-neighbour coupling between spins, and g is the transverse magnetic field along the x axis. The ground state of this model exhibits asecond order phase transition between a paramagnetic phase at large g , and a ferromagneticphase at small g ; the order parameter of this transition is the longitudinal magnetization M z = (cid:104) Ψ GS | σ z | Ψ GS (cid:105) . We show results in units where J = 1.Figure 10(a,b) shows the longitudinal magnetization M z and transverse magnetization M x as a function of transverse field g in the vicinity of phase transition, for different values ofiPEPS bond dimension D . Our implementation of iPEPS reproduces accurately both the SUand FU results reported in Refs. [24–26]. As expected, the SU and FU calculations match wellfar from the critical point where correlations are short-ranged, and the results converge fastwith increasing values of D . Near the critical point FU becomes considerably more accuratethan SU at a given bond dimension. That is, due to the diverging correlation length, a muchhigher value of D is needed for SU to achieve the same level of accuracy as FU with D = 2 , D = 3 predicts the critical point around g = 3 .
05, in good agreement with previous iPEPS results in Refs. [24–26] and Quantum MonteCarlo results in Ref. [56]. g M z (a) SU, D = 2SU, D = 3SU, D = 5FU, D = 2FU, D = 3 g M x (b) n , n n, D = 5n, D = 3n, D = 2n , D = 5n , D = 3n , D = 2 (c) Figure 10: Ground state results benchmarking our iPEPS implementation. Panels (a,b) show(a) the longitudinal magnetization M z and (b) transverse magnetization M x as a function oftransverse field g (in units of J ) in the vicinity of phase transition of the transverse field Isingmodel. These are computed for different iPEPS bond dimensions D using both SU and FUas indicated. Panel (c) shows number density of bosons n per lattice site and the condensatefraction n for the hardcore Bose-Hubbard model as a function of chemical potential µ (inunits of J ), computed for different iPEPS bond dimensions D using SU. Hardcore Bose-Hubbard model
Our second test case is the Bose-Hubbard model (BHM)on an infinite square lattice, H = − J (cid:88) (cid:104) i,j (cid:105) (cid:16) a † i a j + H.c. (cid:17) − µ (cid:88) i a † i a i , (7)where the occupations are restricted to 0 or 1 bosons on each lattice site. Here, a † i , a i arehardcore bosonic creation and annihilation operators at site i , satisfying the commutation16 ciPost Physics Submission relation (cid:104) a † i , a j (cid:105) = (1 − a † i a i ) δ ij . J is the hopping rate between adjacent sites, and µ is theonsite chemical potential. This model undergoes a second order phase transition betweenthe superfluid and the Mott insulator phases at the critical value of µ/J . As before, wepresent results in units where J = 1. Figure 10 shows the number density of bosons n = (cid:104) Ψ GS | a † a | Ψ GS (cid:105) and the condensate fraction n = | (cid:104) Ψ GS | a | Ψ GS (cid:105) | as a function of chemicalpotential µ , for different bond dimensions D of SU iPEPS. Again, our iPEPS calculationsreproduce accurately the ground state results of Refs. [26, 57, 58]. References [1] F. Verstraete, V. Murg and J. I. Cirac,
Matrix product states, projected entangled pairstates, and variational renormalization group methods for quantum spin systems , Adv.Phys. , 143 (2008), doi:10.1080/14789940801912366.[2] U. Schollw¨ock, The density-matrix renormalization group in the age of matrix productstates , Ann. Phys. , 96 (2011), doi:/10.1016/j.aop.2010.09.012.[3] R. Or´us,
A practical introduction to tensor networks: Matrix product states and projectedentangled pair states , Ann. Phys. , 117 (2014), doi:10.1016/j.aop.2014.06.013.[4] I. Cirac, D. Perez-Garcia, N. Schuch and F. Verstraete,
Matrix product states andprojected entangled pair states: Concepts, symmetries, and theorems , URL http://arxiv.org/abs/2011.12127 , Preprint (2020).[5] S. R. White,
Density matrix formulation for quantum renormalization groups , Phys.Rev. Lett. , 2863 (1992), doi:10.1103/PhysRevLett.69.2863.[6] S. R. White, Density-matrix algorithms for quantum renormalization groups , Phys. Rev.B , 10345 (1993), doi:/10.1103/PhysRevB.48.10345.[7] A. E. Feiguin and S. R. White, Finite-temperature density matrix renormal-ization using an enlarged Hilbert space , Phys. Rev. B , 220401 (2005),doi:10.1103/PhysRevB.72.220401.[8] S. R. White, Minimally entangled typical quantum states at finite temperature , Phys.Rev. Lett. , 190601 (2009), doi:10.1103/PhysRevLett.102.190601.[9] A. J. Daley, C. Kollath, U. Schollw¨ock and G. Vidal,
Time-dependent density-matrixrenormalization-group using adaptive effective Hilbert spaces , J. Stat. Mech. TheoryExp. , P04005 (2004), doi:10.1088/1742-5468/2004/04/P04005.[10] M. Zwolak and G. Vidal,
Mixed-state dynamics in one-dimensional quantum latticesystems: a time-dependent superoperator renormalization algorithm , Phys. Rev. Lett. , 207205 (2004), doi:10.1103/PhysRevLett.93.207205.[11] M. J. Hartmann, Polariton crystallization in driven arrays of lossy nonlinear resonators ,Phys. Rev. Lett. , 113601 (2010), doi:10.1103/PhysRevLett.104.113601.[12] C. Joshi, F. Nissen and J. Keeling,
Quantum correlations in the one-dimensional drivendissipative XY model , Phys. Rev. A , 063835 (2013), doi:10.1103/PhysRevA.88.063835.17 ciPost Physics Submission [13] L. Bonnes, D. Charrier and A. M. L¨auchli, Dynamical and steady-state properties ofa bose-hubbard chain with bond dissipation: A study based on matrix product operators ,Phys. Rev. A , 033612 (2014), doi:10.1103/PhysRevA.90.033612.[14] M. Schir´o, C. Joshi, M. Bordyuh, R. Fazio, J. Keeling and H. T¨ureci, Exotic attrac-tors of the nonequilibrium rabi-hubbard model , Phys. Rev. Lett. , 143603 (2016),doi:10.1103/PhysRevLett.116.143603.[15] J. Cui, J. I. Cirac and M. C. Ba˜nuls,
Variational matrix product operators for thesteady state of dissipative quantum systems , Phys. Rev. Lett. , 220601 (2015),doi:10.1103/PhysRevLett.114.220601.[16] E. Mascarenhas, H. Flayac and V. Savona,
Matrix-product-operator approach to thenonequilibrium steady state of driven-dissipative quantum arrays , Phys. Rev.A , 022116(2015), doi:10.1103/PhysRevA.92.022116.[17] A. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert and S. Montangero, Positive tensor network approach for simulating open quantum many-body systems , Phys.Rev. Lett. , 237201 (2016), doi:10.1103/PhysRevLett.116.237201.[18] J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio and D. Rossini,
Clustermean-field approach to the steady-state phase diagram of dissipative spin systems , Phys.Rev. X , 031011 (2016), doi:10.1103/PhysRevX.6.031011.[19] D. Kilda and J. Keeling, Fluorescence spectrum and thermalization in a driven coupledcavity array , Phys. Rev. Lett. , 043602 (2019), doi:10.1103/PhysRevLett.122.043602.[20] E. Gillman, F. Carollo and I. Lesanovsky,
Numerical simulation of critical dissipativenon-equilibrium quantum systems with an absorbing state , New Journal of Physics (9),093064 (2019), doi:10.1088/1367-2630/ab43b0.[21] H. Landa, M. Schir´o and G. Misguich, Multistability of driven-dissipative quantum spins ,Phys. Rev. Lett. , 043601 (2020), doi:10.1103/PhysRevLett.124.043601.[22] E. M. Stoudenmire and S. R. White,
Studying two-dimensional systems with the den-sity matrix renormalization group , Annu. Rev. Condens. Matter Phys. , 111 (2012),doi:10.1146/annurev-conmatphys-020911-125018.[23] F. Verstraete, M. M. Wolf, D. Perez-Garcia and J. I. Cirac, Criticality, the area law, andthe computational power of projected entangled pair states , Phys. Rev. Lett. , 220601(2006).[24] J. Jordan, R. Or´us, G. Vidal, F. Verstraete and J. I. Cirac, Classical simulation ofinfinite-size quantum lattice systems in two spatial dimensions , Phys. Rev. Lett. ,250602 (2008), doi:10.1103/PhysRevLett.101.250602.[25] R. Or´us and G. Vidal,
Simulation of two-dimensional quantum systems on an infinitelattice revisited: Corner transfer matrix for tensor contraction , Phys. Rev. B , 094403(2009), doi:10.1103/PhysRevB.80.094403.[26] J. Jordan, Studies of Infinite Two-Dimensional Quantum Lattice Systems with ProjectedEntangled Pair States , Ph.D. thesis, The University of Queensland, Australia (2011).18 ciPost Physics Submission [27] A. Kshetrimayum, H. Weimer and R. Or´us,
A simple tensor network algorithm fortwo-dimensional steady states , Nat. Commun. , 1291 (2017), doi:10.1038/s41467-017-01511-6.[28] S. Finazzi, A. Le Boit´e, F. Storme, A. Baksic and C. Ciuti, Corner-space renormalizationmethod for driven-dissipative two-dimensional correlated systems , Phys. Rev. Lett. ,080604 (2015), doi:10.1103/PhysRevLett.115.080604.[29] R. Rota, F. Minganti, C. Ciuti and V. Savona,
Quantum critical regime in aquadratically driven nonlinear photonic lattice , Phys. Rev. Lett. , 110405 (2019),doi:10.1103/PhysRevLett.122.110405.[30] G. Carleo and M. Troyer,
Solving the quantum many-body problem with artificial neuralnetworks , Science , 602 (2017), doi:10.1126/science.aag2302.[31] N. Yoshioka and R. Hamazaki,
Constructing neural stationary states for open quantummany-body systems , Phys. Rev. B , 214306 (2019), doi:10.1103/PhysRevB.99.214306.[32] A. Nagy and V. Savona, Variational quantum monte carlo method with a neural-network ansatz for open quantum systems , Phys. Rev. Lett. , 250501 (2019),doi:10.1103/PhysRevLett.122.250501.[33] M. J. Hartmann and G. Carleo,
Neural-network approach to dissipa-tive quantum many-body dynamics , Phys. Rev. Lett. , 250502 (2019),doi:10.1103/PhysRevLett.122.250502.[34] F. Vicentini, A. Biella, N. Regnault and C. Ciuti,
Variational neural-network ansatzfor steady states in open quantum systems , Phys. Rev. Lett. , 250503 (2019),doi:10.1103/PhysRevLett.122.250503.[35] F. Minganti, A. Biella, N. Bartolo and C. Ciuti,
Spectral theory of Liou-villians for dissipative phase transitions , Phys. Rev. A , 042118 (2018),doi:10.1103/PhysRevA.98.042118.[36] D. Kilda, A. Biella, F. Nissen and J. Keeling, iPEPO , https://github.com/The-iPEPO-Project/iPEPO (2019).[37] B. Bruognolo, Tensor network techniques for strongly correlated systems , Ph.D. thesis,Ludwig-Maximilians-Universit¨at M¨unchen (2017).[38] M. Lubasch, J. I. Cirac and M.-C. Banuls,
Algorithms for finite projected entangled pairstates , Phys. Rev. B , 064425 (2014), doi:10.1103/PhysRevB.90.064425.[39] C. Hubig and J. I. Cirac, Time-dependent study of disordered modelswith infinite projected entangled pair states , SciPost Physics (3) (2019),doi:10.21468/SciPostPhys.6.3.031.[40] A. Kshetrimayum, M. Goihl and J. Eisert, Time evolution of many-body localized systemsin two spatial dimensions , URL http://arxiv.org/abs/1910.11359 , Preprint (2019).[41] C. Hubig, A. Bohrdt, M. Knap, F. Grusdt and J. I. Cirac,
Evaluation of time-dependentcorrelators after a local quench in iPEPS: hole motion in the t-J model , SciPost Phys. , 21 (2020), doi:10.21468/SciPostPhys.8.2.021.19 ciPost Physics Submission [42] A. Kshetrimayum, M. Goihl, D. M. Kennes and J. Eisert, Quantum time crystals withprogrammable disorder in higher dimensions , URL http://arxiv.org/abs/2004.07267 ,Preprint (2020).[43] C. Mc Keever and M. H. Szyma´nska,
Dynamics of two-dimensional open quantum latticemodels with tensor networks , Preprint (2020), .[44] H. Weimer, A. Kshetrimayum and R. Or´us,
Simulation methods for open quantum many-body systems , URL http://arxiv.org/abs/1907.07079 , In press (2020).[45] A. A. Gangat, I. Te and Y.-J. Kao,
Steady states of infinite-size dissipative quan-tum chains via imaginary time evolution , Phys. Rev. Lett. , 010501 (2017),doi:10.1103/PhysRevLett.119.010501.[46] L. Vanderstraeten, J. Haegeman, P. Corboz and F. Verstraete,
Gradient methods forvariational optimization of projected entangled-pair states , Phys. Rev. B , 155123(2016), doi:10.1103/PhysRevB.94.155123.[47] P. Corboz, Variational optimization with infinite projected entangled-pair states , Phys.Rev. B , 035133 (2016), doi:10.1103/PhysRevB.94.035133.[48] G. Vidal, Efficient classical simulation of slightly entangled quantum computations , Phys.Rev. Lett. , 147902 (2003), doi:10.1103/PhysRevLett.91.147902.[49] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems , Phys.Rev. Lett. , 040502 (2004), doi:10.1103/PhysRevLett.93.040502.[50] G. Vidal, Classical simulation of infinite-size quantum lattice systems in one spatialdimension , Phys. Rev. Lett. , 070201 (2007), doi:10.1103/PhysRevLett.98.070201.[51] R. Orus and G. Vidal, Infinite time-evolving block decimation algorithm beyond unitaryevolution , Phys. Rev. B , 155117 (2008), doi:10.1103/PhysRevB.78.155117.[52] P. Corboz, R. Or´us, B. Bauer and G. Vidal, Simulation of strongly correlated fermionsin two spatial dimensions with fermionic projected entangled-pair states , Phys. Rev. B , 165104 (2010), doi:10.1103/PhysRevB.81.165104.[53] T. Nishino and K. Okunishi, Corner transfer matrix renormalization group method , J.Phys. Soc. Jpn. , 891 (1996), doi:10.1143/JPSJ.65.891.[54] T. Nishino and K. Okunishi, Corner transfer matrix algorithm for classical renormaliza-tion group , J. Phys. Soc. Jpn. , 3040 (1997), doi:10.1143/JPSJ.66.3040.[55] P. Corboz, J. Jordan and G. Vidal, Simulation of fermionic lattice models in two dimen-sions with projected entangled-pair states: Next-nearest neighbor hamiltonians , Phys.Rev. B , 245119 (2010), doi:10.1103/PhysRevB.82.245119.[56] H. W. Bl¨ote and Y. Deng, Cluster Monte Carlo simulation of the transverse Ising model ,Phys. Rev. E , 066110 (2002), doi:10.1103/PhysRevE.66.066110.[57] A. Kshetrimayum, M. Rizzi, J. Eisert and R. Or´us, Tensor network annealing al-gorithm for two-dimensional thermal states , Phys. Rev. Lett. , 070502 (2019),doi:10.1103/PhysRevLett.122.070502. 20 ciPost Physics Submission [58] J. Jordan, R. Or´us and G. Vidal,
Numerical study of the hard-core Bose-Hubbard model on an infinite square lattice , Phys. Rev. B79