Approach to calculation of long-range Coulomb interaction matrix elements in ion crystals Текст научной статьи по специальности «Физика»

ISSN 2072-5981

Volume ІЗ, No. 2 pages 27-35 2011

http://mrsej.ksu.ru

Established and published by Kazan University Sponsored by International Society of Magnetic

Resonance (ISMAR) Registered by Russian Federation Committee on Press,

August 2, 1996 First Issue was appeared at July 25, 1997

© Kazan Federal University (KFU)

”Magnetic Resonance in Solids. Electronic Journal' (MRSey) is a

peer-reviewed, all electronic journal, publishing articles which meet the highest standards of scientific quality in the field of basic research of a magnetic resonance in solids and related phenomena. MRSey is free for the authors (no page charges) as well as for the readers (no subscription fee). The language of MRSey is English. All exchanges of information will take place via Internet. Articles are submitted in electronic form and the refereeing process uses electronic mail. All accepted articles are immediately published by being made publicly available by Internet (Editor@ksu.ru).

Editors-in-Chief Jean Jeener (Universite Libre de Bruxelles, Brussels) Boris Kochelaev (KFU, Kazan) Raymond Orbach (University of California, Riverside)

Executive Editor

Yurii Proshin (KFU, Kazan) Editor@ ksu.ru

Editors

Vadim Atsarkin (Institute of Radio Engineering and Electronics, Moscow) Detlef Brinkmann (University of Zurich, Zurich) Yurij Bunkov (CNRS, Grenoble) John Drumheller (Montana State University, Bozeman) Mikhail Eremin (KFU, Kazan) Yoshio Kitaoka (Osaka University,

Osaka)

Boris Malkin (KFU, Kazan) Haruhiko Suzuki (Kanazawa University, Kanazava) Murat Tagirov (KFU, Kazan)

*

In Kazan University the Electron Paramagnetic Resonance (EPR) was discovered by Zavoisky E.K. in 1944.

Approach to calculation of long-range Coulomb interaction matrix

elements in ion crystals

O.A. Anikeenok

Kazan Federal University, Kremlevskaya 18, Kazan, 420008, Russia

E-mail: Oleg.Anikeenok@ksu.ru

(Received November 18, 2011; revised December 22, 2011; accepted December 25, 2011)

Expressions for the calculation of matrix elements of the Coulomb interaction of the electron with an infinite crystal lattice in ion crystals have been obtained. In the general case, matrix elements are calculated on the orbitals belonging to ions of different sites. Gaussian-type orbitals are used in the calculations. All expressions are absolutely and rather rapidly converging series in the space of the inverse lattice vectors. At present the value of this interaction being only one-center is estimated by the Madelung constant, in other words, by calculating the electrostatic potential in the lattice site, i.e., in a point.

PACS: 61.50.Ah, 61.72.S

Keywords: localized orbitals, infinite crystal lattice, Madelung energy, orbital energy

1. Introduction

The method of secondary quantization with a partly nonorthogonal one-particle basis was developed in [1-4]. Its application to admixture centers [4], when the overlap integrals of ion orbitals are sufficiently small, allows one to build mathematically correct expressions for the perturbation-theory series within virtual and real charge transfer processes. To calculate physical characteristics of an admixture center ab initio, it is necessary to separate an U-region around the admixture ion where all interactions should be written exactly and the other part of the crystal, which is taken into account in the ion approximation. It is convenient to present the charge of nuclei Z of ions in the U-region as Z = q + n + m [2], where q is the charge of the ion in the host crystal, n is the number of electrons in ions in the considered configuration, m is the deviation of the ion charge in the host crystal for the considered configuration. Such a partition allows one to rather simply group all interactions in the U-region, at least in ion crystals, with respect to their order of magnitude, e.g., when calculating the crystal field parameters or amplitudes of the transition of electrons between ions [2-4]. Note that the value of the ion charge - q in the U-region is a definite number and not a fitting parameter. In this approach it is necessary to calculate two-center matrix elements (TME) and one-center matrix elements (OME) of the long-range (LR) Coulomb interaction of an electron with an infinite crystal lattice.

At present the value of this interaction being only one-center is estimated by the Madelung constant, in other words, by calculating the electrostatic potential in the lattice site, i.e., in a point [5-8]. The exclusion is [9], where the energy of the Coulomb interaction of the 5-electron with the infinite crystal lattice was calculated.

The article below is arranged as follows. General expressions for the calculation of TME of the LR Coulomb interaction of the electron with the infinite crystal lattice are derived in the second section. Formulas necessary for the calculation of matrix elements on 5, p, d, f orbitals of ions are given as well. On the basis of the expressions for TME of the LR Coulomb interaction, the expressions for OME of the LR Coulomb interaction are presented and particular cases are considered.

2. General part

Let the radial part Rnl ion orbital \ynlm (r) where the electron is located has the form of the Gaussian basis expansion (GTO)

Rni = Z exP (-ar2 )• (1)

Let the first ion be located at the site with the radius-vector R0 + r;., R0 = 0 . The second ion is in the site with the radius-vector rb . The charge qp is in the site with the radius-vector Rn + rp . Let us denote r0 = rb - r;., R = Rn -( - rp) ,where Rn is the vector of the n - th crystal unit cell, r;. and rp are the vectors of the ions of the unit cell, g, g are the quantum numbers of ion orbitals.

The matrix element of the Coulomb interaction between the electron and the charge qp calculated on the wave functions of the first and second ions has the form

(Mr - rj)

||Mr - rb )) = (Mr )|-

qp

r - Rn -(rj - rp )

-(R n + rp)

We define the F (n1n2n3) functions as follows

2n F (nn2n3 ) = ~qp Zaibk jexP (-ar2)) yn z"3 |r - R|-1 exP -Pk (r - r0)

n\r -(rb - r j

(2)

z.

(3)

For example, the matrix element calculated on the |pz (r ) functions has the form

(pz (r)-r-pR|pz (r -r0 ft = ^[F(002)-z0F(001)] .

(4)

Let us further present the F(n1n2n3) functions in a form convenient for calculations. To this end, we perform transformations used in [2]

-Lr, = j= jdvexp[-(r - R)

2 v2

(5)

2 (

V =-

3

^2 du

dv =----------

a

(6)

'ik

where aik = a + ak . After transformations (5) and (6), the integration over x, y, z in (3) is reduced to table integrals [10] and the following expression is obtained for the F(n1n2n3) functions

F (n1n2 n3 ) = -qp Z aibk — lj dun

V aik J 0 s=1

(1 - U2 ) ns-Z [((s - Cs )u2 ]

=0 (4a-k )msms! ws=0 w

!(ns - 2ms - ws )!

X exp

-aik (R - c)

22

exp

afa

V aik J

where [ns /2] is the integer part of a number in brackets; x01 = x0, x02 = y0, x03 = z0 denote the

coordinates of the vector r0; cs = pkx0s / aik denote the components of the vector c;

Rj = Rx, R2 = Ry, R3 = Rz denote the coordinates of the vector R . After multiplication of three curly brackets in (7), the multipliers of the form

{Rx -ci) ( -C2 ) (Rz -c3 )3 expf-aik (R -c)

\2 2 ' u

(8)

where ns = ns - 2ms - ws, are presented in each term of the resulting sum. According to the above

notation R = R n - R - rp) .We introduce a vector r = r;. - rp . Since the ion positions in the unit cell

are rather arbitrary, we consider that the vector r is defined in all points of the unit cell and let us introduce the D (r, nj, n2, n3) function whose domain is the unit cell

D (r, nJ, K n3 ) =Z n|(R

-f

n -2m -w

exp

)

(9)

where xj = x, x2 = y, x3 = z are the coordinates of the vector r; Rnj = Rrx., Rn2 = Rny, Rn3 = Rnz denote the coordinates of the vector Rn. The D(r, nj,n2,n3) function is a periodic function of r with the period of the crystal unit cell (see, e.g., [11]) and the same as in [11], when the Fourier coefficients D(g, n1,n2,n3) of the D(r, n1,n2,n3) functions are found, integration over the unit cell can be reduced to the integration over the whole space. Thus

D (i% n^, n2, n3 ) = Z D (g, ^, n2, n3 )exp [i (gr )~\-.

(10)

D (g,«!,n2,n3

xnynzn exp{-|^aikr2u2 + i(gr) }dxdydz exp[-i(gc)],

(11)

where vc is the volume of the unit cell, g is the vector of the inverse lattice, i is the imaginary unit. By performing integration in (11), we obtain

D ( ^, n2, n3 )= —

vc

V aik

3

2 — ~ if 1

ns i| —

2

n

s=1

/ . \ ns -2h / \ ns-h

(-igs ) S S f 1 ^

h i(n, - 2hs )i

Vaiku J

X| U3 Iexp

2

4aku

exp [-i (gc )^,

'iku J

(12)

where g1 = gx, g2 = gy, g3 = gz are the coordinates of the vector g. Let us further introduce the Fjb (nln2n3) functions. To this end, we substitute R = Rn - r into expression (7) and perform summation over vectors Rn with the use of expressions (9)-(12), then we substitute r = r;. - rp into the obtained expression and perform summation over vectors rp of the whole unit cell. These calculations are bulky but rather simple and as a result, we obtain

g

c

Fjb (n1n2n3 ) =-n1 in2 in3 !Za>bk

V aik J

exp

aibk r2 ' *0

v a ,

IIdu V u

XZ fjb (nl, gx )fjb (n2 , gy )fjb (n3, gz )exp

4aiku

V ™iku J

Z qp exp i'g (rj- rp) exp [-4

V p

where

( 2 V

(1 - u ) n-2m 1

fjb (n, gs ) = Z ( \m i Z --j

m=0 ( 4a;-k) m i w=0 i

-igs

V 2aik J

(-1)h

Z

h=0 hi(n - 2m - w - 2h )i

u*/a,

(14)

Then we obtain the explicit expressions of the fjb (n, gs) functions for the values n = 0,1, 2,3, 4

Pkx0S

-i^T~, fb (0, gs ) = 1, fb (1, gs ) = Z, fb (2, gs )=-Zs2 +

aik 2aik

fjb <3-g- )=iZ* + itZ- ■ fb( *)=+ 3k

(15)

It is seen that the fjb (n, gs) functions do not depend on the parameter u . Below we will call the sum over the index p in parentheses in (13) as a structural factor Gj (g). It is presented as follows

Gj (g) = g(1) (g) +iG(2) (g) = Z qp exp ['g (rj - rp)

G(1)(g )= cos (j )F1 (g ) + sin (j )F2 (g ) G^2)(g ) = sin (j )F1(g )-cos (j )F2 (g )

F1 (g ) = Z qpcos (grp), F2 (g ) = Z qpsin (grp).

p p By substituting (16) in (13) we obtain for the Fjb (n1n2n3) functions

(16)

Fjb (n1n2n3 ) =-n1 in2 in3iZaibk

i.k

V aik J

exp

f aibk r 2 A 1 " r0

v a

(17)

XZ fjb ( % gx ) )jb (n2, gy )fjb (^, gz ) exp

2

V 4aiku J

Gj (g)exp[-i (gc)].

By performing integration over the parameter u in (17) we obtain the final expression for the Fjb (n1n2n3) functions

3 2n2

Fjb (n1n2n3 ) =-n1 in2 in3 !Z aibk

vc ~T

V aik J

exp

aib

____L r 2

*0

V aik J

(18)

XZ fjb («1, gx )fjb (n2 , gy )fjb (n3, gz ) exp

2

V 4aik J

exp [-i (gc )].

Now we introduce the operator HLR of the Coulomb interaction between the electron and the infinite crystal lattice

g

5

g

3

g

W (r - rj )HLR | W? (r - rb f = (w? (r - rj )

- Z

n, p

■|W'(r - rb ) .

Thus, the matrix elements of the operator HLR calculated on the orbitals of the lattice ions are expressed in terms of the Fjb (n1n2n3) functions. For example, we obtain for the |pz (rf and \pz (r - r0f orbitals

{pz (r )| Hlr\pz (r - r0 f = - [ Fjb (002)- ^ F]b (001)] .

(19)

One-center matrix elements

One-center matrix elements can be expressed in terms of the Fj (n1n2n3) functions. To find these functions, it is sufficient to take r0 = 0 in (18) and we obtain

3 2n2

Fj (n1n2 n3 )=----------------------------------------------------n1in2in3iZ aibk - Z f (nl, gx )f (n2, gy )f (n3, gz)

Vc Tt Vaik J "

2

-exp

V 4aik J

. (20)

The f (n, gs) functions are obtained from expressions (14) for the fjb (n, gs) functions, if only the terms with the index w = 0 are kept there. To obtain the f (n, gs) functions for the values n = 0,1,2,3,4, it is necessary to assume x0s = 0 in the expression Zs in formulas (15). The f (n, gs) functions for the values n = 5, 6 are given below

zs =-i^, f (5, gs ) = ^ zs5 +-L- zs3 +- 1

2—

5i s 24— s 32— s’

1

1

1

1

(21)

f (6,gs ) = — z6 +--------z^ +------— z2 +

Ss' 6i s 96— s 64— s 384ai

3

Thus, all matrix elements for the s, p, d, f orbitals of ions are determined.

Particular cases

The matrix element of the operator HLR on the |pz (r f wave functions and at F2(g) = 0, r0 = (x0,0,0) according to (18) and (19) is expressed in terms of the function

1

Fb (°02) = T-17TZab

2vn1'

f

V aik J

exp

aA

'iHk v2

V —k

x0pk

\

xj-

V aik J

nx + y,ny + z3nz

Z

ax, ny,

F(

22 n f nz 1 1

nx, ny , nz,

22 / n \

+

V b J

+

(22)

X exp

ak

22 f ny 1

+

V b J

+

3

5

where F1(g) = Fl[nx, ny, nz), g = {lnnx / a,2nny / b, 2nnz / c), a, b, c are the lattice constants. The quantities x., y., zj, x0 in the argument of the cosine are given in relative units.

The one-center diagonal matrix elements calculated with the help of (20) are the energy E. of the Coulomb interaction between the orbital electron located in the site r. and the crystal lattice, including the interaction with the charge q. located in the same site. At the same time, it is natural to refer the

interaction between this electron and its nucleus to the Hartree-Fock energy [2-4]. In this case it is necessary to exclude the interaction between the considered electron and the charge q. from the

energy E., i.e., to introduce the quantity E. (/) defined as follows (in [9] it is denoted as E.)

Ej (/) = Ej - q.Eo, (23)

where E0 is the interaction between the considered electron and the unit positive charge at the same site, / is the orbital. The energy E0 can be expressed in terms of the F0 (nln2n3) functions, which are determined by formula (3), if r0 = 0, R = 0, qp = 1. If it is taken into account that F0 (n1n2 n3 ) 0 only when n1, n2, n3 are even numbers, then, according to (7), we obtain

F0 (n!n2 n3 )=-

n1 !n2!n3!

m!

;Z

( i V+m

a a,,

V aik J

m1!m2!m3! 2m (2m +1)!! i k where n1 = 2m1, n2 = 2m2, n3 = 2m3, m = m1 + m2 + m3.

We write below the expressions F0 (n1n2n3) ,which are used further

(24)

F0 (000) —X

a,ak

i,k

V aik J

, F0 (200) = F0 (020) = F0 (002) = -3Xa<ak

i,k

V aik J

F0 (220) = -15Z

1 V’

a a,,

V aik J

F0 (400) = F, (040) = ^-£

5 i k

( 1 v3

a a

V aik J

(25)

According to (20) and (25), the energy of the interaction between the arbitrary 5-orbital E. (5) and the lattice determined by (23) can be written as

i,k

V aik J

4^Z G.(g) (

—Z^^exP

.2 A

4a.-

+ 2q,

'ik J

a*'2

n

(26)

Let us denote the expression in square brackets in (26) as E() (5) and assume aik = 2a then

4^n G] (g)

-exp

8a

(27)

The expression (27) is the energy E. (5) for the 5-orbital composed of an exponent with the index a [9]. The a. value is determined in [9] analogous to the Madelung constant aM. The

constant aj for the NaCl crystal was calculated in [9] for different a values with the accuracy of up to 25 digits after the decimal point. For example,

3

c g

a = 0.1,

aj = 1.7429785198333593881232629,

a = 1.

aNa = 1.7475645946331821906362119.

a = 10 and a = 100, aNra = 1.7475645946331821906362120.

The Madelung constant aM for NaCl calculated by direct summing in [12] with the accuracy of up to 25 digits after the decimal point is

aM = 1.7475645946331821906362120.

It is seen that one and the same number (with given accuracy) is in the expression (26) for the arbitrary 5-orbital in square brackets if aik > 4 (a> 2) by the order of magnitude. In the case of NaCl, for such

aik the energy E(a(5) = 0.32851544 a.u. given with the accuracy of up to 8 digits after the decimal

point coincides with the Madelung energy. The explanation of this fact is obvious. According to the Gauss theorem, if a spherically symmetric charge distribution does not overlap with a point charge, this charge distribution can be considered as the point charge as well.

Below we consider p-orbitals. According to (20), (23) and (25), the energy of the interaction between the px -orbital E. (px) with the crystal lattice can be written as

1

E. (px )=Z aak

8 i,k

V aik J

4n

Z

2aik

1

VG.(g) (

^exp

2

4ak

1

4 ( ak2

+3 q.l

3 J V n

(28)

The expression in square brackets in (28) is denoted as E(1) (px) and we assume ak = 2a then

e!"(p,Z

4n

4a

--1

V G. (g) (

^exp

§1 v

8a

(29)

The expression (29) is the energy E. (px) for the px -orbital composed of an exponent with the index a. Structural factors G. (g) for the NaCl and KMgF3 crystals are given in [9]. The calculations according to the formula (29) for the Na ion in the case of the NaCl crystal for the values a = 1,2,5 give the value EiN}a(px ) = 0.32851544 a.u., i.e., the Madelung energy. The calculations according to the formula (29) for the Mg ion in the case of the KMgF3 crystal for the values a = 1, 2, 5 give the

value E_M)(px ) = 0.82429794 a.u., i.e., the Madelung energy [9]. The explanation of this fact is also obvious. The crystals are cubic and the quadrupole interaction for the p orbitals is zero. To calculate the energies of the py, pz -orbitals, it is necessary to substitute the explicit variable gx in (28) for the variables gy, gz, respectively.

For example, the interaction between the crystal lattice and the |dxy (r ) and determined following Fj [n1n2n3) functions:

d

2 2 x - y

(r) orbitals is

F} (400) = -nZ

a a,,

yc i,k

V aik J

Z

24ak 2a

- +

1 V g. (g) ( g- V

-exp

4a

(30)

ik J

c g

3

7

Long-range Coulomb interaction in ion crystals The Fj (040) function is obtained by substituting the explicit variable gx in (30) for the variable gy

jr 2 ____

Fj (220) = ---Xa,a„

2v,

c i ,k

X

.if £. -i] m exp f ^

2aik 2aik 22 4ak

(31)

In conclusion of this section, we present the expression for the Fj [n1n2n3) functions in which numerical integration is used. It can be useful when the ak values are rather large

Fj (n!n2n3 ) = - — ni !n2 !n3 !X qp X aibk

Vaik

X-

X f (s , gs )exp

V 4aiku

-f (ni,0)f (n2,0)f (n3,0)

(32)

where Xj. = Xj, Xj2 = yj, x;3 = zj are the coordinates of the vector rj. and analogously xps are the

coordinates of the vector rp. The second term in the formula (32) eliminates the divergence at g = 0 .

If the expression (20) is called the “exact value”, then the formula (32) gives the values which in the most cases coincide with this “exact value” to 8-10 digits after the decimal point already in the calculations with the accuracy set in the programs by default.

3

3

5

3. Discussion

The above TME and OME of the LR Coulomb interaction between the electron and the infinite crystal lattice arise naturally within the approach of [1-4]. The comments about TME follow. For ion crystals or admixture centers in ion crystals, TME is small when the region of the overlap of anion and cation orbitals coincides with the region of the change of the electrostatic potential sign. However, if the average values (r) of the orbitals of ions considerably differ, the value of these matrix elements may

be noticeable and they have to be taken into account in the ab initio calculations. In our opinion, the two-center matrix elements should be also estimated in the calculation of the amplitudes of the charge transfer over the anion or cation sublattice, since the overlap region fits the region of the constant electrostatic potential sign. For example, the amplitudes of the transition over the oxygen sublattice of La2CuO4 were calculated in [13]. The value of these amplitudes is within 0.4 + 0.6 eV but only the

short-range interaction in taken into account. At the same time, for a -NaV2O5 according to (19) and (22) the value of the corrections from the long-range Coulomb interaction over the <j -bond ~ - 0.105 eV, and that over the - -bond ~ 0.174 eV. Thus, in the general case the long-range Coulomb interaction should be estimated when such amplitudes are calculated.

Conclusions

In the general case, the U-region around the admixture ion or the unit cell determined above is several coordination spheres, where the short-range interactions can be written exactly. Formulas obtained in this work allow one to rather exactly take into account the long-range Coulomb interaction as well. Hence, we obtained a method for the ab initio calculation of physical quantities in the U-region, at least at rather small overlap integrals needing no introduction of any parameters.

References

1. Anikeenok O.A. Phys. Solid State 45, 854 (2003)

2. Anikeenok O.A. Phys. Solid State 47, 1100 (2005)

3. Anikeenok O.A. Phys. Solid State 48, 1878 (2006)

4. Falin M.L., Anikeenok O.A., Latypov V.A., Khaidukov N.M., Callens F., Vrielinck H., Hoefstaetter A. Phys. Rev. B 80, 174110 (2009)

5. Ewald P.P. Ann. derPhysik64, 253 (1921)

6. Evjen H.M. Phys. Rev. 39, 675 (1932)

7. Sabry A., Ayadi M., Choukn A. Computational Materials Science 18, 345 (2000)

8. Horiuchi S., Shirakawa T., Ohta Y. Phys. Rev. B 77, 155120 (2008)

9. Anikeenok O.A. Magn. Reson. Solids 11, 1 (2009)

10. Prudnikov A.P., Brychkov Yu.A., Marichev O.I. Integraly i Ryady, Nauka, Moskow (1981) (in

Russian)

11. Ziman J.M. Principles of the Theory of Solids, Cambridge University Press, Cambridge (1972)

12. Hajj F.Y. J. Chem. Phys. 56, 891 (1972)

13. McMahan A.K., Annet J.F., Martin R.M. Phys. Rev. B 42, 6268 (1990)