Density matrix renormalization group method

Density matrix renormalization group method DMRG Swapan K Pati Theoretical Sciences Unit JNCASR, Bangalore

Outline Introduction (DMRG and Quantum many body) Numerical renormalization group (K G Wilson) DMRG method ( S. R. White) Extension and a few applications

What is it? Density matrix renormalization group (DMRG) is a numerical technique for finding accurate approximations of the ground state and the low-energy excited states of strongly interacting quantum systems. Its accuracy is remarkable for one-dimensional systems with very little amount of computational effort. It is however limited by the dimensionality or range of interactions. The method is kind of “iterative method” and is based on the truncation of the Hilbert space used to represent the Hamiltonian in a controlled way, keeping the most probable eigenstates! The physical understanding of quantum many-body systems is hindered by the fact that the number of parameters describing the physical states grows exponentially with the number of particles, or size of the system.

More formally in a nutshell: DMRG method • For large systems • Accuracy comparable to exact results • Variational and non-Perturbative • No problems with frustration or fermions • It can calculate: • All ground state properties (energies, correlation functions, gaps, moments) • Finite temperature properties • Classical systems at finite temperature • Dynamical quantities (frequency dependent) • Time evolution Limitations: Convergence depends on details of the system (dimensionality, boundary conditions, range of interactions) and efficient programming is very complicated.

Quantum Many-Body Problem: System of N quantum mechanical subsystems Examples: Hubbard Model:an effective model for electrons in narrow band (eg, -d or –f electron metal ions). It is applicable for atoms, clusters, molecules, solids,….. One-tight binding band, local Coulomb interactions 4N degrees of freedom – states: |0>, , and Interesting ground state properties: (AFM at half-band filling, n=1), 1D: Luttinger liquid; 2D: d-wave superconductivity (n <1?) Dynamical properties like conductivity and temperature dependence are also quite interesting.

Heisenberg model: Strong coupling limit of the Hubbard model at n=1 Antiferromagnetic exchange J=4t2/U Localized QM spin degrees of freedom: (2S+1)N for N spin-S objects. A model to describe quantum magnetism in most of the oxide materials or any system with localized spin orbitals Bethe-ansatz (closed form exact) solution exist only in 1D A good model for describing the parent phase of high-Tc cuprates Ground state, dynamics and low-temperature properties quite interesting.

How does one study many-body interactions? Analytic: Mean-field theories Strong and weak coupling expansions (perturbative methods) Field theoretical methods Mostly uncontrolled Numerical: Exact diagonalization Configuration Interactions and Coupled Cluster Quantum Monte Carlo Dynamical mean-field theory (DMFT) DMRG Extremely involved, each method has its own difficulties Microscopic understanding of systems for applications in magnetic, optical, electrical, mechanical, transport…..phenomena

What is RG Method? The basic idea behind a renormalization group method is to apply a transformation to the Hamiltonian which eliminates unimportant degrees of freedom for the description of the system within a given energy range. For example, if we are interested in the low-energy states of a system with a energy cut-off E, one integrates out energy modes with energy dE is a small energy interval Then we rescale the parameters of the new system so that it reproduces the previous one. Given a H of a system with N variables, a RG transformation Ra is a mapping in the Hamiltonian space which maps H to H’ : H’=Ra(H). H’ now has N’ variables where N’=N/a which is less than N. RG transformation must be unitary; i.e., it has to preserve Z=Tr exp(-H/kT) so that ZN’[H’] = ZN [H].

However, an exact transformation is not possible. Wilson and others: Work in Fourier space and use a perturbative scheme in order to analytically solve this problem. Extremely successful method for solving Kondo problem and Anderson Impurity problem • K. G. Wilson and F. Kogut, J. Phys. Rep C 12, 75 (1974). • K. G. Wilson, Rev. Mod. Phys. 43, 773 (1975). • H. R. Krishnamurthy, J. W. Wilkins and K. G. Wilson, Phys. Rev. B 21, 1044 (1980).

Numerical Renormalization Group (K. G. Wilson, 1974). Can it be applied for Correlated Lattice problem? Idea behind all lattice renormalization group methods is to enlarge the system iteratively but keeping only a constant number of basis states. Integrate out the degrees of freedom numerically for obtaining low-energy properties. Let H be a Hamiltonian describing an interacting electronic systems on a lattice with L sites. Each site has four states: |0>, |down>, |up> & |2>. The dimension of the Hilbert space for L=100 with Nup=Ndown=50 is 1058, which is not intractable numerically. The idea is to obtain the Low-energy eigen-states of this system keeping only a small number of states, say 100.

REAL Space algorithm Let Bl be a block describing the first l sites for which we only keep m states to describe the H. The same goes for Bl’ block also with m’ states. When we put these two blocks together, the H of the new block Bl+l’ has dimensions mm’. Solve H l+l’ and keep only lowest m energy eigenstates • Isolate a finite system (N) • Diagonalize numerically • Keep m lowest energy eigenstates • Add another finite system (N) • Solve (2N) system and iterate the process.

By iterating this procedure, one obtains recursion relations on the set of coupling constants which define the Hamiltonian and the properties in the thermodynamic limit. The message: Low-energy states are most important for low-energy behavior of larger system However, only for Kondo lattice or Anderson impurity models

Very bad for other quantum lattice models: • Hubbard (Bray, Chui, 1979) • Heisenberg (White, 1992, Xiang and Gehring, 1992) • Anderson localization (Lee, 1979) Kondo impurity problem:Hierarchy in the matrix elements; Boundary conditions seem not important. Just have a look for a tight-binding model: 2tii – ti, i+1 – t i, i-1 If one puts two blocks together: does not represent the full system Another way….

Two same size boxes (of length L, 1D, 1-electron problem) Will putting together the ground states of L-length box give rise to the ground state of box of size 2L? NO Treatment of boundary becomes critical: Gr states of a chain of 16 atoms (open) and two 8 atoms chains (filled)

When boundary becomes critical: (White and Noack, 1992) Use combinations of boundary conditions (BCs): Diagonalize a block, HL with different combinations of BCs. Use orthogonalized set of states as new basis. Fluctuations in additional blocks allow general behavior at boundaries. HL Fixed-Fixed Fixed-Free Free-Free Free-Fixed HL HL HL HL 2L Diagonalize superblock composed of n blocks (each block size L) Project wavefunctions onto size 2L block, orthogonalize Exact results as n becomes large

The total number of states is 228 This amazing accuracy is achieved by just keeping only ~100 states! Relative error vs no of states kept

Use of density matrix Steven White, 1992 Divide the many-body system:how? Density matrix projection Divide the whole system into a subsystem and an environment We know the eigenfunction for the whole system: how to describe the subsystem block best? Reduced density matrix for the subsystem block: Trace over states of the environment block, all many-body states.

Density Matrix “When we solve a quantum-mechanical problem, what we really do is divide the universe into two parts – the system in which we are interested and the rest of the universe.” - Richard P. Feynman (Statistical Mechanics : A set of lectures; Westview press, 1972) When we include the part of the universe outside the system, the motivation of using the density matrices become clear. So what does this mean ?

Let be a complete set of vectors in the vector space describing the system, and let be a complete set for the rest of the universe. Now let A be an operator that acts only on the system, ie A does not act on The most general way to write the wavefunction for the total system is Now we have, Density Matrix

Due to the Hermitian nature of ρ it can be diagonalized with a complete orthonormal set of eigenvectors with real eigenvalues We define the operator ρ to be such that, Note that ρ is Hermitian. Again, = Tr Aρ

So, we have <y|A|y> = Tr Ar and If we let A be 1, we obtain If we let A be , we have Therefore, we have and Orthonormal set of eigenvectors and real eigenvalues.

Any system is described by a density matrix ρ, where ρ is of the form and • The set is a complete orthonormal set of vectors. • Given an operator A, the expectation of A is given by Thus, is the probability that the system is in state . If all but one are zero, we say that the system is in a pure state; otherwise it is in a mixed state . Notice that,

Consider a pure state: 1 Notice thatρ2 = ρ and Tr ρ = 1 Consider a mixed state: 50% and 50% 2 Notice thatρ2 ρ and Tr ρ = 1

3 50% and 50% mixture of & Notice thatρ2 ρ and Tr ρ = 1 Diagonal elements = Populations Off-diagonal elements = Coherence It is important to note that in 2 and 3, both cases we describe a system about which we know nothing, ie, A state of total ignorance.

Density matrix: Eigenstates of density matrix form complete basis for subsystem block Eigenvalues give the weight of a state Keep the m eigenstates corresponding to m highest eigenvalues Eigenstate of the whole system thus given by: Schmidt decomposition The optimal approximation Entanglement states (mutual quantum information): DM can be defined for pure, coherent superposition or statistical averaged states

Density matrix renormalization group: Formulation • Diagonalization of a small finite lattice • Division of system • Reduction of the subsystem block via density matrix • Renormalize the matrix formulation of all the operators • Add one or two sites (few possible degrees of freedom) • Construct the bigger lattice • Repeat all the steps • Environment block: • a) Exact sites only • b) Reflection of subsystem block • c) Stored block from a previous step

Technical Details: For nn-Heisenberg Hamiltonian: Given that the ground state of n-sites are known Let |m> with m=1,2,3,……..l specify the complete set of states for the subsystem L and |n>, with n=1,2,3,……k for the rest of the system (n-L). Generally they are considered to be same (L=R=n-L) L R The wavefunction then can be written as We like to find optimal states to represent the subsystem with less number of states than l  a < l.

We wish to minimize by varying aan and La subjected to <La| La’>= da a’ Without loss of generality, we can write In terms of matrix notation, we thus have And we minimize d over all La, Ra and aa. Use Linear algebra: The solution is produced by singular value decomposition Minimization of d: the m largest-magnitude of diagonal elements of D are aaand the corresponding columns of U and V are Laand Ra

These optimal states Laare also eigenvectors of the reduced density matrix of the subsystem as part of the whole system. SVD is now so that U diagonalizesr Thus, the description of the optimal states, La, is best represented by keeping the eigenvectors corresponding to highest eigenvalues of the reduced density matrix. Each eigenvalue (wa) of the density matrix describes the probability of the subsystem being in La, with Sa wa=1 The error then is the deviation of Sa=1,,m wa from unity because of truncation..

Basis states for Ground state 4 Site Problem Density matrix for 2 sites from 4 sites ψG

Basis states for two sites 4 sites ψG can be written in a general way as and since

Reduced density matrix: How does it look? In Matrix representation, the r for 2 sites within 4 sites is Xab are the nonzero elements

Transformation of operators: L R Density matrix for subsystem(L): Diagonalize the density matrix and obtain a transformation matrix O mm’;a composed of m eigenvectors of density matrix Famm’with a=1,……….m. Any operator A ij;i’j’ is transformed as Dimension of X m times m

Iteration: Add two spin-sites in the middle. Left block now consists of L+1 sites and Right block R+1 sites (generally L=R). L l p p+1 r R The basis for left block is now |i> = |La>*degrees of freedom for the new site (s) Similarly for right block, the basis is |j> In this basis, the total Hamiltonian for the L+R+2 sites is: i and i’ are the basis states of left block with new site p j and j’ are the basis of right block with new site p+1

DMRG Algorithm • Make four initial blocks, each consisting of a single site. Set up relevant matrices representing each block. Form the Hamiltonian matrix for the whole system • Diagonalize the Hamiltonian to get target state: Y (m,s,t,n) Expectation values of various quantities can be calculated at this point • Form the reduced density matrix for two-block subsystem (2 sites) using r( m ,s; m ‘ s ‘)=Stn Y( m, s, t, n) Y ( m ‘, s ‘, t, n) • Diagonalize r to find a set of eigenvalues wa and eigenvectors, u (m,s; a). Discard all but the largest m eigenvalues and the corresponding eigenvectors: O=u (m,s; a) • Form matrix representations of various operators of the two-block subsystem. • Form a new block-1 by changing the basis to the u a and truncating to m states using P’ = O P O+. Its dimension is m X m. • Replace old block-1 with this new block 1. • Replace the old block-4 with the reflection of this new block 1. • Go to step -2.

System size 2 ways of building up superblock; choice of environment block Infinite system algorithm: Environment block: reflection of the subsystem block Superblock grows by 2 lattice sites per iteration Finite system algorithm: Starting point: infinite system method Size of whole system same: environment block size shrinks Zipping back and forth: iterative convergence Environment block: use from previous iteration

Infinite system DMRG System size grows by 2 sites every iteration

Finite system DMRG Start: Infinite DMRG generated system Left block grows Right block shrinks Reflect the whole system Left block grows Right block shrinks First iteration ends

Few Remarks • Most accurate method for interacting systems in low-dimensions. • Higher dimension: if X states are required for obtaining accurate results in one-dimensional systems, X d states are necessary to obtain comparable accuracy in d-dimensional systems. • All kinds of Static, Dynamical, Thermodynamic properties can be calculated. • It is not a RG technique in strict sense. No fixed point or RG flow. • The Hamiltonian matrix that one encounters from iteration to iteration, remains roughly of the same order, but the matrix elements keep changing. In this sense it can be called a renormalization procedure  Coupling constants keep changing while the system size increases, as in the RG procedure in a blocking technique.

While writing the code: Use quantum numbers: Any basis states can be partitioned by various quantum numbers. All operators, in general, consist of blocks connecting definite quantum numbers. Eg. Sz total and N total (fermionic case) Operators can be stored in sparse matrix form. Do all arithmetic's using sparse matrices. Non-abelian quantum numbers, like S2, also possible (though complicated) Disk usage: Information not used in current loop may be written to disk (eg, infinite DMRG generated blocks for finite DMRG)

Applications: • Strongly correlated electronic systems • Nuclear Physics • Quantum information theory • Quantum Chemistry • Classical Statistical Physics • Soft condensed matter Physics Total number of papers published with the string “density matrix renormalization” in their title or abstract from 1993 to 2005 is more than5,000(obtained from ISI database) Distributed Multimedia Research Group Design Methodology Research Group Direct Marketing Resource Group Data Management Resource Group Groupe de Renormalisation de la Matrice Densite (GRMD) !!! All are DMRG indeed !

Original Reference and some review articles S. R. White, “Density matrix formulation for quantum renormalization group” Phys. Rev. Lett. 69, 2863 (1992). S. R. White, “Density matrix algorithm for quantum renormalization group” Phys. Rev. B48, 10345 (1993). Steven White was awarded the 2003 Aneesur Rahman Prize for Computational Physics by the American Physical Society for the development of DMRG method. S. Ramasesha, S K Pati, Z. Shuai, J.L. Bredas, “The DMRG method: application to the low-lying electronic states of conjugated polymers”, Adv. in Quant 38, 120 (2000). S K Pati, S. Ramasesha and D. Sen, “Exact and approximate techniques for quantum magnetism in low-dimensions”, as a chapter in “Magnetism: from molecules to materials”, edited by M. Drillon and J. Miller, 2002. U. Schollwock, “The density matrix renormalization group”, Rev. Mod. Phys. 77, 259 (2005).

Symmetrized DMRG Methods

Try to exploit different symmetries of Model Hamiltonians Some possible symmetries: [H, J] = [H, C2] = [H, P] =0 Electron – Hole Symmetry (J): only for ½-filled fermionic system End – to – End interchange or 180o rotational symmetry ( C2 ) Spin-Parity symmetry ( P ): Only for Sztot=0 How do they operate? Electron – Hole symmetry operator interchanges the creation and annihilation operators at a site with a phase. Example : ai† = ( -1 )i bi Transformation of the Fock space of a single site under J i Ji | 0 > = | ↑↓ > ; Ji | ↑ > = ( -1 )i | ↑ > Ji | ↓ > = (-1 )i | ↓ > ; Ji | ↑↓ > = (-1 ) | 0 >

C2 symmetry operation : • Rotation through an angle 1800 • Interchanges the states of the left and right halves of the system with a phase factor. Consider a state for system+Env: | µ, σ, σ’, µ’ > C2| µ, σ, σ’, µ’ > = (-1 )γ | µ’ , σ’ , σ , µ > γ = Phase= ( nµ + nσ ) ( nµ’ + nσ’ ) Where m (m’) refers to the eigenvectors of the left (right) block DM and s (s’) refers to the Fock space states of the new left and right sites added. nµ , nµ’ , nσ and nσ’ are the occupancies in the states | µ > , | µ’ >, | σ > , | σ’ >

Spin-Parity operation : • Flips the electron spin at a particular site. Transformation under spin-parity operation : Pi | 0 > = | 0 > , Pi | ↑ > = | ↓ > Pi | ↓ > = | ↑ > , Pi| ↑↓ > = | ↓↑ > The full spin-parity operator ( P ) is the Direct Product of single site parity operators (Pi ) : P = ∏ i Pi This operator bifurcates the Sztot=0 spin states into Even and Odd total Spin states.

If we use all the three symmetry operators, together with Identity, they Form an Abelian group with 8 irreducible representations. eA+, eA- , oA+ , oA- , eB+ , eB-, oB+ , oB- Projection operator for a given irrudecible representation Г PГ = 1/h ∑ χГ ( R ) R Where R ‘s are the symmetry operations, χГ ( R ) is the character of R inГand h is the order of the group. Since the symmetries do not exist for “only system” or “only environment” blocks, the implementaion of the symmetries require “symmetrization” and then “unsymmetrization” steps in every DMRG iteration.

1. The Matrix represenatation for symmetry operators of the full system with (2n +2 ) sites in direct product space are obtained as the direct product of corresponding matrices. < µ , σ , σ’ µ’ | R2n+2 | ν,ح , ح ’ , ν’ > = < µ | Rn | ν > < σ | R1 | ح > < σ’ | R1 | ح ’ > < µ ’ | Rn | ν ’ > • Matrix representation in the basis of eigen vectors of new density matrix : Matrix ( R ) • as a direct product of Rn and R1 given as < µ, σ | Rn+1 | ν, ح > = < µ | Rn | ν > < σ | R1 | ح > 3. This matrix is Renormalized by the Transformation Rn+1 = O† Rn+1 O O is the Matrix whose columns are the chosen eigenvectors of the density matrix • The Hamiltonian matrix in the direct product basis can be transformed: where S is the symmetry elements matrix in the direct product basis. H2n+2 = S† H2n+2 S 5. The wave function then is back transformed to unsymmetrized basis: |ysymm> = S-1 |yunsymm>

This symmetrized DMRG scheme has been implemented in Hubbard Model for the calculation of excitation gap. Ref : Phys. Rev. B 54, 7598 (1996) Energy gaps ( measured from the ground state ) of the lowest state in each subspace for chain length varying 40 to 50 , for two different values of U/t. Splitting of the Hubbard bands

Density matrix renormalization group method