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.

In the modelling of these systems, the inclusion of degrees of freedom which are external to the lattice, such as a driving field or a bath of oscillators, requires extending the description from a closed to an open quantum lattice model, as illustrated in Fig. 1.Open quantum systems are often well described by a Lindblad master equation [21] which facilitates the study of a range of collective phenomena including non-equilibrium criticality [22][23][24][25][26][27], quantum chaos [28,29] and time-crystallinity [30][31][32], many of which have no counterparts in closed systems at equilibrium.However, to better understand, control and utilise the dissipative non-equilibrium dynamics of correlated quantum systems, simulation techniques which are scalable to large lattices are still missing, especially in higher dimensions.
The investigation of large many-body quantum systems is hindered by the exponential growth of the Hilbert space.As the size of the system increases, solving the Lindblad master equation exactly using methods such as diagonalization of the Liouvillian or averaging over ensembles of exact quantum trajectories [33][34][35] quickly become infeasible.To simplify the problem, many have resorted to a mean field type approximation [25,[36][37][38][39][40][41][42] in which correlations between small individual subsystems are approximated by an average field.This simplification, however, may often give qualitatively incorrect results in regions where inter-subsystem correlations be-come important -for instance, near criticality.Moreover, key aspects such as entanglement and quantum information cannot be treated at this level.Progressing beyond mean field approximations should therefore involve the systematic inclusion of correlations between subsystems in a controlled and tractable manner.
In this vein, phase space methods such as those based on the Wigner [43], Positive-P [44] and Q [45] representations attempt to find classical stochastic processes for which the hierarchy of couple moments is a good approximation to that of the quantum problem.For highly nonlinear problems, phase space techniques often fail dramatically in important regimes [46][47][48].Cluster based methods [24,49] separate large lattices into small clusters and capture correlations within lattice sites belonging to each cluster, an approach which can become inaccurate when correlation lengths exceed cluster sizes.Variational approaches [50,51] based on the parametrization of the state in terms of a suitable functional and their optimisa-FIG.1.An open quantum lattice model of interacting spins.Nearest neighbour spins are coupled via a hopping J and interact with an external bath via a coherent drive Ω and/or a dissipative process γ.The open system can be modelled by describing the unit cell and its environment using a tensor network.
tion relies on good intuition, which may not be available for some problems.Recently, methods based on neural networks and the variational minimization of an appropriate cost function [52][53][54] have provided an interesting proof of concept, however, like most of these methods, they are restricted to small system sizes or may fail to capture long range correlations.
A different approach is to restrict the growth of the system's Hilbert space by retaining only the most important correlations or most probable states [55].Tensor network (TN) methods [56] belong to this class.Here, truncation of Hilbert space is controlled by the so called bond dimension (usually denoted D or χ) of indices which connect a set of tensors representing the quantum state.In the context of closed quantum many-body systems, the significant success of TN methods is underpinned by an area law in the growth of entanglement entropy possessed by ground states of gapped Hamiltonians [57].For open systems the picture is much less clear.In particular it is not obvious whether transient or steady states can be efficiently represented by a TN.Nevertheless, in the context of dissipative or driven-dissipative systems, we can reasonably expect that in many cases, dissipative processes should curtail the growth of entanglement and limit correlations generated by entangling dynamics.
Despite this expectation, TN algorithms for open systems [58][59][60][61][62] have mostly been restricted to onedimensional lattices where the simple geometry plays a central role in the algorithm.In dimensions greater than one, progress has been limited.The work of [63] introduced the Infinite Projected Entangled Pair Operator (iPEPO) to represent the mixed state of an infinite periodic two-dimensional square lattice and employed the so called simple update (SU) algorithm to apply Lindblad dynamical map evolving the system in real time towards a steady state.Although SU is efficient, in order to integrate the equation of motion it isolates a subsystem -for example one unit cell -from the rest of the lattice and applies the dynamical map to the subsystem in isolation until a steady state is reached.It has been questioned whether this approach can produce accurate results and there are concerns over the convergence of this method in non-mean-field regimes [64].While algorithms going beyond SU exist for closed and finite temperature systems [65][66][67], advancing beyond the SU approach in the driven-dissipative context remains undeveloped.
In this paper we devise a new TN method to accurately simulate time dynamics and steady states of many-body quantum lattice models in two spatial dimensions and directly in the thermodynamic limit.The method uses the iPEPO as an ansatz for the mixed state of the open system and incorporates techniques inspired by those presented in [68] -Full Environment Truncation (FET) and fixing the network to Weighted Trace Gauge (WTG)to calculate accurate time dynamics and steady state solutions of open quantum lattice models.The central step in the algorithm involves finding an optimal truncation of enlarged bonds with respect to an objective function appropriate for mixed quantum states.
The method successfully reproduces numerically exact calculations for both dynamics and steady-states while also agreeing with results obtained using the so called Corner Space Renormalization method of [55].Importantly, it performs well in non-mean field limits, proving able to capture substantial correlations in the presence of dissipation and therefore enabling the study of regimes which are accessible to current experiments but lie well beyond the applicability of existing techniques.
The paper is organised as follows.In section II we describe the algorithm including a brief introduction to the Lindblad master equation and the TN ansatz.As a benchmark we calculate time dynamics of a dissipative transverse quantum Ising model in section III A and find that the systematic inclusion of correlations -controlled by the TN bond dimension -coupled with the incorporation of the unit cell's environment when truncating enlarged bonds yields results which agree very well with the exact dynamics.Furthermore we demonstrate the applicability of the algorithm outside the exactly solvable regime.In section III B we show that the FET method outperforms the SU method by finding more optimal truncations of enlarged bonds by removing redundant internal correlations in the network.Finally in section III C we show that lattice models with drive and dissipation can also be treated using this method and compare steady state results for a driven-dissipative hard core boson model with literature values.In section IV we conclude with a short discussion.

A. Master Equation
The goal of the algorithm is to calculate time dynamics and steady states of translationally invariant twodimensional quantum lattice models, which interact with a bath via a Lindblad master equation (1) ( = 0) where Ĥ governs the coherent dynamics of the system and the dissipator D, which models the coupling of the system to its bath has the form with Lα being the Lindblad operators.We focus on the case of time-independent nearest neighbour Hamiltonians such that H can be decomposed as a sum of Hermitian operators which act non-trivially on at most two nearest neighbour lattice sites.Although the algorithm allows for up to two-local dissipators, for simplicity, we focus only on local coupling to the environment such that each Lindblad operator acts on one site only and respects the translational invariance of the Hamiltonian.

B. TN Ansatz
We represent the system's density matrix ρ (t) as an infinite Projected Entangled Pair Operator (iPEPO).The iPEPO is composed of a network of tensors {A j }, where we associate each node j of the network with one site of the square lattice shown in Fig. 2 (a).To reflect the translational invariance of the system and to simplify the algorithm, we use a pair of independent tensors A j and A l to represent the unit cell.The infinite system is the repetition of this unit cell over the two-dimensional plane.Each sixth-rank tensor A has a pair of physical indices of dimensions d and a set of four bond indices of dimension D, reflecting the coordination number z = 4 associated with a square lattice.The physical dimension d corresponds to the dimension of the local Hilbert space at each lattice site (d = 2 for the two-level spin), whereas D is a variational parameter which controls the accuracy of the ansatz.It is convenient to use the vectorized form of the density operator, which at the level of the iPEPO corresponds to vectorization of the pair of local Hilbert space indices as shown in Fig. 2 (a) and has the effect of transforming the iPEPO into the form of a infinite Projected Entangled Pair State (iPEPS) commonly used in TN algorithms for two-dimensional closed systems [56].Finally, to each unique bond we associate a bond matrix σ.
As with other algorithms based on Matrix Product Operators (MPOs), the PEPO ansatz is not inherently positive and therefore not all PEPOs represent physical states.For the present case of an infinite PEPO (iPEPO) we do not have access to the full spectrum of eigenvalues and it has been shown for the case of MPOs that the problem of deciding whether a given iMPO represents a physical state in the thermodynamic limit is provably undecidable [69].We therefore rely on the positivity of the dynamical map to maintain the physicality of the iPEPO throughout the time evolution and find in practice that, in most cases, the reduced density matrices calculated from the iPEPO are physical.
We refer to all of the spins in the system which are not part of the unit cell as its environment (see Fig. 1.), not to be confused with system's bath which is accounted for in the Lindblad master equation (1).Since the system is infinite, we represent the environment approximately by associating to each tensor in the unit cell an effective environment E. E is itself made up of a set of tensors including four corner transfer matrices C µν and four half row or half column tensors T µ , where the labels µ and ν take the appropriate first letter of left, right, up and down as illustrated in Fig. 8. (a-b).
We consider two distinct types of effective environment.The "trace effective environment" E tr of Fig. 8 (a) is calculated by first tracing over the local Hilbert space dimensions d of the tensors at each node of the network giving the set of fourth-rank tensors {a tr j } as shown in Fig. 7 (a).We use E tr to calculate the reduced density matrices of the system.Secondly, the "Hilbert-Schmidt Environment Truncation (FET) algorithm is used to find the isometries ũ and ṽ which maximize the fidelity between truncated and untruncated bonds.
effective environment" E hs (Fig. 8 (a)) is that formed by first contracting a hs j = A j A † j giving the Hilbert-Schmidt inner product of the tensor A i with itself, where all bond indices {D} are left open as shown in Fig. 7 (c).E hs is used during the algorithm to calculate an optimal truncation of enlarged bond dimensions as discussed in section II D. In both cases we calculate the effective environment using a corner transfer matrix method [70][71][72][73][74][75].In particular, we use a variant of the Corner Transfer Matrix Renormalization Group (CTMRG) algorithm [76] which makes use of an intermediate SVD to improve stability, details of which are given in Appendix A.

C. Time Evolution
To calculate dynamics and find a TN representation of the steady state we use a time evolving block decimation (TEBD) algorithm.The time evolution is obtained by application of the dynamical map ρ t = e tL ρ 0 .In principle it may also be possible to find the steady state directly by searching for the ground state of the Hermitian operator → L † L, for example, via imaginary time evolution.However, in general, L † L is a highly non-local operator and is therefore not straightforward to implement using standard techniques for an infinite systems [77].Finally, access to the transient dynamics is often of direct interest in many physical contexts.
The dynamical map e tL is approximated by a set of Trotter layers as it is common in algorithms based on TEBD.In particular, consider the evolution of the state from a time t to a short time later t + τ , then, in vectorized notation, where we note that the density matrix is vectorized column-by-column, the dynamical map takes the form The Liouvillian superoperator L is two-local and can therefore be written as a sum of superoperators acting on nearest neighbours of the square lattice, where the labels α and β correspond to the coordinates of the lattice site j and l respectively.The full Liouvillian takes the form The Hamiltonian part of the evolution is included in the superoperator H and the dissipative part in the superoperator D each are constructed as shown in equations ( 5) and ( 6) respectively: We then split the vectorized operators in the exponent into those acting on even and odd pairs of lattice sites along both the x and y lattice dimensions, giving four sets of vectorized operators L e x , L o x , L e y and L o y where which allows us to decompose e tL into a set of layers via a Trotter decomposition with τ = t/n where n 1 is the Trotter number with Each dynamical map in the decomposition is applied to pairs of nearest neighbour tensors A j and A l in turn.
We first construct the linear map L (A j A l ) where the linear operator L j ,l j,l acts on the pair of tensors A j and A l such that A j A l behaves as a vector in the linear map as illustrated in Fig. 2. (b).By repeated application of this map, an approximation to the tensor e τ L A j A l (Fig. 2 (c)) is calculated using Krylov subspace methods, eliminating the need for explicit calculation of e τ L , where τ is a real number for the case of real time evolution.
To complete the update, the resulting tensor A j,l = e τ L A j A l needs to be decomposed into a new pair of tensors A j and A l , illustrated in Fig. 2 (c-d).Typically this is done via singular value decomposition (SVD), where in general, the new bond dimension D -equal to the number of singular values associated with the SVD -will be enlarged (D > D) and therefore needs to be truncated in an appropriate way for the algorithm to remain efficient, in particular, we would like to truncate D back to D after each dynamical map.

D. Truncation of Enlarged Bonds
For TNs without closed loops (acyclic), finding an optimal truncation benefits greatly from the ability to efficiently apply a gauge transformation and re-cast a network to a so called canonical form, for details we refer the reader to [56].For TNs with closed loops (cyclic) however, such a canonical form cannot be defined uniquely and truncating the enlarged bond in an optimal way is much less straightforward.Moreover, cyclic TNs can host so called internal correlations which have no influence on properties of the quantum state but can cause computational problems if they are allowed to accumulate [68].
After applying the dynamical map we choose to decompose the tensors using SVD and truncate the bond irrespective of the state of the environment, leaving a new dimension D ≥ D chosen such that only those singular values greater than some small tolerance D 1 are retained.We are then left with a bond matrix σ with the remaining D singular values along its diagonal and the tensors A i and A j as shown in Fig. 2 (d).The final step in the truncation involves replacing σ with the product ũσṽ † where ũ and ṽ are isometries of dimension (D , D) such that ũũ † = ṽṽ † = I and σ is a new D dimensional diagonal bond matrix.The enlarged bond is then truncated by contracting A j and A l with ũ and ṽ as illustrated in Fig. 2

(e).
To calculate the set ũ, σ and ṽ we adapt the Full Environment Truncation (FET) algorithm of [68], which prescribes a method to find the truncation of an internal index of an arbitrary network for closed systems, optimal with respect to a fidelity measure for pure states.In our case, since we are dealing with an open system, we optimize the truncation with respect to an objective function suitable for mixed states.More precisely, we maximize a mixed state fidelity measure between the state ρ in which the enlarged bond dimension is left untruncated and the state φ in which the same bond has been truncated by ũ, σ and ṽ.Supposing that a global maximum is found, this procedure finds the isometries which leave φ as close as possible to ρ with respect to the chosen fidelity measure.
We choose to maximize the fidelity F (ρ, φ), which has the Hilbert-Schmidt inner product of ρ and φ in its numerator and the geometric mean of their purities tr(ρ 2 ) and tr(φ 2 ) in its denominator [78] Since squaring F is convex, the ρ and φ which maximize F 2 (ρ, φ) also maximize F (ρ, φ).We therefore construct F 2 (ρ, φ) tr ρ 2 as a Rayleigh quotient of tensors which can be maximized to find an optimal ũ, σ and ṽ.Details of the optimization procedure are given in Appendix B.
Finally, there exists a gauge freedom across the newly truncated bond which we fix to so called Weighted Trace Gauge (WTG) as described in [68].This allows for the recycling of the environment E hs calculated for use at each FET step of the algorithm as an initial guess for the renormalization procedure (CTMRG in our case) which precedes the following FET step thereby reducing the number of renormalization iterations required at each step.We refer to the algorithm outlined in this section as Full Environment Truncation in Weighted Trace Gauge (WTG+FET ).
It is straightforward to recover a Simple Update (SU) method by bypassing the FET and WTG steps above and instead choosing both ũ → ũsu and ṽ → ṽsu as D × D matrices with all diagonal entries equal to one and all other entries equal to zero and by retaining the D largest singular values of σ in the truncated σsu .In general, the set of ũ, ṽ and σ we find using FET are not equivalent to ũsu , ṽsu and σsu showing that, in the general case, SU does not yield a truncation which is optimal with respect to the objective function we use.A comparison between SU and WTG+FET is made in section III B.

A. Dissipative Transverse Ising Model
As a first benchmark of the algorithm we simulate dynamics of a dissipative transverse quantum Ising model with Hamiltonian where V is the hopping coupling, h x is the strength of a transverse field and z is the lattice coordination number which we set to z = 4 for the square lattice.The spins undergo dissipation at a rate γ described by local Lindblad jump operators Lj = √ γ 1 2 σy j − iσ z j , which are the same at each lattice site.For zero transverse field h x /γ = 0, the purely dissipative dynamics D (ρ dis ) = 0 drive the system towards a steady state ρ dis = |↓ x ↓ x | which does not commute with the Hamiltonian and thus ordered phases of the Hamiltonian can be frustrated by the dissipation.Moreover, in the specific case of h x /γ = 0, this Liouvillian belongs to a family of efficiently solvable dissipative models [79] (see Appendix D for further details) in which correlations remain localized and therefore the Liouvillian admits an efficient exact solution for local observables.We denote this method EXACT and use it as a benchmark.
For all parameters considered, we initialize the lattice spins in a product state ρ 0 = |↑ z ↑ z | and simulate their evolution in time in strongly dissipative (V /γ = 0.2, h x /γ = 0), moderately dissipative (V /γ = 1.2, h x /γ = 1.0) and weakly dissipative (V /γ = 4.0, h x /γ = 0) regimes, as well as in a regime (V /γ = 0.5, h x /γ = 1.0) which does not admit an efficient solution using the EXACT method.For all results pertaining to this model we choose D = 10 −8 and set the convergence criteria for both the CTMRG and FET algorithms to 10 −10 .We choose a time step τ γ = 0.01 in all cases except for the weakly dissipative regime where we choose τ γ = 0.005.
In each regime we calculate reduced density matrices ρ j and ρ l for each lattice site labelled j and l in the twosite unit cell as well as the set of four nearest neighbour reduced density matrices ρ jl and four next nearest neighbour reduced density matrices ρ jj where j and j are at a distance of 2 lattice constants rather than √ 2, i.e. they are in the same row or column.Although we find that all reduced density matrices within each set are equivalent to a high precision, it is convenient to plot expectation values averaged over each set.We therefore calculate the average magnetization m x = 1 2 (tr (σ x ρj ) + tr (σ x ρl )) as well as the average purity of the single site reduced density matrices Π 1 = 1 2 tr ρ2 j + tr ρ2 l as function of time.To compare larger reduced density matrices we calculate S xx 12 and S xx 13 , where S xx jl (t) = tr(σ x j ⊗ σx l ρ t ), again averaged over the four possible choices for j and l.Finally we show the infidelity I(t) = 1 − F(t) of each truncation averaged over the four trotter layers which make up every time step τ where F is the mixed state fidelity equation (9).Results are plotted for a range of bond dimensions D and the environment dimensions χ tr and χ hs , where we choose χ tr = χ hs = χ in each case, and where χ tr and χ hs are associated to the effective environments E tr and E hs , respectively.Finally, we have confirmed the convergence of the results with respect to increasing χ tr and χ hs in all results shown.

Strong Dissipation
In Fig. 3 (a) we plot the results of the strongly dissipative regime, in which the dissipative process dominates and where the spins are strongly damped.The exact dynamics of the system can be summarised as follows.From the initial product state, the average single site expecta- Figure 3 (a.v)plots the infidelity of truncation I(t), the qualitative behaviour of which is similar for all values of D. As the the dynamics progress from the initial product state and correlations begin to deviate from zero, I(t) increases from I 1 where the error introduced by truncation of enlarged bonds is negligible, to a larger finite value which indicates that the truncation causes the state to deviate slightly from the exact dynamics, nevertheless, for D = 4 and D = 5, I(t) remains below ≈ 10 −10 at all times and is an indicator of the accuracy of the results.We note here that, I(t) has a dependence on the time step τ and this should be considered when comparing this parameter across different values of τ .

Moderate and Weak Dissipation
An example of the moderate dissipation regime is presented in Fig. 3 (b).In this case, the hopping strength is comparable to the dissipation and therefore the exact dynamics display some transient oscillations which are quickly damped by the dissipation.Here again, the exact solution contrasts significantly from the D = 1 solution in which the dynamics tend towards a pure steady state with all spins in the |↓ state.We find that WTG+FET reproduces the exact dynamics to good precision for the single site observables for D > 3.While S xx 12 and S xx 13 also show good agreement with EXACT.A weak dissipation case for V /γ = 4.0 and h x /γ = 0.0 is plotted in Fig. 3 (c).The weakly damped oscillations of the EXACT results at early times reflect the dominance of the hopping term in this regime.While the D = 1 solution gives incorrect results, the results for D = 5 and D = 6 reproduce the exact solution early in the transient phase and begin to deviate from the exact dynamics after approximately tγ = 2 − 3 while still retaining the same qualitative behaviour.The fact that a larger bond dimension is required to reproduce the exact results is indicative of the greater role played by correlations in this coherent hopping dominated regime.

Outside Exactly Solvable Regime
For finite transverse field h x , the Lindblad master equation does not fulfil the conditions for an efficient exact solution using the EXACT method and correlations may not remain localised, nevertheless WTG+FET makes no assumption as to extent of correlations and should therefore be applicable for these parameters.As an example, a case for V /γ = 0.5 and h x /γ = 1.0 is presented in Fig. 3 (d).Using WTG+FET we find that the dynamics converge as the iPEPO bond dimension is increased.Results for D ∈ (1, 4, 5, 6) converge very well for D ≥ 5.The behaviour of the system is similar to the efficiently solvable cases; after some transient phase, the initial pure product state tends towards a correlated mixed state which is qualitatively different from the mean field solution.The infidelity of truncation Fig. 3 (d.v)remains below I(t) < 10 −8 for the converged results, which is in line with previous benchmarking results.

B. Comparison with Simple Update
To highlight differences between the WTG+FET and SU truncation methods, we compare the results calculated using each method in the moderate damping regime (V /γ = 1.2, h x /γ = 0) of section III A for a range of bond dimensions.All parameters are the same for both methods; τ = 0.01 and D = 10 −8 and CTRMG and FET convergence criteria set to 10 −10 , with the only difference being in how ũ, ṽ and σ are calculated.
As well as comparing the observables m x (t) and S xx 12 (t), we provide a quantitative measure of the accuracy of each method by calculating the trace distance between the EXACT reduced density matrix at each time step and the corresponding reduced density matrix calculated using the different TN methods.In particular we find the trace distance T 2 (t) of the nearest neighbour reduced density matrices T 2 (ρ jl , φ jl ) = 1  2 tr (ρ jl − φ jl ) † (ρ jl − φ jl ) where T 2 (t) is averaged over the four nearest neighbour reduced density matrices of the two site unit cell.By observing m x (t), and S xx 12 and the trace distance T 2 in Fig. 4 (a-c) it is clear that the SU method does not reproduce the EXACT results to the same accuracy as WTG+FET .Fig. 4 (a) and its inset demonstrates that, while WTG+FET shows clear systematic improvement in accuracy as D is increased, SU shows only minor and not clearly systematic reduction in T 2 (tγ = 10) even if D is increased well beyond that for which WTG+FET demonstrates good convergence.For values of D > 3, T 2 is consistently about an order of magnitude smaller for WTG+FET than for SU, demonstrating the much better compression and greater accuracy of WTG+FET.The observables in Fig. 4 (b-c) calculated using SU deviate from the EXACT dynamics considerably compared to WTG+FET (compare to Fig. 3 (b)), at times tγ 2, the SU method struggles to accurately capture the EX-ACT dynamics for all bond dimensions shown.
Finally we compare how the two algorithms deal with internal correlations in the network and compare the fidelity of truncation at each time step.TNs with closed loops (or cyclic TNs) can suffer from an accumulation of internal correlations, which do not contribute to any property of the quantum state.To achieve an optimal TN representation of the state at each truncation step, it is necessary to remove these internal correlations.Furthermore, a build up of these correlations can lead to problems in computation and breakdown of algorithms [68].The cycle entropy S cycle defined in [68] prescribes a way of quantifying the extent of internal correlations in the network, and is conveniently expressed in terms of the bond environment, details of its calculation in the present case are given in Appendix E. The cycle-entropy S cycle plotted in Fig. 5 (a) shows the extent of internal correlations in the network as a function of time.Initially the network, which represents a product state, has no internal correlations.In time, the extent of internal correlations grows and saturates at a finite value.Importantly, S cycle grows more slowly and saturates at a smaller value for WTG+FET than it does for SU, illustrating that the proper truncation of bonds reduces the extent of internal correlations in the network.Although the growth of S cycle in this case is relatively benign, the failure of SU to curtail the accumulation of internal correlations may contribute to the breakdown of the algorithm in some circumstances.As a final comparison we plot the infidelity of truncation I as a function of time for the two different methods in Fig. 5 (b) and find that the WTG+FET method outperforms SU , decreasing the infidelity between truncated an untruncated bonds by approximately an order of magnitude.Although the variational degree of the ansatz is the same in each case -they have same D and χ -the method by which enlarged bonds are truncated is crucially important in finding an optimal representation, thereby greatly reducing the accumulation of errors do to inadequate truncation and ultimately giving the most accurate results.

C. A Driven-Dissipative Hard Code Boson Model
In driven-dissipative quantum lattice models dissipation to the bath is replenished via a coherent or incoherent drive.Driven-dissipative systems constitute an important class of models with direct relevance to experimental platforms such as driven coupled photon arrays in a variety of architectures [80].In this section we calculate steady state properties of a driven-dissipative hard core boson model which can be mapped to a lattice of interacting spin-1/2 particles.The Hamiltonian is given in the rotating frame by where ∆ = ω p − ω c is the detuning between the pump frequency ω p and the on site energy ω c , F is the pump field strength, J is the hopping coupling and the sum j,l runs over nearest neighbours in the lattice of coordination number z.The spins undergo dissipation at a rate γ described by a Lindblad operator Lj = √ γ σ− j , which is the same at each site and where the spin raising and lowering operators are defined as σ± ≡ 1 2 (σ x ± iσ y ).We compare steady state expectation values with those calculated using the Corner Space Renormalization method [55].To this end we consider an array of TABLE I. Steady state values of a hard core boson model on an infinite square lattice with parameters ∆/γ = 5.0, F/γ = 2.0 and J/γ = 1.0 calculated using WTG+FET .In each case we use a time step of τ γ = 0.0025.For comparison we tabulate results for the same parameters from the corner space renormalization method [55] for different sizes Nx × Ny.
hard core bosons with ∆/γ = 5, F/γ = 2 and J/γ = 1 and calculate the average single site boson density n = 1/2(n j +n l ), the nearest neighbour ( j, l )correlation functions g (2) averaged over all combinations of ( j, l ), where g Finally we calculate the average real part of [tr (σ − ρ ss )] at each lattice site.
Staring from an initial product state, we find the steady state for a set of parameters D, χ and D , where convergence in time is achieved when all expectation values ô up to next nearest neighbour fulfil a convergence criterion of t < 10 −6 where We use the steady state iPEPO calculated for one set of variational parameters as an initial state for the next until convergence to the desired precision is achieved.Results of this procedure are given in TABLE I along with comparable results from [55].
The steady state values converge as the iPEPO variational parameters are increased and are comparable to the results of the Corner Space Renormalization method.Where we might expect increasing the Corner Space Renormalization lattice size N x × N y will give results closer to the WTG+FET method, which represents the thermodynamic limit directly, we find that the opposite is true, with a lattice size of 4 × 4 closer in agreement to WTG+FET than 8 × 8.This discrepancy could be due to finite size effects or spatial symmetry breaking, which may be present in the N x × N y results, and is not observed in the iPEPO solution where we have enforced two-site translational invariance by choosing a two-site unit cell.

D. Anisotropic Dissipative XY Model
Having demonstrated the capabilities of the algorithm, by choosing an interesting example we now show that the method is highly suitable to address physical questions.In particular, two dimensional systems can host a unique set of phenomena, here we explore the stability with respect to fluctuations of a spontaneously symmetry broken staggered-XY (sXY) phase in the steady state of an ansiotropic dissipative XY model.While the mean field theory predicts that the sXY phase is stable in two dimensions, it is not clear whether it remains accessible if fluctuations at the microscopic level are accounted for and if any long range order associated with the sXY phase is present.The anisotropic dissipative XY model has a Hamiltonian of the form with a nearest neighbour hopping J and coordination number z = 4, as well as dissipation described by local Lindblad operators Lj = √ Γσ − j at each lattice site.The (Gutzwiller) mean-field (MF) phase diagram, plotted in Fig. 6. (e), was studied in [81] and shows that, for J/Γ > 1/4, the steady state hosts a staggered-XY symmetry broken phase in which the spins divide into A and B sublattices with angles ±θ relative to the x = y line on the Bloch sphere as depicted in Fig. 6(c).The spontaneous breaking of this continuous U (1) symmetry means that θ can take any value and allows for vortexlike topological defects in the lattice.The question of whether or not the sXY phase is accessible in two dimensions if corrections beyond MF theory are accounted for has previously been addressed using a Keldysh field theory approach [81,82].There, an effective model is constructed by mapping the spins to bosons, an approach which does not capture the microscopic physics of the spin model but addresses the behaviour in the long wavelength limit.In that approximation they found that the steady state physics of the effective model is described by a partition function in the same universality class as the classical XY model and therefore one should expect a Kosterlitz Thouless transition in two dimensions.However it is also predicted in [82], based on a simple MF theory analysis, that the effective temperature of the model will be greater than the Kosterlitz Thouless temperature, such that the ordered phase will not be accessible when quantum fluctuations are included and any long range algebraic order will be absent or at least significantly diminished.We can now use our method to address this question exactly by directly solving the microscopic spin model close to the transition point J/Γ = 1/4, where the MF theory is expected to break down.Moreover, we are able to give a quantitative picture of the system by calculating not only local observables as a function of time, but also spatial correlation functions in the steady state.
We first find the steady state iPEPO representation of the model for a bond dimension D = 1 -equivalent to a MF solution -at J/Γ = 0.3, which lies just within the sXY phase.To do this, we initialize the iPEPO in a state for which the symmetry is explicitly broken σ x A = − σ x B = 1 and calculate the D = 1 steady state with WTG+FET.Then, using the symmetry broken D = 1 iPEPO solution as an initial state, we systematically add quantum fluctuations by calculating steady states for bond dimensions D ∈ [3,4,5,6] until convergence.Results of this procedure are presented in Fig. 6.For bond dimensions D = 3 we find that the system remains in the sXY phase.For D ∈ [4,5,6], however, the spin magnetization m z (t), which is uniform across the lattice, is slightly modified and the magnetizations m x (t) and m y (t) on each sublattice slowly tend towards zero such that the continuous symmetry is no longer broken-depicted in Fig. 6 (d)-and the sXY phase is therefore unstable to fluctuations, corroborating the Keldysh field theory predictions of [82].This proves that long wavelength fluctuations captured by the approximate theory dominate over other microscopic fluctuations.In Fig. 6 (f) we plot the correlation function S xx k,j = σ x j σ x k (note that σ x j σ x k = 0) which shows a staggered structure reminiscent of the sXY phase where correlations at a radii r (see Fig. 6 (e) inset) corresponding to an odd number of steps on the lattice are zero, whereas even step correlations are finite and decay with r.Considering only the even step correlations in Fig. 6.(f) inset, we find that the decay is well approximated by an exponential function of the form S xx r∈even ∝ e −ηr with η ≈ 1.07, any long range algebraic order which may have been associated to the symmetry broken phase is not present in the iPEPO solution suggesting that the system is in the disordered phase.Good convergence is found for D = 6 and τ γ = 0.01 resulting in infidelity of truncation I(t) < 10 −9 .

IV. DISCUSSION
We have developed a new TN algorithm capable of accurately simulating dynamics of dissipative quantum lattice models on a two-dimensional square lattice directly in the thermodynamic limit.The method adapts the Full Environment Truncation (FET) and Weighted Trace Gauge (WTG) fixing techniques of [68] to dealing with the iPEPO TN ansatz for mixed states.Comparisons with exact numerical results demonstrate an excellent accuracy of the method and its performance across different dissipative regimes.Contrasting with the more efficient but much less accurate simple update truncation scheme, we have proven that it is necessary to optimally truncate enlarged bonds to obtain accurate results.We have shown the applicability of the technique for calculating steady state properties of driven-dissipative systems by comparison with literature results.The methods performs well in regimes where mean-field approximation fails, proving able to capture substantial correlations in the presence of dissipation.Finally we have shown that a staggered-XY phase of the dissipative anisotropic XY model predicted by mean field theory is not stable if correlations are included and while a remnant of the staggered structure remains in the correlation function, it's decay is well approximated by an exponential function and no long range order remains.
As with similar algorithms for iPEPS, the principal contribution to the computational complexity of the algorithm comes from the calculation of the effective environment which is updated at each time step (here using CTMRG).The leading cost of the version of CTMRG we use arises from a singular value decomposition of order O(χ 3 hs D 6 ), improvements in performance can therefore be achieved by optimizing this step, for instance, using a fixed point method such as the FPCM [76] or approximating the effective environment by using a boundary matrix product state to represent the boundary of the system.Numerous algorithms have been developed to calculate the fixed point including a time-evolving block decimation (TEBD) [83,84] or variational MPS-tangent space methods (VUMPS) [76,[85][86][87][88] and can lead to significant speed up for TNs which are close to being critical [76].
As well as accurately determining steady state properties such as long range equal-time correlation functions, this work facilitates the calculation of more complex dynamical properties e.g.dynamical correlation functions and fluorescence spectra of strongly correlated driven dissipative quantum lattice models.A significant advantage of both the FET method of truncating enlarged bonds and the WTG method of fixing the TN gauge, is that they can be used in tensors networks of arbitrary geometries, provided the bond environment can be calculated efficiently.In this regard, straightforward adaptations of the method we have presented in this work could be used to treat driven-dissipative models with longer range interactions or those defined on more complicated network structures such as hyperbolic lattices [5] as well as problems related to functional quantum biology [89,90].
Appendix A: Calculating the Effective Environments Given tensors representing the unit cell of the 2D lattice, we calculate the effective environments E tr and E hs using a variant of the Corner Transfer Matrix Renormalization Group (CTMRG) algorithm [70][71][72][73][74][75][76].To improve stability and convergence properties of the CTMRG algorithm as well as the conditioning of the bond environment Υ jl we find it helpful to use the variant of CTMRG presented in [76], which makes use of an intermediate singular value decomposition.Following [76] we refer to Fig. • (e) We define F lu a ≡ U u a Su a where singular values of magnitudes (relative to the largest singular value) less than some small tolerance are truncated to improve stability.
• (f) We next use the so called biorthogonalization procedure (see [76] for further details) to calculate P l and P l − , the first step of which is to contract F lu a with F ld a and perform a SVD to find W l a , Ql a and the diagonal matrix Σl 2 a .
• (g,j) We calculate the projectors P l a = F lu a Ql a Σl +  of the tensors at neighbouring sites j and l.(d) The bond environment Υ j,l is the contraction of E hs j,l and the updated tensors A j and A l with enlarged bonds {D j } ≥ D. (e-g) Using Υ j,l the terms in the fidelity between the truncated (φ) and untruncated (ρ) density matrices are calculated by contracting with the isometries u, v and the bond matrix σ.

FIG. 2 .
FIG. 2. Main steps in the time evolution algorithm.(a) The vectorized form of the iPEPO.(b) The contraction of A j and A l with the dynamical map L. (c) Singular value decomposition (SVD) of e τ L AjA l .(d) The D singular values of relative tolerance greater than D are retained in the diagonal bond matrix σ .(e) The isometries ũ and ṽ truncate the enlarged bond from D to D giving the new bond matrix D. (f) The updated tensors Ãj and Ãl .(g) The FullEnvironment Truncation (FET) algorithm is used to find the isometries ũ and ṽ which maximize the fidelity between truncated and untruncated bonds.

FIG. 3 .
FIG.3.Dynamics of the dissipative Ising model for hx/γ = 0 in (a) strong (V /γ = 0.2), (b) moderate (V /γ = 1.2) and (c) weak (V /γ = 4.0) spin damping regimes calculated using WTG+FET for a range of bond dimensions D and superimposed with results calculated using the EXACT method.(d) Shows results for a regime not applicable to the EXACT method (V /γ = 0.5 and hx/γ = 1.0) but which can be treated with WTG+FET .In each case we plot (i) the magnetization m x (t), (ii) the average purity Π1 of the single site reduced density matrices, (iii) the nearest-neighbour S xx 12 , (iv) the next-nearest neighbour S xx 13

FIG. 4 .
FIG. 4. Comparison between WTG+FET (solid), SU (dotted) and EXACT (dashed line) in the moderately damped regime of the dissipative Ising model V /γ = 1.2, h x /γ = 0.0.(a) Trace distance T2(t) as a function of time and at tγ = 10 (inset) for a range of bond dimensions.(b) Magnetization m x (t) and (c) nearest-neighbour S xx 12 show that WTG+FET outperforms SU.

FIG. 5 .
FIG. 5. (a) The accumulation of internal correlations in time quantified by the cycle entropy S cycle (t) is more effectively curtailed by WTG+FET and (b) the infidelity of truncation I(t) is an order of magnitude smaller than SU at each truncation step.Results form moderately damped regime of the dissipative Ising model V /γ = 1.2, h x /γ = 0.0 with D = 4.

FIG. 6 .
FIG. 6. Fate of staggered-XY phase at J/Γ = 0.3.(ab) 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 ∈ [3, 4, 5, 6].(c-d) Representation of a 2 × 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 transition at J/Γ = 1 4 .Inset: Radii r = | k − j| of an odd (red) and even (blue) number of steps on the lattice.(f) The correlation function S xx j,k = σ x j σ x k versus distance r has a staggered form which 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 decaying with r.Inset: Exponential fit to even step correlations giving exponent η ≈ 1.07.
FIG. 7. Tensor diagrams representing some of the steps involved in performing the left-move component of the CTMRG algorithm used to calculate the effective environment E hs .

1 / 2 ab , F ru a ≡ Su 1 / 2 ab
V u † a , F ld a ≡ U d a Sd 1/2 ab and F rd a ≡ Sd 1/2 ab V d †

a
and P l − a = F ld a W l † a Σl + a with Σl + being the Moore-Penrose pseudoinverse of Σl.• (f-h) We repeat steps (d-j) to calculate P l b and P l − b by replacing a ↔ b in the upper and lower half system transfer matrices.Using these projectors the updated environment tensors T l b , T l a , Clu a , Clu b , Cld a , Cld b are calculated and normalized as shown in Fig.7 (h,j,k).This is one iteration of the left-move component of this CTMRG algorithm.A similar sequence of steps is used to perform the rightmove, up-move and down-move steps in CTMRG.The set of directional moves are repeated in series until the vectors of singular values of the corner transfer matrices converge.It is possible to perform right-move at the same time as left-move by following the biorthogonalization routine starting with F ru b and F rd b calculated in step (b) above, similarly for up-move and down-move.

FIG. 8 .
FIG. 8.The environment of the unit cell.(a) The trace effective environment E tr j of the iPEPO tensor Aj used to calculate the d × d reduced density matrix ρj.(b) The Hilbert-Schmidt effective environment E hs j of the tensor a hs j used in constructing the bond environment.(c) The effective environment E hs j,l