硅的能带结构 - a Monte Carlo view

Pop E. Self-heating and scaling of thin body transistors[M]. stanford university, 2005

---

能带是半导体最重要的性质之一，也是Electron Monte Carlo最重要的输入，能带就是电子能量和动量的关系，也就是色散关系.

对于一个真空中的自由电子，他的色散关系很简单，

$$E=\frac{1}{2}m(v_x^2+v_y^2+v_z^2)=\frac{p_x^2+p_y^2+p_z^2}{2m}$$ 根据德布罗意关系，我们知道粒子的动量和波长满足 $$p_x=\hslash k_x,p_y=\hslash k_y,p_z=\hslash k_z$$ 在实际半导体中，色散关系并不是像自由电子这样这么简单的，半导体导带中的自由电子和内部原子、其他电子均存在着相互作用，即使对研究最充分结构最简单的硅来说，它的电子能带结构也是比较复杂的.

但幸运的是，电子满足的费米-狄拉克分布告诉我们，电子的数目在超过费米能级后，其数目会随着能量的增加而大幅衰减. 因此，对于电子导电的半导体来说，自由电子可能允许的状态基本都集中在导带的底部，相当于被困在了这个势阱里. 因此对于大部分情况来说，我们没有必要给出所有的导带形貌，我们只需要对导带底部做出一些合理的描述即可. 由于导带底部是导带的能量最低点，相当于稳定平衡位置，而在平衡位置附近由于势能的一阶导数（力）为0，所以任何势能分布都可以采用二阶简谐来近似，也就是抛物线，弹簧的胡克定律或者固体物理中推导的晶格相互作用，都是基于这种近似得到的.

所以我们可以认为，对于导带底部分布的电子，它的色散关系的形式应该还是满足

$$E=E_c+\frac{1}{2}m(v_x^2+v_y^2+v_z^2)=\frac{p_x^2+p_y^2+p_z^2}{2m}$$ 这种形式，其中$E_c$是导带底部的能量. 由于半导体内部原子等其他因素的作用，此时的分母会偏离电子的实际质量$m$，我们把内部势场的对电子的影响全部弄到分母上去，使得色散关系还能是抛物线的形状，我们把新的分母定义为电子的有效质量，不同半导体的有效质量是不同的. 同时，由于晶体具有各向异性的性质，不同方向的电子有效质量也是不同的，我们于是可以把导带底部的色散关系写成 $$E=E_c+\frac{{\hslash }^2}{2}[\frac{(k_x-k_{0,x}{)}^2}{m_x^{\ast }}+\frac{(k_y-k_{0,y}{)}^2}{m_y^{\ast }}+\frac{(k_z-k_{0,z}{)}^2}{m_z^{\ast }}]$$ 其中$m_x^{\ast },m_y^{\ast },m_z^{\ast }$分别是x、y、z方向的电子有效质量，$k_{0,x},k_{0,y},k_{0,z}$是导带底的坐标，可以看出这是一个旋转椭球的色散关系. 注意我们在这篇文章中讨论的坐标全部是动量空间（波矢空间）的概念，与实空间分布没有任何关系，动量空间的分布只决定电子的性质. 从这里也可以看出我们用Monte Carlo方法是怎么理解波粒二象性的，我们用波的色散描述粒子的性质，驱动粒子在实空间中进行运动.

所以表达式就是这样了，假如我们想针对一个实际半导体，比如硅，我们要描述它的导带底部的色散关系，我们需要知道什么呢？

首先是$E_c$，如果我们把价带顶的能量定义为0，那么$E_c$就是带隙$E_g$了，它是一个温度相关的量. 根据实验，我们可以得到下面这个关系式，

$$E_g(T)=1.1756-8.3131\times {10}^{-5}T-2.6814\times {10}^{-7}{T}^2$$ 在0K时Si的带隙为1.1756 eV. 我们还需要知道导带底部的位置，要知道对于实际半导体而言，它导带能量最小值的点并不在布里渊区的原点$(0,0,0)$！而且由于晶格的对称性，不一定只存在一个能量最低点，而是可能存在很多个对称的能量最低点。对于硅来说，它存在6个对称的导带最小值点，分别出现在$[100],[\overline{1}00],[010],[0\overline{1}0],[001],[00\overline{1}]$这几个方向上，且导带最低点的位置位于对应方向布里渊区中心到布里渊区边界的0.85倍处的位置.

$[100]$方向布里渊区边界的波矢和动量分别是

$$q_m=2\pi /a,p_m=\hslash q_m$$ a是硅的晶格常数. 于是硅的导带底部的能带结构可以表示成了 $$E^s=E_c+\frac{{\hslash }^2}{2}[\frac{(k_x-k_{0,x}^s{)}^2}{m_x^{\ast }}+\frac{(k_y-k_{0,y}^s{)}^2}{m_y^{\ast }}+\frac{(k_z-k_{0,z}^s{)}^2}{m_z^{\ast }}]$$ 其中$s=1,2,3,4,5,6$，分别代表6个对称的导带底. 另一个需要的量就是等效质量，由于旋转椭球面的transverse方向是对称的（相当于声子有两个T声学支），我们可以把等效质量约化成横向有效质量和纵向有效质量两个有效质量， $$m_l=0.9163\times m_0,m_t=0.1880\times m_0\times \frac{E_g(0)}{E_g(T)}$$ 现在我们基本上可以完全描述硅的能带结构了，但是人们发现这个抛物线近似并不是太好.. 因为它毕竟只是简谐近似，只会在电子能量很低、非常接近导带底的时候才近似成立. 那么只要引入一定非抛物近似，就可以扩展近似成立的范围了. 于是我们就有了non-parabolic Kane band， $$E_k^s(1+\alpha E_k^s)=\frac{{\hslash }^2}{2}[\frac{(k_x-k_{0,x}^s{)}^2}{m_x^{\ast }}+\frac{(k_y-k_{0,y}^s{)}^2}{m_y^{\ast }}+\frac{(k_z-k_{0,z}^s{)}^2}{m_z^{\ast }}]$$ 我们把电子高于导带底部的能量定义为$E_k$，再引入一个非抛物系数$\alpha$，来更好地描述能带结构，当然$\alpha$也是一个经验系数， $$\alpha =0.5\times \frac{E_g(0)}{E_g(T)}$$ 一般取为0.5就足够了. 我们可以看一下在预测态密度上的区别，non-parabolic还是要好很多的. 注意态密度只是电子可能存在的态的数目，并不是电子的数目，电子可以不填充在相应的态上.

态密度的计算公式为

$$\text{DOS}(E)=\mathrm{d}Z/\mathrm{d}E$$ 能量为$E$时的态密度定义为$E+\mathrm{d}E$和$E$这两个等能面之间围成的态的数量. 而根据量子力学，相空间中一个量子态占据的体积是固定的，为${\hslash }^r$，其中$r$为体系的维数. 代入parabolic能带和non-parabolic能带，可以分别得到 $$\begin{array}{rl}& g_{\text{para}}(E_k)=\frac{m_d^{3/2}}{\sqrt{2}{\pi }^2{\hslash }^3}\sqrt{E_k}\\ & g_{\text{non-para}}(E_k)=\frac{m_d^{3/2}}{\sqrt{2}{\pi }^2{\hslash }^3}\sqrt{E_k(1+\alpha E_k)}(1+2\alpha E_k)\end{array}$$ 到这里我们对硅的能带结构的描述已经很比较完全了，那么按照这样的结构我们应该怎么样初始化一些电子出来呢？我们知道处于热平衡状态的电子的平均能量是$\frac{3}{2}{k}_{B}T$，假设现在这只电子处于某个导带最低点附近，那么在对应valley初始化电子的过程就可以描述为 $$\begin{array}{rl}{p}_x& =p_{0,x}+\sqrt{m_x^{\ast }{k}_{b}T}\times \text{Gassuian-Dis}\\ p_x& =p_{0,y}+\sqrt{m_y^{\ast }{k}_{b}T}\times \text{Gassuian-Dis}\\ p_x=p_{0,z}+\sqrt{m_z^{\ast }{k}_{b}T}\times \text{Gassuian-Dis}\end{array}$$ 而后可以代入能带结构获取对应动量Electron的能量， $$ϵ(\mathbf{k})=\frac{-1+\sqrt{1+4\alpha \gamma }}{2\alpha }$$ 其中 $$\gamma =\frac{{\hslash }^2}{2}[\frac{(k_x-k_{0,x}{)}^2}{m_x^{\ast }}+\frac{(k_y-k_{0,y}{)}^2}{m_y^{\ast }}+\frac{(k_z-k_{0,z}{)}^2}{m_z^{\ast }}]$$ 总觉得这里有些奇怪.. 不过先这样把. 用Plotly可以画一下我们得到的300K下Silicon中的电子分布.

import numpy as np
import pandas as pd

pi = 3.14159265
pi2 = 6.28318531
echarge = -1.0  # electron charge
ecoulom = -1.60217653e-19  # electron charge in Coulombs
kb = 8.61734318e-05  # eV/K
hbar = 6.58211915e-16  # eV*s
m0 = 5.68562975e-16  # eV*s^2/cm^2
eps0 = 5.52634972e05  # e/V/cm permittivity of free space

class Valley(object):
    def __init__(self, px0, py0, pz0, mx, my, mz, T, alpha, PQMAX, number_of_electrons):
        self.px0 = px0
        self.py0 = py0
        self.pz0 = pz0
        self.mx = mx
        self.my = my
        self.mz = mz
        self.T = T
        self.alpha = alpha
        self.PQMAX = PQMAX
        self.number_of_electrons = number_of_electrons
        electrons = []
        E = []
        for i in range(self.number_of_electrons):
            electron = self.create_electron(T)
            electrons.append(electron)
            E.append(self.energy(electron))
            
        self.electrons = pd.DataFrame(np.hstack([np.array(electrons), np.array(E).reshape(-1, 1)]), columns=['px', 'py', 'pz', 'E'])
        
    def create_electron(self, T):
        # create an electron in thermal equilibrium
        px = py = pz = 100
        while (px > self.PQMAX or py > self.PQMAX or pz > self.PQMAX):
            px = np.random.normal() * np.sqrt(kb * T * self.mx) + self.px0
            py = np.random.normal() * np.sqrt(kb * T * self.my) + self.py0
            pz = np.random.normal() * np.sqrt(kb * T * self.mz) + self.pz0
            electron = np.array([px, py, pz])
        return electron

    def gamma(self, electron):
        return 0.5 * (
            (electron[0] - self.px0) ** 2 / self.mx
            + (electron[1] - self.py0) ** 2 / self.my
            + (electron[2] - self.pz0) ** 2 / self.mz
        )
        
    def energy(self, electron):
        return (-1 + np.sqrt(1 + 4 * self.alpha * self.gamma(electron))) / (2 * self.alpha)
    
class Silicon(object):
    def __init__(self, T=300, number_of_electrons=600):
        # number of electrons
        electrons_per_valley = int(number_of_electrons / 6)
        self.number_of_electrons = 6 * electrons_per_valley
        
        # temperature
        self.T = T

        # lattice information
        self.asi = 5.411e-8  # lattic constant, cm

        # effective electron mass
        self.ml = 0.9163 * m0
        self.mt = 0.1880 * self.E_gap(0) * m0 / self.E_gap(self.T)

        # Dos electron mass
        self.md = 6 ** (2 / 3) * (self.ml * self.mt**2) ** (1 / 3)

        # brillouin zone information
        self.QMIN = 1.0
        self.QMAX = 2 * pi / self.asi
        self.PQMAX = hbar * self.QMAX  # hbar * q = p

        # Kane band
        self.alpha = 0.5 * self.E_gap(0) / self.E_gap(self.T)
        
        # six valleys
        self.valley_kx0 = Valley(-0.85 * self.PQMAX, 0, 0, self.ml, self.mt, self.mt, self.T, self.alpha, self.PQMAX, electrons_per_valley)
        self.valley_kx1 = Valley(0.85 * self.PQMAX, 0, 0, self.ml, self.mt, self.mt, self.T, self.alpha, self.PQMAX, electrons_per_valley)
        self.valley_ky0 = Valley(0, -0.85 * self.PQMAX, 0, self.mt, self.ml, self.mt, self.T, self.alpha, self.PQMAX, electrons_per_valley)
        self.valley_ky1 = Valley(0, 0.85 * self.PQMAX, 0, self.mt, self.ml, self.mt, self.T, self.alpha, self.PQMAX, electrons_per_valley)
        self.valley_kz0 = Valley(0, 0, -0.85 * self.PQMAX, self.mt, self.mt, self.ml, self.T, self.alpha, self.PQMAX, electrons_per_valley)
        self.valley_kz1 = Valley(0, 0, 0.85 * self.PQMAX, self.mt, self.mt, self.ml, self.T, self.alpha, self.PQMAX, electrons_per_valley)

    def E_gap(self, T):
        return 1.1756 - 8.8131e-5 * T - 2.6814e-7 * T**2

    def DoS(self, E):
        return (
            self.md ** (3 / 2) * np.sqrt(2)
            / (pi**2 * hbar**3)
            * np.sqrt(E * (1 + self.alpha * E))
            * (1 + 2 * self.alpha * E)
        )

    def DoS_para(self, E):
        return (
            self.md ** (3 / 2) * np.sqrt(2)
            / (pi**2 * hbar**3)
            * np.sqrt(E)
        )

si = Silicon()