Criticality and Excitation Gap in Quantum Systems: Applications of Continuous Matrix Product States in Imaginary Time
Abstract
We demonstrate an efficient method that allows for simultaneous determination of the ground state, low energy excitation properties and excitation gap in quantum many body systems. To this aim we first use the abinitio optimization principle of tensor networks (TN), to show that the infinite density matrix renormalization group (iDMRG) in the real space is associated in a natural manner to the infinite timeevolving block decimation (iTEBD) implemented on a continuous matrix product state (MPS), and defined in imaginary time. We illustrate this association showing that the (imaginary) time matrix product state (MPS) in iTEBD reproduces accurately the properties of the twodimensional (2D) classical Ising model, verifying in this way that the time MPS corresponds to a welldefined physical state. We apply then our scheme to the onedimensional (1D) quantum Ising chain, where the time MPS is defined in continuous imaginary time. It is found that the time MPS at or close to the critical point is always maximally entangled, but gives essentially different correlations than those obtained from iDMRG; we argue that these correlations describe the dynamic correlations for the ground state.
pacs:
75.40.Mg, 71.27.+a, 11.25.Hf, 02.70.c=1
I Introduction
Density matrix renormalization group (DMRG) (1); (2); (3) and timeevolving block decimation (TEBD) (4); (5); (6); (7) belong to the most important algorithms in quantum manybody physics. Based on the idea of Wilson’s numeric renormalization group (8), White proposed DMRG in 1992 to simulate the ground states of onedimensional (1D) quantum systems (1); (2). Its basic idea is to add one site to the block and its copy each time, then renormalize the Hilbert space by integrating out the less entangled part of basis. Beyond its original target, DMRG has been extended with great success to simulate finitetemperature thermal states (9); (10); (11) and twodimensional (2D) quantum models of finite size (12), as well as to simulate quantum chemistry problems (13).
Besides the renormalization group (RG) terminology, DMRG is now better understood in the language of matrix product states (MPS) and matrix product operators (MPO) (14). Specifically, the ground state obtained by DMRG was found to be an MPS that indicates a renormalization flow of the Hilbert space. Such a new perspective stimulates us to implement DMRG in a more efficient way and provides new clues for further improvements, e.g., the central matrix technique to speed up the convergence (15), and a perturbation theory of DMRG to improve its accuracy (16). Moreover, MPS has been proven to be extremely powerful as a variational ansatz for 1D quantum models defined on both discrete (14) and continuous space (17); (18); (19).
Different from the RG scheme, TEBD proposed by Vidal in 2004 is a general algorithm based on tensor network (TN) (4); (5); (6); (7). Using TrotterSuzuki decomposition (20); (21), the calculation of the ground state becomes the contraction of a 2D TN, where the TN is contracted layer by layer to an MPS until it is convergent. After each contraction, the dimension of the MPS increases, and an optimal truncation is introduced to reduce the dimension. In this sense, TEBD solves the simulation in a contractionandtruncation (C&T) way.
The C&T scheme is very popular and useful to compute not only ground states but also general TN contraction problems. Compared with the DMRG idea, the C&T scheme is much easier to be generalized to simulate higherdimensional systems. One famous example is tensor renormalization group (TRG) (22); (23); (24); (25); (26); (27); (28); (29); (30); (31); (32); (33); (34). In TRG, the TN is deformed by singular value decomposition (SVD) and then contracted in a certain way so that the resulting TN restores its original geometry. After each contraction, the total number of tensors in the TN decreases and the bond dimensions of the tensors increase exponentially. Then truncations, which can be obtained locally (called simple or cluster update) (22); (23); (24); (27); (29); (31); (32); (33) or nonlocally (called full update) (25); (26); (28); (30); (34), are implemented to bound the dimensions. The main difference among these truncation schemes are the ways of keeping important basis by considering different environments, which lead to different accuracy and computational cost.
Recently, the abinitio optimization principle (AOP) based on TN was proposed (35). Different from the C&T idea, AOP “encodes” an infinite TN into the summation of finite number of degrees of freedom embedded in an entanglement bath that is determined by the fixed point of a set of selfconsistent equations. Interestingly, AOP explicitly shows that for infinite systems, DMRG and TEBD (dubbed as iDMRG (3) and iTEBD (6); (7), respectively) are just two eigensolvers for the same set of selfconsistent eigenvalue equations. Such a unified scheme paves the way to develop new algorithms that inherit the advantages of both and to dig novel concepts to further characterize the physical systems described by the TN.
In this work, we explicitly build an equivalence between iDMRG and iTEBD based on the AOP scheme (35). Considering a 2D TN that can represent a 2D classical partitioning or a (1+1)dimensional quantum theory, we show that when we do iDMRG, which is the renormalization along the spatial direction, in fact we are doing iTEBD in the other direction of the TN simultaneously, and vice versa. Two MPS’s are defined accordingly, which are the wellknown MPS giving a realspace renormalization flow in iDMRG (dubbed as spatial MPS) (14) and a novel translationally invariant MPS in iTEBD (dubbed as time MPS). In a 2D statistical model, we verify our scheme by showing that both MPS’s give accurately the physical properties as expected. In a 1D quantum chain, we find that the time MPS defined in the continuous imaginary time exhibits different properties from the spatial MPS that corresponds to the ground state defined in the discrete space. The time MPS is a maximally entangled state and its correlation is found to be the dynamic correlation of the ground state, from which the excitation gap can be accurately obtained. Our work demonstrates that not only the ground state, but also the properties beyond (such as the excitation gap) can be accurately obtained with our TN scheme.
Ii A brief review of infinite timeevolving block decimation and density matrix renormalization group
In this section, we briefly review the formulation of iDMRG (1); (2); (36) and iTEBD (4); (5). In particular, we explain how to use these two algorithms to simulate the ground state of a 1D quantum system in the thermodynamic limit.
For the iTEBD, it follows the C&T scheme, where the ground state wave function is in an infinite MPS made of a small number of atensors that are repeated indefinitely. To obtain the optimal truncation, one transforms the MPS in its canonical form and truncate the basis according to the canonical Schmidt numbers. In this way, the truncation error is minimized.
Different from iTEBD, the realspace renormalization of the physical Hilbert space is the key in iDMRG. The ground state MPS is not explicitly translationally invariant, and actually represents a renormalization flow.
ii.1 Infinite Matrix Product State
Here, we consider a 1D quantum chain. The physical degrees of freedom on each site give a local dimension Hilbert space , where the local basis are denoted as . A general form of a pure state can is written as
(1) 
with the coefficient matrix. One can see that the size of this matrix increases exponentially with the number of sites. To resolve such an “exponential wall”, it has been proposed to write in an MPS form (37); (14); (38) that reads
(2) 
where is a thirdorder tensor, i.e., a () matrix for each value of ( the bond dimension of the virtual index ). Such an representation can be readily generalized to infinite systems with translationally invariant MPS, i.e. for (see Fig. 1).
Given a quantum state, its MPS representation is in general not unique, but has guage degrees of freedom. For each virtual bond , we can define an invertible square matrix , and rewrite the state Eq. (2) by inserting an identity in each virtual bond as
(3) 
The new MPS is formed by that is again a tensor satisfying
(4) 
The MPS Eq. (2) is then written as
(5) 
which gives exactly the same state as that formed by . Therefore, the gauge of MPS is equivalent to the direct sum of the groups of Isomorphisms of dimensioned complex vector spaces
(6) 
To fix the guage, one can define the leftnormalized form of the MPS, where the local tensor satisfies
(7) 
Similarly, the rightnormalized form is defined as
(8) 
To have the translational invariance, an MPS with a positivedefined diagonal matrix on each virtual bond is introduced as (see Fig. 2)
(9) 
In this way, the canonical form of an MPS can be defined, which is connected with the left and right orthogonal form by
(10) 
In this form, one can prove that is actually the bipartite entanglement spectrum of the state.
ii.2 Infinite density matrix renormalization group
In the DMRG scheme, one typically starts with a short block (of length ) dubbed as the system (S) and its copy as the environment (E). The basis are denoted as and , respectively. Then one grows the chain to add two site between the system and the environment. Then the total state can be written as
(11) 
where and represent the basis of sites added to the system and environment. The aim of the density matrix projection is to determine a subset of states () that optimally approximate the ground state of the enlarged system (environment) block. Accordingly, and are defined by the truncation matrices as
(12)  
(13) 
and are the isometries that realize the truncations of the basis, satisfying Eqs. (7) and (8) (see Fig. 11 (a) and (b)).
Then the ground state with truncated basis are written as
(14) 
To minimizes the truncation error that is indicated by the quadratic cost function
(15) 
the truncation matrices is given by the dominant eigenvectors of the reduced density matrix
(16)  
(17) 
Then with eigenvalue decomposition, one can obtain and as
(18) 
One can easily see that , and in fact give the optimal singular value decomposition (SVD) of in Eq. (11), i.e.
(19) 
showing that gives the entanglement spectrum of the ground state.
It is well known that the ground state obtained by DMRG is actually an MPS formed by , and . From Fig. 3 (c), one can see that the MPS is in an orthogonal form as introduced in the last subsection, and gives a renormalization group (RG) flow of the physical Hilbert space. The direction of the RG flow is determined by the orthogonal conditions shown in Eqs. (7) and (8).
ii.3 Infinite timeevolving block decimation
In this section we will discuss of infinite time evolving block decimation, i.e. we will show how to update the iMPS for the state after an application of operator :
(20) 
where can be real or imaginary. Let us assume that describe an infinite translational invariant system in the form (9), and consists of nearestneighbour interactions only, i.e. , where contains the interaction between sites and . We can then discretize time as with and and use the TrotterSuzuki decomposition (39); (40), which approximates the operator . For example, the first order expansion reads:
(21) 
which contains an error due to the non commutativity of bond Hamiltonians, . The second order expansion similarly reads:
(22) 
where we have to rewrite the Hamiltonian in the following way:
(23) 
All time evolutions on odd and even bonds respectively commute among each other, and can be carried out at the same time.
So we assume is specified by an infinite matrix product operator iMPO (41); (42); (15); (3). This iMPO is represented by a tensor , where and are physical indices and and are bond indices.
(24) 
We can find the MPO of time evolution operator, contracting an MPO for the odd bonds and the MPO for the even bonds. Now we show how construct the MPO for the odd bonds. Let us consider the Trotter step for all odd bonds of a chain:
(25) 
It would therefore be desirable to have in some form containing tensor products , to maintain the MPS form. At this aim we carry out a singular value decomposition:
(26) 
Moreover, the MPO representing the operator in the equation (25) is the following:
(27) 
Then the global MPO can be formed trivially from local MPO. Finally the algorithm occurs in two steps:

Contraction: Firstly we consider the initial state in the form (9). Then we contract the tensor with the tensor of the MPO, producing an evoluted iMPS charactered by the new tensor :
(28) The tensor doesn’t change during the time evolution (because of the translational invariant). The bond dimension of new iMPS became , is larger than the rank of initial iMPS. The computational cost of this step is . Then we brought again the iMPS in canonical form.

Truncation: A final iMPS is obtained from canonical form by truncating all bond indices. In particular, on each bond we preserve the first values of the bond indices, corresponding to the largest Schimidt coefficients.
Iii tensor network encoding and time matrix product state
In Ref. [(35)], the Abinitio optimization principle (AOP) for the ground states of translationally invariant manybody systems was proposed based on TN. The idea is to embed the local subsystem of a unit cell in an entanglement bath that is determined selfconsistently. In the language of TN, the contraction of an infinite TN is encoded into the contraction of the original local tensor with a proper boundary. The boundary tensors are determined by a set of selfconsistent eigenvalue equations.
In this section, we first briefly review the standard AOP where the matrices appearing in the eigenvalue equations are assumed to be Hermitian. Then by employing the idea of DMRG, we modified the eigenvalue equations to get rid of the Hermitian assumption, generalizing AOP for nonHermitian problems. In this modified scheme, we show that two iMPS’s appear along the two directions in the encoding of the 2D TN: one is the iMPS of iDMRG that represents an RG flow, and the other is the translationally invariant iMPS of iTEBD. In other words, the RG transformations in iDMRG are actually the truncations of the iMPS in iTEBD. It means that when iDMRG is implemented in one direction, in fact iTEBD is implemented at the same time in the other direction.
iii.1 Encoding of Hermitian tensor network
The central part of AOP is to solve the selfconsistent eigenvalue equations that determine the boundary tensors (denoted as and ). First, one chooses a cell tensor that is defined as the only inequivalent tensor in the TN. Note there are many choices for a cell tensor, for example, a tensor properly formed by several cell tensors is still a cell tensor.
With the boundary tensors, two local eigenvalue equations are defined by two matrices (assumed to be Hermitian) as
(29)  
(30) 
where is obtained by the SVD of the boundary tensor , i.e. , (Fig. 7). The tensor is introduced from in order to transform the nonlocal generalized eigenvalue problem to a local regular problem [Eq. (31)], where and are required to be Hermitian. The proof can be found in the Appendix in Ref. [(35)]). Then, the boundary tensors and are obtained as the dominant eigenstates of and , respectively, i.e.,
(31)  
(32) 
One can see that such two eigenvalue problems are closely related to each other: one is parametrized by the solution of the other. When the boundary tensors simultaneously solve both equations, i.e. they converge to a selfconsistent fixed point, the whole infinite TN is encoded into the local contraction of the cell tensor with the boundary tensors (Fig. 8).
Specifically speaking, we start with the local contraction (a scalar) given by
(33) 
The eigenvalue equations indicate that is maximized by the boundary tensors. Then, we repeatedly do the substitution with Eq. (31). It should be noted that the normalization of is required here as a constraint. Then Eq. (33) is equivalently transformed as
(34) 
with an infinite MPS formed by the tensor and an infinite MPO formed by .
Then, one utilizes the fact that is the solution of the maximization of . It means that the MPS formed by maximizes Eq. (34), i.e., is the ground state of the MPO. In other words, we have a (nonlocal) eigenvalue equation
(35) 
under the assumption that the ground state of can be effectively represented as an MPS. By substituting Eq. (35) in Eq. (34) repeatedly, the infinite TN formed by can be reconstructed, meaning the whole TN is encoded into the local contraction given by Eq. (33). Note that for the second reconstruction with Eq. (35), the MPS should be normalized. It is a nonlocal constraint that leads to the definition of in Eq. (30).
iii.2 A generalized tensor network encoding
In the formal subsection, the matrices in the eigenvalue equations are required to be Hermitian, so that the infinite TN contraction can be encoded into completely local problems that can be simulated in a more efficient and simple way. In this subsection, we show that iDMRG (3) can be rephrased by the TN encoding scheme, where an equivalence between iDMRG and iTEBD is explicitly built.
Here, we take the onesite iDMRG as example. Comparing with the encoding of Hermitian TN, the general idea of the TN encoding given by iDMRG is the same, where the infinite TN contraction is encoded into local eigenvalue problems. But the details are quite different. As shown in Fig. 9 (a), (b) and (c), there are three local eigenvalue equations, which are given by three matrices
(36)  
(37)  
(38) 
, and are the eigenstate of these three matrices, respectively,
(39)  
(40)  
(41) 
and , which are the left and right orthogonal part of , are obtained by QR decomposition [Fig. 9 (d)] as
(42) 
with and are isometries, satisfying orthogonal conditions as
(43)  
(44) 
The above scheme can be reinterpreted in the iDMRG language. and represent the system and environment superblock. is the effective Hamiltonian with its ground state. By the QR decompositions on , the renormalization of the basis of the system and environment are given by and , and is the center matrix.
Similar to the AOP, the original TN can be reconstructed by repeatedly using the eigenvalue equations. We start from a local contraction (a scalar ) of and the boundary tensors as shown at the bottom of Fig. 10. Then with Eqs. (39) and (40), is transformed to the product of an MPO with an MPS and its conjugate. One can see that like iDMRG, the MPS is formed by ’s on the left side (leftorthogonal part) and ’s on the right side (rightorthogonal part) with the center matrix in the middle. Then by using the fact that such an MPS optimally gives the ground state, we can use the corresponding eigenvalue equation to reconstruct the entire TN (Fig. 10).
In the encoding, there are three constraints: the normalizations of , and the MPS. The first two constraints are obviously local. For the third one, by utilizing the orthogonal conditions of and , the normalization of the MPS is equivalent to that of , which is also local. Thus, all eigenvalue problems are local and regular.
In the above scheme, one can explicitly see that (or ) also gives an MPS, and interestingly, such an MPS is updated in the way of iTEBD contracting along the spatial direction. Let’s pay attention to the eigenvalue equation that satisfies [Eq. (39)]. If one chooses to solve it using a power method, i.e. to find the ground state of by updating with the product . Such a product is equivalent to evolving the local tensor of the MPS in iTEBD with the MPO formed by in the vertical direction. In the evolution, the bond dimension of the local tensor increases exponentially. Then, truncation is implemented by contracting with the isometry .
Note that is obtained from the QR decomposition of . Considering the eigenvalue equation [Eq. (41)], in fact represents one half of the infinite environment of the vertical MPS, thus gives the optimal truncation matrix. All these arguments also apply to and .
Iv Applications
In this section, we apply the TN encoding scheme to 2D classical Ising model and 1D quantum Ising model. There are intrinsic differences between their TN representations. For the 2D classical partitioning, the two dimensions of the TN are both spatial and discrete, and are equivalent to each other. For the 1D quantum chain, one dimension of the TN is spatial and discrete, and the other is the time dimension that is continuous. Though both models give the same central charge that corresponds to a free fermionic field theory in the real space, we find intrinsically different properties in the imaginary time characterized by the temporal correlations and the corresponding “central charge”. It means that the criticality of a ()D theory or a 2D TN should be characterized by two scalings in the two directions of the TN.
iv.1 Twodimensional classical Ising partition function
To demonstrate the validity of our theory, we study the 2D classical model. Its partition function can be directly written in a 2D TN, where the local tensor is the probability distribution of some local Ising spins. Here, we take the Ising model on square lattice as an example, where the local tensor is defined as
(45) 
with the inverse temperature and the spin index . Note that the local tensor of the TN can also be chosen as the contraction of several ’s.
From the exact solution (43), we have the critical temperature . At the critical temperature, the dominant eigenstate of the transfer matrix can be approximated as an MPS. It is well known that each MPS with a finite bond dimension corresponds to a gapped state with a finite correlation length and entanglement entropy . the central charge can be extracted by the scaling behavior (44); (45); (46). Specifically speaking, with different , and satisfy
(46)  
(47) 
The coefficient gives the central charge that corresponds to a free fermionic field theory (44). By combining Eqs. (46) and (47), one can readily have
(48) 
which is independent of the calculation parameters.
In Fig. 11, we show that the MPS from iDMRG has the same correlation length and entanglement entropy as the vertical MPS obtained simultaneous in the iTEBD in the other direction with the difference . In other words, these two MPS’s, though updated within two different schemes and located in two different directions of the TN, are connected by a guage transformation. By choosing different bond dimension cutoff and cells to construct the tensor , the relation between and shows robustly a logarithmic scaling, giving accurately the central charge. The precision increases with the size of the cell tensor.
iv.2 Onedimensional quantum Ising chain
For the ground state of the 1D quantum chain that is essentially a (1+1)D theory, the vertical MPS is defined on the continuous time dimension. Here, we take the 1D quantum Ising chain in a transverse magnetic field as an example, where the Hamiltonian reads
(49) 
with and the spin operators on the x and z directions. The ground state has been exactly solved by fermionization (47) with central charge .
The ground state calculation can be transformed to a TN contraction problem (6) by TrotterSuzuki decomposition (20). The local tensor of the TN is obtained by the imaginarytime evolution operator with the Trotter step an infinitesimal real number. Then, we use the approach explained in Sec. III (B) to solve the TN contraction (similar to iDMRG (3) but applied directly to TN). The ground state is given by the MPS extended in the spatial direction. In this scheme, another welldefined MPS appears simultaneously in the time direction. It is interesting to see the physical properties of the time MPS as well as its scaling behavior at the critical point.
In Fig. 12, the correlation length and entanglement entropy of the spatial MPS converges rapidly to finite values as the Trotter step decreases. At this situation, the Trotter error that is estimated as is negligible.
For the time MPS, its correlation length diverges as decreases, thus is not welldefined. Considering that the (imaginary) time length is determined by both and the Trotter step that characterizes the discretization of time, we redefine the time correlation length as . Fig. 13 (a) shows that converges rapidly when reducing .
It is known that while using MPS to approximate the ground state at the critical point, its entanglement entropy fulfills a logarithmic scaling law versus correlation length and the bond dimension , as shown in Eqs. (46) and (47) (44); (45); (46). As shown in Fig. 14, our results with precisely demonstrate such behaviors.
In Fig. 15, the scaling of against of the spatial MPS with and is given. The results accurately give the central charge with an error . Note that for different choices of , there are small corrections to both and versus , leading together to better results of .
The properties of the time MPS are quite different. For the temporal entanglement (18), the time MPS is always maximally entangled, i.e., the entanglement spectrum is flat. It means the temporal entanglement spectrum exactly satisfies the logarithmic scaling law
(50) 
The reason is that when the Trotter step approaches zero, the evolution of the spatial MPS is nearly identical. Thus as shown in the inset of Fig. 16, the dominant term of the “evolution” of the time MPS is the tensor product with the identity of its “virtual” bonds (which actually are the physical bonds of the spatial MPS). Such temporal entanglement spectrum appears for any magnetic fields, meaning Eq. (50) holds no matter the system is at the critical point or not.
Different from the temporal entanglement, we show that the correlation length of the time MPS is actually the dynamic correlation length of the ground state, whose properties rely on the criticality of the model. Taking the transverse Ising chain as an example, the dimensions of time and space have been proven to be equivalent transforming the model into a 2D classical partition function. Two correlation length satisfy
(51) 
with the rescaling factor in this model. Such a relation is verified at the critical point, as show in Fig. 17.
On the other hand, we find that the time MPS is also critical at , satisfying a scaling law against the bond dimension that is similar to the spatial MPS. In Fig. 16, we show that and satisfy Eqs. (46) with the exponent and (50), respectively. By fitting with Eq. (48), one can still define a “temporal central charge” , which is given in Fig. 18 with different and .
For comparison, we show in Fig. 19 and of the time MPS at when the magnetic field is away from the critical point. One can see that is only and quickly converges to a finite value as increases. On the other hand, still strictly satisfies the logarithmic scaling law given by Eq. (50), meaning the time MPS is maximally entangled just the same as the critical ground state.
It is known that at a noncritical point, the dynamic correlation function decays exponentially as
(52) 
with the excitation gap. Thus, the excitation gap can be given by the dynamic correlation length as
(53) 
As shown in Fig. 20, such a relation is precisely found with our scheme, meaning that the excitation gap can be accurately determined from the time MPS.
Our simulations show that the time MPS’s at or away from the critical point are both maximally entangled, satisfying the logarithmic scaling given by Eq. (50). Though these two time MPS’s have the same entanglement property, it is quite interesting to see that the correlation functions are essentially different. While the entanglement actually gives the properties of the fixed point state (the ground state of the transverse evolution), the correlation length contains the information of the corresponding excitations: the criticality of the evolution determines the scaling behavior of the correlation. But it is still unclear if the scaling coefficient corresponds to a central charge in the conformal field theory.
V Conclusion
In this work, we propose a unified scheme based on the AOP of TN, where the spatial MPS in iDMRG and time MPS in iTEBD are defined in the two directions of the TN. The validity of these two MPS’s are shown by applying to the 2D classical Ising model. Then we demonstrate on the 1D transverse Ising model that the time MPS gives the dynamic correlation of the ground state, from which the excitation gap can be accurately obtained. Our scheme can be readily applied to other 1D models, including those defined in continuous space and time. It can also be generalized to simulate the dynamic correlation as well as the excitation gap of 2D quantum systems.
Acknowledgments
We thank Ian McCulloch for stimulating discussions. This work was supported by ERC ADG OSYRIS, Spanish MINECO (Severo Ochoa grant SEV20150522), FOQUS (FIS201346768), Catalan AGAUR SGR 874, Fundació Cellex, and EU FETPRO QUIC.
References
 S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
 S. R. White, Densitymatrix algorithms for quantum renormalization groups, Phys. Rev. B 48, 10345 (1993).
 I. P. McCulloch, Infinite size density matrix renormalization group, arXiv:0804.2509 .
 G. Vidal, Efficient Classical Simulation of Slightly Entangled Quantum Computations, Phys. Rev. Lett. 91, 147902 (2003).
 G. Vidal, Efficient Simulation of OneDimensional Quantum ManyBody Systems, Phys. Rev. Lett. 93, 040502 (2004).
 G. Vidal, Classical Simulation of InfiniteSize Quantum Lattice Systems in One Spatial Dimension, Phys. Rev. Lett. 98, 070201 (2007).
 R. Orús and G. Vidal, Infinite timeevolving block decimation algorithm beyond unitary evolution, Phys. Rev. B 78, 155117 (2008).
 K. G. Wilson, The renormalization group: Critical phenomena and the Kondo problem, Rev. Mod. Phys. 47, 773 (1975).
 R. J. Bursill and T. Xiang and G. A. Gehring, The density matrix renormalization group for a quantum spin chain at nonzero temperature, J. Phys. Cond. Matter 8, L583 (1996).
 X. Q. Wang and T. Xiang, Transfermatrix densitymatrix renormalizationgroup theory for thermodynamics of onedimensional quantum systems, Phys. Rev. B 56, 5061 (1997).
 A. E. Feiguin and S. R. White, Finitetemperature density matrix renormalization using an enlarged Hilbert space, Phys. Rev. B 72, 220401(R) (2005)
 E. M. Stoudenmire and S. R. White, Studying TwoDimensional Systems with the Density Matrix Renormalization Group., Ann. Rev. Cond. Matter Phys. 3, 111 (2012).
 G. K.L. Chan and A. Keselman and N. Nakatani and Z. D. Li and S. R. White, Matrix Product Operators, Matrix Product States, and ab initio Density Matrix Renormalization Group algorithms, arXiv: 1605.02611.
 U. Schollwök, The densitymatrix renormalization group in the age of matrix product states., Ann. Phys. 96, 326 (2011).
 I. P. McCulloch, From densitymatrix renormalization group to matrix product states., J. Stat. Mech. 2007, P10014 (2007).
 E. Tirrito, S.J. Ran, A. J. Ferris and M. Lewenstein, An efficient perturbation theory of density matrix renormalization group, arXiv: 1605.00940.
 F. Verstraete and J. I. Cirac, Continuous Matrix Product States for Quantum Fields, Phys. Rev. Lett. 104, 190405 (2011).
 M. B. Hastings and R. Mahajan, Connecting entanglement in time and space: Improving the folding algorithm, Phys. Rev. A 91, 032306 (2015).
 V. Stojevic, J. Haegeman, I. P. McCulloch, L. Tagliacozzo and F. Verstraete, Conformal data from finite entanglement scaling, Phys. Rev. B 91, 035120 (2015).
 M. Suzuki and M. Inoue, The STTransformation Approach to Analytic Solutions of Quantum Systems. I General Formulations and Basic Limit Theorems, Prog. Theor. Phys. 78, 787799 (1987).
 M. Inoue and M. Suzuki, The STTransformation Approach to Analytic Solutions of Quantum Systems. II TransferMatrix and Pfaffian Methods, Prog. Theor. Phys. 79, 645 (1998).
 M. Levin and C. P. Nave, Tensor Renormalization Group Approach to TwoDimensional Classical Lattice Models, Phys. Rev. Lett. 99, 120601 (2007).
 H. C. Jiang, Z. Y. Weng, and T. Xiang, Accurate Determination of Tensor Network State of Quantum Lattice Models in Two Dimensions, Phys. Rev. B 101, 090603 (2008).
 Z. C. Gu, M. Levin, and X. G. Wen, Tensorentanglement renormalization group approach as a unified method for symmetry breaking and topological phase transitions, Phys. Rev. B 78, 205116 (2008).
 J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Classical Simulation of InfiniteSize Quantum Lattice Systems in Two Spatial Dimensions, Phys. Rev. Lett. 101, 250602 (2008).
 Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, and T. Xiang, Second Renormalization of TensorNetwork States, Phys. Rev. Lett. 103, 160601 (2009).
 Z. C. Gu and X. G. Wen, Tensorentanglementfiltering renormalization approach and symmetryprotected topological order, Phys. Rev. B 80, 155131 (2009).
 R. Orús and G. Vidal, Simulation of twodimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009).
 L. Wang and F. Verstraete, Cluster update for tensor network states, arXiv:1110.4362.
 Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang and T. Xiang, Coarsegraining renormalization by higherorder singular value decomposition, Phys. Rev. B 86, 045139 (2012).
 W. Li, J. von Delft and T. Xiang, Efficient simulation of infinite tree tensor network states on the Bethe lattice, Phys. Rev. B 86, 195137 (2012).
 P. Czarnik L. Cincio and J. Dziarmaga, Projected entangled pair states at finite temperature: Imaginary time evolution with ancillas, Phys. Rev. B 86, 245101 (2012).
 P. Czarnik and J. Dziarmaga, Projected Entangled Pair States at Finite Temperature: Iterative SelfConsistent Bond Renormalization for Exact Imaginary Time Evolution, Phys. Rev. B 92, 035120 (2015).
 H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz and R. Orús, Infinite projected entangled pair states algorithm improved: Fast full update and gauge fixing, Phys. Rev. B 92 035142 (2015).
 S. J. Ran, Ab initio optimization principle for the ground states of translationally invariant strongly correlated quantum lattice models, Phys. Rev. E 93 053310 (2016).
 A. J. Daley and C. Kollath and U. Schollwöck and G. Vidal, Timedependent densitymatrix renormalizationgroup using adaptive effective Hilbert spaces, J. Stat. Mech.: Theor. Exp., P04005 (2004).
 S. Östlund and S. Rommer, Thermodynamic Limit of Density Matrix Renormalization, Phys. Rev. Lett. 75, 3537 (1995).
 R. Orus, A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States, Ann. Phys. 349, 117 (2014).
 M. Suzuki, Fractal decomposition of exponential operators with applications to manybody theories and Monte Carlo simulations, Phys. Lett. A 146, 319 (1990).
 M. Suzuki, General theory of fractal path integrals with applications to manybody theories and statistical physics, J. Math. Phys. 32, 400 (1991).
 G. M. Crosswhite and A. C. Doherty and G. Vidal, Applying matrix product operators to model systems with longrange interactions, Phys. Rev. B 78, 035116 (2008).
 V. Nebendahl and W. Dür, Improved numerical methods for infinite spin chains with longrange interactions, Phys. Rev. B 87, 075413 (2013).
 L. Onsager, Crystal Statistics. I. A TwoDimensional Model with an OrderDisorder Transition, Phys. Rev. 65, 117 (1944).
 G. Vidal, J. I. Latorre, E. Rico and A. Kitaev, Entanglement in quantum critical phenomena, Phys. Rev. Lett. 90, 227902 (2003).
 L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir and J. I. Latorre, Scaling of entanglement support for matrix product states, Phys. Rev. B 78, 024410 (2008).
 F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore, Theory of FiniteEntanglement Scaling at OneDimensional Quantum Critical Points, Phys. Rev. Lett. 102, 255701 (2009).
 P. Pfeuty, The onedimensional Ising model with a transverse field, Ann. Phys. 57, 7990 (1970).