说明：
1.代入计算时要特别注意量纲
2.由于代码的编写时间不同，导致同一物理量在不同代码文件中所用的符号不一致
3.如果lsqcurvefit函数不能运行，可能是没有安装相关工具包，可以根据提示安装
4.拟合计算的结果特别依赖初值，以及边界条件，甚至有“人为调参”的痕迹
5.如果有问题欢迎讨论，欢迎指出错误
1.俘获态模型(trapping state model)
俘获态模型假设GB主导载流子的传输，忽视了晶粒电阻的贡献(grain resistivity)，在实际的材料中晶界散射一般是在室温或低于室温下才占主导，可以通过电导率随温度的变化判断一般为电导率与温度成正比例关系。在晶界附近有很多不完全的成键，对应的产生能够俘获电子的缺陷态，从而产生势垒(蓝色)，周围产生空间电荷区。
1.1 当晶界散射占主导时
\sigma_{\mathrm{GB}}=e^2dn\left(\frac{1}{2{\pi}m^*k_{\mathrm{B}}T}\right)^{1/2}\exp\left(\frac{-E_{\mathrm{b}}}{k_{\mathrm{B}}T}\right)
d: 晶粒的尺寸可以理解为晶粒的平均尺寸
E_{\mathrm{b}}: 势垒的高度
E_{\mathrm{b}}=\left\{
\begin{aligned}
&\frac{e^2d^2N}{8\varepsilon} &dN<Q_{\mathrm{t}}\\
&\frac{e^2Q_{\mathrm{t}}^2}{8N\varepsilon} & dN>Q_{\mathrm{t}}\\
\end{aligned}
\right.
Q_{\mathrm{t}}: 俘获态在晶界处的面密度 (C \cdot m^{-2})
\varepsilon: 静态介电常数(F \cdot m^{-1})
N: 离子杂质的浓度(通常为掺杂原子的浓度 (m^{-3})
对于重掺杂的热电材料，只有一部分晶粒被耗尽(depleted)，此时 E_{\mathrm{b}}可以通过 dN>Q_{\mathrm{t}} 的公式来表征
1.2 化简
实验上当一种材料是GBS占主导时电导率的表达式可以简化为下式:
\ln({\sigma}T^{1/2})=\frac{-E_{\mathrm{b}}}{k_{\mathrm{B}}T}+{\mathrm{constant}}
可以通过做\ln({\sigma}T^{1/2})-1/k_{\mathrm{B}}T的关系曲线，这时载流子浓度应该与温度无关；或者做\ln({\mu}T^{1/2})-1/k_{\mathrm{B}}T的关系曲线，斜率就是-E_{\mathrm{b}}
1.3 适用条件
在室温以及室温以下一般是以晶界散射为主，而当温度升高时声学声子散射会变为主导，这一规律可以由Matthiessen’s rule来描述，\mu^{-1}=\mu_{\mathrm{APS}}^{-1}+\mu_{\mathrm{GBS}}^{-1}，由晶界散射主导的载流传输特性是载流子迁移率随温度升高而升高，而当声学声子散射占主导时载流子迁移率随温度升高而降低，另外这一rule要求材料是均匀的(晶粒尺寸差不多)，但是并不是所有多晶材料都满足这一条件。
2.两相模型(Two-phase model)
基于俘获态模型的适用条件，当耗尽层的厚度大于晶界自身或者可以与载流子的平均自由程相比较甚至大于载流子的平均自由程，Matthiessen’s rule被打破，kuo et al提出了两相模型，耗尽层作为分离的相与晶粒相连。
2.1 两相模型下的电导率
\sigma_{\mathrm{total}}^{-1}=(1-t_{\mathrm{GB}})\sigma_{\mathrm{G}}^{-1}+t_{\mathrm{GB}}\sigma_{\mathrm{GB}}^{-1}
S=\frac{S_{\mathrm{G}}\frac{1-t_{\mathrm{GB}}}{K_{\mathrm{G}}}+S_{\mathrm{GB}}\frac{t_{\mathrm{GB}}}{K_{\mathrm{GB}}}}
{\frac{1-t_{\mathrm{GB}}}{K_{\mathrm{G}}}+\frac{t_{\mathrm{GB}}}{K_{\mathrm{GB}}}}
K_{\mathrm{G}}：为温度梯度，由于t_{\mathrm{GB}}一般很小，所以对塞贝克系数影响不大(即晶粒内部的塞贝克系数)。
G指的是晶粒，GB指的是晶界，t_{\mathrm{GB}}是晶界的的尺寸分数(晶界的空间占比)，对于重掺杂半导体，在晶粒内部是以声学声子散射占主导，在晶界内是以晶界散射为主。
\sigma_{\mathrm{AP}}=AT^{-3/2}
\sigma_{\mathrm{GB}}=CT^{-1/2}\exp\left(\frac{-E_{\mathrm{b}}}{k_{\mathrm{B}}T}\right)
\sigma_{\mathrm{total}}^{-1}=\left(A_{1}T^{-3/2}\right)^{-1}+\left[C_{1}T^{-1/2}\exp\left(\frac{-E_{\mathrm{b}}}{k_{\mathrm{B}}T}\right)\right]^{-1}
2.2 两相模型的典型例子
2.3 计算代码
1.需要将文件与代码放在同一目录下
2.文献中的数据通过 origin 获取，与实际的数据会略有差别
(Hf0.3Zr0.7)0.98Nb0.02CoSb.txt
温度(K) 电导率
K sigma(104·S·m-1)
Temperature Electrical conductivity
310.771917785248 0.33558841727966
364.533472301409 0.467121501203444
411.67822010789 0.624275209409123
463.372022527276 0.778920169255976
512.584522430531 0.93311064964373
559.729270237011 1.04097231774284
608.941770140267 1.11498378384588
662.703324656428 1.17916332791581
712.329374979039 1.17916332791581
761.541874882294 1.21262542731353
811.167925204905 1.21262542731353
860.380425108161 1.24386883180347
909.592925011416 1.22814778169716
959.218975334027 1.2629999537189
1010.49922733406 1.22814778169716
1057.64397514054 1.2629999537189
1102.72097085024 1.27916713315169
(Hf0.3Zr0.7)0.98Nb0.02CoSb.m
%基于两相模型拟合电导率
clc,clear
%(Hf0.3Zr0.7)0.98Nb0.02CoSb
%原始数据
%第1列: 温度
%第2列: 电导率
fileID = fopen('(Hf0.3Zr0.7)0.98Nb0.02CoSb.txt'); %获取打开文件的文件句柄(指针)
data = textscan(fileID, '%f%f', 'headerlines',3); %只读取原始数据的第1-2列，从第4行开始读取
fclose(fileID);%关闭文件句柄
%将元胞数组转换为矩阵
matrix = cell2mat(data);
Tdata = matrix(:,1); %温度 K
Sigmadata = matrix(:,2)*10^4; %d电导率 S/m
Kb = 1.38066*10^-23; %J/K 1J=1V*A*s 1W = 1J/s
h = 6.62608*10^-34; %J*s
hr = h/(2*pi); %简约普朗克常数
e = 1.60219*10^-19; %1C=1A*s
me = 9.10953*10^-31; %Kg
K = Kb/e;
param0 = [10^9,10^6, 95]; %参数初值 最终结果和初值有很大关系
%param(1): A1; param(2): C1; param(3):Eb
Sigma = @(param, Tdata)((param(1).*Tdata.^(-3/2)).^(-1)+(param(2).*Tdata.^(-1/2).*exp(-param(3)*10^(-3)*e./(Kb*Tdata))).^(-1)).^(-1); %同时拟合三个参数
[param, resnorm] = lsqcurvefit(Sigma, param0, Tdata, Sigmadata); %emf:拟合的得到的态密度有效质量，resnorm残差平方和
display(param); %显示拟合得到的值
Tdatac = 200:1200; %c: curve 用于绘制计算得曲线
SigmaGBSdatac = param(2).*Tdatac.^(-1/2).*exp(-param(3)*10^(-3)*e./(Kb*Tdatac)); %晶界散射占主导
SigmaAPSdatac = param(1).*Tdatac.^(-3/2); %声学声子散射占主导
SigmaGBSAPSdatac = Sigma(param, Tdatac); %同时考虑晶界散射与声学声子散射
%semilogx(x,y) 绘图，x轴对数刻度，y轴线性刻度，
%semilogy(x,y) 绘图，x轴线性刻度，y轴对数刻度，
%loglog(x,y) 绘图，x,y轴都为对数刻度
figure(1)
%仅考虑GBS
semilogy(Tdatac, SigmaGBSdatac*10^(-4), 'b-', 'LineWidth', 2); %以10为底的对数坐标(y轴)
hold on
%仅考虑APS
semilogy(Tdatac, SigmaAPSdatac*10^(-4), 'm-', 'LineWidth', 2); %以10为底的对数坐标(y轴)
hold on
%考虑GBS+APS
semilogy(Tdatac, SigmaGBSAPSdatac*10^(-4), 'r-', 'LineWidth', 2); %以10为底的对数坐标(y轴)
hold on
semilogy(Tdata, Sigmadata*10^(-4), 'ro'); %数据点为空心圆
hold on
xlabel('T [K]');
ylabel('\sigma [10^{4}\cdotS\cdotm^{-1}]');
xlim([200,1200]);
ylim([0.1,10]);Mg2Si0.9875Sb0.0125+4wt.%MgO.txt
温度(K) 电导率
K sigma(104·S·m-1)
Temperature Electrical conductivity
320.986390367611 9.48126337633109
339.936220013887 10.1594771367606
360.382088842764 10.5529574699815
386.064094810745 11.0249707848602
409.00336227729 11.4387927669319
426.457152740965 11.8000136040463
452.139158708945 12.1027592699827
475.327766039257 12.0332784437531
499.513732824637 12.0332784437531
522.702340154949 11.9504277761214
546.888306940328 11.6514926213687
572.570312908308 11.2948182937323
598.252318876288 10.8862048859465
621.4409262066 10.5529574699815
647.12293217458 10.2181385289009
671.308898959959 9.84847623802204
696.990904927939 9.48126337633109
721.426211577085 9.19102388944388
748.604256727666 8.68679738578482Mg2Si0.9875Sb0.0125+4wt.%MgO.m
%基于两相模型拟合电导率
clc,clear
%Mg2Si0.9875Sb0.0125+4wt.%MgO
%原始数据
%第1列: 温度
%第2列: 电导率
fileID = fopen('Mg2Si0.9875Sb0.0125+4wt.%MgO.txt'); %获取打开文件的文件句柄(指针)
data = textscan(fileID, '%f%f', 'headerlines',3); %只读取原始数据的第1-2列，从第4行开始读取
fclose(fileID);%关闭文件句柄
%将元胞数组转换为矩阵
matrix = cell2mat(data);
Tdata = matrix(:,1); %温度 K
Sigmadata = matrix(:,2)*10^4; %d电导率 S/m
Kb = 1.38066*10^-23; %J/K 1J=1V*A*s 1W = 1J/s
h = 6.62608*10^-34; %J*s
hr = h/(2*pi); %简约普朗克常数
e = 1.60219*10^-19; %1C=1A*s
me = 9.10953*10^-31; %Kg
K = Kb/e;
param0 =[2*10^9,5.5*10^6, 95]; %参数初值 最终结果和初值有很大关系
%param(1): A1; param(2): C1; param(3):Eb
Sigma = @(param, Tdata)((param(1).*Tdata.^(-3/2)).^(-1)+(param(2).*Tdata.^(-1/2).*exp(-param(3)*10^(-3)*e./(Kb*Tdata))).^(-1)).^(-1); %同时拟合三个参数
[param, resnorm] = lsqcurvefit(Sigma, param0, Tdata, Sigmadata); %emf:拟合的得到的态密度有效质量，resnorm残差平方和
display(param); %显示拟合得到的值
Tdatac = 250:850; %c: curve 用于绘制计算得曲线
SigmaGBSdatac = param(2).*Tdatac.^(-1/2).*exp(-param(3)*10^(-3)*e./(Kb*Tdatac)); %晶界散射占主导
SigmaAPSdatac = param(1).*Tdatac.^(-3/2); %声学声子散射占主导
SigmaGBSAPSdatac = Sigma(param, Tdatac); %同时考虑晶界散射与声学声子散射
%semilogx(x,y) 绘图，x轴对数刻度，y轴线性刻度，
%semilogy(x,y) 绘图，x轴线性刻度，y轴对数刻度，
%loglog(x,y) 绘图，x,y轴都为对数刻度
figure(1)
%仅考虑GBS
semilogy(Tdatac, SigmaGBSdatac*10^(-4), 'b-', 'LineWidth', 2); %以10为底的对数坐标(y轴)
hold on
%仅考虑APS
semilogy(Tdatac, SigmaAPSdatac*10^(-4), 'm-', 'LineWidth', 2); %以10为底的对数坐标(y轴)
hold on
%考虑GBS+APS
semilogy(Tdatac, SigmaGBSAPSdatac*10^(-4), 'r-', 'LineWidth', 2); %以10为底的对数坐标(y轴)
hold on
semilogy(Tdata, Sigmadata*10^(-4), 'ro'); %数据点为空心圆
hold on
xlabel('T [K]');
ylabel('\sigma [10^{4}\cdotS\cdotm^{-1}]');
xlim([250,850]);
ylim([5,40]);2.4 计算结果
HZNCS计算结果图
MSS+MO计算结果图
3.基于Drude-Sommerfeld自由电子模型
\mu_{\mathrm{w}}=\frac{3h^3\sigma}{8{\pi}e(2m_{\mathrm{e}}k_{\mathrm{B}}T)^{3/2}}\left[\frac{\exp\left[\frac{|S|}{k_\mathrm{B}/e}-2\right]}{1+\exp\left[-5\left(\frac{|S|}{k_\mathrm{B}/e}-1\right)\right]}+\frac{\frac{3}{\pi^2}\frac{|S|}{k_\mathrm{B}/e}}{1+\exp\left[5\left(\frac{|S|}{k_\mathrm{B}/e}-1\right)\right]}\right]
3.1 加权(权重)载流子迁移率与迁移率的关系(用于拟合得到势垒)
\mu_{\mathrm{w}}\approx\mu\left(\frac{m^{\mathrm{*}}}{m_{\mathrm{e}}}\right)^{3/2}
3.2 加权载流子迁移率计算代码
4.参考文献
[2]: Snyder, G.J., et al., Weighted Mobility. Advanced Materials, 2020. 32(25): p. 2001537.
评论区