- Email: [email protected]
Analytic integration of the continued fraction expansion of a density of states
Solid State Communications, Printed in Great Britain.
ANALYTIC
INTEGRATION
Vo1.50,No.5,
pp.401-404,
OF THE CONTINUED
FRACTION
1984.
EXPANSION
003%1098/84 $3.00 + .OO Pergamon Press Ltd.
OF A DENSITY
OF STATES
G. ALLAN Laboratoire
d'6tudes
des Surfaces et Interfaces*,Institut Superieur d'Electronique 3,Rue Franqois Baes, F59046 LILLE Cedex (France)
du Nord
M.C. DESJONQUERES Centre
d'Etudes
Nucleaires de Saclay,Service de Physique des Atomes F91191 GIF SUR YVETTE Cedex (France)
et des Surfaces
D. SPANJAARD Laboratoire
de Physique Received
des Solides *,Universitk
November
de Paris-Sud,F91405
28, 1983, in revised
form February2
ORSAY Cedex (France)
1984byJ.Friedel**
We show that the continued fraction expansion of a density of states can be analytically integrated and expressed in terms of simple elementary integrals.
I INTRODUCTION
II SOME CONTINUED
Since a few years, the continued fraction method [l] has proved to be very efficient to calculate the electron or the phonon density for systems without 3-dimensional translation symmetry. It has been used near bulk defects (vacancies,dislocations,twin boundaries), near surfaces and for bulk disordered alloys. One of the main problem in this method was the termination of the continued fraction as one cannot calculate an infinite number of levels. Solutions have been given for a density of states with gaps [2] or for strongly oscillating recursion coefficients [3]. Here we consider the semi-elliptic termination assuming all the coefficients are equal after the n-th level. But we concentrate on the continued fraction integration. A numerical integration can always be done but sometimes one needs more accurate results. This is the case when one calculates the variation of the cohesive energy with interatomic distances and its first and second order derivatives in order to get the equilibrium distances and the interatomic force constant_s6 In such cases, an accuracy better than 10 must be reached. This is quite difficult using a numerical integration when the density of states shows very sharp peaks. Such peaks arise for example when a bound state is close to the band edge.
continued
FRACTION
PROPERTIES
Let US first consider fraction F,,(z):
a
truncated
1 F,(z)=
(1) bl x-al-
b
n-2
b x-an_l- -
n-l
s-a n It can be written
as[l]:
(2)
where Rn_l(~) and S,(z) are respectively (n-l)th and n-th order polynomials which satisfy the following recursion formulae: Rn= (z-an+l)*Rn-l-bn*Rn-2
(3.a)
Sn= (x-an)~Sn_l-bn_l*Sn_2
(3.b)
with
In the next section, we show that a simple analytical method can be used to obtain the exact result. Let us quote that a similar method has been recently developed by other authors[4] but limited to a continued fraction with few levels. Here we generalise to any order provided one is able to calculate all the real and complex roots of a real polynomial.
RO(z) and So(z) equal
to 1.
One can also easily Rn_l&,-l
- Rn-2*S, =b n-l+.(Rn-2&,_2
*
Laboratoire associe au Centre Nationalde la Recherche Scientifique. ** On behalf of the late Prof. A. Blandin's office 401
show that:
=b n-l bn_2
- Rn_3&.l)
.. b2 (R1*S1 - RO*S2)
=b n-l bn_2 . . .. . b2 bl
(4.a) (4.b) (4.c)
402
THE CONTINUED
FRACTION
EXPANSION
Such a truncated fraction is obtained for a finite cluster. The S,(E) polynomial roots are equal to some of the cluster eigenstates. One often embeds this cluster into a" infinite effective medium i.e one approximates the following levels to get for example the right band limits:
1 G,(z)=
(5) bl
OF A DENSITY
OF STATES
Vol. 50, No. 5
- b" b"_l....b2
bl
Im G"(E)=
JT
(13)
D" where Dn=2b Sz- 2b" S"_l S (E-a)+ 2b2S2 n n-l n
(14)
When E does not belong to the band, Im G,(E) is equal to zero save possibly for some discrete values of E corresponding to bound states.
z-alAs the highest coefficient of S" is equal to 1, the 2n-th coefficient of D" is equal to 2b(l-b"/b). This shows that D" is a 2n-th order polynomial if b is not equal to b and a (2n-l)-th order one if "bn is equal to be
b n-l x--a,_1s-an-b,,f(s) Let us recall here that the density states n(E) we are dealing with is simply
of
If we roots, we get
call
Ei
the
2n
(or
2n-1)
D
n
'i Im G"(E)= n(E)=- 1 Im G"(E) ll
1
J I(E-a)2-4bl
(6)
(15)
E-Ei
i where
From F"(z), to get G"(z), simply to formally replace bn_l by b,,f(z)/(z-a")).
one has b"_l/(l-
- b" b"_l
.. . b2 bl
ci=
(16) Dr!,(Ei)
Then one gets
and DA is the first-order
derivative
of D".
R"_l(x)-b"*R"_2(z)*f(z) (7)
G"(z) = S"(z)-b,,S"_l(z)"f(z)
Here we assume that all the recursion coefficients after the n-th level are equal to their infinite values a and b. In this case, f(z) is given by
Let us notice that only some of the real roots of D" correspond to bound states as some spurious states are roots of the conjugate denominator in (12) and are not poles of G". To integrate over energy the imaginary each part of G,.,(E), one now has to evaluate integral
1 f(z) =
(8) z-a-b,f(z)
EF
&zz
Ii=
dE /
or along the real axis (E-a)-S JR
where Ei is complex.
f(E)=
(9)
2b
One can show [5] that
where
RF Ii
R = )(E-a)2-4bl
s =
+1
S = +i s = -1
(17)
E-Ei
a-Ldb
E < a-2v'b a-2Jb < E < a+2Jc a+26 < E
-
E dE
+
(10)
&-
(11)
+ [4b-(Ei-a)2]
EF
(2a-Ei)
dE K
I EF
dE (18) (E-pi) fi
Within
Then one gets [2b Rn_l - b" R"_2(E-a)l+bn
=
E=a+2 &
R,_,S&
the usual transformation
sin t
(19)
(12)
G"(E)= [2b S" - b" Sn-1 (E-a)]+b"
S,_l SK
one easily gets Ii = 2 fi
This expression has a continuous art when th_e energy E lies within the imaginary band (a-2 Jp b < E < a+Zib) which is given by:
cos tF + (a-Ei)*(tF+v/2)+ tF
[4b-(Ei-a)2]
J
_n,2
dt (20) 26
sin t + (a-Ei)
Vol. 50, No. 5
THE CONTINUED
If we call Ji the integral uF
2
FRACTION
EXPANSION
in (20)
du
Ji = -
(21) a-Ei
I
4JTu u2 +-+1 a-Ei
-1
with u=tg(t/2). If ul and u2 are the roots one gets 2
1
Ji=-a-E.1
of the denominator,
UF-U2 (Log~-Log-
u2+1
UF-U1
u1+1
u2-u1
)
(22)
where Log is the complex logarithm. Let us notice that the two logarithms contraction in (22) can lead to 2v indetermination of the Ji imaginary part. Then adding the weights of the occupied bound states lying below the Fermi level, one can easily get the number of occupied states. The cohesive energy calculation can also be done using the same procedure
=--
I
E n(E) dE
al=0.0024611, bl=0.0798919,
a2=-0.0779382, b=0.0343536
a=-0.0751389,
n (El 16
A
14
10 E Im G,(E) dE
t 1 ci TI.
=
(24)
dE
1
Ii
+ Ei
1
blb2
‘0.
(25)
;
3
6 4
bn_l (26)
=
”
d
8
I
One has also to add the occupied bound if any multiplied by their state energies weights.Let us notice that these weights are the residue of G,(E) on the corresponding poles Eb. They can be easily computed and are given by
wJ.
have shown that a split-off state above the band is found for all investigated band-fillings Nb of the metal. On the other hand, a split-off state below the band exists only when Nb<6. When the hydrogen local density of states Nb'6, exhibits a sharp resonance just at the bottom of the band as seen in Fig. 1 corresponding to an hydrogen atom sitting on a centred site at The parameters of the equilibrium distance. calculation have been given in [6], and lead to the following values of the coefficients of the continued fraction (the Rydberg is taken as energy unit):
(23)
=F
1
II I
a-
403
OF STATES
12
RF
EC =
OF A DENSITY
b,S,$,_
2 0
-0.4
I
I
1
1
-0.3 -0.2-0.1 0 E(Ryd1
L
)
0.1 0.2 0.3
1
q&.1
-q_lsn
+ (Eb -a)b,S,_l-2bS,
S,‘, is polynomials E=Eb.
the derivative all of S,(E), in (26) must be calculated at energy
III ILLUSTRATIVE
Fig.1: Local density of states on a hydrogen atom adsorbed in a centered position on the (100) face of a bee transition metal with 6 d electrons.
EXAMPLE
The present method is much faster than when a direct numerical integration. Moreover, the density of states exhibits sharp resonances as explained in the introduction, we gain also a very large factor on the accuracy of the results, which is required in the determination of equilibrium distances or vibration frequencies for instance. We will illustrate now this point on a specific physical example. Let us consider the adsorption of hydrogen on the (100) face of bee transition metals. Thuault-Cytermann et al.[b], using a tight-binding method and a continued fraction expansion of the adatom density of states exact to the third moment and terminated as above,
The split-off state above the band is found at the energy Eb CO.313960 with a weight ~~~0.277360. The calculation of the zero-th and first order moments of the density of states (respectively N and ~1) should give 1 and al; this is a mean to check the accuracy of our method. In table I, we give the computing time and the absolute errors on N and ~1 obtained from our method and from a numerical integration using the trapezoidal rule with 2000 equal subintervals. The computations have been performed on an Amdahl V7 computer with single precision variables.
THE CONTINUED
404
FRACTION
EXPANSION IV
Computi g time (10% Dur method rapesoidal
340 rule
43500
OF A DENSITY
OF STATES
Vol. 50, No. 5
CONCLUSION
AN
10-5 -3.10-3
Table I. Comparison of computing accuracy of our method and a direct integration.
-3.10-5 10-3 time and numerical
One sees that a considerable gain of time and accuracy is obtained with our The use of double precision v_agriables method. would lead to an accuracy better than 10 .
We have shown that the integration of the imaginary part of the resolvent G(E) approximated by a continued fraction with n exact levels can be easily evaluated in terms of elementary integrals. The whole procedure is not completely tractable by hand as one must calculate the real and complex roots of a with the usual computer polynomial. However, libraries, this is quite easy and the result is much more accurate than the one we get using a direct numerical integration.
REFERENCES 1.
R. HAYDOCK, Sol.Stat.Phys. x,216(1980), M.J.KELLY, Sol.Stat.Phys.g,296(1980) and references therein, H. EHRENREICH, F. SEITZ, D. TURNBULL Editors (Academic Press New York), J.P. GASPARD Thesis, Orsay, France (1975).
2.
P.TURCHI,F. DUCASTELLE,G.TREGLIA, J.Phys.C (Solid State Physics) 2,2891(1982).
3.
A. TRIAS,M. KIWI,M. 2,1859(1983).
WEISMANN,Phys.Rev.
-B-
4.
K. MASUDA,Phys.Rev.
5.
RYZHIK,Table of I.S. GRADSHTEYN,I.W. integrals,series and products,Academic Press(New York London)1965 p.89.
6.
C. THUAULT-CYTERMANN, M.C. DESJONQUERES and D. SPANJAARD, J.Phys.C (Solid State Physits) x,5689(1983).
B-26,5968(1982).
ANALYTIC
INTEGRATION
Vo1.50,No.5,
pp.401-404,
OF THE CONTINUED
FRACTION
1984.
EXPANSION
003%1098/84 $3.00 + .OO Pergamon Press Ltd.
OF A DENSITY
OF STATES
G. ALLAN Laboratoire
d'6tudes
des Surfaces et Interfaces*,Institut Superieur d'Electronique 3,Rue Franqois Baes, F59046 LILLE Cedex (France)
du Nord
M.C. DESJONQUERES Centre
d'Etudes
Nucleaires de Saclay,Service de Physique des Atomes F91191 GIF SUR YVETTE Cedex (France)
et des Surfaces
D. SPANJAARD Laboratoire
de Physique Received
des Solides *,Universitk
November
de Paris-Sud,F91405
28, 1983, in revised
form February2
ORSAY Cedex (France)
1984byJ.Friedel**
We show that the continued fraction expansion of a density of states can be analytically integrated and expressed in terms of simple elementary integrals.
I INTRODUCTION
II SOME CONTINUED
Since a few years, the continued fraction method [l] has proved to be very efficient to calculate the electron or the phonon density for systems without 3-dimensional translation symmetry. It has been used near bulk defects (vacancies,dislocations,twin boundaries), near surfaces and for bulk disordered alloys. One of the main problem in this method was the termination of the continued fraction as one cannot calculate an infinite number of levels. Solutions have been given for a density of states with gaps [2] or for strongly oscillating recursion coefficients [3]. Here we consider the semi-elliptic termination assuming all the coefficients are equal after the n-th level. But we concentrate on the continued fraction integration. A numerical integration can always be done but sometimes one needs more accurate results. This is the case when one calculates the variation of the cohesive energy with interatomic distances and its first and second order derivatives in order to get the equilibrium distances and the interatomic force constant_s6 In such cases, an accuracy better than 10 must be reached. This is quite difficult using a numerical integration when the density of states shows very sharp peaks. Such peaks arise for example when a bound state is close to the band edge.
continued
FRACTION
PROPERTIES
Let US first consider fraction F,,(z):
a
truncated
1 F,(z)=
(1) bl x-al-
b
n-2
b x-an_l- -
n-l
s-a n It can be written
as[l]:
(2)
where Rn_l(~) and S,(z) are respectively (n-l)th and n-th order polynomials which satisfy the following recursion formulae: Rn= (z-an+l)*Rn-l-bn*Rn-2
(3.a)
Sn= (x-an)~Sn_l-bn_l*Sn_2
(3.b)
with
In the next section, we show that a simple analytical method can be used to obtain the exact result. Let us quote that a similar method has been recently developed by other authors[4] but limited to a continued fraction with few levels. Here we generalise to any order provided one is able to calculate all the real and complex roots of a real polynomial.
RO(z) and So(z) equal
to 1.
One can also easily Rn_l&,-l
- Rn-2*S, =b n-l+.(Rn-2&,_2
*
Laboratoire associe au Centre Nationalde la Recherche Scientifique. ** On behalf of the late Prof. A. Blandin's office 401
show that:
=b n-l bn_2
- Rn_3&.l)
.. b2 (R1*S1 - RO*S2)
=b n-l bn_2 . . .. . b2 bl
(4.a) (4.b) (4.c)
402
THE CONTINUED
FRACTION
EXPANSION
Such a truncated fraction is obtained for a finite cluster. The S,(E) polynomial roots are equal to some of the cluster eigenstates. One often embeds this cluster into a" infinite effective medium i.e one approximates the following levels to get for example the right band limits:
1 G,(z)=
(5) bl
OF A DENSITY
OF STATES
Vol. 50, No. 5
- b" b"_l....b2
bl
Im G"(E)=
JT
(13)
D" where Dn=2b Sz- 2b" S"_l S (E-a)+ 2b2S2 n n-l n
(14)
When E does not belong to the band, Im G,(E) is equal to zero save possibly for some discrete values of E corresponding to bound states.
z-alAs the highest coefficient of S" is equal to 1, the 2n-th coefficient of D" is equal to 2b(l-b"/b). This shows that D" is a 2n-th order polynomial if b is not equal to b and a (2n-l)-th order one if "bn is equal to be
b n-l x--a,_1s-an-b,,f(s) Let us recall here that the density states n(E) we are dealing with is simply
of
If we roots, we get
call
Ei
the
2n
(or
2n-1)
D
n
'i Im G"(E)= n(E)=- 1 Im G"(E) ll
1
J I(E-a)2-4bl
(6)
(15)
E-Ei
i where
From F"(z), to get G"(z), simply to formally replace bn_l by b,,f(z)/(z-a")).
one has b"_l/(l-
- b" b"_l
.. . b2 bl
ci=
(16) Dr!,(Ei)
Then one gets
and DA is the first-order
derivative
of D".
R"_l(x)-b"*R"_2(z)*f(z) (7)
G"(z) = S"(z)-b,,S"_l(z)"f(z)
Here we assume that all the recursion coefficients after the n-th level are equal to their infinite values a and b. In this case, f(z) is given by
Let us notice that only some of the real roots of D" correspond to bound states as some spurious states are roots of the conjugate denominator in (12) and are not poles of G". To integrate over energy the imaginary each part of G,.,(E), one now has to evaluate integral
1 f(z) =
(8) z-a-b,f(z)
EF
&zz
Ii=
dE /
or along the real axis (E-a)-S JR
where Ei is complex.
f(E)=
(9)
2b
One can show [5] that
where
RF Ii
R = )(E-a)2-4bl
s =
+1
S = +i s = -1
(17)
E-Ei
a-Ldb
E < a-2v'b a-2Jb < E < a+2Jc a+26 < E
-
E dE
+
(10)
&-
(11)
+ [4b-(Ei-a)2]
EF
(2a-Ei)
dE K
I EF
dE (18) (E-pi) fi
Within
Then one gets [2b Rn_l - b" R"_2(E-a)l+bn
=
E=a+2 &
R,_,S&
the usual transformation
sin t
(19)
(12)
G"(E)= [2b S" - b" Sn-1 (E-a)]+b"
S,_l SK
one easily gets Ii = 2 fi
This expression has a continuous art when th_e energy E lies within the imaginary band (a-2 Jp b < E < a+Zib) which is given by:
cos tF + (a-Ei)*(tF+v/2)+ tF
[4b-(Ei-a)2]
J
_n,2
dt (20) 26
sin t + (a-Ei)
Vol. 50, No. 5
THE CONTINUED
If we call Ji the integral uF
2
FRACTION
EXPANSION
in (20)
du
Ji = -
(21) a-Ei
I
4JTu u2 +-+1 a-Ei
-1
with u=tg(t/2). If ul and u2 are the roots one gets 2
1
Ji=-a-E.1
of the denominator,
UF-U2 (Log~-Log-
u2+1
UF-U1
u1+1
u2-u1
)
(22)
where Log is the complex logarithm. Let us notice that the two logarithms contraction in (22) can lead to 2v indetermination of the Ji imaginary part. Then adding the weights of the occupied bound states lying below the Fermi level, one can easily get the number of occupied states. The cohesive energy calculation can also be done using the same procedure
=--
I
E n(E) dE
al=0.0024611, bl=0.0798919,
a2=-0.0779382, b=0.0343536
a=-0.0751389,
n (El 16
A
14
10 E Im G,(E) dE
t 1 ci TI.
=
(24)
dE
1
Ii
+ Ei
1
blb2
‘0.
(25)
;
3
6 4
bn_l (26)
=
”
d
8
I
One has also to add the occupied bound if any multiplied by their state energies weights.Let us notice that these weights are the residue of G,(E) on the corresponding poles Eb. They can be easily computed and are given by
wJ.
have shown that a split-off state above the band is found for all investigated band-fillings Nb of the metal. On the other hand, a split-off state below the band exists only when Nb<6. When the hydrogen local density of states Nb'6, exhibits a sharp resonance just at the bottom of the band as seen in Fig. 1 corresponding to an hydrogen atom sitting on a centred site at The parameters of the equilibrium distance. calculation have been given in [6], and lead to the following values of the coefficients of the continued fraction (the Rydberg is taken as energy unit):
(23)
=F
1
II I
a-
403
OF STATES
12
RF
EC =
OF A DENSITY
b,S,$,_
2 0
-0.4
I
I
1
1
-0.3 -0.2-0.1 0 E(Ryd1
L
)
0.1 0.2 0.3
1
q&.1
-q_lsn
+ (Eb -a)b,S,_l-2bS,
S,‘, is polynomials E=Eb.
the derivative all of S,(E), in (26) must be calculated at energy
III ILLUSTRATIVE
Fig.1: Local density of states on a hydrogen atom adsorbed in a centered position on the (100) face of a bee transition metal with 6 d electrons.
EXAMPLE
The present method is much faster than when a direct numerical integration. Moreover, the density of states exhibits sharp resonances as explained in the introduction, we gain also a very large factor on the accuracy of the results, which is required in the determination of equilibrium distances or vibration frequencies for instance. We will illustrate now this point on a specific physical example. Let us consider the adsorption of hydrogen on the (100) face of bee transition metals. Thuault-Cytermann et al.[b], using a tight-binding method and a continued fraction expansion of the adatom density of states exact to the third moment and terminated as above,
The split-off state above the band is found at the energy Eb CO.313960 with a weight ~~~0.277360. The calculation of the zero-th and first order moments of the density of states (respectively N and ~1) should give 1 and al; this is a mean to check the accuracy of our method. In table I, we give the computing time and the absolute errors on N and ~1 obtained from our method and from a numerical integration using the trapezoidal rule with 2000 equal subintervals. The computations have been performed on an Amdahl V7 computer with single precision variables.
THE CONTINUED
404
FRACTION
EXPANSION IV
Computi g time (10% Dur method rapesoidal
340 rule
43500
OF A DENSITY
OF STATES
Vol. 50, No. 5
CONCLUSION
AN
10-5 -3.10-3
Table I. Comparison of computing accuracy of our method and a direct integration.
-3.10-5 10-3 time and numerical
One sees that a considerable gain of time and accuracy is obtained with our The use of double precision v_agriables method. would lead to an accuracy better than 10 .
We have shown that the integration of the imaginary part of the resolvent G(E) approximated by a continued fraction with n exact levels can be easily evaluated in terms of elementary integrals. The whole procedure is not completely tractable by hand as one must calculate the real and complex roots of a with the usual computer polynomial. However, libraries, this is quite easy and the result is much more accurate than the one we get using a direct numerical integration.
REFERENCES 1.
R. HAYDOCK, Sol.Stat.Phys. x,216(1980), M.J.KELLY, Sol.Stat.Phys.g,296(1980) and references therein, H. EHRENREICH, F. SEITZ, D. TURNBULL Editors (Academic Press New York), J.P. GASPARD Thesis, Orsay, France (1975).
2.
P.TURCHI,F. DUCASTELLE,G.TREGLIA, J.Phys.C (Solid State Physics) 2,2891(1982).
3.
A. TRIAS,M. KIWI,M. 2,1859(1983).
WEISMANN,Phys.Rev.
-B-
4.
K. MASUDA,Phys.Rev.
5.
RYZHIK,Table of I.S. GRADSHTEYN,I.W. integrals,series and products,Academic Press(New York London)1965 p.89.
6.
C. THUAULT-CYTERMANN, M.C. DESJONQUERES and D. SPANJAARD, J.Phys.C (Solid State Physits) x,5689(1983).
B-26,5968(1982).