{'fraction_non_alphanumeric': 0.031622632686462475, 'fraction_numerical': 0.00917699321954641, 'mean_word_length': 4.211188953193218, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast convergence of imaginary time evolution tensor network algorithms by recycling the environment

We propose an environment recycling scheme to speed up a class of tensor network algorithms that produce an approximation to the ground state of a local Hamiltonian by simulating an evolution in imaginary time. Specifically, we consider the time-evolving block decimation (TEBD) algorithm applied to infinite systems in 1D and 2D, where the ground state is encoded, respectively, in a matrix product state (MPS) and in a projected entangled-pair state (PEPS). An important ingredient of the TEBD algorithm (and a main computational bottleneck, especially with PEPS in 2D) is the computation of the so-called environment, which is used to determine how to optimally truncate the bond indices of the tensor network so that their dimension is kept constant. In current algorithms, the environment is computed at each step of the imaginary time evolution, to account for the changes that the time evolution introduces in the many-body state represented by the tensor network. Our key insight is that close to convergence, most of the changes in the environment are due to a change in the choice of gauge in the bond indices of the tensor network, and not in the many-body state. Indeed, a consistent choice of gauge in the bond indices confirms that the environment is essentially the same over many time steps and can thus be re-used, leading to very substantial computational savings. We demonstrate the resulting approach in 1D and 2D by computing the ground state of the quantum Ising model in a transverse magnetic field.


I. INTRODUCTION
The study of interacting quantum many-body systems remains a central problem in modern physics. Understanding how the microscopic degrees of freedom organize themselves collectively is key to explaining the large variety of complex, emergent phenomena exhibited by quantum many-body systems. The study of such systems, for instance on a lattice, faces a major computational challenge, given the exponential growth of the Hilbert space dimension with the lattice size. There are several numerical approaches available, including exact diagonalization of the Hamiltonian on small lattices or quantum Monte Carlo, based on sampling an imaginary time evolution in Hamiltonians that are free from the sign problem. Among these and other methods, tensor network approaches stand up for their ability to approximating the full ground state wave-function of the system's ground state, which is efficiently stored in a tensor network state.
The matrix product states (MPS) 1,2 is the best known example of a tensor network state. It is the basis of the density matrix renormalization group (DMRG) method 3,4 to efficiently represent the ground state of a 1D quantum system in both finite and infinite [5][6][7][8][9][10][11] lattices. The MPS is also the basis of methods to simulate time evolution, including the time evolving block decimation (TEBD) algorithm 12,13 and its DMRG-like formulation 14,15 which have been recently applied to study dynamical properties of 1D quantum system [16][17][18] , as well as the more recent algorithms based directly on the time-dependent variational principle 19 . Here we will focus on the TEBD algorithm, originally proposed for finite systems and subsequently extended to address systems in the thermodynamic limit, where it is referred to as infinite TEBD algorithm (iTEBD) 20,21 . In the thermodynamic limit, the many-body state is represented by an infinite MPS (iMPS), which consists of a finite unit cell of MPS tensors that is repeated throughout the infinite chain. On the other hand, the projected entangled pair state (PEPS) [22][23][24] or tensor product state (TPS) [25][26][27] was introduced as a natural generalization of the MPS for 2D systems. In the thermodynamic limit, one can again use an infinite PEPS (iPEPS) 28,29 , which has a finite unit cell of tensors repeated throughout the infinite 2D lattice.
Both in 1D and 2D, there are two main strategies to find a good approximation to the ground state of a local Hamiltonian H: one can either minimize the expectation value of H, succinctly referred to as energy minimization, or one can simulate the action, on some initial state, of the evolution operator e −tH in imaginary time in the long time limit, referred to as imaginary time evolution. In a 1D lattice, energy minimization (e.g. through the DMRG algorithm) is preferred over imaginary time evolution, because it is known to converge faster to the ground state. In contrast, in 2D systems most current implementations are based on imaginary time evolution using arXiv:1411.0391v1 [quant-ph] 3 Nov 2014 the TEBD algorithm supplemented with a large variety of algorithms to contract the tensor network, including schemes based on a boundary MPS 28 , the corner transfer matrix [29][30][31][32][33][34] , and the tensor renormalization group and its various important generalizations [35][36][37][38] . Among the most exciting explorations of ground states in challenging 2D systems by means of imaginary time evolution on iPEPS one can find e.g. recent calculations of the t-J model on the square lattice 39,40 and honeycomb lattice 41 , and of the J1-J2 frustrated Heisenberg model on the square lattice 42 .
In order to simulate the action of the imaginary time evolution operator e −tH on some initial state, the TEBD algorithm proceeds by breaking this evolution operator into the product of a large number of operators e −δH , each representing a small time step, δ 1, and by then applying each of these operator e −δH successively. Each application of a small time step e −δH results in the growth of the bond dimension of the tensors in the tensor network, which in turn results in a significant increase of computational costs. In order to keep the computational cost in check, it is thus important to truncate the bond index back to its original size. An optimal truncation of the bond index of a tensor (that is, a truncation that minimizes the error introduced in the many-body state), requires the computation of the so-called environment, which contains information of how the truncation affects the rest of the many-body state.
In a 1D system with open boundary conditions (OBC), the cost of computing the environment is similar to the cost of applying the time evolution e −δH . However, in 2D the cost of computing the environment is much larger than the cost of applying e −δH , and in fact computing the environment is the major bottle-neck of the whole TEBD algorithm. So much so, that it has become a standard practice to bypass the computation of the environment altogether, and use instead just local information as a guide to perform the truncation of the bond indices, in what is known as the simplified update (SU) 43,44 . In spite of its wide use within the PEPS community, the SU is not well-justified, and it is known to sometimes lead to poor PEPS approximations of the ground state. In other words, the SU, with its lower computational cost, allows to use a PEPS with larger bond dimension; however it does not exploit this larger bond dimension in an optimal way, potentially producing worse approximations to the ground state than what one would have obtained with a properly optimized PEPS with smaller bond dimension. Therefore the so-called full update (FU), which requires the computation of the environment, is still the most reliable option in spite of its prohibitive cost.
In this paper we thus focus on simulating imaginary time evolution using TEBD with FU, and propose a strategy to reduce its computational cost. Our strategy aims mainly at speeding up PEPS computations of 2D ground states, but it is also be useful in MPS calculations in 1D systems with periodic boundary conditions (PBC) as opposed to the OBC case alluded to above, where computing the environment is also significantly more expensive than just applying e −δH . The key observation is that while at the beginning of the imaginary time evolution the many-body state is changing rapidly, in later iterations of the TEBD algorithm the changes are much less significant. This may not be directly apparent in the environment, which might still be subject to big changes. However, these big changes correspond mostly to a random choice in the gauge freedom existing in the tensor network. Specifically, on each bond index connecting two tensors of the MPS or PEPS, one always introduce the product of an invertible matrix and its inverse without affecting the represented many-body state. Moreover, it is important to remember that the environment is only used to determine the right choice of truncation of the bond dimension, and that a small change in the environment will not change this choice significantly. Therefore, if we manage to fix the gauge freedom in the environment so that an equivalent gauge is chosen over several time steps, then we may be able to compute the environment only once and then recycle or re-use it to guide the truncation during many time steps. As we will see, this strategy actually leads to very significant computational savings.
In order to be able to recycle the environment, we first need to bring the tensor network into a canonical form that fixes the gauge degrees of freedom in the bond indices. For an MPS with OBC, there is already a natural canonical form as well as a well-understood procedure to obtain it. For an MPS with PBC and a PEPS there is no obvious choice of a canonical form and it is not a priori clear that a canonical form even exists. Here we propose an iterative procedure that, when applied to an MPS with OBC, has the usual canonical form as its fixed point. The same iterative procedure can then be applied to an MPS with PBC and (after suitable generalization) to a PEPS, and we verify numerically that a unique fixed-point is reached in each case. We then define the canonical form of an MPS with PBC and of a PEPS to be this fixed-point. For the purpose of incorporating the environment recycling schemes into the TEBD algorithm, we introduce an iterative method to obtain the canonical form of an MPS and a PEPS 45 . We notice that this type of canonical form of an PEPS is similar to the "quasi-canonical" form used in Refs. 46 and 47 to study the ground state of a 2D system by means of the SU.
For concreteness, in this work we discuss the recycling of the environment in a ground state computation by imaginary time evolution in an infinite system both in 1D and 2D, using iMPS and iPEPS respectively. However, we emphasize that environment recycling can also be used in a more broader class of algorithms. On the one hand, one can recycle the environment when simulating imaginary time evolution on finite systems, even though such algorithms are not currently very popular. On the other hand, the very same ideas can also be applied to the computation of a classical partition function in 2D and 3D, using a scheme where a tensor network (e.g. an iMPS and an iPEPS, respectively) is used to encode the dominant eigenvector of a (1D or 2D, respectively) transfer matrix by a power method, where now this transfer matrix replaces the imaginary time evolution operator e −δH . Indeed, after applying this transfer matrix a few times, the tensor network starts to approach the dominant eigenvector; from that moment on, the corresponding environment will no longer change significantly over iterations and it can be recycled.
Let us next summarize the structure of the paper. In Sec. II, we briefly introduce MPS and PEPS. In Sec. III, we revise the canonical form of an iMPS and introduce an alternative iterative method to obtain it. We also provide methods to fix the remaining gauge degrees of freedom (mostly, complex phases) within a canonical iMPS. Then in Sec. IV we propose a canonical form for iPEPS and a scheme to fix all its gauge degrees of freedom. In Sec. V we present an environment recycling scheme for TEBD on iMPS and some benchmark results for 1D quantum Ising model with transverse magnetic field. Next, in Sec. VI the environment recycling scheme is developed for TEBD on iPEPS and applied to study 2D quantum Ising model with the transverse magnetic field. Finally, Sec. VII contains our conclusions.
Note added: in this paper we use the terms MPS and PEPS to denote tensor networks in 1D and 2D respectively, and the term TEBD to denote an algorithm that simulates a time evolution (in imaginary time) according to e −tH in order to obtain an MPS/PEPS approximation to the ground state of H. In the literature one often finds that the term "iPEPS algorithm"is used to denote what we here call instead the TEBD algorithm for iPEPS. In this paper we find it important to distinguish between the tensor network being used (e.g. MPS, PEPS, or even iMPS, iPEPS) and the algorithm that is employed to optimize it (DMRG, TEBD, etc).

II. TENSOR NETWORK STATES
In this section we briefly revise the matrix product state and the projected entangled pair state employed to represent the many-body state of 1D and 2D quantum lattice systems respectively.

A. Matrix product states
Let us consider a 1D quantum lattice system consisting of N sites. Each lattice site has an internal degree of freedom, for example a spin, and is represented by a local d-dimensional Hilbert space H d = C d . The N -site lattice is then described in the Hilbert space H = (H d ) ⊗N . A general pure state of the system can be written in the local basis as where |s 1 , s 2 , . . . , s N = |s 1 ⊗ |s 2 ⊗ . . . ⊗ |s N and |s i (for i = 1, . . . , N and s i = 1, . . . , d) is the orthonormal basis of the local d-dimensional Hilbert space at site i, and c s1s2...s N is a complex component of a rank-N tensor c. The total number of components c s1s2...s N is d N , and hence proliferates exponentially with the number of lattice sites. We can decompose the rank-N tensor c into a product of lower-rank tensors without changing the state of the system. More precisely, if the lattice has OBC, the pure state in Eq. (1) can be rewritten as follows, where A s1 and A s N corresponding to each value of s 1 and s N are the 1 × χ row and χ × 1 column vectors respectively, and each A s k (for k = 2, . . . , N − 1) corresponding to each value of s k is a χ × χ matrix. The pure state represented by Eq. (2) is called as a matrix product state; see Fig. 1(a). This MPS is parametrised by (N −2)χ 2 d+2χd variational parameters, and accordingly the state is represented by a number of parameters which increases polynomially with N . The coeficient χ ≥ 1 is called the bond dimension and understood as a refinement parameter. More precisely, χ plays an important role in qualifying the MPS representation as the larger the χ, the better the MPS representation. If χ = 1, the MPS represents a product state otherwise the MPS describes an entangled state. We can always choose an upper bound of χ to represent exactly the general state of N -site lattice, say χ = d N/2 , but this will result in poor use of computational resources because the MPS is no longer an efficient representation.
In case that the lattice obeys PBC, the MPS expressed by Eq. (2) can still represent correctly the state of the system. This is straightforward because the tensor c decomposed into a product of tensors {A s k } N k=1 does not impose any specific boundary condition on the lattice. However, it is not a feasible representation in the sense that there is no connection between the tensors representing the first and the last sites of the lattice. Therefore the correlation between these two sites is not transmitted through the link of the system. Alternatively, we can generalize this representation where all the tensors including the first and the last ones have two bond indices. The generalized MPS is then written as follows, Note that the trace (Tr) is taken to ensure that after contracting all the matrices {A s k } N k=1 we have a scalar quantity. A diagrammatic representation of the MPS with PBC is shown in Fig. 1(b), where the trace is taken by connecting the left bond of tensor A s1 with the right bond of tensor A s N . In reality, the representation of an MPS with PBC often causes some computational difficulties in simulation. More explicitly, simulation using MPS with PBC is much more computationally expensive than employing MPS with OBC. As the cyclic structure of the MPS with PBC, it requires larger computational resources to capture the amount of entanglement in the system. The computational cost scale as O(χ 5 ) 48 comparing to O(χ 3 ) when simulating using MPS with OBC. However, this scaling sometimes can be reduced to O(χ 3 ) 49,50 .
Note that MPS can also represent the state of an infinite-size homogeneous system using only a small number of parameters by employing its invariance under translations. This results in what is called infinite MPS (iMPS). An iMPS is normally constructed by repeatedly assigning a unit cell made up of a finite set of tensors along the lattice chain. For example, a one-site iMPS is represented as follows, In this iMPS, the unit cell contains a pair of tensors {Γ, λ} where the tensor Γ is positioned at each lattice site and the diagonal matrix λ lies on the link between each two adjacent sites. The iMPS is respresented by only χ 2 d + χ variational parameters. This is the most striking advantage of using translational invariant iMPS because we only need a small number of parameters to represent the ground state of a lattice system in the thermodynamic limit. Consequently, using an iMPS to represent a state when simulating a 1D infinite system obviously reduces the computational complexity as well as the computer memory.

B. Projected entangled pair states
Consider a 2D quantum system on a lattice of N sites where each lattice site is represented by a local Hilbert space H d ∼ = C d . To represent the system in terms of a TN, one option is to use the so-called projected entangled pair state [22][23][24] which is generally defined as is the coordinate of the lattice site i), which are connected to some set of neighbor sites according to the geometry of the lattice. Each tensor of the PEPS has n (the coordination number, or the number of nearest neighbors of one lattice site) bond indices of dimension D and one physical index of dimension d. The choice of n in the TN depends on the geometry of the lattice, and can in principle be chosen arbitrarily (although this affects the computational complexity of contracting the TN). The function F contracts all the tensors according to this pattern and then performs the trace to obtain a scalar quantity. In Fig. 2, we show diagrammatically several PEPSs for the systems corresponding to different geometries with OBC. In case the geometry of a system is a square lattice pattern with PBC (n = 4), the PEPS consists of tensors which have four bond indices and one physical index. Overall, the PEPS depends on O(N D 4 d) variational parameters. Note that the bond dimension D determines the quality of PEPS, and if D is chosen to be large enough, any state of a 2D lattice system can be well-represented by a PEPS. As a generalization of the MPS representation to 2D, PEPS is a good representation for a class of variational wave functions when studying the ground state of 2D quantum lattice systems. The rationale behind this is that a PEPS is constructed naturally to capture the geometry of the lattice pattern, such that the scaling of entanglement entropy of a subregion is proportional to the surface area of the subregion, it fulfills the area law [53][54][55] . More precisely, the entanglement entropy of an l ×l block scales as the number of bonds that are cut in the TN, which is O(l) 56 . Thus a PEPS can represent efficiently the ground state of a 2D quantum lattice system with a small number of parameters which scales at worst polynomially with the system size. Furthermore, it can also effectively capture the quantum correlations existing in the system. There is some prospect of capturing also dynamical properties, as long as the dynamical growth in entanglement entropy is not too large. Besides, as similar to MPS we can also employ an infinite PEPS (iPEPS) to represent a state of a 2D homogeneous system which is parametrized by a fairly small number of parameters.

III. CANONICAL FORM OF AN IMPS AND GAUGE DEGREES OF FREEDOM FIXING
In this section we revise the canonical form of an iMPS and introduce an iterative method to obtain it. This method as we will see can be generalized to apply to MPS with either OBC or PBC. Besides, some methods applied to fix the gauge degrees of freedom in a canonical iMPS will also be introduced explicitly.

A. Canonical form of an iMPS
For simplicity, let us consider an infinite lattice chain represented by a one-site translationally invariant iMPS in Eq. (4). This iMPS is said to be in its canonical form if λ contains the Schmidt coefficients and therefore we can always write it in the Schmidt decomposition at any bond of the chain as follows, with |φ α L and |φ α R being orthonormal bases. Equivalently, this iMPS must satisfy the orthonormality constraints defined as, where ρ L = ρ R = λ 2 and I are reduced density and identity matrices respectively. To transfer an arbitrary iMPS represented by Eq. (4) into a canonical form, one can apply the scheme introduced in Ref. [21]. Here, as an alternative, we propose an iterative method that can also handle well this task.
The key idea of this iterative method is quite similar to the one in Ref. [21]. The main difference is that we do not need the left-and right-dominant eigenvectors of a transfer matrix to find the resolutions but instead we will find the resolutions locally. More concretely, our scheme consists of the following fundamental steps: (i) Compute the matrices V L and V R which coefficients are defined as follows, which are shown graphically in Fig. 3(i). As V L and V R are Hermitian and positive matrices, we can (ii) Insert two resolutions I = (Y T ) −1 Y T and I = XX −1 into the bonds of the iMPS as illustrated in Fig. 3(ii). Then we get a newλ which contain the singular values from taking the SVD of . Re-assign Γ ≡Γ and λ ≡λ then go back to step (i).
The above steps are repeated until the spectrum consisting of singular valuesλ converge to a fixed point. Then the iMPS represented by {Γ,λ} fulfills the conditions in Eqs. (7,8), and thus is in canonical form. Note that because the canonical form of an iMPS is obtained locally without taking into account the information of the whole system, we can also apply this scheme to find the canonical form for MPS with PBC when it is translationally invariant under shifts of a certain unit cell of sites.

B. Fixing gauge degrees of freedom in an iMPS
We now introduce two possible methods to fix the gauge degrees of freedom between two different canonical iMPSs which represent a same physical state |Ψ . Suppose that the canonical iMPSs are one-site translationally invariant and denoted as |Ψ A and |Ψ B . In this notation, |Ψ A and |Ψ B are characterized by single tensors A and B respectively. For simplicity, let us assume that |Ψ A and |Ψ B are normalized. Due to the canonical form representations in these iMPSs, we come up with the relation A = U † BU where U can be a unitary matrix or a diagonal matrix consisting of complex phases. To fix the gauge degrees of freedom between |Ψ A and |Ψ B , we need to search for the matrix U . In the following, with assumption that U is unitary, we will explain in details how to obtain U with two different methods.

Direct method
Suppose that X is the left-dominant eigenvector of the transfer matrix created by two tensors A and B with corresponding left-dominant eigenvalue µ. We have the following equation Note that for a clear demonstration, we will incorporate graphical representations for some expressions in the text. For example, the Eq. (11) is shown graphically as = .
Let us replace A s αβ = χ γ,δ=1 U * αγ B s γδ U δβ into the above equation we have = . (12) Multiplying both sides of the the Eq. (12) with U † yields which is simplified as We realize that XU † is also an eigenvector of transfer matrix created by tensors B and B † with corresponding eigenvalue µ. However, as |Ψ B is defined in its canonical form, possible solution of Eq. (14) is XU † = I. Thus, we can easily extract U = X and µ = 1. The matrix U needed for fixing the gauge degrees of freedom between the two iMPSs |Ψ A and |Ψ B has been found exactly and this corresponds to the dominant eigenvector of the transfer matrix created by tensors A and B.

Iterative method
We can also determine the matrix U iteratively. To do that, let us define a local fidelity as follows, This equation is visualized as , where we have defined M = AU † B † . Our task is to find the matrix U that maximize the fidelity F . More Specifically, we first initialize an arbitrary U , e.g. U = I and then iterate the following steps: (i) Apply the SVD to decompose the matrix M such that M = W SV † .
(ii) Assign U = V W † and subtitute it back into the Eq. (15).
(iii) Compute the fidelity which is now defined as F = Tr(S).
(iv) Re-compute the matrix M = AU † B † and then go back to step (i).
The above steps are iterated until the fidelity converges to a fixed point, i.e. F = 1.

IV. CANONICAL FORM OF AN IPEPS AND GAUGE DEGREES OF FREEDOM FIXING
Similar to an iMPS, an iPEPS representing a pure state |Ψ of a 2D lattice system is not unique, and there exists a lot of different iPEPSs representing the same state physical |Ψ . These iPEPSs can be distinguished by some gauge degrees of freedom. For instance, for a one-site translationally invariant iPEPS |Ψ A , characterized by tensor A, one can always insert M −1 M = I and K −1 K = I (M and K are arbitrary invertible matrices) into horizontal and vertical links of the iPEPS |Ψ A respectively without changing its physical properties. More concretely, we define a new iPEPS |Ψ B characterized by tensor B given by which describes a same physical state as |Ψ A . The existence of gauge degrees of freedom sometimes causes lots of numerical difficulties in dealing with iPEPS representation. Therefore, we need to fix the gauge degrees of freedom in every iPEPS. This can be done by defining a new type of canonical form for an iPEPS in what follows.

A. Canonical form of an iPEPS
For the sake of simplicity, we consider a one-site translationally invariant iPEPS which is represented by tensors {Γ, λ h , λ v } and is visualised as in Fig. 4. We can also represent the iPEPS with a single tensor A by contracting the matrices λ h and λ v to tensor Γ such that A = λ h λ v Γ. Note that for a TN with closed loops such as an iPEPS, it is impossible to define a canonical form based on the orthonormality of a bipartite Schmidt decomposition.
However, in what follows we define a new type of canonical form for an iPEPS in the sense that the tensors in the unit cell of the iPEPS satisfy some particular constraints. This definition for a canonical form iPEPS is completely different from the one given in Ref. 52 where the canonical form of a finite PEPS is termed from characterizing the existence of symmetries given by the fact that two representations of the same injective PEPS are related by unique invertible matrices. In what follows, we define the canonical form based on the orthogonality constraints of local tensors of the iPEPS.
Definition: An iPEPS is defined in its canonical form if it is represented by the tensors {Γ, λ h , λ v } satisfying the following constraints, where we have defined The conditions defined by Eqs. [17][18][19][20] can be called the left, right, down and up canonical constraints for the canonical form of the iPEPS and are shown diagrammatically in Fig. 5. By generalizing the iterative method applied to transfer an iMPS into its canonical form, we can obtain the canonical form for an arbitrary iPEPS. In particular, we iterate the following steps which are also illustrated diagrammatically in Fig. 6,  λ 1h andλ 1v are then obtained by taking the SVDs (iii) Define a new tensorΓ which is obtained by contracting all the remaining tensors V h , Fig. 6(iii). Re-assign Γ ≡Γ, λ h ≡λ 1h and λ v ≡λ 1v then go back to step (i).
The above steps are repeated until the fixed point of the spectrum described by diagonal tensorsλ h ,λ v is reached. Then the iPEPS represented by tensors {Γ,λ h ,λ v } is in its canonical form and satisfies the conditions defined in Eqs. [17][18][19][20].

B. Fixing gauge degrees of freedom in an iPEPS
One of the most significant benefits of working with the canonical form of an iPEPS is that we can fix its gauge degrees of freedom that can reduce many numerical difficulties. Even if an iPEPS is in canonical form, there is still a gauge freedom to choose the complex phases of each index of the tensors. Here, by assuming that two different iPEPSs describing a same quantum state are in the canonical forms, we suggest an iterative method to fix the gauge degrees of freedom between them. In particular, let us call these two iPEPSs as |Ψ A and |Ψ B which are characterized by tensors {Γ A , λ Ah , λ Av } and {Γ B , λ Bh , λ Bv }, respectively. We define two new tensors A = λ tensors also characterize |Ψ A and |Ψ B and are different from each other by phases or unitary matrices. More precisely, they are related by A = U † h U † v BU h U v where we assume that U h and U v are unitary matrices. Fixing the gauge degrees of freedom between |Ψ A and |Ψ B requires determining matrices U h and U v . By generalizing the iterative method proposed in Sec. III B 2, we can find U h and U v such that they maximize the local fidelity, defined as, which is graphically represented as .
More precisely, U h and U v are obtained iteratively as follows. To start we initialize them arbitrarily, for example we could start from identity matrices (or some random matrices), and repeat following steps: We then obtain (v) Compute the fidelity F and check whether F = 1, up to some acceptable small tolerance. If not, go back to step (i).

V. ENVIRONMENT RECYCLING IN THE TEBD ALGORITHM FOR IMPS (1D)
In this section we explain the environment recycling scheme in the TEBD algorithm for infinite 1D systems. For concreteness, we describe the scheme for the simple quantum Ising model with transverse magnetic field, where the time evolution operator is broken into gates representing the interaction and magnetic field. The infinite system described by this model can be represented by a one-site translationally invariant iMPS. The scheme can be generalized to arbitrary Hamiltonians and larger unit cells.
A. Imaginary-time evolution for the 1D quantum Ising model Consider an infinite spin-1/2 chain described by the Hamiltonian of the quantum ferromagnetic Ising model defined as where σ z , σ x are Pauli matrices, h is the amplitude of the transverse magnetic field. In order to find the ground state of the system employing the iTEBD algorithm, we firstly rewrite the Hamiltonian as where H where we have defined As it is expected to find the ground state of the system which is represented by a one-site translationally invariant iMPS, we need to express the time evolution operator in terms of the matrix product operator (MPO) representation. As illustrated in Fig. 7, at each link between two lattice sites, we employ the SVD to decompose the interaction term as with the same contribution between the interaction and magnetic field terms at each site, we define the MPO as follows, The infinite MPO (iMPO) is then constructed by which is translationally invariant under a single rankfour tensor g of dimensions κ × κ × d × d. To obtain the ground state of the system represented by an iMPS, we apply successively the iMPO to an initial state |Ψ 0 until the convergence is achieved. Mathematically, this is described by where m = T /δ and is the number of iterations in iTEBD algorithm.

B. 1D environment recycling
We now describe the environment recycling scheme and how to incorporate it into the iTEBD algorithm to accelerate the convergence. We first need to determine quantities which characterize the environment in each update of the iTEBD algorithm. We assume that the system is represented by a canonical iMPS |Ψ A0 at time t 0 of the evolution. This iMPS is characterized by a tensor A 0 with bond dimension χ. In the next update, at time t 0 +δ, the iMPS is updated to |Ψ Θ with bond dimension κχ where Θ = gA 0 . The iMPS |Ψ Θ is then transferred into its canonical to perform the bond dimension truncation to obtain a new iMPS |ΨÃ with bond dimension χ. This truncation is equivalent to finding matrices P and Q which dimensions are κχ × χ and χ × χκ respectively such thatÃ = QΘP. (36) Note that the tensors P and Q are the objects that contain the information of the environment that we want to recycle. We assume that the environment characterized by tensors P and Q changes slightly for the next few time steps of the evolution, and thus these tensors can be recycled instead of being re-calculated from the scratch at every iteration. In order to propose the environment recycling scheme, we first need to construct a so-called renormalized gate denoted as G which is obtained by following the fundamental steps: (i) Transfer the iMPS |ΨÃ into its canonical form which is described by tensors {Γ 1 ,λ 1 }. We obtain new iMPS |ΨÃ 1 withÃ 1 = λ 1Γ1 λ 1 . There is also a relation betweenÃ 1 andÃ such that A 1 = Q cÃ P c which is shown graphically as where Q c and P c are χ × χ matrices and can be obtained easily from the relation Q cÃ P c = λ 1Γ1 λ 1 .
(ii) Fix the gauge degrees of freedom between two iMPSs |Ψ A0 and |ΨÃ 1 . This can be done by applying the methods mentioned in Sec. III B to find the complex phases or unitary matrix U . As a result, two tensors A 0 andÃ 1 are related to each other byÃ 1 = U A 0 U † and graphically shown as = (38) (iii) From the tensors {U, Q c , Q, g, P, P c }, we construct a tensor G such that G = U † Q c QgP P c U which is visualized as G is called the renormalized gate.
If we apply the renormalized gate G on the intial iMPS |Ψ A0 , a new iMPS |Ψ A1 is obtained where A 1 has the same bond dimension with A 0 and the gauge degrees of freedom between them are fixed. At this point the environment recycling scheme can be applied by keeping applying the renormalized gate G successively to tensor A 0 for several times to get a new state |Ψ A N before computing a new environment for the next environment recycling process.

C. 1D benchmark results
We have incorporated the environment recycling scheme into the iTEBD algorithm for iMPS to study the ground state of the 1D Ising model with the transverse magnetic field h = 1.05. In Fig. 8, we show the relative errors during the convergence of the ground state energy per link and local magnetization per site versus the number of times that we incorporate the environment recycling scheme into the iTEBD algorithm with respect to different numbers of recyclings N Re 's, denoted as i Re . The total number of time steps is N ts = N Re i Re . Note that unless N Re is extremely large, it may take the same number of time steps N ts to obtain the same ground state of the system for different N Re 's. However, when N Re > 1, the computational effort is already reduced because the number of environment calculations is reduced. Therefore, the larger the N Re , the faster the ground state convergence.
A question naturally arises: can we keep recycling a particular environment as many times as we want? In principle, the answer is yes and the maximum of N Re can be reached when the dominant eigenvector of the renormalized gate is found. Roughly, this can be understood as recycling the environment with an infinite number of times, i.e. N Re → ∞. However, if N Re is chosen to be very large, the scheme can be unstable as the state of the system might fall into the local minima. Therefore, it depends on the problem, and in the general case we need to determine an appropriate value for N Re based on numerical experiments.

VI. ENVIRONMENT RECYCLING IN THE TEBD ALGORITHM FOR IPEPS (2D)
In the previous section we have shown that the environment recycling scheme works quite well for 1D system using iMPS. Although it is not a crucial improvement of the TEBD with iMPS with OBC (recycling the environment in that case does not reduce significantly the computational cost), it is a good start to build the environment recycling scheme for 2D systems. In this section we will present in detail a scheme for 2D infinite system using iPEPS. For concreteness, we propose the scheme to study the ground state of an infinite square lattice system described by the 2D quantum Ising model with a transverse magnetic field.
A. Imaginary-time evolution for the 2D quantum Ising model The ferromagnetic Ising Hamiltonian for a quantum lattice system in 2D has the following form where < r, r > represents the nearest-neighbor sites. For a square lattice, each lattice site has four nearest neighbors. We then group the nearest-neighbor interactions into two terms, for interactions along the horizontal and vertical directions as follows, where we have used indides i and j to represent horizontal and vertical directions in the square lattice respectively, i.e. r = (i, j). The Hamiltonian can be also written as where H In order to obtain the ground state of the system we again employ the imaginary-time evolution. To start we need to express the imaginary-time evolution operator in terms of the Projected Entangled Pair Operator (PEPO). This can be done by decomposing the imaginary-time evolution operator as follows, where we have defined We construct the PEPO in a similar way with what we did above for the MPO. The scheme is explained diagrammatically in Fig. 9. More concretely, we apply SVD to decompose the interaction components in both horizontal and vertical directions of the lattice as follows, and respectively; see Fig. 9(c). With the equal contributions of the interaction terms (from both directions) and magnetic field term, the PEPO at each lattice site is constructed by F , (46) as illustrated in Fig. 9 (d). In Fig. 9 (e), the iPEPO is then defined as a tensor product of identical local terms which is given by The ground state of system will be obtained by acting on the iPEPS with the iPEPO gate as many times as needed until it is converged. We keep in mind that, as a result of acting the iPEPO gate on the iPEPS at each update, the bond dimension of this iPEPS will also increase by κ. Thus, one needs to perform a truncation so that the bond dimension will remain the same after each update and the obtained iPEPS is also the best approximated ground state. At this point we can take advantage of transferring the iPEPS into its canonical form to perform the truncation. This scheme is also known as the SU 43,44,46,47 , applied here for a one-site translationally invariant system. The advantage of the SU is obvious as the computational cost is reduced essentially as compared to the FU. However, in most of the cases, results obtained from this scheme are less accurate, especially when applied to a critical system where the correlations are long-range. In these cases we need the FU, which involves calculating the environment and will be explained next.

B. Full update for a one-site translationally invariant iPEPS
The FU is typically implemented using two-site translationally invariant iPEPS 28,29,33 . In what follows, we will describe the FU for a one-site translationally invariant iPEPS.
At each time step of the evolution, after applying the iPEPO gate to the iPEPS, we get a new state |Ψ g = g ∞ |Ψ A which is represented by tensor Θ with bond dimension κD. We need to find four optimal tensors {P h , Q h , P v , Q v } to insert into four bonds of the tensor Θ simultaneously such that the updated iPEPS |ΨÃ represented byÃ = Q h Q v ΘP h P v with bond dimension D is well-approximated to |Ψ g and this can be achieved by employing the FU.
The FU involves two following fundamental processes: 1. Environment calculation: There are a lot of feasible ways to compute the environment. For instance, one can utilize the coarse-graining tensor renormalization group methods, like TRG/SRG and HOTRG/HOSRG [35][36][37][38] . Alternatively, the environment can be also computed by means of the iTEBD algorithm applied to a boundary iMPS 28 . Here we apply the Corner Transfer Matrix method (CTM) [29][30][31][32][33][34] to compute the environment. This method is known to be more efficient than the MPS schemes when studying the system in the vicinity of the phase transition 29 . In our case, an infinite square lattice L(a), where a (l1l2),(r1r2),(d1d2),(u1u2) = s Θ s l1r1d1u1 (Θ s l2r2d2u2 ) * , is constructed by contracting the physical indices of |Ψ g and Ψ g |. We decompose L(a) into two contiguous regions A and B; see Fig. 10(a). B is chosen to contain four sites where the tensors are connected to each other by the links denoted as r (right), d (down), l (left) and u (up). A plays a role as an environment around B and is required to be computed. Applying the CTM method, this environment is computed approximately and characterized by a set of four corner tensors {C 1 , C 2 , C 3 , C 4 } and eight edge tensors {T a 1 , T a 2 , T a 3 , T a 4 , T b 1 , T b 2 , T b 3 , T b 4 }; see Fig. 10(a). Note that we can further contract the network to obtain the environment represented by only six tensors {E 1 , E 2 , E 3 , E 4 , E 5 , E 6 } around a specific link; see Fig. 10(b).

FIG. 10.
Compute the environment around a specific link using CTM method.
(a) An infinite square lattice L(a) is decomposed into two contiguous parts A and B.
In order to obtain optimal tensors {P h , Q h , P v , Q v }, we need to variationally optimize the pairs of tensors (P r , Q r ), (P d , Q d ), (P l , Q l ) and (P u , Q u ) so that they optimally truncate the bond dimensions of the tensor Θ corresponding to the links r, d, l and u, respectively. This can be done via the standard variational method explained below.

Standard variational method: Let us define a cost
function as follows where ζ = (P r , Q r , P d , Q d , P l , Q l , P u , Q u ). Note that f (ζ) is a quadratic function for every variables in ζ, and hence we can use standard variational method to find ζ that minimizes f (ζ). More concretely, let us fix all tensors in ζ except for the first two tensors (P r , Q r ). We then find P r and Q r that minimize the cost function by following the steps: (i) Fix Q r with some arbitrary initial tensor or the tensor obtained from previous iteration and find P r . Then we can write Eq. (48) as a quadratic scalar expression as follows, where P r is understood as a reshaped vector with D 2 κ components and matrices R, S, T can be obtained from contracting appropriate TNs including the contraction of the whole environment around the link r.
(ii) From Eq. (49), the minimum of f (P r , P † r ) with respect to P † r is obtained only if P r = R −1 S. (iii) Fix the tensor P r and find Q r with the same procedure as steps (i) and (ii) above.
Update the new tensors P r and Q r in ζ and then turn into the links d, l and u to update (P d , Q d ), (P l , Q l ) and (P u , Q u ), respectively. The above process are iterated until all the tensors in ζ converge to fixed tensors. The convergence can be recognized by tracing of the cost function computed in the form of Eq. (48). Note also that instead of applying the optimization scheme described above one could also apply other methods such as the conjugate gradient algorithm to optimize these tensors.
Once the tensors (P r , Q r , P d , Q d , P l , Q l , P u , Q u ) are converged, we can set P h = (P r + P l )/2; Q h = (Q r + Q l )/2;P v = (P d + P u )/2 and Q v = (Q d + Q u )/2. This is due to the links r and l are equivalent and so are the links d and u. The new iPEPS is represented bỹ If the state is evolved by a large number of time steps, where at each time step one apply the iPEPO to the iPEPS combined with the bond dimension truncation using the above FU scheme, a converged iPEPS is obtained which represents approximately the ground state of the system.
Note that the computational cost of this FU is quite expensive, as the cost roughly scales as O(χ 3 D 6 κ 6 ). This dominant cost is due to the environment calculation. Compared to the FU for a two-site translationally invariant iPEPS, where the dominant cost is O(χ 3 D 6 ) as stated in Ref. 29, this one-site translationally invariant scheme is more expensive by a factor κ 6 , and hence, is much less efficient. Fortunately, we can get rid of the factor κ 6 by truncating the tensor Θ before sandwiching it to tensor a to calculate the environment. The truncation can be performed by using tensors {P h , Q h , P v , Q v } which can be obtained from the previous update or the SU. In our simulation we use the latter and it works very well. With this optimization, we reduce the computational cost down to O(χ 3 D 6 ) which is the same as for the case of using two-site translationally invariant iPEPS.

C. 2D environment recycling scheme
The above FU provides a method to study the ground state of a one-site translationally invariant iPEPS, with the same computational cost as the usual two-site algorithm. However, this is still computationally very expensive and more importantly it requires a long time to converge to the ground state. The main aim of this paper is to show how we can further improve this drawback. This can be done by incorporating the environment recycling scheme into the iTEBD algorithm for iPEPS to speed up the convergence.
The 2D recycling environment scheme is constructed similarly with the case of 1D. More precisely, suppose that at time t 0 of the imaginary-time evolution, the iPEPS |Ψ A0 is in the canonical form and characterized by tensor A 0 with bond dimension D. At time t 0 + δ we apply the iPEPO to this iPEPS by contracting tensor g with A 0 . This results to tensor Θ of bond dimension κD. To keep the bond dimension of the iPEPS fixed at D we employ the FU described above to get a new tensor (50) Note that the tensors {Q h , Q v , P h , P v } indirectly contain the environment. As the environment is assumed to change slightly for the next few time steps of the evolution, we can recycle them. Similar to the case 1D, the recycling environment scheme includes two phases: define the renormalized gate G and successively apply G a few time steps on A 0 . In order to determine G we proceed following fundamental steps: (i) Transfer the iPEPS |ΨÃ into its canonical form which is charaterized by tensors {Γ,λ h ,λ v }. We then assignÃ 1 = λ h λ vΓ λ h λ v for the iPEPS |ΨÃ 1 . Next, we can easily find the matrices {Q hc , Q vc , P hc , P vc } such thatÃ 1 = Q hc Q vcÃ P hc P vc which can be seen graphically as, (ii) Find the unitary matrices to fix the gauge between two tensors A 0 andÃ 1 such thatÃ 1 = (iii) Construct the renormalized gate defined as G = The recycling process is implemented by applying successively the renormalized gate G to state |Ψ A0 for several times to get a new state |Ψ A N . Next, we re-assign state |Ψ A0 = |Ψ A N before applying another environment recycling process.

D. 2D benchmark results
In Fig. 11, we plot the local magnetization m z (h) per site which is defined as where |Ψ h is the ground state pertaining to the magnetic field h. The plot also shows the comparisons between different update schemes: SU and FU. We can see that the result obtained from the FU are much better than the results of SU, especially in the vicinity of the criticality. Results of environment recycling scheme for the 2D quantum Ising model are shown in Fig. 12 using iPEPS with bond dimensions D = 3 and χ = 30 is used for computing the environment with CTM. We observe how fast the algorithm drives the system to the ground state by calculating the energy per link E(h) and magnetization per site m z (h), with different numbers of environment recyclings near the critical point h = 3.05. In both plots we realize that the larger the number of environment recyclings the faster the convergence into ground state we obtain. We show extra plots in Fig. 13 to compare the time consumed when the system is evolved with the same amount of computational time T = N ts δ using different  (3,30). The plots show that the larger the N Re , the smaller the time elapsed to evolve the system up to time T . We also fit the data and observe that the time measured almost decays exponentially with the number of environment recyclings.
Note also that as in the case of 1D environment recycling scheme, it is not always allowed to use a freely large number of recyclings although the larger it is the faster convergence to ground state. This is also applied for the case 2D. The reason is that instability appears if the same environment used to update the state that might be significantly different from the initial state, and hence it is contrary to the assumption that the states between successive time evolution update steps are close enough to each other. Therefore, we need to wisely choose an appropriate number of recyclings to stablize the algorithm when dealing with a specific problem.

VII. CONCLUSIONS
In this paper we have first introduced a canonical form for both MPS (with PBC) in 1D and PEPS in 2D, as the fixed point of an iterative procedure, and used this canonical form as the basis to completely fix the gauge freedom in the MPS and PEPS. Then we have introduce a strategy to recycle the environment during a simulation of imaginary time evolution using the TEBD algorithm with iMPS and iPEPS.
By recycling the environment, we have obtained a significant reduction of computational costs while using the FU approach. The speed-up is accomplished because we only need to compute the environment once and can apply it for several time steps of the imaginary time evolution. To validate the schemes, we have studied the ground state of the quantum Ising model with transverse magnetic field in 1D and 2D and shown that the convergence into ground state is greatly accelerated -the computational time required to evolving the system for a given amount of imaginary time decays abruptly. However, the incorporation of environment recycling scheme into the TEBD algorithm with iPEPS sometimes causes numerical instability. For instance, if the initial state on which the time evoution is applied is very different than the ground state of the system, we find it is best to first use a FU scheme without recycling of the environment for a few time steps, and only use recycling once the many-body state has stopped changing dramatically at each time step. Similarly, the number N Re of time steps during which the environment is recycled needs to be adjusted to ensure that the environment does not become too out-dated.
In summary, we have proposed an environment recycling strategy that significantly accelerates ground state computations using iPEPS. We envisage that this strategy will help boost the application of PEPS to challenging 2D problems, such as the resolution of t-J and Hubbard models 39,40 .