# Full text of "Spatial decay of the single-particle density matrix in insulators: analytic results in two and three dimensions"

## See other formats

Spatial decay of the single-particle density matrix in insulators: analytic results in two and three dimensions (N O o (N > cn o I O o S. N. Taraskin,!'! D. A. Drabold,2^[] and S. R. EUiotti^l] 'Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 lEW, UK ^Department of Physics and Astronomy, Ohio University, Athens, Ohio 45701 (Dated: February 1, 2008) Analytic results for the asymptotic decay of the electron density matrix in insulators have been obtained in all three dimensions {D = 1 — 3) for a tight-binding model defined on a simple cubic lattice. The anisotropic decay length is shown to be dependent on the energy parameters of the model. The existence of the power-law prefactor, oc r"^''^, is demonstrated. PACS numbers: 71.15.Ap, 71.20.-b, 71.15.-m In a phenomenological classical approach to atomic dy- namics, a very local picture of interatomic interactions is often highly appropriate. The simplest example is the use of Keating "springs" to describe small atomic oscil- lations in solids [Q. In quantum mechanics (QM), the situation is superficially quite different. If are single- particle wave functions (for example, the Kohn-Sham or- bitals of density-functional theory), then the electronic contributions to the total energy, forces or the dynamical matrix can all be obtained from the single-particle den- sity matrix (DM): p(r,r') = Eocc V'l* (r)V'i(r') §. For a condensed-matter system, virtually all of the Hamilto- nian eigenstates ipi are extended and oscillatory through- out the volume of the system, except some of the states near the band edges in the case of disordered systems. To the extent that p(r,r') is built from objects that are extended in real space, it is a nontrivial feature of the quantum mechanics of extended systems that in fact the DM can be localized, even exponentially so in insulators ii- The possibility of a local formulation of QM goes back at least to Wannier (5 , whose " Wannier functions" decay in real space. Kohn [ 3| extended this work, and empha- sized the principle of " nearsightedness" , that is the de- pendence of local properties, such as forces, on the local environment. The possibility of real-space local formu- lations of QM has inspired vigorous activity on efficient "order-N" methods for electronic structure In early pioneering work, Kohn |Q showed that the DM and Wannier functions for a ID centrosymmctric two- band model decay exponentially in insulators: p(r, 0) cx exp(— A|r|), where A cx E^ip 0, although Ismail-Beigi and Arias 1^ find instead that A cx Eg^p as Eg^p 0. Dcs Cloizeaux proved that the DM decays exponentially quite generally in insulators [^. Recently, He and Vanderbilt pof demonstrated that in ID there is in fact a power- law prefactor cx r~^/^, and that the characteristic inverse decay length. A, of both the DM and Wannier functions is universal, although the prefactor is different in each case. In this Letter, we report analytic results for the spa- tial decay of the DM for a two-band tight-binding (TB) Hamiltonian for all three dimensions, D — 1 — 3. The Hamiltonian we use is similar in spirit to that of Harri- son's bond orbital model [0 and is the minimal model which contains the basic features of an insulator. Both the power-law prefactor and the exponential spatial de- cay of the DM are found in the asymptotic regime. Thus, we reproduce the results of Refs. [0, |o| , and extend them to 2D and 3D. We derive the dependence of the inverse decay length. A, in terms of the fundamental energy pa- rameters of the model. We further show how, in certain cases, the exponential localization arises from the trun- cation of the sum defining p. Precise numerical results support the analytic asymptotic behaviour in all three dimensions. Our model is of centrosymmctric nearest-neighbor TB form on a simple cubic (s.c.) lattice described by the Hamiltonian: (1) Here, each site i is characterized by two bare orthogonal electronic states, (/i = 1,2) with bare energies (Ae = El — £2 > for definiteness) , which together form the orthonormal bare basis set. These may be interpreted as bonding and antibonding states centered at each lat- tice site. The hopping integrals, t^^, between similar states (the same n) on different nearest-neighbour sites j (to i) transform the on-site energy levels into bands, while the hopping integrals, t^^i , between different states (/i 7^ ^') on different sites are responsible for inter-band hybridization. For this model, the eigenstates |k,7) (k is the wavevector and 7 = 1,2 indexes the bands) and dispersion are analytically available and the spectrum exhibits the main features of a semiconductor /insulator (i.e. two bands separated by a gap for a certain range of parameters). We have written the density operator for the Hamil- tonian given by Eq. (|l|), first, in the eigenvector basis, |k, 7), then expanded the eigenstates in the site basis \ifi) and evaluated the density matrix elements, Pifi.j^' 2 -20 ^ -40 -60 3 2 ;^ 1 -1 (a) — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — (b) ^ s 1 \ 1 1 1 1 ' (c) « S » S 9 9 fl ? g g 5 = S ? ^«88iiiSiifiiiii,,. _ □ioiiSSSfifiBSSSHBSnoo o ^ □ □ □ Q g Q a ^ S Q 6 O E3 Q Q S a D Q 10 20 30 40 10 20 30 40 FIG. 1: The dependence on distance of the logarithm of the DM llnp^^l (a,b) and of the effective inverse decay length, A — — ln|p^^|/cr (c,d) in the 2D case (the prefactor is in- cluded in the definition of A for easier comparison of analyti- cal results with numerical ones). In (a) and (c), the DM has been calculated on the lattice sites situated around the main diagonal [11] and characterized by polar angles 4> € [lO^jSO"], while in (b) and (d), Vy — and i-'x varies (in the [10] direc- tion). The circles represent exact numerical results obtained using Eq. ^ or Eq. (^), while the crosses are the results of the analytical expressions (§)-(^ (they are practically indis- tinguishable in (a) and (c)). In (c), the different analytical contributions to the effective inverse decay length are shown by the solid line (the ln(2yl)-term in Eq. (^)), squares (the second term in Eq. (^), triangles (the sum of the last two terms in Eq. (^) and diamonds (the prefactor in Eq. (^). In (d), the solid line is the ln(2AV-contribution to the effective inverse decay length (see Eq. (|6|)), the squares represent the contribution of (— In | J^^ (z/x/LA)|/j^x) to A and the diamonds signify the prefactor in Eq. in the site basis. For fixed site indices i and j, the DM contains four elements, which, in the case of a semicon- ductor at zero temperature, obey the following relations: Pit^,m = ^ii s-'^d piij2 = Pi2ji- Below we investigate the decay properties only for the off-diagonal elements, Pii,j2 (the analysis and results for the diagonal elements Pifj.,jfi are very similar to those for Piij2 presented below). The expression for pa,j2 = p{^ij) can be conveniently represented via an integral in reciprocal space over the first Brillouin zone: '5k (2) where {ai^a;, avy^ avz} — — is the connection vec- tor (with the unit-cell parameter a taken hereafter to be unity, and (e.g. a = x,y, z m 3D) being an inte- ger), S'k = (1/2) ^^-^-^ exp (ik • Tij) is the structure fac- tor and = (Ae -f 2(tii — t22)Sk)/Ati2 is the only en- ergy parameter of the model which, for the symmetric case (^11 = <22 = t), analyzed below for simplicity, is k- independent i.e. Ak = A = Ae/4ti2. In the model without interband hybridization (ti2 = or A — > cxd), the DM given by Eq. (||) is obviously pij = 0. In order to study the spatial decay of the density matrix, we should take the integral in Eq. ^ and investigate its depen- dence on the lattice indices i^a ■ This can be done analyt- ically in the asymptotic regime for all or some i/q, 1, if the energy parameter of the model obeys the inequality, A > D, which is the case for weak interband hybridiza- tion in semiconductors. Eq. (|^) has been derived without making use of the particular symmetry type of the underlying lattice, and is valid for any primitive lattice (one atom per unit cell) . The structure factor, S'k, reflects the differences between various lattices. The DM given by Eq. (^) can be eas- ily calculated numerically for any reasonable choice of parameters, and the results of such calculations support the main conclusions derived below analytically for the s.c. lattice, subject to some additional restrictions (sym- metric case, til = ^22, and weak hybridization, A > D). To proceed with the evaluation of the DM given by Eq. (^, first we expand in a series the denominator of the integrand (S^/A < 1, if A > D), separate the vari- ables and then evaluate all D integrals making use of the orthogonality of the cos(fcct;/Q.)-functions involved. This results in the following general exact form for the DM for all three dimensions: P^c (-ir (4/1)2^- -1 /l^ fe=0 (2fc')' {2k' + i)j:d h^y.{k - m + Vy)!] be- but that for 31?, to)! (to -|- fx)!(m + (4:A)''{k' (3) if <T = h'a is odd, and zero otherwise. Here k' ~ k + V = k + {a ~ l)/2. El = l/[fc!(fc + v^y.] and ^2,3 are hypergeometric finite sums, the one for 21?, ^2'= ELol/H(fc-m)!(m ing expressible in closed form |lj ^3^j:t=oi'^rn + u^ + i.yy./[mlik t'y)!(TO + I'-x + ^'y)!(^ — m + J^z)!], cannot so be written, in principle |l^ . The orthogonality of the cos(fcct^'a)- functions, being a consequence of the Fermi level lying in the gap (only in this case are the upper limits of the integrals in Eq. (H) equal to tt), is an important feature in deriving Eq. (0). The original series is truncated for small indices which results in the appearance of the ex- ponential factor, (4A)~2''-i^ before the sum in Eq. (|^). This factor then appreciably contributes to the energy (A)-dcpendent part of the decay length of the DM (al- though see Fig. |^(e) where this is not the case). Further analysis is based on the use of Stirling's ap- proximation for Va ^ 1 which allows us to evaluate asymptotically both the internal finite sum S3 in the 3-D case and the external infinite sums over k in all dimen- sions. For example, if all ^ 1, then keeping the lead- ing term in the asymptotic expansions, we can write the 3 -10 -20 -30 -40 -50 3 (a). 1 (b) - i 8 is (c) - ** — 1 — 1 — i — 1 — 1 — (d) — ' 1 1 1 ' . » (e) •"i, i - "■■■■18 ••:tiiiiiiii< - 8 a S „ _ □□nanSSfififlfifi 10 20 v.+v„+v. 10 20 v,+v„ 10 20 V, 30 FIG. 2: The dependence on distance of | In p^^ \ (a-c) and of A = — ln|p^^|/(T (d-f) in the 3-D case. In (a) and (d), the DM has been calculated on the lattice sites situated around the main diagonal [111] for solid angles defined by the polar angle (j) G [20^ 70°] and azimuthal angle 6 £ [20°, 70% in (b) and (e), ly^ — and and i^y vary around the plane diag- onal and are characterized by the polar angle (f) G [20°, 70°] (around the [110] direction); and in (c) and (f), = fy = and varies (in the [100] direction). All the symbols have the same meaning as in Fig. |l| except in (e), where the solid circles represent the contribution to A due to the exponential term in Eq. (^, the squares and triangles are due to J^^ and "^fx+fyi respectively, while the diamonds signify the prefactor in Eq. d). following expression for the DM in all three dimensions: p.„^(-l)''(27r)-^/Mn^.J exp[-aA] , (4) where the inverse decay length, X{A, Ua), is given by A(A^^a) = ln(2A) + ^ + ^-lna , (5) with cr_i = Y.'^a^ and CTin = Va\nva. As fol- lows from Eqs. (0)-(||) and Figs. 0(a,c) and §(a,d), the DM in this model of an insulator: (i) exponentially de- cays with the distance-related parameter cr; (ii) has a power-law prefactor which is proportional to v~^l'^\ (iii) is anisotropic; and (iv) the inverse decay length of the DM depends on the energy parameters of the model. The behaviour with effective distance of In \pu^ \ and the effec- tive inverse decay length, A = — In l/ci, are shown in Figs. |l|(a,c) for the 2D case and in Figs. ||(a,d) for the 31? case, respectively. The different points correspond to all the lattice nodes within a relatively wide solid (po- lar in 2D) angle around the main diagonal. Examina- tion of the results shows that the slowest decay of the DM occurs along the main diagonal, i.e. the [111] direc- tion in the 313 case. The analytical expression ^ allows the different contributions to the effective decay length FIG. 3: The dependence of the inverse decay length A (a) and the effective inverse decay length A (b) on the energy- dependent parameter A of the model in 3-D. The data in (a) given by the thick solid and long dashed lines were obtained by linear regression from the slope of In |p| evaluated numerically using Eqs. (^ and (^ for the [111] direction, respectively. The thin solid line in (a) represents Eq. (^) (the difference between the thin and thick solid lines is due to the prefactor in Eq. (^) while the short dashed line shows the ln(274)-contribution to A. The solid line in (b) represents the analytical result ac- cording to Eqs. (^)-(^), while the long dashed line is obtained by exact numerical evaluation of Eq. (0) or Eq. (^. The short dashed line in (b) shows the ln(2A)^contribution, while the dot-dashed line gives the 1/yl '^-contribution. The DM in (b) has been calculated on the particular representative node characterized by v^^ = Vy = = 9. to be estimated. It is clearly seen in Figs. |](c) and §(d) that the major A-dependent contribution comes from the ln(2^)-term in Eq. (||), while the role of the second term, '-^ \/A^, in Eq. (||) is not important at all, especially in comparison with the significant A-independent contribu- tion from the last two terms in Eq. (||) (see Fig. ||(b)). In the case when one of the indices, e.g. ^ 1, but the other ones are small and finite (-D > 2), Mil exp[-i/,ln(2.4)] n^-. (^) , (6) where J stands for the Bessel function, the asymptotic form for which, with large argument | p^ , must be used in Eq. (|6|). Bearing in mind that J^„(i/xM) « {v^/Ay'^/'^ for S> 1, it is easy to see that the power-law prefactor, (X v~^/'^^ is restored in Eq. (^). The results for this particular case are presented in Figs. |l|(b,d) and H(c,f). Again, the contribution of the term ln(2v4) in the inverse decay length is dominant, and the inverse decay length is now larger than for decay along the main diagonal (cf. e.g. Figs, ga) and§(c)). Finally, if t'x, J^y 2> 1 and is finite (in 3-D), 2'KV^Vy exp 4 A (7) where v± = Vy^:^ Vy and, for both Bessel functions, known asymptotic expressions (for the argument in 3i^^ and for the argument and index in J^^ ) must be used [ p2[ . These asymptotics for both Bessel functions restore the v~^/'^ prefactor for the DM in the same manner as for Eq. (6). The energy (A)-dependent contribution to A (see Fig. 2(e)), in this case, comes mainly from the asymp- totic expression for J^^ which is of exponential type if v+/{A,Jir^) < 1 (and of cos- type otherwise) p2| . Another important point to discuss is the dependence of the typical inverse decay length. A, on the energy pa- rameters of the model. The only energy-dependent pa- rameter in this model is the parameter A which is sim- ply related to the spectral characteristics of the model, A = A£/4ii2 = DAEopt/iP^E^), where ASb = ADt is the band width, AE'opt = Ae is the optical (direct) band gap and the ratio of hopping integrals, (3 = ti2/t, which can be connected to the thermal gap, AE'th = [A£;2p^ + {l3AEbfY/^ - AEb. The spectral parameters (Ai?b, Ai?opt and Ai^th) are not independent values be- cause the hopping integrals t^^' in the Hamiltonian (|l|) can depend on the orbital energies and obtaining such a dependence, e.g. ii2(Ae), requires more accurate anal- ysis which is not necessary for our approach (that is why we plot Xvs A but not the direct gap width Ae in Fig. ||). However, the limiting case, Ae — > 0, can be investigated, at least numerically (analytical results are valid only for A > D) under the reasonable assumption that ti2 ^ as Ae ^ 0, and the results are presented in Fig. ||(a), confirming the linear scaling of the inverse decay length, Aoc^cx A£for^->0 |]. This model can be analysed analytically only for a s.c. lattice and thus, strictly speaking, we stick to the atomic orbital representation. However, numerics show that the results (exponential decay of the density ma- trix and variation of the decay length with parameters of the Hamiltonian) are qualitatively the same for differ- ent underlying lattices. This is not surprising because the asymptotic behaviour at large distances hardly de- pends on the local structural details (the local features of the lattice can influence the value of A by a factor of ~ 2; see Figs. Thus, we believe that our results, at least qualitatively, are general. Indeed, the inverse de- cay length decreases with decreasing A (see Fig. ||), thus reflecting the known delocalization tendency of the DM with increasing metallicity (in the bonding/antibonding representation A ~ ^^2/^1, being the ratio of the bonding-antibonding splitting energy and the metallic bandwidth energy). Moreover, the analytical expressions for the decay length give correct order-of-magnitude es- pressions for the spatial decay of the density matrix de- cay in insulators. The results have been obtained for a tight-binding Hamiltonian defined on the s.c. lattice having two orthogonal bare states on each node. Two additional assumptions, namely symmetric bands (equal in-band transfer integrals) and relatively weak intraband hybridization, allow us to derive exact asymptotic results for the density matrix in all three dimensions. The main features of the analytic solution are: exponential spa- tial decay of the density matrix, anisotropy of the decay length, the dependence of the decay length on the energy parameters of the model, and the existence of a power- law prefactor (cx r~-°/^). All the analytic asymptotic re- sults have been supported by precise numerical solutions of the problem. We are grateful to Peter Paule for giving us access to the finite hypergeometric summation codes. We thank L. He and D. Vanderbilt for helpful comments on the manuscript. DAD acknowledges support from the NSF under grants DMR 0081006 and DMR 0074624. Electronic address: sntlOOO@cus.cam.ac.uk Electronic address: irabold@ohio.edu; Work imple- mented while on le ave at Trinity Colleg e, Cambridge [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] Electronic address: srel@cus.cam.ac.uk P. N. Keating, Phys. Rev. 145, 637 (1966). S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999). W. Kohn and R. J. Onffroy, Phys. Rev. B 8, 2485 (1973). J. J. Rehr and W. Kohn, Phys. Rev. B 10, 448 (1974). G. H. Wannier, Phys. Rev. 52, 191 (1937). W. Kohn, Phys. Rev. W. Kohn, Phys. Rev. timates of A for Si and C |14 [15] To conclude, we have derived analytic asymptotic ex- Lett. 76, 3168 (1996). 115, 809 (1959). S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett. 82, 2127 (1999). J. D. Cloizeaux, Phys. Rev. 135, A698 (1964). L. He and D. Vanderbilt, Phys. Rev. Lett. 86, 5341 (2001). W. Harrison, Electronic structure and the properties of solids: the physics of the chemical bond (Freeman, San Francisco, 1980). G. N. Watson, Theory of Bessel Functions (Cambridge University Press, Cambridge, 1966), 2nd ed. M. Petkovsek, H. S. Wilf, and D. Zeilberger, A=B (A. K. Peters, Wellesley, 1997). Using the relation A — V2/V1 (relevant for the bond- ing/antibonding representation) and the values of pa- rameters given on p. 90 of Ref. [Q, we can estimate A « ln{2A) = 1.22 (1.87) for diamond Si (C), which are rather close to the values A ~ 1.15 (2.06) obtain by real- istic calculations in Ref. These estimates are coarse since the diamond lattice does not have bond orbitals on a s.c. lattice and we use only the simplest approximation for A. X. Zhang and D. A. Drabold, Phys. Rev. B 63, 233109 (2001).