Transcript MOBILITY Matlab code解說

Student：光心君
Date：2010/04/15
Proof:

2
d V1
dx
2

dx
2

dx
2
2

x
2


 [ A]
N N
[V ]
N 1
 [ R ho ]
N 1
si
V3
V1
 1

V2
 si
2
V 2  2V 3

2
V 1  V 3  2V 2
x
2
d V3
dx
V 2  2 V1
x
2
d V2

2
d V


2
推廣
2
d V
 si
dx
 3
n 
2
V
n 1
V
n 1
2
x
 2V
n  2
V
V
n  n 1
2
2
x
x
 si
經整理變矩陣程式
 2
 x2

 1
 x2

 0

1
x
2
2
x
2
1
x
A
2

  1
0 
 
V


 1

1     2
V2 
2
x     
 V  
2  3 
 3


2
 
 x 
V








推廣
 2
 x2

 1
 x2

 0




0

1
x
2
0
2
1
x
x
2
1
x
2
2
0
2
1
x
x
2
2
0
0
1
x
2

0 

0 




0 

2 
2
 x 
 V1

V
 2



V
 N







 N 1
NN
Rho
 V0=A\Rho initial guess(V0=A-1R)

Code:
A = zeros(N,N); % Matrix for 2nd differential operator
A(1,1)=1/dx0^2; %boundary condition Vsurface=Vs
A(N,N)=1/dx(N-1)^2;
Explain:
Initial Condition  V(1)=Vs ,V(N)=0
1.V(1)=Vs
2.V(N)=0

A (1,1) 
1
2
dx 0
 AV  V1  Vs
and
Rho (1) 
Vs
2
dx 0
1
A( N , N ) 
dx
2
( N 1 )
 AV  V N  0
and
Rho ( N )  0






for j=2:N-1
avgdx=(dx(j-1)+dx(j))/2;  avgdx=Δx
A(j,j-1) = 1/dx(j-1)/avgdx; A(j,j-1) = 1/ Δx2
A(j,j) = -(1/avgdx)*(1/dx(j-1)+1/dx(j)); A(j,j)=-2/ Δx2
A(j,j+1) = 1/dx(j)/avgdx;  A(j,j+1) = 1/ Δx2
end;
1
 2

 x2

 1
J=2  2
x

 0

J=N-1 


0

x
2
0
2
1
x
x
2
1
x
2
2
0
2
1
x
x
2
2
0
0
1
x
2
0 

0 




0 

2 
2
 x 
 V1

V
 2



V
 N
NN







 N 1
%****************CALCULATED PARAMETERS****************

A(1,1)=1/dx0^2; %boundary condition Vsurface=Vs
A(N,N)=1/dx(N-1)^2;
%**************POISSON EQUATION SETUP*****************************

Rho(1)=Vs/dx0^2; %bondary condition on the surface
Rho(N)=0; %bondary condition V(N)=0
 1
 x2

 1
 x2

 0




0

0
0
2
1
x
x
2
1
x
2
2
0
2
1
x
x
2
0
0
0
2

0 

0 




0 

1 
2
 x 
 V1

V
 2



V
 N
NN







 N 1
Vs

R
ho

1
2

dx 0

R ho 2

 

 R ho N 1

0










From subprogram: shhole01R.m
%**************************scale set up
*****************************
xscale = linspace(xstart,xend,N).'; % New scale cm
dx0=(xend-xstart)/real(N);
dx= dx0/au;
% Mesh separation in a.u.
au = 0.5262E-8; % atomic unit in cm(波耳氫原子半徑)

dd=1/2/(dx^2);
% (a.u.)^-2
x ' 

x end  x start
real ( N )
Xscale xstart~xend中有N個元素的行向量
x 
x '
a .u
dd 
1
2x
2

找由1至N所對應的potential值以帶入薛丁格方程式
%**************potential set up*****************************
V=zeros(N,1);
% Potential in Hr (Hr = 27.212; % 1 Hartree in eV)
for j=2:(N-1)
V(j) = interp1(xscaleI,VI,xscale(j))/Hr;
end
V(1)=20;
%boundary condition
V(N)=20;
%boundary condition


Interp1做一維的內插法
Hr Hartree energy ,the atomic unit of energy.
Time independent equation:

2

d  ( x)
2
2m
dx
2
 V ( x) ( x)  E  ( x)
2

d 1
2
dx
2
d 2

2
dx
2
d 3

2
dx
2

其中H=

 2  2 1
x
2
 1   3  2 2
x
2
 2  2 3
x
2
2m
d  ( x)
2
dx
2
V
  n  En n
,且  0
2
1  
(
2m
2
2  
(
2
3  
(
2m
2
x
2
 2  2 3
x
E    N 1
)  V1
 1   3  2 2
2m
   N  N    N 1 
 0,  N 1  0
 2  2 1
x
pf
2
)  V3
推廣
)  V2
2
N  
(
2m
 N 1  2  N
x
2
V=∞
Ψ0=0
V=∞
ΨN+1=0
xstart
xend
)  VN

矩陣形式
 2
  x 2  V1

2
1


2
2m  x


0

推廣
 A1

B

0


0

1
x
2
x
2
 V2
2
1
x
2
B
0
A2
B
0
B
A3
B
0
0
0

0


1

2

x

2
 V3 
2

x
2

 V1

2
mx

2


2

2mx


0


1 



 2
  3 
0 
 1 



0
2









0 


A N  N  N   N  N 1
2
Aj 
B 
mx

2
2
2mx
2
Vj

2
2mx
2
2
mx

2
 V2
2
2mx
2



2


2

2mx

2

 V3 
2
mx

0
1 



 2
  3 
%******** Schrodinger Equation ***************
H = zeros(N,N);% light hole
for j=2:(N-1)
1
A

2

Vj
j
2
H(j,j) = V(j)+2*dd/m1;
2  x  m1
end
H(N,N)=V(N)+2*dd/m1;
H(1,1)=V(1)+2*dd/m1;
for j=2:N
H(j-1,j) = -dd/m1;
1
H(j,j-1) = -dd/m1;  B 
2
2  x  m1
end

解波函數Ψ 與E

[Y,D]=eig(H); % Eigen vectors(Y) and Eigen values(D)
[lambda1,key1] =sort(diag(D));
%sort:以行為單位,將每一行的向量由小到大排列
Y1 = Y(:,key1); 取key1行的一整列的元素
E1=lambda1*Hr+Ev1;
  N N


   N 1 
E    N 1
Y=eigenvector : [Ψ]N ×1(Ψ以行向量的方式儲存在矩陣Y裡)
D=eigenvalue :  E 1 0
0  (Ei存放在矩陣D的對角元素)

0

 0


E N 
0
0

藉由Ψ去計算n(x)
%**************** Calculating hole densities ******************
for j=1:N  hole density /cm2
p1(j)=Do1*k*T*log(1+exp((Ef-E1(j))/k/T)); %hole den in heavy
p2(j)=Do2*k*T*log(1+exp((Ef-E2(j))/k/T)); %hole den in light
end p ( x )   D ( E ) f ( E ) dE


2D
E  E0
2D
Do1 = md1*m0/3.1415/(hb)^2/6.24146E11; %density of state (#/eV/cm2)
D
Do2 = md2*m0/3.1415/(hb)^2/6.24146E11; %density of state (#/eV/cm2)
2D
(E ) 
m

*
2
for j=1:N jj=1:Nhole density at each valley /cm3
YY1(j,jj)=(Y1(j,jj))^2*p1(jj)/dx0; %hole den. (heavy band in #/cc)
YY2(j,jj)=(Y2(j,jj))^2*p2(jj)/dx0; %hole den. (light band in #/cc)

p3 D ( x ) 


E  E0
i
2
D 2 D ( E ) f ( E ) dE


else
求R: hole density /cm3
R(j) = R(j)+(YY1(j,jj)+YY2(j,jj));
if xscale(j)>=xstart & xscale(j)<=xend;
Nep(j) = interp1(xscaleO,R0,xscale(j));
xstart=0.0;
xend=fregion;
p3 D ( x ) 


E  E0
i
2
D 2 D ( E ) f ( E ) dE
Nep(j)=+Nep0*exp(-beta*V0(j));
p  N c  exp[ 
(Ec  E f )
]
k BT
Nen=+ni^2./Nep(只考慮classical case)