## Abstract

The density-matrix renormalization group (DMRG) method has established itself over the last decade as the leading method for the simulation of the statics and dynamics of one-dimensional strongly correlated quantum lattice systems. The DMRG is a method that shares features of a renormalization group procedure (which here generates a flow in the space of reduced density operators) and of a variational method that operates on a highly interesting class of quantum states, so-called matrix product states (MPSs). The DMRG method is presented here entirely in the MPS language. While the DMRG generally fails in larger two-dimensional systems, the MPS picture suggests a straightforward generalization to higher dimensions in the framework of tensor network states. The resulting algorithms, however, suffer from difficulties absent in one dimension, apart from a much more unfavourable efficiency, such that their ultimate success remains far from clear at the moment.

## 1. Introduction

In the last decade, the density-matrix renormalization group (DMRG), invented by White [1,2], has established itself as the leading method for the simulation of the statics and dynamics of one-dimensional strongly correlated quantum systems, both at zero and finite temperatures (for more recent reviews, see [3–5]). Nevertheless, attempts have been made to extend this method to the simulation of two-dimensional quantum systems and (three-dimensional) small clusters, such as nuclei or molecules. These attempts are currently being advanced owing to the alternative view of the DMRG as a variational method in the space of matrix product states (MPSs) [6–9] because they allow for an easy extension of the ansatz class to higher dimensions [10,11], which might hold interesting promise for such finite clusters of particles. In this short overview, following [5], I will first introduce and discuss MPSs, introduce the most basic version of the DMRG as a variational method and discuss its dual nature as a variational and renormalization group (RG) method, as well as the reasons for its successes and limitations. I will conclude by a short discussion of the generalizations to higher dimensions.

## 2. Density-matrix renormalization group in one dimension

### (a) Matrix product states

Let us assume that our system consists of *L* sites on a one-dimensional chain with local state spaces {|*σ*_{i}〉} on site *i* of dimension *d* each. Then, it can be shown (e.g. [5]) that any quantum state
2.1
where |** σ**〉=|

*σ*

_{1},

*σ*

_{2},…,

*σ*

_{L}〉, can be brought (for example) by a sequence of singular value decompositions (SVDs) exactly into the form of an MPS, namely 2.2 where on each lattice site

*i*we introduce a set of

*d*matrices

*M*

^{σi}. To see this, let us briefly recall the key properties of an SVD (e.g. [12]): for an arbitrary, rectangular matrix

*A*of dimension (

*m*×

*n*), it provides a decomposition

*A*=

*USV*

^{†}, where the matrices

*U*,

*S*and

*V*

^{†}have dimensions (

*m*×

*p*), (

*p*×

*p*) and (

*p*×

*n*), where , the smaller of the two original dimensions. In addition, the three matrices have the following properties:

*U*consists of

*p*orthonormal columns, which can be read as a set of orthonormal vectors (hence,

*U*

^{†}

*U*=

*I*); similarly,

*V*

^{†}consists of

*p*orthonormal rows, which again can be read as a set of orthonormal vectors (and

*V*

^{†}

*V*=

*I*); they are referred to as the left and right singular vectors.

*S*is a non-negative diagonal matrix, with

*r*<

*p*positive singular values on the diagonal;

*r*is the rank of

*A*. In order to carry out the transformation into an MPS, we start by reshaping

*c*

_{σ}into a matrix

*Ψ*

_{σ1,(σ2 … σL)}=

*c*

_{σ1 … σL}of dimension (

*d*×

*d*

^{L−1}). This is now singular value decomposed as 2.3 The

*U*-matrix has dimension

*d*×

*d*. As the next step, we reshape into a matrix

*Ψ*

_{(a1σ2),(σ3…σL)}of dimension (

*d*

^{2}×

*d*

^{L−2}), where the row dimension may be smaller if singular values of the previous SVD happen to be zero. This is again singular value decomposed as 2.4 with the dimension of the

*U*-matrix being

*d*

^{2}×

*d*

^{2}. If we continue this procedure to the end of the chain, we end up with a sequence of

*U*-matrices as

*U*

_{σ1,a1},

*U*

_{(a1σ2),a2},

*U*

_{(a2σ3),a3}and so on, with the last matrix being

*U*

_{(aL−1σL),1}. The total state then reads 2.5 To simplify the notation, we now define the

*d*matrices

*M*

^{σi}per site as and use that the sum over the

*a*

_{i}is then nothing but a sequence of matrix multiplications of

*M*

^{σi}matrices as in equation (2.2).

It is useful to introduce a graphical notation for MPSs (figure 1), where each matrix is represented by a blob with outgoing legs, namely a vertical physical leg and horizontal legs corresponding to the matrix indices. The rule is that all connected legs are contracted, leading to the MPSs as defined above.

If we review the above sequence of decompositions and the sizes of the *U*-matrices, the resulting *M*^{σ} matrix dimensions will be (1×*d*),(*d*×*d*^{2}),…,(*d*^{L/2−1}×*d*^{L/2}),(*d*^{L/2}×*d*^{L/2−1}),…,(*d*^{2}×*d*),(*d*×1): the *M*-matrices are always cut in the row dimension by a factor of *d* by the reshaping from the *U*-matrices; beyond the centre of the chain, the column dimension will be smaller than the row dimension of the singular value decomposed matrix, hence the dimensions of the *U*-matrices will go down again. In systems of interest, where *d*^{L/2} will be exponentially large in system size, this exact representation is obviously useless in numerical practice. But let us assume that we consider an MPS where matrix dimensions do not exceed some number *D*, which will be *O*(100) through *O*(1000), and that these approximate MPSs are a very close approximation to the exact quantum states of interest. Under this assumption, we may ask about quantum mechanics in this restricted set of MPSs of dimension *D*.

Before we consider typical MPS operations, we note that MPS representations of quantum states are not unique: between any matrix pair *M*^{σi}*M*^{σi+1}, we may insert *I*=*XX*^{−1}, where *X* is an invertible matrix of suitable dimensions, changing and . One of the main usages of this is to bring the MPS into *canonical form*: any MPS can be written in left-canonical form
2.6
where the matrices *A*^{σi} are *left-normalized*, , or in right-canonical form
2.7
where the matrices *B*^{σi} are *right-normalized*, , or in mixed-canonical form
2.8
where the matrix *Ψ* is of matching dimensions; one can also contract e.g. , which generates a variety of ansatz classes for the DMRG. To see how this can be done iteratively, consider a state that is already partially in left-canonical form on its left end,
2.9
Then, one way of making *M*^{σi+1} left-canonical is to reshape the set of matrices into one matrix *M*_{(σi+1ai),ai+1}, which upon an SVD reads . *U* is reshaped into a set of matrices , which meet the left-normalization condition owing to the column orthonormality of *U*. *SV*^{†} is multiplied into *M*^{σi+2}, and the next iterative step occurs at site *i*+2. Similarly, one can start from the right to obtain right-normalized matrices. In fact, the proof of the representation in equation (2.2) generated left-normalized matrices right away, owing to the orthonormality properties of the *U*-matrices.

If we form states |*a*_{i}〉_{A} on block *A* formed from sites 1 through to *i* and |*a*_{i}〉_{B} on block *B* formed from sites *i*+1 through to *L* as
2.10
and similarly for |*a*_{i}〉_{B}, then the left- and right-normalization conditions correspond to orthonormality conditions _{A}〈*a*_{i}|*a*′_{i}〉_{A}=*δ*_{ai,a′i} and _{B}〈*a*_{i}|*a*′_{i}〉_{B}=*δ*_{ai,a′i}, respectively. Generally, it is true that one can write
2.11
where . This establishes the connection between MPS notation and block growth , which—taking into account decimation of the states |*a*_{i+1}〉_{A}—makes the connection to RG growth steps.

In the mixed-canonical representation, the nature of the MPS as a systematic low-entanglement approximation to a quantum state becomes most obvious. If we do an SVD on *Ψ*=*USV*^{†} and absorb matrices *U* and *V*^{†} into the adjoining *A*- and *B*-matrices, this does not change the orthonormality of the block states because of the properties of *U* and *V*^{†}, but the state now takes the form
2.12
where the sum runs over the *r* non-vanishing (positive) diagonal elements *s*_{ai} of *S*, linking up states in *A* and *B* in pairwise fashion. This is the so-called Schmidt decomposition of a quantum state, with *r* being the Schmidt rank. If we have to cut the sum at some *D*<*r*, the quality of the approximation will depend on how rapidly the *s*_{ai} decay; the optimal approximation is given by keeping the *D* largest *s*_{ai}. But it is easy to see that at the same time, the reduced density operators for *A* and *B* are given as and analogously . A good approximation is then possible if the eigenspectrum of the reduced density operators decays quickly. It is usually not known, but we may argue that a quickly decaying spectrum corresponds to small entanglement , and MPSs systematically approximate low-entanglement states well. Let us take that for granted and postpone the crucial question when this is actually the case to later, and focus on operations on MPSs.

The most important operations on MPSs are overlaps and expectation values. Let us consider the calculation of 〈*ϕ*|*ψ*〉, the generalization to is completely obvious. We have
2.13
Transposing the scalar formed from the (which is the identity operation), we arrive at adjoints with reversed ordering,
2.14
To represent this graphically, we introduce an inverted matrix-blob to correspond to the complex-conjugate matrix *M*^{σi}* (figure 1). In the resulting pictorial representation of the overlap calculation (figure 2), this calculation then becomes much simpler, if we follow the rule that all bond indices are summed over. To evaluate this expression efficiently, the following order of contractions is advisable:
2.15
which contracts in *O*(*LD*^{3}) operations. In order to calculate an expectation value, say at site *i*, the sum to be evaluated there changes to
2.16
where *E* is the result of all the contractions up to site *i*−1.

If we exploit left-normalization, it is easy to see that for 〈*ψ*|*ψ*〉 matrix contractions collapse trivially to the identity *I* (similarly for right-normalization, if we rearrange the contractions to go from *L* to 1). In mixed-canonical notation, expectation values for operators that sit right next to the point where left- and right-normalized matrices meet become particularly efficient, as essentially all contractions become trivial.

The last important operation on an MPS is its *compression*. Assume that some operation (like the application of a matrix product operator (MPO) to an MPS, see below) generates an MPS of dimension *D*′, but we can only handle dimension *D* efficiently. How do we approximate |*ψ*〉 by , an MPS of dimension *D*, such that the distance between the two states becomes minimal, i.e. the approximation becomes optimal? Essentially two techniques exist, one based on SVD and not necessarily optimal, but frequently very close to it, and a variational technique, that is more complicated, but optimal.

Let us consider an MPS in mixed-canonical representation, with *Λ*^{[ℓ]}=diag(*s*_{1},*s*_{2},…),
2.17
from which we read off the Schmidt decomposition , where the states on *A* and *B* form orthonormal sets, respectively. We now look for the state that approximates |*ψ*〉 best in the 2-norm and can be spanned by a smaller number of *D* states each in *A* and *B*. This is achieved by retaining the *D* largest *s*_{aℓ}, and the compressed state simply reads . If normalization is desired, the retained singular values must be rescaled.

This procedure rests on the orthonormality of the states on *A* and *B*, therefore can only be carried out at one bond. In order to shrink the state on all sites, we have to work our way through all mixed-canonical representations, say from right to left, truncate and shift the boundary between left- and right-normalized matrices by one site to the left, using techniques from canonization.

After the first step of right-canonization of a left-canonical state, it reads
2.18
where I have already reshaped *B*, which is right-normalized and guarantees that states formed as are orthonormal. But so are the states
2.19
as SVD guarantees *U*^{†}*U*=1: we are just doing a basis transformation within the orthonormal basis set constructed from the left-normalized *A*^{σi}. Hence, we have a correct Schmidt decomposition as
2.20
The difference to a right canonization is now the truncation: matrices *U*, *S* and *B*^{σL} are truncated (and singular values possibly renormalized) to , and just as explained before: retain the *D* largest singular values. is still right-normalized. The next *A*^{σL−1} to the left, and are multiplied together to form *M*^{σL−1}. By reshaping and singular value decomposing,
2.21
we obtain right-normalized *B*^{σL−1}, truncate *U*, *S* and *B*^{σL−1} to , and , and the procedure continues. At the end, the compressed MPS is right-normalized and given by -matrices.

The disadvantage of the procedure is that there is a one-sided dependence of truncations on each other: truncations further to the left depend on those to the right, but not vice versa. If the truncation is small, the introduced additional inaccuracy is minor; however, for cases where large truncations may occur, the dependence might become too strong and the truncation far from optimal.

The optimal approach is to start from an ansatz MPS of the desired reduced dimension, and to minimize its distance to the MPS to be approximated iteratively, i.e. by changing its *M*^{σ} matrices iteratively. The matrices play the role of variational parameters.

Minimizing with respect to is a highly nonlinear optimization problem in the . But this can be done iteratively and linearly as follows. Starting with an initial guess for , we sweep through the set of site by site, keeping all other matrices fixed and choosing the new , such that the distance is minimized. Repeating this sweep through the matrices several times will lead to a converged optimal approximation.

The new is found by extremizing with respect to , which only shows up in . We find
The sum over ** σ*** runs over all physical sites except

*i*. Assuming that we keep at each iteration in mixed-canonical form, , the first term simplifies owing to the normalization conditions to . Then, we obtain the equation 2.22 where is given by 2.23 with

*L*and

*R*as indicated in figure 3. In fact, calculating

*L*and

*R*is nothing but carrying out the first steps of an overlap calculation, starting from the left or right. If one sweeps through the system from left to right and back, one can build

*L*and

*R*iteratively from previous steps, which is the most efficient way.

To make this work for an entire chain, we have to shift the boundary between the left- and right-normalized matrices as we move through the chain. But this can be done by applying one step of normalization (by SVD) to the matrices on a single site.

### (b) Matrix product operators

In the previous section, we have considered a local operator , which ties in nicely with the local MPS notation. For operators like the Hamiltonian which consists of sums of (in practice usually) local terms of operator products, the natural question arises whether such operators can also be expressed in a form resembling an MPS. This is indeed the case and leads directly to MPOs. The most general operator acting on our *L*-site system is given as
2.24
Interpreting the scalar coefficients as *c*_{(σ1,σ1′),…,(σL,σL′)}, we can, by analogy to an MPS, conclude that SVDs allow a decomposition as
2.25
where at each lattice site, we introduce a set of *d*^{2} matrices *W*^{σiσi′}. The previous statements about gauge degrees of freedom, orthonormalization types, etc. generalize directly to MPOs.

The graphical representation of an MPO is also an obvious extension from the MPS case, with an ingoing physical leg down and outgoing physical leg going up (figure 4).

The application of an MPO to an MPS is very elegant (figure 5), The elegance stems from the observation that the form of the MPS remains invariant, with an increase of the matrix size: the new MPS matrices, 2.26 have the multiplied dimension of the MPS and the MPO. At the same time, like in the case of an overlap, the operation is only of small polynomial, not exponential complexity. Of course, it remains to show that it is in fact possible to construct the MPO of, e.g. an Hamiltonian ‘by hand’, because the conceptual way of using an SVD might be exponentially complicated. But as we will see, this is no big problem indeed [13,14].

Let us consider the following simple Hamiltonian
2.27
As this is a shorthand for sums of tensor products of operators while omitting all the identity operators on other sites, it is convenient to reconsider the building block combined with its associated projector |*σ*〉〈*σ*′| to become an operator-valued matrix such that the MPO takes the simple form
2.28

In order to construct the MPO, we move through an arbitrary operator ‘string’ appearing in : starting from the right end, the string contains unit operators, until at one point, we encounter one of (in our example) three non-trivial operators. For the field operator, the string part further to the left may only contain unit operators; for the interaction operators, the complementary operator must follow immediately to complete the interaction term, to be continued by unit operators further to the left. For book-keeping, we introduce five corresponding states of the string at some given bond: state 1, only units to the right; states 2,3,4, one , , just to the right; state 5, completed interaction somewhere to the right. Comparing the state of a string left and right of one site, only a few transitions are allowed: by the unit operator , by , by , by . Furthermore by , by and by , to complete the interaction term, and for a completed interaction by the unit operator . Furthermore, all string states must start at 1 to the right of the last site and end at 5 (i.e. the dimension *D*_{W} of the MPO to be) to the left of the first site. This can now be encoded by the following operator-valued matrices:
2.29
and on the first and last sites
2.30
Inserting the explicit operator representations gives the desired *W*^{σσ′}-matrices for the MPO. The MPO in this example has dimension *D*_{W}=5. For longer-ranged Hamiltonians, further ‘intermediate states’ have to be introduced: for
2.31
the bulk operator would read
2.32
as one can verify by working out the operator string paths allowed by this construction. Except for special cases like exponentially decaying interactions, which find a very compact representation [15,16], the dimension of the MPO will grow with interaction range; however, we may always compress an MPO like an MPS to a smaller dimensional object with a minimal loss of information.

### (c) The finite-system density-matrix renormalization group algorithm

The so-called *finite-system DMRG algorithm* is nothing but a variational ground-state search in the ansatz space of *D*-dimensional MPSs; hence, it is also called vMPS. Let me mention that historically, a finite-system DMRG was not framed in the language of an MPS, but rather in terms of the analysis of reduced density operators of subsystems (blocks) of the chain. If, on the other hand, one assumes that MPSs are very efficient at encoding one-dimensional states (as they are, see §2*e*), then the following algorithm for ground-state searches in that state class follows naturally—and it then turns out that White [1] had had the perfect intuition, as the algorithms are identical.

As increasing *D* leads to ansatz spaces that are supersets of the previous ones, the quality of approximation improves monotonically; in view of the exact MPS representation of quantum states, the limit is exact. These two features allow for reliable extrapolation of results obtained for sequences of finite *D*.

Finding the ground state is equivalent to minimizing
2.33
or, using a Lagrange multplier *λ*, minimizing
2.34
Inserting the MPS encoding of |*ψ*〉, this turns into a highly nonlinear problem in the variables which is unsolvable in numerical practice. But in fact, this problem can be turned into a sequence of linear problems, whose solutions will lead to an iterative improvement of the solution (in the sense that energy is monotonically lowered), where, with suitable numerical tricks, it can be ensured (in the sense of actual numerical practice) that the (energetically) optimal state in this ansatz class can be reached. One starts from a guess state, which can be picked either randomly (i.e. *A*^{σ} matrices with random entries) or be obtained from some warm-up algorithm such as the so-called *infinite-system* DMRG algorithm [1,2]. We consider the MPS state provided by the initial guess and pick one site ℓ, while all other sites remain inert. We then optimize energy with respect to the *A*^{σℓ} matrices belonging to this site. Then, the variables appear in equation (2.34) only in a quadratic form, for which the determination of the extremum is a benign linear algebra problem. This will lower the energy, and find a variationally better state, but of course not the optimal one. Now one continues to vary the matrix elements on another site for finding a state again lower in energy, moving through all sites multiple times, until the energy does not improve anymore.

Let us analyse the parts of equation (2.34). If we keep the state suitably mixed normalized, the normalization conditions imply that 2.35

Considering , with in MPO language and taking the extremum of equation (2.34) with respect to , we find
2.36
with the mathematical objects as defined graphically in figure 6. This is in fact a very simple eigenvalue equation; if we introduce a matrix *H* by reshaping and a vector *v* with , we arrive at an eigenvalue problem of matrix dimension (*dD*^{2}×*dD*^{2}),
2.37
Solving for the lowest eigenvalue *λ*_{0} gives us a , which is reshaped back to , *λ*_{0} being the current ground-state energy estimate. In general, *dD*^{2} is too large for an exact diagonalization, but as we are only interested in the lowest eigenvalue and eigenstate, an iterative eigensolver that aims for the ends of the spectrum will do. Typical methods are the Lanczos or Jacobi–Davidson large sparse matrix solvers.

The optimal algorithm then runs as follows.

— Start from some initial guess for |

*ψ*〉, which is right-normalized, i.e. consists of*B*-matrices only.— Calculate the

*R*-expressions iteratively for all site positions*L*−1 through to 1 iteratively.—

*Right sweep*: starting from site ℓ=1 through to site*L*−1, sweep through the lattice to the right as follows: solve the standard eigenproblem by an iterative eigensolver for*M*^{σℓ}, taking its current value as the starting point. Once the solution is obtained, left-normalize*M*^{σℓ}into*A*^{σℓ}by an SVD to maintain the desired normalization structure. The remaining matrices of the SVD are multiplied to the*M*^{σℓ+1}to the right, which will be the starting guess for the eigensolver for the next site. Build iteratively the*L*expression by adding one more site. Move on by one site, , and repeat.—

*Left sweep*: starting from site ℓ=*L*through to site 2, sweep through the lattice to the left as follows: solve the standard eigenproblem by an iterative eigensolver for*M*^{σℓ}, taking its current value as the starting point. Once the solution is obtained, right-normalize*M*^{σℓ}into*B*^{σℓ}by an SVD to maintain the desired normalization structure. The remaining matrices of the SVD are multiplied to the*M*^{σℓ−1}to the left, which will be the starting guess for the eigensolver for the next site. Build iteratively the*R*expression by adding one more site. Move on by one site, , and repeat.— Repeat right and left sweeps, until convergence is achieved. Convergence is achieved if energy converges, but the best test is (using MPOs) to consider to see whether an eigenstate has been reached; this expression should approach 0 as closely as possible.

In this iterative process, the energy can only go down, as we continuously improve by varying the parameters. In fact, this most simple version of the DMRG can be improved in terms of stability and performance but the core algorithm remains unchanged; for details, see e.g. Schollwöck [5].

### (d) Variational or renormalization group method?

In the exposition just given, the DMRG clearly emerges as a variational method, and the link to true RG methods (with scaling and truncating) is quite tenuous. It can, however, be put on a more solid footing if one considers the reduced density operators of subsystems and takes those to the thermodynamic limit [6]. For these, a fixed point relationship holds. The DMRG, therefore, sets up an RG flow in the space of reduced density operators, not in the space of Hamiltonians, as other RG methods would do.

### (e) Why does it work and why does it fail?

The performance of all methods presented here rests on whether a quantum state can be approximated well by an MPS with bonds of a manageable dimension *D*. As we can cut a state into two parts across such a bond, we have to consider bipartite quantum systems *AB*. Any state , where the states |*i*〉_{A} and |*j*〉_{B} form orthonormal bases of dimensions *N*_{A} and *N*_{B}, respectively. Thinking of the *ψ*_{ij} as entries of a rectangular matrix *Ψ* (dimension *N*_{A}×*N*_{B}), the reduced density matrices *ρ*_{A} and *ρ*_{B} take the forms
2.38
If we assume |*ψ*〉 to be normalized, the eigenvalues *w*_{a} of *ρ*_{A} are non-negative, and sum to 1. This allows us to interpret them directly as statistical weights. If we assume that we know |*ψ*〉 exactly, the question of approximability reduces to how much statistical weight rests in the *D* eigenstates with the largest eigenvalues. In the few cases where this is the case, such analyses have been carried out in one and two dimensions starting with Peschel *et al.* [17] and Okunishi *et al.* [18]. They reveal that in one dimension for gapped systems, eigenvalues *w*_{a} generically decay exponentially fast (roughly as ), which explains the success of the DMRG, but in two-dimensional stripe geometries of size *L*×*W*, *L*≫*W*, *c*∝1/*W*, such that with increasing width *W* (increasing two dimensionality) the eigenspectrum decay is so slow as to make DMRG inefficient.

Usually, we have no clear idea about the eigenvalue spectrum; but it turns out that in such cases, entanglement entropies can serve as ‘proxy’ quantities, namely the von Neumann entanglement or entanglement entropy. It is given by the non-vanishing part of the eigenvalue spectrum of *ρ*_{A} as
2.39
If we now consider a bipartitioning *A|B* where *AB* is in the thermodynamic limit and *A* of size *L*^{𝒟}, with 𝒟 the spatial dimension, the so-called *area laws* [19–23] predict that for ground states of short-ranged Hamiltonians with a gap to excitations entanglement entropy are not extensive, but proportional to the surface, i.e. *S*(*A*|*B*)∼*L*^{𝒟−1}, as opposed to thermal entropy. This implies *S*∼*cst*. in one dimension and *S*∼*L* in two dimensions. At criticality, a much richer structure emerges, which usually involves the presence or absence of logarithmic corrections (, in one dimension and two dimensions), see Srednicki [20], Vidal *et al.* [24], Latorre *et al.* [25], Gioev & Klich [26] and Barthel *et al.* [27]. It should be emphasized that these properties of ground states are highly unusual: in the thermodynamic limit, a random state out of Hilbert space will indeed show extensive entanglement entropy with probability 1.

In a mathematically non-rigorous way, one can now make contact between the DMRG and the area laws of quantum entanglement: between two *D*-dimensional state spaces for *A* and *B*, as provided by a *D*-dimensional MPS, the maximal entanglement is in the case where all eigenvalues of *ρ*_{A} are identical and *D*^{−1} (such that *ρ*_{A} is maximally mixed); meaning that one needs a state of dimension 2^{S} and more to encode entanglement *S* properly. This implies that for gapped systems in one dimension, an increase in system size will not lead to a strong increase in *D*; in two dimensions, *D*∼2^{L}, such that the DMRG will fail, even for relatively small system sizes, as resources, i.e. MPS dimensions, have to grow exponentially (this however does not exclude very precise results for small two-dimensional clusters or quite large stripes). Critical systems in one dimension are borderline cases: the logarithmic correction makes them untreatable in the thermodynamic limit, but enough data for excellent finite-size extrapolations can be obtained.

## 3. Generalizations of the density-matrix renormalization group ansatz to higher dimensions

### (a) Tensor network states

Extensive attempts to simulate two-dimensional quantum systems within the conventional paradigm of DMRG as a one-dimensional method have been made over the last 15 years; generically, the two-dimensional system is seen as a one-dimensional snake where interactions that on the true lattice are short-ranged are now fictitiously long-ranged (with a range corresponding to the smaller of the two linear sizes). As discussed above, the entanglement scaling implies a matrix dimension that grows exponentially with system size; nevertheless, impressive results on wide strips have been obtained ([28,29] among the first and most recent ones), but it remains clear that except for astute finite-size extrapolations, there is no way towards the thermodynamic limit.

On the other hand, we may interpret the two matrix indices of *M*^{σi} as bond indices, with the row and column corresponding to bonds going to the left and right neighbour on a one-dimensional chain, as already implied by the graphical representation. It is, therefore, a natural extension to consider *tensors* with four indices that correspond to the four bonds up, down, left and right on a two-dimensional square lattice and similar generalizations on other lattices (triangular, hexagonal and so on). Instead of an MPS, we would consider [10] a tensor product state (TPS) or tensor network states (TNSs)—a popular other name being projected entangled pair states (PEPS; [11])—reading
3.1
where the stands for a contraction over all connected legs in a tensor network; in the one-dimensional case this reduces immediately to the matrix multiplication of an MPS. For a graphical representation, see figure 7.

The attractive feature of this ansatz is the following: imagine a vertical cut through a system of size *L*×*L*. In a PEPS, this will cut *L* legs of dimension *D* each. This allows for *D*^{L} configurations. The maximum amount of entanglement across the cut that can be encoded is given by the logarithm as ; this means that it is proportional to the surface and hence the area law for gapped two-dimensional systems is immediately obeyed, indicating that this is a valid ansatz class in two dimensions.

It should be remarked that the concept of a TNS is much more general than just a set of tensors reflecting the underlying lattice structure. The TNS can also encode hierarchy of scales and make use of disentanglers to keep descriptions as efficient as possible, as is done by the multi-scale entanglement renormalization ansatz (MERA) scheme [30].

### (b) Tensor network algorithms: finding a projected entangled pair states ground state

At first sight, the solution of two-dimensional quantum ground-state problems seems immediately at hand. We may also generalize from MPOs to tensor product operators (TPOs), express a Hamiltonian as a TPO efficiently, and solve the variational ground-state equation iteratively as in one dimension, varying one tensor after another while keeping all others fixed.

The problems that occur are twofold: (i) the contraction of the tensor network can only be approximate, never exact as in the one-dimensional case and (ii) numerical stability is an issue as orthonormalization properties are lost. The first is immediately obvious if one tries to contract two tensors versus two matrices: the product of two matrices is again a matrix, ; however, the contraction of two 4-leg tensors over one leg produces one 6-leg tensor. This can be turned into an exact mathematical statement that the contraction of a two-dimensional tensor network is indeed an non-deterministic polynomial time hard problem, i.e. impossible on a classical computer. Approximation schemes, therefore, have to be involved. One example would be to read a PEPS as a succession of rows of tensors applied to an MPS at the bottom end, to be completed by an MPS at the top end. The tensor rows can be interpreted as MPOs and we would then have a succession of MPO–MPS contractions; these lead to a multiplicative growth of bond dimensions, which is checked by an MPS compression algorithm. But this is obviously only approximate. The second problem is less fundamental, but as annoying. In our exposition of one-dimensional DMRG, we have made extensive use of the left- and right-normalizability of MPSs. But this was only possible because we considered open boundary conditions; for a ring (periodic boundary conditions), there is no separate left- and right-normalization. One immediate consequence is that the eigenproblem of the DMRG turns into a generalized eigenproblem, *Hv*=*λNv*, whose numerical stability depends on the condition number of *N*, which can turn out to be large. In a two-dimensional system, even with open boundary conditions, any tensor singled out for variational optimization generates a ring-like structure of tensors, implying potential numerical instability.

Let me conclude this outlook on two dimensions by a general statement on TNSs. As all approximative methods, they have some built-in bias. In the case of MPSs and TNSs, it is their preference for *low-entanglement* states. If the accuracy is not very high, and the energy resolution therefore not very good, among very close-lying low-energy states those with less entanglement will be preferred, which might bias to the wrong kind of ground state, e.g. a fixed singlet pattern over spin-liquid-like fluctuating ground states. In one-dimensional MPS (DMRG) applications, *D* can nowadays be pushed to sizes of 10 000 or so, resolving energies down to machine precision, and this bias does not seem to create serious problems. In two-dimensional applications, where current dimensions are very small, this issue is potentially much more severe!

## 4. Outlook

Another important extension of the DMRG I have completely omitted is the extension to the time domain, i.e. the simulation of quantum dynamics close and far from equilibrium. There are a variety of methods, all based on Vidal’s insight that time evolutions of MPSs can be done very efficiently [31–35], which have given rise to numerous successful applications. Here, the current effort is mainly to try to extend the range of such simulations. In ground-state calculations, in one dimension, essentially all algorithmic questions of interest have been settled satisfactorily; in two dimensions, a sober assessment would indicate that at the moment, the best PEPS results do not yet significantly exceed the results of other methods, and the road to higher accuracies and higher stability seems very steep, but probably not hopeless.

## Footnotes

One contribution of 11 to a Theme Issue ‘New applications of the renormalization group in nuclear, particle and condensed matter physics’.

- This journal is © 2011 The Royal Society