单层 TMD 材料的三能带紧束缚模型构建

引言

单层过渡金属二硫化物（TMD，如 $\mathrm{MoS_{2}}$、$\mathrm{WS_{2}}$、$\mathrm{MoSe_{2}}$等）是近年来非常受关注的一类二维半导体材料。他们在单层极限下表现为直接带隙体系，并且具有明显的自旋-轨道耦合效应和能谷（valley）自由度，这使得许多新奇的物理都集中在 $K/K'$ 点附近。在这个能量区间内，价带顶与导带底几乎完全由过渡金属原子的 $d$ 轨道组成，它们分别是 $d_{z^{2}}$、$d_{xy}$ 与 $d_{x^{2} - y^{2}}$。基于这一事实，我们可以构建一个仅包含这三个轨道的紧束缚模型，也就是常提到的三能带模型。尽管这个模型的基组很小，但已经足够描述 $K/K'$ 点附近的带边性质。

• 紧束缚模型： 所谓紧束缚模型，是指利用局域的原子轨道来描述晶体电子结构的办法
  1. 它假设电子主要局域在各个原子附近，即波函数由一组正交的原子轨道线性组合而成，并允许电子通过有限距离的跳跃在晶格中传播。
  2. 它的哈密顿量通常由原位（on-site）能级和 近邻/远近邻跳跃项组成，参数可以由对称性约束得到简化，并通过拟合 DFT 能带得到。
  3. 特点： 形式直观、计算高效、适合分析带边、能谷以及自旋-轨道耦合等物理性质。

TMD 的晶体结构与对称性分析

在构建紧束缚模型之前，我们首先要弄明白研究对象的晶体结构，这包含了材料的晶胞形状以及对称性信息。而紧束缚模型中跳跃矩阵的各项参数并非完全独立，受材料对称性的约束，这些参数之间有存在相应的关联，这些关联操作可以大大减少跳跃矩阵中独立参数的数量，使问题简单化。因此，我们有必要弄明白 TMD 材料的对称性信息。

• 注： 为了避免误会，接下来我们所提到的 TMD 材料均为单层材料。而 TMD 的单层材料中有两种常见的相：分别是 1T 相和 2H 相，其中 1T 相为金属， 2H 相为半导体（或绝缘体）。我们下面所讨论的均为 2H 相。

TMD 的实空间结构与布里渊区

单层 TMD 材料的化学式为 $\mathrm{MX_{2}}$ ，其中 $\mathrm{M}$ 为过渡金属元素（如 $\mathrm{Mo}$、$\mathrm{W}$ 和 $\mathrm{V}$ 等），$\mathrm{X}$ 为硫族元素（如 $\mathrm{S}$、$\mathrm{Se}$ 和 $\mathrm{Te}$ 等）。

其结构如图 1 所示，图中灰色原子为过渡金属原子，黄色原子为硫族原子，它们组成了一种类似“三明治”的结构，金属原子形成了二维的六角晶格，硫族原子分别位于金属原子的上下两侧，形成了三棱柱形的配位环境，形状如图 1 (b) 所示

而由于 $\mathrm{MX_{2}}$ 的原胞为六角晶型，他的布里渊区同样为六角形，如图 1 (c) 所示。在这个布里渊区中，除了原点 $\Gamma$ 和边中心 $M$ 以外，我们对位于六边形顶点的 $+K$ 与 $-K$ 点比较感兴趣，它们决定着能谷（valley）中存在许多有趣的物理性质。

TMD 的对称性

单层 TMD 的点群为 $D_{3h}$ ，它主要包括以下几种对称性：

1. 三重旋转对称性： 绕 $z$ 轴旋转 $120\degree$ 后系统保持不变，即

   $$\begin{align} C_{3}:\begin{bmatrix}x\\y\end{bmatrix}\longrightarrow\begin{bmatrix}\cos{(\theta)} & -\sin{(\theta)}\\ \sin{(\theta)} & \cos{(\theta)}\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix},\quad\theta = \frac{2\pi}{3} \end{align}$$

2. 水平镜面反射： 关于 $xy$ 平面反射后系统不变，即

   $$\begin{align} \sigma_{h}:(x,y,z)\longrightarrow(x,y,-z) \end{align}$$

3. 竖直的镜面反射： 分别沿 $\mathrm{M-X}$ 键做反映（以图 1 中沿 $yz$ 反映面为例）

   $$\begin{align} \sigma_{v}:(x,y,z)\longrightarrow(-x,y,z) \end{align}$$

以上三种对称操作可以相互组合，并形成更多的对称操作，下面将部分操作总结在表 1 中

表 1 ：TMD 材料的部分对称性操作

基组的确定及对称性算符

在构建紧束缚的哈密顿量之前，我们需要明确系统的基组。众所周知，“紧束缚模型”也称为“原子轨道线性组合法”（LCAO，linear combination of atomic orbitals），它使用原子轨道作为基组表示系统的波函数，因此我们需要明确具体使用哪些原子轨道。

以 $\mathrm{MoS_{2}}$ 为例，G.B.Liu 等人使用 DFT 软件计算得到它的原子轨道投影能带，如图 2 所示。不难看出 $\mathrm{MoS_{2}}$ 的价带最高（VBM）与导带最低（CBM）均主要由原子的 $d$ 轨道贡献，而 $p$ 轨道占比相比之下就小很多，$s$ 轨道更是几乎没有。由于我们通常只关心 VBM 和 CBM 的带隙与色散，因此只使用 $M$ 原子的这三个 $d$ 轨道来构成紧束缚模型的基组是可行的。为此我们认为 $\mathrm{X}$ 原子的所有轨道以及 $\mathrm{M}$ 原子的 $sp$ 轨道都可以忽略不计，他们不参与紧束缚模型的构建。

在此，我们的紧束缚模型基组确定为 $M$ 原子的三条 $d$ 轨道：

从轨道的基函数表达式不难判断，他们对应了几种对称性：

1. $|d_{z^{2}}\rangle$ 的基函数为 $z^{2}$ ，与 $xy$ 无关，因此它对于 $xy$ 平面内的操作不变，即 $C_{3}$ 操作与 $\sigma_{v}$ 操作对它不起作用：

   $$\begin{align} C_{3}:|d_{z^{2}}\rangle \longrightarrow |d_{z^{2}}\rangle;\quad \sigma_{v}:|d_{z^{2}}\rangle \longrightarrow |d_{z^{2}}\rangle \end{align}$$

2. $|d_{xy}\rangle$ 的基函数为 $xy$ ，它对于沿 $yz$ 平面的反映产生符号改变：

   $$\begin{align} \sigma_{v}:|d_{xy}\rangle \longrightarrow -|d_{z^{2}}\rangle \end{align}$$

3. $|d_{xy}\rangle$ 与 $|d_{x^{2} - y^{2}}\rangle$ 组合起来时，它对三重旋转操作有以下关系

   $$\begin{align} C_{3}:\begin{bmatrix}|d_{xy}\rangle\\|d_{x^{2} - y^{2}}\rangle\end{bmatrix} \longrightarrow \begin{bmatrix}\cos{(2\frac{2\pi}{3})} & -\sin{(2\frac{2\pi}{3})}\\ \sin{(2\frac{2\pi}{3})} & \cos{(2\frac{2\pi}{3})}\end{bmatrix} \begin{bmatrix}|d_{xy}\rangle \\ |d_{x^{2} - y^{2}}\rangle \end{bmatrix} \end{align}$$

   这个式子的意思是，当一个波函数可以用 $|d_{xy}\rangle$ 与 $|d_{x^{2} - y^{2}}\rangle$ 这两种轨道的线性组合表示时，它经过三重旋转后的结果仍可以使用这两个轨道表示，旋转操作只改变了两种轨道的组合系数。需要指出的是，尽管这个式子看起来与公式 $(1)$ 很像，但他们是完全不同的概念。当前对称操作作用的基底是原子轨道，而 $(1)$ 作用的基底则是位置坐标 $xy$ 。

那么，在 $\left\{|d_{z^{2}}\rangle,|d_{xy}\rangle,|d_{x^{2} - y^{2}}\rangle\right\}$ 基组中，上面两种对称操作可以用下面的算符表示：

有了以上两种对称操作后，表 1 中对称操作对应的矩阵形式就都可以得到了。

由于 $X$ 原子不参与到模型中，系统只由 $M$ 原子组成，并且所有 $M$ 原子都在同一平面内， $z$ 方向上的对称操作对这个系统的基函数也不起作用，因此 $D_{3h}$ 点群中其它与 $z$ 方向相关的操作，如 $xy$ 平面反映 $\sigma_{h}$ 与绕 $x$ 轴的二次旋转 $C_{2}$ 等均不需要考虑。

紧束缚哈密顿量的构建

坐标系与近邻原子

在此我们不妨明确一下系统的坐标系，如图 3 所示。

图中 $ab$ 表示两个晶格矢量，他们分别是

$l$ 是相邻两个 $M$ 原子的距离，在本文中它也表示晶格常数。最近邻原子、次近邻原子与第三近邻原子的位置分别用 $R_{n}$ 、$\tilde{R}_{n}$ 以及 $2R_{n}(n = 1,2,\cdots,6)$ 表示，它们对应的分数坐标为

紧束缚模型哈密顿量

由于我们研究的晶胞内只包含一个 $M$ 原子，每个原子有三个 $d$ 轨道，因此三能带模型的紧束缚哈密顿量在动量空间中表示为一个 $3\times3$ 的矩阵

为了方便起见，我们不妨定义以下记号

公式 $(10)$ 中每一项对应的表达式为

其中 $E_{ij}(\boldsymbol{R})$ 即为轨道间的跳跃项，

它表示位于原点的 $i$ 轨道跳跃到位于 $\boldsymbol{R}$ 的 $j$ 轨道的能量积分，这里 $\boldsymbol{R}$ 表示近邻的原子。

这个积分有一个十分重要的性质 $E_{ij}(\boldsymbol{R}) = E_{ji}(-\boldsymbol{R})^{*}$ ，证明如下

此外，当有了其中一个近邻原子轨道的 $E_{ij}(\boldsymbol{R})$ 后，我们可以通过对称性操作推导出其它近邻原子轨道对应的跳跃积分 $E_{ij}(\boldsymbol{R}')$ ，它们可以通过下面的式子推得

$E(R)$ 表示由 $E_{ij}(\boldsymbol{R})$ 组成的矩阵， $\hat{g}_{n}$ 是具体的对称操作， $D(\hat{g}_{n})$ 是 $\hat{g}_{n}$ 对应的矩阵形式，它们已在公式 $(8)$ 中给出。

在 TMD 中，根据考虑近邻原子的多少不同，跳跃项分最近邻项、次近邻项和第三近邻项，下面将以最近邻模型和第三近邻模型为例分别阐述它们的具体形式。

最近邻模型

现在我们考虑最简单的一种情况，即只考虑最近邻原子的跳跃，如图 4 所示

那么哈密顿量包含以下两项

其中 $H^{\mathrm{on-sote}}$ 是原位项，$H^{\mathrm{NN}}$ 则是最近邻跳跃项，它们的形式如下。

原位项

原位项是紧束缚模型里最简单的一项，我们把 $R = R_{0} = 0$ 代入公式 $(13)$ 后，可得

即各个原子轨道的本征值，它其实是比最近邻原子“更加近邻”的那一项。需要指出的是，在 $\Gamma$ 点，导带存在简并现象，因此原位项里 $\phi_{2}$ 与 $\phi_{3}$ 的本征值都是 $\epsilon_{2}$。

最近邻项

最近邻项包括 $R_{1},R_{2},\cdots,R_{6}$ 共 6 个位点，我们不妨设 $R_{0}$ 到 $R_{1}$ 的跳跃积分都为实数

其一共包含 9 个参数，但它们不完全是相互独立的，受对称性约束，有些参数之间存在关联，下面我们来寻找一下。

不难发现， $R_{1}$ 与 $R_{4}$ 刚好互为镜像，即 $R_{1} = -R_{4}$ ，可以应用公式 $(15)$ 以及公式 $(14)$ 得到：

使用下面的 SymPy 脚本，可以解出部分参数之间的关系

1
#!/usr/bin/env python3
2

3
import sympy as sym
4
from sympy import sqrt, S
5

6
t11, t12, t13 = sym.symbols("t11, t12, t13", real=True)
7
t21, t22, t23 = sym.symbols("t21, t22, t23", real=True)
8
t31, t32, t33 = sym.symbols("t31, t32, t33", real=True)
9

10
DSigmaV = sym.Matrix([[1,  0, 0],
11
                      [0, -1, 0],
12
                      [0,  0, 1]])
13

14
t = sym.Matrix([[t11, t12, t13],
15
                [t21, t22, t23],
16
                [t31, t32, t33]])
17

18
equation_matrix = DSigmaV * t * DSigmaV.T - t.T
19
equations = list(equation_matrix)
20
solutions = sym.solve(equations, (t11, t12, t13, t21, t22, t23, t31, t32, t33))
21
print(solutions)
22

23
# Output:
24
# {t12: -t21, t13: t31, t23: -t32}

将 $t_{12} = -t_{21},t_{13} = t_{31},t_{23} = -t_{32}$ 代入，可以得到简化后的跳跃矩阵：

现在我们只需要 $t_{11},t_{12},t_{13},t_{22},t_{23},t_{33}$ 这六个参数即可描述最近邻的跳跃了。

在确定了 $E(R_{1})$ 后，我们可以通过公式 $(15)$ 得到 $R_{2}$ 到 $R_{6}$ 对应的跳跃项

• $R_{2}$ 可以通过 $R_{1}$ 先转动 $\frac{2\pi}{3}$ ，再沿 $yz$ 面反映得到：

  $$\begin{align} E(R_{2}) = D(\sigma_{v}\cdot C_{3})E(R_{1})D(\sigma_{v}\cdot C_{3})^{\dag} \end{align}$$

• $R_{3}$ 可以通过 $R_{1}$ 直接转动 $\frac{2\pi}{3}$ 得到：

  $$\begin{align} E(R_{3}) = D(C_{3})E(R_{1})D(C_{3})^{\dag} \end{align}$$

• $R_{4}$ 刚好是 $-R_{1}$ ：

  $$\begin{align} E(R_{4}) = E(R_{1})^{\dag} \end{align}$$

• $R_{5}$ 是 $R_{1}$ 转动 $\frac{4\pi}{3}$ 得到：

  $$\begin{align} E(R_{5}) = D(C_{3}\cdot C_{3})E(R_{1})D(C_{3}\cdot C_{3})^{\dag} \end{align}$$

• $R_{6}$ 是 $R_{1}$ 转动 $\frac{4\pi}{3}$ 再沿 $yz$ 面反映得到：

  $$\begin{align} E(R_{6}) = D(\sigma_{v} \cdot C_{3} \cdot C_{3})E(R_{1})D(\sigma_{v} \cdot C_{3} \cdot C_{3})^{\dag} \end{align}$$

将以上式子代入公式 $(12)$ 可得到 $\boldsymbol{k}$ 点处完整的最近邻跳跃贡献

将参数 $\epsilon_{1},\cdots,\epsilon_{3},t_{11},\cdots,t_{33}$ 代入公式 $(16)$ 并对角化它，即可得到能带。

能带求解

以 $\mathrm{MoS_{2}}$ 为例，这里借用刘贵斌教授 Phys Rev B 给出的数值：

可得能带图如下：

不难看出，最近邻模型对带边的描述效果不错，比如对于 $K/K'$ 点的 VBM 和 CBM ，以及 $\Gamma$ 点的 VBM ，最近邻模型都给出的十分接近于 DFT 的能带。如果我们只关心 $K/K'$ 点的能带性质，那么最近邻模型已经是一个相当可用的工具。但如果我们想让最近邻模型的能带与 DFT 复合的更好，仅使用最近邻跳跃就不够了，我们需要考虑更多的跳跃项。

第三近邻模型

第三近邻模型在最近邻模型的基础上增加了次近邻以及第三近邻跳跃，这样可以得到更符合 DFT 的结果。他的哈密顿量包括

其中 $H^{\mathrm{NNN}}$ 是次近邻跳跃贡献， $H^{\mathrm{TNN}}$ 是第三近邻跳跃贡献。

第三近邻项

我们不妨从第三近邻跳跃开始，如图所示

第三近邻跳跃的矢量对应于 $2R_{n}$ ，它刚好是最近邻跳跃 $R_{n}$ 矢量的倍数，因此这种跳跃也符合最近邻跳跃的对称性约束，我们直接复用最近邻跳跃矩阵 $E(R_{1})$ 的形式即可。为方便区分， $E(2R_{1})$ 的矩阵元用 $u_{ij}$ 表示，即

因此第三近邻跳跃的哈密顿量形式如下：

次近邻项

不妨令次近邻跳跃 $\tilde{R}_{1}$ 对应的积分为 $E(\tilde{R}_{1})$ ，有

次近邻跳跃如图所示，此时跳跃矢量不再是 $R_{n}$ 的倍数，它的方向有所改变，因此我们需要寻找其他的对称性来减少参数数量。

仔细观察，不难看出，对于 $\tilde{R}_{1}$ ，我们可以先做 $C_{3}$ 转动变换为 $\tilde{R}_{3}$ ，再经过 $yz$ 面反映变换回 $\tilde{R}_{1}$ 。即

同样，使用下面的 SymPy 脚本，可以解出 $E(\tilde{R}_{1}) 中部分参数之间的关系

1
r11, r12, r13 = sym.symbols("r11, r12, r13", real=True)
2
r21, r22, r23 = sym.symbols("r21, r22, r23", real=True)
3
r31, r32, r33 = sym.symbols("r31, r32, r33", real=True)
4

5
DC3 = sym.Matrix([[1,          0,         0],
6
                  [0,    -S(1)/2, sqrt(3)/2],
7
                  [0, -sqrt(3)/2,   -S(1)/2]])
8

9
r = sym.Matrix([[r11, r12, r13],
10
                [r21, r22, r23],
11
                [r31, r32, r33]])
12

13
equation_matrix = (DSigmaV * DC3) * r * (DSigmaV * DC3).t - r
14
equations = list(equation_matrix)
15
solutions = sym.solve(equations, (r11, r12, r13, r21, r22, r23, r31, r32, r33))
16
print(solutions)
17

18
# Output:
19
# {r12: -sqrt(3)*r13, r21: -sqrt(3)*r31, r22: -2*sqrt(3)*r32/3 + r33, r23: r32}

因此可得简化后的跳跃矩阵

有了 $\tilde{R}_{1}$ 后，我们可以通过对称性陆续得到其它的次近邻跳跃项 $\tilde{R}_{2},\cdots,\tilde{R}_{6}$：

• 为方便理解，我们先从 $\tilde{R}_{4}$ 开始，刚好有 $\tilde{R}_{4} = -\tilde{R}_{1}$ ，因此有

  $$\begin{align} E(\tilde{R}_{4}) = E(\tilde{R}_{1})^{\dag} \end{align}$$

• $\tilde{R}_{2}$ 可以使用 $\tilde{R}_{4}$ 经过两次 $C_{3}$ 转动得到

  $$\begin{align} E(\tilde{R}_{2}) = D(C_{3} \cdot C_{3})E(\tilde{R}_{1})^{\dag}D(C_{3} \cdot C_{3})^{\dag} \end{align}$$

• $\tilde{R}_{3}$ 可以使用 $\tilde{R}_{1}$ 经过一次 $C_{3}$ 转动得到

  $$\begin{align} E(\tilde{R}_{3}) = D(C_{3})E(\tilde{R}_{1})D(C_{3})^{\dag} \end{align}$$

• $\tilde{R}_{5}$ 可以使用 $\tilde{R}_{1}$ 经过两次 $C_{3}$ 转动得到

  $$\begin{align} E(\tilde{R}_{5}) = D(C_{3} \cdot C_{3})E(\tilde{R}_{1})D(C_{3} \cdot C_{3})^{\dag} \end{align}$$

• $\tilde{R}_{6}$ 可以使用 $\tilde{R}_{4}$ 经过一次 $C_{3}$ 转动得到

  $$\begin{align} E(\tilde{R}_{6}) = D(C_{3})E(\tilde{R}_{1})^{\dag}D(C_{3})^{\dag} \end{align}$$

此时次近邻跳跃的哈密顿量形式如下：

将参数代入 $H^{\mathrm{NNN}}$ 和 $H^{\mathrm{TNN}}$ 中，即可得到第三近邻模型下的能带结构。需要注意的是，在使用第三近邻模型时， $H^{\mathrm{on-site}}$ 与 $H^{\mathrm{NN}}$ 的参数也需要更新。

能带求解

这里同样使用刘贵斌教授给出的数值：

可得能带图如下：

显然，这次模型所得能带与 DFT 符合的很好，这说明考虑次近邻跳跃与第三近邻跳跃后的三能带模型已经可以很好地满足带边的描述需求。

自旋-轨道耦合的引入

在前面的讨论中，我们构建了一个仅包含轨道自由度的紧束缚模型。然而，对于单层 TMD 材料而言，自旋-轨道耦合（spin-orbit coupling，SOC）不仅仅是一种微小的能量修正，它是决定能带精细结构以及谷激子性质的关键因素，因此我们有必要研究如何将 SOC 引入紧束缚模型。

原子内 SOC 的基本形式

所谓自旋轨道耦合，就是电子的自旋角动量与轨道角动量相互作用的一种狭义相对论效应。当电子在原子核产生的电场中高速运动时，相对论效应会使其感应出有效磁场，这种磁场与电子自身的自旋磁矩产生相互作用，从而导致能级分裂和谱线精细结构的产生。

SOC 的通用形式为

其中：

• $\sigma$ 是泡利矩阵；
• $\boldsymbol{p}$ 是电子的动量；
• $v(\boldsymbol{r})$ 是电子感受到的静电势；

在真实体系中，这个表达中的 $\nabla v(\boldsymbol{r})$ 并不容易得到，因此我们将它与前面的系数合并为 $\lambda$ 作为一个经验参数给出，此时 SOC 可以写成

其中：

• $\lambda = \frac{\nabla v(\boldsymbol{r})}{2mc^{2}}$ 表示 SOC 作用的强度，它决定了 SOC 对能级劈裂的贡献；
• $\boldsymbol{L}$ 是电子的轨道角动量；
• $\boldsymbol{S} = \frac{\hbar}{2}\sigma$ 是电子的自旋角动量。

我们将用这个式子对紧束缚模型中的 SOC 作用进行建模。需要指出的是，一般认为 SOC 仅作用于 on-site 项，即它只作用于同一个原胞内的同一个原子，不涉及跳跃项。

三能带模型中的 SOC 作用

在三能带模型中，我们选取的轨道基组为：

当 $\boldsymbol{L}$ 算符作用到它们时， $L_{x}$ 与 $L_{y}$ 分量所得结果可能会得到 $|d_{yz}\rangle$ 、 $|d_{xz}\rangle$ 分量，而这些分量并不在三能带模型的基组内；相反，当使用 $|L_{z}\rangle$ 算符作用时，所得结果仍可用现有的基组表达。因此在三能带近似下，可以将 $m_{z}$ 视为一个号量子数， SOC 被简化为

下面我们来讨论三能带基组下 SOC 的矩阵元。

$L_{z}$ 算符

对于各个实轨道，它们都可以写成复球谐函数的线性组合：

• $|d_{z^{2}}\rangle = Y_{2,0}$
• $|d_{xy}\rangle = \frac{\mathrm{i}}{\sqrt{2}}(Y_{2,-2} - Y_{2,2})$
• $|d_{x^{2} - y^{2}}\rangle = \frac{1}{\sqrt{2}}(Y_{2,-2} + Y_{2,2})$

对每个球谐函数，有

将 $L_{z}$ 作用于实轨道时，有

• $\hat{L}_{z}|d_{z^{2}}\rangle = 0$
• $\hat{L}_{z}|d_{xy}\rangle = \frac{i}{\sqrt{2}}(-2Y_{2,-2} - 2Y_{2,2}) = -2\mathrm{i}\frac{1}{\sqrt{2}}(Y_{2,-2} + Y_{2,2}) = -2\mathrm{i}|d_{x^{2} - y^{2}}\rangle$
• $\hat{L}_{z}|d_{x^{2} - y^{2}}\rangle = \frac{1}{\sqrt{2}}(-2Y_{2,-2} + 2Y_{2,2}) = 2\mathrm{i}\frac{-\mathrm{i}}{\sqrt{2}}(-Y_{2,-2} + Y_{2,2}) = 2\mathrm{i}|d_{xy}\rangle$

因此， $\hat{L}_{z}$ 的矩阵元 $\langle i|\hat{L}_{z}|j\rangle$ 中有一下关系

• $\langle d_{z^{2}}|\hat{L}_{z}|j\rangle = 0$
• $\langle d_{xy}|\hat{L}_{z}|d_{x^{2} - y^{2}}\rangle = \langle d_{x^{2} - y^{2}}|\hat{L}_{z}|d_{xy}\rangle^{*} = -2\mathrm{i}$
• $\langle d_{xy}|\hat{L}_{z}|d_{xy}\rangle = \langle d_{x^{2} - y^{2}}|\hat{L}_{z}|d_{x^{2} - y^{2}}\rangle = 0$

把上面的关系写成矩阵形式有

$S_{z}$ 算符

由于我们忽略了 $L_{x}$ 与 $L_{y}$ 的贡献，因此电子自旋也只需考虑 $S_{z}$ 的贡献

此时我们可以得到完整的 SOC 矩阵：

这里我们相当于引入了新的自由度，即自旋，因此需要使用克罗内积（Kronecker Product，也称“直积”）符号 $\otimes$ 。

考虑自旋-轨道耦合后的能带

前面我们已经得到了紧束缚模型的哈密顿量 $H^{0}$ 和完整的 SOC 矩阵 $H^{\mathrm{SOC}}$ ，现在可以将它们组合起来

这里 $I_{2} =\begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}$ 表示单位矩阵。

借用刘贵斌教授 Phys. Rev. B 给出的数值， $\mathrm{MoS_{2}}$ 中 $\lambda = 0.073$ ，代入 $H^{\mathrm{TB + SOC}}$ 中，可得能带如图 5 所示

不难看出，在包含 SOC 后，三能带模型与 DFT 计算结果仍然符合地不错。比如 $K$ 谷里 VBM 较大的劈裂和 CBM 微小的劈裂都被良好地复现出来，而 $\Gamma$ 谷 VBM 存在简并的特性也一样得到保留。这说明在只取 $L_{z}$ 近似下的 SOC 就可以描述出较为合理的自旋-轨道劈裂。

参数的拟合

到目前为止，我们的模型在建立好之后一直依赖查文献来获取模型的参数，这对于 $\mathrm{MoS_{2}}$ 这样已经十分熟悉的体系而言是没什么问题的，但对于一个文献鲜有报道的体系，我们就不得不自己来拟合这些参数。

我们来简化一下这个问题，它等同于：

有一个函数 $f$ ，它接受一个参数集 $\left\{p\right\}$ 并输出结果 $E$ ，现在求一个 $\left\{p\right\}$ 使得 $E$ 尽可能接近目标值 $E_{0}$ ，即

$$\begin{align} \min_{p}{\sum{\left[f(\left\{p\right\}) - E_{0}\right]^{2}}} \end{align}$$

• 其中 $f$ 即为紧束缚模型；
• $\left\{p\right\} = \left\{\epsilon_{1},\epsilon_{2},t_{11},\cdots\right\}$；
• $E_{0}$ 为目标能带的具体数值。

这是一个典型的最优化问题，直接交给程序去完成即可。

目标能带

这里我们以 $\mathrm{MoS_{2}}$ 为例演示如何通过拟合 HSE 能带（不带 SOC ）得到三能带紧束缚模型的参数值。图 6 中突出显示的为目标能带。

构造损失函数

损失函数的构造其实相当简单，只需

其中 $E[p](n, \boldsymbol{k})$ 是在当前参数空间 $p$ 中计算所得到的能带值。这个损失函数其实就表示当前能带“距离目标能带有多远”，式中 $|\quad|^{2}$ 就表示求两个量的欧氏距离。

程序设计和计算结果

这里我们使用 Julia 的 Optim.jl 包，大概流程如下

1
# Evaluate bands, pvec is the vector of TB params
2
function forward_bands(pvec, kpath, model, ...)
3
    # Const and diagonalize the Hamiltonian
4
    Hk = ...
5
    E = eigvals!(Hk)
6
    return E
7
end
8

9
# Loss function
10
function make_loss(Etarget, kpath, model, cache)
11
    function loss(pvec)
12
        Etb = forward_bands(pvec, kpath, model, cache)
13
        s = 0.0
14
        @inbounds for b in 1:3, i in 1:size(Etb,2)
15
            d = Etb[b,i] - Etarget[b,i]
16
            s += w * d*d
17
        end
18
        return s
19
    end
20
    return loss
21
end
22

23
# main function
24
loss = make_loss(Etarget, kpath, model, cache; mask=mask, p0=p0, reg=reg, nsegs=nsegs)
25

26
# 1) Coarse: Nelder-Mead (robust, no gradients)
27
res1 = optimize(loss, p0, NelderMead())
28
p1 = Optim.minimizer(res1)
29
@info "After NelderMead: loss = $(Optim.minimum(res1))"
30

31
# 2) Refine: BFGS (fast local)
32
res2 = optimize(loss, p1, BFGS(), Optim.Options(show_trace=true))
33
p_opt = Optim.minimizer(res2)
34
@info "After BFGS: loss = $(Optim.minimum(res2))"