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

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).