半导体中的电子运动及迁移率模拟

电子的运动

在经典的图像里，电子有质量，电子受到电场力产生了加速度，在加速度的作用下，电子速度增加发生了运动，但是这个图像在能带理论里已经摒弃了。实际上，在固体内电场力的作用引起了电子动量的变化，进而引起了电子波矢态\(\mathbf{k}\)的变化， \[ \mathbf{E}q = \frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t} = \hbar \frac{\mathrm{d}\mathbf{k}}{\mathrm{d}t} \] 波矢态\(\mathbf{k}\)的变化相应地使得电子波包群速度发生变化， \[ \mathbf{v}_g(\mathbf{k}) = \frac{1}{\hbar} \frac{\partial \epsilon}{\partial \mathbf{k}} = \frac{\partial \epsilon}{\partial \mathbf{p}} \] 波包群速度就是能带里能量对动量的偏导，波矢态\(\mathbf{k}\)确定了，电子运动的速度就确定了。

电子可以取到的波矢\(\mathbf{k}\)，离散地分布到允许的能带中。所以从净的效果来看，电子之所以能够被电场力加速，就是因为它的波矢态是可以在动量守恒的条件下、在电场力的作用下发生相应变化的。根据泡利不相容原理，不能有两个或两个以上的电子处于完全相同的态，如果所有允许的态都被电子填满了，电子的波矢无法发生相应变化，自然也就无法被电场力加速了，也就没有电流了，所以宏观体现为不导电。当然实际过程可能是说有一个电子被电场力加速，它的态发生了变化，跑到了另一个态上，另一个态的电子无处可去，只能回到这个电子之前的态上，宏观体现为外加了电场，里面电子的状态没有发生任何改变，没有电子被电场力加速，自然也就没有电流了，这里已经完全不再类似于经典的受力加速图像了。电子在半导体中受电场力的驱动在实空间中运动，但是在半导体的某个位置电子所能取到的状态，受到当地的能带结构约束，而并不是像真空的自由电子一样可以随意的被加速和移动，这也是为什么必须借助能带图来理解电子在某个结构中的分布和运动。

迁移率模拟

用蒙特卡洛模拟来计算电子迁移率，可以非常直观地认识迁移率、能带和色散关系的物理意义，以及电子到底是如何在半导体内运动的，郝跃老师的书里对氮化镓中的电子在电场力下运动的过程有一些描述：

GaN电子的漂移速度首先随外电场强度的增大而上升，随后在一个较强的电场强度达到峰值（对于纤锌矿结构，约180kV/cm）；电场继续增强，则出现一个速度逐渐下降的广阔负微分电阻区，直到在足够强的电场下速度达到电子饱和速度（约\(1.5\times 10^7\text{cm/s}\)）。由于具有较高的击穿电场强度（GaN为\(3\times 10^6\text{V/cm}\)以上，AlN为\(12\times 10^6\text{V/cm}\)以上），GaN材料和AlN材料可以承受很强的电场。低场下电子都位于导带中最低的\(\Gamma\)谷，速度和电场呈近似线性关系，其比例系数即为低场迁移率；当电场增强时，电子能量增加，有机会进入能量较高、有效质量也较大的卫星谷，令迁移率总体呈下降趋势，但GaN中\(\Gamma\)谷中的电子有效质量较大、室温下极性光学声子散射左右很强以及\(\Gamma\)谷和卫星谷之间较大的能量间距，使得电子的谷间转移在较高的电场下才能大规模出现，此时电子的漂移速度达到峰值。当电子进入卫星谷后，其动能和速度减小，有效质量变大，令电子漂移能力变弱，形成负微分电阻特性；电场继续增加时，电子继续向卫星谷转移，电子能量仍上升，但有效质量大且新的谷间形变势散射成为有效的散射机制，是的电子能量随电场增加的速度变慢。两者的综合作用导致形成了广阔的负微分电阻区，最后电子速度达到饱和。

电子在给定半导体的能带中受到电场力作用，动量改变，引起波矢\(\mathbf{k}\)的变化，因为电子能够取到的态分布在导带上，这个过程沿着某一个能带进行，此时这个电子的能量和动量关系（即决定了它在电场力下如何运动的关系，色散关系），由它当前所处的能带位置决定（比如下面图中的不同能带即不同的线上，能量动量关系是不一样的）。在运动的过程中，电子可能会发生散射，跳跃到同一个能带中的其他位置（谷内散射，intravalley scattering）或者跳跃到不同的能带中（谷间散射，intervalley scattering）。注意这种跳跃只是电子的状态发生了跳跃，电子在实空间中的位置是不发生改变的。散射后，电子运动的能量动量关系由跳跃后抵达的能带位置决定，在电场力的作用下沿着新的能带运动，直到发生下一次散射，在这个过程中，一般认为导带是足够空的。所以从这个讨论可以看到，和真空中的自由电子完全不同，并不是电子能量越高则电子的速度越高，因为能带不再是一个简单的抛物线了。什么情况下电子的迁移率会大一些，就是当它运动在一个有效质量低（抛物线比较陡峭）且散射强度也较低的谷中，如果电子的能量很高但运动在一个有效质量大（抛物线比较平缓）的谷中，它的迁移率实际上很低。

用蒙特卡洛方法模拟体材料迁移率的程序，总的来说可以为两个部分，

• 模拟前：准备好半导体的能带信息和各散射机制的散射率，以数组形式存储
• 模拟中：初始化电子们，在每一时间步中对每一个电子处理电场加速和散射过程，同时记录当前步的电子信息

体材料迁移率的模拟只发生在动量空间，不涉及实空间的边界条件，不需要自洽地求解泊松方程，这简化了模拟流程，而且在迁移率模拟中虽然本质上叫作Ensemble模拟，实际上各个电子的运动是彼此独立的。原始的程序来自于https://nanohub.org/wiki/BMC，亚利桑那州立大学的Dragica Vasileska老师在2000年左右开源供大家学习的Fortran程序，这里简单地用Python针对硅重新写了一遍，没有采用像numba或者向量化这些加速操作，大概比原始的Fortran程序慢了300倍左右。

模拟准备

能带

可以参考硅的能带结构 - a Monte Carlo view，这里以硅的非抛物Kane带为例，也可以替换成一些DFT计算的结果，实际上原理都是一样的。硅一共有六个等价的导带，导带底部可以用椭球来描述。

对于第s个导带椭球，它的表达式可以写成 \[ E_k^s\left(1+\alpha E_k^s\right)=\frac{\hbar^2}{2}\left[\frac{\left(k_x-k_{0, x}^s\right)^2}{m_x^*}+\frac{\left(k_y-k_{0, y}^s\right)^2}{m_y^*}+\frac{\left(k_z-k_{0, z}^s\right)^2}{m_z^*}\right] = \gamma \] 其中\(k_{0, x}^s, k_{0, y}^s, k_{0, z}^s\)代表这个椭球的中心坐标，\(\alpha\)称为非抛物因子，这里取为\(0.5 \,\mathrm{eV}^{-1}\)，这里在计算的时候要非常小心，因为\(\alpha\)是带能量量纲的，semiconductor people喜欢用电子伏特和厘米作为单位制，如果使用标准单位计算的话，\(\alpha\)的值要除以一个电子伏特。沿椭球方向的等效质量设置为\(0.91m_0\)，垂直椭球方向的等效质量设置为\(0.19m_0\)。当然椭球可能不太好处理，我们可以简单地做一下坐标变换，我们定义一个电导有效质量\(m_c\)，

\[ m_c = \frac{3}{\frac{1}{m_l} + \frac{2}{m_t}} \] 假设对于长轴位于x轴的椭球，可以定义\(k_x = k_x^\prime \sqrt{m_l / m_c}, k_y = k_x^\prime \sqrt{m_t / m_c}, k_z = k_z^\prime \sqrt{m_t / m_c}\)，且由于电子的性质只和偏离导带底部的位置有关，我们可以只考虑偏离导带底的波矢部分，把所有的椭球都平移到布里渊区的零点，同时也不需要真的模拟六个椭球导带，只需要模拟三个互相正交的椭球导带就可以了，因为波矢的绝对值不影响电子性质，只要导带椭球在布里渊区里摆放的方向是一样的，实际上电子在里面运动就是一样的。将上述坐标变换代回能带的表达式中，得到

这会给波矢初始化和能量评估带来一些简便性。假设我们已经得到电子的波矢值了，那么可以解出来能量的表达式为 \[ E = \frac{-1 + \sqrt{1 + 4\alpha \gamma} }{2\alpha} = \frac{2\gamma}{\sqrt{1 + 4\alpha \gamma} + 1} \] 因为\(\sqrt{1+4\alpha\gamma}\)和1的差距很小，我们在数值计算里要避免这种两个相差很小的数相减的情况，因此上面的表达式转换成了最右边的式子以保持数值精度。

群速度表达式为 \[ v_x =\frac{1}{\hbar} \frac{E}{\partial k_x} = \frac{\hbar k_x^\prime}{1+2\alpha E} \sqrt{1/(m_l m_c)} \]

此时假如在x方向的电场力作用下运动了\(\tau\)时间，那么x方向的经坐标变换后的波矢变化就是 \[ \mathrm{d}k_x^\prime = - \frac{1}{\hbar}E_x q \tau \sqrt{m_c / m_l} \] #### 散射率

实际模拟过程中为了简化散射率的计算，一般是按照能量分成小段，把散射率储存成查找表，而不是每对应一个状态就重新计算一次散射率。比如我们假设电子在模拟过程中最高可以取到的能量是1eV，把\(0\sim 1\,\mathrm{eV}\)离散成n_lev=1000个小段，把每一段的各个散射机制对应的散射率都储存成一个数组，散射率这个数组中抽样，超出1eV的电子散射率也都按照1eV去处理。可以参考Electron Monte Carlo (1) - 散射，这里只选取了五种散射机制，包括g-intervalley吸收散射，g-intervalley发射散射，f-intervalley吸收散射，f-intervalley发射散射，以及谷内声学声子散射，没有考虑离子库伦散射。这五个散射都可以认为是各向同性散射，这里的各向同性指的是相对椭球来说的各向同性，也就是坐标变换后的\(k_x^\prime, k_y^\prime, k_z^\prime\)是各向同性的，可以按照各向同性抽样去抽取散射后新的波矢大小，这里可以参考蒙特卡洛中的采样。

下面时各种散射率随着电子能量变化的示意图，这里的能量均以导带底作为原点。从各种机制下散射率的大小可以清楚地看到电子的运动在电场力和晶格散射的作用下是如何达到平衡的。起初，电子能量较低，散射率较低，电子自由飞行，从电场中获得能量，能量增加。随着能量的增高，散射率逐渐增加，在各种散射机制中，发射声子的g/f的谷间散射占据了绝对优势，因此每发生一次散射，电子大概率发射一颗光学声子，此时电子的能量降低。当两种机制达到平衡时，即从电场中获得能量收益和发射声子的能量损失平衡时，电子的运动也就达到了平衡，此时就可以统计电子的平均漂移速度来计算迁移率了。

模拟过程

初始化

处于热平衡状态的电子的平均能量是\(\frac{3}{2}k_BT\)，假设系统处于300K，那么先按照高斯分布抽样电子的能量，能量确定后，依据能量抽样电子波矢，最后随机抽样电子所处的valley。初始化完所有的电子以后，对每一个电子抽样初次的自由飞行时间。

漂移-散射

整个Ensemble蒙特卡洛模拟的逻辑是这样的，模拟沿着时间推进，每一个时间步的间隔是\(\mathrm{d}t\)，在每一个时间步里模拟这个电子的行为。首先，对于每一个电子，我们要抽取它在发生散射前的飞行时间，由下面这个式子确定， \[ \tau = -\ln(r) /\Gamma_m \] 其中\(\Gamma_m\)是一个人为选取的固定值，这个数值大于模拟体系中可能出现的最大散射率，这里引出了电子蒙特卡洛模拟中重要的自散射（self-scattering）的概念，相关的讨论可以参考Electron Monte Carlo (1) - 散射。抽样出了\(\tau\)以后，这时候根据\(\tau\)的数值存在两种可能性，如果\(\tau>\mathrm{d}t\)，这说明在当前这个时间步里这颗电子没有被散射，它一直在被电场加速，此时，新的\(\tau = \mathrm{d}t - \tau\)。如果\(\tau < \mathrm{d}t\)，这说明这个电子被电场加速了\(\tau\)的时间后，发生了散射，这时在模拟中，我们让这个自由电子自由飞行\(\tau\)的时间，然后根据它飞行后的能量，抽样散射的机制，根据对应的机制，更新散射后这颗电子的波矢、能量、所处的valley信息，之后，再次抽样它的自由飞行时间\(\tau_\text{new}\)，如果\(\tau_\text{new} + \tau > \mathrm{d}t\)，那么此时的情况就和第一种情况一样了，这颗电子在电场的作用下自由飞行完剩下的\(\mathrm{d}t - \tau\)的时间，在当前时间步结束后距离下一次散射的时间就变成了\(\tau = \tau_\text{new} + \tau - \mathrm{d}t\)。如果还是\(\tau_\text{new}<\mathrm{d}t-\tau\)，那么就重复这个抽样-飞行-散射的过程，直到\(\tau_\text{new} > \mathrm{d}t - \tau\)。

Code

MC class

import os
import sys

import tqdm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import constants as c

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.style.use('science')

class MC(object):
    def __init__(
        self,
        number_of_electrons=1000, # 模拟电子数
        E_max=1*c.e, # 最大电子能量，用于确定散射率查找表
        number_of_bands=1000, # 电子能量离散的段数，用于确定散射率查找表
        simulation_time=5e-12, # 模拟总时间
        dt=1e-14, # 步长
        electric_x=0, # x方向电场强度
        electric_y=1e5, # y方向电场强度
        electric_z=0, # z方向电场强度
    ):
        # 全部采用SI国际单位制
        self.t = 0 
        self.number_of_electrons = number_of_electrons
        self.E_max = E_max
        self.number_of_bands = number_of_bands
        self.E_array = np.linspace(0, self.E_max, self.number_of_bands)
        self.simulation_time = simulation_time
        self.dt = dt
        self.electric_x = electric_x
        self.electric_y = electric_y
        self.electric_z = electric_z
        self.valley_indexes = np.array([0, 1, 2]) # 用3个正交的valley替代原始的6个valley

        self.rho = 2329 # 密度
        self.ml = 0.91 * c.m_e # 椭球纵向有效质量
        self.mt = 0.19 * c.m_e # 椭球横向后效质量
        self.md = (self.ml * self.mt * self.mt) ** (1 / 3) # 态密度有效质量
        self.mc = 3 / (1 / self.ml + 2 / self.mt) # 电导有效质量
        self.alpha = 0.5 / c.e # 非抛物系数
        self.T = 300 # 温度
        self.vs = 9040 # 声速

        # Coordinate transformation
        self.factor_l = np.sqrt(self.mc / self.ml)
        self.factor_t = np.sqrt(self.mc / self.mt)

        # scattering
        self.Deformation_potential_g = 5.23e10 * c.e
        self.Deformation_potential_f = 5.23e10 * c.e
        self.phonon_energy_g = 63e-3 * c.e
        self.phonon_energy_f = 59e-3 * c.e
        self.number_of_g_valleys = 1
        self.number_of_f_valleys = 4

        # intra, g_em, g_ab, f_em, f_ab分别是各种机制的散射率数组
        intra = self.intravalley_acoustic(self.E_array)
        g_em, g_ab = self.intervalley(
            self.E_array,
            self.Deformation_potential_g,
            self.phonon_energy_g,
            self.number_of_g_valleys,
        )
        f_em, f_ab = self.intervalley(
            self.E_array,
            self.Deformation_potential_f,
            self.phonon_energy_f,
            self.number_of_f_valleys,
        )

        self.stacked_scattering_rate = np.vstack((intra, g_em, g_ab, f_em, f_ab)).T
        
        # 计算累积分布，用于抽样
        self.cumsum_scattering_rate = np.cumsum(self.stacked_scattering_rate, axis=1)
        
        # 确定最大总散射率，用于确定抽样时间
        self.max_cumsum_scattering_rate = np.max(self.cumsum_scattering_rate)
        
        # 用最大总散射率归一化累积分布，用于抽样
        self.normalized_cumsum_scattering_rate = (
            self.cumsum_scattering_rate / self.max_cumsum_scattering_rate
        )

    # 计算态密度
    def DOS(self, E):
        dos = (
            self.md ** (3 / 2)
            * np.sqrt(2)
            / (np.pi**2 * c.hbar**3)
            * np.sqrt(E * (1 + self.alpha * E))
            * (1 + 2 * self.alpha * E)
        )
        dos[np.isnan(dos)] = 0
        return dos

    # 玻色分布，用于确定散射率
    def bose(self, E_phonon):
        return 1 / (np.exp(E_phonon / (c.k * self.T)) - 1)

    # 计算谷间声子散射的散射率
    def intervalley(self, E, D, E_phonon, Z_iv):
        scatOconst = np.pi / (2 * self.rho / c.hbar)
        Nq = self.bose(E_phonon)
        rate_em = (
            Z_iv * D**2 * (Nq + 1) * self.DOS(E - E_phonon) * scatOconst / E_phonon
        )
        rate_abs = Z_iv * D**2 * (Nq) * self.DOS(E + E_phonon) * scatOconst / E_phonon
        return rate_em, rate_abs

    # 计算谷内声学声子散射的散射率
    def intravalley_acoustic(self, E):
        D = 6.55 * c.e
        return (
            np.pi
            * c.k
            * self.T
            / self.rho
            / c.hbar
            / self.vs**2
            * D**2
            * self.DOS(E)
        )

    # 给定电子波矢，根据Kane带计算电子能量
    def get_electron_energy(self, kx, ky, kz):
        gamma = c.hbar**2 / 2 * (kx**2 + ky**2 + kz**2) / self.mc
        energy = (4 * self.alpha * gamma) / (
            2 * self.alpha * (np.sqrt(1 + 4 * self.alpha * gamma) + 1)
        )
        return energy

    # 给定电子波矢、能量、所处的valley，计算电子速度
    def get_electron_velocity(self, kx, ky, kz, energy, valley_index):

        if valley_index == 0:
            vx = (
                c.hbar
                * kx
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.ml))
            )
            vy = (
                c.hbar
                * ky
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.mt))
            )
            vz = (
                c.hbar
                * kz
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.mt))
            )
        elif valley_index == 1:
            vx = (
                c.hbar
                * kx
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.mt))
            )
            vy = (
                c.hbar
                * ky
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.ml))
            )
            vz = (
                c.hbar
                * kz
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.mt))
            )
        else:
            vx = (
                c.hbar
                * kx
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.mt))
            )
            vy = (
                c.hbar
                * ky
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.mt))
            )
            vz = (
                c.hbar
                * kz
                / (1 + 2 * self.alpha * energy)
                * np.sqrt(1 / (self.mc * self.ml))
            )

        return vx, vy, vz

    # 根据给定的能量，各向同性地抽样电子波矢
    def generate_electron_vector(self, energy):
        gamma = energy * (1 + self.alpha * energy)
        k = np.sqrt(2 * gamma / c.hbar**2 * self.mc)
        phi = 2 * np.pi * np.random.random()
        cos_theta = 1 - 2 * np.random.random()
        sin_theta = np.sqrt(1 - cos_theta**2)
        kx = k * sin_theta * np.cos(phi)
        ky = k * sin_theta * np.sin(phi)
        kz = k * cos_theta
        return kx, ky, kz

    # 漂移过程，给定当前时刻的电子波矢，所处valley，计算电场加速acceleration_time后电子的性质
    def drift(self, kx, ky, kz, valley_index, accleration_time):

        valley_index = int(valley_index)

        if valley_index == 0:
            dkx = -self.electric_x * c.e * accleration_time * self.factor_l / c.hbar
            dky = -self.electric_y * c.e * accleration_time * self.factor_t / c.hbar
            dkz = -self.electric_z * c.e * accleration_time * self.factor_t / c.hbar

        elif valley_index == 1:
            dkx = -self.electric_x * c.e * accleration_time * self.factor_t / c.hbar
            dky = -self.electric_y * c.e * accleration_time * self.factor_l / c.hbar
            dkz = -self.electric_z * c.e * accleration_time * self.factor_t / c.hbar

        else:
            dkx = -self.electric_x * c.e * accleration_time * self.factor_t / c.hbar
            dky = -self.electric_y * c.e * accleration_time * self.factor_t / c.hbar
            dkz = -self.electric_z * c.e * accleration_time * self.factor_l / c.hbar

        kx = kx + dkx
        ky = ky + dky
        kz = kz + dkz

        energy = self.get_electron_energy(kx, ky, kz)
        return kx, ky, kz, energy

    # 抽样自由飞行时间
    def sample_flight_time(self):
        return -np.log(np.random.random()) / self.max_cumsum_scattering_rate

    # 根据给定的电子能量，确定它的能量区间，用于确定散射率
    def find_electron_index(self, energy):
        return max(0, np.searchsorted(self.E_array, energy) - 1)

    # 散射过程，输入当前电子的状态，输出散射后电子的状态
    def scatter(self, kx, ky, kz, energy, valley_index):
        valley_index = int(valley_index)
        index = min(self.find_electron_index(energy), int(self.number_of_bands - 1))
        cumsum_scattering_rate = self.normalized_cumsum_scattering_rate[index]
        band = np.searchsorted(cumsum_scattering_rate, np.random.random())
        if band == 0:
            # intravalley_acoustic
            new_energy = energy
            new_valley_index = valley_index
        elif band == 1:
            # intervalley_g_emission
            new_energy = energy - self.phonon_energy_g
            new_valley_index = valley_index
        elif band == 2:
            # intervalley_g_absorption
            new_energy = energy + self.phonon_energy_g
            new_valley_index = valley_index
        elif band == 3:
            # intervalley_f_emission
            new_energy = energy - self.phonon_energy_f
            new_valley_index = np.random.choice(
                np.delete(self.valley_indexes, valley_index)
            )
        elif band == 4:
            # intervalley_f_absorption
            new_energy = energy + self.phonon_energy_f
            new_valley_index = np.random.choice(
                np.delete(self.valley_indexes, valley_index)
            )
        else:
            # self scatter, nothing happens
            return kx, ky, kz, energy, valley_index

        return *self.generate_electron_vector(new_energy), new_energy, new_valley_index

    # 漂移-散射
    def flight_scatter(self, kx, ky, kz, energy, valley_index, tau):
        if tau > self.dt:

            # 此时在整个dt时间内没有发生散射，一直被电场力加速
            kx, ky, kz, energy = self.drift(kx, ky, kz, valley_index, self.dt)
            tau = tau - self.dt
            return kx, ky, kz, energy, valley_index, tau
        else:
            # 直到tau_new > dt - tau, 让电场力加速完整个dt周期

            kx, ky, kz, energy = self.drift(kx, ky, kz, valley_index, tau)
            kx, ky, kz, energy, valley_index = self.scatter(
                kx, ky, kz, energy, valley_index
            )
            tau_new = self.sample_flight_time()

            while (tau_new + tau) < self.dt:
                tau += tau_new
                kx, ky, kz, energy = self.drift(kx, ky, kz, valley_index, tau_new)
                kx, ky, kz, energy, valley_index = self.scatter(
                    kx, ky, kz, energy, valley_index
                )
                tau_new = self.sample_flight_time()

            left_acceleration_time = self.dt - tau
            tau = tau + tau_new - self.dt
            kx, ky, kz, energy = self.drift(
                kx, ky, kz, valley_index, left_acceleration_time
            )
            return kx, ky, kz, energy, valley_index, tau

    # 初始化电子
    def initialize_electron(self):
        energy = -1.5 * c.k * self.T * np.log(np.random.random())
        valley_index = np.random.choice(self.valley_indexes)
        kx, ky, kz = self.generate_electron_vector(energy)
        vx, vy, vz = self.get_electron_velocity(kx, ky, kz, energy, valley_index)
        tau = self.sample_flight_time()
        return kx, kz, kz, energy, valley_index, tau, vx, vy, vz

    # 开始模拟
    def run(self):
        self.t = 0
        self.electrons = []
        self.hists = []
        
        # 初始化生成self.number_of_electrons个电子
        for i in range(self.number_of_electrons):
            self.electrons.append(self.initialize_electron())
        self.electrons = np.array(self.electrons)

        # 开始模拟
        print('begin simulation... (∪.∪ )...zzz')
        for i_t in tqdm.tqdm(range(int(2 + self.simulation_time / self.dt))):
            for i in range(self.number_of_electrons):
                # flight-scatter process
                self.electrons[i, :6] = self.flight_scatter(*self.electrons[i, :6])
                self.electrons[i, 6:9] = self.get_electron_velocity(
                    *self.electrons[i, 0:5]
                )
            # save current distribution
            self.hists.append(
                (
                    self.t,
                    np.average(self.electrons[:, 3]), # energy
                    np.average(self.electrons[:, 6]), # vx
                    np.average(self.electrons[:, 7]), # vy
                    np.average(self.electrons[:, 8]), # vz
                )
            )
            self.t += self.dt
        print('simulation ends... （＾∀＾●）ﾉｼ')
        self.hists = np.array(self.hists)

a single run

mc = MC(simulation_time=5e-12, electric_y=5e6)
mc.run()

begin simulation... (∪.∪ )...zzz

100%|██████████| 501/501 [00:16<00:00, 30.79it/s]

simulation ends... （＾∀＾●）ﾉｼ

# scattering calculation

plt.plot(mc.E_array / c.e, mc.stacked_scattering_rate[:, 0])
plt.plot(mc.E_array / c.e, mc.stacked_scattering_rate[:, 1])
plt.plot(mc.E_array / c.e, mc.stacked_scattering_rate[:, 2])
plt.plot(mc.E_array / c.e, mc.stacked_scattering_rate[:, 3])
plt.plot(mc.E_array / c.e, mc.stacked_scattering_rate[:, 4])

plt.text(0.5, 3e13, r'$\mathbf{f-em}$', color='purple')
plt.text(0.5, 9.9e12, r'$\mathbf{g-em}$', color='red')
plt.text(0.3, 4e12, r'$\mathbf{intravalley \; acoustic}$', color='blue')
plt.text(0.5, 1.2e12, r'$\mathbf{f-ab}$', color="#FF7F00" )
plt.text(0.5, 2.3e11, r'$\mathbf{g-ab}$', color="green" )

plt.xlabel(r'$e\, (\mathrm{eV})$')

plt.ylim(1e10, 1e14)
plt.yscale('log')
plt.ylabel(r'$\Gamma\, (\mathrm{s}^{-1})$')

plt.savefig('scattering.png')

能量随时间的演化：

plt.plot(mc.hists[:, 0]/1e-12, mc.hists[:, 1]/c.e)

plt.xticks([0, 1, 2, 3, 4, 5])
plt.xlabel(r'$t\, (\mathrm{ps})$')

plt.ylim(0, 0.3)
plt.yticks([0, 0.1, 0.2, 0.3])
plt.ylabel(r'$e\, (\mathrm{eV})$')

plt.savefig('energy.png')

漂移速度随时间的演化：

plt.plot(mc.hists[:, 0]/1e-12, mc.hists[:, 2]/1e5, label=r'$v_x$')
plt.plot(mc.hists[:, 0]/1e-12, mc.hists[:, 3]/1e5, label=r'$v_y$')
plt.plot(mc.hists[:, 0]/1e-12, mc.hists[:, 4]/1e5, label=r'$v_z$')

plt.xticks([0, 1, 2, 3, 4, 5])
plt.xlabel(r'$t\, (\mathrm{ps})$')

plt.ylim(-3, 1)
plt.ylabel(r'$v \, (10^5 \, \mathrm{m/s})$')

plt.legend()

plt.savefig('velocity.png')

迁移率模拟

fys = [2e4, 5e4, 1e5, 5e5, 1e6, 5e6, 7e6]
mcs = [MC(simulation_time=5e-12, electric_y=fy) for fy in fys]
for mci in mcs:
    mci.run()
    
velocities = [- np.average(mci.hists[200:, 3]) for mci in mcs]

plt.plot(fys, velocities, marker='s', markerfacecolor='None', markeredgecolor='black', linestyle='--', color='black')

plt.xlim(1e4, 1e7)
plt.xlabel(r'electric field$\, (\mathrm{V/m}) $')
plt.xscale('log')

plt.ylabel('drift velocity$\, (\mathrm{m/s})$')
plt.ylim(1e3, 2e5)
plt.yscale('log')

plt.savefig('mobility.png')