目录

载流子迁移率

关键词：迁移率，电声耦合（电子-声子相互作用），温度相关载流子输运性质，电导率（conductivity）

在本教程中，我们将介绍如何计算石墨烯（graphene）中的声子限制的电子迁移率。结合电子态和玻尔兹曼输运方程，同时考虑声子模式与电子-声子相互作用，使用 DFT 方法计算电子迁移率。我们将使用 QuantumATK 新版中提供的分析模块 ElectronPhononCoupling 和 Mobility 功能来完成。

计算迁移率 $\mu$ 的关键是计算电子的弛豫时间 $\tau$。弛豫时间有两种计算方法：

  1. 全角度(k,q)-相关：这种方法完整的考虑 $\tau$ 对电子和声子的波矢 $\mathbf{k}$ 和 $\mathbf{q}$ 的相关，即 $\tau=\tau(\mathbf{k},\mathbf{q})$。下面称为 (k,q)-相关方法。
  2. 各向同性散射：这种方法仅考虑能量有关的 $\tau=\tau(E)$，假定 $\tau(E)$ 在动量空间中的变化是各向同性的。下面称为 E-相关方法。

下面的计算表明，两种方法得到的结果相差无几，但是第二种方法比第一种方法计算速度大大提高。

文末将介绍理论背景知识和更多相关信息，已发表结果参见[GMSB16]。

提示

本教程使用特定版本的QuantumATK创建，因此涉及的截图和脚本参数可能与您实际使用的版本略有区别，请在学习时务必注意。

结构弛豫与电子态计算

在 Builder 中添加一个石墨烯结构：

设置晶格常数C方向为20 Å：

将结构在视图中居中（Coordinate Tools → Center → Apply），发送至脚本生成器Script Generator。

要进行精确的弛豫时间和迁移率计算，首先要精确的优化A、B方向的晶格常数并计算能带。

计算结束后，点击 LabFloor 中包含在 Graphene_relax.hdf5 的 Bandstructure 图标，使用右侧的 Bandstructure Analyzer 可以查看能带图。

将鼠标防止在一条能带上就可以看到相关的信息。比如上图中高亮的一条价带的编号为3，导带为4。这两条能带是与迁移率有关的关键能带，接下来我们将专注于这两条能带。

石墨烯中的声子

下面我们将计算石墨烯的动力学矩阵（dynamical matrix）。由动力学矩阵，可以得到声子谱和振动模式，从而可以测验计算的精确度。

步骤如下：

在 Sciprt Generator 主面板中，将 Default output file 改为：Graphene_dynmat.hdf5，脚本发送到 Job Manager，保存输入文件Graphene.py，完成计算。

为确保石墨烯电子态单胞计算时的k点采样与此处声子计算时的超胞一致，有必要设置 $\mathrm{k_i^{unit cell}} = \mathrm{k_i^{supercell}} \times \mathrm{N_i^{supercell}}$ ($\mathrm{i} = \mathrm{A},\mathrm{B},\mathrm{C}$)。

由于$11 \times 11$超胞尺寸较大（242原子），因此计算动力学矩阵非常耗时。不过，计算可以将原子位移并行化。由于单胞包含两个原子，每个原子在x，y，z三个方向上位移，因此有6个计算要完成。当总的计算数目乘以计算所用参数Proecesses per displacement等于计算所用的总核数时，计算效率最高。比如此例子中，如果设置Proecesses per displacement = 4，计算在24核心上完成只需要 30 分钟。

计算结束后，使用 Bandstructure Analyzer 查看声子谱结构，如下：

图34. 声子谱结构。有 3 个声学模式：最下方一条，平面外模（ZA）在长波极限下∝q2，另外两条在小波矢附近与 q 线性相关。横声学模（TA）的能量低于纵声学模（LA）。

石墨烯的迁移率

接下来，我们来看如何计算石墨烯中的电子迁移率。在已经计算的电子态结构和动力学矩阵的基础上，要计算迁移率还需要三步：

  1. 计算哈密顿量导数
  2. 计算电子-声子耦合
  3. 计算迁移率

其中第2、第3步分别有(k,q)-相关方法和E-相关方法两种，下文分别介绍。

哈密顿量微分

为了计算电子-声子耦合，需要首先计算哈密顿量的一阶偏微分。哈密顿量对第 $i$ 个原子的笛卡尔坐标$\alpha$的一阶偏微分，$\partial \hat{H}/ \partial R_{i,\alpha}$，可以使用 HamiltonianDerivatives 对象计算得到。

现在，按下面步骤操作：

注意

此处的超胞重复次数需要和计算动力学矩阵的一致。

发送任务到任务管理器 Job manager，保存输入文件graphene_dHdR.py，运行计算任务。

注意

与动力学矩阵类似，哈密顿量导数计算相当耗时，但是可以将原子位移在多核上并行。此例子在24核上并行的计算时间大约1小时。

(k,q)-相关方法

要计算弛豫时间 $\tau(\mathbf{k},\mathbf{q})$ 和迁移率 $\mu$，需要在 $\mathbf{k}$- 和 $\mathbf{q}$-点的精细网格上计算电子-声子耦合。

电声耦合

单击Script Generator，添加如下计算单元：

以下两个计算单元自动添加：

在本例里，由于动力学矩阵和哈密顿量微分都已经计算过了，因此可以将这两个计算删去，添加两个Analysis from File，直接使用已有的结果。

打开第二个Analysis from File，选择Graphene_dynmat.hdf5，加载DynamicMatrix数据对象，如下图：

打开第三个Analysis from File，选择Graphene_dHdR.hdf5，加载HamiltonianDerivatives数据对象，如下图：

最后在ElectronPhononCoupling里设置如下参数：

k和q的采样范围选择为仅在布里渊区的相关区域。对于前者是在$K$-点（[1.7063,0.0,0.0]）附近，后者是在$\Gamma$-点（[0.0,0.0,0.0]）附近。

在Script Generator主面板里设置输出文件为Graphene_eph.hdf5，将脚本发送至Editor，在下面一段按如下进行修改： 

 electron_phonon_coupling = ElectronPhononCoupling(
     configuration=bulk_configuration,
     dynamical_matrix=dynamical_matrix,
     hamiltonian_derivatives=hamiltonian_derivatives,
     kpoints_cartesian=kpoints,
     qpoints_cartesian=qpoints,
     electron_bands=[3,4],
     phonon_modes=All,
     energy_tolerance=0.01*eV,
     initial_state_energy_range=[-0.5,0.5]*eV,  //修改这一行
     store_dense_coupling_matrices=False,       //修改这一行
     )

对$k$- 和 $q$-空间进行高密度采样意味着计算非常耗时。合理设置energy_tolerance和initial_state_energy_range两个参数可以让计算仅考虑相关的电子-声子耦合矩阵元， 这样可以大大减少计算时间和对内存的需求。本例的计算在24核服务器上大约需要1.5小时。有关这两个参数的详细说明参见手册：ElectronPhononCoupling

将脚本发送到Job Manager，保存输入文件为Graphene_eph.py，完成计算。

迁移率

在LabFloor里，选择Graphene_relax.hdf5文件，将其中的Object Id ‘gID001’发送至Script Generator，添加如下计算单元，并修改参数：

Fermi shift 0.13 eV相当于载流子浓度$n=10^{12}\mathrm{cm}^{-2}$

在Script Generator主面板里，将 Default output file 设为Graphene_mu-full.hdf5，将脚本发送至Job Manager，保存输入文件Graphene_mu-full.py，完成计算。

计算一旦完成，在LabFloor里Graphene_mu-full.hdf5，选择$\mu_{e,h}$图标，点击Text Representation。 

 +------------------------------------------------------------------------------+
 | Mobility Report                                                              |
 | ---------------------------------------------------------------------------- |
 | Input parameters:                                                            |
 |     Temperature       = 300.00 K                                             |
 |     Fermi level shift = 0.13 eV                                              |
 |     Energy broadening = 0.0030 eV                                            |
 |     q-grid refinement = 1                                                    |
 |                                                                              |
 +------------------------------------------------------------------------------+
 | Trace of linear responce tensors:                                            |
 +------------------------------------------------------------------------------+
 |                                                                              |
 | Electrons:                                                                   |
 |                                                                              |
 |     Mobility                 =  2.56e+05 cm^2/(V*s)                          |
 |     Conductivity             =  7.71e+06 S/m                                 |
 |     Seebeck coefficient      = -6.11e-05 V/K                                 |
 |     Thermal conductivity     =  4.46e+01 W/(m*K)                             |
 |     Carrier density (2D, xy) =  3.76e+11 cm^-2                               |
 |                                                                              |
 |                                                                              |
 | Holes:                                                                       |
 |                                                                              |
 |     Mobility                 =  2.56e+06 cm^2/(V*s)                          |
 |     Conductivity             =  4.27e+04 S/m                                 |
 |     Seebeck coefficient      =  1.81e-01 V/K                                 |
 |     Thermal conductivity     =  4.02e+00 W/(m*K)                             |
 |     Carrier density (2D, xy) =  2.09e+08 cm^-2                               |
 |                                                                              |
 +------------------------------------------------------------------------------+

计算得到的电子迁移率为$2.56 \cdot 10^{5}\ \mathrm{cm^{2} V^{-1} s^{-1}}$，与过去文献中在载流子浓度为$n = 10^{12} \mathrm{cm^{-2}}$和室温时的结果一致。

E-相关方法

电声耦合

下面我们用计算能量相关的弛豫时间的$\tau(E)$的方法来计算迁移率，为此我们需要重新计算电声耦合。我们假定弛豫时间在$\mathrm{k}$-空间是各向同性的。这意味着我们只需要从Dirac点（K点）向外计算一条线段上的电声耦合就足够了，很显然这样可以大大降低计算量。

单击Script Generator，添加如下计算单元：

与上一节一样从Graphene_dynmat.hdf读入动力学矩阵，从Graphene_dHdR.hdf5，读入哈密顿量微分。

在ElectronPhononCoupling单元，将参数作如下修改：

在Script Generator主面板里，将Default output file设置为Graphene_eph-line.hdf5，脚本发送至Editor，在ElectronPhononCoupling一段，作如下修改： 

 electron_phonon_coupling = ElectronPhononCoupling(
     configuration=bulk_configuration,
     dynamical_matrix=dynamical_matrix,
     hamiltonian_derivatives=hamiltonian_derivatives,
     kpoints_cartesian=kpoints,
     qpoints_cartesian=qpoints,
     electron_bands=[3,4],
     phonon_modes=All,
     energy_tolerance=0.01*eV,                 //修改这一行 
     initial_state_energy_range=[-0.5,0.5]*eV, //修改这一行     store_dense_coupling_matrices=False,
     )

脚本发送至Job Manager，保存输入文件为Graphene_eph-line.py，完成计算。

由于计算只对k空间的一个线段采样，计算比(k,q)-全相关方法大大提速，这个计算在24核服务器上计算大约只需要5分钟。

迁移率

现在我们可以用两步计算来得到基于能量相关的弛豫时间$\tau(E)$的室温下的迁移率$\mu$。

单击Script Generator，添加以下计算单元：

打开Analysis from File，选择文件Graphene_relax.hdf5，加载Object Id为‘BulkConfiguration_1’的BulkConfiguration。

另外有三个计算单元自动添加：

在此例中，这三个数据都已经计算过了，因此可以将其删除，选择从相应的文件读取。添加三个‘Analysis from File’，按如下排序：

单击第二个‘Analysis from File’，选择Graphene_dynmat.hdf5，加载文件中的DynamicalMatrix数据。

单击第三个‘Analysis from File’，选择Graphene_dHdR.hdf5，加载文件中的HamiltonianDerivatives 数据。

单击第四个‘Analysis from File’，选择Graphene_eph-line.hdf5，加载文件中的 ElectronPhononCoupling 数据。

最后单击打开 Mobility，按如下图设置参数：

在Script Generator主窗口中，将Default output file设为 Graphene_mu_line_full.hdf5，脚本发送到Job Manager，保存输入文件为Graphene_mu_line_full.py，单击开始按钮完成计算。

计算结束后，就可以使用各向同性方法计算迁移率了。在本例中，我们使用前面计算得到的Graphene_relax.hdf5和Graphene_mu_line_full.hdf5中的bulk configuration 结果和 (k,q)-相关迁移率。

单击打开Script Generator，添加：

在 Script Generator 主面板中将 Default output file 设为 Graphene_mu-line-iso.hdf5，将脚本发送至 Editor。

在 Editor 里，修改脚本，读入全(k,q)-相关的迁移率数据，并修改mobility的计算： 

 # -------------------------------------------------------------
 # Analysis from File
 # -------------------------------------------------------------
 path = u'Graphene_relax.hdf5'
 configuration = nlread(path, object_id='BulkConfiguration_1')[0]
 
 # 增加这一段
 # -------------------------------------------------------------
 # Mobility
 # -------------------------------------------------------------
 mobility_full = nlread('Graphene_mu-line-full.hdf5', Mobility)[-1]
 
 # -------------------------------------------------------------
 # Mobility
 # -------------------------------------------------------------
 kpoint_grid = MonkhorstPackGrid(
     na=99,
     nb=99,
     )
 
 mobility = Mobility(
     configuration=configuration,
     method=Isotropic,
     energies=numpy.linspace(-0.24, 0.24, 100)*eV,
     inverse_relaxation_time=numpy.linspace(0, 1e+12, 100)*Second**-1,
     mobility_object=None,
     electron_bands=[3,4],    # 修改这一行
     kpoints=kpoint_grid,
     temperature=300*Kelvin,
     fermi_shift=0.13*eV,
     energy_broadening=0.003*eV,
     calculate_hall_coefficients=False,
     )
 nlsave('Graphene_mu-line-iso.hdf5', mobility)

将脚本发送至 Job Manager，保存 Graphene_mu-line-iso.hdf5，选择 Mobility，点击 Text Representation。 

 +------------------------------------------------------------------------------+
 | Mobility Report                                                              |
 | ---------------------------------------------------------------------------- |
 | Input parameters:                                                            |
 |     Temperature       = 300.00 K                                             |
 |     Fermi level shift = 0.13 eV                                              |
 |     Energy broadening = 0.0030 eV                                            |
 |     q-grid refinement = 1                                                    |
 |                                                                              |
 +------------------------------------------------------------------------------+
 | Trace of linear responce tensors:                                            |
 +------------------------------------------------------------------------------+
 |                                                                              |
 | Electrons:                                                                   |
 |                                                                              |
 |     Mobility                 =  2.58e+05 cm^2/(V*s)                          |
 |     Conductivity             =  1.79e+07 S/m                                 |
 |     Seebeck coefficient      =  7.91e-05 V/K                                 |
 |     Thermal conductivity     =  1.16e+02 W/(m*K)                             |
 |     Carrier density (2D, xy) =  8.63e+11 cm^-2                               |
 |                                                                              |
 |                                                                              |
 | Holes:                                                                       |
 |                                                                              |
 |     Mobility                 =  1.68e+06 cm^2/(V*s)                          |
 |     Conductivity             =  4.94e+04 S/m                                 |
 |     Seebeck coefficient      = -1.25e-03 V/K                                 |
 |     Thermal conductivity     =  5.98e+00 W/(m*K)                             |
 |     Carrier density (2D, xy) =  3.66e+08 cm^-2                               |
 |                                                                              |
 +------------------------------------------------------------------------------+

计算得到的迁移率为 $2.58 \cdot 10^{5}\ \mathrm{cm^{2} V^{-1} s^{-1}}$，与(k,q)-相关方法得到结果完全一致。

迁移率与温度的关系：(k,q)-相关方法和能量相关方法的结果

对能量相关方法更严格的测试是计算迁移率在室温范围对温度的依赖关系。

从下图可以看出，在 $100 \mathrm{K} \leq T \leq 300 \mathrm{K}$ 范围里，两种方法计算的迁移率对温度的依赖关系都比较好的再现了$1/T$特性。

理论知识

迁移率 $\mu$ 与电导率 $\sigma$ 有如下关系： $$ \mu = \frac{\sigma}{ne} \tag{1}\label{1}$$

$n$ 表示载流子浓度，$e$ 表示电子电荷。

本文中，我们将采用半经典玻尔兹曼输运方程（BTE）计算电导率。在零磁场和恒定电场下，定态玻尔兹曼方程简化为： $$\frac{q\mathbf{E}}{\hbar} \cdot \nabla_{\mathbf{k}} f_{\mathbf{k} n} \approx \frac{q\mathbf{E}}{\hbar} \cdot \nabla_{\mathbf{k}} f^0_{\mathbf{k} n} = \left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} \tag{2}\label{2}$$

方程($\ref{2}$)右侧，常被称为碰撞积分（collision integral）,用来描述不同的散射机制。由多种机制散射后，系统可耗散至稳态。方程左侧与电场强度线性相关，比例系数由分布函数决定。

电声子散射由碰撞积分描述，通常采用弛豫时间近似： $$\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} = \frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k}n}}$$

描述声学声子的准弹性散射（弛豫时间近似）。用 $\lambda$ 标记声子模式，$n$ 标记电子能带，$k$ 标记电子波矢。由 $|\mathbf{k}n\rangle$ 态到 $|\mathbf{k}'n'\rangle$ 态的跃迁几率可以由 费米黄金规则（FGR）得到：

$$\begin{split}P_{\mathbf{k}\mathbf{k'}}^{\lambda nn'} &= \frac{2\pi}{\hbar} |g_{\mathbf{k}\mathbf{k'}}^{\lambda n n'}|^2 [ n_{\mathbf{q}}^{\lambda} \delta \left(\epsilon_{\mathbf{k}'n'}-\epsilon_{\mathbf{k}n}-\hbar \omega_{q \lambda} \right) \delta_{\mathbf{k}',\mathbf{k}+\mathbf{q}} \\[.5cm] &+ (n_{\mathbf{-q}}^{\lambda}+1) \delta \left( \epsilon_{\mathbf{k}'n'}-\epsilon_{\mathbf{k} n}+\hbar \omega_{-q \lambda} \right) \delta_{\mathbf{k}',\mathbf{k}-\mathbf{q}} ],\end{split}$$

方括号中第一项和第二项分别表示吸收一个声子和放出一个声子。在 QuantumATK 中，使用 ElectronPhononCoupling 分析模块计算电声子耦合矩阵元 $ |g_{\mathbf{k}\mathbf{k'}}^{\lambda n n'}|$。从而得到电声子碰撞积分：

$$\begin{split}\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} &= - \sum_{\mathbf{k}'n'} \left[f_{\mathbf{k}n}\left(1-f_{\mathbf{k}'n'}\right)P_{\mathbf{k}\mathbf{k}'}^{nn'}-f_{\mathbf{k}'n'}\left(1-f_{\mathbf{k}n}\right)P_{\mathbf{k}'\mathbf{k}}^{n'n}\right]\end{split} \tag{3}\label{3}$$

弛豫时间近似

弛豫时间近似，是采用非常简单的形式来表示碰撞积分： $$\begin{split}\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} &=& -\frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k}n}}\end{split} \tag{4}\label{4}$$

其中 $\delta f_{\mathbf{k}n} = f_{\mathbf{k}n}-f^0_{\mathbf{k}n}$。因此，线性玻尔兹曼输运方程可写为： $$q \mathbf{E} \cdot \mathbf{v}_{\mathbf{k} n} \frac{\partial f^0_{\mathbf{k}n}}{\partial \epsilon_{\mathbf{k}n}} = -\frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k} n}}.$$

于是有： $$\delta f_{\mathbf{k}n} = -q \mathbf{E} \cdot \mathbf{v}_{\mathbf{k} n} \frac{\partial f^0_{\mathbf{k}n}}{\partial \epsilon_{\mathbf{k}n}} \tau_{\mathbf{k} n}.$$

弛豫时间可由（$\ref{3}$）式碰撞积分的一般形式计算得到。在准弹性散射近似下，($\omega_{q \lambda}\rightarrow 0$)，($\ref{4}$）式的特殊形式成立，（$\ref{3}$）式中的碰撞积分可简化为：

$$\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll}^{RTA} = \sum_{\mathbf{k}'n'} \left(f_{\mathbf{k}'n'}-f_{\mathbf{k}n}\right) P_{\mathbf{k}\mathbf{k'}}^{nn'} \equiv -\frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k} n}}$$,

因为在该近似下，$P_{\mathbf{k}\mathbf{k'}}^{nn'} = P_{\mathbf{k}'\mathbf{k}}^{n'n}$。

由此，我们得到散射几率（弛豫时间倒数）的表达式： $$\frac{1}{\tau_{\mathbf{k} n}} = \sum_{\mathbf{k}'n'} \left(1-\frac{\delta f_{\mathbf{k}'n'}}{\delta f_{\mathbf{k}n}}\right) P_{\mathbf{k}\mathbf{k'}}^{nn'} \tag{5}\label{5}$$

这里，应用了细致平衡方程 $P_{\mathbf{k}\mathbf{k'}}^{nn'}(f^0_{\mathbf{k}'n'}-f^0_{\mathbf{k}n})=0$，并假设散射各向同性。

为计算弛豫时间，我们引入散射角的概念 $$cos(\theta_{\mathbf{k}\mathbf{k'}}) = \frac{\mathbf{v}_{\mathbf{k}'n'} \cdot \mathbf{v}_{\mathbf{k}n}}{|\mathbf{v}_{\mathbf{k}'n'}| |\mathbf{v}_{\mathbf{k}n}|}$$

其中，$\mathbf{v}_{\mathbf{k}n}$ 表示电子速度。由文献【GMSB16】可得： $$\frac{1}{\tau_{\mathbf{k} n}} =\sum_{\mathbf{k}'n'} \frac{(1-f^0_{\mathbf{k}'n'})}{(1-f^0_{\mathbf{k}n})}\left(1-cos(\theta_{\mathbf{k}\mathbf{k'}})\right) P_{\mathbf{k}\mathbf{k'}}^{nn'}$$

迁移率

由弛豫时间，我们可以计算迁移率： $$\mu = -2q \frac{\sum_{\mathbf{k}n} v_{\mathbf{k} n}^2 \frac{\partial f^0_{\mathbf{k}n}}{\partial \epsilon_{\mathbf{k}n}} \tau_{\mathbf{k} n}}{\sum_{\mathbf{k}n}f^0_{\mathbf{k}n}} \tag{6}\label{6}$$

其中因子 2 表示自旋简并。

简化表达式和形变势

由上述公式推导可以看到，迁移率计算是一项大型数值计算任务，这也是采用多个近似的进行简化处理的原因。一个非常常用的近似是形变势近似。电声子耦合矩阵元用解析公式表示，从而直接计算上述公式的积分值。

在形变势近似中，电子-声子耦合项表示为： $$ | \langle \mathbf{k + q}|\delta V_\mathbf{q} | \mathbf{k} \rangle| = |M_{\mathbf{k},\mathbf{q}}| = D^{(1)}\cdot q $$

其中 $D^{(1)}$ 表示一阶形变势。一些声子模式（特别是光学模）零阶形变势已足够准确，即 $|M_{\mathbf{k},\mathbf{q}}|=D^{(0)}$ 与 $q$ 无关。

考虑仅有一个有效声速为$v_{eff}$的声学声子模，产生各向同性一阶形变势$D_{eff}$，迁移率在高温区的解析表达式可写为： $$\mu = \frac{4e\,\hbar\,\rho \,v_F^2 v_{eff}^2}{\pi\, n \, D_{eff}^2 k_B T} \tag{7}\label{7}$$

其中 $\rho$ 表示质量密度，$n$ 为载流子浓度。石墨烯中的 $n$ 与费米能量 $E_F$ 相关，$E_F = \alpha \sqrt{n}$。在高温区，迁移率与温度成反比。在低温极限下，迁移率服从 Bloch-Gruneisen 规律。在 Bloch-Gruneisen 温度 ($T_{BG}=2\hbar\,k_F\,v_s\,/k_B$) 以下，由于反散射 (back-scattering) 的相空间缩减，载流子寿命随温度降低而迅速增加，迁移率比 $T^{-1}$ 更快的增加。【文献KTJ12b】

参考文献