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