A New Expansion For Energy Spectra Of Solids

A new cluster expansion (CE) is proposed to calculate band structure and adiabatic potential of an insulator at high pressures. This kind of the CE is developed to permit using Lovdin orthogonalization of atomic orbitals for innnite crystals. The strong binding approach (LCAO method) and OPW one are modiied. The technique suggested also ts for the crystals having impuirities, d-electrons, etc.


Introduction
The aim of our study is to create the theory allowing one to describe the properties of an insulator at high pressures.The method of calculating energy structure and adiabatic potential of a crystal is proposed in this work.
Available methods of band structure calculation are cumbersome in the case of high pressure and require much machine time.Moreover, we cannot be sure for these methods to give correct results at high compressions.So, at present time we don't know regular ab-initio methods for band structure calculation of insulators under extremal conditions.The aim of the work suggested is to ll in this gap.
For high compressed crystals (or having d-electrons, impuirities, etc.) one has to take care for changing electron wave functions of atom in crystal as it is being compressed.In view of the fact, two methods of band structure calculation, namely, the LCAO and OPW are modi ed by Lovdin orthogonalization and cluster expansion (CE) developed in 7] for the sake of building Bloch functions.The modi cations are called CELO and CEOPW respectively (chapters 2 and 3).In the chapter 4 of this work new equations for the adiabatic potential (AP) of a crystal are obtained in reciprocal space using CE.
In conventional LCAO method trial Bloch functions are built as for isolated atom orbitals neglecting their overlapping in crystal 1].In the OPW method, core states are considered by this approach.Beside, there is a number of works 2-5] taking account for the overlapping by the rst order in overlapping integrals (OI).However, for the cases mentioned above overlapping among orbitals of neighbour atoms seems to be too large to be treated in those 1-5] ways.Therefore, trial Bloch functions are required to be constructed by exactly orthogonal Wannier functions instead of isolated atom ones.
The Wannier functions are believed to be approximated by linear combinations of atomic functions orthogonalizated each to others after Lovdin 6].However, unlike of molecules, in extended solids exact Lovdin orthogonalization is an extremely di cult calculational problem.It is the CE that permits practical work.Using the CE, the problem of Lovdin orthogonalization for localized orbitals of whole crystal is reduced to similar problems for respective clusters.These problems may be solved easily in low orders of the CE (the CE for band structure calculation is shown 7] to converge rapidly with n (as a power law in the OI)).
The CE developed by the authors for band structure calculation is analogous to Abarenkov-Antonova's one suggested for electron density of a crystal 8-10] (the CE of 8-10] is used for calculating the AP).
The approach suggested are applicable to calculate the same bands as their "nonorthogonalized" origins but in more wide scale of pressures.They also t without limitation for the substances having d-electrons, impurities, etc.

CELO method
This method is a modi cation of the LCAO one.It is disscussed by example of lled bands of insulator.The OI among basis orbitals of the LCAO method may be large as such in the case of high pressures (impuirities, delectrons).Moreover, they are included in LCAO equations being multiplied by the nearest neighbour number z ~10.So, high-order OI corrections should be taken into account 4].To permit any further approximation we build, before all, functions ls exactly (after Lovdin 6]) orthogonalized each to others (we use those as a basis to construct exactly orthogonal set of trial Bloch functions), ls = X l 0 X s 0 ' l 0 s 0 (r) T ?1=2 l 0 s 0 ; ls : ( Here ' ls is an isolated atom function, l is a lattice vector, N is the number of unit cells (the crystals with one atom per cell are considered), s = (n; l; m) denotes occupied states of atom.Each one-particle state (all the quantities are equal Z=2; Z being the nuclear charge) is supposed to be occupied by two electrons with spins directed "up" and "down".There is no division between core and valence electrons states in this work excluding its last chapter.The orbitals ' ls belonging to the same atom are regarded to satisfy the Shr odinger equation.In equation (1) T is the so-called 8-10] metric matrix with dimension L L (L = N Z=2).Its elements are the OI between respective atomic orbitals ' ls T l 0 s 0 ;ls = Z ' l 0 s 0 (r)' ls (r)d 3 r; l; l 0 = 1 : : : N; s; s 0 = 1 : : : Z=2: ( T is a hermitian matrix.In equation ( 1) the matrix T 1=2 such that T 1=2 T ?1=2 = I; T 1=2 T 1=2 = T and T ?1=2 !I as T !I is used.
The next step is building a trial Bloch function of CELO method using equation (1) k (r) = X s ks (r)C s ks (r) = N ?1=2X l exp(ikl) ls (r); (3) with C s being variational parameters.Index stands for lled bands.In order that the trial functions equation (3) be orthonormal it is necessary and su ciently for matrix of coe cients C s to be unitary.
Altering the equation E k = h k jHj k i=h k j k i for band energy on the C s we obtain the following set of equations for C s X s 0 fh ks j H j ks 0 i ?E k ss 0 g C s 0 = 0: (4) Here H = T + V (r) is the Hamiltonian of an electron, with V (r) including interaction with all nuclei (not with cores) and self-consistent interaction with all the rest electrons in the crystal.Conventionally the V (r) is presented by the superposition V (r) = P l v l , with v l v(j r ?l j) being self-consistent potential centered on the l-th lattice site.The energy E k of a band electron is found by solving secular equation for the equation ( 4) with coe cients h ks j H j ks 0 i.
The clusters which di er each of other only by transposition of indices in writing l 1 : : : l ] are considered to be equal.Here we need in matrices of A t = I ?T ?t with t = 1; 1=2; A 1=2 = P equation ( 6), A 1 = P equation (15) (see below on): A t m] = P m]; t = 1=2 P m]; t = 1 = I ?(T m]) ?t ; A t N] A t : (8) All matrices in equation (7) have dimension L L. The T m] denotes cluster metric matrix.According to the de nition 8-10] it is T m] ls;l 0 s 0 = T ls;l 0 s 0 ; l and l 0 2 m]; ll 0 ss 0 in other cases: One can see that each (e.g.n-th) order term of the CE is constructed to take account of the overlapping integrals among LO of atoms belonging only to the same n-particle cluster.For every n such clusters comprise whole crystal.Beside, each term of the cluster expansion is formed to avoid multifold taking into account the contribution coming from the same cluster.Indeed, the n-th order term includes (due to the multiplier (?1) m+n ) contribution of the nparticle cluster minus analogous contributions of the n ? 1, n ?2-particle, etc. clusters contained in this n-particle one (these contributions were already accounted in the previous orders of the CE).
Although the basic formula of CE (equation ( 7)) had been obtained in 10] for an arbitrary matrix A t , the theorems which permit practical calculations had been proved there only for speci c case of matrix P (equation ( 8) with t = 1).The equations for Bloch functions necessary for band structure calculation contain the matrix P (equation ( 8) with t = 1=2).Therefore, these theorems should be rederived for the P.Here we present only the formulations, with the arguments given in 7].
Theorem I.All elements Q 1=2 n] ls;l 0 s 0 with l and(or) l 0 = 2 n] are equal to zero.
Theorem II.The Q 1=2 n] is equal to zero even if orbitals of one atom of the cluster n] are not overlapped with all the rest orbitals of this cluster.
This property allows to develop the nearest (or the second, the third, etc.) neighbours approximation for the Q 1=2 n].
Theorem III.Extracting square root of cluster matrix T m] it is enough to solve characteristic equation det k I ?T m] k having an order m Z=2 instead N Z=2 as for whole metric matrix T, with m usually being much smaller than N.
Using translational symmetry, it may be shown 7] that n] ks (r + a) = n;?a] ks (r) exp(ika); (10) where n; ?a] l 1 ?a : : : l n ?a] is a cluster shifted by the lattice vector a with respect to the original cluster n].Nevertheless, the function n ks (r) from equation ( 9) is a Bloch function i.e. n ks (r + a) = n ks (r) exp(ika).
Substituting the equation ( 9) into coe cients h ks j H j ks 0 i of equation ( 4) we obtain the CE for them.Solving secular equation for the equation (4) using two-, three-particle, etc. clusters we obtain band structure of a crystal in respective approximation.

CEOPW method
Let's go on discussing the CEOPW method as to conduction bands of an insulator.In this method basis plane waves are orthogonalized to wave functions k (r) of lled bands like they are done in ordinary OPW one.A trial Bloch function kc (r) of a conduction electron is kc (r) = kc (r) ?P ;k 0 k 0 (r) R k 0 (r 0 ) kc (r 0 )dr 0 ; kc = V ?1=2P g a k (g) exp f?i(k + g)rg : (11) Here a k (g) are variational parameters, g is reciprocal lattice vector.
The rigid core approximation is used for describing occupied states in conventional OPW.In the CEOPW method the lled band functions k (r) may be orthogonalized each to other as preciously as we desire being calculated by the CELO method.Owing to that the orthogonalization of equation ( 11) is of the same accuracy.In the conventional OPW method this orthogonality is only in the rst order in the OI 1].
Substituting equation (1) into equation ( 4) and writing ' (r ?l 0 ) = X s ' s (r ?l 0 )C s ; we bring the k into the form suitable for further work (see equations ( 3), ( 6)) k (r) = N ?1=2P l exp(ikl) ' (r ?l) ?P l 0 0 ' 0 (r ?l 0 ) l 0 0 ;l ; l 0 0 ;l = ll 0 0 ?P ss C s 0 0 (T ?1=2 ) l 0 s 0 ;ls C s ; (12) with C s being variational parameters of the CELO method.They are regarded to be already calculated.The functions k in the equation (12) depend on pressure both explicitly (via the OI contained in the ) and implicitly (via the C s contained in the ' and ).
Using equation ( 12) we obtain for the a k (g) the following set of equations P g 0 fH gg 0 (k) ?E kc gg 0 g a k (g 0 ) = 0; H gg 0 (k) = h 2 2m (k + g) 2 gg 0 +V C (g ?g 0 ) + V ex (k + g; k + g 0 ) + V PK (k + g; k + g 0 ): (13) Here V C (q) = 4 e 2 = ?q 2 0 ( (q) ?Z) ( 14) is the form-factor of the Coulomb interaction of exceeding electron with all the rest ones, (q) is the Fourier transform of their density (r) = 2 X ls ( ' s (r ?l)' s (r ?l) ?X l 0 s 0 ' s (r ?l)P ls;l 0 s 0 ' s 0 (r ?l 0 ) ) : (15) In equation (15) 0 is the unit cell volume, sums go over atomic states instead of because = P k k k = P ks ks ks since the matrix of CELO variational parameters is unitary.

Adiabatic potential
Adiabatic potential (AP) is the energy of electrons in crystal, with the ions being xed in arbitrary positions R = l + u l close to equilibrium ones l.
In the Hartree-Fock (HF) approximation the AP is expressed only via oneelectron density matrix (r j r 0 ; fRg) depending on ion coordinates fRg as on parameters.In this approximation the AP consists of electron kinetic energy T e , electron-electron Coulomb and exchange interactions, U C and U ex respectively, electron-ion interaction U ei , and ion-ion one, U ii .In this work the (r j r 0 ; fRg) is expressed in terms of s (r ?R).
Using CE, the equation for AP had been obtained in the 8-10] in rspace.This result is very useful for the imperfect crystals with impurities chaotically spread in the crystal.For perfect crystals under pressure (or for the cases when impurities form ordered superlattices) it is convenient to work with the CE in reciprocal space.Expanding the (r j r 0 ; fRg) as a Fourier series we express the AP in reciprocal space in the following way T e = h 2 =(2m) X q q 2 f(q j q); ( U ex = ?e 2 =(4N 0 ) X pq v C (q) j S(p) j 2 f( j +p)f ( +q j +p+q); (20) Here the f( j +p) = 2 X ss 0 I s ( )V ss 0 ( )I s 0 ( +p); F(p) = X f( j +p); (23) respectively, are Fourier transforms of the nondiagonal and diagonal elements of electron density matrix of an atom a ected by the electrons of all the rest atoms of crystal through the orthogonalization.It may be shown that f(q j q 0 ) = f (q 0 j q), F(k) = F (?k).The c (k) is the Fourier transform of core electron density.In the equations ( 20)-(23) V ss 0 (q) = ss 0 ?P ss 0 (q); P ss 0 (q) = X h P hs;0s 0 exp(iqh); where h = l+u l ?u 0 is the distance between two atoms, one of them being placed near zero lattice vector (u 0 is its displacement).Using the CE for matrix P (equation (7) with t = 1) developed in 8-10] we obtain the CE for HF part of the AP in reciprocal space.
Beside HF part (equation ( 20)-( 23)) the AP contains many-electron contributions, e.g.Van-der-Vaals interaction.There is a number of problems, with many-particle e ects being dramatic.Speci cally, for rare gas crystals the Van-der-Vaals term is a single negative part of the AP providing lattice stability.This term may be obtained, for example, by the second order of perturbation theory in the di erence between exact and HF potentials 11] E (2) = ?X kk 0 0 0 jhk k 0 0 jv C (1 ?p)j 0 c 0 cij 2 E c + E 0 c 0 ?E k ?E k 0 0 : (24) Here v C = e 2 =r is pair the Coulomb interaction potential, p is the operator transposing primed and unprimed indices.The indices v and C are used here for valence and conduction bands, E k and E c being respective band energies.
Unlike of 11], the functions j k i and j ci are proposed to be calculated by the CEOPW and CELO methods respectively (instead of strong binding approach).

Discussion
So, the technique is developed to permit the calculation of band structure and adiabatic potential in the same approach under high pressures.Beside, the methods suggested are necessary in the case of impurities, d-electrons even at normal pressure.
Another application of these methods is phase transitions under pressure, speci cally, the insulator-metal transitions in rare gas crystals.However, one delicate moment appears in this theory: The results obtained here are valid only in adiabatic approximation, with the adiabatic parameter = h!D =E g (! D is Debye frequency) containing forbidden gap width E g in the denominator 12].As the pressure increases this gap decreases and the adiabatic parameter may become larger at the transition point.Then, phonons play an active part in forming band structure.Fortunately, insulator-metal transition pressure is estimated 13] to be nonsensitive to the electron-phonon contributions.Thus, the approach suggested does not lose its validity for nding the transition pressure.
Finally, the question about quantitative expedience of using the theory developed occurs.We discuss the behaviour of matrices S = T ?I, P and P in two-particle cluster approximation 7] and the nearest neighbour one for Ne and Ar.In the Tabl.I and II the largest elements of these matrices are given as a functions of relative compression V=V 0 (the matrix S is the rst order in the OI of the matrices 2P and P).It is seen from the tables that the "cluster" responsible for the overlapping di ers of the rst order in overlapping integrals insigni cantly ( and 0:01) until V=V 0 0:7 (about at this value of compression the transition "insulator-metal" is predicted to occur 4,5]).So, the theory for band structure calculation 2,3] and for describing insulator-metal transition 4,5] under pressure seems to be quantitatively precious for Ne and Ar.In the case of hard elements (for example Kr, Xe) this deviation is larger due to the presence of delectrons.
So that the next orders in overlapping integrals are required.In other words, cluster expansion should be used.
Table 1.Ne.The overlapping of the nearest neighbour orbitals.