# Steady States of Infinite-Size Dissipative Quantum Chains via Imaginary Time Evolution

###### Abstract

Directly in the thermodynamic limit, we show how to combine imaginary and real time evolution of tensor networks to efficiently and accurately find the nonequilibrium steady states (NESS) of one-dimensional dissipative quantum lattices governed by the Lindblad master equation. The imaginary time evolution first bypasses any highly correlated portions of the real-time evolution trajectory by directly converging to the weakly correlated subspace of the NESS, after which real time evolution completes the convergence to the NESS with high accuracy. We demonstrate the power of the method with the dissipative transverse field quantum Ising chain. We show that a crossover of an order parameter shown to be smooth in previous finite-size studies remains smooth in the thermodynamic limit.

Introduction–The out-of-equilibrium behavior of dissipative many-body systems is of relevance to experimental platforms such as trapped ionsSchindler et al. (2013), cold atoms Schwager et al. (2013), superconducting circuitsHouck et al. (2012); Schmidt and Koch (2013); Hur et al. (2016), and nanoelectromechanical systems Lozada-Vera et al. (2016). Theoretical activity has recently increased around developing numerical methods for determining the non-equilibrium steady states (NESS) of such dissipative lattices Werner et al. (2016); Weimer (2015a); Cui et al. (2015); Finazzi et al. (2015); Mascarenhas et al. (2015); Dorda et al. (2015); Raghu et al. (2016); Jin et al. (2016); Kshetrimayum et al. (2016). This includes tensor network methods Schollwöck (2011); Verstraete et al. (2008); Cirac and Verstraete (2009); Orús (2014), which serve as an efficient numerical ansatz for states that obey an area lawEisert et al. (2010) (i.e. have small correlations between their bipartitions). It is proved that the mutual information of the NESS of a local dissipative quantum system satisfies an area law Brandao et al. (2015); Mahajan et al. (2016), assuring that their tensor network representation, if found, will be computationally efficient. Also, a proof for the stabilityCubitt et al. (2015) of NESS against local system perturbations assures that theoretically determined NESS of translationally invariant systems are of relevance to experiments, where translational invariance can only be approximate.

The focus of the present work is on the problem of numerically finding the tensor network representation of the NESS of translationally invariant one-dimensional systems in the thermodynamic limit. The infinite-size version of the time-evolving block decimation (iTEBD) algorithm Vidal (2007) enables a (local) time-evolution directly in the thermodynamic limit using a matrix product state (MPS) that spans only a single unit cell. Also, the iTEBD algorithm can be used for both imaginary and real time evolution with local operators.

When imaginary time evolution is performed with a Hamiltonian, the fixed point lies in the ground state manifold of the given Hamiltonian. An efficient tensor network can accommodate the entire imaginary time evolution trajectory if the fixed point obeys an area law. Furthermore, the convergence toward the fixed point is exponentially fast with a rate proportional to the energy gap of the Hamiltonian. Therefore, imaginary time evolution with iTEBD is a very efficient way to obtain ground states of local Hamiltonians in the thermodynamic limit.

For real time evolution of dissipative systems, the fixed point lies in the NESS manifold. However, real time evolution with iTEBD Orús and Vidal (2008) is not always an efficient way to obtain the NESS since portions of the real time evolution trajectory from the initial state to the NESS may not obey an area law even when the NESS does so itself Cui et al. (2015). Further, the rate of convergence to the NESS may be very slow Mascarenhas et al. (2015); Cai and Barthel (2013). For finite-size chains, recent worksCui et al. (2015); Mascarenhas et al. (2015) tackle these problems using variational methods. These approaches, though accurate, can not be easily generalized to infinite-size systems.

All prior approaches to finding groundstates and NESS with tensor networks have exploited the area law property of the target state to achieve efficiency. By additionally exploiting a closely related but distinct physical property, exponential decay of two-point correlators, we here show how to alleviate the inefficiencies of real time evolution to the NESS for infinite-size one-dimensional tensor networks: we construct an auxiliary local Hamiltonian such that its ground state approximates the NESS, and perform iTEBD imaginary time evolution with the auxiliary Hamiltonian. In this way, we are able to bypass any highly correlated portions of the real time evolution trajectory and to arrive in the area law-obeying neighborhood of the NESS exponentially quickly. We further improve the convergence to the NESS using a real time evolution of the Lindbladian.

The auxiliary Hamiltonian approach that we present holds potential for the thermodynamic limit of higher dimensional dissipative systems as well, since both imaginary time evolution and variational optimization Corboz (2016) have been demonstrated to obtain ground states with the higher dimensional tensor network ansatz of infinite-size projected entangled pair states (iPEPS) Jordan et al. (2008).

As a demonstration of our method, we use it to probe the thermodynamic limit of a crossover in an order parameter of the one-dimensional dissipative transverse field Ising model. To our knowledge, non-mean-field studies of this crossover have only previously been performed in the finite-size limit. Our result shows that this crossover remains smooth in the thermodynamic limit.

Method–The equation of motion for a (discrete) quantum system coupled to a Markovian environment is given by the Lindblad master equation (LME) Breuer and Petruccione (2002). Under the Choi isomorphism (), the LME is

(1) |

with

(2) |

where is the system Hamiltonian and are dissipative operators.

A matrix product density operator (MPDO)Verstraete et al. (2004) is a tensor network that serves as an efficient ansatz for a mixed state density matrix of a one-dimensional lattice when the correlations between real space bipartitions of the state are small Zwolak and Vidal (2004); Verstraete et al. (2004). Under the Choi isomorphism, the MPDO Verstraete et al. (2004) can be written as a matrix product state (MPS) Perez-Garcia et al. (2007): , where the are tensors of dimension and is referred to as the “bond dimension”. The chief advantage of such an ansatz is that while needs to be exponentially large in the system size for the MPS to be exact, a much smaller (i.e. computationally tractable) value of yields extremely high accuracy for states that obey an area law Eisert et al. (2010). Another advantage of this ansatz is that the entanglement spectrum (denoted below) between subblocks of any bipartition of the lattice is readily calculated Vidal (2004). In the case of infinite-size systems with translational invariance, the MPS or MPDO is referred to as an “iMPS” or “iMPDO”, and it can span a Hilbert space as small as a single unit cell of the target physical state Vidal (2007).

The NESS (denoted ) of dissipative systems described by the LME are defined by , with the constraint that the are vectorized forms of positive operators with unit trace (i.e. physically valid density matrices; see below for further discussion). The authors of Ref. Cui et al. (2015) observe that will also be the ground state of the nonlocal Hamiltonian . For finite-size chains, they present a variational method for finding the (non-degenerate) ground state of as a way of determining . However, their method does not apply directly (i.e. without costly extrapolation from finite-size scaling) in the thermodynamic limit. Since is nonlocal, imaginary time evolution with iTEBD also can not be used to find its ground state. Here we show that it is possible to construct a local auxiliary Hamiltonian such that the iTEBD algorithm may be used to approach the infinite-size NESS via imaginary time evolution:

(3) |

where is any vectorized density matrix such that .

As an initial motivation we observe that if is a local Hamiltonian with positive eigenvalues, will be a nonlocal Hamiltonian with the same ground and excited states as the local . This suggests the possibility of finding a local Hamiltonian whose ground state is at least a good approximation to the ground state of the nonlocal .

We assume that is a translationally invariant local operator; this corresponds to the case of a translationally invariant local system Hamiltonian and translationally invariant local dissipation operators . can therefore be expressed as a sum of translationally invariant local terms () and we may write

(4) |

We may use as a benchmark, and as a figure of merit. It has been analytically proven in Ref. Brandao et al. (2015) that the two-point correlator shows an exponential decay in NESS. Therefore the long-range couplings in do not play a significant role in determining its ground state and we may truncate the second sum in Eq. (4) by setting . Further taking the root of each remaining term we arrive at the proposed local auxiliary Hamiltonian for imaginary time evolution to the NESS:

(5) |

If the (numerical) gap between the lowest two eigenvalues of is less than one, will increase the gap since is positive semi-definite. This will yield faster convergence toward the ground state of (for imaginary time evolution the rate of convergence is proportional to the gap). While for will not generally commute with for , we find that the advantage gained (see next section) outweighs this drawback.

The form of can change as long as . with larger support leads to longer range couplings in ; for of infinite support, becomes equal to (when ). We may therefore decrease the distance between and the ground state of by increasing the support of , albeit with an increased computational cost.

Rather than relying on imaginary time evolution alone, the following hybrid method can be more efficient: with small , relatively large timestep, and of small support, imaginary time evolution can be used to rapidly converge to the area law-obeying neighborhood of the NESS, after which real time evolution with iTEBD can minimize with successively smaller timesteps and successively larger . The Trotter error of iTEBD vanishes exponentially in the timestep size, and the error associated with a finite vanishes exponentially in when there is an area law.

As mentioned, for the converged solution to be physically valid it must correspond to a positive operator. If the NESS is unique, it is shown in Ref. Cui et al. (2015) that this requirement is automatically satisfied when . We may therefore assume a unique NESS and assure (near-)positivity by attaining very small and converging the spectrum. This is similar to what is done in Ref. Cui et al. (2015), which assures (near-)positivity by variationally reducing to below a chosen threshold and converging various observables, and Ref. Mascarenhas et al. (2015), which assures (near-)positivity by variationally minimizing . As an alternative, which we do not implement here, positivity may be enforced even in the case of non-unique NESS by projecting onto the physical subspace. We outline this procedure in the Appendix.

Numerical Results–There have been various numerical investigations of the 1D quantum Ising model with nearest-neighbor coupling, uniform transverse magnetic field, and on-site dissipation Cui et al. (2015); Werner et al. (2005a, b); Yin et al. (2014); Orth et al. (2008); Schwager et al. (2013); Cai and Barthel (2013); Hu et al. (2013); Ates et al. (2012); Lee et al. (2011); Joshi et al. (2013); Marcuzzi et al. (2014); Weimer (2015b); Rose et al. (2016). Here we look at the following incarnation: the system is governed by Eqs. (1) and (2) with the Hamiltonian given by , where ()

(6) |

is the lattice site index, and is the transverse field strength; the local dissipation terms are given by , where is a decay rate from spin up to spin down.

We choose a tensor network that is structured such that two physical sites are associated to each of the numerical sites of a two-site iMPDO. For both imaginary and real time evolution, a fourth order Suzuki-Trotter expansion of the evolution operators is used. For the imaginary time evolution, we choose the following 4-local form () for :

(7) |

where the local dissipation term is given as . The overlap between and in this form is two sites; forms resulting in three site overlap would yield more intersite couplings in and therefore a better match between the ground state of and , but the above form is numerically more efficient and sufficient for our demonstration. It is straightforward to verify that this form satisfies .

For the finite-size version of this model, a previous numerical study Hu et al. (2013) shows that a smooth crossover occurs in the value of the up-spin density as is increased from small to large values. To demonstrate the theoretical validity of the imaginary time method put forth in the previous section, we probe this crossover (for ) in the thermodynamic limit with imaginary time evolution and real time evolution independently. The purpose here is to demonstrate that the imaginary time method can come close to the NESS very efficiently and robustly (i.e. with very crude parameters and convergence). In the Appendix, we check the convergence of the imaginary time evolution vs. both and imaginary time step size at the middle of the crossover. From there we determine that fixing and the imaginary timestep size at (in units of the inverse Ising interaction strength squared) is sufficient for this purpose. For real time evolution we also use a timestep size of (in units of inverse Ising interaction strength). The simulations are run until the spectrum is converged according to the following criterion: , where is the largest singular value in the iMPDO. At each value of the same initial state is used for both imaginary time and real time evolution. The results are illustrated in Fig. 1. The spectra in the middle panel shows that the imaginary time evolution () and the real time evolution have converged to a weakly correlated subspace. The bottom panel shows that is much smaller than the smallest energy scale in for real time evolution and imaginary time evolution (when ). The results in these two panels indicate proximity to the NESS. Because is smallest for real time evolution, the data in the top panel for the real time evolution can be considered the most accurate, and the increasing overall match between the imaginary time data and real time data for with higher shows that the accuracy of the imaginary time evolution improves with larger . Though we do not investigate it here, it is possible that the significant enhancement in accuracy with larger is due to the removal of metastability Rose (); Macieszczak et al. (2016). Taken together, these results demonstrate the conceptual validity of the imaginary time evolution method as an alternative to real time evolution for approaching the state space neighborhood of the NESS. The level of accuracy achieved with the imaginary time method is remarkable given the crude convergence scheme and enormous truncation of in forming . The physical explanation for the success of the truncation of is the exponential decay of the two point correlators in the NESS. For comparison, we show simulation results for a 2-local form of in the Appendix. The 2-local form also permits successful imaginary time evolution toward the NESS, but with poorer accuracy.

We note that our results indicate that in this one-dimensional case the crossover in remains smooth in the thermodynamic limit; this is in contrast to the sharp transition predicted for the two-dimensional version of the model Weimer (2015a, b). This is the first study of this crossover in the thermodynamic limit beyond the mean-field theory.

We also demonstrate the functionality of the full hybrid method (imaginary time evolution followed by real time evolution). For the same system parameters as the other data, the top panel of Fig. 1 shows the converged data from the hybrid method when the state at the end of the imaginary time evolution for is converged to with real time evolution.

To demonstrate that the hybrid method remains efficient where real time evolution becomes inefficient, with and , we evolve a highly correlated initial state with the two methods separately and find that the hybrid method converges rapidly in several minutes of computation time with while real time evolution over the same number of timesteps becomes trapped in an area law-violating sector of the Hilbert space even with as high as . The time evolution of is shown in Fig. 2, while the evolution of the spectra and further details are given in the Appendix. It is clear that the initial imaginary time evolution bypasses the highly correlated region of the real-time trajectory, and brings the iMPDO close to the neighborhood of the NESS; thus, an efficient real time evolution becomes possible, and very high accuracy can be achieved.

Discussion– Imaginary time evolution of tensor networks with iTEBD is an efficient way of approximating the ground states of many-body quantum systems in the thermodynamic limit when an area law holds for the ground state; the area law permits high accuracy with a computationally tractable bond dimension. Here we have shown that imaginary time evolution of an iMPDO with iTEBD is an efficient way of approximating the NESS of dissipative quantum chains when the NESS has both an area law and exponential decay of two-point correlators; the area law permits the iMPDO to approximate the NESS with high accuracy with a computationally tractable bond dimension, and the exponential decay of two-point correlators permits the construction of a local auxiliary Hamiltonian with which to perform imaginary time evolution with iTEBD. We have demonstrated that imaginary time evolution can very efficiently reach the area law-obeying state space neighborhood of the NESS. Real time evolution can be employed after the imaginary time evolution to efficiently obtain the NESS with very high accuracy. Alternatively, but perhaps less efficiently, imaginary time evolution may be used alone to reach the NESS by converging not only in the timestep size and bond dimension, but also the support size of the local terms in the auxiliary Hamiltonian.

In the course of demonstrating our method with the dissipative transverse field quantum Ising chain, we have shown that a crossover of an order parameter that is smooth in the finite-size limit remains smooth in the thermodynamic limit, in contrast to the two-dimensional case.

With the higher dimensional tensor network of iPEPS, both variational optimization and imaginary time evolution have been shown to work Jordan et al. (2008); Corboz (2016), and real time evolution toward NESS was also recently demonstrated Kshetrimayum et al. (2016); the method presented here may thereby be extended to higher dimensions.

We note that the imaginary time evolution method presented here can also apply to finite-size chains by using TEBD instead of iTEBD, and it would be interesting to compare the performance of the hybrid method in this paper with the variational methodsCui et al. (2015); Mascarenhas et al. (2015) for finding the NESS of finite-size chains.

YJK and AAG acknowledge discussions with Ian McCulloch, Roman Orús and Augustine Kshetrimayum, and the hospitality of the Institut für Physik at Johannes Gutenberg-Universität. AAG acknowledges early pedagogical interactions at the University of Queensland with Guifré Vidal and Ho N. Phien regarding MPS and TEBD. We thank Dominic C. Rose for insightful remarks regarding metastability. Some of the simulations were coded using the Uni10 tensor library Kao et al. (2015). This work is supported by the MOST in Taiwan through Grants No. 104-2112-M-002 -022 -MY3, 105-2112-M-002 -023 -MY3.

## References

- Schindler et al. (2013) P. Schindler, M. Müller, D. Nigg, J. Barreiro, E. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller, and R. Blatt, Nat. Phys. 9, 361 (2013).
- Schwager et al. (2013) H. Schwager, J. I. Cirac, and G. Giedke, Phys. Rev. A 87, 022110 (2013).
- Houck et al. (2012) A. A. Houck, H. E. Türeci, and J. Koch, Nat. Phys. 8, 292 (2012).
- Schmidt and Koch (2013) S. Schmidt and J. Koch, Annalen der Physik 525, 395 (2013).
- Hur et al. (2016) K. L. Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux, and M. Schiró, Comptes Rendus Physique 17, 808 (2016).
- Lozada-Vera et al. (2016) J. Lozada-Vera, A. Carrillo, O. P. S. Neto, J. K. Moqadam, M. D. LaHaye, and M. C. Oliveira, EPJ Quantum Technology 3, 1 (2016).
- Werner et al. (2016) A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Phys. Rev. Lett. 116, 237201 (2016).
- Weimer (2015a) H. Weimer, Phys. Rev. Lett. 114, 040402 (2015a).
- Cui et al. (2015) J. Cui, J. I. Cirac, and M. C. Bañuls, Phys. Rev. Lett. 114, 220601 (2015).
- Finazzi et al. (2015) S. Finazzi, A. Le Boité, F. Storme, A. Baksic, and C. Ciuti, Phys. Rev. Lett. 115, 080604 (2015).
- Mascarenhas et al. (2015) E. Mascarenhas, H. Flayac, and V. Savona, Phys. Rev. A 92, 022116 (2015).
- Dorda et al. (2015) A. Dorda, M. Ganahl, H. G. Evertz, W. von der Linden, and E. Arrigoni, Phys. Rev. B 92, 125145 (2015).
- Raghu et al. (2016) M. Raghu, C. D. Freeman, S. Mumford, N. Tubman, and B. Swingle, (2016), arXiv:1608.05074 [cond-mat] .
- Jin et al. (2016) J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio, and D. Rossini, Phys. Rev. X 6, 031011 (2016).
- Kshetrimayum et al. (2016) A. Kshetrimayum, H. Weimer, and R. Orús, arXiv preprint arXiv:1612.00656 (2016).
- Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
- Verstraete et al. (2008) F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57, 143 (2008).
- Cirac and Verstraete (2009) J. I. Cirac and F. Verstraete, J. Phys. A: Math. Theor. 42, 504004 (2009).
- Orús (2014) R. Orús, Ann. Phys. 349, 117 (2014).
- Eisert et al. (2010) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
- Brandao et al. (2015) F. G. Brandao, T. S. Cubitt, A. Lucia, S. Michalakis, and D. Perez-Garcia, Journal of Mathematical Physics 56, 102202 (2015).
- Mahajan et al. (2016) R. Mahajan, C. D. Freeman, S. Mumford, N. Tubman, and B. Swingle, arXiv preprint arXiv:1608.05074 (2016).
- Cubitt et al. (2015) T. S. Cubitt, A. Lucia, S. Michalakis, and D. Perez-Garcia, Commun. Math. Phys. 337, 1275 (2015).
- Vidal (2007) G. Vidal, Phys. Rev. Lett. 98, 070201 (2007).
- Orús and Vidal (2008) R. Orús and G. Vidal, Phys. Rev. B 78, 155117 (2008).
- Cai and Barthel (2013) Z. Cai and T. Barthel, Phys. Rev. Lett. 111, 150403 (2013).
- Corboz (2016) P. Corboz, Phys. Rev. B 94, 035133 (2016).
- Jordan et al. (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
- Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, 2002).
- Verstraete et al. (2004) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
- Zwolak and Vidal (2004) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
- Perez-Garcia et al. (2007) D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum Inf. Comput. 7, 401 (2007).
- Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
- (34) See Supplemental Material at [URL will be inserted by publisher] for convergence test of bond dimension and time steps, and a comparsion of real tiem evolution and hybrid method.
- Werner et al. (2005a) P. Werner, K. Völker, M. Troyer, and S. Chakravarty, Phys. Rev. Lett. 94, 047201 (2005a).
- Werner et al. (2005b) P. Werner, M. Troyer, and S. Sachdev, Journal of the Phys. Soc. of Japan 74, 67 (2005b).
- Yin et al. (2014) S. Yin, P. Mai, and F. Zhong, Physical Review B 89, 094108 (2014).
- Orth et al. (2008) P. P. Orth, I. Stanic, and K. Le Hur, Phys. Rev. A 77, 051601 (2008).
- Hu et al. (2013) A. Hu, T. E. Lee, and C. W. Clark, Phys. Rev. A 88, 053627 (2013).
- Ates et al. (2012) C. Ates, B. Olmos, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. A 85, 043620 (2012).
- Lee et al. (2011) T. E. Lee, H. Häffner, and M. C. Cross, Phys. Rev. A 84, 031402 (2011).
- Joshi et al. (2013) C. Joshi, F. Nissen, and J. Keeling, Phys. Rev. A 88, 063835 (2013).
- Marcuzzi et al. (2014) M. Marcuzzi, E. Levi, S. Diehl, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. Lett. 113, 210401 (2014).
- Weimer (2015b) H. Weimer, Phys. Rev. A 91, 063401 (2015b).
- Rose et al. (2016) D. C. Rose, K. Macieszczak, I. Lesanovsky, and J. P. Garrahan, (2016), arXiv:1607.06780 [cond-mat.stat-mech] .
- (46) D. C. Rose, (private communication).
- Macieszczak et al. (2016) K. Macieszczak, M. Guţă, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. Lett. 116, 240404 (2016).
- Kao et al. (2015) Y.-J. Kao, Y.-D. Hsieh, and P. Chen, Journal of Physics: Conference Series 640, 012040 (2015).

## Appendix A Appendix

### a.1 Positivity enforcement via projection onto physical subspace

We outline here a more rigorous method to ensure positivity of the state vector . A basis of physical state vectors can be constructed by performing the Choi isomorphism on a set of random physical density matrices (which may or may not be orthonormalized): . (The constructed physical basis may be optimized by restricting to basis states that obey an area law.) Then will correspond to a positive operator. Normalization, if needed, can be performed as mentioned in Ref. [Cui et al. (2015)]. Such a projection may, however, increase , especially since the constructed physical basis will not necessarily span the entire physical subspace. If the increase in after the projection is too large, the minimization of may be resumed with , or different random physical bases may be tried until one is found that gives a satisfactory result for .

### a.2 Convergence test of bond dimension and timestep size for imaginary time evolution

Here we show results for convergence tests of the imaginary time evolution method with the 4-local form of (see main text). For , the middle of the crossover in occurs at about , and we perform the tests with these parameters fixed at these values.

The first test assesses convergence vs. the bond dimension with a fixed timestep size of . A random initial state is chosen as the initial condition for the simulation with the smallest value of . When the simulation is stopped, is increased and the simulation is again started with the initial state as the final state at the previous value of . At each value of , the simulation is stopped when . The results are shown in Table 1.

D | ||
---|---|---|

5 | 0.251167 | 0.081437 |

7 | 0.251116 | 0.081345 |

9 | 0.251163 | 0.081198 |

11 | 0.251245 | 0.081079 |

13 | 0.251319 | 0.081008 |

15 | 0.251363 | 0.080971 |

17 | 0.251392 | 0.080945 |

19 | 0.251410 | 0.080928 |

21 | 0.251419 | 0.080920 |

23 | 0.251423 | 0.080916 |

25 | 0.251425 | 0.080913 |

The second test assesses convergence vs. the imaginary timestep size with a fixed bond dimension of . The test is conducted in the same manner as the first test. The results are shown in Table 2.

0.1 | 0.251516 | 0.080544 |

0.01 | 0.251380 | 0.080964 |

0.001 | 0.251367 | 0.081004 |

0.0001 | 0.251365 | 0.081008 |

### a.3 Simulation results for 2-local form of

For the dissipative Ising model specified in the main text, we here show a demonstration of the imaginary time method with a 2-local form of :

(8) |

It is straightforward to verify that this form satisfies , where is defined in the main text. Fig. 3 compares converged values of , the entanglement spectrum, and from real time evolution to those from imaginary time evolution. We use the same bond dimension, timestep sizes, and convergence criterion as for the 4-local (see main text). For each value of , the same random initial state is used for the real and imaginary time evolutions (for all ), which are both performed on a two-site iMPDO. Inside the crossover region, the data shows rough qualitative agreement between the different evolutions. Comparing with the results in the main text, it is clear that accuracy is greatly enhanced by using the 4-local form of instead of the 2-local form. This is because the 4-local form results in more inter-site couplings in , making the groundstate of a better approximation for the groundstate of .

### a.4 Comparison of hybrid method and real time evolution

In this section we give further details related to Fig. (2) in the main text, where the hybrid method is shown to be efficient where real time evolution becomes inefficient. The 4-local form of is used (see main text) and the simulation parameters are , , , and . The same hardware (a single CPU laptop) is used for all simulations.

For the pure real time evolution simulation, a highly entangled initial state is created with but before the first timestep the bond dimension is increased to with the additional 18 singular values () initialized at zero. The simulation is then run for timesteps so that the total physical runtime is the same as the largest timescale () in the system. Fig. (4) shows the time evolution of the spectrum. It is clear that the system state becomes trapped in an area-law violating portion of the statespace and does not reach the NESS. The computational time for the simulation is on the order of hours.

For the hybrid method simulation, the initial state is the same as the one used for the pure real time simulation but the bond dimension is kept at during both imaginary and real time evolution. The imaginary time evolution is run for timesteps, which takes about one minute, after which real time evolution is performed until , which takes about another eight minutes. The hybrid method thereby achieves in less than ten minutes of computation time, while pure real time evolution remains stuck at for more than ten hours of computation time. Fig. (5) shows the evolution of the spectrum during the full hybrid method evolution.