Electron correlation in 1D electron gas Yan Jun July 26, 2008

1

The equation of motion [1]

The equation of motion for the classical one-particle distribution function f (x, p, t) under an external potential Vext (x, t) can be described using Boltzmann’s equation [2], ∂f (x, p, t)/∂t + p/m · ∇x f (x, p, t) − ∇x Vext (x, t) · ∇p f (x, p, t) Z − ∇x υ(x − x0 ) · ∇p f (x, p, x0 , p0 , t)dx0 dp0 = 0,

(1)

where υ(x) is the Coulomb interaction potential and f (x, p, x0 , p0 , t) is the two-particle distribution function. We define, f (x, p, x0 , p0 , t) = f (x, p, t)f (x0 , p0 , t)g(x − x0 ),

(2)

where g(x) is the static pair-correlation function. It is a finite value at the origin and one at long distances. In Hartree approximation, where the system is uncorrelated, g(x) = 1 everywhere. In Hartree-Fock approximation, where the exchange interaction is considered, g(x) interpolate from to 1 at large distance [3, 4]. 1

1 2

at the origin

If one writes f (x, p, t) = f0 (p) + δf (x, p, t),

(3)

where δf (x, p, t) denotes the deviation from the equilibrium distribution function f0 (p) induced by the external potential. By Linearizing Eq. (1), one has the following equation: [∂/∂t + p/m · ∇x ] δf (x, p, t)
R − ∇x Vext (x, t) + ∇x ψ(x − x0 )δf (x0 , p0 , t)dx0 dp0 · ∇p f0 (p) = 0,

(4)

where ∇x ψ(x) = g(x)∇x υ(x).

(5)

The effective field felt by the particle is then, Eef f (x, t) = −∇x Vef f (x, t) Z

= −∇x Vext (x, t) − ∇x υ(x − x0 )δf (x0 , p0 , t)dx0 dp0 Z − [g(x − x0 ) − 1]∇x υ(x − x0 )δf (x0 , p0 , t)dx0 dp0 .

(6)

The first two terms on the right-hand side of Eef f (x, t) are the usual macroscopic electric field, and the third term corresponds to the local-field corrections. Only the first two terms are present in the RPA.

2

Deduction of G(q)

Start from Eq. (5), we will obtain ψ(q),υ(q), and G(q) using Fourier transform. Eq. (5) can be written as, ∇x ψ(x) = ∇x υ(x) + [g(x) − 1]∇x υ(x). 2

The Fourier transform of ψ(x) and υ(x) is given by, X

ψ(x) =

ψ(q)eiq·x ,

q

υ(x) =

X

0

υ(q0 )eiq ·x .

q0

The pair correlation function g(x) is related to the static structure factor through [3] g(x) − 1 =

1 X [S(k) − 1] eik·x N k

(7)

Then Eq. (5) becomes, X

X 1 X 00 [S(k) − 1] eik·x iq00 υ(q00 )eiq ·x N k q0 q00 X 1 XX 00 0 [S(k) − 1] iq00 υ(q00 )ei(k+q )·x . = iq0 υ(q0 )eiq ·x + N k q00 q0

iqψ(q)eiq·x =

q

X

0

iq0 υ(q0 )eiq ·x +

Under the condition, q = q0 = k + q00 , the above equation must satisfy that, qψ(q) = qυ(q) +

1 X [S(q − q00 ) − 1] q00 υ(q00 ). N q00

ψ(q) = υ(q) +

1 X q0 υ(q0 ) [S(q − q0 ) − 1] N q0 q

As a result,

= υ(q)[1 − G(q)],

(8)

where G(q) = −

1 X q0 υ(q0 ) [S(q − q0 ) − 1] N q0 q υ(q) 3

(9)

G(q) is different for 3D, 2D and 1D case:  Z q · q0 1 d3 q0  0  − [S(q − q ) − 1] ,    n q02 (2π)3        Z  q · q0 d2 q0 1 G(q) = [S(q − q0 ) − 1] , −  n q · q0 (2π)2          1 Z q 0 υ(q 0 ) dq 0   [S(q − q 0 ) − 1] ,  − n qυ(q) 2π

(3D) [1]

(2D) [5]

(10)

(1D) [6]

where n = N/Ω is the electron density. We have used in the above equation, υ(q) = 4πe2 /q 2 and υ(q) = 2πe2 /q for the case of 3D and 2D, respectively.

3

Deduction of Vef f (q, ω)

Starting from Eq. (6), we will obtain Vef f (q, ω) by Fourier transform. Eq. (6) is written as, Z ∇x Vef f (x, t) = ∇x Vext (x, t) +

∇x ψ(x − x0 )δf (x0 , p0 , t)dx0 dp0 .

The Fourier transform of each quantity in the above equation is (discard the coefficient denoting the volume), Vef f (x, t) =

X

Vef f (q, ω)eiq·x e−iωt ,

q,ω

Vext (x, t) =

X

Vext (q, ω)eiq·x e−iωt ,

q,ω

ψ(x − x0 ) =

X

0

ψ(q)eiq·(x−x ) ,

q,ω

δf (x0 , p0 , t) =

X

0

δf (q, ω, p0 )eiq·x e−iωt .

q,ω

4

Finally, Vef f (q, ω) can be expressed as, Vef f (q, ω) = Vext (q, ω) + ψ(q)δf (q, ω, p0 )Ωdp0 = Vext (q, ω) + υ(q)[1 − G(q)]δn(q, ω),

(11)

where δn(q, ω) is the induced charge density.

4

Deduction of (q, ω), χ(q, ω)

The induced charge density can be expressed as, δn(q, ω) = χ0 (q, ω)Vef f (q, ω)

(12)

= χ(q, ω)Vext (q, ω).

(13)

Combine Eq. (11)-(13), we can obtain, χ(q, ω) =

χ0 (q, ω) , 1 − υ(q)[1 − G(q)]χ0 (q, ω)

(14)

where G(q) = 0 reduces to the case of RPA. The dielectric function is defined as, Vext (q, ω) VH (q, ω) + Vext (q, ω) 1 = 1 + υ(q)χ(q, ω) υ(q)χ0 (q, ω) = 1− 1 + υ(q)χ0 (q, ω)G(q)

(q, ω) =

where we have used VH (q, ω) = υ(q)δn(q, ω) [5]

5

(15) (16) (17)

5

The self-consistent solution

Equations (9), (14) and (17), together with the relation between the structure factor and the imaginary part of the dielectric function, provide a set of equations which have to be solved self-consistently. They are written in the following for the 1D’s case [6], χ(q, ω) =

χ0 (q, ω) , 1 − υ(q)[1 − G(q)]χ0 (q, ω) υ(q)χ0 (q, ω) , 1 + υ(q)χ0 (q, ω)G(q)

(19)

q 0 υ(q 0 ) dq 0 [S(q − q 0 ) − 1] , qυ(q) 2π

(20)

(q, ω) = 1 − 1 G(q) = − n

Z

(18)

~ S(q) = − πn

Z

∞

Im[χ(q, ω)]dω   Z ∞ 1 ~ Im = − dω. πnυ(q) 0 (q, ω)

(21)

0

We have obtained in the last report that, (in atomic units)  2  2 ω − ω− (q) 1 Reχ0 (q, ω) = ln 2 , 2 πq ω − ω+ (q)   −1/q (ω (q) < ω < ω (q)) − + Imχ0 (q, ω) =  0 Otherwise,

(22)

(23) (24)

where ω± (q) = |kF q ± 21 q 2 |.

References [1] K. S. Singwi, M. P. Tosi, R. H. Land, and A. Sj¨olander, Electron correlations at metallic densities, Phys. Rev. 176, 589 (1968). 6

[2] http : //en.wikipedia.org/wiki/Boltzmann equation [3] David Pines, Elementary excitations in solids, p74-76, W. A. Benjamin, Inc. 1964. [4] Jorge Kohanoff, Electronic structure calculations for solids and molecules - Theory and computational methods, p26, Cambridge, 2006. [5] T. Inaoka, T. Nagao, S. Hasegawa, T. Hildebrandt, and M. Henzler, Twodimensional plasmon in a metallic monolayer on a semiconductor surface: Exchange-correlation effects, Phys. Rev. B 66, 245320 (2002). [6] W. I. Friesen and B. Bergersen, Dielectric response of a one-dimensional electron gas, J. Phys. C 13, 6627 (1980).

7

% ----------Common Parameters, v(q) and w+(q), w-(q) ---------kf=0.41*0.529177;a=1.0/0.529177;m0=0.6; n=kf/pi;epsilon=6.25; % single-spin component % Be cautious to choose Ni, Nj, dq, dw % Finally, only 1:Ni/2*dq is valid, due to the calculation of G(q) Ni=400;Nj=2000; %w(Nj) should larger than w1(Ni) dq=0.001;dw=0.0003; % atomic unit 0.005/a.u. 0.0003Ry = 0.0082eV % Calculate w+(q) and w-(q) in atomic unit for i=1:Ni q(i)=i*dq; w1(i)=(kf*q(i)+0.5*q(i)^2)/m0;

% w_plus

w2(i)=abs(kf*q(i)-0.5*q(i)^2)/m0; % w_minus end %figure;plot(q,w1,'ok',q,w2','or');hold on; % error report if(w1(Ni) > Nj*dw) warning('w is smaller than w+, … please reset Nj or Ni');pause; end % Calculate Coulomb potential v(q), Bergesen's paper, % weekly report 2008-8-1 for i=1:Ni z=(a*q(i))^2; v(i)=quadl(@(x)exp(-x)./x,z,50,1e-6); v(i)=v(i)*exp(z)/epsilon; end %figure;plot(q,v,'ok');hold on;pause % ------------------------- STLS ------------------------% Calculate chi0(q,w), weekly report 2008-7-26, also Bergesen's paper chi0(1:Ni,1:Nj)=0.+0.i; for i=1:Ni for j=1:Nj w(j)=j*dw; temp1=m0/(pi*q(i))*log(abs((w(j)*w(j)-w2(i)*w2(i)+1.d-50)/… (w(j)*w(j)-w1(i)*w1(i)+1.d-50))); % Real part if(w2(i)
% Self-consistent iternaration. G(1:Ni)=0.;chi(1:Ni,1:Nj)=chi0(:,:);% Initialize l=['-k';'-r';'-y';'-r';'-y';'-r';'-y';'-r';'-y';'-r']; for iter=1:15 chiold(1:Ni,1:Nj)=chi(1:Ni,1:Nj); % Update chi for i=1:Ni for j=1:Nj chi(i,j)=chi0(i,j)/(1.-v(i)*(1.-G(i))*chi0(i,j)); end end % Whether converged? conv=0; for i=1:Ni/2 for j=1:Ni/2 conv=conv+abs(chi(i,j)-chiold(i,j)); end end conv if(conv <1d-3) break;end % Calculate S(q) for i=1:Ni S(i)=-sum(imag(chi(i,:)))*dw/(pi*n); end %figure;plot(q,S(:)); % Calculate G(q) G(1:Ni)=0.; for i=1:Ni for j=-Ni:Ni % !!! Especially note here if (j~=0 && abs(i-j)0) if(j>0) G(i)=G(i)-q(abs(j))*v(abs(j))/… (q(i)*v(i))*(S(abs(i-j))-1)*dq/(2.*pi*n); elseif(j<0) G(i)=G(i)+q(abs(j))*v(abs(j))/… (q(i)*v(i))*(S(abs(i-j))-1)*dq/(2.*pi*n); end end end

end %plot(q(1:Ni/2)/kf,G(1:Ni/2),l(iter,:));hold on; end %figure;plot(q/kf,G(:)); % find out plasmon energy, % the maximum of chi(q,w) at each q corresponding to wp, % Note, we consider only half Ni/2*dq, because the calculation % of G(q) for q>Ni/2*dq is invalid. for i=1:Ni/2 temp=0; for j=1:Nj if(temp