An introduction to numerical methods in low-dimensional quantum systems

An introduction to numerical methods in low-dimensional quantum systems

André Luiz Malvezzi

Departamento de Física, Faculdade de Ciências Universidade Estadual Paulista, Caixa Postal 473, 17015-970, Bauru, SP, Brasil

ABSTRACT

This is an introductory course to the Lanczos Method and Density Matrix Renormalization Group Algorithms (DMRG), two among the leading numerical techniques applied in studies of low-dimensional quantum models. The idea of studying the models on clusters of a finite size in order to extract their physical properties is briefly discussed. The important role played by the model symmetries is also examined. Special emphasis is given to the DMRG.

1 Introduction

From a mathematical point of view, quantum mechanics can be regarded in many respects as an eigenvalue problem. Typically, one has to to calculate eigenvalues and eigenvectors of Hamiltonians. Very often, especially in low temperature problems, the knowledge of the ground state and the first few excited states yields considerable insight into the physics of a given system. Symmetry properties and the known quantum numbers of the desired state can be used to reduce the Hilbert space. However, even a large reduction factor will eventually be overcome by the exponential growth. This means that requirements for memory and running time of relatively small systems can be prohibitive, even for the best computers, when a full diagonalization of a finite cluster is attempted. Typically, many degrees of freedom have to be integrated out of the original problem to make it accessible to present day computer capacities.

These notes consist of an introductory course to the Lanczos Method and Density Matrix Renormalization Group Algorithms(DMRG), two among the leading numerical techniques applied in studies of low-dimensional quantum models. Another important numerical technique, the Quantum Monte Carlo, is discussed by R. R. dos Santos. The numerical approach allows for a direct and unbiased calculation of physical properties for finite clusters, from which a phase diagram can be constructed.

For those outside the field of one- or two-dimensional quantum systems it may seem that working on problems at such low dimensionalities is completely irrelevant to three-dimensional reality. However, there are in fact many reasons for working in low-dimensional physics, in both statistical mechanics and other fields.

First, restricting the study to one dimension continues to demonstrate itself as an efficient and effective laboratory for the development and consideration of new theoretical ideas, many of which are intended for application to real higher-dimensional systems. Both thermal and quantum fluctuation effects are larger in lower dimensions.

Secondly, there are in fact real systems which display a high degree of low-dimensionality, at least in condensed matter physics.

Thirdly, one-dimensional systems share the possibility with higher-dimensional systems of undergoing quantum phase transitions as some parameter is varied while holding the temperature fixed at T = 0. In experimental situations, this parameter could be pressure or doping, for example, such as in the case of high-temperature superconductors where, upon the variation of doping, the ground state changes from antiferromagnetic to superconducting.

This manuscript is organized in the following way. In Section 2 the Hamiltonians for the Hubbard, t-J and Heisenberg models are defined. These models and their generalizations are aimed to describe the microscopic features of matter. Besides that, they present rich mathematical aspects. In Section 3, the idea of studying the models on clusters of a finite size in order to extract their physical properties is discussed. The important role played by the model symmetries is also examined. In Section 4 the Lanczos method is presented and some of its most popular generalizations discussed. In Section 5 we start by examining the basic ingredients of the DMRG method and presenting the two kinds of DMRG algorithms. Then, the way measurements of observables are performed in DMRG is discussed and, finally, the most important additions to the original method are briefly mentioned. Section 6 contains some final comments and the acknowledgments.

2 Models

In this Section we present some representative quantum models. These models and their generalizations have wide applicability in statistical mechanics as well as in condensed matter physics.

2.1 The Hubbard Model

A starting point for the microscopic approach to electron motion in crystals may be obtained by analyzing the energy levels of the atoms involved. Lattice models are used, because the atomic structure in the physical systems determines the possible places at which the electron can be found. If two neighboring atoms have overlapping orbitals with very similar energies, the orbitals can hybridize and allow electrons to travel from one atom to the other. On the other hand, the repulsion between electrons is very strong due to their charge. The simplest approximation for the interaction between the electrons is to restrict it to the case when both electrons are on the same site (or atom). On-site Coulomb repulsion and nearest-neighbour (n.n.) hopping are already the terms of the Hubbard model[4,5], which can be characterized by the Hamiltonian

where (c_js) is a creation (annihilation) operator for an electron of spin s = ­,¯ in a Wannier orbital at lattice site i and áijñ denotes n.n. pairs. In the sum, bonds áijñ are summed over only once each. Here t is the matrix element for tunneling from one lattice site to a neighbouring one, i.e. the overlap of n.n. electron wavefuctions. The letters h.c. denote the Hermitian conjugate of the immediately preceding term. The fermionic operators obey the anticommutation relations

and

Also, n_is = c_is is the number operator for electrons of spin s at site i, and U the Coulomb repulsion. The charge or number operator at site i is

The states of the model are given by specifying the four possible configurations of each site (its level can either be empty, contain one electron with either of two spins, or two electrons of opposite spins) on a lattice made of L sites.

However, the knowledge of the parameters t and U is not enough to characterize the system. One also needs to know:

• The dimensionality of the system, i.e., whether it is a one-dimensional (1D) chain, a two-dimensional (2D) plane or a full three-dimensional (3D) system.

• The geometry of the lattice. The lattice structure can introduce frustration, and its symmetry properties greatly influence the behaviour of traveling electrons.

• The boundary conditions (BC), i.e., whether the lattice is open, closed or has some special condition at its surface.

• The filling or density,

where á...ñ denotes a ground state expectation value.

• The temperature T which is the third energy scale along with t and U.

The first term in Hamiltonian Eq. (2.1) is called hopping term and the second one Coulomb term. The hopping term alone can be shown to lead to a conventional band spectrum and one-electron Bloch levels in which each electron is distributed throughout the entire crystal (a metal). The Coulomb term alone would favor local magnetic moments, since it suppress the possibility of a second electron (with oppositely directed spin) at singly occupied sites (an insulator). When both terms are present the competition between them brings about a transition between the metallic phase and the Mott insulating phase[6].

2.2 The t – J Model

If the Hubbard model is considered in the limit where U/t is large, the strong coupling limit, the number of doubly occupied sites is small. This leads to the derivation of effective models, most prominently the t - J model[7]. In the t– J model, the Hubbarb model complexity is reduced by projecting out the states with double occupancy. Using this procedure one gets the t - J Hamiltonian

where

is the electron spin operator at site i, with

being the Pauli matrices. The tilde on the c-operators in the hopping term refers to the fact that a creation operator cannot introduce an electron on a site where another electron, even of the opposite spin, is already located. Formally, one would express this as

In the derivation of the t - J model from the Hubbard model, one takes into account intermediate double occupancies between states which are suppressed by their high energy but help to lower the kinetic energy by making the electrons more mobile. In such processes, an electron hops on an already occupied site, then either electron can hop back to the original site. Double occupancies are forbidden in the t - J model, but the physics can be included by adding an additional term. Since the intermediate processes in the Hubbard model are only possible if the electrons on the adjacent sites have opposite spin, this term can be described by an interaction that favors the singlet state compared to the triplet state. This is precisely what the second term in Hamiltonian Eq.(2.6) does. To show this, we note that _i ·_j has two eigenvalues: if the two electrons are in a singlet state, the eigenvalue is -3/4, and if the electrons form a triplet, the energy is 1/4. If the effect of the n_in_j-term is included, one gets the following energies from the J-term:

one or both sites unoccupied

0

spins on both sites forming a triplet state

0

spins on both sites forming a singlet state

– J.

Therefore, the term effectively lowers the energy for states in which two electrons with opposite spin are on adjacent sites. From second order perturbation theory of the Hubbard model we obtain J = 4t²/U[2].

2.3 The Heisenberg Model

In the case of half-filling (n = 1), when n_i º 1 for all sites, hopping becomes impossible in the t - J model. In addition, the n_in_j-term is reduced to a constant. Thus, both terms can be neglected and the t-J model at half filling is a spin 1/2 Heisenberg model,

It is important to remember that this is an effective model, i.e., although the spins in the model are localized they are meant to describe a system of mobile electrons. The _i ·_j interaction is called a spin exchange interaction.

All these models were originally designed to describe physical systems. The Heisenberg model was introduced before the Hubbard model and is not restricted to just spin 1/2 systems. The restriction in the above case resulted from the fact that electrons are spin 1/2 fermions.

A possible generalization for these models is to extend the n.n. summation to more distant neighbours.

3 The study of finite-size systems

Having set up Hamiltonians which are believed to contain the physics we are interested on, we are left with the formidable task of calculating measurable quantities to gain some understanding of the models.

Here we choose to work with two unbiased numerical techniques, namely, Lanczos and DMRG. These methods are unbiased (i.e., in contrast to mean-field based approaches) in the sense that they do not make an initial assumption on the nature of the ground state of the system. However, as we will see later, each of these numerical methods has its limitations.

In studying finite-size systems the general idea is to construct a matrix representation of the Hamiltonian for a given number of sites (system size). The Hamiltonian is then diagonalized in order to obtain the spectra and calculate measurable quantities, e.g., spin and charge correlation functions. We repeat this calculation for systems of different sizes and extrapolate the results toward an infinite size system (i.e., towards the thermodynamic limit). This agenda is usually well succeeded provided that we have enough data for a reliable extrapolation. The number of lattices with different sizes needed for the extrapolation procedure to converge depends heavily on the model being studied and even on the set of parameters being used for a given model. The Finite Size Scaling theory[8] is specially useful when critical behaviour is present since, in this case, a strong dependence of the physical quantities on the system size is expected[9]. However, in many cases a fairly good ideia of the properties in the thermodynamic limit can be achieved by examining just a few system sizes. The models we are dealing with involve many parameters (e.g., n, t, U, and J in the models of Section 2). In order to build up a phase diagram we have to systematically cover the model parameter space to see how the physical quantities depend on each of these parameters.

We can roughly divide the above task in two steps:

• Build up a representation for the Hamiltonian, diagonalize it, and calculate the measurable relevant quantities.

• Interpret the results and construct the phase diagram.

The second step strongly depends on the model being studied and on what kind of physics we are looking at. If the subject being studied is an advanced topic in physics (as it usually is), then a considerable experience in that field of research might be necessary to carry that step out. On the other hand, the first one is, in principle, much simpler and involves the knowledge of basic concepts in quantum mechanics and numerical analysis; topics covered in most undergraduate courses in physics. The Lanczos and DMRG are tools for carrying out the first step. To understand why such special techniques are needed let us suppose the Heisenberg model is to be studied on a lattice of L sites. Each site has two possible states: spin up and down. A lattice with L sites has 2^L states and this is the dimension of the Hamiltonian matrix. Similarly, for the t - J and Hubbard models we have 3^L and 4^L, respectively. Due to this exponential growth with L even small lattices (typically 10 sites or so) generate Hamiltonians too big to be handled by present-day computers using standard diagonalization algorithms.

3.1 Symmetries

In the process of constructing a representation for the Hamiltonian it is very useful to take advantage of the model symmetries. Many models, including those presented in Section 2, exhibit conservation of total spin number, total spin in the z-direction and total charge, i.e.,

where H is the model Hamiltonian and

In addition, these operators also commute with one another, i.e.,

so that the eigenvalues of H, , S^z, and N are simultaneous good quantum numbers, which we will denote simply by E, S(S+1), S^z, and N respectively.¹ In the numerical treatment of a given model it is possible to consider eigenstates which simultaneously diagonalize H and all operators associated with its symmetries. We do this by choosing to work in a representation in which the symmetry operators are always diagonal, selecting a subspace or sector of Hilbert space with particular eigenvalues of those operators, and diagonalizing H in this particular sector. As exemplified below, total spin in the z-direction and total charge are easily implemented. However, total spin quantum number ()² is much more difficult to specify, and also quite difficult to measure, since in terms of the fundamental operators c_is it is extremely non-local:

where i and j both run throughout the entire lattice, i.e., they do not just represent n.n. terms. However, in cases in which choosing the direction of quantization in the z-axis is arbitrary, i.e., in which there is full rotational (or SU(2)) symmetry, the eigenvalues of H, N, and ² will be independent of the eigenvalues of S^z, which may assume any value ranging from – S to S. This can be quite useful in determining the total spin quantum number S of the ground state if it is only possible to specify the projection S^z in the method of investigation.

Let us denote by E(S^z) the ground state in the sector specified by the spin projection S^z. If E(S^z) < E(S^z+1) and E(S^z) = E(S^z-1) = ... = E(), where

then the absolute ground state, (i.e., that in the Hilbert space unrestricted by the specification of S^z) contains a state of spin S^z and no states of any higher spin.

Strongly correlated electron models also often exhibit a particle-hole symmetry. This symmetry relates the creation of an electron to its destruction in the following way. Consider the transformation

Under this transformation, the n.n. hopping terms transform according to

remaining unchanged under this transformation. However, the number operators transform according to

and similarly the conduction electron spin operator transforms according to

Here R represents a spin rotation by p about the y-axis. If this particle-hole symmetry applies to the full Hamiltonian, i.e. not only to its non-interacting kinetic part as shown here but also its interaction terms, then there will be a one-to-one correspondence between eigenstates with quantum number (E,S,S^z,N) and (E,S,-S^z,2L-N). In particular, if we wish to determine the properties of a system away from half-filling, we need only do this below half-filling and the physics above half-filling will be identical.

In order to reduce finite size effects it is common to work with periodic boundary conditions (PBC) to eliminate the lattice surface in given a direction. In this case, the system gains translational invariance in that direction². This symmetry expresses itself through the commutation of H with the translational operator T, which can be defined as

where

and d_i is the state at site i (e.g., for the Hubbard model d_i can be either one electron with spin up or down, two electrons with opposite spins, or an empty site). In Eq. (3.15) we have assumed, for the sake of simplicity, a one-dimensional lattice of size L. In higher dimensions translational invariance in each direction can be treated separately. The eigenvalues of T are eik, k = 0,1,2,...,L-1, which we will label by the number k. An eigenstate f_k of T with eigenvalue k is given by

where U = ei^kT and C a normalization constant. If f_k = 0 then is not possible to construct a state with momentum k from |d₁d₂...d_L-1d_Lñ. Note that S, N, S^z, and k are simultaneous good quantum numbers.

Another symmetry shared by all the models introduced and Section 2 and many others is the charge-conjugation symmetry. This symmetry implies a one to one correspondence between eigenstates with quantum number (E,S,S^z,N) and (E,S,-S^z,N) and is present if the Hamiltonian commutes with the parity operator, whose effect is to flip the particle spin (|­¯00­ñ® |¯­00¯ñ). In the presence of charge-conjugation symmetry we can choose to work only with S^z £ 0 or only with S^z ³ 0.

Let us consider the t - J model on a chain with four sites under PBC. Each site can be in one of the following states: one electron with spin up (|­ñ) or down (|¯ñ), or empty (|0ñ). The dimension of total Hilbert space is 3^L = 3⁴ = 81.. We divide the Hilbert space in sectors labeled by the quantum numbers N, S^z, and k. Below, we denote each sector by [N,S^z,k,(dimension of the sector)] and write down the states for some illustrative cases.

• [0,0,0,1]

• [1,0.5,0,1] [1,0.5,1,1] [1,0.5,2,1] [1,0.5,3,1]

• [2,0,0,3]

• [2,0,1,3]

• [2,0,2,3]

• [2,0,3,3]

• [2,1,0,2]

• [2,1,1,1]

• [2,1,2,2]

• [2,1,3,1]

• [3,0.5,0,3] [3,0.5,1,3] [3,0.5,2,3] [3,0.5,3,3]

• [3,1.5,0,1] [3,1.5,1,1] [3,1.5,2,1] [3,1.5,3,1]

• [4,0,0,2] [4,0,1,1] [4,0,2,2] [4,0,3,1]

• [4,1,0,1] [4,1,1,1] [4,1,2,1] [4,1,3,1]

• [4,2,0,1]

The above example gives an ideia of how helpful symmetry implementation can be. Indeed, a further reduction in the Hilbert space can be achieved in some cases by considering a lattice reflection symmetry, namely [W,H_{t J}] = 0, where

The eigenvalues of W are +1 and -1. For instance, the sector [2,0,0,3] can be broken into two subspaces, indexed by W = -1 ([2,0,0,-1,1]) and +1 ([2,0,0,1,2]), and given by:

• [2,0,0,-1,1]

• ([2,0,0,1,2])

In higher dimensions, besides reflection with respect to different axes, rotations by several angles are also available. Note that these operations might or might not commute with the translational operator depending on the value of k. For instance, reflection symmetry can not be employed to break sector [2,0,1,3] in sub-sectors.

By implementing the model symmetries we can significantly reduce the computational effort required for diagonalizing the Hamiltonian. In some cases we can study lattices large enough to reveal bulk properties using standard routines to diagonalize each sector of the Hamiltonian. For instance, a 4×4 Heisenberg system can be fully diagonalized if most symmetries discussed above are implemented[10]. Since this approach yields the full spectrum we can construct the exact partition function for the system and, therefore, obtain the thermal behaviour exactly at arbitrary temperature.³

4 Lanczos Method

In this Section, an algorithm is presented which allows us to determine numerically the ground state and some excited states for Hamiltonian operators on finite clusters. The basic idea of the Lanczos method[11,12] is that a special basis can be constructed where the Hamiltonian has a tridiagonal representation. Once in this form the ground state of the matrix can be found easily using standard library subroutines such as Numerical Recipes or IMSL.

The tridiagonal matrix is constructed iteratively. First, it is necessary to select an arbitrary normalized vector |y₀ñ in the Hilbert space of the model being studied. The overlap between the actual ground state |Y₀ñ, and the initial state |y₀ñ should be nonzero. If no a priori information about the ground state is known, this requirement is usually easily satisfied by selecting an initial state with randomly chosen coefficients in the working basis that is being used. If some other information of the ground state is known, like its total momentum and spin, then it is convenient to initiate the iterations with a state already belonging to the subspace having those quantum numbers (and still with random coefficients within this subspace).

After |y₀ñ is selected, we define a new vector by applying the Hamiltonian H to the initial state. Subtracting the projection over |y₀ñ, we obtain

that satisfies áy₀|y₁ñ = 0. Now, we can construct a new state that is is orthogonal to the previous two as

It can be easily checked that áy₀|y₂ñ = áy₁|y₂ñ = 0. The procedure can be generalized by defining an orthogonal basis recursively as

where i = 0, 1, 2, ..., and the coefficients are given by

supplemented by b₀ = 0, |y-1ñ = 0. In this basis, it can be shown that the Hamiltonian matrix becomes,

i.e., it is tridiagonal as expected. Once in this form the matrix can be diagonalized using standard library subroutines. However, note that to diagonalize completely the model being studied on a finite cluster a number of iterations equal to the size of the Hilbert space (or of the subspace under consideration) is needed. In practice this would demand a considerable amount of CPU time. However, one of the advantages of this technique is that accurate enough information about the ground state can be obtained after a small number of iterations (typically of the order of ~ 100 or less). The ideia behind Lanczos method is a systematic improvement of a given variational state that is used to represent the ground state [14]. The procedure just described assumes that H is a an hermitian matrix. If that is not the case then a generalized algorithm is needed[15].

The eigenvalues of (4.5) steadily approach the lowest eigenvalues of H and its eigenstates are expanded in the Lanczos basis |y_iñ. Each state |y_iñ is represented by a large set of coefficients, when it is itself expanded in the basis selected to carry out the problem (e.g, the basis used in example of Section 3). Thus, in practice, it is not convenient to store each one of the |y_iñ vectors individually, since such a procedure would demand a memory requirement equal to the size of Hilbert space sector multiplied by the number of Lanczos steps. A simple solution to this problem consists of running the Lanczos subroutine twice. For instance, if |Y₀ñ = å_i f_i |y_iñ, then in the first run the coefficients f_i are obtained, and in the second the vectors |y_iñ are systematically reconstructed one by one and used to build up |Y₀ñ in the original basis.

A common difficulty with the Lanczos method is that finite precision arithmetic causes the |y_iñ to lose their orthogonality. A consequence of that is the appearance of spurious eigenvalues (and eigenstates) in the spectra. One fix is to repeatedly reorthogonalize, which is too costly since, in this case, all |y_iñ have to be stored. Another is to partially reorthogonalize, and a third option is to ignore the problem. We can choose the latter if just the extremal (lowest or highest few) eigenvalues are needed, which is often the case. A final way of avoiding the problem is to stop the Lanczos whenever the problem starts to appear, calculate the ground state, and use it as an initial state |y₀ñ to restart the Lanczos procedure. This resets the basis but keeps the information from the previous Lanczos running. By pushing this ideia to its limit, we can always perform just two Lanczos steps (i.e., work just with |y₀ñ and |y₁ñ), diagonalize a 2×2 matrix, and use its lowest eigenstate as a new |y₀ñ. This is the so-called modified Lanczos[16,17]. The modified Lanczos converges more slowly than the original Lanczos but has the convenience of having the ground state always at hand. An even more pedestrian technique is the power method[19], which consists of successively applying the Hamiltonian to the initial state until all excited states are filtered out and only the ground state remains. This procedure is the slowest in speed of convergence, but in simple problems is enough and easy to program.

One of the greatest appeals of the Lanczos method is the possibility of calculating dynamical properties of a given Hamiltonian in finite clusters[18,12]. The technique permits accurate calculations of energy and momentum dependent dynamic correlation functions which are observable in scattering experiments, such as Neutron Scattering (spin dynamics) and Photoemission Spectroscopy which measures the spectral function of the system[13].

Since Lanczos yields a number of excited states another interesting possibility is the calculation of finite-temperature quantities. Examples of successful attempts in this front are Refs. [20] and [21].

The main limitation of Lanczos technique is the size of the clusters that can be studied. Recently, attempts have been made to reach larger cluster. The basic ideia is the following. If |f_iñ is a complete basis we are working with, then the ground state can be formally represented as

In general, all |f_iñ contribute significatively for the sum but, in some cases, it might happen that several states have very small weight g_i. In fact, for most models studied a very small percentage of the states in the basis dominate the sum in Eq. (4.6). This suggests that useful results can still be obtained if a large part of the Hilbert space is simply discarded. The truncated Lanczos[22] implements this idea by systematically enlarging and reducing the working Hilbert space in a controlled way and by keeping a fixed number of states (typically a million or so) in the basis. The number of states kept represent a very small percentage of the total Hilbert space but the algorithm is designed to search for and keep the dominant states in sum Eq. (4.6). It is also important to choose a basis |f_iñ which uses as much information about the system as possible. For every physical systems we should use a different basis depending on physical insight we have about the system. For instance, if we know that there is a tendency to dimerization then it is convenient to construct a basis of dimers. A well chosen basis leads to a better representation of |Y₀ñ for a given number of states kept in the truncation procedure.

The efficiency of the truncation technique depends heavily on the model being studied. In particular, it seems suitable for problems with gaps in the spectrum (like a spingap)[23].

5 Density Matrix Renormalization Group Algorithms

5.1 Overview

The basic agenda to overcome the system size limitations is to use a basis in which the ground state can be represented by only a few base states. In other words, a procedure must be found to identify or construct the important states and neglect or discard all others so that the piece of the Hilbert space one operates on remains small. The truncated Lanczos, briefly discussed in the previous Section, is a possible approach to this agenda. It has the advantage of working with a basis formed by states that have an intuitive meaning so that the results can be easily interpreted. In addition, dynamical information can be obtained without difficulty. The DMRG technique we are about to discuss is an alternative approach to that agenda. The innovation of the DMRG is that it does not hold on to a specific basis, but optimizes the basis it uses in the steps leading to the calculation of the ground state. A disadvantage of the DMRG method is that the base states chosen by the algorithm are not intuitive, and the description of the state requires the measurement of observables. For the measurement process, one needs a representation of the operators in the current basis. Consequently each operator that needs to be measured must be stored, and every time the basis is changed all of them have to be transformed. This is expensive in time and memory. Another disadvantage is that dynamical information cannot be easily obtained.

Historically, DMRG has its roots in the renormalization group approach pioneered by Wilson[24]. Considering a real-space blocking version for lattice systems of Wilson's original approach, the basic idea is to start with a small system that can be handled exactly. The system size is then increased without increasing the size of the Hilbert space until the desired system size is reached.

Increasing the system without increasing the Hilbert space is typically done in two steps:

• The system size is increased, and as a consequence, the Hilbert space grows at the same time.

• The Hilbert space is truncated to its original size keeping the system size constant.

To characterize such a renormalization procedure two basic questions have to be answered:

• How is the enlargement done ?

• Which criterion to apply in the second step to differentiate between the states we will keep from those we will discard ?

In Wilson's approach, we start with blocks of a certain (small) size. In the first step, two such blocks are linked to form a block which is twice as large. The Hamiltonian of this larger block is then exactly diagonalized and its eigenstates are used as base states. The criterion for keeping states is their energy, and only eigenstates whose energy lies below a certain threshold are kept. The states which are kept characterize the new block that is again linked to an identical block, and the process is iterated.

This approach proved to be very effective for the Kondo model[24]. However, for other strongly correlated systems like those in Section 2 it was not successful[25]. The main reason for this failure lies in choosing the block eigenstates as the states to be kept. Since the block was not previously connected to the rest of the system ( another identical block in the case above ) its eigenstates have inappropriate features at the block ends, making them a poor choice as a basis to represent the ground state of a larger system, formed by putting together two (or more) blocks. This problem was pinpointed by White and Noack in Ref. [1] and an attempt was made to fix it by combining eigenstates from several different blocks under various BC. Let us see how DMRG fixes this problem.

5.2 Enlargement and Reduction in the DMRG Procedure

As mentioned in Section 5, the two most important characteristics of a renormalization procedure are the way the system grows and how the decision is made on which states are kept in the Hilbert space by the truncation step. In this Section, these elements of the DMRG procedure will be discussed. How these elements are used in the global DMRG algorithm will be the subject of the subsequent Sections.

Figure 5.1 shows the most important structures used in the DMRG algorithm. The elementary unit is a site, and is described by the states d_i(i = 1,¼,D), in which the site can be found⁴. A block B(l,m) consists of a number of sites l and its Hamiltonian H_B contains only terms involving the sites inside the block. To represent B(l,m) and H_B we associate an m-dimensional basis where m is in general smaller than the full Hilbert space of the block. The states in the basis are grouped in symmetry sectors labeled by a set of quantum numbers (e.g., S^z and N), which makes H_B a block-diagonal matrix. We also store the matrix elements of H_B between these states. The block is grown by adding a site to it, and together they form the enlarged block B^e. If |b₁ñ¼|b_mñ and |d₁ñ¼|d_Dñ represent, respectively, the basis of block and a site then the basis of the enlarged block can be constructed from the direct product

The dimension of the Hilbert space for B^e is the product of the dimensions of the Hilbert space of B(l,m) and a site, i.e, m×D. A possible mapping of i and j-values onto k-values in Eq. (5.1) is k = (i - 1)D + j.

The next step in the DMRG method is the formation of the superblock Hamiltonian (Fig. 5.2). The superblock consists of two enlarged blocks connected to each other. In Fig. 5.2 open boundary conditions (OBC) are applied to the superblock. These BC are the most widely used in DMRG for it yields the best results for a given computational effort. We will discuss PBC latter.

The DMRG method focus on a single eigenstate of the superblock Hamiltonian (usually the ground state), called the target state, which is used to construct the density matrix⁵. The ground state of the superblock is calculated (using Lanczos or any other method). We then eliminate the states from the basis of the enlarged block that contribute the least to the ground state of the superblock. To calculate those, the density matrix is used.

The concept of the density matrix was developed in statistical mechanics[26] by considering the problem of a system in contact with a much larger bath. The ground state of the universe, i.e. system and bath, is known, and the question is which states of the system contribute the most to this ground state. This is what the density matrix can tell us. One can express the ground state of the universe (the superblock) in a basis that is the tensor product of the basis vectors of the system (one of the enlarged blocks) and the bath (the other enlarged block),

Hence many of the eigenstates of the system contribute to the one ground state of the universe. The density matrix of the system is given by

We show an actual example of such a calculation below. The density matrix has the same dimension and block-diagonal structure of the Hamiltonian H_e, for the enlarged block. If we denote by |u_a ñ (a = 1,¼,m×D) the eigenstates of r and by w_a its eigenvalues then å_a wa = 1 and w_a is the probability of the system being in the state |u_a ñ given that the universe is in the state |Y₀ñ.

This is the information we need to decide which states to keep in a renormalization group approach. In order to make an optimal decision of which states to discard and which to keep, it is a good criterion to consider the weight w_a of the states in the ground state of a larger system, which we eventually want to describe. We must order the |u_a ñ by their eigenvalues in a decreasing order and use the first m of those states with largest eigenvalues to form a new basis for the enlarged block B^e, which will then become B(l+1,m). In symbols,

where the rows of the m×(m×D) matrix O are formed by the |u_a ñ previously selected. The change of basis in Eq. (5.4) renormalizes the Hilbert space, cutting its size back to m. Constructed in this way, the blocks are being prepared to be connected to another block in the next step, when a new superblock will be formed. By using the density matrix states we somehow look into the future and adapt the block for it. Besides H_B we also need to store other operators representing the sites at the border of the block. These operators are necessary to construct the interaction between the block and the site, when forming the enlarged block and also needed to be transformed according to Eq. (5.4).

To illustrate the DMRG steps of enlargement and truncation, a full DMRG step for the antiferromagnetic spin 1/2 Heisenberg chain will be performed[27]. The starting point is a block B(1,2) of a single site. The possible states of the single site are

For convenience the up/down basis is chosen. The basis itself is not stored. The only data that is stored are the operators needed to progress the algorithm namely, the operators needed to build the Hamiltonians for the enlarged block and the superblock.

For one isolated site without external fields the Hamiltonian is zero. Since the up/down basis was chosen, the other operators are the spin matrices given by

To build the enlarged system, another site is added. In this case the basis of the block is the same as the basis of the added site,

and the operators look the same as those of the block. Thus, the basis of the enlarged block is, by Eq. (5.1),

The Hamiltonian H_e for the enlarged block B(2,4) has non-zero elements, and describes the interactions of the sites in B(2,4). H_e consists of the H_B, describing the interactions within the block, and the interactions between the rightmost spin of the block and the new site. In the above basis the Heisenberg Hamiltonian of the enlarged block is

and it looks as follows:

In Eq. (5.9) the indices b and d refer to the operators acting on the Hilbert space of the block and the site, respectively, and I is the unit matrix. In this first step we have m = D = 2 but, as the block grows in size in the following steps we will have m > D. Note that only representations for the Hamiltonian of the block and for the operators S⁺, S^-, and S^z of the rightmost site of the block and the new site are needed to construct the enlarged block.

The superblock is constructed by taking the enlarged block as the left block and connecting it to another enlarged block on the right (Fig. 5.2). In the so-called infinite size method, discussed in the next Section, the right block is the same as the left block, only spatially reflected so that the site last added to the left block is connected with the site added last to the right block. The rightmost site of the left block becomes the leftmost of the right block.

In addition to the Hamiltonians of the enlarged blocks, one needs representations of the spin operators of the rightmost site of the enlarged block. In order to construct representations for S⁺, S^-, and S^z in the basis of the enlarged block we have to calculate the tensor product of the unit matrix of the block Hilbert space and the operator in the representation of the basis of the rightmost site. For instance, the ()_e matrix, the S⁺-operator of the spin on the rightmost site in the basis of the enlarged block, is given by

Representations for ()_e and ()_e are obtained in a similar way. The basis for the superblock is the tensor product of the bases from the two enlarged blocks being connected:

In general, |ñ and |ñ are distinct basis. Assuming that we want to calculate the ground state properties, it is possible to exploit S^z conservation and the fact that the ground state belongs to the subspace with S^z = 0. Therefore, we can concentrate only on states in this symmetry sector:

The Hamiltonian of the superblock consists of three parts: the internal Hamiltonians of the two enlarged blocks and the exchange arising from the spin interacting at the connection between them:

where the prime refers to operators in the Hilbert space of the second enlarged block forming the superblock. We now build a representation for H_s in the basis |⁽⁰⁾ñ:

The ground state energy of is E₀ = (1/4)(3+2) and the corresponding eigenvector is

In order to decide which states of the left enlarged block are the most important for the ground state of the superblock one uses the density matrix given by Eq. (5.3). Combining Eqs. (5.12) and (5.13) we obtain

This allows us to identify the coefficients a_ij in Eq. (5.2) (and in Eq. (5.3)) and they are all zero except for a₁₄, a₂₂, a₂₃, a₃₂, a₃₃, a₄₁. For the density matrix we get

Note that r and H_e (Eq. (5.10)) share the same block-diagonal structure. The eigenvalues of r are (1/12)(2 + ) » 0.02 for each of the triplet states and (21 + 12)(12(2 + )) » 0.94 for the singlet state. The base states are then ordered according to the size of the respective eigenvalues, with the singlet state (largest eigenvalue) coming first. The transformation matrix O in Eq. (5.4) is given by

After determining the basis and the transformation, the representations of all operators used to describe the enlarged block are changed to the new basis. Applying the transformation to the H_e (Eq. (5.10)) leads to

In the present simple case H_e and r have the same eigenvectors. That is the reason why the above transformation diagonalizes H_e. The same transformation is done with the other operators that will be needed for future calculations. One example is the S⁺-operator, which has the following representation in the new basis

Note that, even though a site has been added to block B(1,2) to form block B(2,2), the dimension of Hilbert space did not change, due to the truncation performed. The states kept in the truncation are those with higher probability to be found in the ground states of the superblock system.

In this example we have performed the truncation in the Hilbert space in order to illustrate the procedure. In a practical calculation the system in our example would be too small to already start the truncation. Usually, we know from the beginning how many states will be kept. Thus, in the first steps the block is grown (sites are added) without truncation until the number of states needed to describe the block becomes larger than the number of states we want to keep⁶. If we had, for instance, decided to keep m = 20 states in the computation, the chain would be grown to a size of 5 sites with 2⁵ = 32 states, which is the first block size with Hilbert space dimension exceeding 20. Then, a truncation would create a block B(5,20) from an enlarged block with 5 sites and in the following steps all blocks would have dimension 20, even though they represented a increasing number of sites.

We often use the sum of the density matrix eigenvalues of the discarded states (1 -

_{= 1} w_a) as a measure for the severity of the truncation. The goal is to keep this number as small as possible. In many cases it has been found that this number is roughly proportional to the error in the energy[28]. The proportionality factor is of course model-dependent. In doped fermionic models, we need to keep more states to achieve a good accuracy than in a spin model. Even for a given model the accuracy for a given truncation may depend on the parameters being used in the calculation (e.g., couplings and symmetry sector). For instance, close to a phase transition line or inside a critical (massless) phase the strong quantum fluctuations tend to reduce accuracy. In the example above we discarded two of the triplet states leading to a truncation error of 0.04, which is unacceptably high. Truncation errors in actual calculations are usually kept smaller than 10^-4.

In this Section the focus was on one DMRG step. Enlargement of the block by adding a site, the formation and diagonalization of the superblock, the calculation of the density matrix, and the truncation procedure were discussed. In the next Sections we describe how several DMRG steps are combined to calculate the properties of a given model.

5.3 The infinite system algorithm

The first implementation of the DMRG method was the infinite system algorithm[2,3]. The goal was to use DMRG's advantage to decouple the system size and the size of the Hilbert space and calculate ground state energies of large systems, i.e., system sizes that are unreachable for exact diagonalizations, eventually converging to the thermodynamic limit.

In the infinite system algorithm, the left enlarged block is connected to its own mirror image on the right side, so that the number of sites of the superblock is increased by two in each step. Growing the block and truncating the Hilbert space is done as explained in Section 5. Schematically the algorithm can be described as follows:

1. Grow the chain to a size in which its Hilbert space dimension is just larger than m, the number of states to be kept. This is the first enlarged block.

2. Form the superblock by adding an identical enlarged block on the right such that the sites, which were added last, are next to each other.

3. Diagonalize the superblock, calculate the density matrix with respect to the left enlarged block.

4. Change basis of the left enlarged block to the eigenbasis of the density matrix, keep only the m states with the largest density matrix eigenvalues. The transformed, truncated left enlarged block becomes the block for the next iteration.

5. This block is enlarged, i.e., a site is added on the right.

6. Continue with step 2 until convergence is reached.

In Fig. 5.3, two successive iterations of the infinite system algorithm are shown. The starting point is a block with l sites that is described by a basis with m states. Enlargement and construction of the superblock (step 2), leads to the situation portrayed in Fig. 5.3.a.

After diagonalizing the superblock, finding the transformation, and truncating the enlarged block according to its density matrix (steps 3 and 4), one arrives at the situation depicted in Fig. 5.3.b. The new block describes a chain with l + 1 sites, but uses a basis with only m states. Enlarging the block (step 6) and building the superblock (step 2) leads to the situation in 5.3.c, and the procedure is repeated.

From the computational point of view the most difficult part is the calculation and the subsequent diagonalization of the superblock Hamiltonian. The diagonalization can be done with the Lanczos method but any other method (e.g., the Davidson algorithm[29]) can be used. The computational effort depends on the size of the Hamiltonian matrix and accuracy needed for the ground state. Since the superblock Hamiltonian is block-diagonal, diagonalizing only the sector that has the proper quantum numbers reduces the matrix by a factor that depends on the superblock, block, and model. Typically it is of order 10-20. Since the total Hilbert space for the superblock has dimension (D×m)², the most important determining factor for the size of the Hamiltonian is the number m of states kept in the block. In actual calculations m is typically a few hundred of states and is limited due to restrictions in computer memory and CPU-time. As we work with different models, the size D of the Hilbert space for a single site also affects the computational effort needed to reach a given accuracy in the results.

To summarize: in the infinite system algorithm, the system size is increased in each step while the number of states kept to describe the blocks is constant. The goal is to grow the chain to a long-enough length, so that the energy and short range correlations around the center have converged. The convergence is checked by keeping track of the difference DE₀ between the ground state energy of the superblocks in two sucessive steps⁷.

5.4 The Finite System Algorithm

In the finite system algorithm the goal is no longer to reaches the thermodynamic limit, but rather to restrict ourselves to a finite system size L. In the beginning, until the superblock size reaches the system size, the algorithm is identical to the infinite size algorithm. When the system size reaches L, i.e. the enlarged block has L/2 sites, the left block is grown further but on the right an enlarged block with a smaller number of sites is used in order to keep constant the number of sites in the superblock. As soon as the decreasing size of the right block reachs a single site the procedure is stopped. One such iteration in which the left block has been calculated for all possible sizes, nearly up to L, is called a sweep over the system. After one sweep is done, we start all over with a small left block. However, from now on, the information about the best representation of the right block that complements the left block to the desired system size L is now present since it has been calculated in the previous sweep. When the optimal basis for a specific size of the left block is determined with DMRG, the result is stored and used in the next sweep as best possible guess for the optimal states describing the right block. The steps of the finite system algorithm are compiled in the following list:

1. In the first sweep use the infinite size algorithm until the superblock size reaches the chain size L under investigation. After every truncation save all operators of the reduced block to disk.

2. Enlarge the left block size l + 1.

3. Read a block of size L - l - 2 from disk; this is the right block.

4. Enlarge the right block to the size L - l - 1.

5. Form the superblock from right and left enlarged blocks.

6. Diagonalize the superblock, calculate the density matrix with respect to the left enlarged block.

7. Change the basis of the left enlarged block to the eigenbasis of the density matrix, keep only the m states with the largest density matrix eigenvalues, save the block with the basis to disk. The transformed, truncated left enlarged block becomes the left block for the next step.

8. Continue with step 2 until the right block becomes a single site

9. If the right block is a single site, begin a new sweep over the system, i.e. construct a superblock with a left enlarged block containing two sites. Continue with step 3 until convergence is reached.

To illustrate the progress of the algorithm, two general steps are portrayed in Fig. 5.4.

At first glance Fig.5.4 looks like Fig.5.3. In fact the algorithms are very similar. The left bock is enlarged, the superblock is constructed and diagonalized (Fig.5.4.a). Then the density matrix for the left enlarged block is constructed and the block is reduced (Fig.5.4.b). This reduced block is then again enlarged, the superblock is built up and so on. The difference between the two algorithms is that the right block in Fig. 5.3 has the same number of sites l as the left side. In the finite system algorithm the right enlarged block always complements the left one to the target size L and thus becomes shorter, while the later grows.

As an example, let us consider a calculation for the Heisenberg chain with L = 16 sites and a truncation of m = 24 states. The following superblocks have to be formed:

• Initial sweep:

[B(1,2),B(1,2)] [B(2,4),B(2,4)] [B(3,8),B(3,8)] [B(4,16),B(4,16)]

[B(5,24),B(5,24)] [B(6,24),B(6,24)] [B(7,24),B(7,24)] [B(8,24),B(6,24)]

[B(9,24),B(5,24)] [B(10,24),B(4,16)] [B(11,24),B(3,8)] [B(12,24),B(2,4)]

• Following sweeps:

[B(1,2),B(13,24)] [B(2,4),B(12,24)] [B(3,8),B(11,24)]

[B(4,16),B(10,24)] [B(5,24),B(9,24)] [B(6,24),B(8,24)]

[B(7,24),B(7,24)] [B(8,24),B(6,24)] [B(9,24),B(5,24)]

[B(10,24),B(4,16)] [B(11,24),B(3,8)] [B(12,24),B(2,4)]

If the ground state of the superblock is even under reflection symmetry, left and right side can be interchanged. That means that the reduced left blocks stored during the first half of a sweep (when the size of the left block is smaller than L/2) can be already used as right blocks in the second half of the same sweep. However, if the ground state of the superblock does not show reflection symmetry, the reduced blocks from the left and right side have to be constructed and stored independently. Therefore, in the absence of reflection symmetry stored space and CPU-time are doubled.

The finite system algorithm is terminated when convergence is reached, i.e. when the energy in succeeding sweeps does not improve (decrease) any more. It happens, however, that the energy stays on a certain level for two or three sweeps only to further decrease afterwards. Therefore, one cannot just compare the energies of the last two sweeps performed. The convergence behaviour of the model should be taken in account. This behaviour can be investigated by looking at small systems, where calculations are inexpensive. The number of sweeps necessary for convergence depends strongly on the system size L, the number of states kept m, and the model itself. Fermionic models need more sweeps than spin models, specially when they are doped. Typically, spin models converge in less than 10 sweeps, even for fairly large chains (L = 100 or larger), while calculations for the doped t - J model with the same m usually need roughly twice as many sweeps, even for smaller system sizes.

The first sweeps do not yield very accurate results. Their purpose is to generate a good set of blocks of different sizes. Therefore, the first sweeps are normally done with a small number of states. When convergence is approached the number of states kept in the truncation can be increased in order to improve accuracy. This helps saving CPU-time, specially when L is large.

5.5 Measurement of observables

The ground state energy of the superblock is determined every time the superblock is diagonalized. The value is used to determine whether convergence is reached. It turns out that in the finite system algorithm the energy is lowest when the two blocks forming the superblock have the same size. Therefore, this symmetric configuration is used to measure also all other observables in which one is interested.

Unfortunately the values for the other observables, such as the S^z value of a certain spin or the spin correlation between spins on different sites, are not as easily obtainable as the energy.

This is caused by the change of basis that is performed in every step. Even if we start out with a base where the demanded properties of the basis states are known, or could easily be calculated, this knowledge fades fast with the repeated linear combination of basis states of new base systems. Of course we could keep track, for instance, of áS^zñ for each site in each state but the computational effort would be enormous.

The way that is chosen is to carry out the measurement by actually evaluating the operator in the ground state. The expectation value of S^z on site i is calculated as

This is the textbook formula, but the difficulty of applying it here is not visible from this equation. The problem is the basis. The ground state is expanded in a basis that evolved in every step of the algorithm and could not have been anticipated at the beginning.

When the site i is added to the block, the representation of the S^z-operator in the basis of the enlarged block is known (in Eq. (5.11) the calculation is done for ). But in general this it too early to measure á ñ, because the symmetric configuration was not yet reached. In order to still have the right representation for in the symmetric configuration, the matrix has to be updated and stored every time the basis changes. Updating means that the basis change has to be performed on . If () denotes the representation of the S^z-operator on site i in the Hilbert space of the enlarged block with j sites (j ³ i) and O_j is the matrix that transforms and cuts the basis before adding site j+1, the representation of after the truncation is

When another site, the (j+1)th, is added, the representation of the operator also has to be adjusted according to

and truncated with O_j+1. Following this procedure we always have a representation of the operator in the current basis.

When it is time to perform measurements, i.e. when both blocks forming the superblock have the same size, we only have to find the representation of the operator in the Hilbert space of the superblock. If, e.g., the operator is acting on a site inside the left block this means tensorizing it with the unit element acting on all remaining spaces, namely, the two central sites and right block.

If i is small, i.e., the site is close to the end of the block, a lot of basis changes have to be performed before the measurement is carried out. Due to the truncations that go with this procedure, the accuracy is decreased. We expect a greater accuracy from observables on sites close to the middle of the chain.

The issue is somewhat more complicated for nonlocal operators, e.g. spin correlations like C_s(i,j) = áñ. In general one could just take the representations of the involved S^z-operators and multiply them, when the symmetric configuration is reached. However, there is a more accurate way to proceed. As an example, we consider a spin-spin correlation where j = i + 1 and the symmetric configuration is reached at L/2 = i + 2. The representation of the two operators at the chain size i+2 are

Therefore, one gets for the spin correlation

Another way of calculating the matrix is to multiply the two operators as soon as possible. In this case the operation can be done when the enlarged block size is i + 1, which gives

Then, from now on, (C_s(i,i+1))_i+1 is tranformed as a whole. Its representation for the enlarged block with i+2 sites is

Comparing Eq. (5.26) with Eq. (5.28) the difference is only a a

₊₁O_i+1-factor between the two S^z-operators. Without the truncation of states it would have been the product of two unitary matrices and thus be a unit matrix. The two ways of calculating (C_s(i,i+1))₊₂ would be equivalent. With the truncation, however, this is no longer true, as one can immediately see calculating O^†O with the projector in Eq. (5.19). The factor ₊₁O_i+1 leads to a loss of accuracy in the matrix multiplication since instead of multiplying the matrices, only their projections are multiplied. This error becomes worse the further the sites i and j in C_s(i,j) are apart, because every separating site adds another pair O^†O.

For DMRG procedure this means that we have to use the latter approach, multiplying the operators as soon as possible. Prior to the calculation, a list of the observables that we are interested in measuring has to be made. When growing the chain, the on-site operators are stored as soon as they are generated and updated every time the basis is changed. The products of two-site operators are formed and stored as soon as on has a representation for both operators, then they are also updated. In the case of operators that involve more single site operators, e.g. pairing operators, we have to proceed in the same way.

If we are measuring correlations between sites located on different blocks we can not multiply them before the symmetric configuration is reached. This means that measurements of correlations across the center will always have larger error than correlations where both sites are located in the same block.

From these explanations it has hopefully become clear that measuring observables involve additional storage and calculations. In order to save computer time the measurement process is started as late as possible in the calculations, after convergence has been reached. In the infinite size algorithm we restrict ourselves to the sites close to the middle of the chain, which were the last to be added. In the finite size algorithm measurements are not performed in the initial sweeps when the energy has not converged yet.

It is generally difficult to estimate the error on the values for the observables other than the energy for which it has been established that the error is proportional to the truncation error as discussed before. Unfortunately there is no known method to calculate the error of the observables from any other measured quantity. However, by checking how stable the results are as we change the number m of states kept in the truncation, we can have a qualitative control of the error. The error on the energy is generally smaller than that of the observables because it is a quantity that averages over the sum of many terms.

5.6 General remarks

When working with fermionic models such as the Hubbard and t - J models we have to implement the anticommutation of the fermionic operators Eqs. (2.2) and (2.3)[30].

In order to implement PBC in DMRG a specific superblock configuration is required[3]. DMRG performs poorly under PBC. Typically, if a given accuracy is obtained under FBC with m states kept in the truncation, then m² states will be needed to achieve the same accuracy under PBC[3].

An important improvement to DMRG involves keeping track of the wave function from step to step and perform a transformation into the basis corresponding to the current superblock. Since a good initial guess speeds up the Lanczos or Davidson convergence, this saves time in the diagonalization of the superblock[31].

When DMRG procedure converges to a fixed point the superblock ground state can be simply written as a matrix-product form and also be rederived through a simple variational ansatz making no reference to the DMRG construction. These very interesting analytical results are obtained in Ref. [32] and give some insight into the mechanisms working behind DMRG algorithms.

6 Conclusions and acknowledgments

The field of numerical simulations can stand by itself as a third way of doing science and its interaction with experiment and theory is very fruitful. We believe that scientific research in nearly every area in physics can benefit from it. As the power of nowadays computers increases rapidly the relatively new field of numerical simulations gains more and more prominence. To keep up with the advances in the hardware new methods and algorithms are being developed and traditional ones are being improved. The two techniques examined here are typical examples of such methods and algorithms.

The author is thankful to those (too many to name here!) who carefully read this manuscript helping to improve it in many ways. The author also acknowledges hospitality at the Instituto de Física de São Carlos - USP and the financial support from Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP and Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq - Brasil.

References

[1] S. R. White and R. M. Noack, Phys. Rev. Lett. 68, 3486 (1992).

[2] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).

[3] S. R. White, Phys. Rev. B 48, 10345 (1993).

[4] P. W. Anderson, Phys. Rev. 115, 2 (1959).

[5] J. Hubbard, Proc. Roy. Soc. A276, 238 (1963); A277, 237 (1964); A281, 401 (1964).

[6] For recent review on metal-insulator transitions see M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).

[7] P. W. Anderson, Science 235, 1196 (1987); for the derivation of the t -J model from the Hubbard model see, e.g., P. Fulde, Electron Correlations in Molecules and Solids, Vol. 100 of Solid-State Sciences (Springer, 1991).

[8] For a review see, e.g., M. N. Barber, Vol. 8 of Phase Transitions and Critical Phenomena (Academic Press, London, 1984).

[9] S-K Ma, Modern Theory of Critical Phenomena (Addison-Wesley Pub., 1997).

[10] M. Kikuchi and Y. Okabe, J. Phys. Soc. Japan 58, 679 (1989).

[11] C. Lanczos, J. Res. Nat. Bur. Stand. 45, 255 (1950).

[12] D. G. Pettifor e D. L. Weaire, The Recursion Method and Its Applications, Springer Series in Solid-State Sciences, Vol. 58 (Springer, Berlin/Heidelberg, 1985).

[13] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).

[14] For a more formal discussion of the Lanczos method see, e.g., H. Q. Lin and J. E. Gubernatis, Computers in Phys. 7, 400 (1993) and references therein.

[15] See, e.g., A. L. Malvezzi, M.A. Thesis, Univ. Federal de São Carlos - Programa de Pós-graduação em Física (1991).

[16] E. Dagotto and A. Moreo, Phys. Rev. D 31, 865 (1985).

[17] E. R. Gagliano, E. Dagotto, A. Moreo, and F. C. Alcaraz, Phys. Rev. B 34, 1677 (1986).

[18] E. R. Gagliano and C. A. Balseiro, Phys. Rev. Lett. 59, 2999 (1987).

[19] A. Jennings and J. J. Mckeown, Matrix Computation, (John Wiley & Sons/New York 1983).

[20] F. C. Alcaraz and R. R. P. Singh, Phys. Rev. B 47, 8298 (1993).

[21] J. Jaklic and P. Prelovøsek, Adv. in Phys. 49, 1 (2000).

[22] J. Riera and E. Dagotto, Phys. Rev. B 47, 15346 (1993); J. Pekariwicz and J. R. Shepard, Phys. Rev. B 56, 5366 (1997).

[23] E. Dagotto, G. B. Martins, J. Riera, A. L. Malvezzi, and C. Gazza, Phys. Rev. B 58, 12063 (1998).

[24] K. Wilson, Rev. Mod. Phys. 47, 773 (1975).

[25] J. W. Bray and S. T. Chui, Phys. Rev. B 19, 4876 (1979); C. Y. Pan and X. Chen, Phys. Rev. B 36, 8600 (1987); M. D. Kovarik, Phys. Rev. B 41, 6889 (1990); T. Xiang and G. A. Gehring, J. Magn. Magn. Mater. 104-107, 861 (1992).

[26] R. P. Feynman, Statistical Physics: A Set of Lectures (Benjamim, Reading, MA, 1972).

[27] M. Laukamp, Computational Study of ZN-doped Quantum Spin Chains and Ladders, PhD thesis (1998).

[28] O. Legeza and G. Fáth, Phys. Rev. B 53, 14349 (1996).

[29] E. Davidson, Computers in Physics 7, 519 (1993); J. Comput. Phys. 17, 87 (1975).

[30] S. Caprara and A. Rosengren, Nuclear Phys. B B493, 640 (1997); R. M. Noack, S. R. White, and D. J. Scalapino, in Computer Simulation Studies in Condensed-Matter Physics VII - Springer Proceedings in Physics, Vol. 78 (Edited by D. P. Landau, K. K. Mon, and H.-B Schüttler), 85 (1994).

[31] S. R. White, Phys. Rev. Lett. 77, 3633 (1996).

[32] S. Rommer and S. Östlund, Phys. Rev. Lett. 75, 3537 (1995); Phys. Rev. B 55, 2164 (1997).

Received on 27 August, 2002

1

It will be clear from the context whether by

S^z or

N we mean the operator or its eigenvalue.

2

In fact, a cyclic BC is enough to gain translational invariance.

3

As discussed below, this is a commodity we will not have when working with Lanczos or DMRG methods.

4

Here, the index

i is not labeling a site in a lattice but the states accessible to a given site. For instance,

D = 4 and 2 for the Hubbard and Heisenberg model, respectively.

5

It is possible to target several eigenstates simultaneously but, for a given computational effort, the accuracy decreases rapidly.

6

In this initial steps with no truncation the matrix

O has dimension (

m×

D)×(

m×

D) and Eq. (5.4) becomes a simple change of basis.

7

(D

E

₀)/2 converges to the ground state energy per site in the thermodynamic limit.

Publication Dates

• Publication in this collection
  23 Apr 2003
• Date of issue
  Mar 2003

History

• Received
  27 Aug 2002