The Recursive Solution of the Schrodinger Equation

SOLID STATE PHYSICS. VOLUME

35

The Recursive Solution of the Schrodinger Equation ROGERHAYDOCK Caomdish Laboratory, University of Cambridge, Cambridqe, England

...........

1. Basics.. . . . . . . . . . . . . . . . . . . . . .

........... .................................................... ........................... a Chain.. . . . . .

1 1. Irregular Solutions on a Chain

14. Symmetry.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15. Vector Chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11. Green Functions and the Local Density of States.

.

17. Green Functions of a Finite Chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18. Green Functions of an Infinite Chain and Their Integrals.. . . . . . . . . . . . . . . . . 19. Examples of Green Functions of Infinite Chain Models. . . . . . . . . . . . . . . . . . . . 111. Perturbation Theories, Energy Differences, and Computing the Density of States . 20. Perturbation Theory for the Local Density of States, Van Hove Singularities, and Chain Asymptotics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21. Singularities in a Band of Finite Width 22. Computing the Local Densities of States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23. The Average (Total) Density of States. . 24. Differences in Energy and Other Integra 25. Binding Energy.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .... 26. Perturbation Theory for the Chain 27. Perturbation of the Green Function nts . . . . . . . . 28. Perturbation Theory Applied to Integrated Quantities. . . . . . . . . . . . . . . . . . . . .

216 216 218 220 222 225 228 229 23 1 234 236 240 242 243 244 247 247 2 50 250 255 259 263 263 270 270 215 216 217 281 283 286

215 Copyright @ 1980 by Academic Press, Inc. All rights of reproduction in any form reserved. ISBN 0-12-607735-5

216

ROGER HAYDOCK

IV. Relation to Other Methods and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29. Relation to Other Methods.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30. The Method of Moments .......................... 3 I . Transfer Matrices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32. Feenberg Perturbation The 33. The Cluster Bethe Method 34. Equation of Motion Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

287 287 289 290 29 1 29 1 292 292

I. Basics 1. INTRODUCTION

The local environment approach to the electronic structure of solids requires an alternative to band theory for solving the Schrodinger equation. In the limit of weak interaction between electrons and atoms, for example in simple metals, the electronic structure is determined by the long-range periodicity of the atomic potentials. Band theory exploits this aspect of the physics to express electronic properties of the solid as a coherent superposition of the electronic properties of all the atoms. When the electrons interact strongly with the atoms, this picture breaks down and properties no longer depend on long-range periodicity, but rather on only the first few shells of neighbors of each atom. The d electrons in transition metals are prime examples of this regime. While band theory is still a valid, formal solution to the Schrodinger equation, the physics is better understood by means of a solution that explicitly accounts for the role of local environment. It is just such a method of solution of the Schrodinger equation that we present in this chapter. Electronic orbitals that are localized near each atom or bond can provide a basis for the electronic states of the solid. The atomic or bond properties that are independent of environment determine the local basis orbitals. Their interactions, expressed by the tight-binding model, determine the way in which electrons couple throughout the solid to produce its properties. The author of Chapter 2 discusses the various ways of constructing such a basis, while in this chapter we concentrate on determining the properties of the solid from those of the basis orbitals. The local density of states describes the effect of the rest of a solid on one region. Physically, the local density of states is the intensity of each eigenstate on a particular atom or bond. Mathematically, it is the magnitude squared of the projection of each eigenstate on a local orbital, so that the local density of states gives the coupling of an orbital to eigenstates at each

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

217

energy. Most properties of solids including binding energy, photoemission spectra, and magnetism can be related to the local density of states. The basic problem is to find the local density of states from a tightbinding or localized orbital model. The physics demands that it be done in such a way that the local orbital itself has the greatest effect, and that successively more distant orbitals have lesser effects. The solution must define a hierarchy of environments so that their relative influence is explicitly displayed in the local density of states. The solution of this problem is the determination of linear combinations of the local orbitals such that the electron must pass through each state as it diffuses away from the local orbital from which it starts. The local density of states for a basis orbital is determined by the sequence of environments described by orthogonal states I U,), where I U,) is a linear combination of basis orbitals, and I U , ) is 1 4,). The 1 U , ) must satisfy the condition that

%lU,>

= a,IU,)

+ b,+llU,+l) + ~ , l U , - l ) 7

(1.1)

where % is the Hamiltonian of the model. In practice the states I U,) are localized on the shell of atoms n hops from the atom accommodating 1 40). The parameters a, and b, describe the coupling of each environment to itself and its neighbors. The powerful result that justifies this chapter is that any model can be transformed into this pseudo-one-dimensional chain model where the importance on I 4,) of each I U,) is unambiguous. This is difficult to believe until one realizes that any matrix can be transformed into a tridiagonal or Jacobi matrix and that this is all the method involves. However, the numerous applications of this idea and the physical implications of a chain model for any system have only just begun to be exploited. The generality of this result suggests a new way of looking at quantum mechanics. Since any quantum system can be transformed into a chain model, we need only investigate such chain models in order to see the varieties of quantum phenomena that are possible. To understand a particular physical system, we need only find the chain model appropriate to it. There is a large body of mathematics that bears on the analysis of the chain model. The subjects of orthogonal polynomials, classical moments, and continued fractions have numerous useful results. For example, the local density of states for I 4,) can be written in terms of the parameters of the chain model as the continued fraction no(E) = -Im(l/{E - a, - b:/[E- a , - h:/(E - a , - ...)I})/n, (1.2) where the energy E has an infinitesimal positive imaginary part. A broad and well-developed mathematical base has greatly speeded the application of this approach as well as giving confidence in the results.

218

ROGER HAYDOCK

The recursive solution of the Schrodinger equation is well adapted to the digital computer. Nowadays, computations of experimental quantities are produced along with any theory as a matter of course. Being discrete, the methods of this chapter are easily implemented on digital computers and a library of Fortran programs is available.’ To cover this new approach to quantum mechanics completely would take a whole volume, and so in this chapter we develop the topics of most immediate use to solid state physicists. We assume that the reader is familiar with the basic ideas of wave and matrix mechanics. The results are stated as precisely as is consistent with making their physical applications clear, but we have substituted short outlines of the arguments instead of formal proofs. There are still many equations needed for those who would wish to apply the methods to specific problems. We hope these will not discourage the general reader from finding out what the methods are about. Formal proofs for these results, where they are available, are to be found in Akhiezer,’ Shohat and T a m a r k i ~ and ~ , ~ S z e g ~along , ~ with a great many other useful results. The material on asymptotic approximations (cf. Section II1,20) has not, so far as we know, been developed formally. The chapter is organized in order of sophistication of the applications. The elementary part finishes with Section 11,17, on Green functions. Thereafter the topics apply to particular kinds of calculations such as local density of states, energy differences, and photocurrents. Material on statistical averaging of models has been left out altogether due to lack of space. The chapter finishes with a summary of fruitful directions for further work. 2. NOTATION

There are a number of reasons for adopting a notation slightly out of the ordinary. The methods to be described are particularly adapted to computations on a digital computer and our notation makes that implementation straightforward. It is an economical and thus abstract notation, and for this reason we devote this section to explaining it and relating it to more familiar ways of expressing quantum mechanics. The computer is a discrete machine and hence not suitable for calculations on continuous quantities. Quantum mechanics and many other physical

’ C. M. M. Nex, The Cambridge Recursion Library.

N. I. Akhiezer, “The Classical Moment Problem.” Oliver & Boyd, London, 1965. J. A. Shohat and J. D. Tamarkin, “The Problem of Moments,” Math. Surv. I, rev. ed. Am. Math. SOC.,Providence, Rhode Island. 1950. G. Szego, “Orthogonal Polynomials,” Colloq. F’ubl. XXIII, rev. ed. Am. Math. SOC.,Providence, Rhode Island, 1957.

’

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

219

theories involving linear equations are usually expressed in terms of continuous mathematics : wave functions, differential equations, and integrals. These theories can be formulated in terms of an infinite, discrete set of basis functions rather than the continuous bases of position or momentum (plane waves). Formulated discretely, the mathematical representation is that of infinite matrices and the notation is that of infinite-dimensional matrix algebra. As well as being discrete, the computer is also a finite machine and while the formulation of the theory is in terms of an infinite basis, the implementation must always be finite. The crucial physical approximation in every calculation is the choice of finite basis to be used. One of the advantages of the approach we take here is that this physical approximation is clearly distinguished from the various numerical approximations that are usually made in a computation. The choice of discrete basis has been covered in Chapters 1 and 2 and we say no more about it here other than to define it formally as a set of wave functions, not necessarily orthogonal, {40(r), 41(r), &(r), . . . , q$,(r), . . .} where r represents the coordinates of the system. We refer to this basis when relating the matrix mechanics to wave mechanics. The first constituent of the notation is the states represented by column vectors. A state u is represented as an infinite column vector with elements uo, ul,u 2 , . . . . The basis relates it to a wave function u(r) by

The second constituent of the notation is the linear operators, which are represented by matrices. A linear operator A is represented as an infinite square matrix with elements {amn:m, n = 0, 1,. . .}. This is related to an operator d from wave mechanics by amn =

(2.2)

J 4 z ( r ~ 4 n ( r )dz,

where the integral is over all of coordinate space. The most important operator for our purposes is the Hamiltonian. The time-independent Schrodinger equation in wave mechanics,

{C - h2/(2m)Iv2 + J‘(r)>ua(r)

=

Eaua(r)

(2.3)

becomes, in matrix notation, Hu, = E,u,,

(2.4)

where H is the Hamiltonian matrix and CI an index containing all the quantum numbers.

220

ROGER HAYDOCK

The third constituent of the notation is the adjoint operation, denoted by a superscript dagger (t).A state u has an adjoint ut, which represents a row vector whose elements are the complex conjugates of the elements of u. A linear operator A has an adjoint A t whose matrix is the complex conjugate transpose of the matrix A . Treating scalars as 1 x 1 matrices, the adjoint of a scalar is its complex conjugate. Row and column vectors interchange, and square matrices go to other square matrices. Taking the adjoint of the adjoint gets one back to the original matrix, vector, or scalar. We combine states and linear operators according to the rules of matrix algebra. For example, two states u and v combine to give a scalar representing their overlap s by s = V+U, (2.5) i.e., multiplication of a row vector and column vector. If u(r) and u(r) are the wave functions corresponding to the states, then s=

Jv(r)*u(r) dz.

Another property is that the adjoint of a product is the product of the adjoints in reverse order. We occasionally take a liberty with the notation so far defined,and that is to refer to states and linear operators, expressed in terms of a continuous basis such as position or momentum, just as if the basis were discrete. In this case the summations implied by matrix multiplication go over to integrals with the appropriate measure. This completes the brief description of the notation. Further information can be found in any textbook on matrix algebra. 3. THEHAMILTONIAN AND ORTHOGONALITY

The purpose of this section is to define in familiar terms the matrices and vectors that make the recursive computational methods efficient. The author of Chapter 2 discussed a number of methods for constructing a local orbital basis for the Hamiltonian. In general these methods produce a set of electron orbitals {q5,,(r)} that are not orthogonal, in the sense that

/q5:(r)q5,,(r)dT # 0,

for m # n.

(3.1)

To construct an orthogonal basis such as Wannier states out of linear combinations of the {&(r)} would, first of all, obscure the chemical distinction of an orbital that is localized on one atom or bond. Second, it is a long computation, which is completely unnecessary.

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

22 1

The first matrix needed is the overlap matrix S for the basis. Its elements are defined by

CSI,,

=

/ M r ) M ) dz.

When we write two states u and v in terms of the {$,(r)},

The physical overlap or inner product of these states is given by

vtsu.

(3.4)

This expression relates the overlap of the two states to the overlaps of the basis orbitals. The expression vtu is not a physical overlap. It is just a formal summation. Throughout this chapter, when we refer to states u and v as orthogonal, we mean physically orthogonal, and hence that

v+Su = 0.

(3.5)

The second matrix required for computation is that of the Hamiltonian operator YF. When this acts on a basis orbital 4,,(r), it produces a state that we approximate or in some cases write exactly as a combination of basis orbitals 00

#4n(r) = C Hrnn$rn(r)n=O

(3.6)

By the Hamiltonian matrix H,we mean the matrix These parameters H,, are not to be confused with matrix elements of the Hamiltonian

When the {@,(r)} are not orthogonal, these do not give the expansion of X&(r). The YFmnare related to H,, by the overlap matrix S,

222

ROGER HAYDOCK

The matrix S must be self-adjoint, because (3.10)

and similarly,

Z,,= A?,*,.

(3.1 1)

In general H is not a self-adjoint matrix. In the case where the {q!+,(r)} are orthogonal, then S is the identity and the H,, and Z . , are the same. Many of the methods in Chapter 2 give H directly. and for the methods of this chapter, In some cases, one knows S and [Z,,] one needs to compute H.The relationship is

H = S-'[Z,,].

(3.12)

In many cases S has properties that make this computation easy. If the basis {$,,(r)} is normalized, then the diagonal elements of S are unity, and in general the off-diagonal elements are relatively small. If the {q$,(r)} are well localized, S- falls off rapidly in real space.5Whatever approximation is used for H,it must have the property that SH is self-adjoint. A reasonable way of doing this is to approximate S- by a self-adjoint matrix T. We ensure that SH is self-adjoint by approximating H by

H = T[iW,,]TS.

(3.13)

Methods for doing this computation are part of standard libraries.

4. THECHAINMODEL Any quantum-mechanical model can be transformed into one or more chain models, which can then be solved. In this section we introduce and discuss the chain model before going on in other sections to develop the transformation that makes it universal. A quantum-mechanical model consists of a Hamiltonian, a set of orbitals on which that Hamiltonian acts, and a physical interpretation of the orbitals. The Hamiltonian describes the motion of the system among the orbitals of the model. The orbitals and their interpretation relate this motiotl to the physics. We assume that the Hamiltonian is hermitian, conserving probability, and that the overlap of orbitals is defined by some hermitian inner product. The chain model is specified by a sequence of orthonormal orbitals {uo, u,, . . .} and two sets of real parameters {ao, a l , . . .} and { b l ,b Z , .. .}, P.W.Anderson, Phys. Ref.. 181, 25 (1969).

RECURSIVE SOLUTfON OF SCHRODINGER'S EQUATION

223

FIG.1 . A graphic representation of the chain model including the states {un}and the recurrence parameters {an}and {&}.

which describe the action of the Hamiltonian H on the orbitals by a symmetric three-term recurrence relation,

Hun = a,u,

+ b,+lu,+l + b,u,-l.

(4.1)

The recurrence is called symmetric because the component of u , + ~in Hun is the same as that of u, in Hun+1. This is partly a result of the hermiticity of the Hamiltonian. The overlap of orbitals is defined in terms of a matrix S so that u;su,

(4.2)

= h,,.

If the basis for the matrices and vectors is an orthogonal one, then S is the identity and H a self-adjoint matrix. However, as indicated in Chapter 2, many methods of constructing local orbitals produce nonorthogonal bases. In these cases S is self-adjoint, reflecting hermiticity of overlap, and while the matrix of H is not self-adjoint,SH is, reflecting hermiticity of the Hamiltonian. The sequence of orbitals {uo, . . . ,u,, . . .} may be finite or infinite, and for the first and last recurrence relations, orbitals u- and uN+ are taken as zero. The chain model gets its name from its graphic representation shown in Fig. 1.The orbitals are represented by vertices (nodes) and the Hamiltonian is represented by the edges (connections between vertices). Graphical representations of models are interesting as a separate subject. Any Hamiltonian may be represented in terms of a set of orbitals, the vertices, and the nonzero elements of the Hamiltonian, the edges. For instance, a one-electron treatment of a square lattice of atoms with interactions only between neighboring s orbitals may be represented by a graph that divides the plane into equal squares. The associated recurrence relation is no longer three term, but the methods we apply to the chain can be generalized and applied to more complicated models. However, for computational purposes, the chain model is best. The chain model is equivalent to expressing the matrix of H as a Jacobi matrix (a J matrix or tridiagonal matrix). To see this, we relate u, to a column

,

,

224

ROGER HAYDOCK

vector filled with zeros except that the nth element is unity. The matrix representation of H in this basis follows from Eq. (2.1) and gives

a. bl H= 0 0

.. .

0 0 b2 a2 b, 0 b3 a, b,

.. .

0 b2

.. .

.. .

... ... ...

...

(4.3)

Thus the chain model is a symmetric representation of H . There is also a physical interpretation of the chain model. The orbital u, represents the initial state of the system, for example, an electron on a particular atom in a solid. The chain model permits the system to hop from u, to u,, which could represent a linear combination of orbitals on atoms neighboring u,. Similarly, having reached u,, the system could either hop back to u, or on to u2, representing (say) second nearest neighbors of the original atom. The chain establishes an ordering of states of the basis in terms of how many intermediate states must be passed to reach a given state. Thus the system must pass through u, on the way to u2 and so forth. Another way of expressing this is to say that the behavior of the system, initially at u,, is most strongly influenced by u, and less so by each succeeding orbital in the chain. The stronger influence on motion of the system comes from orbitals nearer in the sense of being accessible via fewer intermediate orbitals. The chain model is the mathematical expression of the physical concept of local environment, introduced in Chapter 1. Each orbital on the chain represents a more distant part of the environment of the initial orbital, and the parameters of the recurrence specify the effect of that environment on the motion of the system. The chain model describes the evolution of the system from a given initial state. Orbitals that do not couple to this initial state do not appear in the chain. The orthonormality of the orbitals in the chain simply ensures that the recurrence will be symmetric for a hermitian Hamiltonian. There is a more general chain model6 that does not depend on any inner product or definition of overlap. It can be used for nonhermitian models or to describe transitions between states of a system. In this chapter we deal only with the symmetric chain model because of its numerical stability and because there are several outstanding problems concerning unsymmetric chains. A chain is a generalized eigenstate. The simplest chain contains one orbital uo, and in this case Eq. (4.1) becomes

Hu,

' R. Haydock, J . Phys. A 7,2120 (1974)

= aouo,

(4.4)

RECURSIVE SOLUTION OF SCHR~DINGER'SEQUATION

225

which means that uo is an eigenstate of H.A longer chain represents an invariant subspace, a set of states such that H acting on any state gives another state in the subspace. The chain model corresponds to the smallest invariant subspace containing the initial state, or for that matter, any other state in the chain. An example of a hypothetical physical system that leads directly to a chain model is a string of hydrogen atoms. In a one-electron treatment of this problem, we allow the electron to hop from one s orbital to another on a neighboring atom. In this case all the {a,} are the s orbital energy and the {b,} are the hopping element from one atom to the next.

5. TRANSFORMATION TO A CHAIN The transformation of quantum-mechanical models to chain models can be done analytically in some cases and numerically in other cases. While the basic idea of the transformation is simple, starting and ending it involve some complications. In this section we discuss initiating, recurring, and terminating the transformation. The transformation to a chain has been known in one form or another for over 100 years. During this time it has been discovered independently by people working in many different areas. It is not always obvious that the underlying mathematics is the same, but these methods have been applied by Whitehead in nuclear physics,' Mori in statistical mechanics,* as well as many others. In solid state physics, it was suggested by Haydock, Heine, and Kelly' as a means of calculating electronic densities of states. Earlier Lanczos" suggested it as a means of diagonalizing matrices by computer. The three-term recurrence is the method of solving the classical moment p r ~ b l e m ,and ~ . ~Chebyshev" applied it to the uniform approximation of a function on an interval. The starting point of the transformation is some general quantum model that we wish to investigate by converting it to a chain model. The general model consists of a set of states, their physical interpretation, and a Hamiltonian that determines the motion of the system among its states. The model may be specifiedby wave functions and a differential operator, or by a set of orbitals and matrix elements. Unless the physics is one dimensional in some fundamental way, the Hamiltonian will connect each state with

' R. R. Whitehead, A. Watt, B. J. Cole, and D. Morrison, Adi,. Nutl. Phys. 9, 123 (1977). * H. Mori, Jpn. Prog. Theor. Phys. 34, 399 (1965). R. Haydock, V. Heine, and M. J. Kelly, J . Phys. C8, 2591 (1975). C. Lanczos, J . Res. Natl. Bur. Stand. 45, 255 (1950). " P. L. Chebyshev, Bull. Acad. Imp. Sci. St. Petersbourg 1, 193 (1859). lo

226

ROGER HAYDOCK

several or even an infinite number of other states rather than the two of the chain model. We express the general model in terms of matrices. The states are represented by column vectors, and in order to relate overlap of states to matrix multiplication, there is an overlap matrix S whose elements are the overlaps of the various Cartesian basis vectors (vectors with one 1 and the other elements zero). S must be self-adjoint and positive definite to represent the hermitian inner product of two states u and v,

(v, u) = vtsu.

(5.1)

The Hamiltonian is represented by a matrix H so that Hu is the state resulting from acting with the Hamiltonian on u. The elements of H are to be distinguished from the “matrix elements of the Hamiltonian,” which in this notation are given by SH (cf. Section 1,3). The hermiticity of the Hamiltonian makes this product self-adjoint although H alone may not be. The first step of the transformation is the choice of a starting state u,. It determines what physics the chain represents. For example, if we wish to explore the electronic structure at a metal surface, then u, should be a state localized at the metal surface. If we want the electronic bond order between two atoms, then u, should be a linear combination of orbitals on the two atoms. Each choice of u, produces a different chain, and so we choose u, to give a chain containing the information we want. The second property of uo is that u&SH”u,must be finite for all finite n. For any finite matrices S and H,this will be true, but in the case of analytic chain models, this condition can be important. What is says is that the high-energy components of uo die off rapidly. For convenience we assume that u, is normalized, so that u,su,

=

1.

(5.2)

The transformation proceeds by constructing ul, the nearest environment of uo, by using the recurrence Eq. (4.1),

Hu,

=

aouo

+ blU,.

(5.3)

The left-hand side is known and a,, bl, and u1are to be determined. Applying orthonormality determines a. :

a, = u&SHu,. Subtracting a. u, from both sides isolates blu,,

b,u,

= (H - a0)uo.

(5.4) (5.5)

Orthonormality fixes b:,

b: = [(H- ~ o ) u o ] ~ s [ (-Hao)uo].

(5.6)

227

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

b: is greater than or equal to zero because the inner product is positive definite. Suppose 6: is greater than zero. It is convenient to take b, as the positive root of b:, which gives

u1 = ( H - ao)u,/b,. (5.7) We observe that by construction u, is normalized to one and orthogonal to uo . Thus we have calculated the first recurrence in the chain. For the general step of the transformation, suppose we have constructed orthonormal states {uo, ul, . . . ,u,} along with real parameters {ao,a,, . . . , a,-,} and {b,, b 2 , .. . , b,}, so that no b is zero. In this step we construct u,+ I , a,, and b,+ by again using the recurrence Eq. (4.1):

+

+

Hu, = a , ~ , b,+lu,+l b,u,-l. (5.8) We calculate Hu, and now project out its components of u, and u , , - ~ . Using some matrix algebra and the self-adjointness of S,

uL-

,SHu,

=

(Hun- l ) t S ~ = , (u~SHU,-1),

(5.9)

which is just b, as constructed in the previous step. Likewise, a, = uiSHu,.

(5.10)

Solving for the new state, we get

bn+lun+l = ( H Normalizing the new vector gives b;+,

=

[(H

-

a,)u,

-

-

(5.11)

a,)u, - b,U,rl.

b,u,- ,ItS[(H - a,)u,

-

b,u,-

,I.

(5.12)

This is positive (assuming the new vector is not zero) because Eq. (5.12) represents the inner product of a vector with itself. Again, take b,+, as the positive root, leaving

u,+1 = C(H

-

a,)% - bnun-ll/bn+l.

(5.13)

By construction u,+, is normalized to one and orthogonal to u, and u , - ~ . However, to extend the chain it must also be orthogonal to u , , - ~ ,. . . ,uo. We show this by considering the expression

,

(5.14) b,+ l ~ L S ~ , = + uASHU,, obtained from Eq. (5.11) with m less than n - 1. Using the self-adjointness of S we get b, + ,u: sun+

, = u: HtSU, = (HU,)tSU,.

(5.15)

Now H operates on urnand the result is a linear combination of urn-1, urn, and urn+ all of which have zero overlap with u,, and so u,+ is orthogonal

,

228

ROGER HAYDOCK

to all the previous states in the chain. It is clear that the construction of the chain is very subtle compared to, say, Schmidt orthogonalization, which uses all states at each stage. The transformation terminates when the new state has norm zero. In Eqs. (5.6) and (5.12) we bypassed the possibility that b;,, could be zero. Because S is positive definite, this can only occur if (H - a,)u, - b , ~ , - is ~ zero. Termination of the chain means that the smallest invariant subspace containing the initial state is finite dimensional, and spanned by the chain states. In most solid-state systems, the chain will not terminate. 6. A FREEELECTRON A simple analytic example of the transformation to a chain model is the free electron in three dimensions. We take a liberty with the notation in this example to refer to a continuous basis, namely, the position. The Hamiltonian is the usual one: = -h2V2/2m.

H

(6.1)

Suppose we are interested in the behavior of an electron initially localized about the origin. A suitable choice for initial state is

(n/Jlr>3’2 exp[ -(nr)’/2],

(6.2) where A- is a localization length and r the position. uo has been normalized to 1 and is a simple Gaussian. Since uo is spherically symmetric and H preserves spherical symmetry, we can proceed in spherical coordinates : uo =

H

=

~

(-

d2 + z”) dr2 r dr ’

-h’ 2m

and An appropriate unit of energy is E =

h2A2/2m.

Proceeding with the transformation, a. = u&Huo= 4 ~ ( n / & ) ~

I

a,

0

which gives

-c[(Ar)’ - 31 e ~ p [ - ( k ) ~ ] dr, r ~ (6.6)

229

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

Further results are

bl

=

~1

= - (3/2)”’[(h)’

;)!.!I”’(

(3/2)’/’~,

(6.8) -

3/2]u0.

(6.9)

The general step of the transformation is

u,

=

[(;)!(.

+

m=O

(-l)m(Lr)’m/[(n

-

m)!(m

+ ;)!m!]]uo, (6.10)

+ a, = ~ ( 2 + n 3/2). b,

=

~ [ n ( n 1/2)]”’,

(6.11) (6.12)

The expression for u, is proportional to the original Gaussian modulated by the Laguerre’ polynomial L,,!/’(L’r’). The transformation does not terminate because there are an infinite number of states available to the electron. However, the weight of u, moves progressively outward from the origin. This example demonstrates a number of properties of the transformation. The transformation is economical because it only introduces states that couple to the original state. The coupling between the states in the chain grows larger as we move outwards, reflecting the fact that a localized electron decays rapidly in free space.

’

7. NUMERICAL TRANSFORMATION TO A CHAIN

The numerical transformation to a chain is extremely stable and economical of storage. In this section we discuss the implementation of the transformation on a digital computer and show why the results are so much better than one might otherwise expect. Most quantum models cannot be solved analytically, and so we must resort to computation. Systems with high symmetry can often be reduced to a large number of few-dimensional problems as with band theory. However when the symmetry is low or nonexistent, a large number of degrees of freedom (dimensions)may be strongly coupled. For the standard matrix diagonalization techniques, the amount of store goes up with the square of the number of dimensions. If the Hamiltonian can be stored in a compact form-as is often the case because of its simplicity -then the chain transformation (Lanczos method) is highly advantageous because the store needed is only twice the number of dimensions. Because the chain contains successively less important environments of the initial state, l2

S. 1. Gradshteyn and I. M. Ryzhik, “Table of Integrals, Series, and Products.” Academic Press, New York. 1965.

230

ROGER HAYDOCK

the physics is often contained in only the first few states and so the whole transformation is unnecessary. Until 1972 this transformation was believed to be numerically unstable. Calculation of d band densities of states’ showed no signs of this instability, and a study of the symmetric Lanczos method by PaigeI4 showed that there was an algorithm that gave extremely good results for eigenvalues of matrices reduced to Jacobi form by this transformation. His analysis of the propagation of rounding error showed that computation of b,+ by taking the square root of the right-hand side of Eq. (5.12) and then using the same value in the calculation of u,+~ (each b,+ appears in two recurrences and one could conceivably recalculate b, + the second time) was the most stable procedure. This is also the fastest form of the algorithm. Even with the stable algorithm, there is a cumulative decay of orthogonality in the chain, due to rounding error. Suppose H is an N-dimensional matrix with no degeneracy. If uo contains a bit of all the eigenstates, transformation to a chain should terminate after u N - , because all degrees of freedom are exhausted. What happens in practice is one finds that bN is small but nonzero and the chain continues. This produces a section of chain modeling the system and weakly coupled to the first part of the chain. This is called a “ghost” because it gives rise to multiple eigenvalues for each eigenvalue of the original matrix. These numerical ghosts continue to be produced as long as the transformation continues. The cumulative loss of orthogonality does not make the chain a less accurate model of the original system until linear independence of chain states is lost. Although nonorthogonality accumulates, eigenvalues, etc., are accurate to a single rounding error. As long as {uo, . . . ,uN- are linearly independent, we can construct a linearly independent set of vectors {vo, . . . ,vN- 1 } with the property

,

’

VJU,

=

.,a

(7-1)

There may be rounding errors in the computation of {ao; a,, . . . , a N - 2 }and {bl, b 2 , .. . , b N - , } but the components on both sides of each recurrence

Hun = a,u,

+

bn+lu,+l

+ b,u,-1,

0 I n I N - 2,

(7.2)

are equal to within the rounding error. The states may not be quite orthogonal, but Eq. (7.2) does not depend on their orthogonality. Using the fact that v: Hun = (Htv,)tu,, l3

l4

0 I n IN - 2,

R. Haydock, V. Heine, and M. J. Kelly, J . Phys. CJ, 2845 (1972). C . C . Paige. J . Inst. Math. Appl. 10, 373 (1972).

(7.3)

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

23 1

we see that the v obey a recurrence for Ht, Thus

HtVn= u,V, v:HuN-,

+ b,+ Iv,,+ + b,v,= 0,

= (Htv,JtuN_1

0I nI N - 2.

(7.4)

for 0 I H 5 N - 3,

(7.5)

1,

and SO H U N - can only have components on vN- and vNlate :

uN-1

=

’,which we calcu-

~lt-1HuN-l.

(7.7)

This now determines the last recurrence,

(7.8) which is unsyrnmetric6.l5 in that c N - , may not equal b N - , . It shows that despite an accumulation of rounding error in the orthogonality of the states, the chain can be terminated so as to be accurate with a simple rounding error rather than accumulated rounding error. The reason for this atonishing result is that the orthogonality of the states in the chain model is not really necessary, but only a convenience for some other uses. The real requirement is that the states be linearly independent. A general relaxation of orthogonality would involve us in the unsymmetric chain, which is beyond the scope of this chapter. Making an unsymmetric correction at the end of a chain is a completely practical way of eliminating cumulative rounding errors although in most applications one does not need that accuracy. Wilkinson’ has suggested reorthogonalization of each new state with respect to the preceding states as a way of preventing the accumulation of rounding error. From the preceding argument we see that there is nothing to be gained from this and much to be lost if it corrupts the three-term recurrence relation, which is the real reason for stability. The symmetric chain transformation has been implemented without end corrections by Nex’ and is available from him as part of a library of Fortran programs relating to chain models. HUN-,

8.

= uN-IUN-1

CLUSTER SIZE, CHAIN LENGTH, AND

f

CN-luN-2,

BOUNDARY CONDITIONS

There are several simple rules of thumb for the most economical use of the numerical chain transformation. By allocating various amounts of storage, one controls the size of the cluster or number of dimensions M l5

J. H. Wilkinson, “The Algebraic Eigenvalue Problem.” Oxford Univ. Press (Clarendon), London and New York, 1965.

232

ROGER HAYDOCK

included in the Hamiltonian. By allocating various amounts of time for the computation, one controls the length L of the chain model produced. The basic relation between resources is that the time T required is roughly

T = TML,

(8.1)

where the constant z depends on particular Hamiltonian, program, and machine characteristics. Various overheads account for minor deviations from this. A rough measure of the quality of a chain transformation is the agreement between eigenvalues of the chain model and of the original model. For this purpose we suppose that the original model has a continuum of eigenstates spread over a range W of energy. The resolution of the calculation is thus the average spacing of levels in the chain, ( A E ) , divided by W, r = (AE)/W. The number of distinct levels contained in an initial state uo depends on the symmetry of that state and the symmetry of the Hamiltonian. We count how many degrees of freedom S of the Hamiltonian are equivalent under symmetries of the Hamiltonian and uo. For example, suppose the Hamiltonian describes hopping among s orbitals in a cubic lattice. If uo is one of the s orbitals, then the cubic rotations in general make the other orbitals equivalent in groups of 48. For some it may be less, but for most S = 48. Thus for a cluster of M degrees of freedom and ignoring the splitting that could be induced by an asymmetric boundary, which would be an error anyway, the number of eigenstates in uo is about

M o = MIS.

(8.3)

Thus if we solved the cluster exactly as we could not do better than a resolution of about r

=

l/Mo = S / M .

(8.4)

The number of distinct levels represented in a chain model is simply L. There is no point in computing more levels in the chain that M o , because the extra levels will only be ghosts produced either by loss of orthogonality or a weak splitting of degeneracy caused by an asymmetric cluster boundary. Thus if we wish to compute a chain of length L, we must include at least LS dimensions in the cluster. Whether we use a cluster larger than LS depends on the purpose of the chain model. For total energy or other integrated quantities LS is big enough. For differences in energy and other differences in integrated quantities the cluster should be bigger than LS. For smooth

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

233

local densities of states, the clusters should be larger, something like 2LS or more. The reader is referred to Sections 11,18, III,22, and III,24. Thus to obtain a given resolution r with minimal resources,

L

=

l/r,

T = zS/r2.

These are only estimates, but experience has shown them to be, if anything, very conservative. The resolution is usually much better than estimated because of subtleties in the convergence of eigenvalues in chain models. Intuition often misleads the user into using much larger M than necessary for a given L. For continuous systems the energy resolution is the basic number from which most other errors can be estimated. In the case where one seeks convergence of a particular level rather than average convergence over a band, the error estimation is more complicated and is discussed later. It is sometimes suggested that when approximating an infinite periodic system by a finite cluster, imposing periodic boundary conditions on the cluster improves the model. The point of this suggestion is that by taking a finite cluster with a surface, one destroys the translational symmetry of the Hamiltonian. Adding elements to the cluster Hamiltonian that allow an electron to hop out of the surface of a cluster at one point and back into it at another can restore the translational symmetry. This is called a periodic boundary condition as opposed to a cluster boundary condition, where the electron is reflected at the surface. Periodic boundary conditions do not in practice give any improvement in results. One possible reason for this is that even for a Hamiltonian that possesses translation symmetry, if the initial state uo does not also have this symmetry then the chain model will not have it either. The symmetries of the chain model are those both of the Hamiltonian and the initial state uo . We have performed calculations on similar clusters with both periodic and cluster boundary conditions. The results are always slightly different but one cannot be said to be consistently better than the other. Thus we recommend whichever is most convenient. The shape of the cluster boundary affects some calculations, particularly local densities of states. The result of experience is that it is best to make the boundary of the cluster as rough as possible while retaining the symmetry of the Hamiltonian. This avoids the chain model picking up strong interference from coherent reflections from the boundary. It is much better that the boundary should scatter diffusely.

234

ROGER HAYDOCK

9. d ELECTRONS ON AN fcc LATTICE As a numerical example, let us consider the one-electron treatment of the d bands of an fcc metal, say nickel. The states of the system are the five atomiclike orbitals with d symmetry on each atom. They are exactly orthogonal to other orbitals on the same site and approximately orthogonal to orbitals on other sites. The Hamiltonian describes hopping of an electron from atom to atom induced by overlap of orbitals on one site with the potential of a nearest-neighbor site. This hopping is described by the three Slater-Koster l6 parameters: dda when the two orbitals have s symmetry about the axis joining atoms, ddn for p symmetry about that axis, and dd6 for d symmetry about that axis. The elements of H,hopping parameters, between orbitals on nearest-neighbor atoms, are determined by performing the appropriate rotation, obtaining values that in general depend on dda, ddn, and ddd. In this example we approximate the crystal by a cluster of four unit cells of fcc lattice, containing 62 atoms. We wish to model the behavior of a d electron in an infinite crystal, and so we choose initial states in the center of the cluster as having an environment most like the crystal. In an infinite fcc crystal, d orbitals can be chosen so that there are only two that are nonequivalent; one coupling to the E, states and the other coupling to the T2,states. The cluster selected does not possess the full point group symmetry of the fcc lattice. The T2, symmetry is split by the cluster boundary into a doublet and a singlet, while the E, doublet is also split. Table I lists parameters for the first 11 levels for chains initiated by each of the five orbitals on the central atom of the cluster. There are a total of 310 degrees of freedom, which is very modest for a calculation of this type but serves for illustration. The chains have not terminated, but as we shall see later, considerable information can be derived from them. The first few levelsof the chains display the triply degenerate T2,and the doubly degenerate E, symmetries. This is because until the chain states involve atoms on the boundary of the cluster, they have the full point group symmetry of the lattice. At the third level, the boundary splits the two lattice symmetries and this enables us to see the effect of the boundary. There is little to learn from the chain states except that because the Hamiltonian involves only nearestneighbor hopping, the weight of the states moves outward with each level until the states hit the boundary and scatter off it. The behavior of the states is analogous to the free-electron example, although the medium is discrete rather than continuous.

l6

J . C. Slater and G. F. Koster, Phvs. Reo. 94, 1498 (1957).

235

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

TABLE 1. CHAIN PARAMETERS FOR d ELECTRONS IN A N fcc CRYSTAL’”

n

a”

Orbital 0.0 -0.736282B - 02 -0.2856738 - 01 -0.160718B - 01 -0.134218B - 01 -0.112190E - 01 -0.131124E - 01 -0.671958E - 02 -0.1651808 - 01 -0.123618B - 01 0.0

0 1 2 3 4 5 6 7 8 9 10

1 0.100000E 01 0.295643B - 02 0.235335B - 02 0.179840B - 02 0.162690B - 02 0,1562546 - 02 0.174530B - 02 0.153999E - 02 0.195674B - 02 0.173710B - 02 0.168461B - 02

+

Orbital 2 0 3 4 5 6 7 8 9 10

0.0 -0.736282B - 02 -0.285674B - 01 -0.203862B - 01 -0.1347126 - 01 0.9070248 - 03 -0.596882B - 02 -0.4637698 - 02 -0.128520B - 01 0.1280766 - 02 0.0

0 1 2 3 4 5 6 7 8 9 10

0.0 -0.736282B - 02 -0.285673B - 01 -0.203862B - 01 -0.3347128 - 01 0.907012B - 03 -0.596881B - 02 -0.463773B - 02 -0.128520B - 01 0,1280756 - 02 0.0

1

2

b,Z

a“

b,2

+

0.100000E 01 0.2956436 - 02 0,2353358 - 02 0.205787B - 02 0.1529316 - 02 0.164361E - 02 0.155990B - 02 0.131672B - 02 0.134450E - 02 0.160917B - 02 0.158227B - 02

Orbital 0.0 -0.262667E - 02 - 0.833 132E - 02 -0.962267B - 02 -0.5401956 - 02 - 0.472927B - 02 -0,3988386 - 02 -0.611713B - 02 -0.219311B - 02 -0.107605B - 01 0.0

4 0.100000E + 01 0.21 1155E - 02 0.158428B - 02 0.121493E - 02 0.6018818 - 03 0.13 1840E - 02 0,8699938 - 03 0.8958456 - 03 0.854633B - 03 0.126003B - 02 0.125162B - 02

Orbital 5 0.0 -0.262665E -0.833131E -0.l16066B -0.128611B -0,7368588 -0.194090B -0,5327688 0,4531856 -0.123153E 0.0

02 02 01 02 02 02 02 04 01

0.100000E + 01 0.211 155E - 02 0.158428B - 02 0.121205B - 02 0.610331B - 03 0.135928B - 02 0,9850646 - 03 0.829151B -,03 0.9559986 - 03 0.127546E - 02 0.147460B - 02

Orbital 3

a

dd6

”

0.100000E + 01 0.295643E - 02 0,2353358 - 02 0.205787B - 02 0.152932B - 02 0.164361E - 02 0.155990B - 02 0.131672E - 02 0,1344506 - 02 0.160917B - 02 0.158227B - 02

bi is the normalization of the starting orbital. ddu = -0.027784 Ry, ddn = 0.012535 Ry, -0.001554 Ry.

Parameters for chains initiated with each of the five d orbitals on the atom of the center of a cluster of 62 atoms. Hopping parameters were determined from the three Slater-Koster parameters ddu, ddx, and dd6.

236

ROGER HAYDOCK

We leave this example for the moment, but return to it for a demonstration of various techniques of analysis of chain models, discussed in further sections. The demonstration program for this example is available as part of the recursion library. 10. EIGENVALUES AND EIGENSTATES OF A CHAIN

The physics of the chain model reflects the physics of the original system and thus enables us to easily calculate many important quantities. In this section we study the eigenvalues and eigenstates of the original models by means of the chain. We interpret the chain as a set of states of the system that are coupled by the Hamiltonian in a particularly simple way. The Hamiltonian is completely specified on the chain states by the recurrence parameters. Since all the physics is contained in the knowledge of eigenstates of a system, let us derive them for the chain. An eigenstate is some linear combination of the states {uo, uI, . . .}, which when operated on by the Hamiltonian gives a state that is proportional to the original linear combination. Thus we wish to find a linear combination Pnun such that

I."=

(10.1) which is just the time-independent Schrodinger equation. In matrix form, this equation is E - UO -bl

-b1 E -al

0 -bz

0 0

...

0 0

0

-b3

E - ~3

...

...

The standard approach to this matrix equation is to find the zeros of det(E - H),which are the eigenvalues, and then find the eigenvectors. The Jacobi form of the matrix in Eq. (10.2) gives rise to a simple recurrence relation for its determinant. Let An(E)represent the determinant of the first n rows and columns of E - H.Define = 0,

(10.3)

Ao(E) = 1.

(10.4)

237

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

Suppose we know A,(E) and A n - l ( E ) ;let us expand An+, ( E ) in the elements of the ( n 1)th row as follows:

+

det

E

...

E

bn-1 0

-

...

(E

=

- bn-1

- an-2

- a,)

+ b, det

-

det

I

an-l b, E

-0h, 1 - a,

... E - a n P 3 ... - bnP2 E ... 0

bn-2 0 an-2 - b,-, - b,-l E - a,-

... E - a , - 3 - bn-2 0 ... - bn-2 E - a n P 2 0 ... 0 - bn-l - b,

, (10.5)

On the right-hand side, the first determinant is just A,(E), and the second determinant we expand in terms of the elements in the last column, to get

An+ ,(El

=

(E

- a,)

An(E) - b,Z An-

1 m .

(10.6)

In the case where the chain is finite and there is an unsymmetric termination, then by the same argument the final recurrence is

AN(E) = ( E

- a N - l )A N - l ( E ) -

b N - I c N - lA N P 2 ( E ) .

(10.7)

The determinants are polynomials in E with the same order as the subscript. The eigenvalues E, are the zeros of AN(E). Let us now find the eigenstates. Operate with H on the left-hand side of Eq. (10.1) and collect coefficients of u,: m

m

n=O

n=O

1 (a, P?) + b,, ,P f i + b, PF! ,)u, = 1 E, PjP)u,.

(10.8)

We take P?), to be zero. Because the {u,} are linearly independent, we may equate coefficientson the two sides to get a recurrence for the eigenstates,

E,P?)

=

a,P?)

+ b,+lP?il + b,P?l1.

(10.9)

We see that this is adjoint to the chain recurrence, by looking at the case of an unsymmetric termination,

EaP$L2

=

aN-2P$L2

+ cN-,P$Ll + b N - 2 P $ L 3 .

(10.10)

238

ROGER HAYDOCK

From this recurrence we see that the eigenstates are related to the determinants. Since the linear chain can be defined independently of any inner product, there is no natural normalization for the eigenstates. However, it is convenient to make the eigenstates have coefficient one on uo.This makes the coefficientsa set of polynomials in E , and we now include E explicitly in their definition: P-,(E) = 0,

(10.11)

Po@) = 1,

(10.12)

EPn(E) = anPn(E) + bn+IPn+l(E) + bnPn-,(E).

(10.13)

For an unsymmetric termination to a chain of length N ,

+

+

E P N - ~ ( E=) ~ N - ~ P N - ~ (CEN)- ~ P N - ~ ( E~) N - ~ P N - ~ ( (10.14) E). The {P,(E)} are related to the determinants {A,@)} by

A,(E) = bl b2 b3 . . * bnPn(E).

(10.15)

For a finite chain, the condition that for eigenvalues E,, AN(E,,.) is zero, implies that (Ea

- aN-

1IPN-

l(Ea)

- bN-

1PN-2(Ea)

= O*

(10.16)

We can interpret this as saying that if there were an extra level in the chain, no matter what the coupling to it was, P,(E,) would be zero. In future, we write PN(E)as if a bN could be calculated although it is arbitrary. This is much simpler, and bN will always cancel out of formulas. In the case of an infinite-chain model, we determine the eigenvalues differently. The determinants diverge in the sense that the values of the determinant for energies between eigenvalues goes to infinity. However, if we plot the zeros of successive determinants, in the infinite limit, these converge to the eigenvalues of the chain. An equivalent way of looking at the infinite chain is to say eigenvalues are those energies for which the P,(E,) are bounded as n goes to infinity. This is analogous to the condition that the wave functions must be bounded at infinity for a one-dimensional continuous system. Thus for an infinite chain we have the concept of localized states for which P,(Ea) dies off exponentially and extended states for which P,(E,) oscillates with constant amplitude. The zeros of successive {P,(E)) form a Sturm4 sequence in that they separate each other. This and many other useful properties of polynomials are proved by S ~ e g oThis . ~ is an analog of the properties of nodes of the wave functions of Liouville systems. This property facilitates computing eigenvalues for chains. It also implies that there are no degeneracies in the chain.

RECURSIVE SOLUTION OF SCHR~DINGER'SEQUATION

239

The completeness and orthogonality properties of these eigenstates are particularly useful. The polynomials {P,(E)} are orthogonal in the sense that there is a nonnegative weight function no(E), which we later see is the local density of states, such that (10.17)

The expression for completeness is

no(E)

1 Pn(E)Pn(E') = 6(E -

(10.18)

n

where E and E are both eigenvalues. The polynomials define a transform analogous to, but more general than, the Fourier transform. A state

v

1 unun W

=

(10.19)

n=O

and its transform u(E) are related by (10.20)

The inverse transform is (10.21)

The analogy to the Fourier transform is that the chain {u,} replaces Y and the spectrum {E,} replaces k. The transform defined by the polynomials suggests an alternative formulation of a nondegenerate quantum-mechanical model. Beginning with a spectral density no(E), nonzero on {E,}, we define states by their amplitude on each eigenstate, u = u(E),

(10.22)

so that u(E,) is the amplitude of u on the eigenstate with energy E, (remember our system is nondegenerate). Thus

s_, a,

vtSu =

u(E)+u(E)no(E)dE,

HU = Eu(E),

(10.23)

(10.24)

and so the spectral representation of a system can also be transformed to a chain.

,

240

ROGER HAYDOCK

Other useful properties of the eigenstates (polynomials) follow from the Darboux-Christoffel formula4:

This can be proved by applying the recurrence relations to the sum on the left-hand side. In the case as E , goes to E 2 , we get

where the prime indicates the derivative with respect to E . The subject of orthogonal polynomials provides a wealth of theorems that can be usefully applied to the eigenstates of the chain model. If we wish to evaluate the sum (10.27)

we may do so by operating the recurrence backwards. Each polynomial P,(E) is related to the determinant A,(E) by Eq. (10.15). We reverse the order of all rows and columns so that each determinant starts with row n and column n in the upper left-hand corner and works down to row zero and column zero in the lower right-hand corner. Now let

and continue according to

until ao(E) is found. This method is quick, accurate, and economical. 1 1. IRREGULAR SOLUTIONS ON A CHAIN Just as with a second-order differential equation, there are solutions to the chain model that are irregular at uo . In the previous section we discussed the eigenstates, which are the regular solutions of the three-term recurrence. For a given E , a solution satisfying Rn + 1( E ) = C(E - an)Rn(E) - b n Rn - 1(E)l/bn + 1

(11.1)

RECURSIVE SOLUTION OF SCHRODINGER'S EQUATION

24 1

is determined by the choice of two values, R o ( E ) and R , ( E ) . Hence therq are two linearly independent solutions. The eigenstates {P,(E)} result from the initial conditions

P d E ) = 1,

(11.2)

Pl(E) = ( E - aoM1.

(11.3)

A linearly independent solution { Q,(E)}, called the irregular polynomials or polynomials of the second kind, results from the initial conditions

QoW

= 0,

(11.4)

Qi(E) = 1-

(11.5)

Alternatively we can think of them as the regular polynomials for a chain with uo deleted. The {Qn(E)}are also orthogonal polynomials, but with a different weight function from the {Pn(E)}.They have the Sturm sequence property and they have the added separation property that the zeros of P,(E) and Q,(E) alternate along the E axis. Qn(E)is of order n - 1 and so has one fewer zero than Pn(E). In analogy with the second-order differential equation, any pair of independent solutions possess a relationship similar to that of the Wronskian theorem. Writing the recurrences for Pn+,(E) and Qn+,(E) and multiplying by Q,(E) and Pn(E),respectively, gives an Qn(E)Pn(E) + bn + 1 €?,(Elf', + 1 (El

+ bn Qn(E)Pn-

EPn(E)Qn(E) = anPn(E)Qn(E) + bn+ ,Pn(E)Qn+ , ( E )

+ hnPn(E)Qn-

EQn(E)Pn(E)

=

(11.6) (11.7)

Subtraction results in the identity bn+ 1 CQn + l(E)Pn(E) - Qn(E)Pn + I (El1 =

C Q n ( W n - ,(E)

-

Qn

- l(@Pn(QI.

( 1 1.8)

We define the analog of the Wronskian W,C{Pm(E)I, {Qm(E)>I

=

b n + I C Q n + ,(E)Pn(E) - Q n ( W n + 1(E)19 (11.9)

and by chaining the equations like (11.8) together, we show that it is a constant,

KC{Pm(E)} {Qrn(E)11 = b 1. 9

(11.10)

An elegant proof of the Darboux-Christoffel identity follows from doing the same manipulations with { P n ( E l ) }and { P n ( E 2 ) } .

242

ROGER HAYDOCK

The associated functions, {A?,(E)}are useful in calculating changes in the local density of states (cf. Section 111,20). If the {P,(E)} are orthonormal with respect to no@), then the associated function A?,(E) is the Hilbert transforms of P,(E) weighted by no(E),

}

%(E) = Re{ J~mCP.(Eo)no(Eo)/(E - E0)1 dE0 .

(11.11)

One can demonstrate that the associated functions are also solutions of the three-term recurrence for n > 1, although in general they are not polynomials. They are orthogonal to the { P , ( E ) } in the sense that where no(E) > 0,

2 Pn(E)2,(E)

=

(11.12)

0.

n=O

12. THECONSTANT CHAIN The chain model with constant parameters turns up again and again and has a number of useful properties. It is defined by a, = a,

(12.1)

b, = b.

(12.2)

Thus the recurrence relation for its polynomials is

EP,(E)

=

aP,(E)

+ bP,+ i ( E ) + bP,- i(E).

(12.3)

This is satisfied by, P,(E)

=

sin{(n

+ 1) cos-’[(E - a)/(2b)]}/sin{cos-’[(E

-

a)/(2b)]}. (12.4)

These are Chebyshev polynomials of the second kind” and the first few are (12.5) (12.6) (12.7) The polynomials are proportional to determinants of various length chains, and so we can find the eigenvalues of a chain of length n by determining the zeros of P,(E),

E:) = a

+ 2b cos[rnx/(n + l)],

(12.8)

where m = 1,2,. . . ,n. As expected, P,(E) has n zeros and they separate those of P,- , ( E ) and P,, , ( E ) according to the Sturm property. They are located

RECURSIVE SOLUTION OF SCHRODINGER'S EQUATION

243

between a - 2b and a + 2b, and in the limit as n goes to infinity they fill the interval a - 2b < E < a + 2b. The eigenstates are sinusoidal sums of the states in the chain. In the band of eigenvalues, 8 = cos-'[(E - a)/(2b)] is always real and between 0 and n. Thus for a given E, the eigenstate is proportional to

1 sin[(n + i)e]u,. W

(12.9)

n=O

This is just the Bloch state for a one-dimensional semi-infinite crystal. Clearly the constant chain is the one related to the Fourier transform. To illustrate some of the other properties of orthogonal polynomials we can determine the spectral function no(E). The orthogonality condition is J-Wm[sin(n

Jo

+ 1)O sin(m + l)O/[(sin2 O)]no(E) dE = dnm.

(12.10)

The eigenvalues correspond to 8 between 0 and n, and so changing the variable to 8 gives

+

sin(n + l)Osin(m

+ 1)Ono(E)(2b/sin8)dO = dnm.

From the orthogonality of sin(n + 1)0 and sin(m

+ l)O, we see that (12.12)

no(E) = sin cos-'[(E - a)/(2b)]/(2b), or no@) =

(12.11)

+

a - 2b < E < a 2b, ( 12.13) E < a - 2b,E > a 2b.

[4bZ - (E - a)2]'/2/(4b2),

+

This is in fact the local density of states at site zero on the chain. The constant chain is important because it is related to an isolated band. The band extends from a - 2b to a 2b and includes all phase factors between 0 and n. Models with bounded spectra transform to chains that are asymptotically constant.

+

13. THECHAIN MODELOF A FREE ELECTRON Let us return to the chain model for a free electron and derive its eigenvalues and eigenstates. The chain model was defined by a, = ~ ( 2 n+ 3/2),

(13.1)

b, = E(n)(n

(13.2)

+ 1/2)'12.

244

ROGER HAYDOCK

The polynomials determined from these chain parameters are related to the Laguerre polynomials,'2 P,(E)

= (-

P,(E)

=

1)"[(1/2)!n!/(n + 1/2)!]'/2L,'/2(E/&)

(- 1)"[(1/2)!n!(n

+ 1/2)!]""

n

x

(13.3)

1 (-l)"(E/&)'"/[(n

-

m)!(m

m=O

+ 1/2)!m!].

(13.4)

They are orthogonal, with spectral weight

(13.5) The zeros of P,(E) are all at E > 0 and as n goes to infinity they fill the positive real-energy axis. Using the chain eigenstates and the chain states, we can write expressions for the eigenstates of the free electron, (13.6) Substituting the expressions for P,(E) and u, in terms of Laguerre polynomials, +&) = (A/&)""

m

1 [(1/2)! n ! / ( n + 1/2)!]L;/2[(Ar)2]

n=O

x L,"~(E/E) e~p[-Ar)~/2].

(13.7)

The summation can be done and gives

+E(r) = (1/2) ( ~ / & ) ~ / ~ ( 1 / 2~) !X ~ [ E / ( ~ E ) ~ J , , , [ A ~ ( E / E ) ~ ~ ~'12.I / [ A ~ ( E / E ) ~ (13.8) Thus the chain model for the free electron leads to the spherical Bessel functions as it should. 14. SYMMETRY Taking advantage of symmetry in a Hamiltonian can save time and storage in computation of a chain model or considerably simplify the analytic transformation to a chain. The only difficulty is that a full analysis of symmetries may be time consuming and complicated in itself. In this

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

245

section we briefly go through the application of symmetry to the transformation to a chain model, and show with an example how it can simplify a problem. Symmetry of the Hamiltonian can lead to degeneracy in its spectrum. Let H be a Hamiltonian with eigenvalues { E , } , possibly degenerate, corresponding to an orthonormal choice of eigenvectors {w,}. If in terms of the {w,}, T is the unitary representation of some symmetry of H , then H and T commute, T H = HT.

(14.1)

Hw, = E,w,

(14.2)

H T w , = E,Tw,.

(14.3)

Hence, implies Thus T w , is also an eigenstate of H with the same energy as w,. Unless exp(i+,)w, = T w , for all a, T w , is linearly independent of w, for some 1and E , is thus degenerate. Conversely, suppose E , = E,, for 1# p, and Z is the identity; we can construct a unitary transformation T = Z - (w, - w,)(w~ - w , ) ~ .

(14.4)

T commutes with H and hence represents a symmetry of H . It is not always possible to give a simple physical interpretation to transformations that commute with H , but we treat them all in the same way. The chain model has no symmetries. This follows from the nondegenerate nature of the chain’s eigenvalue spectrum. For the purposes of this section we may add to our other pictures of the chain model that of a model containing states of only a single symmetry, all the other symmetries having been “factored out.” The most efficient use of the chain transformation results when the starting state belongs to an irreducible representation of the symmetries of H . When this is the case, the only states coupled to uo and hence appearing in the chain model are those belonging to the same row of the same irreducible representation of the symmetries of H . This means that only such states need be retained, and so in computations the same resolution can be achieved with much less storage and time because the dimensionality of the matrix of H is reduced. As an example consider a one-electron treatment of a square array of hydrogen atoms. We construct a tight-binding model with nearest-neighbor hopping only. The basis states are {I,,,”}, representing an s orbital at lattice

246

ROGER HAYDOCK

site (m,n) in Cartesian coordinates with a lattice spacing of one. If we take lo,, for the starting state and wish to obtain a chain that agrees to N levels with the result for the infinite lattice, then we must include in the cluster for calculation 2 N 2 + 2N + 1 sites so that we can always take N - 1 hops on the lattice without passing outside the cluster. However, loo belongs to the one-dimensional representation of the square lattice point group (the cluster has no translation symmetry because it is finite). Thus lo, only couples to the combinations of states that are symmetric under square rotations. We construct them as follows:

,

N - 1 of the Ig), Counting the number of basis states, we need lo,,, and N 2 / 4 - N / 2 - 1 of the if N is even or N 2 / 4 - N / 2 - 314 of the lE)n if N is odd. Thus making use of the symmetry, we need only about 118 as much storage to do the same calculation. There is a corresponding saving in time. The chains produced by starting states corresponding to different rows of the same irreducible representation are similar. Suppose T is a unitary representation of a symmetry that takes u, into Tu,: then because T commutes with H,the nth state u, in the chain generated by u, goes to Tun, the nth state in the chain generated by Tu,, and the chain parameters are identical. This means that a chain model need be calculated only once for each irreducible representation. Note also that states belonging to different irreducible representations or different rows of the same irreducible representation do not interact.' This allows us to construct starting states with special properties. If we wish to calculate chains for each row of several different irreducible representations, we can combine this into one calculation by constructing a grand starting state that is a linear combination of one starting state from each irreducible representation weighted by the square root of the dimensionality of that irreducible representation. Because none of these components interact, we show later that various quantities calculated for the single grand chain are equal to the sums of such quantities calculated for the individual chains.

'' J. M. Gallagher. Ph.D. Thesis, Uiversity of Cambridge, 1978 (unpublished).

RECURSIVE SOLUTION OF SCHR~DINGER'SEQUATION

uO

247

"1 "2 u3 FIG.2. A graph representation of a ladder model.

While convenient in the sense of requiring only one chain transformation, the above method of combining several chains representing different symmetries into one chain can reduce resolution. If the eigenvalues of the different symmetries are at very different energies, the band width of the sum will be much greater than any of the separate chains. Thus, for the same length of chain, the resolution can be greatly reduced in the sum. For good resolution it is best to divide a calculation into many chains, each covering a small range of energy. 15. VECTORCHAINS There is a generalization of the chain model where the chain states become multidimensional subspaces and the parameters become operators between subspaces. The simplest example of this generalization is a ladder model whose graph representation is shown in Fig. 2. The pairs of orthogonal states u, and v, define a two-dimensional space orthogonal to the spaces defined by the other pairs of states. The {a,}, {b,}, and {c,} determine how the Hamiltonian acts on each subspace. In terms of matrices these generalized chains correspond to matrices that are block-tridiagonal. The recurrence relations for polynomials become recurrences for vectors of polynomials, with matrix coefficients. The properties of the ordinary chains and polynomials mostly generalize, but care must be taken to account for the noncommutivity of matrices. The greatest benefit of this generalization is that some physical relations can be made explicit in a vector chain when they are obscure in the usual chain model. This is paid for in computations by the complexity of the matrix algebra. 16.

COMPUTATION OF

EIGENVALUES AND EIGENSTATES

In those cases where the eigenvalues and eigenstate of the chain cannot be found analytically, they can be computed numerically. In this section we discuss the computation of these quantities that, in principle, completely

248

ROGER HAYDOCK

describe the physics of a system. The fundamental computation is that of the polynomials associated with the chain. They are defined by P-,(E) = 0, P,(E) = 1, and

C(E - an)Pn(E) - bn Pn - 1(E)l/bn + 1 * If the chain is of length N with an unsymmetric termination, then

(16.1)

C(E - Q N - 2 ) P N - m - bN-ZPN-3(E)I/CN-l.

(16.2)

Pn +

PN-l(E)

=

=

The last polynomial is calculated with an arbitrary b, since we are only interested in its zeros, PN(E)

=

[(E

-

uN-

l)PN-

- bN-

lPN-

2(E)1/bN.

(16.3)

The N zeros of P N ( E )are the N eigenvalues of the chain. We wish to determine the zeros E, of P N ( E )and the values of the polynomials P,(E,), which are the coefficients of u, in the eigenstates of the chain. In the construction of the chain we were able to demonstrate that cumulative loss of orthogonality of the chain states did not make the model inaccurate. Unfortunately, there is no corresponding result for the polynomials. Rounding errors in the computation of the recurrence Eq. (16.1) do accumulate in the higher polynomials. This is balanced by the ease with which the polynomials can be computed to high precision. Whereas transformation to a chain can press the limits of time and storage available, the corresponding eigenvalues and eigenvectors can be calculated in multiple precision with only a fraction of the resources. The accumulation of rounding error in computation can be avoided by using matrix methods. The three-term recurrence defines the matrix in Eq. (4.3). Whereas the full H was too big to use matrix methods for finding eigenvalues and eigenvectors, the J matrix requires much less storage. Matrix methods operate on all states simultaneously, and so there is no accumulation of rounding error. We do not discuss these further here, but refer the reader to Wilkinson.’ Rounding error does not affect the Sturm sequence or separation property of the roots of the polynomial. A method for finding the roots of P N ( E )is to calculate P N ( E )from the recurrence, keeping track of the number of sign changes between P,(E) and P,, , ( E ) as the recurrence proceeds. From the Sturm sequence property, the number of sign changes is the number of roots of PN greater than E. One thus homes in on the roots without danger of skipping any. The accuracy desired in the roots of P N ( E ) must be carried at each stage of the calculation. The larger rounding error occurs when P,(E) and P,_,(E) differ greatly in magnitude. This error is minimized by renormalizing P,(E) and P,- , ( E ) to order one before computing P,+ This is possible because the normalization of a polynomial does not change

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

5

0

10

15

249

”

FIG.3. Schematic behavior of the roots of Pn(.E),{ E g ’ } ,as a function of n.

its roots. Of course, within the recurrence relation, normalizations must be the same. As a function of n, the extremal eigenvalues of a chain converge fastest and the central ones slowest. If we calculate the roots {E:)} of P,(E) we find that the largest and smallest roots vary little with n after the first few levels. This convergence moves in toward the center of the spectrum from the extrema. A qualitative explanation for the phenomena is that (H - a,)u, contains proportionally more of the eigenvalues furthest from a,. This is counterbalanced by subtracting b,u,- For each eigenvalue there is a value of n for which the subtractions overcome the multiplications and the root corresponding to that eigenvalues ceases to change. This is shown schematically in Fig. 3. Once the eigenvalues have been found, the eigenstates follow simply. If the eigenstatesare eventually to be computed, it is wise,during the transformation to the chain, to keep the components of u, in backing storage as they are computed. The sum w(E.2) =

c P,(Ea)u,

N-

1

(16.4)

fl=O

can be accumulated using the recurrence in Eqs. (10.28) and (10.29). If the {u,} are not stored, then we can compute them by applying the same trick as in Eqs. (10.28) and (10.29) to the {ufl}.Suppose N- 1

wo

=

1 c,u,.

(16.5)

fl=O

Define WN = W N + 1 =

(16.6)

0,

and proceed according to

w, = cnuo

+ C(H - a n + l ) W f l + I

-

bn+2Wfl+*I/bn+1

(16.7)

250

ROGER HAYDOCK

until wo has been calculated. Thus the computation is once up the chain at an eigenvalue and once back down to get an eigenstate. Fortran subroutines for calculating the polynomials and the roots are available.'

II. Green Functions and the Local Density of States 17. GREEN FUNCTIONS OF A FINITE CHAIN

A chain model leads to a simple form for the Green function, an economical way of summarizing information about a system's eigenvalues and eigenstates, where their enumeration would be too lengthy to be of any use. Formally the local Green function or matrix element of the Greenian between two states u and v is GO(E) = v'S(E - H)-'u,

(17.1)

where E is a complexenergy, not equal to an eigenvalue of H.Go@) describes the propagation of the system from state u to state v at energy E , hence local because it is relative to two specific states of the system. This is clearer in a spectral decomposition of the Green function. Let {w,} be an orthonormal set of right-hand eigenstates of H with eigenvalues { E , } . The spectral decomposition is

1 w,(E - E,)-'W;

Go(E) = v'S

[ a

I

SU.

(17.2)

Here we see that the Green function has a simple pole, at each eigenvalue, with residue equal to the product of the projection of that eigenstate on u and the complex conjugate of its projection on v. A contour integral of the Green function picks up the eigenstates in the contour that couple to u and v. Any element of the greenian can be expressed as a sum ofdiagonal elements, which have a particularly simple relation to the chain model. For u and v normalized, four independent linear combinations of u and v are

uo = u

+ v,

(17.3)

u

- v,

(17.4)

+ iv,

(17.5)

zo = u - iv.

(17.6)

vo

=

wo = u

25 1

RECURSIVE SOLUTION OF SCHR~DINGER'SEQUATION

In terms of these four normalized states, vV(E - H)-'u

=

[u&S(E- H)-'uo - V&S(E- H)-'vo -i(w&S(E - H)-'w~- z&S(E - H ) - ' z , J ] / ~ . (17.7)

The diagonal elements of the greenian have a number of useful properties as functions of complex E . Consider the general diagonal Green function (17.8)

Go(E) = uXS(E - N)-'u,.

As E goes to infinity in the complex plane in any direction, Go@) goes to zero like 1/E. Second, we can relate G,(E) in the upper and lower halves of the complex E plane by Go(E)t = u&(Et- Ht)-'Stu0 = u&S(Et- H ) - ' u ~= Go(Et). (17.9) This depends on the hermiticity of the model, which implies that S and SH are self-adjoint. Third is the Herglotz3 property. From the spectral decomposition of the diagonal Green function, Eq. (17.2), we see that the residue of each eigenvalue in G,(E) is positive. This forces the imaginary part of Go@) to be negative for E in the upper half of the complex E plane. By hermiticity, the imaginary part of Go(E) is positive for E in the lower halfplane. Thus the poles of G(E), which occur at eigenvalues of H,are real and the zeros of Go@)are also real and separate the poles along the real E axis. A function with the above properties can be written as a continued fraction," Go(E) = 1/(E - bf/{E x ( E - a d . . .I}),

-

b$/[E -

-

.. . - b i -

1/

(17.10)

where we show that the {a,} and {b,} are the same parameters as for the corresponding chain model. The continued fraction is the functional composition of fractional linear transformations (17.1 1) For a, real and b,: positive, one can see that for E in the upper half of the complex plane and G,, in the lower half-plane, G, is also in the lower halfplane. If G,+ , ( E ) goes like 1/E as E goes to infinity, so does G,(E), and implies that

G,+ 1(Et) = G,+ l(EY

G,(Et) l9

=

G,(E)t.

(17.12a) (17.12b)

H. S. Wall, "Analytic Theory of Continued Fractions." Van Nostrand-Reinhold, Princeton, New Jersey, 1948.

252

ROGER HAYDOCK

FIG.4. Schematic behavior of Re[bi+ lGn+l(E)] and E - a, for real E.

If G , + , ( E ) has m poles on the real axis, then G,(E) has m + 1 poles along the real axis. A function with only simple poles and satisfying the analytic properties of a diagonal element of the Green function looks schematically like that in Fig. 4.If we draw a line E - a,, we see that the line and curve have one more intersection than G , , ,(E) has poles and that these intersections, which are the poles of G,(E), separate the poles of G n + l ( E ) .Thus the fractional linear transformation Eq. (17.1 l), transforms one “Green function” into another. Using a finite-chain model we can calculate directly the diagonal element of the greenian or local Green function corresponding to uo.Equation (17.8) defines a diagonal element of the greenian. Equation (4.3) defines the matrix for H in terms of the orthogonal {u,,}. Hence the Go(E) we seek is just

... - b2 ... - b l E - al - b3 ... - b2 E - a 2 0 0 - b3 E - a3 ... 0

E-UO

-bl

0

0 0

-1

( 17.13a) 0

Using the standard formula for an inverse element of a matrix, the zerozero element of the inverse is the ratio of the cofactor of E - a. to the determinant. Let us define the determinant of the matrix with the first n rows and columns deleted as D,(E). (Note that these are defined differently from the A,(E) in Section 1,lO.) Then

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

253

Use the same expansion of determinants as in Eq. (10.5) to get Go(E) = Dl(E)/C(E - ao>Dl(E>- b:D,(E)I.

(1 7 . 1 3 ~ )

Now divide top and bottom by D , ( E ) to get Go@)

=

1/CE - a0 - b:D,(E)/D,(E)l.

(17.13d)

We can define D,(E)/D,(E) as G , ( E ) and continue the procedure, producing a continued fraction and demonstrating that the continued fraction parameters are the same as the chain parameters. Note that the above expansion of Go(E)depends on the orthogonality of {u,}. We discuss later the errors introduced by nonorthogonality in computations. The local density of states for uo is related to the residues of Go@) at its singularities. The Green function Go@) has simple poles on the real axis, at energies {El},and each pole has positive residue p l . The spectrum of Go(E) or local density of states or orbital uo is just (17.14) 1

From the spectral decomposition of the Green function, Eq. (17.2), we see that pa is just the magnitude squared of the overlap of eigenstate w, on uo. Thus the local density of states is just the total density of eigenvalues weighted by the magnitude squared of the overlap between each eigenstate and uo. If we consider the imaginary part of Go(E)along the lines parallel to the real E axis, but slightly above it, we see that no@)

=

1 lim - - Im[Go(E + ie)],

(17.15)

71

&+O

where E goes to zero from above. We abbreviate this by 1 no@) = - Im[Go(E)l, 71

(17.16)

where the limit is implied. Each pole on the real E axis produces a delta distribution in the limit. The real part of Go(E)along the real E axis is the Hilbert transform of the density of states, Re[Go(E)]

r m

=

J

[no(E’)/(E - E’)] dE‘.

(17.17)

-m

It is related physically to the resistance of the system to excitation at other than the eigenenergies. Another way of putting it is that the real part of Go@) describes the response of the system to being driven at a given energy.

254

ROGER HAYDOCK

The Green function can also be expressed as the ratio of two polynomials. The determinants D,(E) are polynomials in E of degree N - n, where N is the number of states in the chain model. Hence from Eq. (17.13b), we see that the two polynomials are Do@) and D,(E). Do@) is proportional to PN(E) defined in Section 1,lO. If we relate D l ( E ) , a chain with uo deleted, we can define polynomials { Q,(E)} as follows : (17.18) (17.19)

Q, + 1(El = C(E - QJQAE) - 6, Q. - 1( E ) l / h+ 1 9

(17.20)

which are the irregular polynomials of Section IJl. In this case Q,(E) is of degree n - 1. Relating these polynomials to the determinants, we get

Go(E) = QN(E)/C~~P&)I.

(17.21)

[Both P,(E) and QN(E)are calculated using the same arbitrary b N ,which cancels out of the above expression.] The zeros of P,(E) are the poles of Go(E) and the zeros of QN(E)are the zeros of Go@), which separate the poles. This is also a Pade approximant.20 We can now relate the polynomial expression for the Green function to the eigenstates of the chain model. As shown in Section IJO, the eigenstates of the finite-chain model are N- 1 wa

a

1Pn(EA*

n=O

(17.22)

When the {u,} are orthonormal, we normalize the eigenstates to get (17.23) Now, substitute this into Eq. (17.2), the spectral decomposition of the Green function : N-

1

(17.24)

Using Eq. (10.26), this simplifiesto Go(E) =

1 CbNPN- I(Ea)PXEa)(E - Ea)I-’*

(17.25)

U

*’

G . A. Baker, Jr. and J. L. Gammel, ed., “The Padt Approximant in Theoretical Physics.” Academic Press, New York. 1970.

255

RECURSIVESOLUTION OF SCHR~DINGER’SEQUATION

This formula gives the normalized spectral weight of each eigenstate on uo and so we can write the density of states directly as

6(E - E a ) / [ b N P N -

=

l(Ea)PN(Ea)I.

(1 7.26)

a

The zeros of P,(E) are E a , and so we obtain the identity QdEa)

=

b~C

~ PN N - l(Ea)I -

‘9

(17.27)

which, of course, is not true for E , other than the zeros of PN(E).In general, this identity relates the values of Q,(E) and P,- ,(E) at the zeros of P,(E). This completes a survey of the theory of Green functions on a finite chain. More information about orthogonal polynomials, continued fractions, and Pade approximants may be obtained from Shohat and Tamarkiq3 Wall,’ Akhiezer,’ S ~ e g oand , ~ Baker.” 18. GREEN FUNCTIONS OF AN

INFINITE CHAIN AND

THEIR INTEGRALS

The continuous eigenvalue spectra that emerge in solid state physics characterize infinite-chain models. In the previous section we established a number of properties of Green functions for finite chains. In this section we extend those results to infinite chains and discuss approximations to Green functions of infinite chains. The analytic structure of the Greenian for a finite system was simple in that it consisted entirely of first-order poles and zeros. As we go to the infinite limit a number of complications occur. For an infinite chain, the simple poles may either remain discrete or become dense in certain ranges of energy, which we loosely call bands. The simplest form of a band is when Go@) acquires a pair of branch points, say E , and E B . The analytic properties of Go(@ force us to draw the branch cut on the real energy axis between EA and E B . What was the residue of poles becomes a discontinuity in the imaginary part of G,(E) across this cut. The local density of states no@) ceases to consist of isolated delta distributions and in the band they coalesce into a continuous spectrum. The band can be more complicated with extra branch points corresponding to Van Hove singularities in the density of states. There can also be essential singularities on the real E axis. In infinite random chains, the local density of states at energies in the localized regions is discontinuous, but dense. The variety of analytic structure possible in the Green function of an infinite chain poses a difficult problem for approximations. The obvious approach is to consider the infinite chain as a limit of finite chains. Let

Gg’(E, Z ) = 1 / { E -

- b:/[E -

-

. . . - b,2/(E - a, - Z ) . .

a]},

(18.1)

256

ROGER HAYDOCK

where Z represents the remainder of the fraction. From Eq. (17.21) this may be written as G6"'(E, Z ) = C(E - an - Z)Qn(E) - bnQn-i(E)I/{C(E - an - Z ) x Pfl(E) - ~ n ~ n - l ( ~ ) l ~ l l . (18.2)

We now use the definition of P,+ ,(E) and Qn+,(E) by means of an arbitrary b,+ to give Gg)(E,Z ) = [Q. + 1( E ) - ZQn(E)/bn+ 1 I/Cb 1(Pn + 1 ( E ) - ZPn(E)/bn+ 1 )I* (18.3)

For a given value of E , we may look upon this as a fractional linear transformation of Z. If E is in the lower half-plane, then Z may be anywhere in the upper half-plane, and the transformation gives the resulting value for G,(E). The real Z axis bounds the upper half-plane and the transformation takes it into a circle. Thus the points within the circle are the possible values of G,(E). Applying some complex analysis, we find that the center of this circle of values for G,(E) is GP'(E)

=

(Qn(E)/Pn(E) + Pn(Et) C Q n + i(E)Pn(E) - Qn(E)Pn + 1(E)I/ x {Pn(E)CP,+ l(E)P,(Et) - pfl+l(~t>Pn(E>ll)/bl, (18.4)

and its radius is

@"'(El

=

I CQn+ i(E)Pn(E) - Qn(E)Pn+i(E)I/CPn+i(E>Pn(Et) - p,+ l(Et)P,(E)ll/h. (18.5)

We may simplify these expressions by considering some general properties of a three-term recurrence. The {P,(E)} and {Q,(E)} are two linearly independent solutions of the three-term recurrence, but with different initial namely, P,(E) is one and Q o ( E ) is zero. In analogy with a secondorder differential equation (of which the three-term recurrence is a generalization), the analog of the Wronskian (see Section A l l ) w,(E)

=

CQ, + ,(E)Pn(E) - Qn(E)Pn+1(E)Ibn+1

(1 8.6)

is independent of n. Hence bn + 1 C Q n + 1(E)Pn(E) - Qn(E)Pn+ 1 (Ell

=

bi.

(1 8.7)

This allows us to write G?'(E)

=

IQ,(E)/Cb iPn(E)I x

+ Pn(Et)/Cbn + iPn(E)I

CPn+l(E)Pn(Et) - P n ( J W n + l(Et>ll

(18.8)

and =

Cbn+iIPn+i(E)P,(Et)- ~n(E)f',+l(Et)ll-l.

(18.9)

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

257

From the Darboux-Christoffel identity, Eq. (10.25), we see that m=O

Remember that we have taken E in the lower half-plane. The convergence of finite-chain approximations to the infinite chain depends on the divergence of Eq. (18.10), which, as n increases, must shrink the circle containing possible values of G,(E) to a point. That this takes place for complex E is the condition for the chain defining a unique Green function, but for real energies the error is always infinite. Nonconvergence of finite approximants to Green functions for real energies afflicts every computational method. The reason for this is that an arbitrarily small perturbation can shift an eigenvalue, which, in moving a pole of the Green function, causes an infinite change in the Green function near that pole. Physically, this is as it should be, since adding atoms to a cluster does slightly shift all the levels of the cluster. The solution to this problem is to examine what we mean by the convergence of the physical properties of a finite system to those of an infinite system. The physical properties are all averages of some sort. At the very least they are densities and they usually involve averages over a number of eigenstates. It is in the “mean” that the computation of Green functions converges for real energies. In precise terms, the convergence that takes place for real energies is that of indefinite integrals of the local density of states

J-J(E)no(E) dE, Eo

F(E0) =

(18.11)

where f ( E ) is a function whose derivatives are all less than or equal to zero over the energies where no(E) is nonzero. no(E)is physical in the sense that there is a finite energy below which no(E) is zero. Given {P,(E), . . . , P,(E)} determined by { u o , . . . , a,- I } and {bl, . . . , bn},choose so that

A,

=

Eo - bflPfl- I(EO)/Pfl(EO)

Rfl+l(E)= ( E - Afl)Pfl(E)- bflP,-l(E)

(18.12) (18.13)

has a zero at E o . Let { E l , . . . , E n } be its other zeros. Now construct polynomials of order 2n, P + ( E )and P-(E), with the following properties: P + ( E m )= {,(Em)’

P - ( E m ) = {,(Em)’ P’+(Em) = f’(Em),

for for for for for

Em I Eo, Em> E,, Em < E,, Em 2 E , , Em# E,.

(18.14) (18.15) (18.16) ( 18.1 7) (18.18)

258

ROGER HAYDOCK

FIG.5. The bounding properties of P + ( E )and P _ ( E ) onf(E), shown schematically.

+

This can be done because a polynomial of degree 2n can satisfy 2n 1 constraints. For E IE , , P + ( E ) and P - ( E ) , respectively, bound , f ( E )from above and below. That Pi@) are bounds depends on the derivatives of f(E)being nonpositive. For E > E , the polynomials bound zero above and below. This is shown schematically in Fig. 5. Thus upper and lower bounds Fi(Eo) of F(E,) are given by

Fi(E0)

J-/ * ( E ) n , ( E ) Eo

=

dE.

(18.19)

This integral can be done exactly2*20*21 and gives (1 8.20)

(18.21) where Prn =

Pn(Ern)Pk + 1(Ern)]-

Cbn+ 1

' 9

(18.22)

which are the same as the residues of G,(E) in Eq. (17.25). The difference between the bounds is pof(E,). Using Eq. (10.26),the limit of the DarbouxChristoffel identity, we get (18.23) If Eo is the energy of a bound state, p o has a finite limit as n increases and the value of the integral depends on whether we include E , or not. Otherwise p o goes to zero as n increases and the bounds on F(E,) come together. C. M. M. Nex, J . Phys. A 1 1 , 653 (1978).

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

259

From the discussion of these two kinds of bounds on the Green function a condition on the chain model emerges. This condition is expressed in terms of (18.24) For the approximants to the infinite Green function to converge, which we shall interpret as a condition for the chain model to represent a physical system, p(E) must be zero everywhere in the complex E plane except at a countable number of points on the E axis. There is no equivalent criterion in terms of the chain model parameters, but for a chain to model a physical system it is sufficient that m

1 l/bn = co

(18.25)

n=l

or (18.26) Thus we have established the way in which finite-chain models approximate infinite-chain models. We have found that the Green function for the finite chain approaches that of the infinite chain everywhere except on the real energy axis. For real energies the density of states of finite chains approaches no limit, yet its indefinite integral over a restricted class of functions does. This says that for real energies we can determine the integrated density of states, but not its derivative, a curious, but as we have shown, a physically sensible conclusion. 19. EXAMPLES OF GREEN FUNCTIONS OF INFINITE CHAIN MODELS

Chain models whose Green functions can be found analytically give insight into the physics of various systems. The physics that can be investigated this way includes the behavior of the Green function inside a band, near various kinds of Van Hove singularities, at various kinds of band edges, near a bound state, when a bound state is near various band edges, for a resonance in a band, and for a resonance near various kinds of Van Hove singularities. The analytically solvable chains are those related to the classical orthogonal polynomials such as Chebyshev polynomials, Legendre polynomials, and Laguerre polynomials. Each of these sets of polynomials satisfies an infinite three-term recurrence that determines for what chain the polynomials are eigenstates. These polynomials also satisfy second-order

260

ROGER HAYDOCK

differential equations that, by the theory of Sturm-Liouville systems, determine their weight functions. The first example is that of the constant linear chain. By choosing the appropriate zero of energy, we may make the (a,} all zero. The continued fraction form of the Green function, Eq. (17.10), gives Go(E) = 1/{E - b2/[E - b2/(E - . . .)I}.

(19.1)

In this case we can proceed to find the Green function by solving the continued fraction directly. This method works whenever some solution of part of the fraction can be related to the solution of the whole. We see that (19.2)

Go(E) = [ E - b2Go(E)]-’,

which is a solvable self-consistent equation. In general such a self-consistent equation for a continued fraction introduces spurious roots, as this does : G,(E)

=

[ E f ( E 2 - 4b2)”2]/(2b2).

(19.3)

We can eliminate the false solution by means of the analytic properties of Go(E). For E in the upper half of the complex plane and the usual definition of Z”’ where the cut is made just below the negative real Z axis, we get [ E - ( E 2 - 4b2)’/2]/(2b2), GdE) = [ E ( E 2 - 4b2)”2]/(2b2),

+

Re{E} > 0, Re{E} < 0.

(19.4)

Hermiticity determines the solution in the lower half of the complex plane. Go@) has square root branch points at f2b, and a branch cut along the real E axis between them. Figure 6 shows the real and imaginary parts of Go(E) just above the real axis. Note that at the band edges both the real and imaginary parts of Go come in with infinite slope.

FIG.6 . The real and imaginary parts of G,(E

+ ie) for a constant chain.

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

26 1

We can investigate resonances and bound states by allowing a, and bl to be chosen arbitrarily, but making the rest of the chain constant. This sort of model is characteristic of an atom adsorbed on a surface or a defect in the bulk of a solid. Remaining in the upper half of the complex E plane, G,(E)

=

{ E - a, - b:[E f ( E 2 - 4b2)’/2]/(2b2)}-1,

Re{E} >< 0. (19.5)

If we divide top and bottom by b:, we see that the denominator is the difference between ( E - ao)/b: and the Green function for the constant chain, = {b:C(E - ao)/b: -

GO(J!~II-~.

(19.6)

Using Fig. 6, we can see the shape of G,(E). First of all, G,(E) only has an imaginary part where Go@) does. Second, if I b , I is greater than 2 1 b 1, there must be a pole of G,(E) outside the band, indicating a bound state. If 1 bl I is less than 2 1 b 1, there may or may not be a bound state, depending on a, . In the band, the local density of states is enhanced where ( E - ao)/b: is between zero and 2 Re{G,(E)}. A resonance occurs when bl is very small and a, is in the band. The density of states is depressed everywhere except the region where E is near a,. This gives rise to a Lorentzian density of states near E = a,. When a, is outside the band, there is a bound state that is shifted from a. by interaction with the band. Physically these cases correspond to various interactions and matching conditions between two systems represented by u, and the constant chain. When bl is large and a, is in the band, we have strong coupling between well: matched systems. This gives rise to bound and antibound states outside the band. As a, moves away from the band, the matching deteriorates until there is just one bound state, only slightly perturbed from a, by the band. In the weak-coupling case, bl is small and when a, is in the band, the systems are well matched. The resonance represents the system slowly leaking away from u, due to coupling with the band. In the weak-coupling and badmatching case, the bound state has very little weight in the constant part of the chain. From the infinite slope of Re{G,(E)} near the band edge, we can see that in the intermediate-coupling, intermediate-matching case, the interaction can change greatly over a small range of a. and b,. The Legendre polynomials, when orthonormalized, satisfy the recurrence

EP,(E)

=

(n

+ 1)[4(n + 1)’ - 1]-1’2Pn+1(E)+ n(4n2 - l)-1’2Pn-l(E).

(19.7)

They are orthogonal with respect to the weight function, nom =

{i;

-1IEI1, otherwise.

(19.8)

262

ROGER HAYDOCK

I----

FIG.7. The real and imaginary parts of G,(E polynomials.

+ k ) for the chain model related to Legendre

In this case it is easiest to find the Green function from the spectrum Go(E) =

+1

2 -

which gives Go@)

=

In[@

( E - Eo)- d E o ,

(19.9)

+ l)/(E - 1)]/2.

(19.10)

Figure 7 shows the real and imaginary parts of this Green function for E just above the real axis. By changing the zero and scale of energy, the band can be given a desired width and position. This chain gives rise to logarithmic singularities at the band edges and a constant density of states, characteristic of two-dimensional systems. The real part of Go(E)goes to infinity and the imaginary part has a discontinuity at the band edges. This chain has a fundamentally different behavior when coupled to another system. Suppose we consider a chain with a. and bl arbitrary and then for n > 0,

b,

1)-

a,

=

0,

G,(E)

=

(b:C(E - a0MJ: - GO(E)l}-

= n(4n2 -

(19.11)

'1'.

As before, l.

(1 9.12)

We see that for b: positive, G,(E) must always have two poles outside the band, corresponding to bound and antibound states. A system with this kind of spectrum couples strongly to anything in the sense of creating bound states. The structure of the interaction is more complicated in this example in that resonances in the band also occur for any coupling.

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

263

The last example is a generalization of the chain model for a free electron. When normalized, the Laquerre polynomials satisfy the recurrence EPr’(E) = (2n + 1 + a)P?’(E) - (n + 1)1’2(n+ 1 + a)’” x Prjl(E) - n’/2(n + L X ) ” ~ P ~ ? ~ ( E ) ,

(19.13)

for LX > - 1 . The polynomials are orthogonal with respect to the weight function N‘,“)(E)=

:{

exp{ -E},

E > 0, otherwise.

(19.14)

Again we find the Green function by means of the spectral theorem, GO@) =

r*

J

0

E$ exp{ -Eo} (E - Eo)- dE,.

(19.15)

Evaluation of this integral involves the incomplete gamma function T(u, u). The result is Go(E) = - r ( l

+ a)T( -LX,

- E ) ( -E)” exp{ - E } .

(19.16)

Depending on LX,the band edge at zero energy takes on the character of the branch point associated with the power LX.The density of states varies like E“ and the real part of the Green function behaves like its Hilbert transform. This Green function represents a band with no upper bound, such as for free particles. In fact, the infinite extent of the band has little effect. The behavior of bound states and resonances near the band edge is determined by the band edge rather than any other parts of the band. 111. Perturbation Theories, Energy Differences, and Computing the Density of States

20.

THEORY FOR THE LOCAL DENSITY OF STATES, VANHOVESINGULARITIES, AND CHAIN ASYMPTOTICS PERTURBATION

We can associate various features of the local density of states of an infinite system with different parts of its chain model. Since all finite chains result in discrete local densities of states, the existence of a continuum and the features related to it must depend on the infinite nature of the chain. These features are the essential singularities of the spectrum, the band edges, and other Van Hove singularities. The asymptotic (large-n) behavior of the chain model determines the distribution of spectral points. In contrast, the low-n behavior of the chain determines the distribution of weight among the

264

ROGER HAYDOCK

spectral points in the local density of states. This makes physical sense in that the shape of a wave function is determined by the local environment, the low-n part of the chain, but if it is extended its energy is determined by the boundary conditions, the asymptotic part of the chain. We can approach the problem of an infinite-chain model as that of a cluster embedded in an effectivemedium, the asymptotic chain. Suppose we have an infinite chain with parameters {a,} and {b,}. We define

GN(E) 1/{E - a, - b,$+,/[E - a N + 1 - b,$+,/(E - u N + ~ - * * * ) I } (20.1) . Now, using Eq. (18.3), we can relate Go(E), the local Green function, to GN(E) by GO(E) = CQN(E)- bN

QN-

l(E)GN(E)l/{blCPN(E) - bN p N - l(E)GN(E)l}. (20.2)

The function G,(E) specifies the influence of the effective medium on the cluster, namely the finite-chain model {ao, . . . ,u ~ - ~{ b} l, ,. . ., b N } .If N is sufficiently large, we get a good approximation by replacing G N ( E ) by F,(E), an asymptotic approximation. The larger N , the more exactly FN(E) approximates GN(E). The asymptotic chain is determined entirely by the essential singularities of the spectrum and vice versa. As we showed in Section 11,17, a finite composition of fractional linear transformations, Eq. (20.2), can only add simple poles to a spectrum. Therefore, any essential singularities in Go(E) must also be in GN(E).On the other hand, the fractional linear transformations can smoothly change the distribution of spectral weight so that the only aspects of GN(E) that are preserved by the fractional linear transformations are the essential singularities. Of course, a smooth match between the two systems prevents any localized states from forming at the interface. Information about the asymptotic chain comes from the original quantummechanical model and from analysis of the finite part of the chain. The essential singularities of a spectrum can often be found from the long-range (asymptotic) properties of the original quantum-mechanical model, such as order or average order. This can determine the asymptotic chain completely. Also one can discover regularities in the computed parameters, which are characteristic of the asymptotic chain. The relation between small variations in the chain and small variations in the local density of states allows us to relate essential singularities to the asymptotic chain. Suppose we have an infinite-chain model with parameters {a,} and {b,}. The regular polynomials {P,(E)} are defined by Eqs. (10.1 1)-(10.13) and the irregular polynomials { Q J E ) } are defined in Eqs. (17.22)-(17.24). Let us introduce a small change in one of the chain parameters. In the first instance, let us change a, by 60,.

265

RECURSIVE SOLUTION OF SCHRODINGER'S EQUATION

The effect of the change in the chain is to mix the old chain's regular and irregular solutions at a particular energy to form the new solutions. This is just the scattering effect of ban. Let {p,(E)> be the regular polynomial solutions of the new recurrence. In terms of the polynomials for the unperturbed chain, the new polynomials are unchanged up to p,(E) because 6a, does not appear in the recurrence for either. For the next polynomial, pfl

+ 1(El = Pfl+ l(E) - 6%P,(E)Ib,+

(20.3)

1.

Except for the above recurrence, the equations for the new polynomials are the same as for the old. Since at a given energy, there are only two linearly independent solutions to a three-term recurrence, for m greater than n, we can write p,,,(E) = a(E)P,(E) + B(E)Q,(E). (20.4)

Taking n + 1 and n + 2 for rn in the above equation, we get two equations for a(E) and P(E) that when solved give

4E)=

Cpn+ 1 (E)Qn+ z(E) - p n + z(E)Qn+ I (E)I/CPn + 1 (@Qn

+ z(E)

-f ' , + z(E)Qn+

P(E) = CPn + 1 (E)pn+ z(E) - Pn + z(E)pn + 1 (E)lICPn+ 1 (QQn + 2(E) - Pn + z(E)Qn+ 1(Ell*

(20.5) (20.6)

The Wronskian theorem, Section IJl, for the {P,(E)} and {Q,(E)} eliminates the denominators and leaves

4 E ) = h+z[Pn+i(E)Qn+2(E) - pn+z(E)Q,+i(E)l/bi,

(20.7)

P(E) = b, + z CP,+ 1(EP, + 2 ( E ) - P, + Z(E)Pn + 1 (E)l/b 1 *

(20.8)

The Wronskian theorem can be used on each of these to obtain 19

(20.9)

B(E) = b, + 1 C P f l W f+l ,(El - P, + 1(E)P,(E)lIbl.

(20.10)

4 E ) = b, + 1CpAE>Qn+

- p n + 1( E ) Q A J ~ I / . ~

Now substitute Eq. (20.3) for p,+ l(E), F,(E) is unchanged,

a(E) = 1

B(E)

+ %PAE)Q,(E)/bi,

(20.1 1) (20.12)

= - ha, P , ( E ) P f l ( W b l .

We can perform the same manipulations for a variation in b, of 6b,. This gives the following results: a(E) = 1 + 6bnCPn- l(E)Qn(E) + bnPn(E)Qn- l(E)(bn + 6 b J 'l/b19

P(E) = -6b,[P,-

,(E)P,(E)

+ bnP,(E)P,-l(E)(b,+ 6bJ

(20.13)

' ] / b i . (20.14)

266

ROGER HAYDOCK

The change in the local density of states can be calculated from the change in parameters. Let no@) be the local density of states for the unperturbed chain and n(E) the total density of states of the chain, namely, (20.15) where { E A } are the eigenvalues of the chain. Then, from Eqs. (17.14) and ( 1 8.23), we get 1

(20.16)

The perturbed density of states is (20.17) In the continuous part of the spectrum, the sum of squares of the polynomials must diverge, and so we may replace p m ( E )for all m greater than n, using Eq. (20.4). As long as n is finite,

We may use the same total density of states for the continuous part of the perturbed and unperturbed chains, because the two-chains' parameters go to the same limit. Define (20.19) the local density of states for the chain with its first state removed, and row

=

Re[E o:;J

no(E)I(E - t o ) ] .

(20.20)

The {Qm(E)}are related to the associated functions {bZm(E)} of Eq. (11.11): Qm(E) = bl rdE)Pm(E) - 9 m ( E ) 7

(20.21)

and so we can use the orthogonality of {P,(E)} amd {bZm(E)}to get 2o(E) = {.(E)2Cno(E)1- + B(E)2Cnl(E)1+ 2b 1 a(E)B(E)ro(E)Cno(E)1- 1- l .

(20.22)

As an example, let us consider the constant chain with a defect at orbital zero. Let us take the band to be of width 4b and centered at zero energy.

RECURSIVE SOLUTION OF SCHRODINGER'S EQUATION

267

This gives the following results for a(E) and p(E) if the defect is represented by nonzero ao: .

a(E) = 1,

(20.23)

B(E) = - ao/b.

(20.24)

From Eq. (19.4) we get no(E) = n,(E) = (4b2 - E2)1/2/(2nb2)

ro(E) =

i

[ E - ( E 2 - 4b2)]'I2/(2b2), E > 2b, E/(2b2), -2b 5 E 5 2b, [ E ( E 2 - 4b2)]ly2/(2b2), E -2b.

-=

+

(20.25 ) (20.26)

Substituting these results in Eq. (20.22),and considering only the energies in the band, we get iio(E) = (4b2 - E 2 ) 1 1 2 [ ( 2 ~ + ) ( b2 ~ 2- a o E ) ] - ' ,

-2b IE I 2 b . (20.27) The reader can verify that this agrees with the result from Eq. (19.5).This form of the new local density of states shows more clearly the effect of the perturbation. When E and a. are of the same sign the new local density of states is enhanced, while it is depressed when they are of opposite sign. Equation (20.22) applies only to energies within a band. We see that the denominator of Eq. (20.29) can go to zero outside the band, where the unperturbed density of states is zero. These zeros represent either bound states or virtual bound states (simple poles for real E on the second sheet of the Green functions), but in this form we cannot tell which, although it is clear from the local Green function, Eq. (19.9, that the bound state occurs if I a. I > 2 1 b I and the virtual bound state otherwise. When we have a large number of small deviations from a solvable chain model we construct a theory analogous to that of multiple scattering. Let us consider only the first-order effectof each scatterer. The first-order change in the density of states 6no(E)is related to the first-order change of a(E) from one, 6a(E), and of p ( E ) from zero, Sfi(E), &o(E) = -2[6a(E)

+ b,ro(E) Gp(E)]no(E).

(20.28)

Relating 6a(E) and Sfl(E) to {dun} and {6bn},we get m

a(E) = 1 {ha,Pn(E)Qn(E)+ dbn + 1CPn(E)Qn+ 1 (El n=O

+ pn + 1(E)Qn(E)l l/bi

c C%Pn(E)Pn(E) + 2 '6bn+,Pn(E)Pn+l(J5)I/b2.

9

(20.29)

m

P(E) = -

n=O

(20.30)

268

ROGER HAYDOCK

We can now combine these equations with Eq. (20.21) to give an expression for 6no(E)in terms of the {P,(E)} and the associated functions {A?,@)}:

(20.31) Equation (20.31) relates the change in the local density of states to a sum of linearly independent functions whose coefficients are the perturbations in the chain model, thus enabling one to either compute the change in local density of states from a change in the chain or vice versa. A special case of this relationship was proposed by Hodges.22Equation (20.31) corresponds to the first approximation to the T matrix of an array of scatterers as the sum of the T matrices of the individual scatterers. It does not include successive scatterings by different scatterers although it does properly account for multiple scattering from the same scatterer if we distinguish between 6b,u,- ,u; and Gb,u,u,t- in the perturbed-chain Hamiltonian. The convergence of this approximation is the same as the first-order convergence of multiple-scattering theory. Singularities can arise due to the interaction of the scatterers and since this is a second-order or higher effect, it is not included in Eq. (20.31). Experience shows it to happen at band edges. The asymptotic form of {a,} and {b,} determine the external points of the continuous spectrum in the local density of states. The energies for which the {P,(E)} remain bounded as n increases but at which P,(E)’ diverges form the continuous part of the local density of states. The recurrence for the P,(E) indicates when this can happen:

cz=o

pn + I (El = C(E/bn+ I 1 - (an/bn + 1)IPn(E) - (b,/b,+ 1)p, - i ( ~ 9 . (20.32)

We see that for a continuous spectrum to exist, the ratio bn/bn+lmust go to one. The values of (E/b,+l) - (a,,/b,,+l) for which P,(E) is bounded are those between minus two and plus two. Some examples are (1) if the {a,} and {b,} go to constants a and b for large n, then the band edges are a f 2b; (2) if the parameters go like an and bn for large n, then the band edges are the limits of (a - 2b)n and (a 2b)n as n goes to infinity. In the case of Laguerre polynomials, these band limits are zero and plus infinity. Equation (20.3 1) enables us to relate essential singularities within the continuum to the asymptotic behavior of the chain model. The position and nature of essential singularities within the bands do not affect the limits of a, and b, because they do not alter the upper and lower limits of the continuous spectrum. These internal singularities do alter the way in which

+

22

C. L. Hodges, J . Phys. Left.38. 187 (1977).

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

269

the chain parameters approach their limits. By going sufficiently far down a chain, the ha, and 66, that give rise to internal singularities can be made as small as necessary. Thus in order to study internal singularities, we need only consider the asymptotic behavior of the sum in Eq. (20.31). From the other point of view, any change in the local density of states involving a particular internal singularity gives rise to an expansion in {f‘m(E)%n(E)}, {Pm(E)%- 1(E)l, and {pm(E)%n+ whose asymptotic properties are unique to the position and nature of that singularity. A particular asymptotic form of ban and 6b, combines with the {P,(E)} and {2?,,,(E)} to give a particular singularity in the change in the local density of states, and conversely a change in the local density of states with a particular singularity has a particular asymptotic form for its expansion in products of polynomials and associated functions. For studying singularities by means of Eq. (20.31) it is often convenient to use asymptotic forms for the { P , ( E ) ) and {2?,,,(E)}.This works because only the asymptotic part of the sum matters anyway. Because the classical orthogonal polynomials all satisfy second-order differential equations, asymptotic forms are particularly easy to find. Even for nonclassical recurrences, asymptotic solutions can be found. A simple generalization of the WKB approximation produces asymptotic forms for a three-term recurrence. Let us think of P,(E) in Eq. (20.31) as a traveling-wave solution to the recurrence, not as a polynomial that is a standing-wave solution. Suppose the coefficients of P,(E) and P,- , ( E ) in Eq. (20.32) vary slowly with n. Let f n ( E ) be the ratio between the wave at successive sites on the chain, Jl(E) = Pfl(E)/P,- l ( J 3 In terms of this ratio, Eq. (20.32) becomes

Jl+ ,(a= C(E - an)/b,+ 11 - b,/Cbfl+ 1 Jl(E)I.

(20.33) (20.34)

If ( E - a,)/b,+ and b,/b,+ vary slowly with n, then so must f,(E). In fact, in the constant chain Jl(E) is a constant in n. Thus for a zeroth approximation we take f:O’(E)

=

f2

I@),

(20.35)

and solve Eq. (20.34) to get

f:Y1(E)

=

{E

- a,

f [ ( E - u , ) ~- 4b,b,+1]”2}/(2b,+1). (20.36)

In the spirit of WKB, we can improve this by correcting the relationship between f,(E) and f,+, ( E ) so that the (m 1)th approximation is

+

f F+W) = f:Y’(E)f:m’(E)/f:T ,(El

(20.37)

270

ROGER HAYDOCK

in terms of the mth approximation. The final result is

f!,?<’)(E)

=

E - a, f {C(E - a,)2

-

46,6,+lf!,~l(E)/f~m’(E)I”Z}/(2~n+l). (20.38)

For large n this gives the form of the traveling-wave solutions. The plus and minus roots in the approximation apply to solutions traveling up or down the chain. The polynomials are the standing waves, which are the symmetric and antisymmetric combinations of the traveling waves. IN A BANDOF FINITE WIDTH 21. SINGULARITIES

In this section we apply the methods of Section III,20 to show how in a finite-width band the singularities are related to Fourier series. The constant chain is convenientlydescribed by an angular variable 8. Suppose the diagonal chain parameters are all a and the off-diagonal ones are all 6. Define

’

8 = cos - [(E - a)/(26)]

(21.1)

so that as we move in energy across the band, 8 varies from zero at the top to 7c at the bottom. In terms of 8, P,(E)

=

sin[(n

+ 1)8]/sin 8.

(21.2)

In order to study singularities, we need the associated functions defined in Eq. ( l l . l l ) , which turn out to be

g , ( ~= )c o s ~ +( ~lie].

In terms of 8,

no@) = sin 8/(7cb). Applying the product formulas for sine and cosine, we get m

6no(E) =

n=O

{6a, sin[(2n

+ 2)8] + 2 6b,+ , sin[(2n + 3)8]}/(nb)’.

(2 1.3) (21.4)

(21.5)

The change in local density of states is given by the Fourier series with coefficients 6a, and 26b,, and the changes in a chain are just the Fourier coefficients of the change in local density of states, hu, = 26’ IOn6n(E)sinC(2n

6b,+

=

6’

rn

J

0

+ 2)8] do,

6n(E) sinC(2n + 3)8] do.

(21.6) (21.7)

22. COMPUTING THE LOCALDENSITY OF STATES The density of states is the first picture a reader looks at to understand a theory of the electronic structure of some complicated material. A band

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

27 1

structure, if one exists, is usually too complicated to take in at a glance, but the density of states gives a summary of theoretical electronic structure as well as directly relating to various spectroscopic properties of a material. For this reason, computation of the density of states from a chain model is the single most important result of the recursive solution of the Schrodinger equation. There are two approaches to the density-of-states computation. The first is to compute a sufficiently long chain model so that most of the levels of a cluster are represented. The local density of states is then smoothed by some method so that it resembles the local density of states of an infinite system. The second method is to use a much larger cluster to calculate the initial part of the chain, which is the same for the cluster as for the infinite system. A theory of the asymptotic behavior of the chain permits computation of the local density of states. A discrete local density of states results directly from the eigenvalues and eigenstates of a finite chain. Section I,16 discussed the computation of these quantities and Eq. (17.18), where pa is computed from the normalized eigenstates, expresses the discrete density of states as a sum of delta distributions. The resulting spectrum is very difficult to compare with the density of states from a band structure unless the delta distributions are very closely spaced indeed. The main method for smoothing such a discrete spectrum uses the convergence of the indefinite integral of the local density of states.21 Equations (18.20)-(18.22) give bounds to the indefinite integral of certain functions over the local density of states. Let us take f(E)in those equations to be one. The smoothed density of states is then taken to be the derivative of the mean of the bounds, (22.1) where the quantities on the right-hand side are defined in Section II,18. The differentiation indicated here can be performed either numerically or analytically. Numerical differentiation can achieve some further smoothing because there is a degree of instability in the analytic form for ns(Eo): ~ s ( E o= )

c‘

En,5 Eo

{pi+I(Em)Qn+ ,(E)Pb(Em) - [Pi+1(Em)I2Qn(Em)

+ P:,+ I (Em)Q:n

+ 1 (Em)Pn(Em) - Pi + 1 (Em)Qn+ l(Em)Pn(Em)}/ (22.2)

c‘

in the notation of Eqs. (18.12)-(18.19), and with the sum including only half the term for Em = E o . The first and second derivatives of the polynomials Pi(E), Qi(E), and Pi(E) can all be calculated recursively by differentiating their recurrence relations.

272

ROGER HAYDOCK

Because the bounds form an envelope for the possible integrated local densities of states, any actual function would vary more rapidly and thus have more structure. Hence this approximation has the minimum structure consistent with the finite chain from which it comes. In general, the chain one might calculate from ns(Eo) would be similar but not in exact agreement with any part of the chain it approximates. However, as the bounds converge, then the chain calculated from ns(E) converges to the initial chain. Another important point about this approximation is that in general the approximation

d (22.3) + tPOf(E0)l dE0 in the notation of Section 11,18, is different from f(Eo)ns(Eo). The above method of smoothing a discrete local density of states has been implemented in Fortran programs in the recursion library.' There are two versions ; one using analytic differentiation and the other using numerical differentiation in the library. Once E , has been chosen, the bounds on the integral at El, . . . , En can also be calculated at the same time. There is a program that takes advantage of this. While it is very fast, the local density of states is mostly plotted at points the user does not choose, which can be inconvenient. Figure 8(a)-(f) shows the sum of the local densities of states, in the above approximation using analytic differentiation, for the d orbitals on an atom at the center of various sized fcc clusters and various length chain models. This example was used in Section I,9. The method and parameters are the same as explained there, but the cluster sizes and chain lengths have been altered. Figure 8 shows how certain features of the local density of states develop as more atoms are included in the cluster. It also shows that for many levels of a small cluster ns(E) still contains resolved levels, delta distributions. Note that the resolved levels are toward the edge of the band in accordance with the property that resolution of individual levels moves in from the extremes of the spectrum as the chain length increases. The size and shape of the cluster affect how much the smoothed local density of states resembles the band structure result. The number of orbitals that are not equivalent under symmetries of the Hamiltonian needs to be at least twice the number of levels to be calculated in the chain. Also the cluster should have no symmetries beyond those of the Hamiltonian for the infinite system. This means choosing clusters shaped like cubes rather than spheres that have additional symmetry. Another way of smoothing a discrete spectrum is to take for the local density of states, the imaginary part of the Green function at points off the real energy axis. If we look at the imaginary part of 1/(E - ie) for real E fs(Eo)ns(Eo)

= -CF-(Eo)

A,, ,-1.1 ,.--‘.1 ,/I;, ,A, ,A,, 273

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

n(E)

n(E)

20r

,

- 0.1

0

E

0.1

0

E

0.1

10 LEVELS

62 ATOMS

5 LEVELS

62 ATOMS

- 0.1

(b)

(a)

n(E)

,

,F

- 0.1

0

- 0.1

0.1

5 LEVELS

364 ATOMS

0 0.1 1687 ATOMS 10 LEVELS

E

(d)

(C)

n(E)

n(E)

,~

- 0.1

0

0.1

2456 ATOMS 15 LEVELS

(el

-0 1

0

2 4 5 6 ATOMS

0.1

E

20 LEVELS

(f)

FIG.8. Smoothed local densities of states of d orbitals of fcc clusters of various sizes and for various lengths of chain model.

274

ROGER HAYDOCK

x-poles of GJE) 0 - z e r o s of GJEI FIG.9. A contour giving a smooth approximant to the density of states of G,(E).

with a fixed positive E, we find that it is peaked about E = 0 and has width of about E. If we find a contour such as that shown in Fig. 9, which stays away from the real energy axis by a distance about the same as the spacing of poles in G,(E), the imaginary part of Go@) along this contour'will give a smooth approximation to the local density of states,

no@) x

1

-

71

G [ E - ie(E)].

(22.4)

Because a contour can be distorted without changing integrals over it, the above approximation gives integrals of any entire function over n,(E) exactly. This method has been used by Jones23 for densities of states of defects in silicon and germanium. There are several other methods of smoothing a discrete spectrum including making a histogram and convolution. A histogram of the total spectral weight in an interval loses some information because not only is the spectral weight significant, but so is the number of levels within the interval-more levels indicating a flatter spectrum than one level containing all the weight. Convolution of the discrete spectrum with a smooth function is an arbitrary, though effective, procedure. One would often like to change the width of the convoluting function for various parts of the spectrum, which brings this method back to the complex contour idea. In all these smoothing procedures the basic problem remains that, quite independent of the chain model, there is no convergence of the local density of states for a cluster to that of an infinite system. This fact indicates that we should only plot the bounds on the integrated local density of states, which 23

R. Jones, Philos. Ma,q. 35,57 (1977).

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

275

does converge. However, this loses all the advantages of physical insight gained by some attempt at a local density of states. An approach that avoids this difficulty is the use of an asymptotic approximation for the chain. If, by use of a finite cluster, one can compute a number of chain parameters that are the same as for the infinite system, then an asymptotic approximation with the correct essential singularities gives an approximation to the local density of states, which converges pointwise. This can only be done if the mat.rix H in the model is sparse, i.e., it has only a finite number of nonzero elements in each row and column. If H is not sparse, then no part of the chain model for a cluster agrees with that of the infinite system. Experience shows that for complicated systems, the rough idea of the spectrum one gets from a quick and simple smoothing of a discrete spectrum is better than a long and difficult analysis of the many Van Hove singularities. One has to use clusters as big as is practical for the latter approach, while the former approach requires smaller clusters but more time. However, asymptotic analysis could become the superior method with improvements in theory. 23.

THEAVERAGE (TOTAL) DENSITY OF STATES

The average density of states is the volume density of the energy density of eigenstates of the whole system. In finite systems the average density of states is proportional to the total density of states, while in infinite systems the total density of states is infinite. In terms of a complete (thus not of a single-chain model) orthonormal basis for the system {Z,}, and its volume V, the average density of states is

n ( ~=)

CC na(E>I/V,

(23.1)

where n,(E) is the local density of states for Z,. In a crystal we can write the average density of states as a finite sum. Two orbitals Z, and Z, are equivalent if there is some symmetry of the Hamiltonian that takes Z, into Z,. This implies that equivalent orbitals have identical local densities of states. Let {Z,} be the set of nonequivalent orbitals. Each occurs with a number density pa given by the number of such orbitals per unit volume. The average density of states is thus

n(E) =

1Ppn,(E).

(23.2)

B

For crystals the average density of states can be computed either by computing the n,(E) separately and then adding them or by constructing a chain model for the average density of states. If the spectra for different Z, are widely separated in energy, the first method gives the best resolution.

276

ROGER HAYDOCK

If the spectra for different Z, are all in the same energy range, then the symmetry properties of Section I,14 enable one to compute the average density of states from only a few starting orbitals. For systems with an infinite number of significant,nonequivalent orbitals, the average density of states can only be approximated by averages over finite subsets of the nonequivalent orbitals. A crystalline surface has an infinite number of nonequivalent sites, but those other than bulk sites form an infinitesimal fraction and hence do not contribute significantly to the average. It is in random systems that the average density of states must be approximately computed by averaging over a finite number of orbitals. Experience shows that this average converges very rapidly. For a system described by a nonorthogonal basis, the calculation of the average density of states involves the sum of local densities of states over orbitals different from the basis orbitals. The problem is that if the basis set {4a(r)}of Section I,3 has been chosen so that the diagonal elements of S, defined in Eq. (3.2), are all one, then because the orbitals overlap each other, adding up local densities of states overcounts each state to give too large a total. A solution to this problem is to consider a complete set of orbitals {Za}such that ZJZ, = d@. (23.3) This is not any physical orthonormality, but is a kind of formal orthonormality with respect to the nonorthogonal basis {q5a(r)}.The. average density of states for a finite system is n ( ~ =)

1

(-.)-I

Im[(S-'Z,)+S(E - H ) - ' Z a ] / V .

(23.4)

a

Approximation of S - ' was discussed in Section I,3. The individual terms of the sum are not local densities of states because they are not of the symmetric form of Eq. (17.8). As a result the off-diagonal Greenian element in each term must be calculated by resolving it into a sum of diagonal elements as in Eq. (17.7). The same ideas as for orthogonal bases apply to equivalent orbitals and the average density of states in infinite crystals and disordered systems. 24. DIFFERENCES IN ENERGY AND

OTHER INTEGRATED

QUANTITIES

The Fermi level as a function of electron density and the contribution to the total energy of occupied states are both examples of integrated quantities that may change if the model Hamiltonian is changed slightly. Under certain conditions such as some structural transitions where atomic volume remains constant, the change in these quantities comes mainly from integrals of the independent-electron density of states (see Chapter 4).

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

277

These differences can be accurately computed by the methods described in previous sections, even when the difference of the integrated quantities are smaller than the absolute error in the integrals, Eqs. (18.20)-(18.21). The absolute error limits of an integral apply to comparison of different clusters. For instance, if we were comparing the Fermi level and electronic contribution to structural energy of hcp and fcc nickel, the absolute errors would apply to the difference between an hcp cluster of 100 atoms and a differently shaped fcc cluster of lo00 atoms. However, if we compare integrated quantities between the same sized clusters with similar overall size and shape, then the errors due to cluster size and shape cancel out and the residue is due to the different stucture. Another way of looking at this is that we can perform the calculation exactly for the finite cluster-the errors bounds for finite clusters are relative to the infinite system-and as the clusters get bigger, the differences must converge to their values for the infinite systems. More will be said about errors in small differences in Section III,27 on perturbation theory. 25. BINDINGENERGY

The energy of interaction or binding between two subsystems, say an adsorbate and a substrate or an impurity and a host lattice, can be calculated recursively. Conceptually, this binding energy is the difference between the total energy with the interactions turned on and the total energy with them turned off. Even for a finite system, this is the difference between two large numbers, but for infinite systems the total energy is infinite and thus it is impossible to calculate binding energies straightforwardly. This is related to Section III,24, where we discussed changes in total energy that were reflected in a change in the average density of states. A change in the average density of states gives rise to a change in energy density (binding energy per atom or volume), so that the change in total energy for an infinite system is infinite. Here we are concerned with changes in total energy that are finite and thus zero when averaged over an infinite system. The binding energy can be calculated directly rather than as the difference of large numbers.24 Under certain circumstances, the binding energy comes mainly from the difference in the sum of energies of all the occupied independent-electron levels, (25.1) occupied

24

occupied

T. L. Einstein and J. R. Schrieffer, Phys. Rev. B 7,3629 (1973).

278

ROGER HAYDOCK

where {E&)}are the independent electron energies in the interacting case, and {ELD)}the energies in the decoupled case. The factor of 2 comes from double occupancy of each level. These energies are the eigenvalues of two hamiltonians HI for the interacting system and HDfor the decoupled system. By decoupled, we mean that HDcan be written as H o a,, where the orbitals acted on by H, and H, are orthogonal. To calculate AU directly, we group the terms so that the cancellation of the infinite parts of the total energy takes place term by term, leaving a single convergent sum,

+

AU

=

2

1

(25.2)

(E:) - ELD)),

occupied

where we have indexed the eigenvalues so as to achieve this. We now calculate AU by a discrete procedure for turning on the intera c t i o n ~ Let . ~ ~us choose some basis (not necessaily orthogonal) {Z,} for the smaller of the decoupled subsystems, say R, . Define intermediate systems with Hamiltonians H, R,, where H, is HI with the rows and columns corresponding to Z,, 1, Z n + 2 ,. . . excluded and R, is H o with the rows and column corresponding to Z,, Z,, . . . , Z, excluded. Although we have not yet mentioned them, we assume that there SD,S o , and so forth related to each Hamiltonian. are overlap matrices SI, The overlap matrices for H,and R,,are defined similarly to the Hamiltonians. Thus H o and 17, are the Hamiltonians for the two parts of the decoupled system. If 8, has N dimensions, then 8, is zero and HN is the same as HI. Physically, the discrete turning on of the interactions is that of successively transferring each degree of freedom from R0to H o . Let { E r ) } be the eigenvalues of 17, and ( E f ) }be the eigenvalues of H,.In terms of these, the change in total energy of transferring Z, from one system to the other is

+

so,

AU,

=

2

[1

(Er)- Er-')) -

occupied

occupied

(Et-') - Er))].

(25.3)

If in this process electrons are transferred between the decoupled subsystems, then there may be terms in AU, in addition to those in Eq. (25.3), which we ignore at the moment. Thus the total energy difference is

AU

N

=

AU,,

(25.4)

n= 1

for 8, of N dimensions. The change in total energy AU, is determined by the local Green functions for Z, in H,-and 8,- First we establish the relationships between { E t ) } and { E r - ' ) } and between {Et)} and {EF-')}. H, is H,-l with one extra 25

N. R. Burke, Surf. Sci. 58, 349 (1976).

RECURSIVE SOLUTION OF SCHR~DINGER'SEQUATION

279

row and column, just as 8,- is i7, with one extra row and column. From the separation theorem of Section I,lO, or the general theory of eigenvalues of matrices, we find that adding an extra row and column to a matrix leaves some eigenvalues unchanged and changes those to which the extra orbital Z, couples. Considering only those that change, the new and old separate each other. The states that couple to the new orbital are just those spanned by the which start with Z,. In fact, orbitals in the chain models of H , and from Eq. (17.13b) we see that {EF)} are the poles and { E F - ' ) } the zeros of ZJS,(E - H,)-'Z,, and similarly {@-')} are the poles and {Eg)} the zeros ofZ,tS,-,(E - B,-J~Z,. If several of the orbitals {Z,} have the same symmetry, the calculation of AU can be simplified still further. Suppose Z, and Z,+, have the same symmetry. Then they both couple to the same states of H,and H , - and the poles of Z,'S,(E - H,)-'Z, and ZL+,S,(E - i7,)-'Zfl+, are, respectively, the same as the zeros of Z ~ + , S , + , ( E- Hn+l)-lZn+l and Z,'S,-,(E - H , - I ) - 'Z,. Thus in any sequence of Z,, . . . , Z, of states of the same symmetry, we need only calculate the Green functions for the first and last orbitals in the sequence to find how the eigenvalues have shifted. However, if each orbital has a different symmetry, it couples to different states and we need to calculate each AU, . In the case of infinite systems, the eigenvalues can be distributed continuously, and so we must write the AU, as integrals over Green functions. Let

a,-,,

G$'(E)

=

ZLS,(E - H , ) - ' Z , ,

Gb")(E)= Z,tS,-,(E

-

(25.5)

H,-,)-'Z,.

(25.6)

For a Green function G,(E) containing only simple poles and zeros,

d dE

-ln[bG,(E)]

( E - EZ)-' -

=

zeros

1 ( E - Ep)-',

(25.7)

poles

where { E , } are the zeros and {Ep}the poles of G,(E). Suppose that we only consider the ground states of the various systems, and that we have a common Fermi level EF for both subsystems. If the subsystem H, is very large compared to a,, then the Fermi level will differ infinitesimally between the coupled and decoupled cases and we can define an energy contour CF that contains the poles and zeros of G$)(E) and G((;)(E)corresponding to the occupied levels. Using Eqs. (25.3) and (25.5-25.7) we write, AU,

=

{$

CF E

d ln[Gg'(E)/Gg'(E)]]/(ni).

(25.8)

280

ROGER HAYDOCK

I

x - pole o - zero --cut

FIG. 10. A schematic picture of the analytic structure of In[@'(E)/G$"(E)]

The Green functions have zeros and poles that introduce ambiguity into the definition of ln{@)(E)/G$"(E)}. However, all the zeros and poles lie on the real axis and so in keeping with the analytic properties of the Green functions we cut the logarithm just below the real axis along the positive real axis at each pole or zero as shown schematically in Fig. 10. We can now generalize to infinite systems where the Green functions have cuts as well as simple poles. We continue to define the logarithm as cut just below the real axis from zero to plus infinity. CF becomes an open contour with endpoints at EF Oi and EF - Oi and, in the case of ground state binding energies, CFextends around all the singularities below EF,as shown schematically in Fig. 11. Equation (25.8) now applies to all cases.

+

Im[EI

I

x -pole o - zero w-cut

Flc. 11. The contour C , shown with the singularities of @"(E)/Gg'(E).

RECURSIVE SOLUTION OF

b b .----. uo

Ua

SCHR~DINGER’S EQUATION

b u1

28 1

b u2

FIG.12. A constant chain with adsorbate.

As an example, let us consider the binding energy of an extra state on a constant-chain model. This is a simple model for an adsorbate on a surface. Figure 12 shows the chain with a dotted connection for the extra state u,. From Eq. (19.4), In [Ghl)(E)/Gbl)@)I =

ln{E[E - ( E z - 4b2)’/’]/(2b2)}, ln{E[E ( E 2 - 4b2)’I2]/(2b2)},

+

Re{E} > 0, Re{E} < 0, (25.9)

for E in the upper half-plane. We suppose that u, was initially doubly occupied if E , > 0 and unoccupied if E , < 0. Substituting into Eq. (25.8) and integrating, we get

AU

=

-2(4b‘ - E $ ) ” 2 / ~ .

(25.10)

This shows that the extra state always binds to the constant chain and that the binding is greatest when the band is half full. If we compared energies for the case of a singly occupied extra orbital, this would introduce an extra term into the binding, namely the energy necessary to bring the extra electron to the Fermi level. FOR 26. PERTURBATION THEORY

THE

CHAIN

Suppose we make a small change in the model Hamiltonian H.It is sensible to ask how this changes the resulting chain model. The answer to this question is a perturbation theory26 of the transformation to the chain model. In this section we develop this theory and in later sections apply it to the calculation of transition matrix elements and other physical quantities. The basic idea of the perturbation theory is to expand the basis states and parameters of the chain model in powers of the perturbing potential. Suppose we have a chain model for H consisting of orbitals Cubo),ui0),. . .} and parameters {aho),u(1O),. . .} and {b\O), bLo),. . .}. Consider now the Hamiltonian H + AV, where A is a dimensionless scalar and Y some perturbing potential. For each value of A, we could transform to a chain model and we find orbitals R. Haydock, J. Phys. A 10,461 (1977).

282

ROGER HAYDOCK

{u,,(1)} and parameters {a,,(A)} and {b,(1)} that all depend on 1. We also know that u p = U,,(O), (26.1) .$O'

=

(26.2)

an(O),

bA0) = b,,(O),

(26.3)

If the orbitals and parameters are analytic functions of 1,we can expand them in a power series about A = 0. The condition for this is essentially that Vdoes not break any symmetry of H. The coefficients of the various powers of 1 in the chain orbitals and parameters can be calculated recursively in accord with the philosophy of this chapter. Let us define (26.4) (26.5) (26.6) By applying the chain transformation to H same order, we obtain

.If)

1

= m

1= O u;m)tSHujt-m)

+ 1Vand collecting terms of the

I- 1

+

1,,;m)tS~"-m-

m=O

1 )9

(26.7)

1

(26.10) where we have introduced the orbitals {x:"')} to simplify the expressions. Equations (26.7)-(26.10) display a number of useful properties of a per, or turbation theory. In order to compute the Ith order change in t ~ , , + ~a,,, l ~ , + one ~ , needs only the lower order corrections for u , + ~ u,,, , I I , , - ~ , b,,+l, b,,, and a,,. Thus the above equations can be used to calculate successive improvements in an approximation. The parameters can only go singular for finite 1 and n if b:? is zero. This is the case when V breaks a symmetry of H, when the chain model for H is finite. Because V breaks a symmetry, there are more degrees of freedom in the perturbed chain, and so blp! will be zero at the end of the unperturbed chain while b:? will not.

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

283

Infinite systems provide far more interesting behavior when perturbed. Let us consider how the orbitals and parameters change as an infinite chain is perturbed by a V that breaks some symmetry. Application of V to u, introduces components not in As n increases, successive multiplications of V tend to make these new components larger and larger. This in Eq. (26.10). The parameter blp!’ process is offset by division by b!,Yp!, represents the strength H. When Y is stronger than H, this leads to an instability characterized by V dominating the asymptotic behavior of the chain rather than H. Physically this represents a phase transition for some tinite L,namely the strength needed for V to dominate the asymptotic chain. An example27of this is a random potential Vperturbing a crystal Hamiltonian H. For sufficientlystrong Ythe extended nature of the eigenstates vanishes and they become localized. A subroutine for computing the quantities in Eqs. (26.7)-(26.10) is available as part of the recursion library.’

uF).

27. PERTURBATION OF THE GREEN FUNCTION AND TRANSITION MATRIX ELEMENTS Small perturbations of a system, such as the oscillating electric field of a photon in a metal, can induce transitions. The transition rate depends on densities of initial and final states, as well as the transition matrix elements, and is the quantity that many experiments measure directly. Thus a simple method for calculating the transition matrix elements is necessary for direct comparison of theory and experiment. The transition matrix elements are the matrix elements of the interaction between the initial and final states of the system. As an example consider the interaction of electrons and photons that leads to photoemission. The initial states are those where the electrons and photons exist separately and the final states are those where the photon has been absorbed and an electron excited. Thus the unperturbed system is described by H, which includes independently the initial and final states, and the perturbation V is the interaction that connects initial and final states. In keeping with the philosophy of the time-independent Schrodinger equation, we view the transitions as between states of the same total energy rather than between states of different energy due to a time-dependent potential. The transition matrix elements are related to the derivative of the Green function with respect to the strength of the perturbing potential. Suppose we want to calculate the transition rate between an initial normalized state

’’ R. Haydock, Philos. Ma.9. [Part] B 3 1 , 97 (1978)

284

ROGER HAYDOCK

xi and a final normalized state x f . Let {vp)} be the initial eigenstates of the system with energies {E:)} and {vi?} the final eigenstates with energies {E:')}, then the transition matrix element at energy E is

1

T(E) = (xftsvp')(V~?+SVvj;'))(v~~+SXi), E = Eif) = EW.

(27.1)

B

Using the spectral representation of the Green function, Eq. (17.2), we can express T ( E)by,

T ( E ) = R,[x!S(E

-

H)-'V(E - H ) - ' x ~ ] ,

(27.2)

where R , [ f ( E ) ] means the residue of the second-order pole of f at E . This is the coefficient of 1/(Z - E)' in a Laurent expansion of f(Z)about the point E. We may now relate Eq. (27.2) to diagonal elements of the Greenian by defining ~0

= (xi

+ x,>/JZ,

(27.3)

and considering the Green function Go(E) = u ~ S ( E- H - AV)-'UO.

(27.4)

Assume that S, H , and V are all real. Denote the first derivations of Go@) with respect to A at A = 0 by Gb(E), then T(E) = R,[Gb(E)l.

(27.5)

Gb(E) can be expressed simply in terms of the perturbed chain described in Section III,26. In Section II,17 we showed that G,(E) is the zero-zero element of ( E - J ) - ' , where J is the tridiagonal matrix of chain parameters for H AV. Using the perturbation theory we determine J , and J , so that

+

J

=

J,

+ A J , + higher orders.

(27.6)

J , is just the chain parameters for H alone. Using this we find T(E) = R , { [ ( E - J J ' J , ( E

- J0)-qoO}.

(27.7)

Define J g ) as the matrix made from the first n rows and columns of .lo. In order to compute T(E)from Eq. (27.7) we need more than the zero-zero element of ( E - Jo)-', which is G,(E). We define polynomials Q!,'"'(E) that are proportional to the (1, rn) cofactor of ( E - 56"))so that [ ( E - J g ' ) - l ] l m = QP*"'(E)/[b,P,(E)].

(27.8)

285

RECURSIVE SOLUTION OF SCHR~DINGER'SEQUATION

These new polynomials are defined by the following equations:

Qlf*")(E)= 0,

for n Imax(1, m),

(27.9)

Qlf*"-')(E)= b$O)PI(E)/bi0), for 1 < n,

Qlf,"'(E)

=

Qim3 "(E).

-

a$")Qlf,"'(E)

(27.10) (27.11)

They obey the recurrence

Qlf;";'(E) = [ ( E

-

(27.12)

b!,')Qi'$)(E)]/b!,?

From this it is clear that they are just a generalization of the Q,(E) defined in Section 1,ll so that Qn(E) = Q!,090'(E). (27.13) In terms of these polynomials,

(27.14) The nth approximation to T ( E ) , T,,(E), is T,(E) =

n- 1

1 [ U ~ " Q ~ ~ * "+( Ebil)Q$o,l-l)(E) )

I=O

+ bi:'lQ!,ov

I+

')(E)]Qil."(E)/[b,Pn(E)]',

(27.15)

where P,(E) = 0, and P&E) is the derivative of P,(E). For an infinite system we define a transition density z,(E), T,(E) =

1 T,(Et') 6(E - @)),

(27.16)

a

where { E t ) } are the zeros of P,(E). The local density of states of H depends on A and from Eq. (17.18) is n6"'(E) =

c &'(A)

6 [ E - Ef'(A)].

+ LV

(27.17)

a

Differentiating and taking the second-order residue gives z,(E) =

1&)(O)E$''(O) a

6 [ E - Et)(O)].

(27.18)

The convergence properties of the transition density are thus related to those of the density of states. This method has not been applied yet in any calculation of transition densities although the method should make the problem almost as simple as densities of states.

286

ROGER HAYDOCK

28. PERTURBATION THEORY APPLIED TO INTEGRATEDQUANTITIES In Sections III,24 and III,25, we discussed the effects of small changes in the Hamiltonian of a system. It is clear that the results of the perturbation theory of the chain and the Green functions apply to those same problems. Previously we stated that the errors in calculating binding energies and changes in integrated quantities were much smaller than the absolute errors of integrals of Green functions. Using the perturbation theory we can now see the explanation for this. Let us consider the weights pJ1) and eigenenergies E,(1) as defined in Eq. (17.18) as a function of the coupling parameter in the Hamiltonian H 1V.If V represents some small change in the system, then the changes in p a ( l ) and &(A) give the change in the projection of eigenstates on the initial orbital uo and the change in eigenvalues as A varies. In terms of these, the spectral decomposition of the Green function is

+

(28.1)

By considering the changes in pa and E, with 1 we can connect changes in various integrated quantities with the perturbation theory. The change in Go@) is, to first order in AA, AGo(E) x

1 {ph(J-)/[E - Ea(1)I + pa(1)%(1)/CE a

-

Ea(1)12I A1* (28.2)

The integrals discussed in Section II,18 can be expressed as contour integrals so that Eq. (18.19) becomes (28.3) where Co is a contour that encloses all singularities of Go(E) with energy less than E , . To first order in Ad, (28.4) The contour only picks up the first-order poles, and so p&(1) is the part of the change in Go@) that produces changes in this kind of integral. At 1= 0,

+ bj Qk0* :)1

I

+

')(&)IQlf. O)(E,)/[b E(E,)] '.

(28.5)

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

287

In Section III,25 we were concerned with the changes in position of the poles, and so Ea(A + AA) - Ea(A) x Ea(A)AA,

and at A

=

(28.6)

0, n- 1

1 [ U { ” Q ~ ~ ~ ~+) (bj”QAo*’E,) ‘)(Ea) + b):) Qko7 ‘)(EX)]Qk1,‘)(E,)/[b C(E,)]’.

pXEL =

I=O

I+

(28.7)

It should be noted that if Vbreaks a symmetry of H then PL and EL will not be defined for the eigenvalues that are split by V. The methods in Section III,24 compare identical clusters with and without V and amount to a numerical differentiation of the Green function with respect to A. The perturbation theory allows analytic differentiation of the Green function with respect to A, and shows that this derivative is a quantity that depends on cluster size independently of the Green function itself. Thus an error theory of the derivative might be constructed analogously to that for the Green function, but the errors are independent. It may be quite easy but it has not yet been done. Some ideas of the errors come from relating ph and Eb to matrix elements of V between eigenstates of H.If {w,} are the eigenstates of H,then (28.8)

EL = w ~ S V W , .

Let uo represent the initial orbital of the chain model, then pi = 2 Re

r

1

1

1 (~~S~,)(W~SVW~)(~~SU~)/(E~ - ED) . (28.9)

L * a

Thus the errors due to calculating perturbations in a cluster result from the difference between the matrix elements of Vbetween eigenstates, of the cluster and of the infinite system. If V is localized, this can be quite small and even if V is not localized, these errors decrease as the cluster gets larger. The convergence of the transition density is in terms of indefinite integrals because of the discrete spectrum of clusters. IV. Relation to Other Methods and Conclusion

29. RELATION TO OTHER METHODS There are a number of other methods for solving the Schrodinger equation in circumstances similar to those discussed in this chapter. In the next few sections we present briefly the advantages and disadvantages of several of the

288

ROGER HAYDOCK

more popular methods. But first, in this section, we examine the general relationships between methods. All methods reduce to the calculation of the product of the Hamiltonian on a state. A Hamiltonian may already be small enough to store in a computer as a matrix or, by astute use of symmetry, it may be reduced to a blockdiagonal form so that each block can be stored in a computer as a matrix. These Schrodinger equations can all be solved efficientlyby standard computer matrix methods. The problems to which we address ourselvesare typical of those resulting from the local atomic environment approach to solidstate physics and have Hamiltonians that are too large to be stored as a matrix but whose matrix elements are simple to compute. Thus we can compute with vectors from the space operated on by the Hamiltonian but not matrices. Hence all methods must compute the vectors or parts of vectors resulting from multiplication of a vector by the Hamiltonian. The information produced by different methods is equivalent. We can reduce the problem still further by considering the action of the Hamiltonian on one vector at a time. If H is the Hamiltonian, and u our chosen vector, then the various vectors H"u can be calculated. We might wish to find vtSH"u, but this can be written as sums of diagonal elements exactly as in Eqs. (17.3)-( 17.6). This reduces the problem to calculating udSH'u,, which is also called the nth moment p,. This is because p, is the nth statistical moment of the local density of states, p, = u&SH"u,=

J-

03

E"no(E)dE.

(29.1)

m

This can be proved by writing the integral over the density of states as a contour around the singularities of the Green function. Thus every method can be described by how it uses moments to calculate physical properties. The moments serve to relate one method to another. Because there is no other information than the moments, all methods are analytically but not necessarily computationally equivalent to moments and hence to each other unless one method loses or ignores some of the possible information. The recursion methods presented in this chapter are analytically equivalent to computing moments. The construction of a chain model is based on a starting vector u, on which we repeatedly operate with H . Since the chain model is just the representation of H in terms of a different basis, we can easily relate the moments to the chain parameters, (29.2) where J is the tridiagonal matrix of chain parameters in Eq. (4.3). It is possible to invert this relationship to find the chain parameters in terms of the moments. These formulas involve determinants of matrices of moments

RECURSIVE SOLUTION OF SCHR~DINGER’SEQUATION

289

and are contained in Shohat and Tamarken3 and Akhiezer.2 Thus n-chain parameters are equivalent to n moments although one should note that chains may be finite, but the set of moments is always infinite. 30. THEMETHOD OF MOMENTS

A widely used method of approximating the density of states is based on calculating its moments.’* A finite number of moments are straightforward to calculate, but then there is the problem of finding a density of states from them. One way is to choose a form with adjustable parameters and fit the moments. This is fine for a few moments, but for large numbers of moments it only works if the moments are linearly related to the parameters. This brings up a fundamental difficulty. The density of states is a nonnegative distribution and there are no functions, all linear combinations of which are nonnegative. Thus, one either gives up the physical property related to hermiticity or gets into a large nonlinear problem. In fact, the continued fraction-related to the chain model-is the classic solution to this problem but, of course, the continued fraction or chain parameters are not related to the moments linearly; see Eq. (29.2). However, most users of moments convert them to a continued fraction by the inverse of Eq. (29.2). The practical difficulty with the moment method is that the numerical calculation of a density of states from moments is unstable. If we have a band with limits fE,, its even moments are of the order of E;. Any feature within the band, say at El, makes a contribution of the order of Ey to the nth moment. As n increases, the fraction of the nth moment due to this feature decreases exponentially. If we are trying to construct a density of states from its moments, the accuracy of the computation must increase exponentially with the number of moments in order to use the information in the higher moments. Methods that partially avoid this problem use polynomial moments (30.1)

where {P,(E)} is some set of polynomials, usually orthogonal with respect to a given weight function. The growth of the A, is different from the ,u, and if the polynomials are “well suited” to no(E), the numerical instability is greatly reduced. Here, “well suited” means that the weight function for the polynomials and no(,?) have the same essential singularities in the sense of Section III,20.

’’ F. Cyrot Lackmann, J . Phvs. (Paris),Syppl. CI,67 (1970).

290

ROGER HAYDOCK

The polynomials obtained from the chain model are optimal for polynomial moments. The orthogonal polynomials that are the wave functions of the chain model are all orthogonal to no(E) except for P,(E), and so for this set the polynomial moments are all zero except for A,. Thus starting from moments and eliminating their numerical instabilities one comes back to the chain model. This is an alternative development of the recursive methods of this chapter. 3 1. TRANSFER MATRICES

Graphs more general than the chain lead to recurrence relations for the Green functions, which when expressed in matrix form, involve transfer matrices.29 In Sections I,4 and I,15 we mentioned that graphs other than chains led to recurrences of more than three terms or vector recurrences for the Hamiltonian. Each recurrence simply specifiesthe effectof the Hamiltonian on a given basis element. These recurrences lead to nonlinear recurrences for the Green function, analogous to Eq. (17.11). The general Green function recurrences involve Green functions that are restricted to certain parts of the Hamiltonian just as the G,(E) in Eq. (17.11) are restricted to the chain with the first n orbitals deleted. In its most general form, the transfer matrix method solves the Green function recurrences to get an approximation to the desired Green function. The philosophy of the transfer matrix method is not to change basis as with the chain model but rather to display the properties of the basis explicitly in the recurrences for the Green functions. The connections in the graph of the Hamiltonian give similar connections between Green functions in the recurrences, while no matter what the Hamiltonian is, its chain model gives a three-term recurrence. The relationships in the Hamiltonian are reflected in the chain parameters rather than the number of terms in the recurrence. In problems that are solvable analytically, the transfer matrix method leads straightforwardly to the solution. Simply terminating Green function recurrences can lead to wrong analytic structure of Green functions in the transfer matrix method. The general recurrences connect Green functions for different parts of the Hamiltonian. To get a sensible answer the recurrences must relate them consistently. If a subset of the recurrences is used to approximate a Green function, the subset may not be consistent. For example, termination of the recurrence in Eq. (17.11) is the same as cutting the chain and always gives an approximation with correct analytic properties. More complicated recurrences cannot be terminated so easily. 29

L. M. Falicov and F. Yndurain, J . Phvs. C 8, 147 (1975)

RECURSIVE SOLUTION OF

SCHR~DINGER’S EQUATION

29 1

The transfer matrix method is only useful for computation when the recurrences have a constant number of terms. This corresponds to the vector chains in Section IJ5, where the transfer matrices are the matrix generalization of the chain parameters. Otherwise the computations become very complicated. 32. FEENBURG PERTURBATION THEORY Feenburg perturbation theory3’ is an example of the transfer matrix approach. We write all the recurrence relations defining the Hamiltonian in terms of basis orbitals {Z,} and then solve for the local Green function of Z, by successive elimination of orbitals from all the equations to get

c

c HOIHlOG:O1(E)]

-1

Go@) = E - Hoo -

2 1

+zo

2 2

+Z0,Zl

1

E - Hll -

9

I-’

H12HzlG(zO~’1)(E),

(32.1) (32.2)

Each recurrence involves a Green function with one more orbital eliminated from the Hamiltonian. The equations give a sort of multiply branching continued fraction for Go@). The branching of the recurrences, Eqs. (32.1) and (32.2), explicitly display the connective properties of the Hamiltonian, but they are awkward for computation. They also suffer from the general difficulty of terminating the transfer matrix equations. Using any subset of the recurrences usually gives an approximation for G,(E) with wrong analytic properties. 33. THECLUSTER BETHEMETHOD The cluster Bethe method3’ is a form of the transfer matrix method that avoids some of the difficulties in approximation. A small part of the Hamiltonian, usually of the order of ten degrees of freedom, is treated exactly while the rest is replaced by a Bethe lattice or Cayley tree. The approximate Hamiltonian represents a small cluster imbedded in a medium described by a Bethe lattice. When the transfer matrix approach is applied to the approximate Hamiltonian, the Green function recurrences can be solved exactly. P. M. Morse and H. Feshbach, “Methods of Theoretical Physics,” Vol. 11. McGraw-Hill, New York, 1953. 3 1 J . D. Joannopaulos and F. Yndurain, Phys. Reii. B 10, 5164 (1974). 30

292

ROGER HAYDOCK

This is because all elements of the Greenian for a Bethe lattice can be found analytically. The recurrences express Green functions for the cluster in terms of those for the Bethe lattice. While this method always gets the analytic properties right, it only gets a few (usually less than four) moments exactly right. The number of moments that are exact is the same as the highest power of the Hamiltonian whose diagonal elements do not have a contribution from the Bethe lattice. Higher moments are correct to the extent which the Bethe lattice mimics the proper environment of the cluster. The Bethe lattice acts as an asymptotic approximation. It has an adjustable bandwidth, but no internal singularities, and so as an asymptotic approximation it is very limited. 34. EQUATION OF MOTION METHOD

The equation of motion method32 approximates solutions to the timedependent Schrodinger equation rather than the time-independent approach we have used so far. The idea is to solve numerically the equation -ih

d dt

-u(t)

=

Hu(t),

(34.1)

where u(t) is a time-dependent state of the system, and u(0) is known. There are a number of ways of doing this. One is to write the formal solution u(t) = exp(itH/h)u(O),

(34.2)

expand the exponential in powers of H,and thereby return to moments. Another is to break time into discrete intervals and use forward and reverse iterations to find a solution. Both of these solutions directly relate to moments and hence are equivalent to other methods. Error theories for the timeindependent methods can be Fourier transformed to get error estimates for this method. As yet this has not been done, nor is there the elegant mathematical apparatus of the time-independent approach.

35. CONCLUSION The main result of this chapter is the transformation of an arbitrary quantum model into a one-dimensional chain model and its subsequent analysis. The chain expresses mathematically the physical concept of local environment. The recursive transformation yields analytic chains for some systems, but it is also convenient and efficient for constructing numerical 32

D. Weaire and A. R. Williams, J . Phys. C 10, 1239 (1977).

RECURSIVE SOLUTION OF SCHRODINGER’S EQUATION

293

chain models enabling the solution of problems that are too big for numerical matrix methods. Eigenvalues, eigenstates, Green functions, and a host of other physical properties of the chain model can be accurately and efficiently computed by recursive methods. Apart from its computational value, the chain suggests a new approach to quantum-mechanical models. Because any quantum model can be transformed into a chain model, chain models exhibit all possible quantum phenomena. Hence physical problems can be qualitatively related directly to chain models. Because of the simple solution of chain models, the qualitative behavior of different physical properties can be determined. The chain model is soluble in a strong sense. Because it is related to the classical moments problem and orthogonal polynomials, a large number of the results from these subjects apply to the chain model. Unlike many methods for solving quantum models, one has rigorous results about the convergence of approximations. Because they are defined recursively, the approximations are ideally suited to computation. There are a number of open problems concerning the recursive solution of the Schrodinger equation. The most important one is that of determining the asymptotic form of the chain model of an infinite Hamiltonian. One can easily compute exact chain parameters for the first few levels of the chain, and rather than proceeding beyond that point with a chain model of the finite cluster, which deviates from that for the infinite system, one would like to determine an asymptotic form for the chain. The difficulty of this appears to relate to determining an asymptotic form for the chain states u,, a subject we did not discuss in the sections on asymptotic methods. In an infinite tight-binding model, u, is contained within the “sphere” of orbitals that are within n hops of u,, . Depending on the symmetry of the model, orthogonality localizes u, on the surface of this sphere for the highly symmetric Cayley tree or allows u, to fill the sphere as in disordered materials. The u, transform to the polynomials P,(E) when H goes to E, and so a simpler form of the problem can be stated as that of finding an asymptotic form for P,(E) given no(E). However, at some point one must grapple with the problem of an asymptotic form for u, in order to relate the asymptotic chain to an arbitrary model. A second open problem is to relate the asymptotic form of a chain to essential singularities in the Green function. The theory of Sections III,20 and III,21 makes a start in this direction, but there is much more to do. In general, the problem of asymptotic approximations is a rich one and particularly so in the sense in which a continued fraction is nonlinear. We have not discussed statistical properties of chain models at all in this chapter. One of the outstanding problems is that of the Anderson model of a disordered solid. Here the Hamiltonian is defined statistically and one would

294

ROGER HAYDOCK

like to transform to a corresponding statistically defined chain where one could discuss the distributions of various physical quantities. The first steps have been made in this direction.27 The application of these methods to some of the most interesting problems in disordered solids will depend on some ability to determine distributions of chain models. Finally, the reader may have noticed that the approach of this chapter has been entirely directed toward independent particle models. However, many-body Hamiltonians must be similarly transformable to chain models. Aside from a small amount of thought and continuing interest, nothing has been done about this and it remains a most interesting problem concerning recursive solution of the Schrodinger equation.