Abstract
Quantum magnetism in low dimensions has been one of the central areas of theoretical research for many decades now. One of the key reasons for the long standing interest in this field has been the existence of simplified models, which serve as paradigms for understanding the role of strong interactions in manyelectron systems. Although simple, these models quite often can not be solved exactly. In this review, we discuss a variety of analytical and numerical methods, which treat the system in a systematic and controlled manner. The central method employed in all the studies is the density matrix renormalization group (DMRG) method. This is supported by small scale numerical exact calculations, and by analytical methods using field theoretic techniques. We have considered a number of magnetic systems, from magnetic clusters to extended lattices, and have found some novel quantum ground states and lowenergy elementary excitations. In some cases, we have also employed the finitetemperature DMRG method to accurately compute the lowtemperature thermodynamic properties such as specific heat and magnetic susceptibility.
Exact and Approximate Theoretical Techniques for Quantum Magnetism in Low Dimensions
Swapan K. Pati, S. Ramasesha, Diptiman Sen
Jawaharlal Nehru Centre for Advanced Scientific Research,
Jakkur, Bangalore 560064, India
Solid State and Structural Chemistry Unit,
Indian Institute of Science, Bangalore 560012, India
Centre for Theoretical Studies,
Indian Institute of Science, Bangalore 560012, India
1 Introduction
The nonrelativistic Schrödinger equation of a system of electrons is spin independent. It therefore appears at first glance that the solutions of the Schrödinger equation should also be spin independent. However, the indistinguishability of the electrons forces the total wave function which is a product of the spin wave function and the spatial wave function to be antisymmetric. This in turn implies that for two electrons, a spatially symmetric wave function should be associated with an antisymmetric spin wave function and vice versa. The different charge distribution in the spatially symmetric and antisymmetric wave functions leads to different coulomb repulsions by virtue of which the spin states which are symmetric and antisymmetric have different energies [1]. Dirac represented the splitting between the energies of the two spin states by the spin operator , where is the exchange integral involving the two spatial orbitals in which the two electrons are singly occupied. In most open shell atomic systems the exchange integral is large enough to force the total spin of the ground state configuration to be the largest permissible value. This in essence is the Hund’s rule of maximum multiplicity and is also the reason why we find transition and rare earth metal ions in the high spin states in nature.
In solids containing transition metal or rare earth ions surrounded by ligands, the relative alignment of the unpaired spins at the metal site is not at all obvious. In order to understand this, we should examine the possible pathways for the delocalization of the valence electrons in the system. If the favorable delocalization pathways involve antiparallel alignment of the metal ion spins then the nature of the exchange interaction between the metal ion spins is antiferromagnetic and ferromagnetic otherwise. This follows from the fact that delocalization of electrons lowers kinetic energy of the electrons and therefore the ground state corresponds to an alignment of spins that allows for maximum delocalization. This is indeed the reason why the ground state of a hydrogen molecule is a spin singlet. In a system with degenerate partially occupied orbitals, the Hund’s coupling favors highspin alignment of the electrons on an ion. If delocalization pathways exist that allow for these highspin states in the process of delocalization, then the alignment of the spins on two neighboring centers will be ferromagnetic. If delocalization pathways exist only when these highspin states are aligned antiparallel, one would have an antiferromagnetic alignment of the spin. Thus, the overall nature of the spin alignment is governed by a competition between the Hund’s coupling and electron delocalization [2, 3].
Given a collection of spins, the exchange Hamiltonian for the system is written as
(1) 
where is the effective exchange integral for the interaction between the spins at sites ’i’ and ’j’. Since the spins in the cluster arise from unpaired electrons of a transition metal ion in a crystal field, it is natural to expect that the spinorbit as well as spinspin interactions of the electrons in the ion could alter the nature of the total spin on the ion by giving the net spin a preferred direction of orientation. Such a situation can be easily handled by treating each as a vector and generalizing the exchange Hamiltonian as
(2) 
Such a model is often referred to as the XYZ spin model [2]. Two extreme cases are often studied, (i) the spin is assumed to have no projection on the XY plane, in which case the resulting model is the Ising model and corresponds to scalar spins, and (ii) the spin is assumed to have no projection on the zaxis, in which case we have an XY model or a planar spin model. The Ising model is a discrete classical model as it consists of no noncommuting operators in its Hamiltonian while the XY model could be classical or quantum mechanical. Usually, while dealing with large site spin systems, it is not uncommon to assume that the spins are classical, in the spirit of Bohr’s correspondence principle.
In the crystalline state, the spins in the solid would be arranged on a lattice. If the exchange interaction is predominant between spins along a single crystalline direction, the model could be treated as a one dimensional array of spins. There are many examples of solids in which this holds [4]. Likewise, it is also possible that the interactions amongst spins is large along two crystallographic directions and weak along a third direction and this would result in twodimensional spin system [5].
In this review article, we will mainly concern ourselves with the study of isotropic spin clusters and onedimensional spin systems, sometimes in the presence of an external magnetic field. We will be mostly interested in properties of the ground state and the lowlying excitations, since these are the states which govern the lowtemperature properties of systems such as the specific heat and magnetic susceptibility. The symmetries of a system often enable us to characterize the energy eigenstates in terms of quantum numbers such as the total spin , the component of the total spin along some particular direction, say, , spin parity (which is a symmetry for states with ), the wave number for a translation invariant system, and possibly other spatial symmetries depending on the structure. We will discuss below how the use of symmetries can help to lessen the numerical effort required to study the lowenergy states.
In the next two sections we introduce some exact numerical methods, and describe the application of these methods to magnetic clusters. In section 4, we discuss two analytical methods which use field theoretic approximations. In section 5, we describe an innovative way of solving the spin Hamiltonians, which goes beyond the conventional techniques and is based on the density matrix renormalization group (DMRG) theory. Various applications of DMRG to the properties of lowdimensional extended chains are described in sections 6 to 8.
2 Exact Calculations
The properties of a spin Hamiltonian can be computed from the eigenstates of the Hamiltonian which are in turn obtained by setting up the Hamiltonian matrix in a suitable basis and diagonalizing it thereafter. While the procedure itself is quite straightforward, the space spanned by the Hamiltonian rapidly increases with the number of the spins in the system. The Fock space dimensionality of a system of spins with spin is given by
(3) 
The Hamiltonian matrix is blockdiagonal in structure with each block corresponding to specified values of the quantities conserved by the Hamiltonian. Thus, for an isotropic spin system, the component of the total spin, , and the total spin are conserved. Restricting the Fock space to specified values of and gives Hilbert spaces whose dimensionalities are smaller than the Fock space dimensionality.
While constructing spin basis functions which are eigenstates of the total operator is quite simple, construction of spin adapted functions (SAF, eigenstates of the total operator) is not direct. Perhaps the simplest and chemically appealing way of constructing SAF is by the valence bond (VB) method which uses the RumerPauling rules. This method is best illustrated by applying it to a system of spins, each possessing a spin of half, in the total spin sector. A total spin singlet can be formed by choosing pairs of sites and spin coupling each of them to obtain a singlet. The product of these singlet pairs will be a spin eigenstate of the operator . This is illustrated in Fig. 1. However, there are more ways to spin couple in pairs than the number of linearly independent singlet states, e.g., the state in Fig. 1 can be expressed as a linear combination of the states and . The overcompleteness can be avoided by resorting to the RumerPauling rules. To implement this rule, we arrange the spins at the vertices of a regular gon. We draw lines between pairs of sites that are singlet paired. According to the RumerPauling rule, the subset of these, encompassing all diagrams (to be called ’legal’ diagrams) with no crossing lines, forms a complete and linearly independent set of states [6].
The RumerPauling rules can be easily extended to construct complete and linearly independent basis sets in higher spin Hilbert spaces involving spin objects. This is done most easily with the help of phantom sites. If we wish to construct VB diagrams for total spin subspace involving spin objects, then we introduce additional sites to be called phantom sites. Besides imposing the RumerPauling rules on the diagrams with sites, we impose the additional constraint that there should be no singlet lines amongst the phantom sites. In Fig. 2, we show a few examples of VB diagrams with higher total spin.
It is also quite simple to extend the VB rules to spin clusters made up of different site spins [7]. If the spin at a site is , then we replace this site by a set of sites, each with spin. We then proceed with constructing the VB basis, as though the system is made up entirely of spin objects, with one difference, namely, we impose the additional constraint that there should be no singlet lines within the subset of sites which replace the spin at site . The VB diagrams with total spin are constructed as before with the help of phantom sites. An example of a legal VB diagram involving higher site spins is shown in Fig. 2.
Generating and storing the VB diagrams on a computer is also quite simple. We associate one bit with every site. The state of the bit is ’1’, if in the VB diagram a line begins at the site and the state is ’0’ if a line ends at the site. Thus, we can associate an integer of bits with every VB diagram involving spin objects. This association is unique if we decipher the bit pattern of the integer corresponding to the diagram from insideout, much like expanding an algebraic expression with multiple parentheses. In Fig. 2, we have also shown the bit pattern and the associated integer for each VB diagram. The VB diagrams are generated on a computer by checking the bit pattern in all bit integers to see if they satisfy the criterion for representing the desired VB diagram. This also allows us to generate the VB diagrams as an ordered sequence of the integers that represent them, a fact that helps in rapid generation of the Hamiltonian matrix.
The Hamiltonian matrix in the VB basis can be easily constructed by knowing the action of the operator for spin particles (i) on a singlet line joining sites and (ii) on the pair of sites singlet paired to two different sites (Fig. 3). In the case of Hilbert spaces with nonzero total , the Hamiltonian involves spin exchange between the real sites only. However, these exchange operators could lead to VB diagrams in which the phantom sites are interconnected. In this event, simply neglecting these resultant states is sufficient to ensure that we are dealing exactly with the spin Hilbert space. The Hamiltonian for spin clusters with arbitrary spins can be treated as consisting of operators with only spin objects. This is done by replacing the spin exchange operator between sites and , by the operator , where the operators and are the usual spin operators.
The matrix representing the Hamiltonian in the VB basis is in general nonsymmetric since the VB basis is nonorthogonal. However, the matrix itself is sparse. There exist efficient numerical algorithms [8] for obtaining the lowlying eigenstates of sparse nonsymmetric matrices, and it is possible to solve a nonsymmetric matrix eigenvalue problem for a million by million matrix with about a 100 million nonzero matrix elements on a powerful PC based workstation.
While VB theory guarantees spin purity in the computed eigenstates of the spin conserving Hamiltonians, it has several drawbacks. Computation of quantities such as spinspin correlation functions and spin densities are not easy because a sitespin operator operating on a VB diagram spoils the spin purity of the diagram. This could be overcome by converting the eigenstate of the Hamiltonian in the VB basis to the constant basis. Another difficulty with the VB procedure is exploiting the spatial symmetries of the problem. Operation by a spatial symmetry operator on a legal VB diagram could lead to illegal VB diagrams. Disentangling these illegal diagrams into legal diagrams can be computationally prohibitive [9].
In more general spin problems, it is often advantageous to use the constant basis and exploit all the spatial symmetries. Partial spin symmetry adaptation in these cases is also possible by using the spin parity operator. The effect of the spin parity operator on a basis state is to flip all the spins in the state. In the sector, it is possible to factor the Hilbert space into odd and even parity Hilbert spaces. The odd (even) parity Hilbert space is spanned by basis vectors with odd (even) total spin. This also has the effect of reducing the dimensionality of the Hilbert space besides providing partial spin symmetry adaptation. It is rather simple to set up the Hamiltonian matrix in the symmetry adapted basis. The Hamiltonian matrix is symmetric and usually very sparse. The lowest few eigenstates can be easily computed by employing the Davidson algorithm. Given these eigenstates, the computation of properties can proceed by converting an eigenstate in the symmetrized basis into that in the unsymmetrized basis. The orthogonality of the basis states as well as the simple rules involved in obtaining the resultant when a basis state is operated upon by any type of spin operators in any combination affords easy computation of a variety of properties of a magnetic system.
The exact diagonalization techniques discussed above are in general applicable to systems whose Hilbert space dimensionality is about 10 million. The major problem with exact diagonalization methods is the exponential increase in dimensionality of the Hilbert space with the increase in the system size. Thus, the study of larger systems becomes not only CPU intensive but also memory intensive as the number of nonzero elements of the matrix increases rapidly with system size. With increasing power of the computers, slightly larger problems have been solved every few years. To illustrate this trend, we consider the case of the spin Heisenberg chain. In 1973, ten years before the Haldane conjecture, De Neef [10] used the exact diagonalization procedure to solve a site spin chain. In 1977, Blote [11] diagonalized the Hamiltonian of a chain of 10 sites. In 1982, Botet and Jullien [12] increased this to 12 sites. In 1984, Parkinson and Bonner [13] solved the site spin problem. In the same year, Moreo [14] solved the site spin chain. In 1990, Takahashi [15] pushed this up to sites. And in 1994, Golinelli et al. [16] have solved for the lowlying states of a site spin chain. The growth in chain length of the longest spin chain solved is almost linear with time, roughly increasing by sites in every three years. Just to remind ourselves, the Fock space dimensionality in this case increases as with chain length N. The size of the matrix also increases similarly and the CPU and storage scales quadratically with the size of the matrix, if we are targeting only a few eigenstates. For this reason, for systems which span much larger spaces, the focus has shifted to approximate techniques.
3 Applications to Spin Clusters
In recent years, some of the magnetic clusters that have been studied extensively are the Mn [17], Fe [18] and V [19] clusters. These clusters show many interesting phenomena such as quantum resonant tunneling and quantum interference [20]. Basic to a proper understanding of these phenomena is a knowledge of the lowenergy excitation spectrum in these systems. The methods discussed under exact diagonalization schemes allow us to calculate the lowenergy excitation spectrum, given a set of exchange constants. However, the exchange constants themselves are not known with any certainty. Therefore, it is all the more important to be able to carry out exact diagonalization studies of lowlying states to infer the possible sign and magnitude of the exchange constants [7].
For the Mn cluster, we show the geometry and the exchange parameters in Fig. 4. The crystal structure suggests that the exchange constant is largest and antiferromagnetic in nature [21]. Based on magnetic measurements, it has been suggested that has a magnitude of 215K. The magnitude and sign of the other exchange constants are based on comparisons with manganese systems in smaller clusters. It has been suggested that the exchange constant and are antiferromagnetic and have a magnitude of about 85K. However, for the exchange constant , there is no concrete estimate, either of the sign or of the magnitude. In an earlier study, the Mn  Mn pair with the strongest antiferromagnetic exchange constant was replaced by a composite spin object [22] and the exchange Hamiltonian of the cluster solved for three different sets of parameters. It was found that the ordering of the energy levels were very sensitive to the relative strengths of the exchange constants. In these studies, was set to zero and the lowlying excited states were computed. Besides, only states with spin up to ten could be obtained because of the replacement of the higher spin ion pairs by the composite spin objects.
The technique described earlier, however, allows an exact computation of the lowlying states of Mn. The results of the exact calculations are presented in Table 1. We note that none of the three sets of parameters studied using an effective Hamiltonian, gives the correct ground and excited states, when an exact calculation is performed. It appears that setting the exchange constant to zero, cannot yield an ground state (Table 1, cases A, B and C). When is equal to or slightly larger than (cases A and B, Table 1), we find a singlet ground state, unlike the result of the effective Hamiltonian in which the ground state has and respectively. The ground state has spin , when is slightly smaller than (case C, Table 1). In all these case, the first few lowlying states are found to lie within 20K of the ground state.
When we use the parameters suggested by Chudnovsky [17] (case D, Table 1), we obtain an ground state separated from an first excited state by 223K. This is followed by another excited state at 421K. Only when the exchange constant is sufficiently strongly ferromagnetic (case E, Table 1), do we find an ground state with an excited state separated from it by a gap of 35K, which is close to the experimental value [23]. The second higher excited state has , and is separated from the ground state by 62K.
In Fig. 5, we show the spin density [24] for the Mn cluster in the ground state for the state. While the manganese ions connected by the strong antiferromagnetic exchange show opposite spin densities, it is worth noting that the total spin density on these two ions is , well away from a value of expected if these ions were indeed to form a spin object.
The Fe cluster is shown in Fig. 6. Each of the Fe ions has a spin of and the ground state of the system has a total spin , with excited state separated from it by about 20K. All the exchange interactions in this system are expected to be antiferromagnetic. While the structure of the complex dictates that the exchange interaction along the back of the butterfly should be small in comparison with the interaction across the wing [25], in earlier studies it was reported that such a choice of interaction parameters would not provide a ground state [26].
Results from the exact calculations of the eigenstates of the Fe cluster using three sets of parameters is shown in Table 2. In two of these cases, is very much smaller than . We find that in all these cases, the ground state has a spin and the lowest excited state has spin, . One of the main difference we find amongst the three sets of parameters is in the energy gap to the lowest excited state (Table 2). For the set of parameters used in the earlier study, this gap is the lowest at 3.4K. For the parameter sets 1 and 3 [27], this gap is respectively 13.1K and 39.6K. While in cases 1 and 2, the second excited state is an state, in case 3, this state also has spin 9.
The spin densities in all the three cases for the ground state are shown in Fig. 7. The spin densities in all cases are positive at the corners. In cases 1 and 2, the spin density on the Fe ion on the backbone is positive and negative on the remaining two Fe sites [28]. However, in case 3, the negative and positive spin density sites for Fe ions in the middle of the edges is interchanged. This is perhaps due to the fact that in cases 1 and 2, the exchange constant is less than , while in case 3, this is reversed. Thus, a spin density measurement can provide relative strengths of these two exchange constants. In all the three case, the difference between the spin densities in the ground and excited states is that the decrease in the spin density in the excited state is mainly confined to the corner Fe sites.
4 Field Theoretic Studies of Spin Chains
Onedimensional and quasionedimensional quantum spin systems have been studied extensively in recent years for several reasons. Many such systems have been realized experimentally, and a variety of theoretical techniques, both analytical and numerical, are available to study the relevant models. Due to large quantum fluctuations in low dimensions, such systems often have unusual properties such as a gap between the ground state and the excited states. The most famous example of this is the Haldane gap which was predicted theoretically in integer spin Heisenberg antiferromagnetic chains [29], and then observed experimentally in a spin system [30]. Other examples include the spin ladder systems in which a small number of onedimensional spin chains interact amongst each other [31]. It has been observed that if the number of chains is even, i.e., if each rung of the ladder (which is the unit cell for the system) contains an even number of spin sites, then the system effectively behaves like an integer spin chain with a gap in the lowenergy spectrum. Some twochain ladders which show a gap are [32], [33] and [34]. Conversely, a threechain ladder which effectively behaves like a halfoddinteger spin chain and does not exhibit a gap is [33]. A related observation is that the quasionedimensional system spontaneously dimerizes below a spinPeierls transition temperature [35]; then the unit cell contains two spin sites and the system is gapped.
The results for gaps quoted above are all in the absence of an external magnetic field. The situation becomes more interesting in the presence of a magnetic field [36]. Then it is possible for an integer spin chain to be gapless and a halfoddinteger spin chain to show a gap above the ground state for appropriate values of the field [3745]. This has been demonstrated in several models using a variety of methods such as exact diagonalization of small systems and bosonization [46, 47]. In particular, it has been shown that the magnetization of the system can exhibit plateaus at certain nonzero values for some finite ranges of the magnetic field. Further, for a Hamiltonian which is invariant under translation by one unit cell, the value of the magnetization per unit cell is quantized to be a rational number at each plateau [37]. In section 8, we will study the magnetization plateau which can occur in a threechain ladder.
In the next two subsections, we will discuss some field theoretic methods which can be used for studying spin chains and ladders. These methods rely on the idea that the lowenergy and longwavelength modes of a system (i.e., wavelengths much longer than the lattice spacing if the system is defined on a lattice at the microscopic level) can often be described by a continuum field theory.
4.1 Nonlinear model
The nonlinear model (NLSM) analysis of antiferromagnetic spin chains with the inclusion of (nextnearest neighbor coupling) and (dimerization) proceeds as follows [48]. The Hamiltonian for the frustrated and dimerized spin chain can be written as
(4) 
The interactions are schematically shown in Fig. 8. The region of interest is defined by and . We first do a classical analysis in the limit to find the ground state configuration of the spins. Let us make the general ansatz that the ground state is a coplanar configuration of the spins with the energy per spin being equal to
(5) 
where is the angle between the spins and and is the angle between the spins and . Minimization of the classical energy with respect to the yields the following three phases.
(i) Neel phase: This phase has ; hence all the spins point along the same line and they go as along the chain. This phase is stable for .
(ii) Spiral phase: Here, the angles and are given by
(6) 
where and . Thus the spins lie on a plane. This phase is stable for .
(iii) Colinear phase: This phase (which needs both dimerization and frustration) is defined to have and ; hence all the spins again point along the same line and they go as along the chain. It is stable for .
These phases along with their boundaries are depicted in Fig. 9. Thus even in the classical limit , the system has a rich ground state ‘phase diagram’ [49].
We can now go to the next order in , and study the spin wave spectrum about the ground state in each of the phases. The main results are as follows. In the Neel phase, we find two zero modes, i.e., modes for which the energy vanishes linearly at certain values of the momentum , with the slope at those points (the velocity) being the same for the two modes. In the spiral phase, we have three zero modes, two with the same velocity describing outofplane fluctuations, and one with a higher velocity describing inplane fluctuations. In the colinear phase, we get two zero modes with equal velocities just as in the Neel phase. The three phases also differ in the behavior of the spinspin correlation function in the classical limit. is peaked at , i.e., at in the Neel phase, at in the spiral phase and at in the colinear phase.
To study the interactions between the spin waves, it is convenient to derive a semiclassical NLSM field theory which can describe the lowenergy and longwavelength excitations. The field theory in the Neel phase is given by a NLSM with a topological term [29, 47]. The field variable is a unit vector with the Lagrangian density
(7) 
where is the spin wave velocity, is the coupling constant (which describes the strength of the interactions between the spin waves), and is the coefficient of the topological term (the integral of this term is an integer which defines the winding number of a field configuration ). Note that is independent of in the NLSM. (Time and space derivatives are denoted by a dot and a prime respectively). For and less than a critical value, it is known that the system is gapless [47, 50]. For any other value of , the system is gapped. For , one therefore expects that integer spin chains should have a gap while halfoddinteger spin chains should be gapless. This is known to be true even for small values of like (analytically) and (numerically) although the field theory is only derived for large . In the presence of dimerization, one expects a gapless system at certain special values of . For , the special value is predicted to be . We see that the existence of a gapless point is correctly predicted by the NLSM. However, as we will see later, according to reliable numerical results from DMRG, is for [51] and it decreases with as shown in Fig. 10. These deviations from field theory are probably due to higher order corrections in which have not been studied analytically so far.
In the spiral phase, it is necessary to use a different NLSM which is known for [52, 53]. The field variable is now an matrix . The Lagrangian density is
(8) 
where , with , and and are diagonal matrices with diagonal elements and respectively. Note that there is no topological term. Hence there is no apparent difference between integer and halfoddinteger spin chains in the spiral phase. A oneloop renormalization group [52] and large analysis [53] indicate that the system should have a gap for all values of and , and that there is no reason for a particularly small gap at any special value of .
Finally, in the colinear phase, the NLSM is known for , i.e., for the spin ladder. The Lagrangian is the same as in (7), with , and . There is no topological term for any value of , and the model is therefore gapped.
The field theories for general in both the spiral and colinear phases are still not known. Although the results are qualitatively expected to be similar to the case in the spiral phase and the case in the colinear phase, quantitative features such as the dependence of the gap on the coupling strengths will require the explicit form of the field theory.
4.2 Bosonization
Another field theoretic method for studying spin systems in one dimension is the technique of bosonization [46, 47, 54, 55, 56]. This technique consists of mapping bosonic operators into fermionic ones, and then using whichever set of operators is easier to compute with. For instance, consider a model with a single species of fermion with a linear dispersion relation , where the denotes the right and leftmoving fermions respectively (with the corresponding fields being denoted by and ), and denotes the velocity. Similarly, consider a model with a single species of boson with the dispersion relation ; the right and leftmoving fields are denoted by and respectively. Then it can be shown that these operators are related to each other as
(9) 
The length parameter is a cutoff which is required to ensure that the contribution from highmomentum modes do not produce divergences when computing correlation functions. It is convenient to define the two bosonic fields
(10) 
Then the fermionic density is given by
(11) 
where is the background density; fluctuations around this density are described by the quantum fields or .
Although the dispersion relation is generally not linear for all the modes of a given system, it often happens that the lowenergy and longwavelength modes can be studied using bosonization. For a fermionic system in one dimension, these modes are usually the ones lying close to the two Fermi points with momenta respectively. One can define right and leftmoving fields and which vary slowly on the scale length ,
(12) 
Quantities such as the density will generally contain terms which vary slowly as well as terms varying rapidly on the scale of ,
(13)  
One can compute various correlation functions in the bosonic language. Consider an operator of the form
(14) 
We find the following result for the twopoint equaltime correlation function at spatial separations which are much larger than the microscopic lattice spacing ,
(15) 
where denotes an interaction parameter which will be described below. Note that the correlation function decays as a power law. In the language of the renormalization group, the scaling dimension of is given by .
We can now study a quantum spin chain using bosonization. To be specific, let us consider a spin chain described by the anisotropic Hamiltonian
(16) 
where the interactions are only between nearest neighbor spins, and . and are the spin raising and lowering operators, and denotes a magnetic field. Note that the model has a invariance, namely, rotations about the axis. When and , the invariance is enhanced to an invariance, because at this point the model can be written simply as . Although the model in (16) can be exactly solved using the Bethe ansatz, and one has the explicit result that the model is gapless for a certain range of values of and (see Ref. [39]), it is not easy to compute explicit correlation functions in that approach. We therefore use bosonization to study this model.
We first use the JordanWigner transformation to map the spin model to a model of spinless fermions. We map an spin or a spin at any site to the presence or absence of a fermion at that site. We introduce a fermion annihilation operator at each site, and write the spin at the site as
(17) 
where the sum runs from one boundary of the chain up to the site (we assume open boundary conditions here for convenience), or is the fermion occupation number at site , and the expression for is obtained by taking the hermitian conjugate of , The string factor in the definition of is added in order to ensure the correct statistics for different sites; the fermion operators at different sites anticommute, whereas the spin operators commute.
We now find that
(18) 
We see that the spinflip operators lead to hopping terms in the fermion Hamiltonian, whereas the interaction term leads to an interaction between fermions on adjacent sites.
Let us first consider the noninteracting case given by . By Fourier transforming the fermions, , where is the lattice spacing and the momentum lies in the first Brillouin zone , we find that the Hamiltonian is given by
(19) 
where
(20) 
The noninteracting ground state is the one in which all the singleparticle states with are occupied, and all the states with are empty. If we set the magnetic field , the magnetization per site will be zero in the ground state; equivalently, in the fermionic language, the ground state is precisely halffilled. Thus, for , the Fermi points () lie at . Let us now add the magnetic field term. In the fermionic language, this is equivalent to adding a chemical potential term (which couples to or ). In that case, the ground state no longer has and the fermion model is no longer halffilled. The Fermi points are then given by , where
(21) 
It turns out that this relation between (which governs the oscillations in the correlation functions as discussed below) and the magnetization continues to hold even if we turn on the interaction , although the simple picture of the ground state (with states filled below some energy and empty above some energy) no longer holds in that case.
In the linearized approximation, the modes near the two Fermi points have the velocities , where is some function of , and . Next, we introduce the slowly varying fermionic fields and as indicated above; these are functions of a coordinate which must be an integer multiple of . Finally, we bosonize these fields. The spin fields can be written in terms of either the fermionic or the bosonic fields. For instance, is given by the fermion density as in Eq. (17) which then has a bosonized form given in Eq. (13). Similarly,
(22)  
where since is an integer. This can now be written entirely in the bosonic language; the term in the exponential is given by
(23) 
where we have ignored the contribution from the lower limit at .
We can now use these bosonic expressions to compute the twospin, equaltime correlation functions . We find that
(24) 
where are some constants. The parameters and are functions of and ; the exact dependence can be found from the web site given in Ref. [39]. For , is given by the analytical expression
(25) 
Note that at the invariant point and , we have , and the two correlations and have the same forms.
In addition to providing a convenient way of computing correlation functions, bosonization also allows us to study the effects of small perturbations. For instance, a physically important perturbation is a dimerizing term
(26) 
where is the strength of the perturbation. Upon bosonizing, we find that the scaling dimension of this term is . Hence it is relevant if ; in that case, it produces an energy gap in the system which scales with as
(27) 
This kind of phenomenon occurs in spinPeierls systems such as ; below a transition temperature , they go into a dimerized phase which has a gap [57].
5 Density Matrix Renormalization Group Method
The method which held the promise of overcoming the difficulty of exploding dimensionalities is the renormalization group (RG) technique, in which one systematically throws out the degrees of freedom of a manybody system. While this technique found dramatic success in the Kondo problem [58], its straightforward extension to interacting lattice models was quite inaccurate [59].
In early 1992, the key problems associated with the failure of the old RG method were identified and a different renormalization procedure based on the eigenvalues of the manybody density matrix of proper subsystems was developed [60, 61]. This method has come to be known as the density matrix renormalization group (DMRG) method and has found dramatic success in solving quasionedimensional manybody Hamiltonians.
In a realspace RG approach, one begins by subdividing the total system into several blocks and proceeds to iteratively build effective blocks so that at each iteration, each effective block represents two or more blocks of the previous iteration, without increasing the Fock space dimensionality of the blocks from what existed at the previous iteration. Usually, one starts with each consisting of a single site. Since the Hilbert space grows exponentially with the increase in system size, one truncates the number of states kept at each iteration.
The main reason for the failure of the old RG methods is the choice of the states retained at each stage of the iteration [60]. White [61], recognized that the weakness of the old RG procedure was in the truncation of the Fock space of a block based on the eigenvalues of the block Hamiltonian being renormalized. He replaced this choice by introducing a completely different truncation scheme than the one that was used in the old quantum RG procedures. The choice is the eigenvalues of the reduced density matrix of the block constructed from the desired state of the full Hamiltonian. The truncated Fock space is now spanned by the eigenvectors of the reduced density matrix of order () corresponding to the highest eigenvalues of the reduced manybody density matrix. The reason for choosing the eigenvalues of the reduced density matrix as a criterion for implementing a cutoff is that, the larger the density matrix eigenvalue, the larger is the weight of the eigenstate of the density matrix in the expectation value of any property of the system. This result becomes evident when all the dynamical operators are expressed as matrices in the basis of the eigenvectors of the density matrix. The expectation value of any operator is simply
(28) 
where is the density matrix eigenvalue. The larger the value of a particular , the larger is its contribution to the expectation value, for a physically reasonable spread in the diagonal matrix elements .
The manybody density matrix of a part of the system can be easily constructed as follows. Let us begin with given state of , which is called the universe or superblock, consisting of the system (which we call a block) and its environment . Let us assume that the Fock space of and are known, and can be labeled as and respectively. The representation of in the product basis of and can be written as
(29) 
where we assume the coefficients to be real, without loss of generality. Then the reduced manybody density matrix for block , is defined as
(30) 
The eigenvalue of the density matrix gives the probability of finding the corresponding eigenstate in the projection of on block . It therefore follows that the eigenvectors with highest eigenvalues of the density matrix of , are the optimal or most probable states to be retained while the system is augmented.
In the early literature in quantum chemistry, the eigenvectors corresponding to large eigenvalues of oneparticle density matrices were employed as the orbital basis for carrying out a configuration interaction (CI) calculation. The eigenvectors of the density matrix were called the ‘natural’ orbitals, and it was observed that the CI procedure converged rapidly when the ‘natural’ orbitals were employed in setting up the Slater determinants [62].
The DMRG scheme differs from the ‘natural’ orbital scheme in two important respects: (i) the reduced density matrices are manybody density matrices, and (ii) the size of the system in terms of the number of sites being studied at each iteration is usually augmented by two sites. However, the Hamiltonian matrix that one encounters from iteration to iteration, remains roughly of the same order while the matrix elements keep changing (renormalized). In this sense, the procedure can be called a renormalization procedure. The coupling constants (the Hamiltonian matrix elements) keep changing while the system size increases, as in the RG procedure carried out within a blocking technique.
5.1 Implementation of the DMRG method
We now describe the procedure to carry out the computations. One starts the computation with a small size system, , which can be exactly solved, , depending on the degree of freedom at each site. By exact diagonalization, one gets the desired eigenstate of that system. The density matrices of the left and right blocks, each consisting of sites (in principle it is not necessary to have the same number of sites for the two blocks, although in practice this is what is most generally used), are obtained from the desired eigenstate. The density matrices are diagonalized and at the first iteration usually all the density matrix eigenvectors (DMEV) are retained. The Hamiltonian matrix of the left and right blocks (denoted by and ) obtained in any convenient basis are transformed into the density matrix eigenvector basis. So also are the matrices corresponding to the relevant site operators in both blocks. Now, the iterative procedure proceeds as follows.

Construct a superblock , consisting of the block , two additional sites and and the block . Thus, at the first iteration, the system has sites.

Set up the matrices for the total Hamiltonian of the superblock in the direct product basis of the DMEV of the blocks and and the Fock space states of the new sites. Considering that the new sites are spin sites with Fock states each, the order of the total Hamiltonian matrix will be , where is the dimension of the block DMEV basis.

Diagonalize the Hamiltonian of the superblock to find the desired eigenstate . Using the state , evaluate all properties of the superblock of interest.

Construct the reduced manybody density matrix, , for the new block . If the system does not possess reflection symmetry, construct the density matrix, , for the new right block as well.

Diagonalize the density matrix, , and if necessary . Usually, the density matrix is blockdiagonal in the component of the total spin of the block, and it becomes computationally efficient to exploit such quantum numbers. Construct a nonsquare matrix , with columns, each column being an eigenvector of the density matrix corresponding to one of the largest eigenvalues. The number of rows in the matrix corresponds to the order of the density matrix.

Construct the matrices corresponding to the Hamiltonian, , of the new left block , and the site spin operators ( and ) of all the necessary sites. The operators are simply the adjoints of the operators.

Renormalize all the matrices corresponding to the block and site operators by using the RG transformation matrix , e. g., . The resulting renormalized matrices are of order and the procedure amounts to a simultaneous change of basis and a truncation.

Replace the block by . If the system does not possess reflection symmetry, replace by .

Go to step 1.
Using the blockdiagonal nature of the density matrix, besides reducing the requirement in CPU time, also allows one to label the DMEV by the appropriate component of the total spin of the block (). The Fock space of the individual sites that are added at each iteration are eigenstates of the site spin and number operators. This allows us to target a definite projected spin () state of the total system.
We now briefly describe the mathematical notations that we have used so far for various states. A state of is given by the tensor product of a state of with quantum number and an index , and a state of the additional site. Thus,
(31) 
A state of a superblock is given by
(32) 
The eigenstate of the Hamiltonian of the superblock can be written as
(33) 
The density matrix for then will have a block structure and can be expressed as
(34) 
The above algorithm is called the infinite lattice DMRG algorithm because this procedure is best suited for the system in the thermodynamic limit, i.e., when the properties of the system are extrapolated to the infinite system size limit.
5.2 Finite size DMRG algorithm
If we are interested in accurate properties of the system at a required size, then it is possible to improve upon the accuracies obtainable from the infinite DMRG procedure. This involves recognizing that the reduced manybody density matrices at each iteration correspond to a different system size. For example, when we carry out the DMRG procedure to obtain the properties of a system of sites, at an iteration corresponding to sites (), the reduced density matrix we construct is that of a block of sites in a system of sites. However, if our interest is in the site system, we should employ the density matrix of the block of sites in a site system. It is possible to construct, iteratively, the site reduced density matrix of the site system. This is achieved by the so called finitesize algorithm [61]. This method provides highly accurate solutions even when the states of the full Hamiltonian have inhomogeneous (symmetry breaking) properties.
To obtain the site result, we should perform the infinite lattice algorithm up to sites first storing all operators in each iteration. Now the algorithm for finite lattices with reflection symmetry (left block = right block), proceeds as follows.

On reaching a system size of sites, obtain the density matrix of the block of sites.

Use the density matrix of sites on the left and that of sites on the right, add two new sites as in the infinite DMRG procedure and obtain the desired eigenstate of the system.

Now obtain the reduced density matrix of the sites from the eigenstate of the previous iteration obtained in the direct product basis of the DMEVs of the site, site density matrices, and the Fock space states of the individual sites.

Go back to step 2, replacing by and by and iterate until a single site results on the right and sites result on the left.

Since the system has reflection symmetry, use the density matrix of the sites on the right and construct the system as builtup from three individual sites on the left and sites on the right. Obtain the desired eigenstate of the system in this basis.

Now obtain the new site density matrix on the left and site density matrix on the right. Replace the singlesite on the left by two sites and sites on the right by sites in step 5.

Repeat steps 5 and 6 until sites are obtained both on the left and right. The properties of the system obtained from the eigenstates at this stage corresponds to the first iteration of the finitesize algorithm. We can now go back to step 1 and carry through the steps to obtain properties at later iterations of the finitesize DMRG algorithm.
In systems without reflection symmetry, the DMEVs of the right and left parts are not identical even if the sizes of the reduced systems are the same. The finiteDMRG algorithm in this case involves first constructing the density matrices of the left part for sizes greater than and on reaching the density matrix of sites, reducing the size of the leftpart and increasing that of the right, from one site to . This will result in the refined density matrices of both the right and the left block of the total system, for block sizes of . At this stage, we can compute all the properties and continue the reverse sweep until the right block is of size and the left block is of size . The forward sweep that follows will increase the block size on the left and decrease that of the right. We would have completed the second iteration when the two block sizes are equal. The forward and reverse sweeps can be continued until we reach the desired convergence in the properties of the whole system.
5.3 Calculation of properties in the DMRG basis
At the end of each iteration, one can calculate the properties of the targeted state [63]. The reduced manybody density matrix computed at each iteration can be used to calculate the static expectation values of any site operator or their products. Care should be taken to use the density matrices appropriate to the iteration. The expectation value of a site property corresponding to the operator can be written as
(35) 
is the density matrix of the block in which the site is situated and is the matrix of the renormalized site operator at site . For calculating correlation functions, one can use a similar equation. The correlation function between two site operators belonging to separate blocks can be written as
(36) 
However, the accuracy of this procedure turns out be very poor if the sites and belong to the same block [61]. The reason is that a feature implicit in the above procedure is the resolution of identity by expansion in terms of the complete basis. Unfortunately, the basis in which the site operators are represented is incomplete and such an expansion is therefore error prone. To circumvent this difficulty, it has been suggested [61] that, one obtains the matrix representation of the products of the site operators from the first occurrence of the product pair and, by renormalizing the product operator , at every subsequent iteration until the end of the RG procedure. Then, the correlation function between and (where and belong to the same block) can be evaluated as
(37) 
This procedure is found to be more accurate in most cases.
5.4 Remarks on the applications of DMRG
The DMRG method is currently the most accurate method for large quantum lattice models in one dimension. It can be applied to interacting bosonic, fermionic or spin models as well as to models which have interactions amongst them. The overall accuracy of the DMRG method is exceptionally high for onedimensional systems with only nearest neighbor interactions. For a spin chain where exact Bethe ansatz groundstate energy is available, the DMRG ground state energy per site in units of the exchange constant , is found to be accurate to seven decimal places with a cutoff [61]. The method is found to be almost as accurate for the onedimensional Hubbard model, where again it is possible to compare the DMRG results with exact results obtained from the exact analytic Bethe ansatz solution [64].
Since higher dimensionality is equivalent to longerrange interactions within one dimension, the model also restricts the range of interactions in one dimension. It has been noted that the number of DMEVs that should be retained in a calculation on higher dimensional systems, for accuracies comparable to accuracy in one dimension, scales exponentially with dimensionality. Thus, to obtain accuracy comparable to that obtained in a chain of sites for a cutoff , in a square lattice, the number of DMEVs needed to be retained for the corresponding 2dimensional lattice is .
Extending the range of interactions to nextnearest neighbors does not significantly deteriorate the accuracy [65]. However, inclusion of cyclic boundary conditions reduces the accuracy of the method significantly, although in onedimension, the DMRG method still would outperform any other method for the same system size. In the DMRG procedure, the most accurate quantity computed is the total energy. In dealing with other quantities such as correlation functions, caution must be exercised in interpreting the results.
The density matrix eigenvalues sum to unity and the truncation error, which is defined as the sum of the density matrix eigenvalues corresponding to the discarded DMEV, gives a qualitative estimate as to the accuracy of the calculation as well as providing a framework for the extrapolation to the limit. The accuracy of the results obtained in this way is unprecedented [66, 67]. The accuracy of the ground state energy per site for the spin chain is limited by the precession of machine arithmetic, i.e., . Similarly, the accuracy persists even while calculating for the Haldane gap, e.g., the gap is evaluated to be .
Another aspect of the DMRG technique worth noting is that the method is best suited for targeting one eigenstate at a time. However, it is possible to obtain reasonable results for a set of states by using an average manybody reduced density matrix constructed as a weighted sum of the density matrices corresponding to each of the states in question. One way of constructing the average density matrix is by using a statistical weight for the chosen set of states; the averaged density matrix in this instance is given by
(38) 
where with and being the Boltzmann constant and temperature respectively. One can thus extend the DMRG method to finite temperatures.
Finite size algorithms have been used extensively to study the edge states and systems with impurities, where substantial improvement of the accuracy is needed to characterize the various properties of a finite system. The DMRG method has been applied to diverse problems in magnetism: study of spin chains with [68], chains with dimerization and/or frustration [51, 65, 69, 70], coupled spin chains [65, 71, 72], to list a few. Highly accurate studies of the structure factor and string order parameter (topological long range order) [66] as well as edge states in Haldane phase systems [73] have been performed. Dynamical properties for both spin and fermionic systems with DMRG have also been reported within the Maximum Entropy Method [74] as well as the continued fraction [75] and correction vector [76] approaches. Also, DMRG has been successfully formulated to obtain lowtemperature thermodynamic properties for various spin systems [77, 78], and the solution of models of spin chains dynamically coupled to dispersionless phonons [79]. Additionally, Nishino and Okunishi have derived two reformulations of DMRG, namely, the product wave function renormalization group (PWFRG) [80], and the corner transfer matrix renormalization group (CTMRG) [81] methods. These methods offer the means of calculating dynamical correlation functions in spin chains as well as highly accurate results for the dimensional Ising model at criticality.
6 Frustrated and Dimerized Spin Chains
It is well known that the onedimensional XY chain can be mapped on to a onedimensional noninteracting spinless fermion model. The isotropic spin chain will then map on to a chain of interacting spinless fermions. According to Peierls’ theorem, a partly filled onedimensional band of noninteracting fermions is unstable with respect to a lattice distortion that results in an insulating ground state. It has been shown that introduction of interactions in the Peierls system leads to a stronger instability. The mapping between the Heisenberg spin chains with equal nearest neighbor exchange interactions (uniform spin chain) and the spinless fermion model suggests that such a spin chain is also unstable with respect to a lattice distortion leading to alternately strong and weak nearest neighbor exchange constants, i.e., a dimerized spin chain. What is of importance is the fact that such a dimerization is unconditional: no matter how strong the lattice is, the lattice dimerizes, since the exchange energy gained due to dimerization always exceeds the strain energy. This result follows from the fact that the gain in exchange energy varies as while the strain energy loss varies as , where is the magnitude of dimerization that leads to the nearest neighbor exchange constants alternating as .
In recent years, many systems which closely approximate the onedimensional spin chain have been synthesized. What has been observed in these experimental systems is that besides the nearest neighbor antiferromagnetic exchange, there also exists a second neighbor exchange of the same sign and comparable magnitude. Such a second neighbor interaction has the effect of frustrating the spin alignment favored by the nearest neighbor interaction. Therefore, a realistic study of these systems require modeling them using both dimerization and frustration. Theoretically, spin chains with only frustration ( model) were studied by Majumdar and Ghosh. Interestingly, they showed that for , the ground state is doubly degenerate and is spanned by the two possible Kekule structures (Fig. 11). It is quite gratifying to note that a century after the Kekule structure for benzene was proposed, there actually exists a Hamiltonian for which Kekule structure happens to be the ground state!
While most of the discussion above is restricted to spin chains, there has been considerable interest in the higherspin chains following the conjecture of Haldane which predicts that for uniform spin chains, the excitation spectrum of integer spin chains is qualitatively different from that of halfoddinteger spin chains. The latter have a gapless excitation spectrum while the excitation spectrum of the former are gapped. The synthesis and study of integer spin chains have indeed confirmed this conjecture.
Notwithstanding many interesting exact analytic solutions for spin chains, there still exist a large number of situations for which such solutions have been elusive. The exact solutions are basically confined to the uniform Heisenberg model and the frustrated and dimerized model along the line plane, with . Reliable numerical study of these models therefore requires developing techniques which are highly accurate so that the results of large finite systems can be scaled or extrapolated to the thermodynamic limit. As has already been discussed, the DMRG technique is ideally suited due to its high accuracy for quasionedimensional systems. in the
The Hamiltonian for the frustrated and dimerized spin chain is given in Eq. (4) and is schematically shown in Fig. 8. A few lowlying states in a sector with a given value of the total spin component, are obtained in representative points in the plane, using the DMRG method. The ground state is always the first (lowest energy) state in the sector. The accuracy of the DMRG method depends crucially on the number of eigenstates of the density matrix, , which are retained. Working with to over the entire plane gives accurate results. This can be verified by comparing the DMRG results for these values with exact numerical diagonalizations of chains with up to sites for spin systems [82] and sites for spin systems [83]. The chain lengths studied vary from sites for to sites for . The DMRG results are also tracked as a function of , the chain length, to verify that convergence is reached well before sites in all cases. The numerical results are much better convergent for open chains than for periodic chains, a feature generic to the DMRG technique [61, 67].
The quantum phase diagrams obtained for a spin chain is shown in Fig. 12. The system is gapless on the line running from to for , and is gapped everywhere else in the plane. There is a disorder line given by ; the peak in the structure factor