(*Программа «Two-temperature Thermal
Conductivity in Solid Metal due to only Electron-Phonon Scattering» -
сокращенное название суть “Electron-Phonon 2T-sigma”*)
(*Вычисляет частоту столкновений s-электрона с импульсом p с фононами и затем интегрированием по p коэффициент теплопроводности*)
(*температура электронов TeK, ионов TiK*)
(*алюминий*)
TeK=800. ;
(*температура электронов в К*)
mu= -0.0000129462;
(*химпотенциал электронов в а.е.*)
TiK=800.
(*температура ионов в К*)
Ti=TiK/11605./27.2 ;
(*температура электронов в а.е.*)
Te=TeK/11605./27.2 ;
(*температура электронов в а.е.*)
rho=2.6989;
(*плотность
в г/см^3*)
A=26.982 ;
(*молярная масса *)
z=3.
(*заряд иона в а.е.*)
M=1836.*A ;
(*масса атома в а.е.*)
na= rho/A*6.022*0.148*0.1;
(*концентрация атомов в а.е.*)
re=M*na ;
(*плотность в атомных единицах*)
qD=(6*Pi*Pi*na)^(1./3)
(*квазиимпульс Дебая в атомных единицах*)
si=6.260 ;
(*скорость продольного звука в км/сек*)
vat=3.*(10^5)/137 ;
(*а.е. единица скорости*)
s1=si/vat
s=s1/5.1
(*скорость звука в атомных единицах*)
ens= -11.1/27.2;
(*ens - дно s-зоны в а.е.*)
mu0=0;
(*положение химпотенциала электронов при Te=0*)
m=((3.*z*na*Pi*Pi)^(2./3))/(2*(mu0-ens))
(*масса s-электрона в а.е.*)
(*вычисление константы томас-фермиевского
экранирования кулоновского взаимодействия за счет s-электронов. Это kap*)
Xx=100
te=(ens-mu)/Te ;
hs[tk_]=
Sqrt[tk]*Exp[te-tk]/((Exp[te]+Exp[-tk])^2);
UK=
NIntegrate[hs[tk],{tk,0.,Infinity}];
dnsdm=
Sqrt[2.]/(Pi*Pi)*((Sqrt[m*Te])^3)*UK/Te ;
kap=Sqrt[4*Pi*dnsdm]
(*обратная длина томас-фермиевского экранирования в а.е.*)
(*вычисление производной химпотенциала
по электронной температуре dmu/dTe*)
mus=mu-ens ;
(*химпотенциал, отсчитанный от дна зоны проводимости*)
ef1[x_]=Exp[-x*x/Te]/(
Exp[-x*x/Te]+Exp[-mus/Te]) ;
is1= NIntegrate[ef1[x],{x,0.,Infinity}] ;
ef2[x_]=(3*x*x-mus)*ef1[x]
;
is2= NIntegrate[ef2[x],{x,0.,Infinity}]
;
mut= -is2/is1/Te
(*производная
химпотенциала по электронной температуре dmu/dTe в а.е.*)
xx=11111
(*вычисление теплоемкости s-электронов на единицу объема cvs*)
fcv[x_]=(3*x*x*(mut-mus/Te)+5*x*x*x*x/Te)*ef[x] ;
icv= NIntegrate[fcv[x],{x,0.,Infinity}];
cvs=Sqrt[2.]/Pi/Pi*((m)^(3/2.))*icv
(*теплоемкость единицы объема s-электронов в а.е.)
(*вычисление среднего квадрата скорости s-электронов vs2*)
ezv=Exp[-xv] ;
hes=Exp[(ens-mu)/Te]
;
hs[xv_]=Sqrt[xv]*ezv/(hes+ezv);
ihs=NIntegrate[hs[xv],{xv,0.,Infinity}];
hvs[xv_]=xv*hs[xv] ;
ihvs=NIntegrate[hvs[xv],{xv,0.,Infinity}];
vs21=2*Te/m*ihvs/ihs
(*средний квадрат скорости s-электронов vs2 в а.е.*)
pF=Sqrt[2*m*(mu0-ens)]
(*Ферми-импульс в а.е.*)
en[p_]=ens+p*p/(2*m) ;
(*закон дисперсии электронов*)
(*вычисление частот столкновений электрона с импульсом p с
фононами*)
cnu=Pi/(re*s*s)*s*((4*Pi*z*na)^2)*2*Pi*m/((2*Pi)^3) ;
(*постоянный множитель для этих частот*)
w[q_]=cnu*q*q/((q*q+kap*kap)^2) ;
(*функция импульса фонона q, содержащая квадрат матричного элемента и множитель от фазового объема фонона*)
(*вычисляет случай 1 с испусканием фонона*)
p1=Sqrt[p*p] ;
(*модуль импульса электрона после испускания фонона. Здесь пренебрегаем изменением абсолютной величины импульса электрона в результате испускания фонона.*)
(*Выражение для частоты столкновений 
электрона с импульсом 
с фононами при их испускании (индекс i означает irradiation). Это выражение представляет собой интеграл.
Причем пределы интегрирования зависят от значения импульса p электрона. Как
эти пределы зависят от импульса, поясняется в сопроводительном тексте
«Пояснения к программе “Electron-Phonon 2T-sigma & kappa”».
Кроме того пояснения приведены ниже при описании рис. 1. В
программе подынтегральная функция в выражении для записывается в
следующем виде:*)
wp[p_,q_]=w[q]/p*(q*q-(p-p1)*(p-p1))/(2*p*p1)*Exp[s*q/Ti] /(Exp[s*q/Ti]-1) ;
(*Комбинация из нескольких постоянных сомножителей, входящая
в интеграл , была обозначена выше как «cnu» (constanta nu)*)
Рис. 1. На рисунке изображена область
интегрирования в интеграле в плоскости 
. Это плоскость электронного импульса 
и импульса фонона q. Метка 
на оси q на
рисунке соответствует дебаевскому импульсу. Метка 
– фермиевский импульс, где z –
число электронов в зоне проводимости в расчете на один атом. Прямая
соответствует точному
рассеянию электрона назад. Квазиимпульс фононов ограничен значением 
. Поэтому область интегрирования ограничивают две прямые: и 
.
(* 0<p<0.5*qD *)
nup1[p_]=
NIntegrate[wp[p,q],{q,0.0001, 2*p }] ;
(*В этой строке находится частота
столкновений при испускании фононов электрона с импульсом 0<p<0.5*qD. См. рис. 1.
Из-за излома на границе области интегрирования мы разбиваем интеграл на два – один с 0<p<0.5*qD и другой с 0.5*qD <p<Infinity *)
(* 0.5*qD <p<Infinity *)
nup3[p_]= NIntegrate[wp[p,q],{q, 0.0001, qD}] ;
(*частота столкновений при испускании фононов электрона с импульсом 0.5*qD <p<Infinity *)
xx=123
(*вычисляет случай 2 с поглощением фонона*)
pm=Sqrt[p*p] ;
(*импульс электрона после поглощения фонона*)
wm[p_,q_]=w[q]/p*(q*q-(p-pm)*(p-pm))/(2*p*pm)/(Exp[s*q/Ti]-1) ;
(*подынтегральная функция для вычисления частот столкновения при поглощении фонона. Эта строка соответствует интегралу *)
(* 0<p<0.5*qD *)
num2[p_]=
NIntegrate[wm[p,q],{q, 0.0001, 2*p }] ;
(*частота столкновений при поглощении фононов электрона с импульсом 0<p<0.5*qD *)
(* 0.5*qD <p<Infinity *)
num4[p_]= NIntegrate[wm[p,q],{q, 0.0001, qD}] ;
(*частота столкновений при поглощении фононов электрона с импульсом 0.5*qD<p<Infinity *)
xx=12345
nuF=nup3[pF]+num4[pF]
;
(*частота столкновений с фононами электрона с импульсом Ферми в а.е.*)
nuFa=nuF*41
(*частота столкновений с фононами электрона с импульсом Ферми в единицах 10^15 1/сек*)
(*вычисление коэффициента теплопроводности*)
faw[p_]
= 1./(( Exp[(en[p]-mu)/(2*Te)]+ Exp[-(en[p]-mu)/(2*Te)])^2)*
((en[p]-mu)/Te+mut)/Te*(en[p]-mu)*
p*p/m/m*2*Pi*p*p*2/((2*Pi)^3)*(2/3.);
(*faw - часть подынтегральной функции для вычисления коэффициента электропроводности. Коэффициент электропроводности представляет собой интеграл
.
faw равно подынтегральному выражению без
множителя с частотой 
, где 
*)
iaw2=NIntegrate[faw[p]/(nup1[p]+num2[p]),{p, 0, 0.5*qD}]
(*полный вклад в теплопроводность электронов с импульсом 0<p<0.5*qD *)
iaw4=NIntegrate[faw[p]/(nup3[p]+num4[p]),{p, 0.5*qD,
50. }]
(*полный вклад в теплопроводность электронов с импульсом 0.5*qD <p<Infinity*)
kappa= iaw2+ iaw4
(*коэффициент теплопроводности в а.е.*)
kappaSI=kappa1*1.090523*(10^4)
(*это – коэффициент теплопроводности в Вт/(м сек К) (система СИ) *)
nuat=(1/3.)*cvs*vs2/iaw ;
(*это – частота s-i столкновений в атомных единицах*)
nus=nuat*41
(*это – частота в единицах 10^15 1/сек*)