{'fraction_non_alphanumeric': 0.035400994676730066, 'fraction_numerical': 0.0122210281658461, 'mean_word_length': 4.461848211848212, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exploring corner transfer matrices and corner tensors for the classical simulation of quantum lattice systems

In this paper we explore the practical use of the corner transfer matrix and its higher-dimensional generalization, the corner tensor, to develop tensor network algorithms for the classical simulation of quantum lattice systems of infinite size. This exploration is done mainly in one and two spatial dimensions (1d and 2d). We describe a number of numerical algorithms based on corner matri- ces and tensors to approximate different ground state properties of these systems. The proposed methods make also use of matrix product operators and projected entangled pair operators, and naturally preserve spatial symmetries of the system such as translation invariance. In order to assess the validity of our algorithms, we provide preliminary benchmarking calculations for the spin-1/2 quantum Ising model in a transverse field in both 1d and 2d. Our methods are a plausible alternative to other well-established tensor network approaches such as iDMRG and iTEBD in 1d, and iPEPS and TERG in 2d. The computational complexity of the proposed algorithms is also considered and, in 2d, important differences are found depending on the chosen simulation scheme. We also discuss further possibilities, such as 3d quantum lattice systems, periodic boundary conditions, and real time evolution. This discussion leads us to reinterpret the standard iTEBD and iPEPS algorithms in terms of corner transfer matrices and corner tensors. Our paper also offers a perspective on many properties of the corner transfer matrix and its higher-dimensional generalizations in the light of novel tensor network methods.


I. INTRODUCTION
It is commonly accepted that understanding the properties of quantum systems of many particles is one of the most important and challenging problems in condensed matter physics. In this sense, much remains yet to be understood. For instance, and in spite of many efforts, it is still unclear what is the mechanism responsible for high-T c superconductivity in cuprates 1 , not to mention iron-based 2 and organic 3 superconductors, for which this mechanism is also a mistery to a great extent. Many other condensed matter phenomena beyond Landau's paradigm of phase transitions have also proven non-trivial to understand. In this respect, a great interest in exotic (i.e. beyond Landau's paradigm) phases of matter has arised recently. Examples of this are, to name a few, topologically ordered phases (where a pattern of long-range entanglement prevades over the whole system) 4 , quantum spin liquids (phases of matter that do not break any symmetry) 5 , and deconfined quantum criticality (quantum critical points between phases of fundamentally-different symmetries) 6,7 .
The standard approach to understand these systems is based on proposing some simplified, physical model that is believed to mimic the relevant interactions responsible for the observed physics. This is the case of e.g. the Hubbard and t − J models for high-T c cuprate superconductors, as well as quantum Heisenberg antiferromagnets with frustrating interactions for some magnetic materials with a spin liquid ground state 8 . Sometimes we are lucky, and these models are exactly solvable. In practice, this means that one can compute some (if not all) relevant properties analitically. But this is not usually the case and, in spite of their apparent simplicity, these models turn out to have such an outstandingly complex behavior that one needs to rely on alternative approaches. Quantum simulations, as proposed by Feyman 9 , are certainly a possibility. Recent experimental results in this direction using e.g. ultracold atoms in optical lattices are in fact really promising 10 . However, the current technological status does not allow yet to fully understand many interesting systems. Thus, one needs to rely on faithful methods to implement numerical (classical) simulations.
From the point of view of numerical simulation algorithms, there has been increasing interest in recent years in the so-called tensor network methods to simulate strongly correlated systems 11 . In these methods the wave function of the system is described in terms of a network of interconnected tensors (a tensor network). As such, tensor network techniques offer efficient descriptions of quantum many-body states of the system that are based on the amount of entanglement in the wave function. The amount and structure of entanglement is a consequence of the chosen network pattern and the number of parameters in the tensors. The most famous example of a tensor network method is probably the Density Matrix Renormalization Group (DMRG) [12][13][14][15] , introduced by White in 1992, and which has been the technique of reference for the last 20 years to compute ground state properties of 1d quantum lattice systems. Recently, though, many important insights coming from quantum information science have motivated the appearance of a host of new tensor network methods. Nowdays it is easy to get lost in the zoo of names for all these methods, e.g. Time-Evolving Block Decimation (TEBD) 16 , Folding Algorithms 17 , Projected Entagled Pair States (PEPS) 18 , Tensor-Entanglement Renormalization Group (TERG) 19 , Tensor Product Variational Approach 20 , Weighted Graph States 21 , Entanglement Renormalization (ER) 22,23 , String-Bond States 24,25 , Entangled-Plaquette States 26 , Monte Carlo Matrix Product States 27 , Tree Tensor Networks [28][29][30][31] and Continuous Matrix Product States 32 , just to name some of the most recent ones. Each particular method has its own advantages and disadvantages, as well as optimal range of applicability.
The possible variations of all these methods are, in practice, unending. In this respect, we feel that a useful tool still relatively unexplored in this context (albeit with exceptions) is the so-called Corner Transfer Matrix (CTM) 71,72 . This was originally introduced by Baxter in 1968 in the context of classical statistical mechanics and exactly solvable models 73,74 . The CTM has very nice properties, specially regarding its spectrum of eigenvalues, and has been a standard tool to find the exact solution of many classical 2d models such as the hardhexagon and related models 71,72 . However, its practical use goes well beyond analytical solutions, and numerical algorithms can be developed to approximate the properties of 2d classical lattice models based on clever truncations in the eigenvalue spectrum of CTMs. Baxter himself was the first to explore this possibility by means of a variational method 73,74 that was an extension of the so-called Krammers-Wannier approximation 75 . In fact, Baxter's method can be understood as a precursor of DMRG but in the context of classical lattice models (this statement will be made more precise later on). The formal combination of CTMs and DMRG was later on put forward by Nishino and Okunishi in the so-called Corner Transfer Matrix Renormalzation Group (CTMRG) 76,77 . Since then, numerical algorithms using CTMs have been widely used, mostly focusing on classical lattice models. From the point of view of quantum lattice systems, Ref. 76,77 already discussed the possibility of simulating 1d quantum systems by using CTMs and a Suzuki-Trotter decomposition of the evolution operator 78,79 . Nevertheless, we believe that CTMs have not yet been fully exploited as a tool in the context of the novel tensor network methods that are being developed for quantum lattice systems.
Our aim in this paper is to cover in part this gap by exploring the applicability of the CTM to develop algorithms for the simulation of quantum lattice systems. We do this mainly in 1d and 2d, which naturally leads us to consider the generalization of the CTM to higher dimensional systems. Following the convention from previous works 73,74,80 , we call this generalization corner tensor. More specifically, in this paper we describe a number of numerical algorithms based on CTMs and corner tensors to approximate ground state properties of quantum lattice systems of infinite size. The methods that we present rely also on the use of Matrix Product Operators in 1d and Projected Entangled Pair Operators in 2d, and naturally preserve the spatial symmetries of the system, including invariance under translations. This, in turn, is a significant difference with respect to some previous approaches for infinite systems that slightly break translational invariance [35][36][37]70 . Also, we will see that these methods produce, in a natural way, individual tensors for the 'bra' and 'ket' parts of local expectation values, in a way similar to the so-called single-layer picture 81 . This does not seem to have remarkable consequences in 1d but, as we shall discuss, it has some interesting implications for the 2d algorithms when compared to other methods such as e.g. iPEPS 37,70 .
In order to prove the validity of our algorithms we provide preliminary benchmarking calculations for the spin-1/2 quantum Ising model in 1d and 2d. Our methods are roughly comparable in accuracy to other well-established approaches such as the iDMRG 33,34 and iTEBD 35,36 methods in 1d and the iPEPS method and TERG in 2d 19,37,70 , thus offering a possible alternative to them. Moreover, we will also discuss that it is possible in principle to use some of the 2d results for the development of a 3d algorithm in combination with some tensor updates, in the same spirit as the algorithm in Ref. 70 for 2d systems. The computational complexity of all the proposed algorithms is also analyzed and, as we shall see, important differences are found depending on the dimensionality and chosen simulation scheme.
For completeness, we also discuss briefly in Appendix B the case of periodic boundary conditions, as well as real time evolution. This discussion will lead us to a beautiful interpretation of the standard iTEBD and iPEPS algorithms in terms of CTMs and corner tensors. We also believe that the algorithms described in this paper will be useful for the pracical implementation and development of further tensor network methods in the future.
This work is organised as follows: In Sec.II we present several preliminary concepts and generalities. These include notions on tensor networks (Sec.II.A), corner transfer matrices (Sec.II.B,C), Matrix Product States and Projected Entangled Pair States (Sec.II.D), and time evo- lution with Matrix Product Operators and Projected Entangled Pair Operators (Sec.II.E). In Sec.III we present two simple algorithms based on CTMs and corner tensors to approximate ground state properties of infinitesize quantum lattice systems. More elaborated versions of these algorithms are presented in Appendix A. After discussing the general approach in Sec.III.A, we consider 1d systems in Sec.III.B and 2d systems in Sec.III.C. A summary of the proposed methods and their complexities is done in Sec.III.E. In Sec.IV we present preliminary benchmarking numerical calculations for 1d and 2d systems. Sec. V contains our conclusions and final remarks. Moreover, in Appendix B we discuss further possibilities, such as 3d quantum lattice systems, periodic boundary conditions, and real time evolution.

II. PRELIMINARY CONCEPTS
The goal of this section is to introduce some preliminary concepts that we will need in the presentation of the algorithms of the forthcoming sections, as well as to provide some background on relevant topics. We start by quickly reminding some basic concepts on tensor networks and tensor network diagrams. Then, we review the basics of the CTM and its properties, as well as several aspects regarding possible gneralisations. After this, we quickly remind the basic definition of Matrix Product States (MPS) and Projected Entangled Pair States (PEPS), and then review the implementation from Ref. 82 of time evolution operators (both in real and imaginary time) using Matrix Product Operators (MPO) and Projected Entangled Pair Operators (PEPO).

A. Tensor networks and diagrams
For our purposes, a tensor is a multidimensional array of complex numbers. A Tensor Network (TN) is a network of tensors whose indices are connected according to some pattern. This connection of indices is done by summing over all the possible values of common indices between tensors. Summing over an index is also called contracting the index. Summing over all the possible indices of a given TN is called contracting the TN.
Instead of using equations, tensors and TNs are more easily handled by using a diagrammatic notation in terms of tensor network diagrams, see Fig.(1). In these diagrams tensors are represented by shapes, and indices in the tensors are represented by lines emerging from the shapes, see Fig.(1.a,b). A TN is thus represented by a set of shapes interconnected by lines. The lines connecting tensors between each other correspond to contracted indices, whereas lines that do not go from one tensor to another correspond to open indices in the TN, see  a scalar. Notice, though, that the total number of operations that must be done in order to obtain the final result of the contraction depends strongly on the order in which tensors in the TN are contracted, see Fig.(2). To minimize the computational cost of a TN contraction, one must thus optimize over the different possible contraction orderings. There are famous examples of TN in the context of many-body physics. For instance, the partition function of a d-dimensional classical lattice model with nearestneighbour interactions is a TN in d dimensions (see Fig.(3)). Also, for quantum lattice systems, the classes of MPS and PEPS are suitable to describe quantum states of 1d and 2d systems respectively (see Fig.(4.a,c)). Other examples make use of an extra 'holographic' dimension accounting for some renormalization group scale 83 , such as Tree Tensor Networks (TTN) [28][29][30][31] (Fig.(5.a)) and the so-called Multi Scale Entanglement Renormalization Ansatz (MERA) (Fig.(5.b)), which is at the basis of Entanglement Renormalization 22,23 . TNs can also be used to describe operators, such as MPOs in 1d and PEPOs in 2d (Fig.(4.b,d)).

B. Corner Transfer Matrix: fundamentals
The concept of Corner Transfer Matrix (CTM) was introduced by Baxter in the context of exactly solvable models in 2d classical statistical mechanics 73,74 . These can be defined for any planar TN. However, for the sake of simplicity, here we shall assume the case of a 2d TN on a square lattice as in Fig.(6). This TN could represent, for instance, the partition function of some 2d classical lattice model as in Fig.(3) or, as we will consider in Sec.III, the imaginary time evolution of a quantum 1d system. In order to define the CTM one makes the following observation: the contraction of the TN can be obtained by multplying four matrices C 1 , C 2 , C 3 and C 4 , one for each corner (see Fig.(5.a)). Thus, where Z is the scalar resulting from the contraction. Matrices C 1 , C 2 , C 3 and C 4 are the Corner Transfer Matrices of the system. They correspond to the contraction of all the tensors in each one of the four corners of the 2d lattice of tensors. The nomenclature 'transfer matrix' is appropriate, since the CTM 'transfers' a vector in the angular direction around the center of the lattice by an angle of π/2 in our case. To further simplify our explanation, let us assume that the four CTMs are equal, Sometimes it is convenient to define diagonal CTMs C d = P CP −1 . Depending on the symmetries of the system (and thus of C), matrix P may be arbitrary, unitary or orthogonal. Let us assume that there are χ different eigenvalues ν α , with α = 1, 2, . . . , χ. Then, the contraction of the full TN adopts the very simple expression The fact that we choose to name the number of different eigenvalues as χ is made on purpose. There is in fact a direct relation between this parameter and the so-called bond dimension of an MPS. The relation between these two seemingly different parameters will be made clearer later on in Sec.III.B, when considering 1d quantum lattice systems. CTMs are of paramount importance in the context of classical statistical mechanics. They have been used to solve the hard hexagon model, as well as many others 71,72 . But they have also been useful in the context of quantum information, since it is known from long ago that their eigenvalue spectrum can be related to the entanglement spectrum of 1d quantum lattice systems, see e.g. Ref. 84,85 . As we will see, this is a key property in the algorithms that will be further described in Sec.III.B. Also, from a numerical perspective, a variational method to approximate the partition function per site of a 2d classical lattice model was developed by Baxter by truncating in the eigenvalue spectrum of the CTM 73,74 . This idea was later on refined by Nishino and Okunishi, who developed the so-called Corner Transfer Matrix Renormalization Group (CTMRG) 76,77 . As such, CTMRG is an algorithm to approximate properties of 2d classical lattice models with isotropic interactions (and thus a high degree of symmetry), and runs by truncating in the largest eigenvalues in magnitude of the spectrum of matrix C 4 , which is a reduced density matrix of the system (see Fig.(6.b)). This spectrum, in turn, is given by the numbers ν 4 α in Eq. (2), and is thus in one-to-one cor- respondence with the spectrum of the CTM C.
As we shall explain later on, from the point of view of quantum states of 1d quantum lattice systems the numbers ν 4 α are, in fact, the spectrum of eigevalues of the reduced density matrix of half an infinite chain. Or what is the same, the spectrum of Schmidt coefficients of half an infinite chain is given by λ α = ν 2 α . Thus, in a way, we can think of Baxter's variational method as sort of a precursor of the truncation scheme used in DMRG and TEBD, but in the context of 2d classical statistical mechanics 114 . Also, it is well known that these CTM methods work assymptotically in the limit of a system of infinite size. In fact, the closer the method is to convergence, the more faithful are the truncations in the eigenvalues of the corresponding reduced density matrix. The convergence and speed of CTMRG for classical 2d models is also remarkably better than other approaches such as the so-called Transfer Matrix Renormalization Group (TMRG) 86,104,105 . This is understandable, because in many situations away from criticality the largest eigenvalue of the CTM is non-degenerate, and there is a rather big gap to the next lower eigenvalues (as compared to the gap of the row-to-row transfer matrix) 115 . This ensures a fast numerical convergence of CTM methods.

C. Corner Transfer Matrix: generalizations
Since the foundational works by Baxter, there have been several attempts to generalize the CTM in a number of different ways. Here we mention some basic possibilities.
The first obvious generalization is for 2d lattices different from the square one. Of course, the CTM can be defined for any planar lattice as long as corners can be defined. This includes the usual lattices, but also more exotic constructions such as lattices with hyperbolic geometry (or negative curvature), e.g. lattice discretizations of Anti de Sitter (AdS) manifolds 87 .
An attempt to improve the efficiency of the truncations involved in CTMRG was also proposed in Ref. 70 , where a directional version of the method was put forward in the context of 2d iPEPS calculations, and where tensors were no longer real and broke rotational symmetry. In Sec.III.B of this paper we will also propose another plausible way of doing this generalization to non-symmetric tensors in the context of the algorithms that we explain here.
Another natural generalization is the case of periodic boundary conditions. Even if this sounds counterintuitive (since a periodic system does not have any corners), it turns out that the ideas from CTMRG can be generalized to deal with this type of systems as well 88 . We will briefly comment on how one can perform this generalization in Sec.III.D.
Moreover, the CTMRG technique can also be generalized to deal with stochastic models, see Ref. 89 for details. Also, a quantum counterpart of the CTM is the so-called Corner Hamiltonian, which essentially is the logarithm of the CTM. Numerical techniques for Corner Hamiltonians have also been developed, see e.g. Ref. 90 .
Finally, the CTM can also be generalized to higherdimensional systems. The natural corner object now is not a matrix, but a tensor with three indices (for a cubic lattice) which we call the corner tensor 116 , see Fig. (7). In terms of corner tensors C 1 , C 2 , . . . , C 8 , the contraction of a 3d lattice of tensors reads where f is a function that performs the corresponding TN contraction (see Fig.(7.a)). This higher-dimensional generalization is quite obvious but, yet, it has some nontrivial consequences since one is dealing now with tensors instead of matrices. In particular, the corner tensors can no longer be diagonalized (or more precisely, there is no unique eigenvalue decomposition since it depends on how one chooses to split the indices of the tensor). Thus, no expression like the one in Eq.(2) for the 2d case can be obtained in general for 3d. Nevertheless, corner tensors still have interesting spectral properties with respect to eigenvalue/singular value decomposition of bipartitions of their indices. The behavior of these singular values will be the key to define simplified numerical approaches for higher-dimensional systems, as we shall do in Sec.III.C.

D. Matrix Product States and Projected Entangled Pair States
Let us now revise briefly the concepts of Matrix Product States (MPS) and Projected Entangled Pair States (PEPS). There is a vast amount of literature on these two families of states, and we refer the interested reader to it for further details (see e.g. Ref. 11 and references therein).
Consider a quantum many-body system of N particles. The quantum state of the system is is the total Hilbert space of the system and H [r] is the individual Hilbert space of particle r. Let us assume that each particle is modelled by a q-level system, so that dim(H [r] ) = q. Given a local basis |i r for each site r, with i r = 1, 2, . . . , q, the quantum state of the system reads The coefficient c i1i2···i N can be understood as a tensor of N indices with O(q N ) complex coefficients. Thus, this is clearly an inneficient representation of the quantum state because the required number of parameters scales exponentially with the size of the system. In order to find an efficient description, one can consider a decomposition of the above tensor into a MPS for 1d, or a PEPS in 2d. These decompositions are shown in the tensor network diagrams of Fig Both MPS and PEPS offer an efficient description of the quantum state |Ψ for 1d and 2d respectively. Moreover, MPS and PEPS are known to have many interesting (and important) properties. For instance, both of them satisfy the so-called area law scaling of the entanglement entropy [91][92][93][94] , which is a key property of low-energy states of most quantum many-body systems (albeit with notable exceptions, see Refs. 95,96 ). Moreover, it is known that ground states of 1d gapped local Hamiltonians can be efficiently approximated with exponential accuracy by an MPS 97 , and the same holds for thermal states in 2d with PEPS 98 . From the numerical point of view, MPS is the relevant class of states at the heart of efficient methods for 1d systems such as DMRG and TEBD. PEPS is also at the basis of simulation methods for 2d systems such as the finite and infinite PEPS methods 18,37 , as well as TERG 19 .
We wish to remark here a couple of properties of MPS and PEPS. First, for a system of size N , the number of parameters in both families of states scales linearly with N , and polynomially in the bond dimension of the tensors (that is, the number of different values for the connecting indices in the MPS or PEPS tensor network, see Fig.(4)). We call this bond dimension χ in the case of MPS, and D in the case of PEPS. Importantly, this bond dimension can be regarded as a measure of the number of parameters in the TN, but also as a measure of the maximum amount of entanglement that can be handled by the wavefunction (see Ref. 18 ). Second, for systems invariant under translations, it is possible to take the thermodynamic limit N → ∞ and consider a MPS or PEPS for a system of infinite size. This is done by repeating the same unit cell of tensors across the whole lattice [33][34][35][36][37] . This trick is at the basis of methods to study infinitesize systems such as iDMRG, iTEBD and iPEPS. The methods that we shall propose in this paper will be for infinite-size systems, and thus rely on this idea as well.

E. Time evolution with Matrix Product Operators and Projected Entangled Pair Operators
Let us now consider the problem of time evolution. We assume that this evolution is generated by a Hamiltonian H that is the sum of local interaction terms on lattice. For simplicity of the explanation, let us also assume a time-independent Hamiltonian with nearest-neighbour interactions, namely although time-dependent cases and more generic types of intereactions could also be considered (including longrange ones 33,34 ). The real time evolution of a given state |Ψ(0) reads It is equally possible to consider the evolution in imaginary time τ in order to get the ground state |Ψ gs of H, namely where we have assumed that the initial state |Ψ(0) has a non-zero overlap with the ground state |Ψ gs , and the appropriate normalization of the state has been included. As is well known, these two different types of evolutions can be approximated using a TN approach for 1d and 2d systems, both for finite and infinite systems. What is important for us, though, is that the evolution operators e −iHt and e −Hτ can usually be approximated by a sequence of well-behaved MPOs 117 in 1d or PEPOs in 2d, as explained in Ref. 82 . This MPO and PEPO approximation will turn out to be a key step in the algorithms that will be proposed in Sec.III.
Let us now remind the basic steps in obtaining MPO and PEPO descriptions for the evolutions generated by H. Here we simply review some of the results from Ref. 82 , but which are important for our purposes. For concreteness of the explanation, let us imagine that H corresponds to the spin-1/2 ferromagnetic quantum Ising model in a homogeneous transverse field, where σ z and σ x are the usual Pauli matrices, and the sum in the first term is over nearest-neighbours. Let us now assume, for simplicity, the case of a 1d system of N particles with periodic boundary conditions. Without loss of generality, we also consider the evolution in imaginary time under this Hamiltonian for a total time T . The first thing to do, and as also done also in other approaches [35][36][37]70 , is to approximate the whole evolution by breaking the total evolution time T into smaller steps of size δτ 1. Then, we have that the evolution operator U ≡ exp(−HT ) can be written as U = [U (δτ )] m = [exp(−Hδτ )] m , with m ≡ T /δτ . Next, we would like to describe the term U (δτ ) = exp(−Hδτ ) by an MPO. We wish to make it in such a way that the resultant MPO is invariant under translations and symmetric under space inversion, since these are symmetries of the original Hamiltonian, and also favors stabiltity in numerical manipulations. As explained in Ref. 82 , this can be achieved in the following way: first, split H = H z + H x , so that we can perform a Suzuki-Trotter decomsposition 78,79 118 . In this splitting, H z contains all the terms with σ z operators and H x the terms with σ x operators. The term exp(−H x δτ ) is just a tensor product of one-body operators, and its MPO representation is trivial (see Fig.(8.a)). However, the term exp(−H z δτ ) requires a little bit more consideration. In order to evaluate the MPO for this term, we first remember that for any ω. Using this property, and a little bit of algebra, it is easy to arrive at an expression for the MPO. If the operator is expressed as where |i r is the basis of eigenstates of σ z at site r (and same for the j's), then the coefficients c i1i2···i N j1j2···j N are given by the MPO where the non-zero components as explained in Ref. 82 . This MPO is represented in Fig.(8.b). Then, in order to obtain a complete MPO combining the H x and H z terms, we just need to combine the tensors for the two pieces as shown in Fig.(8.b). This MPO description is particularly convenient for numerics, since it is symmetric with respect to space inversion, the tensors are real numbers, and it is also invariant under translations. Notice also that one can easily take the thermodynamic limit N → ∞.
It is a good idea to make the tensors in the MPO description real and as symmetric as possible. This can usually be achieved by a variety of tricks depending on the system considered. For instance, as explained in Ref. 82 , for the antiferromagnetic Heisenberg models on bipartite lattices it may be convenient to perform a sublattice rotation prior to finding the MPO representation in order to make the MPO tensors real. It is a good idea to apply tricks like this also to other models, whenever this is possible.
In 2d it is equally simple to find a PEPO description for the evolution operator. The construction is based on that of the MPO for 1d. Let us imagine that this time we have the same Hamiltonian as in Eq.(8) but in a square lattice. Then, we can decompose the evolution term generated by H z into a sum of two contributions: one for interactions along the rows, and one for interactions along the columns, where H row/col z contains the interactions alog the rows/columns. It is clear that, for each individual row or column, one can find MPO descriptions as the ones described before for the 1d case. Therefore, it is just a matter of combining together all these MPOs in a appropriate way in order to obtain a PEPO for the desired evolution. This is represented in the diagram of Fig.(8.c). This PEPO has again nice properties: it is real, and also is symmetric with respect to space inversions in the two lattice directions.
The construction of MPOs and PEPOs that we have reviewed here can be further generalized to a variety of other models and interactions. Notice that it is of course possible to use other alternative approaches to build MPOs and PEPOs for evolution operators. However, these constructions may not always guarantee the symmetry conditions of the obtained TN such as space inversion and translational invariance. As we shall see, it is important for our purposes that the obtained MPOs and PEPOs have these nice symmetries, since this improves the efficiency of the algorithms that we shall propose in Sec.III.

III. ALGORITHMS
In this section we propose a number of numerical simulation algorithms that allow to compute ground state properties of quantum lattice systems of infinite size. Our main focus is on 1d and 2d, but we also discuss briefly the possibility of 3d systems, periodic boundary conditions, and real time evolution. Our algorithms are valid to deal with systems with a different amount of spatial symmetry. As we shall see, more spatial symmetry means simpler and more efficient simulation methods. In 1d this may not be too relevant, since the number of operations in the proposed methods differ only in subleading and multiplicative constant terms. But in 2d this is crucial since the difference in the leading number of operations turns out to be huge, in fact several orders of magnitude, depending on the amount of symmetry. The algorithms that we present here allow to compute ground state properties. As we shall discuss more quantitatively in Sec.IV, these methods are a possible alternative to other methods for systems of infinite size such as iDMRG, iTEBD, iPEPS and TERG 19,[33][34][35][36][37]70 . Importantly, all the algorithms of this section preserve at every step the invariance under translations.
First we discuss the general approach, and then we describe the details of each algorithm. The cases of 3d and periodic boundary conditions are briefly discussed in Sec.III.D. For an overall view of the complexity of the different the methods, one can jump directly to Sec.III.E.

A. General approach
The general idea is quite simple, and can be understood from the diagrams in Fig.(9) and Fig.(10). The goal is to compute the evolution of some quantum state as driven by some Hamiltonian with local interactions for a total time T . As we have explained in the previous section the time evolution operator, both in real and imaginary time, can be decomposed as the action of MPOs or PEPOs for a given time interval δτ on a given quantum state described by an MPS or PEPS. Thus, the whole evolution of the system can be represented by some initial MPS or PEPS, to which one applies the MPO or PEPO driving the evolution for as many steps m = T /δτ as needed, see Fig.(9). In the case of a ground state calculation the number of steps m is infinite or, in practice, very large until convergence of some relevant quantity is achieved.
The main goal of our algorithms is the efficient approximation of expectation values of local obsevables. For instance, let us consider a one-body operator o 1 . Let us also consider the evolved state |Ψ(T ) at some given time T , as represented in Fig.(10) for an MPS (the case of a PEPS is a straightforward generalization). The expec- where U is the corresponding evolution operator that can be decomposed as U = [U (δτ )] m , with U (δτ ) being approximated by an MPO or PEPO. Thus, e.g. the numerator in Eq. (14) can be represented diagramatically as in Fig.(10). The actual expectation value is the ratio of two contracted TNs like the one in the figure, one with the obserble o 1 for the numerator, and one without the observable for the denominator. For a quantum lattice system in d spatial dimensions, these are (d + 1)dimensional TNs, where the extra vertical dimension is time. Therefore, for an MPS one has to deal with the contraction of a (1 + 1)d TN, whereas in 2d one has the contraction of a (2 + 1)d TN.
All the existing TN methods that compute real and imaginary-time evolutions approximate, in one way or another, a contraction like the one described above or similar. For instance, iTEBD and iPEPS methods com- pute an MPS or PEPS approximation to the dominant eigenvector of some transfer matrix operator defined as the MPO or PEPO driving the evolution [35][36][37]70 . But there are more ways of contracting these TNs. For instance, one could think of 'folding' the TN across the red dashed line in Fig.(10) accompanied by a transversal contraction, as has been done in 1d in Ref. 17 .
Our approach in this paper is to use CTMs and corner tensors to approximately contract TNs like the one in Fig.(10) for different types of systems. This is actually a direct generalization of the ideas of CTMRG in Ref. 76,77 . In our case, though, it is worth stressing a peculiarity inherent to the case of quantum lattice systems: any expectation value is by construction a hermitian expresison. This means that the coresponding TN is symmetric with respect to transposition of the vertical indices and complex conjugation, see Fig.(10). Even if quite obvious, this property needs to be explicitely taken into account in all our algorithms.
The methods that we propose are different ways of approximating these expectation values by using CTMs and corner tensors. This is done in two steps. First, one finds a set of renormalized tensors (e.g. renormalized CTMs) accounting effectivly for the most important correlations in the TN. Second, expectation values are approximated by using these renormalized tensors. In what follows we show how to do this in 1d and 2d.

B. 1d quantum lattice systems
In this paper we consider three algorithms using CTMs to approximate expectation values as in Fig.(10). The first algorithm is called 'simplified one-directional 1d method', and is very efficient. This approach is de-signed for systems such that the MPO is also hermitian with respect to the horizontal indices and is explained here. The other two algorithms are considered in Appendix A.a,b. The second algortithm is the 'full onedirectional 1d method', valid for an MPO without the previous symmetry requirement. The third algorithm is the 'full two-directional 1d method', which is equivalent to the CTMRG algorithm generalized to deal with the MPO considered here without extra symmetry requirements. As we shall see, the notion of 'one-directional' or 'two-directional' for each approach refers to the number of directions in which the lattice is simultaneously expanded at every step.
The relevant parameter in these methods will be χ, which is the rank of the renormalized CTMs at every step. From Eqs.(1,2) and Figs. (12.d,13.d,25.c,26.c,28), one can see that this is also the rank of the the renormalized reduced density matrices of the system. Thus, some of the truncations performed by the algorithms in this section are really truncations in the entanglement spectrum (or, equivalently, the spectrum of Schmidt coefficients) of half an infinite chain 119 The leading scaling of the running time (or complexity) of the three methods is the same, namely O(χ 3 ). This is exactly the same complexity as iDMRG and iTEBD methods. However, the multiplicative corrections to this overall complexity are very different for each case, which implies that the total number of operations is also quite different. The simplfied one-directional 1d method is the less time-consuming, whereas the full two-directional 1d method is the one that takes more time.
The procedure in these algorithms is similar to the directional CTM algorithm from Ref. 70 . Namely, we (i) insert rows and columns in the TN in order to expand its structure, (ii) absorb the inserted tensors by contractions towards the horizontal (x) and vertical (y) directions, or towards the corners, and (iii) renormalize the resultant tensors in some proper way.

Simplified one-directional 1d method
Here we assume that the MPO is invariant under transposition of the horizontal indices plus complex conjugation, see Fig.(11.a). This symmetry is taken into account in order to produce a very efficient algorithm. Thus, the TN must be invariant under the index transposition plus complex conjugation both in the horizontal (x) and vertical (y) directions, see Fig.(11.b). This requirement must be satisfied at all steps in the algorithm. In order to achieve this, we consider a renormalized TN that is specified at every step by one CTM C, a half-column transfer matrix T 1 , a half-row transfer matrix T 2 , and the MPO tensor R, see Fig.(11.b). Notice that tensors T 1 and T 2 can also be interpreted as the tensors of two MPSs respectively in the horizontal and vertical directions (see e.g. Fig.(14)). The algorithm implements what we call x-and y-moves. These are as follows: 1. x-move: (a) Insertion. We insert one new column in the TN, as shown in Fig.(12.a).
(b) Absorption. We absorb the new column towards the left, as indicated in Fig.(12.a,b). At this step we produce the new (unrenormalized) tensors C and T 2 .
(c) Renormalization. This is done by means of the isommetry U as shown in Fig.(12.e). Importantly, this isommetry is obtained from the singular value decomposition ot the CTM C as in Fig.(12.c). As a result of this, we obtain the renormalized CTM C and half-row transfer matrix T 2 , see Fig.(12.e). Notice that U is, in fact, the unitary matrix that diagonalizes the reduced density matrix from Fig.(12.d). Finally, since the system is symmetric with respect to horizontal transposition plus complex conjugation, we can just use the new tensors to obtain the complete TN as shown in the diagram of Fig.(11.b).
2. y-move: (a) Insertion. We insert one new row in the TN, as shown in Fig.(13.a).
(b) Absorption. We absorb the new row towards up, as indicated in Fig.(13.a,b). At this tep we produce the new (unrenormalized) tensors C and T 1 .
(c) Renormalization. This is done by means of the isommetry V as shown in Fig.(13.e). This isommetry is obtained again from the singular value decomposition of the CTM C, see Fig.(13.c). As in the previous case, V is the unitary operator that diagonalizes the reduced density matrix in Fig.(13.d). Finally, we use the new tensors to obtain once again the complete TN as in Fig.(11.b).
The whole algorithm now proceeds by iterating these two moves until convergence of some relevant quantity (e.g. the spectrum of singular values of the CTM). This approach is valid for computing imaginary time evolution (and thus ground states) under an MPO with the required symmetry properties, e.g. the quantum Ising or antiferromagnetic Heisenberg and XX models described in Sec.II.E. The leading number of operations of the metod is O(χ 3 ), where χ is the rank of the CTM C. This approach has a number of advantages: it is extremely efficient in finding ground states (the total number of operations is very low), and also preserves the spatial symmetries of the evolution operator at every step. Let us also stress the fact that the net result of the y-move is, in fact, equivalent to (i) a rotation of the 2d TN by π/2, followed by (ii) an x-move, and then followed by (iii) a rotation of the TN by −π/2. This relation is possible because of the symmetries of the TN in Fig.(11.b).
Let us stress another important and intriguing fact about this algorithm, which is that there are no explicit truncations. That is, the rank χ of the CTM is specified from the very beggining in the initial TN, and does not grow at all throughout the evolution. This is a key difference with other methods such as CTMRG and iTEBD, where the analogous to this rank grows at every step in the algorithm and thus needs to be truncated at every step. Here, though, the situation is more subtle. One can think of the truncation begin implemented implicitely. The initial choice of rank χ for the CTM defines already this implicit truncation, and then χ is preserved all along the algorithm. Even if this is conceptually strange, it is really not a problem for the method: every time we insert a column or row, the number of indices of the CTM proliferates in one direction while it is kept constant in the other, which means that the rank of the CTM can not grow. It is this property of the algorithm what implements an automatic implicit truncation and, thus, we do not need to truncate explicitely at any step. However, in Sec.III.D we will see that this property no longer holds in 2d since the indices of the corresponding corner tensors always proliferate in, at least, two directions, so explicit truncations will always be needed in that case.

Expectation values
Once the algorithm has converged, one is in position of computing expectation values of local operators. For this, one needs to consider the contraction of the TN together with the operator for the observable that one wishes to calculate. The CTM method explained before allows to approximate such a calculation. We show this for the general case in Fig.(14) (up to normalization of the wavefunction 120 ) for one-site operators, two-site operators between nearest-neighbours, and also for twopoint correlation functions. As is well known from MPS methods, all these calculations can also be done in O(χ 3 ) operations.

C. 2d quantum lattice systems
In this section we discuss how to generalize the ideas from the previous section to the 2d case. This generalization is conceptually straightforward, but there are a number of significant differences with respect to the 1d case that are worth stressing. Let us summarize these differences: (i) Since more spatial dimensions come into play, we have more possibilities to effectively expand the relevant TN. Thus, we discuss four possible corner algorithms, as opposed to the three discussed in 1d.
(ii) In any of the approaches the indices in the different tensors proliferate at every step in, at least, two different directions. Thus, this growth needs to be always truncated explicitely, regardless of the approach.
(iii) The spectrums of the corner tensors are no longer directly associated to the spectrums of the reduced density matrices of the system. Nevertheless, these spectrums still carry important information about the correlations in the system.
(iv) The complexity of each approach depends dramatically on the level of symmetry and chosen renormalization scheme. Unlike in the 1d case, where all the methods had the same complexity (namely O(χ 3 )), in 2d there are differences in several orders of magnitude (see Table (I) in Sec.III.E).
(v) Our algorithms produce separate tensors for the 'bra' and 'ket' parts of the TN, see Fig.(15). Thus, the TN is always a positive-definite object by construction. This was also the case in 1d, but was not found to be of special relevance there. However, in 2d this is important since the truncation approach differs from the one used in e.g iPEPS and TERG algorithms. In those algorithms, the 'bra' and 'ket' parts of the TN are approximated simultaneously by some common tensors, which in practice breaks the positivity requirement and may lead to numerical instabilities (see e.g. the analysis in Ref. 99 ). In a way, the approach in this paper is similar to the single-layer methods in Ref. 81 .
Before entering into details of the methods, let us mention that another alternative for a 2d algorithm based on CTMs is actually possible. Namely, one could use some 1d method with CTMs to compute effective environments of a 2d iPEPS, and use it in the context of some tensor update to simulate time evolution (e.g. 'simplified' or 'full' updates 37,100 ). This approach is also valid, and has already been explored 70 . However, the algorithms that we explain in this section are based on a completely different approach.
Here, we first present a 'simplified one-directional 2d method'. As in the 1d case, this simplified method is designed for systems such that the tensor defining the PEPO is hermitian with respect to the two spatial directions independently. The implementation of this algorithm is quite efficient and the steps will be explained in detail. In Appendix A.c, we discuss three different 'full' 2d methods which are valid for PEPO tensors without the previous symmetry requirement. In increasing level of complexity, these approaches are called 'full onedirectional, full two-directional, and full three-directional 2d methods'. We only sketch briefly the main idea behind them, and from here the interested reader can infer very easily their step-by-step implementation. As we shall see in the Appendix, these algorithms are less efficient than the simplified directional 2d method, and the complexity for each one of them is also very different. In fact, sometimes the complexity may be higher than the calculation of expectation values itself, as we will see in Sec.III.E. Still, we believe that it is important (at least from a conceptual perspective) to be aware of the existence of all these possibilities.
For 2d quantum lattice systems, the CTMs are generalized to corner tensors (see Fig.(7)), which are tensors with at least three indices, one for each direction of the lattice. Thus, different singular value decompositions of corner tenors may have different ranks. In practice, one can work with an effective rank χ for all the decompositions, and this is the relevant parameter for the 2d methods. However, it is also possible to work with separate ranks, e.g. χ for the vertical direction and D for the horizontal directions. This choice of ranks would be analogous to the iPEPS algorithm, where we have an iPEPS of bond dimension D for which an effective environment with tensors of bond dimension χ is computed (in our approach, we have that χ ≈ χ 2 very roughly, since we produce 'bra' and 'ket' tensors independently). For simplicity, in this paper we choose always the same rank when explaining the details of the methods (including complexity issues). However, we sometimes choose different ranks for the different directions in the numerical calculations of Sec.IV. Let us also remark that, again and as in the 1d case, the truncations in the vertical indices for the full methods are in fact in the entanglement spectrum of a corner of the 2d quantum lattice system, e.g. see again Fig.(7).
As in the 1d case, the procedure for these methods is always the same: (i) insertion of tensors, (ii) absorption of the inserted tensors towards some direction or towards some corner, and (iii) renormalization by means of some isommetry or rectangular matrix. Also, in this section we use the following notation: the spatial (horizontal) directions are called x and y, whereas the temporal (vertical) direction is called z (see Fig.(15)). Also, x-indices refer to indices that connect the tensors in the TN in the x direction (and similar definitions apply for y and z-indices).

Simplified one-directional 2d method
Here we assume that the PEPO is invariant under transposition of the x and y indices independently plus complex conjugation, see Fig.(15.a). As in the 1d case, this extra symmetry is taken into account in order to produce a very efficient algorithm. Thus, the TN must be hermitian in the three directions x, y and z, see Fig.(15.b). This requirement will be satisfied at all steps in the algorithm. Therefore, the renormalized TN is specified at every step by one corner tensor C (analogue to the CTM in 1d), three tensors T 1 , T 2 and T 3 (analogue to the half-row and half-column transfer matrices in 1d), and three tensors X, Y and Z (coresponding to some renormalized iPEPS for each one of the three planes yz, zx and xy), see Fig.(15.b) The algorithm now implements what we call x-, y-and z-moves. As in 1d, these moves are equivalent up to rotations of the whole lattice by π/2. Thus, for simplicity we explain in detail e.g. the x-move, and then explain how the y-and z-moves can be related to the x-move by rotations.
1. x-move: (a) Insertion. We insert one new plane of tensors in y and z directions of the TN, as shown in Fig.(16.a).
(b) Absorption. We absorb the tensors from the new plane towards the left hand side in the x direction, as indicated in Fig.(16.a). At this step we produce the new (unrenormalized) tensors C, T 1 , T 3 and X, see Fig.(16.b).
(c) Renormalization. This is done by the isommetries U y , U z , W y and W z as shown in Fig.(16.d). These are found as follows: first, we find the isommetries U y , U z , W y and W z from the singular value decompositions of tensors C, T 1 and T 3 that are shown in Fig.(16.b). Then, we perform an explicit truncation in the χ largest singular values respectively of all these decompositions, and find the isom- metries U y , U z , W y and W z . With this, we obtain the renormalized tensors C , T 1 , T 3 and X , see 2. y-move: rotate the TN by π/2 in the xy plane as shown in Fig.(17), and do an x-move. Then, rotate the TN back to the original position.
3. z-move: rotate the TN by π/2 in the zx plane as shown in Fig.(17), and do an x-move. Then, rotate the TN back to the original position.
The algorithm works again by iteration of the above three steps until convergence of e.g. the singular value spectrums of the corner tensors. As in the 1d case, this approach is valid for computing ground states by doing imaginary time evolution driven by a PEPO with the required symmetry properties, e.g. the quantum Ising, Heisenberg and XX models described in Sec.II.E. The leading number of operations in this method is O(χ 7 ) if the different contractions are done by making use of the geometric structure of the TN (remember the example from Fig.(2)). As in 1d, this approach has the advantage of being quite efficient, and also keeps the translational invariance of the evolution operator at every step.
Several remarks are in order. First notice that, as hinted previously, the indices in the tensors proliferate in two different directions at each move, and thus one needs explicit truncations in χ. But second, and also unlike in 1d, this time the different truncations are not associated to the truncations in any reduced density matrix of the system. We thus can think of this method as an over-simplified algorithm that implements some 'acceptable' truncation scheme in a very efficient way. We believe that the isommetries found in this way, despite not being the best possible, are still good choices as long as the amount of entanglement in the system is not too large (in a way similar to the simplified update 100 ). In the end, whether this approach is useful or not can only be assessed by numerical simulations, and this is what we will do in Sec.IV. Notice also that these truncations do not correspond to truncations in any entanglement spectrum of the system. For this, we should construct explicitely the reduced density matrices of the system, thus leading to less efficient but probably more accurate algorithms. This is precisely what we discuss in the algorithms of Appendix A.c.

Expectation values
As in 1d, once the algorithm has converged it is possible to compute expectation values of local operators. This is shown for the general case in Fig.(18) (again up to normalization of the wavefunction) for one-site and twosite operators between nearest neighbours as well as for a two-point function along the x direction. This time the required calculations can all be done in O(χ 11 ) operations by choosing the appropriate order in the TN contraction of the TN. The fact that the complexity is O(χ 11 ) means that, for the simplfied one-directional 2d method, this is actually the bottleneck of the calculation, whereas for the rest of the methods explained in Appendix A.c the bottleneck is the calculation of the reduced density matrices in Fig.(30).

D. Summary of methods and complexities
In the previous sections we have discussed two different approaches to simulate quantum lattice systems with corner transfer matrices and corner tensors. The first approach is the one from Ref. 70 : use a (d − 1)-dimensional corner method to approximate effective environents in the context of an iPEPS algorithm in d dimensions. This is the approach that we discussed in Sec.III.D for 3d quantum lattice systems. The other approach, which is the one that we presented in detail here, is to implement directly a corner method in (d + 1) dimensions to approximate the contraction involved in the calculation of expectation values of local observables for quantum lattice systems in d dimensions, assuming an evolution driven by some suitable MPO or PEPO. This has been done in Sec.III.B and Sec.III.C, as well as Appendix A.
The complexity in the algorithms presented in Sec.III.B, Sec.III.C and Appendix A is summarized in Table I. In 1d, all the methods that we have studied here have the same complexity, which is also the same as in iTEBD and iDMRG. However, in 2d the complexity depends on the chosen method. As expected, 2d methods are harder to implement numerically than 1d methods. Nevertheless, these methods can also be implemented in practice within some limitations 121 .

Method
1d 2d and Appendix A, defined as the leding term in the number of operations. The calculation of expectation values is also added for comparison. In 1d the three complexities are the same, but the total running time for each algorithm is different because of subleading corrections and constant multiplicative terms.

IV. BENCHMARK: 1d AND 2d
In what follows we present preliminary numerical results for some of the algorithms discussed previously. Specifically, we have considered the calculation of some of the ground state properties of the ferromagnetic quantum Ising model in transverse magnetic field h from Eq.(8), both in 1d and 2d. The properties of this model are well known 109,110 , and hence it is useful as a first benchmark  for our methods.

A. 1d
In 1d the algorithm in Sec.III.B and the ones in Appendix A.a,b seem to produce similar results for this model. In Fig.(19) we show a calculation of the entanglement spectrum of half an infinite chain for different values of h (focusing on the first 130 spectral values). This is computed from the converged spectrum of eigenvalues of the CTM as in Eq.(2). The calculated spectrums coincide with remarkable accuracy with the ones in the literature (see e.g. Ref. 35,36 ), and were obtained in just a few minutes with very modest computational resources.
We have also computed the two-point correlation function C xx (l), defined as where we have substracted the long-distance term m 2 x , with m x the expectation value of σ x at one site. This correlation function can be computed exactly 109 , and at criticality (h = 1) it tends to decay algebraically with the separation distance l as for some constant a and large l. In Fig.(20) we plot the calculated value of this correlator for χ = 50 and h = 1. We see that our simulation is indeed able to reproduce the assymptotic behavior in Eq.(16) for a large number of sites even for this small value of χ.

B. 2d
In 2d we have computed ground state properties of the Hamitonian in Eq.(8) by using the simplified onedirectional 2d method from Sec.III.C.1. In particular, we have computed the ground state energy per site e 0 , as well as the magnetizations per site m x and m z , defined respectively as the expectation values of the σ x and σ z operators at a given site. For completeness and comparison to the 1d case, we have also computed the entanglement spectrum of a corner of the 2d lattice, as well as of half an infinite plane. In our numerical calculations we employed different truncation parameters depending on the direction. Hence, we used a truncation parameter D for the indices of tensor Z in Fig.(15.b) in the xy plane, and a parameter χ everywhere else.
In Fig.(21) we plot the entanglement spectrum for different values of h both for a corner of the infinite plane, as well as for half an infinite plane, as computed with (D, χ) = (4, 4). These have been calculated from the corner tensors as in e.g. Fig.(7). The finite value of χ used in our simulations truncates these spectrums in 4 and 16 values respectively. Our calculations show that the entanglement spectrum flattens as the critical point is approached, thus involving a much larger amount of entanglement. Interestingly, we also see that long tails in the spectrums tend to be produced close to criticality. When moving slightly away from criticality, we find a crossover region to a regime wehere the spectrums tend to decay very quickly.
In Fig.(22) we show the energy per site e 0 and the magnetization m x as a function of the magnetic field h for (D, χ) = (4, 4). They follow the usual behavior for this 2d quantum system found by other methods (see e.g. Fig.(4) in Ref. 37 ) with good accuracy. In Fig.(23) we show the behavior of the magnetization m z , which is the order parameter, for (D, χ) = (4, 2), (4,4) and (4,6). As expected, close to criticality the order parameter goes to zero according to a critical exponent β as for some constant b. The estimated values of β as well as our estimations of the critical point for these values of the truncation parameters are shown in Fig.(24). Notice that for (D, χ) = (4, 2) the value of the critical exponent is very close to the mean field one β M F = 1/2, but this gets closer to the known Montecarlo result β M C ∼ 0.327 110 as higher χ is considered. A comparison of the critical exponents and critical points as computed by different methods is provided in Table II. We see that this approach, even if quite simple, already obtains a value for the critical exponent that is compatible with the iPEPS approach using CTMs from Ref. 70 , and a value for the critical point already better than the Vertical Density Matrix Approach (VDMA) from Ref. 38,39 . Let us stress that our 2d results are based on the simplest possible 2d algorithm from this paper. At this point, we wish to remind that this algorithm is sort of an over-simplified method because the truncating isommetries are not computed from reduced density matrices at all. Thus, much better accuracies in the critical properties of the system are expected if instead one implements some of the 'full' methods from Appendix A.c. Therefore, the numerical results in this paper should be considered The calculation is for the simplified one-directional 2d method with truncation parameters (D, χ) = (4, 2), (4,4) and (4,6). The inset shows the behavior around the critical region. only as indicative, and not as representative of the best possible performance achievable by corner methods for 2d problems. Nevertheless, we feel that it is quite encouraging that such an over-simplified approach is already able to capture the essential properties of the system within some accuracy. The numerical performance of some of the full 2d methods from Appendix A.c will be considered in a future work 101 .

V. CONCLUSIONS AND FINAL REMARKS
In this paper we have explored the practical use of CTMs and corner tensors to develop tensor network algorithms for the classical simulation of quantum lattice systems of infinite size. At every renormalization step, these methods try to minimize the trucation error by either (i) keeping the largest singular values of a CTM  or corner tensor, or (ii) keeping the largest magnitude eigenvalues of some reduced density matrix. The second scheme maximizes the fidelity between the unrenormalized and renormalized reduced density matrices. In some cases, e.g. 1d systems with some extra symmetries, we have seen that the truncation in (i) is equivalent to that in (ii).
We have focused mainly on ground states of 1d and 2d systems, although we have also discussed briefly other possibilities (3d systems, periodic boundary conditions, and real time evolution). The methods that have been proposed preserve the spatial symmetries of the system, including invariance under translations. We have benchmarked some of this methods by numerically computing several ground state properties of the ferromagnetic spin-1/2 quantum Ising model in 1d and 2d. These numerics, which should be regarded as preliminary, are already quite encouraging. The algorithms of this paper could be a possible alternative to other well-established ways to compute ground state properties of quantum lattice systems in the thermodynamic limit, such as iDMRG and iTEBD in 1d, and iPEPS and TERG in 2d. The computational complexity of the proposed algorithms has also been analized, and we have seen that in 2d it depends strongly on the simulation scheme, leading to differences in several orders of magnitude.
From a broader perspective, it would be interesting to understand the differences between explicit and implicit truncations in the rank of the CTMs. In particular, it would be good to know whether one type of truncation is, by construction, more accurate than the other. In a way, this could be similar to the difference between conventional DMRG and single-site DMRG 112 , where singlesite DMRG seems to produce more accurate data with the same computational resources.
The methods explored in this paper can also be generalized easily to deal with invariance under translations every two (or more) lattice sites. In such a case, one just needs to choose the tensors accordingly in such a way that everything is compatible with the unit cell of the lattice. The different renormalizations need also to be implemented in accordance with the lattice periodicity. Finally, let us mention that these algorithms can also be used in the context of further tensor network generalizations in order to study e.g. systems with internal symetries [42][43][44][45][46][47][48][49][50] and fermionic quantum lattice systems [53][54][55][56][57][58][59][60] . We believe that the methods of this paper will be useful for the practical implementation and development of further tensor network algorithms in the future.
Acknowledgements R. O. acknowledges T. Nishino for a critical reading of this manuscript as well as for providing many insightful comments and important references. EU is also acknowledged for funding through a Marie Curie IIF. In this algorithm the only existing symmetry in the TN is the vertical transposition of the tensors plus complex conjugate, and also translation invariance. Thus, this time the relevant TN is represented by two CTMs C 1 and C 2 , two half-row transfer matrices T 2 and T 3 , and one half-column transfer matrix T 1 , see Fig.(11.c). Notice that tensors T 1 , T 2 and T 3 can again be interpreted as the tensors of three MPS respectively in the horizontal and vertical directions (in the vertical direction there is one for the left hand side, and one for the right hand side). The general idea now is the same as in the previous algorithm: insert, absorb and renormalize. However, since there are less symmetries than in the previous method, the renormalization needs to be done in a different way.
Here are the main steps of the algorithm: 1. x-move: (b) Absorption. We absorb one new column towards the left, and one new column towards the right, as shown in Fig.(25.a,b). At this step we produce four new (unrenormalized) tensors C 1 , C 2 , T 2 and T 3 .
(c) Renormalization. This is done by means of two isommetries U and V as shown in Fig.(25.d). This time, these isommetries are computed by calculating the eigenvalue decomposition of the reduced density matrices in Fig.(25.c). We have that U renormalizes tensors C 1 and T 2 , whereas V renormalizes C 2 and T 3 , see Fig.(  (b) Absorption. We absorb the row towards up, as in Fig.(26.a,b). At this tep we produce three new (unrenormalized) tensors C 1 , C 2 and T 1 .
(c) Renormalization. This is done as in Fig.(26.d) by means of a matrix P . This matrix is computed by finding the eigenvalue decomposition of the reduced density matrix in Fig.(26.c). Importantly, this time P P † is not equal to the χ × χ identity matrix, so we must do the renormalization using P and its inverse P −1 . Notice, however, that the rank of the reduced density matrix in Fig.(26.c) is always χ. Thus, P is a rectangular matrix, and P −1 is the pseudoinverse of P . That is, if the singular value decomposition of P is given by P = F SG † , then its pseudoinverse is defined as The method proceeds again by iterating the above xand y-moves for as long as required until e.g. convergence of the eigenvalue spectrums of the reduced density matrices. Importantly, we can apply this algorithm to compute evolutions by MPOs that are not symmetric. Regarding complexity, this algoithm is again O(χ 3 ), with χ the rank of the reduced density matrices. Also, similar conclusions as in the simplified one-directional 1d method apply for parameter χ in this algorithm: it does not grow throughout the evolution and is fixed from the beginning of the algorithm. Finally, in the case of having an MPO invariant under horizontal transposition plus complex conjugation, this algorithm is just equivalent to the simplfied directional method, although with a less efficient implementation. The following algorithm does no longer use a directional approach with x-and y-moves. Instead, it uses a complete 'radial'-move, expanding both directions of the TN at the same time. It is in fact the CTMRG method 76,77 but adapted to our case (time evolution of a 1d quantum lattice system driven by an MPO). As in the previous section, the system is assumed not to have any spatial symmetry in the horizontal direction (apart from translation invariance). Thus, the relevant TN is described again by two CTMs C 1 and C 2 , two half-row TM T 2 and T 3 , and one half-column TM T 1 , as in Fig.(11.c). Again, the half-row and half-column TMs can be interpreted as defining three different MPSs. The details of the method are as follows: (a) Insertion. We insert two new rows and two new columns, as shown in Fig.(27.a).
(b) Absorption. We absorb the new tensors towards the corners by computing the two (unrenormalized) corner transfer matrices C 1 and C 2 , half-row transfer matrices T 2 and T 3 and half-column transfer matrix T 1 as in Fig.(27.b). (c) Renormalization. Compute the reduced density matrices from Fig.(28), and perform their eigenvalue decomposition. From these decompositions, obtain the two unitary matrices U , V and the matrix P as well as their inverses U † , V † and P −1 . Truncate explicitely in the χ largest eigenvalues in magnitude for each one of these decompositions, and obtain the matrices U, V, P and their pseudoinverses. Then, renormalize the CTMs, half-row transfer matrices and halfcolumn transfer matrix as shown in Fig.(29).
These steps are again repeated until convergence of the spectrum of the reduced density matrices. Notice that this time the rank of the CTM grows at each iteration of the method and hence needs to be truncated explicitely in χ at every step. The leading number of operations for this algorithm is again O(χ 3 ). However, the subleading and multiplicative corrections are higher than in the previous two methods, and hence the total running time is also higher.
c. Full one-, two-and three-directional 2d methods From a generic perspective, the 'full' methods in 2d have several things in common. First, the PEPO does not need to be hermitian in the spatial directions, so we will deal with TNs like the one in Fig.(15.c). Second, the renormalizations are implemented by tensors found from the eigenvalue decomposition of the relevant reduced density matrix of the indices that one wishes to truncate. Third, truncations in the z direction are implemented by isommetries, whereas in the x and y directions one uses rectangular matrices and their pseudoinverses. And fourth, and unlike in 1d, now explicit truncations are always needed.
There are many different ways of puting the above ideas in practice, see Fig.(30). For instance, one could think of a full one-directional approach, similar to the simplified one, but where the tensors that implement the truncations are computed from the reduced density matrices as in Fig.(30.a). Another approach is a full two-directional approach, which keeps the idea of the full two-directional method in 1d but being restricted to the planes in the (2 + 1)d lattice. In this approach, the relevant reduced density matrices are as in Fig.(30.b). In fact, one could even think of combining the one-directional approach in one direction with the twodirectional approach in the other two directions, giving rise to another possible way of expanding the system. Finally, one can think of a full three-directional approach, where the reduced density matrices look as in Fig.(30.c). This last approach is the most direct generalization of the CTMRG algorithm in Ref. 76,77 to higher dimensions. However, this is also the most inneficient alternative. Notice also that in Fig.(30) we consider only reduced density matrices of two indices. However, it would also be possible to consider reduced density matrices of four neighbouring indices in the same direction within a corner in the full two-and three-directional 2d methods (not shown). Such an approach would of course be less efficient, but also more accurate.
We expect stability and numerical accuracy of these methods to increase as the reduced density matrix captures more relevant correlations. But the price that one has to pay is the growth in complexity which, in practice, makes it very difficult to implement numerical simulations in some cases. Assuming calculations of reduced density matrices for two indices (and thus twoindex truncations), we see that the full one-directional approach has complexity O(χ 11 ), the full two-directional approach O(χ 13 ), and the full three-directional approach O(χ 17 ). This is to be compared with the simplified one-directional method from the previous section, which has complexity O(χ 7 ).
These approaches are expected to produce more accurate results than the simplified one-directional approach. However, in this paper we have just implemented numerically the simplfied one-directional 2d approach and seen that this already produces sensible results, see Sec.IV. Thus, we expect these full methods to work even better, despite one has to pay a price in the number of operations required for the different contractions. The numerical exploration of some of these methods will be presented elsewhere 101 .

Appendix B: Further posibilities
We now discuss briefly some other possibilities for corner-inspired methods, namely, 3d quantum lattice systems, periodic boundary conditions, and real time evolution.
d. 3d quantum lattice systems 3d quantum lattice systems (or equivalently, 4d classical lattice systems) are also of relevance for the study of important physical phenomena, e.g. the emergence of fermions and gauge bosons in 3d string-net models 102 , or confinement of quarks in Quantum Chromodynamics. In principle, it should be possible to generalize all the approaches discussed so far in this paper to the 3d case. However, from our experience with the 2d case, we expect that the complexity of the algorithms will be quite large and thus they may become unpractical. Yet, another possibility to deal with a 3d quantum lattice system would be to follow the same idea as in Ref. 70 for the 2d case, namely, to use some 2d method with corner tensors such as the one described in Appendix A.c to approximate the environment of a 3d iPEPS in the context of a 3d iPEPS algorithm. In fact, this would possibly lead to a quite efficient algorithm for studying 3d quantum lattice systems if combined with a simplified update of the tensors 100 . For instance, if the full one-directional 2d method is used, then we would have a 3d method of complexity O(χ 11 ), which is certainly doable. As in Ref. 70 , such an approach may break the traslational invariance by one lattice site of the lattice, but it would be straightforward to generalize our methods accordingly. The numerical exploration of this approach will also be presented elsewhere 101 .

e. Periodic boundary conditions
Indeed, the case of periodic boundary conditions can be considered as well in the context of full algorithms, where the renormalizing matrices are computed by eigenvalue decompositions of the relevant reduced density matrices. This is so since, in the end, periodic boundary FIG. 31: (color online) Examples reduced density matrices for e.g. full one-directional methods with periodic boundary conditions in (a) 1d and (b) 2d. In the 2d case periodic boundary conditions are assumed in one spatial direction only, but these could also be in both spatial directions. If necessary, further periodic boundary conditions in the temporal (vertical) direction could also be considered. Open indices are shown in red.
conditions just add some extra indices to the TN diagram that needs to be contracted, see e.g. Fig.(31). This means that the complexity of the corresponding algorithm will increase (such an increase in complexity is well-known in the DMRG community, see e.g. Ref. 14,15 and also Ref. 40,41 ). From the point of view of the contraction of a 2d TN this provides an alternative to the method in Ref. 103 . Let us also remark that the case of periodic boundary conditions has already been considered in the literature in the context of CTMRG, see e.g. Ref. 88 . Periodic boundary conditions in the temporal direction may also be useful to compute thermal properties, see e.g. Ref. 104-108 . f. Real time evolution: iTEBD and iPEPS reinterpreted So far we have discussed the problem of finding ground states, which is done by imaginary time evolution. However, one may wonder whether the techniques and methods explained here are of any use to compute also the real time evolution of a quantum lattice system, and hence its dynamical properties. We discuss this in what follows.
In a real time evolution, the MPO or PEPO driving the evolution is usually only hermitian with respect to the vertical (time) indices. Thus, if one wished to implement a time-evolution method based on CTMs and corner tensors, 'full' methods (or variations thereof) should be considered. Moreover, in a real time evolution the total running time is finite (unlike in an imaginary time evolution for finding ground states). Therefore, we need to deal with a 'full' approach such that the time direction can grow independently of the spatial directions. The natural options for algorithms are thus the full one-directional approaches from Appendix A.a,b, or in 2d also a combination of the full one-directional 2D approach in the time direction and the full two-directional 2d approach in the space directions (Appendix A.c).
For concreteness let us focus on the 1d case, for which some variation of the full one-directional approach from Appendix A.a may be used. Based on the above considerations, one may be tempted to propose the followiing real-time evolution algorithm: at a given step, (i) repeat many x-moves until convergence (effectively making the size of the system infinite), and (ii) perform one y-move. This can be understood as approximating the action of an MPO over an MPS by means of three tensors, see Fig.(32.a). Then, just repeat these operations for as many (finite) number of steps as required.
The above approach is indeed a possibility. However, numerical simulations show that its accuracy at long times is not as good as the one that can be obtained by using e.g. the standard iTEBD method (results not shown). This is understandable, since the truncations implemented in this approach are more stringent than the ones in iTEBD.
Let us now think of a different approach: we could think again on the action of the MPO over the MPS in terms of CTMs. Then, the aim is to truncate in the horizontal indices by means of some tensor. Following the spirit of the full algorithms described earlier, this can done by computing the reduced density matrix of the indices by means of the CTMs C 1 and C 2 , see Fig.(32.b). Then, we just find its eigenvalue decomposition. With this procedure we obtain a matrix P and its inverse P −1 . A further truncation in the χ largest eigenvalues in magnitude produces the rectangular matrices P and P −1 , which we can use to renormalize the horizontal indices of the tensors.
Remarkably, this procedure is nothing but a reformulation of the the standard iTEBD algorithm for non-unitary evolutions 35,36 . To see this, notice that iTEBD proceeds by (i) finding the canonical form of the resultant MPS after the MPO has been applied, and (ii) truncating its bond index in the largest χ Schmidt coefficients. It is easy to prove that for an MPS in canonical form, the reduced density matrix in Fig.(32.b) is already diagonal. In our laguage, this diagonalization is actually implemented by the matrices P and P −1 . Introducing these tensors in the TN actually orthonormalizes and truncates the bond indices of the evolved MPS. And this is exactly what the iTEBD algorithm does. Also, it is not difficult to imagine that if we generalize the same procedure to 2d systems, then we obtain nothing but a reformulation of the usual iPEPS algortihm but using PEPOs for the time evolution. In such a case, the relevant reduced density matrices can not be computed exactly (since they amount to the contraction of a 2d TN) and must be approximated somehow. Simple approximations lead to efficient algorithms, whereas more elaborate approximations are less efficient but also more accurate.
Let us remark that the reinterpretation that we just found of the iTEBD algorithm in terms of CTMs is in fact similar to the interpretation of DMRG also in terms of CTMs, see Ref. 76,77 .