# Density-matrix renormalization using three classes of block states

## I Introduction

The density matrix renormalization group (DMRG) method introduced by S. R. White in 1992  has rapidly become one of the leading numerical tools for the study of one-dimensional (1D) and quasi-1D strongly correlated systems. Accuracy and flexibility are two main qualities of this approach. Indeed, the DMRG method allows to study large many-body systems containing hundreds of sites and to determine several ground-state and low-lying excited-state properties with a precision nearly comparable to that of conventional exact diagonalization methods. The success of this technique is readily demonstrated by its wide range of applications. These concern Heisenberg spin Hamiltonians , the $`t`$-$`J`$ model , Hubbard-like Hamiltonians on various lattices (e.g., dimerized polymer chains and Bethe lattices) , the Kondo impurity problem , acoustic phonons , etc.

In its simplest form, the method can be regarded as an iterative projection technique, which allows to include the most relevant part of the ground-state wave-function on a limited number of many-body states and which leads to renormalized block interactions between different regions in real space. The accuracy of the calculation is mainly controlled by the number $`m`$ of basis states which are retained in each block after a renormalization iteration. When combined with efficient diagonalization techniques for sparse matrices (e.g., Lanczos or Davidson methods) the ground-state energy of an infinite 1D system can be extrapolated with a remarkable precision at a much lower computational cost than conventional exact diagonalizations (e.g., straightforward Lanczos). However, as the systems under study become more complex (e.g., itinerant electrons in rings, ladders or multi-band wires) and as more delicate properties are considered (e.g, charge and spin gaps, dynamical response functions , etc.) the DMRG method faces its own limitations, particularly since accuracy improves slowly with increasing $`m`$. In such cases the value of $`m`$ and thus the size of the projected Hilbert space required for an appropriate accuracy become excessively large. The involved computer resources exceed reasonable limits and in many relevant situations the calculations are simply non-feasible. It is therefore of considerable interest to improve the efficiency of the DMRG method in order to reduce the computational effort without significant loss of accuracy.

In this paper an extension of DMRG method is presented which takes better profit from the hierarchy among the states spanning the renormalized blocks. As proposed by White in his original paper  each block $`b`$ is best described in terms of the eigenstates $`|i_b`$ of the reduced ground-state density matrix. The corresponding eigenvalues $`\lambda _i`$ represent the population of the block-state $`|i_b`$ and thus provide a natural hierarchy among them. Two classes of states are considered in the original DMRG procedure: the $`m`$ states having the largest $`\lambda _i`$ which are retained and the remaining states which are completely neglected. In this paper a larger flexibility is introduced by distinguishing an additional group of $`p`$ states ($`\lambda _1\mathrm{}\lambda _m\lambda _{m+1}\mathrm{}\lambda _{m+p}`$) which are primarily important for the description of the ground-state wave function at block $`b`$ and which can be neglected when they participate to long-range interblock fluctuations. It is shown that in this way the computer resources involved in the calculations can be appreciably reduced without introducing significant additional inaccuracies. Moreover, it is also interesting to vary the values of $`m`$ and $`p`$ and to quantify how this reflects in the physical properties in order to gain further insight into the internal many-body structure of the renormalized superblocks.

The rest of the paper is organized as follows. In the next section the 3-classes renormalization method is presented and compared with White’s original approach . In Sec. III the method is applied to a few relevant models for which exact or good numerical results are available. Improvements and limitations of the 3-classes method are discussed in particular by analyzing the results as a function of the sizes $`m`$ and $`p`$ of class-1 and class-2 subspaces. A model independent scheme for optimizing the choice of $`m`$ and $`p`$ and for extrapolating to the limit where $`m`$ and $`p`$ are infinite is also considered. Finally, Sec. IV summarizes our conclusions and outlines a self-consistent extension of the present approach.

## II 3 classes DMRG method

In this section we present a density-matrix renormalization method using 3 classes of block states. First of all, White’s original DMRG algorithm for infinite one-dimensional systems  is briefly recalled for the sake of clarity and in order to facilitate comparisons.

The DMRG method is an iterative procedure. In the case of infinite open chains the ground-state properties of the infinite system are obtained from the extrapolation of a succession of calculations on finite-length chains. At each renormalization group (RG) iteration $`N`$, the number of sites $`L`$ is increased by 2 ($`L=2N+2`$). In spite of the fact that $`L`$ increases with $`N`$, the dimension of the Hilbert space is kept constant by means of the following approximation. At each RG iteration the Hilbert space of the $`L`$-site chain is projected onto a subspace $``$ spanned by a limited number of many-body states $`\mathrm{\Psi }_{ijkl}`$ which are intended to take into account the most relevant part of the ground-state wave function. The $`\mathrm{\Psi }_{ijkl}`$ are constructed as the antisymmetrized direct product of many-body states, which are usually referred to as “block states.” These blocks correspond to different regions in real space as illustrated in Fig. 1. Let $`\psi _i^L`$ ($`\psi _i^R`$) stand for $`i`$th many-body state of the left (right) block ($`1im`$) and let $`\varphi _j`$ be the $`j`$th state of a central (unrenormalized) site ($`1jm_1`$). Then

$$\mathrm{\Psi }_{ijkl}=\psi _i^L\varphi _j\varphi _k\psi _l^R.$$
(1)   
White showed that the error introduced by the truncation of the Hilbert space is minimized when the states spanning the subspace of each renormalized block $`b`$ are the eigenvectors corresponding to the $`m`$ largest eigenvalues of the ground-state density matrix reduced to $`b`$. In other words, only the block states having the largest occupations are kept.

The DMRG algorithm can be summarized by four main steps involved in going from iteration $`N`$ to iteration $`N+1`$: i) The ground state of the $`L`$-site system ($`L=2N+2`$), which corresponds to the Hilbert space $`(N)`$ of the iteration $`N`$, is determined. This involves a calculation in a space of limited dimension $`(mm_1)^2`$. Sparse-matrix diagonalization methods (e.g., Lanczos or Davidson algorithm) are normally used. ii) The reduced ground-state density matrices $`\rho ^L`$ and $`\rho ^R`$ of the superblocks are calculated. As illustrated in Fig. 1, the superblocks are formed by one of the side blocks ($`L`$ or $`R`$) and its neighboring site. The eigenvalues $`\lambda _i^L`$ ($`\lambda _i^R`$) and eigenvectors $`\psi _i^L`$ ($`\psi _i^R`$) of $`\rho ^L`$ ($`\rho ^R`$) are determined. iii) The superblock subspace of dimension $`mm_1`$ is projected onto the $`m`$ most populated states derived in the previous step ($`\lambda _1\mathrm{}\lambda _m`$). These are the left and right block-states of the next iteration. In this way the central sites have been projected or renormalized into the side blocks. iv) The Hilbert space $`(N+1)`$ of the iteration $`N+1`$, which corresponds to a system with $`2N+4`$ sites, is constructed according to Eq. (1). The renormalized Hamiltonian is determined by performing the appropriate unitary transformation and projection. Finally, the procedure loops back to step i) until the physical properties extrapolated to $`L=\mathrm{}`$ converge.

Two groups or classes of block states are distinguished in White’s DMRG procedure. The first one is given by the $`m`$ most populated states which are retained, and the second one by the remaining $`mm_1m`$ least populated states which are projected out. This criterion certainly respects the hierarchy of states given by the eigenvalues or populations of the different block states. However, it results in a rather rigid procedure, since the contributions of a block state must be either kept in full or completely neglected. Except in cases where the number of states that can be kept is large enough to ensure the desired accuracy or when a systematic extrapolation to $`m=\mathrm{}`$ can be performed, the choice of $`m`$ is a delicate and drastic one. In general there is no clear gap in the populations of the block states which could justify where the truncation should take place. It is our purpose to render the method more flexible and efficient by the introduction of a new intermediate class of block states which are intended to contribute in a restricted way to the construction of the Hilbert space of the total system. Physically, one expects that the description of the internal degrees of freedom of the superblocks and the renormalization of block-site interactions should be more demanding and more relevant to the ground-state properties than the description of fluctuations between left and right block-states which are farther apart. The intermediate class of states should take this into account. Moreover, with three classes of block states, the transition between kept and discarded many-body processes should be smoother and the actual choice of $`m`$ less critical.

We consider 3 groups or classes of eigenstates of the block density matrix: i) The first group (class 1) is given by the $`m`$ most populated states $`\psi _i`$ where $`1im`$ and $`\lambda _1\mathrm{}\lambda _m`$. ii) The second group (class 2) contains the eigenstates $`\psi _i`$ associated to the following $`p`$ eigenvalues, i.e., $`m+1im+p`$ with $`\lambda _m\lambda _{m+1}\mathrm{}\lambda _{m+p}`$. iii) Finally, the third group (class 3) refers to the remaining $`mm_1(m+p)`$ states. The states of the first and third classes are treated just as in White’s algorithm, i.e., all direct-product states or configurations of the complete system $`\mathrm{\Psi }_{ijkl}`$ which involve only class-1 block-states ($`i,lm`$) are taken into account while configurations involving a class-3 state ($`i`$ or $`l>m+p`$) are all projected out \[see Eq. (1)\]. An intermediate criterion is used for the second class, namely, class-2 states of the left superblock are included only in combination with class-1 states of the right superblock and vice versa. In other words, direct-product wave-functions $`\mathrm{\Psi }_{ijkl}`$ with $`m+1im+p`$ and $`m+1lm+p`$ are neglected. This is illustrated in Fig. 2.

In analogy with configurations interaction methods, the subspace spanned by the set of $`\mathrm{\Psi }_{ijkl}`$ with $`i,lm`$ can be regarded as a complete active space (CAS), the subspace spanned by the $`\mathrm{\Psi }_{ijkl}`$ with $`im`$ and $`m+1lm+p`$ or $`lm`$ and $`m+1im+p`$ defines then a kind of mono-excited space, and the discarded states obtained from $`\mathrm{\Psi }_{ijkl}`$ with $`m+1i,lm+p`$ can be viewed as double many-body excitations. Class-2 states have a relatively small weight in the ground-state wave-function. One therefore expects that with an appropriate choice of $`m`$ and $`p`$ the contributions of $`\mathrm{\Psi }_{ijkl}`$ obtained by the tensor product of two such block states should be particularly small, eventually even less important than the states involving products of a class-1 and class-3 states which are already discarded in White’s DMRG method. Therefore, the accuracy of a 3-classes renormalization should not differ significantly from that of a standard 2-class DMRG-calculations retaining $`m+p`$ states.

Depending on the values of $`m`$ and $`p`$, the 3-classes procedure may allow an appreciable reduction of the dimension of the many-body Hilbert space, namely, from $`D_0=m_1^2(m+p)^2`$ to $`D=m_1^2(m^2+2mp)`$. This is of considerable interest in many practical applications which are often limited by computer-memory needs and most critically by the computer time required for the determination of the ground-state wave-function at each RG iteration. The 3-classes approach can also be applied straightforwardly to more complex lattices such as Bethe lattices, which involve renormalizations at more than two blocks . In these cases the reduction of computational effort is even more important, since $`D_0/D=(m+p)^{n_b}/(m^{n_b}+n_bpm^{n_b1})`$ where $`n_b`$ is the number of renormalized blocks (e.g., $`n_b=z`$ for a Bethe lattice with $`z`$ nearest neighbors).

The eigenvalues $`\lambda _i`$ of the density matrix reduced to block $`b`$ represent the probability of block $`b`$ being in state $`\psi _i^b`$ ($`b=L`$ or $`R`$). Therefore, the traces $`\mathrm{\Sigma }_1=_{i=1}^m\lambda _i`$, $`\mathrm{\Sigma }_2=_{i=m+1}^{m+p}\lambda _i`$ and $`\mathrm{\Sigma }_3=_{i=m+p+1}^{m_1m}\lambda _i`$, measure the weight of each of the 3 classes of block states on the ground-state wave function ($`\mathrm{\Sigma }_1+\mathrm{\Sigma }_2+\mathrm{\Sigma }_3=1`$). Approximating for simplicity the norm of the retained part of the system wave-function by the product of the traces $`\mathrm{\Sigma }_i`$ on the different blocks and keeping only the leading terms in $`\mathrm{\Sigma }_2`$ and $`\mathrm{\Sigma }_3`$ ($`\mathrm{\Sigma }_2,\mathrm{\Sigma }_31`$), one estimates the truncation error $`ϵ(m,p)`$ involved in a 3-classes renormalization as $`ϵ(m,p)=n_b\mathrm{\Sigma }_3+(\genfrac{}{}{0pt}{}{n_b}{2})\mathrm{\Sigma }_2^2`$, where $`n_b`$ refers to the number of renormalized blocks ($`n_b=2`$ in the present calculations). If $`p=0`$, $`\mathrm{\Sigma }_2=0`$ and we recover White’s estimation of the truncation error $`1P_m`$ which is proportional to $`\mathrm{\Sigma }_3`$  ($`\mathrm{\Sigma }_3=1\mathrm{\Sigma }_1=1P_m`$ if $`\mathrm{\Sigma }_2=0`$). As discussed in Sec. III, comparison with exact results show that $`ϵ(m,p)`$ gives very good estimation of the errors involved in 3-classes calculations. In practice, $`ϵ(m,p)`$ can be used either to control the accuracy of the results, which often depend on the model Hamiltonian and its parameters, or to extrapolate to the limit $`ϵ(m,p)0`$ ($`m,p\mathrm{}`$) in order to further reduce the errors. A similar procedure has been already applied in the standard DMRG . Moreover, $`ϵ(m,p)`$ provides a practical criterion for optimizing the value of $`m`$ and $`p`$, for instance, when the size of the Hilbert space $`D(m,p)`$ is fixed. In general it is a good choice to keep $`\mathrm{\Sigma }_2^2`$ and $`2\mathrm{\Sigma }_3`$ of the same order of magnitude. In fact, for a fixed $`D(m,p)`$ a reduction of $`\mathrm{\Sigma }_2^2`$ (e.g, by decreasing $`p`$) implies an increase of $`\mathrm{\Sigma }_3`$ and vice versa.

## III Results

The density-matrix renormalization using 3-classes of states has been applied to several benchmark problems in order to check the accuracy of the method and to analyze the quality of the results as a function of the subspace dimensions $`m`$ and $`p`$. In particular comparison is made with the standard DMRG algorithm . The considered systems include the spin-$`1/2`$ Heisenberg chain, the periodic Hubbard chain as a function of correlation strength $`U/t`$, and the dimerization of polyacetylene using a distance-dependent tight-binding Hamiltonian. In the following our main results are discussed.

### A Spin-$`1/2`$ Heisenberg chain

In this section we consider the infinite spin-$`1/2`$ Heisenberg chain as a representative application of the 3-classes DMRG method to spin systems. This problem is particularly interesting from a methodological point of view since the exact Bethe-Ansatz solution as well as extremely good standard DMRG calculations are available . In Table I our results for the ground-state energy $`E`$ are compared with the exact value $`E_{ex}=J\mathrm{ln}2+J/4`$ and with DMRG calculations using White’s method which corresponds to $`p=0`$. In all cases, $`N=100`$ renormalization iterations are performed. This ensures that the extrapolated ground-state energy of the infinite-length chain is converged at least within $`10^6J`$. Small values of $`m`$ and $`p`$ are considered in Table I in order to compare the convergence and performance of the different methods. Notice, however, that more accurate calculations involving a larger number of states (e.g., $`m=m+p=44`$ ) are not at all demanding with present computer facilities. For $`m`$ and $`p`$ larger than the values shown in Table I, the 2-classes and 3-classes DMRG results converge very well to the exact ground-state energy and are almost indistinguishable from each other.

The 3-classes results recover accurately the exact ground-state energy even if as few as $`5`$ states are kept in the class-1 subspace or complete active space (CAS). For $`m=5`$ and $`m+p=21`$, $`EE_{ex}10^4J`$, which seems an acceptable accuracy taken into account that the size of the Hilbert space is reduced significantly ($`D_0/D=2.4`$). Our results show that the intrablock fluctuations between the $`m`$ class-1 states and the $`p`$ class-2 states play a major role in the description of the ground-state wave function. In fact, neglecting them yields so poor results that for $`m=m+p=5`$ the DMRG procedure diverges after a few iterations. In contrast, the effects of quantum fluctuations between class-2 states of different blocks are much less important for the determination for the ground-state energy. These contributions are neglected in a 3-classes DMRG calculation with $`m=5`$ ($`m=7`$) and $`m+p=21`$ which yields $`EE_{ex}=1.07\times 10^4J`$ ($`EE_{ex}=4.47\times 10^5J`$) and are included in the calculation with $`m=m+p=21`$ which yields $`EE_{ex}=2.55\times 10^5J`$. Similar conclusions are derived for other values of $`m`$ and $`p`$.

The calculations for different $`m`$ and fixed $`m+p`$ provide further insight on the convergence properties of the method. For very small $`m`$ the difference between 3-classes and 2-classes results decreases very rapidly for increasing $`m`$ (see Table I for $`m=5`$ and $`7`$). This indicates that a minimum size of the class-1 subspace is indispensable for obtaining good accuracy. Further increase of $`m`$ yields almost the same accuracy as the calculation with $`m=m+p=21`$. Nevertheless, the reduction of the dimension of the Hilbert space is still appreciable (e.g., $`D_0/D=1.4`$ for $`m=10`$). Let us recall that the error in $`E`$ decreases slowly with increasing $`D`$ already in the standard DMRG method . As in the case of very small $`m`$, the intrablock transitions involving class-2 states are very important for improving the accuracy ($`EE_{ex}=2.06\times 10^4J`$ for $`m=m+p=10`$). We conclude that the 3-classes DMRG procedure improves efficiency without significant loss of accuracy, provided that the size $`m`$ of the CAS is not too small ($`m7`$ in the present case).

The truncation error $`ϵ(m,p)=\mathrm{\Sigma }_2^2+2\mathrm{\Sigma }_3`$ provides a very good estimate of the accuracy of the calculations for all considered values of $`m`$ and $`p`$. In fact, $`EE_{ex}`$ follows very closely a linear relation with $`ϵ(m,p)`$. It is therefore simple to extrapolate to the limit $`ϵ(m,p)0`$ in order to further improve the results (see Table I). A similar procedure is used in standard DMRG . Moreover, the results suggest that the dimensions of class-1 and class-2 subspaces could be chosen in order to optimize accuracy by minimizing $`ϵ(m,p)`$ for a given computational cost.

### B Hubbard chain

The DMRG method converges very rapidly and accurately when applied to problems involving strong electron correlations such as the Heisenberg model. However, the performance of the method is very sensitive to the strength of the Coulomb interactions and to the degree of electron delocalization. For instance, in the case of the Hubbard chain at half-band filling and for a given number of block states $`m`$, the relative error $`\mathrm{\Delta }\epsilon `$ in the ground-state energy increases by an order of magnitude as we move from the strongly correlated to the uncorrelated limit. For example for $`m=80`$ ($`p=0`$), $`\mathrm{\Delta }\epsilon =8.1\times 10^5`$ for $`U/t=76`$, $`\mathrm{\Delta }\epsilon =4.9\times 10^5`$ for $`U/t=4`$, $`\mathrm{\Delta }\epsilon =3.1\times 10^4`$ for $`U/t=1`$ and $`\mathrm{\Delta }\epsilon =3.3\times 10^4`$ for $`U/t=0`$. Moderate and weak electron interactions are thus much more difficult to describe with the DMRG method. In addition, for small $`U/t`$ the algorithm requires a quite larger number $`N`$ of RG steps to converge. Typically, $`N`$ is about $`2`$$`3`$ times larger for $`U=0`$ than for very large $`U/t`$ . It is therefore considerably interesting to analyze the accuracy and convergence properties of the 3-classes DMRG method in applications to itinerant-electron models such as the Hubbard model.

In Fig. 3 results are given for the ground-state energy of the half-filled one-dimensional Hubbard chain as function of the Coulomb repulsion strength $`U/t`$. These were obtained by using the 3-classes DMRG method with various representative values of $`m`$ and $`p`$. The standard DMRG algorithm corresponds to $`p=0`$. In all cases, the number of renormalization iterations $`N`$ is increased until the infinite-length extrapolated energy has converged within an error smaller than $`10^6t`$. Starting from the uncorrelated limit, the accuracy improves considerably with increasing $`U/t`$. Very precise results are obtained particularly for $`U/t>4`$. However, for very large $`U/t`$ ($`U/t>16`$) the relative error may increase slightly. For example, for $`U/t=76`$ ($`U/t=16`$) the relative errors $`\mathrm{\Delta }\epsilon =(EE_{ex})/E_{ex}`$ are $`\mathrm{\Delta }\epsilon =4.2\times 10^4`$ ($`\mathrm{\Delta }\epsilon =1.1\times 10^4`$) for $`m=m+p=35`$, $`\mathrm{\Delta }\epsilon =8.0\times 10^5`$ ($`\mathrm{\Delta }\epsilon =3.5\times 10^5`$) for $`m=35`$ and $`m+p=80`$, and $`\mathrm{\Delta }\epsilon =8.0\times 10^5`$ ($`\mathrm{\Delta }\epsilon =3.0\times 10^5`$) for $`m=m+p=80`$. Notice that, although the differences between the calculation with $`m=m+p=35`$ and with $`m=m+p=80`$ are significant, the 3-classes results follow very closely the $`m=m+p=80`$ curve. This implies that the largest part of the improvements made by increasing $`m`$ from $`35`$ to $`80`$ are recovered by taking $`m=35`$ and $`m+p=80`$, i.e., by including the intrablock fluctuations. Transitions involving class-2 states on two different superblocks are not very important in this case. This allows to reduce the dimension of the Hilbert space by about 30%. The results seem remarkable, particularly if we recall that the 3-classes method is a controlled approximation in the sense that it respects the variational principle. At each RG iteration the energy obtained with $`m=35`$ and $`m+p=80`$ is an upper bound for the energy of a standard DMRG calculation with $`m=80`$ and $`p=0`$.

The same trends hold for weaker interactions. For $`U<4t`$, one observes that the differences between the relative errors in the 3-classes and in the standard DMRG calculations increase. At the same time, however, the $`m=m+p=80`$ curve starts to deviate more significantly from the exact result showing the that block states $`i`$ with smaller eigenvalues of the density matrix (i.e., weaker occupations $`\lambda _i<\lambda _{80}`$) should be taken into account. Quantitatively, neglecting the intrablock fluctuations involving states with $`\lambda _i<\lambda _{80}`$ is a much more important source of error than the interblock fluctuations neglected in the 3-classes method \[see Fig. 3(a)\]. The additional error introduced by the use of the 3-classes procedure ($`m=35`$ and $`m+p=80`$) is only 9% of the total error ($`m=m+p=80`$). Therefore, in order to further improve accuracy for $`U/t4`$ it should be more efficient to increase the size $`p`$ of the class-2 subspace instead of increasing $`m`$ with $`p=0`$ as in the standard DMRG method. This is consistent with our estimation of truncation error which is quadratic in the sum of the class-2 eigenvalues $`\mathrm{\Sigma }_2`$ and linear in $`\mathrm{\Sigma }_3`$.

In order to investigate this question we compare in Fig. 3(b) the accuracy of the 3-classes approach as a function of the size of class-1 and class-2 subspaces. Several values of $`m`$ and $`p`$ are considered which all correspond to approximately the same dimension of the Hilbert space of the complete system: $`D(m,m+p)=16[(m+p)^2p^2]D(80,80)=102400`$. Therefore, all these calculations demand essentially the same computer resources. One observes that there is an optimum choice for $`m`$ and $`p`$ for which the relative error $`(EE_{ex})/E_{ex}`$ is reduced by about 20–30% with respect to the standard calculation with $`m=m+p=80`$. As in the Heisenberg model, a minimum number of class-1 states is indispensable for describing the ground-state of the entire system with a reasonable accuracy. If $`m`$ is too small, for instance $`m=15`$, the results are poorer even if $`p`$ is very large. Once the most important configurations are taken into account as class-1 states, it is more useful to maximize $`p`$ in order to increase the flexibility of the wave-function within each superblock rather than to increase the number of class-1 states. This improves the treatment of intrablock correlations and allows a more precise renormalization of the interactions as the number of sites grows. The results shown in Fig. 3(b) for $`m=15`$ ($`p=221`$), $`m=20`$ ($`p=170`$) and $`m=35`$ ($`p=109`$) illustrate this behavior.

Taking into account the strong dependence of $`EE_{ex}`$ on $`(m,p)`$ and on $`U/t`$ (even for $`p=0`$) it is of considerable practical interest to have a model-independent criterion for choosing $`m`$ and $`p`$ and for controlling the accuracy of 3-classes DMRG calculations. Comparison between Figs. 3(b) and 3(c) shows that the estimation of the truncation error $`ϵ(m,p)=\mathrm{\Sigma }_2^2+2\mathrm{\Sigma }_3`$ follows the trends in $`EE_{ex}`$. Therefore, the minimization of $`ϵ(m,p)`$ provides an appropriate means of choosing $`m`$ and $`p`$ for a given dimension of the Hilbert space. Alternatively, one may use $`m`$ and $`p`$ such that $`ϵ(m,p)`$ is constant in order to ensure approximately the same accuracy in different calculations (e.g., as a function of $`U/t`$ or of the dimerization $`\delta `$). Notice, that $`ϵ(m,p)`$ can be easily computed at each RG iteration.

It is should be also noted that the minimum $`m`$, which yields better accuracy than the standard DMRG method, as well as the optimal choice for $`m`$ and $`p`$ for a given $`D`$ depend on the Coulomb interaction strength $`U/t`$. For example, $`m=20`$ and $`p=170`$ yields the best results for $`U/t1`$, but for $`U/t2`$ this is somewhat less accurate than the standard DMRG calculation. In the range $`0U/t4`$ the optimum $`m`$ ($`p`$) increases (decreases) as $`U/t`$ increases. This is the result of two contributions to the minimization of $`\mathrm{\Delta }\epsilon `$, which compete when the dimension of the Hilbert space is fixed. Increasing $`p`$ improves the description of the wave function within each block as well as the superblock renormalizations since it reduces $`\mathrm{\Sigma }_3`$. However, this is done at the expense of decreasing $`m`$ (or increasing $`\mathrm{\Sigma }_2`$) and therefore tends to worsen the description of the interactions between different blocks. For $`U/t>4`$ the results for $`m=35`$$`80`$ ($`p=109`$$`80`$) are very similar. Close to the optimum $`(m,p)`$ the energy differences remain small. Let us finally recall that the overall convergence and accuracy of DMRG methods improve as local charge fluctuations are reduced. The quantitative differences between different renormalization or projection strategies become therefore less important at large $`U/t`$.

### C Dimerized chains

The purpose of this section is to discuss the application of the 3-classes DMRG method to dimerized polymer chains. The dimerization of polyacetylene is particularly interesting since it depends sensitively on the details of the wave function and thus on the number of states $`m`$ used in DMRG calculations . A distance dependent tight-binding Hamiltonian is considered for the $`\pi `$ valence electrons in polyacetylene. Repulsive interactions between $`\sigma `$ electrons are modeled by a pairwise potential $`E_\sigma (r_{ij})`$ . The hopping integrals $`t(r_{ij})`$ and the repulsive potential $`E_\sigma (r_{ij})`$ are obtained from ab-initio calculations on the ground-state energy and first singlet-triplet excitation energy of the ethylene molecule . The main interest of the uncorrelated case is methodological, since this is the most difficult limit to study using the DMRG method. DMRG calculations for realistic finite values of the Coulomb repulsion $`U`$ are reported in Ref. .

In Fig. 4 results are given for the ground-state energy $`E(\delta )`$ of an infinite polyacetylene chain vs. the dimerization $`\delta =|r_{i,i+1}r_{i1,i}|/2`$, as obtained using the model Hamiltonian described above. The exact tight-binding result is also given for the sake of comparison. As in the undimerized Hubbard chain, the 3-classes calculations with $`m=35`$ and $`m+p=150`$ remove the largest part of the discrepancies between the $`m=m+p=35`$ and $`m=m+p=150`$ standard calculations. In particular the position of the minimum $`\delta _{min}`$ is improved considerably. However, notice that the error $`\mathrm{\Delta }E=EE_{ex}`$ decreases with increasing $`\delta `$ and therefore, $`\delta _{min}`$ remains overestimated. For $`m=m+p=35`$ ($`m=m+p=150`$) one obtains $`\delta _{min}=0.389a_0`$ ($`\delta _{min}=0.306a_0`$). Using $`m=35`$ and $`m+p=150`$ one finds $`\delta _{min}=0.310a_0`$, while the exact result is $`\delta _{min}=0.303a_0`$. As expected, for a given $`m`$ and $`p`$, the 3-classes method converges better towards the more demanding 2-classes calculation ($`m+p`$ states per block) as the dimerization increases.

## IV Summary and Outlook

In this paper, an extension of the the density matrix renormalization group (DMRG) method has been presented, which achieves larger flexibility and improved performances without increasing the computational effort. An intermediate class of block states is introduced between the states which are fully taken account and those which are discarded at each renormalization step. These states contribute only as far as intrablock fluctuations are concerned and provide a smoother criterion for discerning between retained and discarded many-body states. Applications of the method to one-dimensional models (Heisenberg, Hubbard and dimerized tight-binding) show that in this way the involved computer resources can be reduced without significant loss of accuracy. Varying the number $`p`$ of intermediate states relative to the number $`m`$ of block states which are fully retained one obtains further insight into the DMRG method, particularly concerning the internal structure of the renormalized blocks and the relative importance of intrablock and interblock fluctuations. For example, for a given dimension $`D(m,p)`$ of the Hilbert space of the renormalized system, optimum values of $`m`$ and $`p`$ can be determined which minimize the error in the ground-state energy. In these cases the discrepancies with exact results are 20–30% smaller than in standard DMRG calculations. Using a simple estimation of the truncation error involved in a 3-classes renormalization, a model-independent strategy is derived for choosing $`m`$ and $`p`$ and for extrapolating to the limit where $`m`$ and $`p`$ are infinite. In this way reliability and accuracy are further improved.

In order to pursue the development of the 3-classes approach, the fluctuations involving class-2 states of different blocks should be taken into account without increasing the dimension of the renormalized Hilbert space. This can be done by using a self-consistent intermediate Hamiltonian method based on a coupled-cluster approximation. Similar techniques are very successful in large-scale ab initio configurations interaction (CI) calculations on molecular systems . The model space $`𝒮`$ of the effective Hamiltonian coincides with the Hilbert space of the 3-classes DMRG approach (i.e., multiple class-2 states on different blocks are excluded). From the calculation of the ground-state wave function $`|\mathrm{\Psi }`$ within $`𝒮`$ one obtains relations between the wave-function coefficients corresponding to class-1 and class-2 states. Assuming a coupled-cluster expansion for $`|\mathrm{\Psi }`$, these relations allow to infer approximately the coefficients of the wave function beyond $`𝒮`$, i.e., when class-2 states are present on different blocks. These fluctuation effects can be then included as a dressing of the block interactions within $`𝒮`$. Since the dressed or intermediate Hamiltonian depends on $`|\mathrm{\Psi }`$, a set of self-consistent equations must be solved at each renormalization iteration. Notice, however, that this could be done within the Lanczos or Davidson iterative procedure in order to avoid a too large increase of computer time. The detailed formulation of this self-consistent approach is rather lengthy and beyond the scope of this paper. A summary may be found in the Appendix. While an efficient implementation of such an extension should further improve performance, it is also true that its complexity risks to limit possible applications to very difficult cases where a sufficiently large number of states $`m`$ and $`p`$ is not practicable or where a systematic extrapolation to the limit where $`m`$ and $`p`$ are infinite is not possible.

###### Acknowledgements.

Computer resources provided by IDRIS under project 960806 are gratefully acknowledged.    

## Appendix

The aim of this section is to outline a method for taking into account the effect of the neglected fluctuations between the class-2 block-states,

$$\mathrm{\Psi }_{i^{}jkl^{}}=\psi _i^{}^L\varphi _j\varphi _k\psi _l^{}^R|i^{}jkl^{}$$
(2)   
($`m<i^{},l^{}m+p`$) without increasing the dimension of the Hilbert space. In analogy with the theory of intermediate Hamiltonians , we intend to map the Hilbert space $``$ of the standard DMRG procedure ($`m+p`$ states per block) onto the 3-classes DMRG Hilbert space $`𝒮`$ (with $`m`$ class-1 states and $`p`$ class-2 states). An effective Hamiltonian is thereby defined acting on $`𝒮`$ as model space, which takes into account the effects of the outer space $`𝒬=𝒮`$.

Let us consider the secular equation of a state $`\mathrm{\Psi }_{abcd}𝒮`$

$$\underset{|ijkl𝒮}{}C_{ijkl}H_{ijkl,abcd}+\underset{|\alpha \beta \gamma \delta 𝒬}{}C_{\alpha \beta \gamma \delta }H_{\alpha \beta \gamma \delta ,abcd}=C_{abcd}E$$
(3)   
One would like to rewrite this equation in the following effective Hamiltonian form

$$\underset{|ijkl𝒮}{}C_{ijkl}H_{ijkl,abcd}+C_{abcd}\mathrm{\Delta }_{abcd,abcd}=C_{abcd}E$$
(4)   
since in this way the contribution of the states belonging to $`𝒬`$ is taken into account by a diagonal dressing $`\mathrm{\Delta }`$ of the Hamiltonian. $`\mathrm{\Delta }`$ should therefore satisfy the equations

$$\mathrm{\Delta }_{abcd,abcd}=\underset{|\alpha \beta \gamma \delta 𝒬}{}\frac{C_{\alpha \beta \gamma \delta }}{C_{abcd}}H_{\alpha \beta \gamma \delta ,abcd}$$
(5)   
Notice however that the coefficients $`C_{\alpha \beta \gamma \delta }`$ cannot be obtained from the diagonalization of the effective Hamiltonian $`H_{eff}=H+\mathrm{\Delta }`$ since the latter acts only on the model space $`𝒮`$. The main difficulty is therefore to find a sound evaluation of the coefficients $`C_{\alpha \beta \gamma \delta }`$. This problem can be solved using a coupled-cluster approximation.

Let $`S^L`$ and $`S^R`$ be the cluster operators associated to the left and right renormalized blocks and let $`\mathrm{\Omega }=\mathrm{exp}(S^L+S^R)`$ be the wave operator going from $`𝒮_0`$ to $``$, where $`𝒮_0`$ is generated by the class-1 blocks states only. $`𝒮`$ can then be decomposed as $`𝒮=𝒮_0+𝒮_R+𝒮_L`$, where $`𝒮_R`$ ($`𝒮_L`$) is the Hilbert space generated out of the product of class-2 states in the right (left) renormalized block and class-1 states in the left (right) renormalized block. The ground-state $`\mathrm{\Psi }`$ can be expanded as

$`\mathrm{\Psi }`$ $`=`$ $`{\displaystyle \underset{|ijkl𝒮_0}{}}C_{ijkl}|ijkl+`$ (8)   
$`{\displaystyle \underset{|i^{}jkl𝒮_R}{}}C_{i^{}jkl}|i^{}jkl+{\displaystyle \underset{|ijkl^{}𝒮_L}{}}C_{ijkl^{}}|ijkl^{}+`$   
$`{\displaystyle \underset{|i^{}jkl^{}𝒬}{}}C_{i^{}jkl^{}}|i^{}jkl^{},`$   
while the coupled-cluster approximation gives

$`\mathrm{\Psi }`$ $`=`$ $`{\displaystyle \underset{|ijkl𝒮_0}{}}C_{ijkl}|ijkl+`$ (11)   
$`{\displaystyle \underset{\genfrac{}{}{0pt}{}{i^{}}{m<i^{}m+p}}{}}{\displaystyle \underset{|ijkl𝒮_0}{}}S_{i^{},i}^LC_{ijkl}|i^{}jkl+{\displaystyle \underset{\genfrac{}{}{0pt}{}{l^{}}{m<l^{}m+p}}{}}{\displaystyle \underset{|ijkl𝒮_0}{}}S_{l^{},l}^RC_{ijkl}|ijkl^{}+`$   
$`{\displaystyle \underset{\genfrac{}{}{0pt}{}{i^{}l^{}}{m<i^{}m+p,m<l^{}m+p}}{}}{\displaystyle \underset{|ijkl𝒮_0}{}}S_{i^{},i}^LS_{l^{},l}^RC_{ijkl}|i^{}jkl^{}.`$   
Comparing equations 8 and 11 one obtains the definition of the cluster operators

$`C_{i^{}jkl}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{ijkl}{im,lm}}{}}S_{i^{},i}^LC_{ijkl},`$   
$`C_{ijkl^{}}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{ijkl}{im,lm}}{}}S_{l^{},l}^RC_{ijkl}`$   
and the evaluation of the unknown coefficients

$$C_{i^{}jkl^{}}=\underset{\genfrac{}{}{0pt}{}{il}{im,lm}}{}S_{i^{},i}^LS_{l^{},l}^RC_{ijkl}.$$
(12)   
These systems of equations are unfortunately over-defined in most practical cases ($`p>m`$). An optimal choice for $`S^L`$ is obtained by minimizing the square deviation

$`^L`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{i^{}jkl}{m<im+p,lm}}{}}(C_{i^{}jkl}{\displaystyle \underset{\genfrac{}{}{0pt}{}{i}{i<m}}{}}S_{i^{},i}^LC_{ijkl})^2.`$   
Thus,

$`{\displaystyle \frac{^L}{S_{i^{},i}^L}}`$ $`=`$ $`2\left(C_{ijkl}C_{i^{}jkl}{\displaystyle \underset{\genfrac{}{}{0pt}{}{i^{}}{i^{}<m}}{}}S_{i^{},i^{}}^LC_{ijkl}C_{i^{}jkl}\right)`$   
$`=`$ $`2\left(\rho ^L(i^{},i){\displaystyle \underset{\genfrac{}{}{0pt}{}{i^{}}{i^{}<m}}{}}S_{i^{},i^{}}^L\rho ^L(i^{},i)\right)=\mathrm{\hspace{0.33em}0},`$   
where $`\rho ^L`$ refers to the ground-state density matrix reduced to the left renormalized block. Finally, one obtains

$$S^L=\left(P_2^L\rho ^LP_1^L\right)\left(P_1^L\rho ^LP_1^L\right)^1$$
(13)   
where $`P_1^L`$ ($`P_2^L`$) is the projection operator onto the class-1 (class-2) left block-states. Similar equations hold for the right block.

It should be noted that the clusters operators depend on the ground-state of the dressed Hamiltonian, which depends itself on the clusters operators. The whole process should then be iterated up to self-consistency. Let us finally recall that self-consistent approximations are usually very efficient but do not necessarily satisfy the variational principle.