交换项和关联项均会导致额外的能量下降。其中，由于 Hartree–Fock 方法显式构造了 Slater 行列式，因此其计算的交换作用是精确的；主要难以准确描述的是关联项。

上述为一次量子化。由于需要构造巨大的 Slater 行列式，且无法处理粒子数变化的情况，因此需要进行二次量子化来处理这种情况。

二次量子化

在进入二次量子化之前，我们需要引入一些工具来描述二次量子化。

产生与湮灭算符

由于我们需要研究激子，而激子不是不变的——它们可以通过吸收光子而产生，也能通过辐射（发光）或非辐射（声子）过程湮灭。产生和湮灭算符可以很好地描述这种动态的粒子数变化情况。

位置算符与动量算符

在二次量子化的框架下，产生与湮灭算符并不是基础的算符，而是由位置算符和动量算符派生出来的。因此我们先来看位置算符和动量算符是怎么得来的。（以下过程可能和格里菲斯量子力学讲义中略有不同）

对于一个处于 $x$ 位置的粒子[^c3]，其位置期望值应该是：

我们可以看出，一个"算符"是对其后函数做某种操作的指令。位置算符告诉你乘以 $x$，作用在波函数上，就能得到波函数的位置期望值。因此，我们将位置算符定义为 $\hat{x}$。

同样地，我们可能对粒子的运动速度感兴趣，对 $\langle x \rangle$ 关于时间 $t$ 求导可得[^c4]：

注意： 这里没有对 $x$ 求 $dt$ 是因为 $x$ 的积分上下限是负无穷到正无穷，$x$ 为积分变量而与 $t$ 无关，即 $\frac{\partial x}{\partial t} = 0$，因此求导仅对波函数求导。

使用含时薛定谔方程替换 $\frac{\partial \Psi}{\partial t}$，可以得到：

仅在左侧保留 $\frac{\partial \Psi}{\partial t}$ 可得：

其复共轭可以写作：

将式 (9) 和式 (10) 代入式 (7’)，可以消去势能项，仅保留动能项：

提取其中公因子 $-\frac{i\hbar}{2m}$：

对其中 $\int \frac{\partial^2 \Psi^*}{\partial x^2} x \Psi , dx$ 进行分部积分计算，就可以得到：

一般我们使用动量来描述粒子的速度，我们将其转化为动量：

上式又能写成：

可以看出，括号内就是描述粒子动量的算符。它的意思是：对 $x$ 求导并将结果乘以 $-i\hbar$，作用在波函数上，就能得到波函数的动量期望值。因此，我们将动量算符定义为：

我们对位置算符和动量算符还有一个重要的推广，即二者的对易关系为 $[\hat{x}, \hat{p}_x] = i\hbar$，其中：

分别计算两项：

至此，我们完成了对位置算符和动量算符的推导。接下来我们需要在此基础上推导产生和湮灭算符，而其起源则是一维谐振子。

一维谐振子

对于一个一维谐振子（例如弹簧上的小球）来说，它的经典能量是动能和势能之和（在经典力学中，哈密顿量为动能和势能的总和）：

其中，$p$ 为粒子动量，$m$ 为粒子质量，$x$ 为粒子位置，$\omega$ 为角频率。

在一次量子化中，我们把 $p$ 和 $x$ 变成算符 $\hat{p}$ 和 $\hat{x}$，满足 $[\hat{x}, \hat{p}] = i\hbar$。此时哈密顿算符可以写作：

我们求解的思路是将哈密顿算符分解。如果是普通代数量 $u$ 和 $v$，比较容易看出 $u^2 + v^2 = (u - iv)(u + iv)$，但是由于 $\hat{p}$ 和 $\hat{x}$ 是算符，因此不能随意交换顺序。我们先引入约化后的过渡算符，去掉常数的干扰：

此时哈密顿量可以写成更简洁的形式：

这时，我们可以借用平方和分解的思想，将 $u^2 + v^2$ 分解为 $(u - iv)(u + iv)$。表现在算符上可以写作：

注意：

• 我们将两个算符定义为 $(\hat{X} + i\hat{P})$ 和 $(\hat{X} - i\hat{P})$。
• 这二者是不能交换顺序的，不能写作 $(\hat{P} + i\hat{X})$，因为它们依然是非对易的。
• 证明如下：
• 考虑一个连在弹簧上的小球（质量 $m$，频率为 $\omega$）。
• 根据牛顿定律：$F = ma \implies -kx = m \ddot{x}$。我们更喜欢用频率代替 $k$，即 $k = m\omega^2$。
• 因此上式可以写作：$\ddot{x} + \omega^2 x = 0$。
• 这个方程的通解是正弦或者余弦函数，粒子的位置随时间变化的解为：$x(t) = A \cos(\omega t)$。
• 此处的 $A$ 为振幅。
• 而动量则为 $p = mv = m\dot{x}$，需要对时间求导，即：

$$\dot{x}(t) = \frac{d}{dt} [A \cos(\omega t)] = -A\omega \sin(\omega t) \tag{23}$$

• 因此，动量表示为：

$$p(t) = m \cdot (-A\omega \sin(\omega t)) = -m\omega A \sin(\omega t) \tag{24}$$

• 如果我们将位置和动量画在同一个坐标系中，横轴是位置，纵轴是动量，我们可以得到：

$$\frac{x(t)}{A} = \cos(\omega t) \tag{25}$$$$\frac{p(t)}{-m\omega A} = \sin(\omega t) \tag{26}$$

• 它们的平方和为 1，即：

$$\left( \frac{x}{A} \right)^2 + \left( \frac{p}{m\omega A} \right)^2 = 1 \tag{27}$$

• 它代表了：在 $x$-$p$ 相平面上，速度最大的时候位置在 0 点；位置最大（$A$）的时候速度为 0。
• 如果我们使用欧拉公式描述，即 $e^{-i\omega t} = \cos(\omega t) - i\sin(\omega t)$，那么实部对应位置，虚部对应动量，组合即 $\hat{X} + i\hat{P}$。它代表了：在复平面上，位置和动量有 90° 的相位差，即动量和位置是不同步的。
• 如果我们将算符写成 $\hat{P} + i\hat{X}$（这里先借用一个性质：湮灭算符作用在基态上必须等于 0，即 $\hat{a} |0\rangle = 0$），那么作用在波函数上，我们就得到 $\left( -i\frac{d}{dX} + iX \right)\psi = 0$，对应的波函数为 $\psi \sim e^{+X^2/2}$，这是一个发散的函数，归一化积分为无穷大，物理上不存在。因此，我们定义的算符顺序必须是 $(\hat{X} + i\hat{P})$ 和 $(\hat{X} - i\hat{P})$。

产生湮灭算符的性质及物理意义

其中，式 (21b) 定义为湮灭算符 $\hat{a}-$，式 (22) 定义为产生算符 $\hat{a}+$。

接下来我们考虑它们的对易子 $[\hat{a}-, \hat{a}+]$：

其中，为简化计算令 $C = \sqrt{\frac{m\omega}{2\hbar}}$。

对方括号中每一项求对易，可得：

根据算符关系，算符与自身对易，即 $[\hat{x}, \hat{x}] = 0$，因此第一项和第四项为 0。第二项 $-\frac{i}{m\omega}[\hat{x}, \hat{p}]$ 由对易关系 $[\hat{x},\hat{p}]=i\hbar$ 得到 $\frac{\hbar}{m\omega}$，同理第三项也为 $\frac{\hbar}{m\omega}$。

因此，最终可以得到：

因此，式 (21a) 的哈密顿量最终可以分解成[^c5]：

注意到由于算符的次序问题，导致后面是 $+\frac{1}{2}$。利用式 (30) 的对易关系 $[\hat{a}-, \hat{a}+] = 1$，更一般的情况可以写作：

或者：

式 (32) 和式 (33) 是等价的。

下面我们说明为什么式 (21b) 被叫做湮灭算符、式 (22) 被叫做产生算符。

从薛定谔方程 $\hat H \psi = E \psi$ 出发，如果将上升算符 $\hat{a}_+$ 作用在波函数上：

其中对易关系为 $[\hat{H}, \hat{a}+] = \hbar\omega , \hat{a}+$，则：

这说明，$\hat{a}_+$ 作用在波函数上后，新的波函数依然是 $\hat{H}$ 的本征态，本征能量为 $E + \hbar\omega$。这意味着，我们可以使用上升/下降算符（产生/湮灭算符）直接生成解。

然而如果对基态 $\psi_0$（能量 $\frac{1}{2}\hbar\omega$）作用湮灭算符：$\hat{a}_- \psi_0 = 0$，它代表波函数消失。这保证谐振子存在一个最低能量（基态）。

注意： 我们用 $\hat{x}, \hat{p}$ 表示出了产生湮灭算符，反过来也可以用产生湮灭算符表示位置算符和动量算符：

$$\hat{x} = \sqrt{\frac{\hbar}{2m\omega}} (\hat{a}_+ + \hat{a}_-) \tag{36}$$$$\hat{p} = i\sqrt{\frac{m\hbar\omega}{2}} (\hat{a}_+ - \hat{a}_-) \tag{37}$$

这时，$\hat{x}$ 代表着：测量粒子位置时，对粒子数状态进行了扰动，系统同时处于"少一个粒子"和"多一个粒子"的叠加态。同理动量算符也是如此。将位置算符作用在有 $n$ 个粒子的 $|n\rangle$ 态上：

$$\begin{aligned} \hat{x} |n\rangle &\propto (\hat{a}_- + \hat{a}_+) |n\rangle \\ &= \hat{a}_- |n\rangle + \hat{a}_+ |n\rangle \\ &= \sqrt{n} |n-1\rangle + \sqrt{n+1} |n+1\rangle \end{aligned} \tag{38}$$

可以看到，位置算符作用在 $|n\rangle$ 态后，系统处在 $|n-1\rangle$ 和 $|n+1\rangle$ 的叠加态。

注意：

• 产生算符 $\hat{c}^\dagger$ 和湮灭算符 $\hat{c}$ 互为厄米共轭。（在二次量子化中，我们使用厄米共轭符号 $\dagger$ 表示产生算符，不加 $\dagger$ 的代表湮灭算符，作为一般的记号规则。）
• 湮灭算符我们使用之前的无量纲化形式：

$$\hat{c} = C \left( \hat{x} + \frac{i}{m\omega} \hat{p} \right) \tag{39}$$

• 其中 $C = \sqrt{\frac{m\omega}{2\hbar}}$。
• 对其取厄米共轭 $\dagger$：

$$\hat{c}^\dagger = \left[ C \left( \hat{x} + \frac{i}{m\omega} \hat{p} \right) \right]^\dagger \tag{40}$$

• 其中实数 $C$ 保持不变，即 $C^*=C$。又因为力学量算符均为厄米算符，故 $\hat{x}^\dagger = \hat{x}$、$\hat{p}^\dagger = \hat{p}$，中间仅 $i$ 变号，整理得：

$$\hat{c}^\dagger = C \left( \hat{x} - \frac{i}{m\omega} \hat{p} \right) \tag{41}$$

• 我们就得到了产生算符，证毕。
• 因此，我们在后续约定使用 $\hat{c}^\dagger$ 代表产生算符，$\hat{c}$ 代表湮灭算符。

向量计算

我们考虑一个包含三个态的系统，把它写作列向量：第一行定义为真空态，第二行定义为态 1，第三行[^c6]定义为态 2，并作如下约定：

它代表 100% 概率是空的，0% 在轨道（态）1，0% 在轨道（态）2。

它代表电子 100% 在第一个轨道（态）里。

它代表电子 100% 在第二个轨道（态）里。

产生算符的含义是将真空态变为某个态中有电子的情况。以态 1 为例，要把 $\begin{pmatrix} 1 \ 0 \ 0 \end{pmatrix}$ 变成 $\begin{pmatrix} 0 \ 1 \ 0 \end{pmatrix}$，那么态 1 的产生算符应该写作：

将这个算符作用在真空态：

这时电子就产生在态 1 上。同理，态 2 的产生算符和湮灭算符也是如此，这里不再赘述。

如果要将电子从态 1 跳跃到态 2，我们可以写作：

因为算符的作用顺序是从右向左：首先湮灭算符作用在态 1 消灭电子，随后产生算符在态 2 产生电子，它们整体可以写作跳跃算符 $\hat{t}_{1\to2}$。

场算符

前面我们推导得到的哈密顿量（即式 33）是从非相互作用的谐振子得到的，也就是最理想的情况。然而实际更复杂的量子场是由无数个耦合在一起的谐振子组成的，我们并不清楚实际的产生湮灭算符是什么样子的。

现在，我们需要从一次量子化的形式转换为二次量子化。

在一次量子化中，希尔伯特空间内任意一个量子态 $|\psi\rangle$ 都可以用一组完备的基底 ${|n\rangle}$ 展开：

其中 $C_n$ 是投影系数。

注意：

• 如果我们使用哈密顿量作用在态 $|\Psi\rangle$ 上，得到的本征方程为：

$$\hat{H} |\Psi\rangle = E |\Psi\rangle \tag{49}$$

• 它代表了哈密顿算符 $\hat{H}$ 作用在波函数 $|\Psi\rangle$ 上，得到了一个能量值 $E$。表明 $|\Psi\rangle$ 是一个能量定态：使用哈密顿算符测量能量，会得到一个 100% 的确定值 $E$，且测量后系统状态不变。上式为能量的本征方程。
• 类似地，如果有一个位置算符 $\hat{\vec{r}}$ 作用在状态 $|\Psi\rangle$ 上，得到的本征方程为：

$$\hat{\vec{r}} |\Psi\rangle = \vec{R} |\Psi\rangle \tag{50}$$

• 它表明 $|\Psi\rangle$ 是一个位置定态。如果测量位置，那么粒子 100% 处于位置 $\vec{R}$ 处。上式为位置的本征方程。
• 重要： 根据波函数的定义——处于状态 $|n\rangle$ 的粒子在位置 $\vec{r}$ 被发现的概率幅——可以将波函数 $\phi_n(\vec{r})$ 写作：

$$\phi_n(\vec{r}) = \langle \vec{r} | n\rangle \tag{51}$$

• 其中，只要 ${|n\rangle}$ 是完备的，就具备以下性质：

$$\sum_n |n\rangle \langle n| = \hat{1} \tag{52}$$

• 这个算符为单位算符，可以将任何态投影上去且不丢失任何信息。上式叫做完备性关系。

现在我们需要了解位置态。使用 ${|n\rangle}$ 将 $|\vec{r}\rangle$ 展开：

其中，$\langle n | \vec{r}\rangle$ 是位置态在 $|n\rangle$ 基底下的展开系数。由波函数的定义 $\phi_n(\vec{r}) = \langle \vec{r} | n\rangle$，对它求复共轭可得：

因此位置态 $|\vec{r}\rangle$ 可以用能级态 $|n\rangle$ 表示：

在二次量子化中，每个态都对应一个产生算符。根据之前的推导，产生算符 $\hat{c}_n^\dagger$ 作用在真空上，会创造出状态 $|n\rangle$：

同样地，如果有一个产生算符 $\hat{\psi}^\dagger(\vec{r})$，作用在真空上，能够创造出位置在 $\vec{r}$ 的粒子：

我们将式 (56) 和 (57) 代入式 (55)，就可以得到：

由于两个算符作用在真空态上得到相同的态，因此我们可以将两边的真空态去掉，就得到了：

上面这个算符我们将其定义为场产生算符。根据之前的推导，产生算符和湮灭算符互为厄米共轭，因此将上式整体取厄米共轭就能得到场湮灭算符：

场产生算符的物理意义是：在位置 $\vec{r}$ 处创造一个以 $\hat{c}_n^\dagger$ 为基底、贡献为 $\phi_n^*(\vec{r})$ 的粒子。同理，场湮灭算符则是消灭这样的粒子。

注意：

• 在式 (58b) 中，不要将左边的 $\hat{\psi}^\dagger(\vec{r})$ 跟波函数 $\psi$ 搞混，它是一个算符，是用波函数作为系数展开的。
• 右边的 $\hat{c}_n^\dagger$ 物理含义是在第 $n$ 个量子态（轨道）里制造一个粒子。
• $\phi_n^*(\vec{r})$ 代表了第 $n$ 个轨道在 $\vec{r}$ 处的权重，即 $\langle n | \vec{r}\rangle$。
• 例如，在 $r = 0$ 处制造一个粒子，而现在有两个波：
  • 波 1（$\hat{c}_1^\dagger$）：$\cos(x)$，在 $x=0$ 处值为 1。
  • 波 2（$\hat{c}_2^\dagger$）：$\sin(x)$，在 $x=0$ 处值为 0。
• 我们使用场产生算符：$\hat{\psi}^\dagger(0) = \sum_n \phi_n^*(0) \hat{c}_n^\dagger$
  • 波 1：$\phi_1^*(0) = \cos(0) = 1$
  • 波 2：$\phi_2^*(0) = \sin(0) = 0$
• 最后：$\hat{\psi}^\dagger(0) = 1 \cdot \hat{c}_1^\dagger + 0 \cdot \hat{c}_2^\dagger$
• 它代表了：如果要在 $r = 0$ 处创造一个粒子，余弦波占比是 1，而正弦波则不贡献。

场算符的性质及物理意义

1. 反对易关系

场算符（费米子）满足如下反对易关系（${A, B} \equiv AB + BA$）：

它说明了：在同一个位置发生产生湮灭操作会"撞上"（出现 $\delta$ 函数贡献），如果位置不同则互不影响。

两个湮灭算符反对易[^c7]——即 $\hat{\psi}(\vec{r})\hat{\psi}(\vec{r}’) = -\hat{\psi}(\vec{r}’)\hat{\psi}(\vec{r})$。它体现了费米子波函数对粒子交换的反对称性：交换两次湮灭操作的次序会带来一个负号，但物理观测量（密度、能量等）不变。

两个产生算符也满足反对易关系。特别地，取 $\vec r=\vec r’$ 可得 $[\hat\psi^\dagger(\vec r)]^2 = 0$，即不能在同一个位置产生两个相同自旋的电子——这正是泡利不相容原理的体现。

以下我们做式 (59b) 的推导：

我们将求和符号提出：

其中我们有以下基本反对易关系（公理）：

因此，式 (63) 末尾的反对易子变为 $\delta_{nm}$，消掉一个求和符号：

由于基底 ${|n\rangle}$ 是完备的，因此：

证毕。

2. 电子密度算符

我们在 DFT 中关心的电子密度，在一次量子化中可以写成 $\rho = |\psi|^2 = \psi^* \psi$。那么在二次量子化中，可以将其表示为：

上式被称为**（局域）密度算符**或粒子数密度算符，其物理意义是 $\vec{r}$ 处的电子密度，与一次量子化相比只是变成了算符形式。

上式较为重要，且其产生湮灭算符顺序不能更换。对于费米子，由于泡利不相容原理，在位置 $\vec{r}$ 处（给定自旋下）只能存在一个粒子，因此使用粒子数算符只能得到两个结果：

• 原本在 $\vec{r}$ 处的态是 $|0\rangle$（没有粒子），湮灭算符先作用得到 0，说明这里粒子数为 0，无需还原。
• 原本在 $\vec{r}$ 处的态是 $|1\rangle$（存在一个粒子），湮灭算符先作用得到 $|0\rangle$，这时产生算符再作用上去，把态还原成 $|1\rangle$，最终知道这里粒子数为 1。

3. 总粒子数算符

我们将上式的电子密度在全空间积分，就能得到总的粒子数算符：

这个算符的本征值只能为非负整数。

4. 场算符的非厄米性

我们已经知道：

因此，场算符不是厄米的，它本身不是一个可观测量。

激子模型

二次量子化表示多体项

上面我们所做的一切都是为了将式 (4) 最后一项的双体相互作用表达成二次量子化形式。在求得了场算符以后，我们重新回到上面的双体相互作用项：

我们将其写为关于 $V$ 的函数，$V$ 代表两体相互作用势函数，方便后续研究。

我们将 $V$ 从一次量子化（区分粒子编号）翻译成二次量子化（只看占据数）。对于任意算符 $\hat{O}$，我们都可以使用式 (52) $\sum_n |n\rangle \langle n| = \hat{1}$ 插入完备性表示：

其中 $|n\rangle$ 为初态，$\langle m|$ 为末态，$\langle m| \hat{O}|n\rangle$ 为算符 $\hat{O}$ 在该基底下的矩阵元。

因此，我们可以将相互作用势写作如下形式：

其中 $\langle kl | \hat{V} | mn \rangle$ 是矩阵元 $V_{klmn}$，它描述了这个过程发生的概率幅。我们可以将其展开为：

最右侧的 $\phi_m \phi_n$ 为初态，$\phi_k^* \phi_l^*$ 为末态，中间为相互作用。上式的物理含义是：初态两个粒子的概率幅 × 相互作用强度 × 末态两个粒子的概率幅。其中，相互作用强度代表权重——如果两个粒子距离近，权重就大。

算符 $|kl\rangle \langle mn|$ 的物理意义是：如果系统处于 $m, n$ 态，就将状态改为 $k, l$，写成二次量子化的形式为：

式 (72) 代表了：$m, n$ 态上的粒子在发生相互作用以后，被湮灭算符 $\hat{c}_n \hat{c}_m$ 消灭掉，然后使用产生算符 $\hat{c}_k^\dagger \hat{c}_l^\dagger$ 在 $k, l$ 态上产生粒子。

注意： 这里下标不能随意替换，需要保证和矩阵元一致。最外层对应粒子 1 从 $m$ 态到 $k$ 态的转变，内层代表粒子 2 从 $n$ 态到 $l$ 态的转变。

我们将式 (71a) 重新写成二次量子化的形式：

并且我们了解了中间 $\langle kl | \hat{V} | mn \rangle$ 可以用积分写出来。因此，整个双体相互作用可以写作：

我们将带有相同空间变量的项分组放在一起：

根据之前推导的场算符的定义（式 58b 和式 59a），我们可以对应上面方括号内的四项：

我们就得到了最终的二次量子化表述的双体相互作用公式：

或者，为了保证 $\vec{r}’$ 处和 $\vec{r}$ 处的算符是配对的，可以写作：

双带近似

在实际的晶体中存在大量能级，但是为了研究激子，我们仅关注激子涉及的轨道，例如包含价带顶和导带底的基态激子。因此我们先考虑最简单的情况——价带顶的一个电子被激发到导带底（单激子态）。

场算符的完整定义是所有能级态的叠加：

其中 $n$ 是能带索引，$\vec{k}$ 是波矢。$\phi_{n,\vec{k}}(\vec{r})$ 表示：处于第 $n$ 个能带、动量为 $\vec{k}$ 的那个电子，在位置 $\vec{r}$ 处的形状（权重）。在 DFT 计算中，它就是 Kohn–Sham (KS) 轨道。

这里我们使用双带近似——仅保留价带 $n=v$ 和导带 $n=c$，将其限制在这两个带上，且仅包含单个 $\vec k$ 点：

将式 (81) 代入式 (78)，可以得到仅包含价带、导带的两体相互作用：

上式展开后总共 16 项，代表了仅包含价带、导带所有可能的过程。

注意： 此时我们已经认为价带顶有一个电子被激发到导带底，已经形成了激子。我们现在要做的是研究已经形成的电子-空穴对之间的相互作用。

然而，这 16 项中并非所有项都有意义。例如：

1. $\hat{\psi}_c^\dagger \hat{\psi}_c^\dagger \hat{\psi}_c \hat{\psi}_c$：全是导带电子之间的相互作用，但我们考虑的单激子态此时导带仅有 1 个电子，双湮灭作用得 0，因此这一项在单激子态等于 0。
2. $\hat{\psi}_v^\dagger \hat{\psi}_v^\dagger \hat{\psi}_v \hat{\psi}_v$：代表激发完一个电子以后，价带中 $N-1$ 个价带电子两两之间的相互作用——这是价带背景的自相互作用，通常被吸收进基态能量或单粒子项中，不直接影响激子的结合能。
3. $\hat{\psi}_c^\dagger \hat{\psi}_c^\dagger \hat{\psi}_v \hat{\psi}_v$：同时消灭两个价带电子，产生两个导带电子（双激发），不属于单激子过程。

因此，我们需要的是在 $\vec{r}$ 处理导带电子（$\hat{\psi}_c$ 与 $\hat{\psi}_c^\dagger$）以及在 $\vec{r}’$ 处理价带电子（$\hat{\psi}_v$ 与 $\hat{\psi}_v^\dagger$）的项。由于 $\vec{r}$ 与 $\vec{r}’$ 是等价的（积分对称性），可抵消掉前面的 $1/2$ 系数，因此保留的项仅为：

算符从右向左作用，其含义依次是：

1. $\hat{\psi}_c(\vec{r})$：在 $\vec{r}$ 处消灭那个唯一的导带电子。
2. $\hat{\psi}_v(\vec{r}’)$：在 $\vec{r}’$ 处消灭一个价带背景电子。
3. $V(\vec{r}-\vec{r}’)$：计算这两个电子之间的库仑相互作用。
4. $\hat{\psi}_v^\dagger(\vec{r}’)$：把那个价带电子重新放回 $\vec{r}’$。
5. $\hat{\psi}_c^\dagger(\vec{r})$：把导带电子重新放回 $\vec{r}$。

双电子到电子-空穴对

在上面的式 (82b) 中，我们描述的全是电子和电子之间的相互作用。然而，空穴本质就是电子的缺失，因此我们可以改用空穴产生/湮灭场算符代替价带电子算符。对应关系为：消灭价带电子 $\hat{\psi}_v$ 等于产生空穴 $\hat{\psi}_h^\dagger$；产生价带电子 $\hat{\psi}_v^\dagger$ 等于消灭空穴 $\hat{\psi}_h$。直接替换可得：

但是目前的顺序并不是"左产生右湮灭"的标准正规序，因此我们需要交换算符顺序。由于是费米子，交换顺序要带一个负号。根据反对易关系 ${ \hat{\psi}_h(\vec{r}’), \hat{\psi}_h^\dagger(\vec{r}’) } = \delta(\vec{r}’-\vec{r}’)$ 可得：

因为算符坐标都是 $\vec{r}’$，所以 $\delta(\vec{r}’ - \vec{r}’) = \delta(\mathbf 0)$，这代表无穷大的局域贡献，需要通过正则化吸收到单粒子能量中。具体地，把这一项代入式 (83) 的积分可得：

其中 $\hat{\psi}_c^\dagger(\vec{r}) \hat{\psi}_c(\vec{r})$ 是导带电子的密度算符，中间积分描述了导带电子与所有价带电子之间的库仑相互作用。因此整体的物理意义是：导带的那个电子与所有价带背景电子的库仑作用。我们可以将这一项并入到单粒子哈密顿量的本征值中去，双体相互作用部分仅计算激子相对于基态的能量变化。

因此，我们最终得到了：

上式中，本来应该是 $V$ 的电子斥力，由于交换算符带来变号，变成了电子-空穴吸引相互作用，这就是激子能够结合在一起的原因。

注意：

• 哈密顿量 $H_{\text{exc}}$ 中仅显式包含双体相互作用项，但在数值计算上已经隐含了背景电荷的修正。例如计算带隙得到 1.5 eV，这个数值已经扣除了背景电荷、核作用等。
• $H_{\text{exc}}$ 在激子理论中等于电子-空穴的有效双体相互作用，是已经包含了所有背景电子与离子影响后的、被屏蔽的双体相互作用。当我们选择价带、导带的那一刻，它已经包含了所有价带电子彼此的排斥、激子与离子背景的平衡效果——这是一个我们经过重整化后的"真空能级"。

激子哈密顿量和激子态

我们前面推导出来的只是相互作用势能项。现在我们写出激子的整体哈密顿量，同时包含动能项和势能项：

第一项为电子动能项，其中 $m_e^$ 为电子有效质量；第二项为空穴动能项，其中 $m_h^$ 为空穴有效质量[^c8]；第三项为我们以上推导出的电子-空穴吸引相互作用能。在电子动能项中还包含一项带隙值 $E_g$，这是电子跃迁到导带必须跨过的能量。

我们现在得到了激子哈密顿量，需要将其作用在激子波函数上。激子波函数我们用如下形式表示：

从右往左看：$| G \rangle$ 是基态（价带满电子、导带全空）；$\hat{\psi}_c^\dagger(\vec{r}_e) \hat{\psi}_h^\dagger(\vec{r}_h)$ 在 $\vec{r}_h$ 位置创造一个价带空穴，并在 $\vec{r}_e$ 位置创造一个导带电子；$\Phi(\vec{r}_e, \vec{r}_h)$ 是电子出现在 $\vec{r}_e$ 且空穴出现在 $\vec{r}_h$ 的概率幅（在均匀晶体中仅与相对距离 $\vec{r}_e - \vec{r}_h$ 有关）。

说明： 激子波函数的完整图像可以理解为一个弥散在整个晶体中的氢原子（电子-空穴束缚态）。

Wannier 方程

我们将上面得到的激子哈密顿量作用在激子态上。由于项数较多且每一项较复杂，我们将每一项分开看。

(1) 电子动能项作用：

为了简化，记 $\hat{O}_e = -\frac{\hbar^2 \nabla^2}{2m_e^*} + E_g$。我们主要关注的是算符部分如何相互作用：

• 在电子动能哈密顿量 $\hat{T}_e$ 中，最右侧的电子湮灭算符 $\hat{\psi}_c(\vec{r})$ 会越过数 $\Phi(\vec{r}_e, \vec{r}_h)$，因为它是一个复数（波函数的幅值），算符可以穿透它。因此求和与积分均可往前提。
• 后面三个算符作用在基态上：

前两个都是导带相关的算符，第三个是价带（空穴）产生算符，因此前两个算符可以整体与第三个算符互换（它们之间反对易子为零，互不干扰），但需注意每次互换带来的负号。将 $\hat{\psi}_h^\dagger(\vec{r}_h)$ 向前提，可得：

注意： 这里交换算符位置时，因为反对易关系需要相应地加负号；最终偶数次交换符号回到正号。

依旧看电子湮灭算符 $\hat{\psi}_c(\vec{r})$，它和紧跟在后面的 $\hat{\psi}_c^\dagger(\vec{r}_e)$ 满足反对易关系（式 59b）：${ \hat{\psi}_c(\vec{r}), \hat{\psi}_c^\dagger(\vec{r}_e) } = \delta(\vec{r} - \vec{r}_e)$。可得：

这个式子作用在 $|G\rangle$ 上：

• 第一项 $\delta(\vec{r} - \vec{r}_e) |G\rangle$ 保留。
• 第二项 $-\hat{\psi}_c^\dagger(\vec{r}_e) \hat{\psi}_c(\vec{r}) |G\rangle$ 中，湮灭算符先作用在基态上，但基态的导带本身就没有电子，因此这项直接等于 0。

整理可得：

注意： 这里微分算符 $\hat{O}_e(\vec{r})$ 作用在 $\delta$ 函数上，对 $\delta$ 函数求导。根据 $\delta$ 函数的性质[^c9]：

$$\int d\vec{r} \, f(\vec{r}) \, \nabla^2 \delta(\vec{r}-\vec{r}_e) = \nabla^2 f(\vec{r}_e) \tag{93}$$

积分 $\int d\vec{r}$ 和 $\delta(\vec{r} - \vec{r}_e)$ 相抵消——所有 $\vec{r}$ 变成 $\vec{r}_e$：$\hat{\psi}_c^\dagger(\vec{r})$ 变成 $\hat{\psi}_c^\dagger(\vec{r}_e)$，微分算符 $\hat{O}_e(\vec{r})$ 变成 $\hat{O}_e(\vec{r}_e)$；并且因为 $\Phi$ 是关于 $\vec{r}_e$ 的函数，这个微分算符现在直接对 $\Phi$ 里的 $\vec{r}_e$ 求导。

整理可得：

(2) 空穴动能项作用：

类似电子动能项，这里不再做详细推导。同样地，所有的 $\vec{r}$ 都变成了 $\vec{r}_h$，可得：

(3) 相互作用项作用：

将系数和积分全部提到前面：

我们只关注两组算符之间的乘积：

其中电子和空穴湮灭算符分别和右边的对应产生算符运用反对易关系：与"反对易子" $\delta$ 函数贡献结合，另一项中湮灭算符先作用于基态而得 0。因此式 (97) 可以写作：

把原本的式子拼回去：

双重积分被两个 $\delta$ 函数分别抵消——$\delta$ 函数强制让 $d\vec{r}$ 积分塌缩到 $\vec{r} = \vec{r}_e$，同理 $d\vec{r}’$ 塌缩到 $\vec{r}’ = \vec{r}_h$。最后整理可得：

(4) 合并得本征方程

最后，我们将上面所有项加起来，让它等于激子总能量 $E |\Psi_{\text{exc}}\rangle$：

由于左右两边的基底向量相同且对所有项均成立，大括号 ${\dots}$ 内的内容必须相等。于是我们仅保留系数方程：

上式被称为 Wannier 激子方程。

我们将 Wannier 方程和氢原子薛定谔方程比较一下：

上式为氢原子方程。可以看出，Wannier 方程实际上描述的是晶体中的"氢原子模型"，二者的数学形式完全一致。

GW + BSE

我们上面推导的 Wannier 方程实际上包含了三个关键的近似，导致它适用的是大半径激子，即 Wannier–Mott 激子（或直接称为 Wannier 激子）：

1. 有效质量近似——假设能带的带边附近是完美的抛物线：

   $$E_{c\mathbf{k}} - E_{v\mathbf{k}} \approx E_g + \frac{\hbar^2 k^2}{2\mu^*} \tag{104}$$

   我们将复杂的能带信息简化成了有效质量。

2. 介质屏蔽近似——相互作用是被极度简化的：

   $$W \approx \frac{1}{\epsilon_r} \frac{e^2}{|\mathbf{r}_e - \mathbf{r}_h|} \tag{105}$$

   在这里我们忽略了微观原子结构，把背景当作连续介质。

3. 忽略交换项——因为 Wannier 激子有较大的半径（约 $10^1\sim10^2$ Å），电子和空穴的短程交换作用很弱，因此在 Wannier 方程中通常忽略这一项。这会导致使用 Wannier 方程计算的单重态/三重态劈裂不准确。

对于一般情形，钙钛矿等半导体中的激子大都是大半径激子，使用 Wannier 方程已经足够好。但是如果要最精确计算，必须使用 GW + BSE。

GW 近似

我们前面推导的 Wannier 方程是基于多粒子表象的，而 GW 近似处理的是单粒子。它的目的是把 DFT 里假想的 Kohn–Sham 轨道修正为"真实可观测的准粒子能量"。换句话说，GW 是用来"把电子当成真正带着自能的粒子"，而不是 DFT 里的平均场电子。

为了替换掉 DFT 中那个近似的 $V_{xc}$，我们需要得到一个非局域的、含时的自能算符 $\Sigma$。

单粒子格林函数

为了能够解出包含自能 $\Sigma$ 的方程（Dyson 方程），我们首先需要建立"传播"和"散射"的概念。

我们定义单粒子格林函数：

它的每一项含义如下：

• $1, 2$：代表时空坐标 $(\vec{r}_1, t_1)$ 和 $(\vec{r}_2, t_2)$。
• $\hat{\psi}^\dagger(2)$：产生算符，意思是"在 $t_2$ 时刻、$\vec{r}_2$ 位置创造一个电子"。
• $\hat{\psi}(1)$：湮灭算符，意思是"在 $t_1$ 时刻、$\vec{r}_1$ 位置消灭一个电子"。
• $T$：编时算符（time-ordering）。它强制要求：时间晚的算符放在左边，时间早的放在右边（对费米子算符每交换一次还要带一个负号）。
• $\langle \Psi_0 | \dots | \Psi_0 \rangle$：期望值，表示这个过程是在系统基态 $|\Psi_0\rangle$ 上发生的。
• $|G|^2$ 代表电子从 $2$ 运动到 $1$ 的概率密度。

格林函数包含两个物理过程，分别是电子传播和空穴传播：

1. 电子传播（$t_1 > t_2$）：因为 $t_1$ 较晚，产生算符先作用，在 $\vec r_2$ 处加一个电子，然后电子运动，经过 $t_1 - t_2$ 时间后湮灭算符把电子拿走。这段过程描述了导带电子的运动。
2. 空穴传播（$t_1 < t_2$）：因为 $T$ 算符的作用，优先让湮灭算符在 $\vec r_1$ 处拿走一个电子（创造空穴），然后空穴运动，经过 $t_2 - t_1$ 时间后产生算符把电子放回。这段过程描述了价带空穴的运动。

自由电子格林函数把时域和频域联系起来。我们先写出时域格林函数 $G_0(t)$：根据量子力学，如果一个粒子的能量是 $\epsilon_{\mathbf{k}}$，它的波函数随时间的演化遵循：

注意： 为了方便，下面取原子单位制 $\hbar=1$。

而格林函数 $G_0(t)$ 描述的是：在 $t=0$ 时刻产生一个粒子，它在 $t$ 时刻还存在的概率幅。粒子必须是先产生后传播，因此它应该只在 $t>0$ 时存在。我们引入阶跃函数 $\theta(t)$，$G_0(t)$ 就可以写作：

其中：

• $\theta(t)$ 在 $t>0$ 时为 1、在 $t<0$ 时为 0，保证因果律。
• $e^{-i \epsilon_{\mathbf{k}} t}$：相位随时间旋转，这一项代表粒子有能量 $\epsilon_{\mathbf{k}}$。
• 前面的 $-i$ 是因为薛定谔方程左边有一个 $i$，求解的时候到了分母上，就是 $-i$。

下面我们将格林函数转换到频域 $\omega$，对其做傅里叶变换：

代入式 (108) 可得：

注意： 因为 $\theta(t)$ 的存在，积分下限从 $-\infty$ 变成了 $0$。

这里出现数学问题：当 $t \to +\infty$ 时，指数函数是振荡且不衰减的，这个积分严格不收敛。因此我们引入一个假设——粒子不是永生的，会随时间极其缓慢地消失。这样在指数上加入一个无穷小的衰减因子 $-\eta t$（其中 $\eta \to 0^+$），上式变成：

经过数学推导后，就得到了：

• $+i\eta$ 这一项在物理上代表的边界条件是：保证了粒子只能在 $t>0$ 传播。
• 当分母为 0 时（$\omega = \epsilon_{\mathbf{k}}$），即在 $\epsilon_{\mathbf{k}}$ 处出现一个 $\delta$ 函数式的尖峰——这是电子的本征能量。
• 从时域上看，$G_0(t)$ 的模长为 1（$\eta\to 0$ 极限下），因此粒子寿命无限。

说明：

• 在了解 $G(t)$ 后，我们也要了解它的逆 $G^{-1}$。
• 格林函数在数学上的定义是某个微分算符的逆，因此 $G^{-1}$ 本质上就是薛定谔方程里的微分算符。
• 含时薛定谔方程可以写作：

$$i \frac{\partial}{\partial t} \psi(t) = H \psi(t) \tag{113}$$

• 代入格林函数并整理可得：

$$\left( i \frac{\partial}{\partial t} - H \right) G(t) = \delta(t) \tag{114}$$

• 在数学上，把微分算符作用在 $G$ 上会得到一个 $\delta$ 函数。
• 把 $\left( i \frac{\partial}{\partial t} - H_0 \right)$ 记为 $\hat A$。
• 那么 $\hat A G = I$，因此 $G$ 只能是 $\hat A$ 的逆。
• 因此，$G^{-1}$ 等于 $\left( i \frac{\partial}{\partial t} - H_0 \right)$，也就是整体的微分算符。
• 因此，$G(t)$ 描述粒子的运动状态（解），$G^{-1}$ 描述粒子的运动方程（动力学）。
• 频域下的 $G^{-1}$ 为：

$$G^{-1}(\mathbf{k}, \omega) = \omega - \epsilon_{\mathbf{k}} + i\eta \tag{115}$$

• 在频域的表达下，它变成了简单的代数，如果再令 $G^{-1}=0$，那么就得到了本征关系：$\omega = \epsilon_{\mathbf{k}}$。

单粒子格林函数满足如下算符方程（后面我们称之为 Dyson 方程）：

上面我们已经了解了格林函数的逆相当于哈密顿（运动方程）算符。因此，可以写出上式各项的含义：

• $G^{-1}(1,2)$：相互作用系统格林函数的逆。
• $G_0^{-1}(1,2)$：无相互作用（含外势）系统的格林函数逆，完整展开为 $G_0^{-1}(1,2) = \left[i\hbar \frac{\partial}{\partial t_1} + \frac{\hbar^2 \nabla_1^2}{2m} - V_{ext}(1)\right] \delta(1,2)$。
• $V_H(1)$：经典平均静电势（Hartree 势）。
• $\delta(1,2) = \delta(\vec{r}_1-\vec{r}_2)\delta(t_1-t_2)$：时空 $\delta$ 函数。
• $\Sigma(1,2)$：自能函数，包含了所有超越平均场的多体效应（交换与关联）。

说明：

• 准粒子总有效哈密顿量算符可写作：

$$G^{-1}(1,2) = \left[ i\hbar \frac{\partial}{\partial t_1} - H_{\text{total}}^{\text{eff}} \right] \delta(1,2) \tag{117}$$

上式代表了真实[^d1]准粒子遵守的运动方程算符。若 $G^{-1} \psi = 0$，那就是真实准粒子的薛定谔方程。

• 自由格林函数的逆：

$$G_0^{-1}(1,2) = \left[ i\hbar \frac{\partial}{\partial t_1} + \frac{\hbar^2 \nabla_1^2}{2m} - V_{ext}(1) \right] \delta(1,2) \tag{118}$$

上式表示单电子在仅存在核作用下的运动方程算符，这部分没有任何电子之间的关联。

• Hartree 项：

$$V_H(1) \delta(1,2) = \int d3 \, v(1,3) \rho(3) = \int d\mathbf{r}_3 \frac{e^2 \rho(\mathbf{r}_3)}{|\mathbf{r}_1 - \mathbf{r}_3|} \, \delta(1,2) \tag{119}$$

上式表示电子在 $\mathbf{r}_1$ 感受到的 Hartree 势，来源是其他所有电子在 $\mathbf{r}_3$ 处的密度。它仅包含电子之间的经典排斥相互作用，不含量子的交换-关联作用。

• 最后一项自能 $\Sigma(1,2)$ 没有显式的封闭表达式。它包含了所有经典平均场无法描述的量子效应。
• 因此，经过全部加和我们可以写出准粒子波动方程：

$$\left[ i\hbar \frac{\partial}{\partial t_1} + \frac{\hbar^2 \nabla_1^2}{2m} - V_{ext}(1) - V_H(1) \right] \psi(1) - \int d2 \, \Sigma(1,2) \psi(2) = 0 \tag{120}$$

• 在上式中，我们对原算符方程两边同乘 $\psi(2)$ 并对 $2$ 积分。因为 $G^{-1}(1,2)$ 将状态从位置 $2$ 映射到位置 $1$，为了得到位置 $1$ 的状态，必须把位置 $2$ 的波函数作为输入。
• 矩阵形式上，这相当于 $y_i = \sum_j A_{ij} x_j$：我们为了得到 $y_i$，输入必须是 $x_j$。
• 式 (120) 中除了最后一项，所有的 $\delta(1,2)$ 都把对 $2$ 的积分塌缩掉了；只有最后一项不是对角矩阵，因此积分必须保留。

最终，我们可以将全相互作用系统的逆格林函数 $G^{-1}(\mathbf{r}_1, t_1; \mathbf{r}_2, t_2)$ 写成如下形式：

上式方括号内的项都是局域的，因为 $\delta(\mathbf{r}_1 - \mathbf{r}_2)\delta(t_1 - t_2)$ 的存在使仅有对角项起作用。最后一项为非局域项，我们将其写成矩阵形式：

其中，矩阵中每一项都是能量 $\omega$ 的函数。对角元和非对角元具有不同的含义：

(1) 对角元——轨道自己对自己的影响

它包含实部和虚部：

• 实部 $\text{Re}[\Sigma_{ii}(\omega)]$ 代表能量的重整化——电子在轨道 $i$ 上时极化了周围的环境，导致自己的能量发生变化。因此这项修正了 DFT 计算出的能级：

  $$E_i^{QP} \approx \epsilon_i^{DFT} + \text{Re}[\Sigma_{ii}(E_i)] - V_{xc,ii} \tag{123b}$$

  上式代表准粒子能量 = DFT 计算的能量 + GW 计算出的自能 − DFT 中近似的 $V_{xc}$（因为 DFT 计算已经包含了 $V_{xc,ii}$，要去掉以避免重复计数）。

• 虚部 $\text{Im}[\Sigma_{ii}(\omega)]$ 代表寿命的衰减——电子在轨道 $i$ 上会因为散射（电子-电子散射、电声耦合等）而跃迁到其他态，虚部越大寿命越短。

(2) 非对角元——两个轨道之间的相互作用修正

在一般情况下，$H_{ij}$ 非对角元代表了电子从 $j$ 跳到 $i$ 的概率。但在自能中——正如我们上面讨论的——电子会自己极化周围的环境导致能量变化，电子跳跃时也伴随周围被极化的环境一起跳跃。因此这一项包含了电子从 $j$ 态出发、受到扰动到了 $i$ 态、然后到达的复杂过程。

同样地，这一项也包含实部和虚部：

• 实部 $\text{Re}[\Sigma_{ij}]$ 代表跃迁的耦合强度。
• 虚部 $\text{Im}[\Sigma_{ij}]$ 代表跳跃过程中受到的阻力，即发生散射的强度。

在绝大多数半导体中，非对角项通常被忽略（称为 $G_0W_0$），此时对角元占主导——我们认为 DFT 计算出的波函数已经足够好，不太需要非对角元修正。但在强关联体系中，这是不可忽略的。

Dyson 方程

我们将格林函数满足的算符方程 (116) 写成最简单的形式：

上面我们已经了解了它的物理含义：

真实运动方程 = 自由运动方程 − 修正项

上式两边同时从左边乘以自由格林函数 $G_0$ 可得：

其中 $G_0 \cdot G_0^{-1}$ 为单位矩阵 $I$，整理可得：

等式两边右乘全格林函数 $G$，消掉左边的 $G^{-1}$：

整理可得：

上式就是 Dyson 方程，也是我们前面提到的算符方程（式 116）的等价形式。等式两边都有 $G$，看起来是递归的。我们把 $G$ 代入到上式右边的 $G$，整理可得：

如果一直代入，就会得到无限的迭代级数：

在 Dyson 方程中，我们把 $G_0$ 叫做传播子。上式的物理含义是：

粒子从 $i$ 到 $j$ 的总传播概率幅等于所有路径之和，包含粒子所有可能散射的路径叠加：

• 第一项 $G_0$ 代表粒子没有发生散射，直接从 $i$ 到 $j$；
• 第二项 $G_0 \Sigma G_0$ 代表粒子传播 $G_0$，散射了一次 $\Sigma$，然后又传播 $G_0$，到了 $j$；
• 第三项则是散射两次……

我们切换到频域空间来看式 (128) 的最简单情形（动量守恒下 $\Sigma$ 在 $\mathbf k$ 表象对角）：

将含有 $G$ 的项移到一边并代入 $G_0 = \frac{1}{\omega - \epsilon_k}$，可得：

可以看到，原本在 $\epsilon_k$ 的 DFT 能级现在移动到了 $\epsilon_k + \Sigma$ 的 GW 能级。它对应了能带图中的 $y$ 轴；$\mathbf k$ 代表了 $\mathbf k$ 空间中不同的 $\mathbf k$ 点。

当 $\omega = \epsilon_{\mathbf{k}} + \text{Re},\Sigma(\mathbf{k}, \omega)$ 时，分母实部为 0，$|G(\mathbf{k}, \omega)|$ 出现峰值，它代表了这里出现了一条能带。这个时刻的 $\omega$ 值就是准粒子能量。

双粒子格林函数

我们需要算格林函数 $G$，因为它包含了能带信息。但是仅仅从 $G$ 的定义来看，我们无法直接解出来具体的能带结构，需要找到一个 $G$ 满足的方程，然后解这个方程。

上式是格林函数的定义。我们了解算符的定义，但不清楚 $|\Psi_0\rangle$ 的具体形式——因为我们面对的是一个真实多体系统。

我们要找的能带本征能量，本质上是微分方程的特征值（对应时间导数 $i \frac{\partial}{\partial t}$）。因此，我们需要对格林函数关于时间求导；但格林函数中的 $T$ 算符需要先写成数学上的阶跃函数形式才能求导：

前面同样提到过，当 $t_1 > t_2$ 时仅保留右侧第一项，反之保留第二项；负号来自费米子算符交换。

因此，格林函数 $G$ 可以写作：

接下来对时间 $t_1$ 求导。求导会同时作用在 $\theta$ 函数和 $\hat\psi(1)$ 算符上。我们将对阶跃函数求导得到的部分称为边界项，将对 $\hat\psi(1)$ 求导得到的部分称为动力学项：

(1) 边界项： 对 $\theta$ 函数的求导是 $\delta$ 函数，即 $\frac{\partial}{\partial t_1} \theta(t_1 - t_2) = \delta(t_1 - t_2)$。因此：

利用反对易关系 ${\hat\psi(\mathbf{r}_1), \hat\psi^\dagger(\mathbf{r}_2)} = \delta(\mathbf{r}_1 - \mathbf{r}_2)$，边界项直接变成：

(2) 动力学项：

我们直接得到了海森堡运动方程的形式：

说明： 算符的时间演化属于海森堡绘景，必须遵守海森堡运动方程。

因此可以写作：

为了计算对易子，我们将哈密顿量拆分为动能和相互作用两部分：$\hat H = \hat H_0 + \hat H_{\text{int}}$[^d2]，其中 $\hat H_0$ 包含动能和外势能，$\hat H_{\text{int}}$ 包含电子间的相互作用能。

动能部分：

利用对易公式 $[A, BC] = {A, B}C - B{A, C}$：

其中，右侧第二项两个湮灭算符的反对易子为 0；第一项产生湮灭算符的反对易子为 $\delta(\mathbf r_1 - \mathbf x)$，与积分符号抵消：

相互作用部分：

需要计算对易子 $[\hat\psi(\mathbf{r}_1), \hat\psi^\dagger(\mathbf{x}) \hat\psi^\dagger(\mathbf{y}) \hat\psi(\mathbf{y}) \hat\psi(\mathbf{x})]$。看起来繁琐，但产生算符之间反对易子为 0，产生-湮灭反对易子是 $\delta$ 函数，所以最终化简得到：

这里 $3$ 代表空间中任何一个位置，由于积分符号，效果上是把所有位置 $3$ 处的电子贡献加起来。

把两部分合并整理：

我们可以看到上式右侧出现了 $\hat\psi^\dagger \hat\psi \hat\psi \hat\psi^\dagger$（四个场算符），因此被称为双粒子格林函数 $G_2$。

它说明了：如果我们要求单粒子 $G$，就要先求双粒子 $G_2$；要求 $G_2$，得先求 $G_3$……这变成了 BBGKY 等级，似乎陷入了死循环。

Hedin 方程组

Schwinger 导数

为了打破上式的死循环，Schwinger 引入了一个随时间变化的虚拟外场 $\phi$。经过数学推导，双粒子格林函数可以表示为：

代回到运动方程 (145)：

其中多出来的 $V_H(1)$ 是做 Schwinger 替换后吸收掉式 (146) 第二项得到的 Hartree 势。我们比较一下 Dyson 方程：

我们就能得到 $\Sigma$ 的形式：

上式中，$2$ 是起始点，$1$ 是终点，$3$ 是粒子发生散射的中间点。

这个方程的物理含义是：自能 $\Sigma$ 是电子通过库仑力 $v$ 和它自己引发的环境变化 $\frac{\delta G}{\delta \phi}$ 之间的耦合的结果——

• 左边 $G(3,2)$ 是传播子，代表电子从起点 $2$ 跑到中间点 $3$。$\Sigma(1,3)$ 则是电子在 $3$ 受到散射、到达位置 $1$。
• 右边 $v(1,3)$ 是位置 $1$ 和位置 $3$ 之间的库仑相互作用，$\frac{\delta G(1,2)}{\delta \phi(3)}$ 是电子从 $2$ 到 $1$ 的路途中，在位置 $3$ 发生散射的概率幅。

整体描述了：电子从 $2$ 到 $3$，经过 $3$ 时由于排斥周围电子在 $3$ 产生了一个局部带正电的环境；电子随后跑到 $1$，在 $1$ 受到刚刚 $3$ 形成的正电环境影响——受到 $v(1,3)$ 的库仑相互作用。两边描述的物理过程是一致的。

逆介电函数和顶点函数

我们有了式 (149)，为了解出自能 $\Sigma$，需要使用链式法则处理 $\frac{\delta G}{\delta \phi}$：

其中 $V_{\text{tot}}$ 是包括外场扰动后的总有效势：$V_{\text{tot}} = V_{ext} + V_H + \phi$（关心变化时可写 $\delta V_{\text{tot}} = \delta\phi + \delta V_H$）。

这里出现了第四个位置 $4$，因为我们在位置 $3$ 处施加了一个虚拟外场 $\phi(3)$，这会导致不仅 $3$ 位置的电势改变，全空间其他所有位置的电势都会变化——因此我们要将所有受影响的位置 $4$ 加起来。

我们把 $\frac{\delta V_{\text{tot}}}{\delta \phi}$ 定义为逆介电函数 $\epsilon^{-1}$：

它描述了外场 $\phi$ 如何引起总势场的变化（产生屏蔽效应）。

我们把式 (149) 右侧库仑势 $v$ 与逆介电函数的卷积定义为 $W(1,4)$：

这就是 GW 近似中的 $W$——屏蔽库仑相互作用。

我们再看链式法则得到的另一项。根据数学恒等式 $\frac{\delta G}{\delta V} = -G \frac{\delta G^{-1}}{\delta V} G$，我们定义：

因此可得：

$\Gamma$ 我们称之为顶点函数（vertex function）。它的物理意义是：在 $3$ 处的势场变化对电子从 $1$ 到 $2$ 的逆传播有什么影响——即电子在 $3$ 点感受到外场时发生的"非平凡散射"。

我们将得到的项全部代回式 (149)，最后可得（紧凑算符记号下）：

我们将顶点函数中的 $G^{-1}$ 替换成 Dyson 方程中的形式：

其中：

• 第一项：自由格林函数 $G_0$ 只包含动能和外势，与总势场无关，因此为 0。

• 第二项：

  $$-\frac{\delta [-V(1)\delta(1,2)]}{\delta V(3)} = \frac{\delta V(1)}{\delta V(3)} \delta(1,2) = \delta(1,3)\delta(1,2) \tag{158}$$

• 第三项含 $\Sigma$，需要进一步链式法则。

因此，整体可以写作：

对其中 $\frac{\delta \Sigma(1,2)}{\delta V(3)}$ 再使用链式法则引入 $G$（因为 $\Sigma$ 是 $G$ 的泛函，而 $G$ 依赖于 $V$）：

其中 $\frac{\delta G(4,5)}{\delta V(3)}$ 这项根据前面的式 (155) 可写为：

将其代回式 (159)，最终可以得到顶点方程：

库仑屏蔽和极化率

根据式 (152) 以及总有效势的定义，$\frac{\delta V(4)}{\delta \phi(3)}$ 可以写作：

其中 $\frac{\delta \phi(4)}{\delta \phi(3)} = \delta(4,3)$。Hartree 势的定义是：

它的导数可以写成：

我们依然可以使用链式法则，将 $\frac{\delta \rho}{\delta \phi}$ 展开：

其中根据极化率（不可约极化）的定义 $P(X, Y) = \frac{\delta \rho(X)}{\delta V(Y)}$，上式可以写作：

因此式 (153) 右侧第二项变成：

将上式代回 (153)，那么 $W$ 可以写作：

• 第一项为裸库仑相互作用，代表两个电子之间的纯库仑相互作用：

  $$\int d3 \, v(1,3) \delta(4,3) = v(1,4) \tag{170}$$

• 第二项为屏蔽修正：

  $$\int dX \, dY \, v(4, X) P(X, Y) \left[ \int d3 \, v(1,3) \frac{\delta V(Y)}{\delta \phi(3)} \right] \tag{171}$$

  屏蔽修正代表了：电子 A 排斥了周围的电子并"吸引"空穴，导致电子 A 周围形成了一圈带正电的极化云；电子 B 原本仅感受到 A 的排斥作用，但因为正极化云的存在，被抵消了一部分——这就是屏蔽。

最终可以得到 $W$ 的 Dyson 形式：

如果我们希望用格林函数 $G$ 把极化率 $P$ 表达出来，就能把微观量（$G$）和宏观响应（$P$）联系起来。

在前面（式 68），我们推导出过密度算符 $\hat\rho(1) = \hat\psi^\dagger(1) \hat\psi(1)$。根据格林函数的定义 $G(1,2) = -i \langle T \hat\psi(1)\hat\psi^\dagger(2) \rangle$，如果 $2$ 趋近于 $1$（取时间无穷小正向极限 $1^+$），那么粒子数密度可以用格林函数写出：

这样就能根据极化率的定义写出格林函数的表达：

由式 (155) $\frac{\delta G}{\delta V} = G\Gamma G$，直接代回式 (175)，可以写出极化率：

封闭循环与 $G_0W_0$

根据上面的推导，我们得到了五个相互依赖的方程——Hedin 方程组。重新写出并看一下它们的物理意义：

(1) 极化率方程：

它描述了介质对外场的响应能力。尽管 $G \times G$ 看起来仅描述了一对粒子，但格林函数本身已经包括了所有电子能带和所有量子态信息的加和，因此 $P$ 表示的是整个电子海的集体响应。顶点 $\Gamma$ 描述了这一产生过程中的相互作用修正。

(2) 屏蔽相互作用方程：

它描述了：裸库仑力 $v$ 通过被介质极化云 $P$ 反复散射，最终变成屏蔽相互作用 $W$。它表明粒子产生的极化场会导致库仑力变弱。

(3) 自能方程：

它描述了电子传播过程中：通过与自身形成的极化场耦合所产生的能量修正，以及与所有其他电子发生非经典、非局域的纠缠。中间的耦合强度由 $\Gamma$ 决定。

(4) 顶点方程：

它描述了粒子间相互作用的顶角修正。真实的电子不仅仅是一个点电荷，它周围笼罩着一团复杂的极化云，使得它与外界场的相互作用不再是简单的点接触（$\delta$ 函数），而是复杂的散射过程。

(5) 格林函数 Dyson 方程：

它把所有电子间的相互作用全部封装到 $\Sigma$ 中，加到自由电子 $G_0$ 上，从而得到相互作用系统下真实的传播子 $G$。

上面五个方程构成一个封闭的自洽循环：

如果我们知道初始的 $G$ 和 $\Gamma$（假设 $\Gamma = \delta\delta$），我们就能算出介质的极化 $P$；极化又能算出屏蔽相互作用 $W$（相互作用不是常数，而是由介质状态决定的）；我们就能根据 $W, G, \Gamma$ 算出自能 $\Sigma$（如果不知道电子在什么样的介质中运动，就不知道它会被修正多少）；这时根据新的自能又能算出新的顶点修正 $\Gamma$[^d3]（当电子略微改变运动方式时，环境立刻重新组织响应它）；自能一旦改变，电子的传播 $G$ 就必须改变——回到循环开始。

上面的方程组意味着：电子不是在一个固定环境中运动，而是在一个被它自己不断重塑的环境中运动。构成循环的原因是：电子的传播、介质的响应和相互作用的有效形式并非独立物理量，而是由同一批电子共同决定的集体结果；任何一个量的改变都会通过其他量反馈回来，迫使它必须以循环、自洽的方式求解。

那么我们能不能像 DFT 那样直接迭代求解 GW 方程组？

回到 $\Sigma(\omega)$ 中，前面推导过 $\Sigma(\omega)$ 是存在虚部的，它代表了粒子的寿命（能量-时间不确定关系）：

这意味着 GW 里的电子不再是一个锐利的粒子，而是一个有宽度的波包。GW 近似中电子一定会因为散射而到达别的状态，因此 $\Delta t$ 是有限的，$\Delta E$ 必须是非零的，这会导致粒子在能谱中出现能量展宽。但是 DFT 里电子寿命是无限[^d4]的，因此能带图每条线都是无限细的。当我们把有宽度的波包粒子代回去重新计算新的 $W$，得到的屏蔽效应会变弱、相互作用变强，粒子的寿命就会更短，波包会更展宽，粒子更模糊。因此迭代后的能带不再清晰，而是糊成一团。

更重要的是，GW 不是变分理论。我们使用自洽迭代求解 KS 方程的时候，理论基础是变分原理：存在一个明确的泛函，其真实解是该泛函的极值（通常是最小值）。任意试探波函数的能量期望值一定大于等于基态能量，即 $E[n] \geq E[n_0]$，唯一变量是电子密度。但 GW 的方程不是最小化某个物理量的结果，而仅仅是满足某些自洽关系，因此迭代不能保证收敛到物理解。

为了打开这组方程组的循环，我们直接忽略顶点修正，强行让 $\Gamma(1,2;3) = \delta(1,2)\delta(1,3)$（即纯粹的点接触），代回后所有的 $\Gamma$ 积分都化简。这时，极化率变成了：

注意： 这里通过 $\Gamma$ 的两个 $\delta$ 强制 $3 = 4 = 2$，得到的 $\delta$ 函数消掉了两个积分变量。

此时极化率变成了电子和空穴的简单乘积。自能方程通过同样的化简变成：

这就是为什么它叫做 GW 近似——电子自能只取决于传播 $G$ 与屏蔽相互作用 $W$ 的乘积。屏蔽方程由于极化率的简化而被简化，但形式没有发生改变。

之所以叫准粒子近似：这里的电子是有寿命的，最终会散射消失；并且它不仅是一个电子，还包括了它周围被它极化的整体电子云。

虽然我们现在知道这组方程怎么求解，但是 GW 的第一步是构建无相互作用格林函数 $G_0$，即：

注意：

• 这是实空间表象下的格林函数，而我们前面推导的式 (112) 则是本征态表象下的。我们可以使用完备性关系：

$$\hat{1} = \sum_{n\mathbf{k}} |n\mathbf{k}\rangle \langle n\mathbf{k}| \tag{184}$$

• 将格林函数转换为：

$$\begin{aligned} G(\mathbf{r}, \mathbf{r}') &= \langle \mathbf{r} | \hat{G} | \mathbf{r}' \rangle \\ &= \langle \mathbf{r} | \left( \sum_{n\mathbf{k}} |n\mathbf{k}\rangle \langle n\mathbf{k}| \right) \hat{G} \left( \sum_{m\mathbf{q}} |m\mathbf{q}\rangle \langle m\mathbf{q}| \right) | \mathbf{r}' \rangle \end{aligned} \tag{185}$$

• 因为 $\hat{G}$ 在本征基底下是对角的（$\langle n\mathbf{k} | \hat{G} | m\mathbf{q} \rangle = G_{n\mathbf{k}} \delta_{nm}\delta_{\mathbf{kq}}$），因此双重求和退化为单重求和：

$$\begin{aligned} G(\mathbf{r}, \mathbf{r}') &= \sum_{n\mathbf{k}} \underbrace{\langle \mathbf{r} | n\mathbf{k}\rangle}_{\phi_{n\mathbf{k}}(\mathbf{r})} \cdot \underbrace{\langle n\mathbf{k} | \hat{G} | n\mathbf{k} \rangle}_{\frac{1}{\omega - \epsilon_{n\mathbf{k}}}} \cdot \underbrace{\langle n\mathbf{k} | \mathbf{r}' \rangle}_{\phi_{n\mathbf{k}}^*(\mathbf{r}')} \\ &= \sum_{n\mathbf{k}} \frac{\phi_{n\mathbf{k}}(\mathbf{r}) \phi_{n\mathbf{k}}^*(\mathbf{r}')}{\omega - \epsilon_{n\mathbf{k}}} \end{aligned} \tag{186}$$

对于 $G_0$ 来说，我们需要上式中的波函数和能级。DFT 的能量虽然不准确，但波函数通常是可靠的。因此 GW 近似不再去修正波函数，而是去修正能量。

我们现在要求解准粒子哈密顿量 $H_{QP}$，已有 DFT 的哈密顿量 $H_{DFT}$：

因此我们需要解出来：

虽然 $\Sigma$ 是非对角的，但 $\boldsymbol V_{xc}$ 在 KS 轨道基底下也是近似非对角占优的，且二者非对角元的大小往往可比，因此微扰项 $(\Sigma - V_{xc})$ 的非对角元会大量抵消，仅保留对角占优结构。因此我们仅计算对角元：

上式中 $\epsilon_{n}^{DFT}$ 是 DFT 计算的本征值；$\Sigma(\epsilon_n) - V_{xc}$ 是用自能替换掉 DFT 中近似的交换关联势；$Z_n$ 是准粒子重整化因子（$Z_n = \left(1 - \partial \Sigma / \partial \omega \big|_{\omega=\epsilon_n}\right)^{-1}$），代表准粒子峰的相干谱权重。$Z_n$ 介于 0 和 1 之间——典型半导体中为 $0.7\sim0.8$，剩下的权重转移到非相干"卫星峰"上。

因此，在得到 DFT 计算的 $\epsilon_0$ 和 $\phi_0$ 后，构建格林函数 $G_0$ 和屏蔽 $W_0$，算出自能 $\Sigma = i G_0 W_0$，进而得到准粒子能量：

这时直接输出准粒子能量 $E_{QP}$，不再进行循环——此时得到的就是 $G_0 W_0$ 近似。

Bethe–Salpeter 方程

BSE 哈密顿量

尽管我们在 GW 过程中涉及了"电子-空穴对"型的极化率 $P = -i GG$，但 GW 关心的本质上仍是单个电子——极化率的作用只是改变了单个电子感受到的有效势，我们关心的也只是单个电子被屏蔽后的样子。

因此，我们需要使用 Bethe–Salpeter 方程（BSE） 去解电子-空穴对的薛定谔方程，才能得到电子和空穴是如何相互作用、形成束缚激子的。我们的目的很明确：得到双粒子格林函数。

为了得到双粒子格林函数，我们需要解以下方程：

其中 $U(2’, 2)$ 是外加势场（非局域），它代表在 $2$ 消灭一个粒子、在 $2’$ 产生一个粒子。$\delta U$ 代表这只是一个微扰，$\delta G$ 是电子传播受到的相应变化。

物理上：在 $(2’, 2)$ 对系统施加一个微扰，那么在 $(1, 1’)$ 的电子传播状态会发生什么变化。左侧的 $L$ 本质上是电子-空穴对的传播子。

说明： 为什么对电子传播施加一个微扰就变成了双粒子的传播子？我们做以下推导。

• 在量子力学中，系统的状态随时间演化由时间演化算符 $S$ 控制：

$$S = \exp\left( -i \int H(t) \, dt \right) \tag{193a}$$

• 我们在哈密顿量 $H$ 里施加一个外场 $U(\mathbf{r}, t)$，它会和电子密度 $\hat{\rho} = \hat{\psi}^\dagger \hat{\psi}$ 耦合：

$$H_{\text{total}} = H_0 + \int U(\mathbf{r}, t) \, \hat{\psi}^\dagger(\mathbf{r}) \hat{\psi}(\mathbf{r}) \, d\mathbf{r} \tag{193b}$$

• 因此，时间演化算符变为：

$$S[U] = \exp\left( -i \int (H_0 + U \cdot \hat{\psi}^\dagger \hat{\psi}) \, dt \right) \tag{194a}$$

• 如果我们对 $U$ 求导，那么 $U$ 的系数（$\hat{\psi}^\dagger \hat{\psi}$）会出现在前面：

$$\frac{\delta S}{\delta U} = \left( -i \hat{\psi}^\dagger \hat{\psi} \right) \times \exp(\dots) \tag{194b}$$

• 但我们要求 $\frac{\delta G(1,1’)}{\delta U(2’,2)}$，需要使用链式法则。在此之前，要写出格林函数包含时间演化算符 $S$ 的完整定义：

$$G(1, 1') = -i \frac{ \langle \Phi_0 | T [ S \, \hat{\psi}(1) \hat{\psi}^\dagger(1') ] | \Phi_0 \rangle }{ \langle \Phi_0 | S | \Phi_0 \rangle } \tag{195}$$

• 因此可以写出：

$$\frac{\delta G}{\delta U} = \frac{-i \langle T \, S \, \hat{\psi}(1) \hat{\psi}^\dagger(1') \hat{\psi}^\dagger(2) \hat{\psi}(2') \rangle}{\langle S \rangle} - G(1, 1') \cdot \frac{-i \langle T \, S \, \hat{\psi}^\dagger(2) \hat{\psi}(2') \rangle}{\langle S \rangle} \tag{196}$$

• 第一项就是双粒子格林函数 $G_2(1, 2’; 1’, 2)$，它包含了四个场算符。
• 第二项是单粒子格林函数乘以系统平均密度，即 $G(1, 1’) \times \langle \rho(2’,2) \rangle$。它代表了非关联部分，因为这是两个独立事件的乘积，二者没有发生任何相互作用。
• 如果我们希望知道外场 $U$ 对电子传播 $G$ 的纯关联贡献，必须减去这个非关联部分——剩下的才是关联 $L$：

$$L(1,2;1',2') = \frac{\delta G(1,1')}{\delta U(2',2)} \tag{197}$$

我们现在得到了双粒子的传播子定义，但是我们需要找到双粒子传播子的动力学方程，即 $L$ 满足什么动力学规律。

式 (197) 是 $L$ 的定义式，而动力学传播要从这里出发再推导。我们前面知道，格林函数和它的逆满足：

我们对上式两边对 $U$ 求导，然后右乘 $G$ 并移项可得：

它代表了：粒子先传播 $G$，然后受到一个微扰的变化 $\delta G^{-1}$，再继续传播 $G$ 的过程。

注意： 这里我们需要计算的是扰动让 $G$ 改变了多少，使用现在的基准来预测扰动后 $G$ 的变化 $\delta G$。这类似于一阶微扰：

$$\delta \psi_n = \sum_{m \neq n} \frac{\langle m | V | n \rangle}{E_n - E_m} \psi_m \tag{200}$$

我们求一阶微扰时使用的是没有微扰时候的能级。上式类似于式 (128) 中的第二项，仅保留了一阶微扰。

我们现在需要得到 $\frac{\delta G^{-1}}{\delta U}$ 这一项。回顾式 (118) Dyson 方程，$G^{-1}$ 包含无相互作用部分加上电子之间的 Hartree 项 $V_H$ 和自能项 $\Sigma$。引入外场 $U$ 后：

然后我们对 $U$ 求导：

其中 $-\frac{\delta U}{\delta U} = -\mathbb{I}$，而 $V_H$ 和 $\Sigma$ 都是通过 $G$ 间接依赖于 $U$ 的，因此可以写作：

我们把括号中前面那项定义为相互作用核 $K$。重新代回式 (199)：

其中：

• 第一项 $-G(-\mathbb{I})G = GG$，我们把它叫做无相互作用电子-空穴传播子 $L_0$。
• 第二项 $-G(-K \cdot L)G = GGKL$，写作 $L_0 K L$。

最终，双粒子传播子动力学方程为：

我们发现双粒子传播子方程很像单粒子 Dyson 方程（式 128），其中：

因此式 (205) 也可以做迭代展开：

它表明：双粒子的传播是两个分别与环境有无穷次耦合的粒子传播 $G$，再加上二者之间无穷次相互作用 $K$ 的卷积。

为了把 $K$ 的形式写出来，我们结合它的定义：

第一项 $\frac{\delta V_H}{\delta G}$ 源自 Hartree 势的变分。习惯上在 BSE 中把它称为交换核 $K^x$。这里需要小心术语：这个"交换"指的是 BSE 中的两种相互作用核之一（直接 vs. 交换），来源是 Hartree 势的变分，对应虚光子在电子-空穴对中的瞬时复合-重生过程；这跟 Hartree–Fock 中源于泡利不相容原理的"交换作用"是不同的概念[^d5]。

说明：

• 根据 Hartree 势定义：

$$V_H(1) = \int v(1, 2) \rho(2) \, d2 = -i \int v(1, 2) G(2, 2^+) \, d2 \tag{209}$$

即在点 $1$ 测得的所有其他位置 $2$ 处的电荷 $\rho(2)$ 对 $1$ 产生的库仑力的总和。

• 对其求导得到相互作用核 $K^x$：

$$K^x(1, 2; 3, 4) = \frac{\delta V_H(1)}{\delta G(3, 4)} \delta(1,2) \tag{210}$$

在上式中，位置 $1, 2$ 已经被作为观测点，因此需要使用 $3, 4$ 作为扰动源来求导。这个扰动对应于在 $3$ 处注入一个电子、在 $4$ 处注入一个空穴。

• 代入 Hartree 势的表达式：

$$\begin{aligned} K^x &= -i \int v(1, 2) \frac{\delta G(2, 2^+)}{\delta G(3, 4)} \, d2 \\ &= -i \int v(1, 2) \delta(2, 3) \delta(2, 4) \, d2 \\ &= -i \, v(1, 3) \, \delta(3, 4) \end{aligned} \tag{211}$$

• 上式包含 $\delta(3, 4)$，即只有当我们在同一个地方放一个电子-空穴对时才起作用。它说明：只有当我们注入的电子-空穴复合时，才会有一个非零的响应场 $v$ 发射出去，产生交换核作用。复合产生的能量通过 $v$ 扩散到位置 $1$。
• 量子场论中，没有超距作用，所有的力都是通过交换粒子来实现的。 因此可以认为这里的电子-空穴复合产生的是虚光子；$v$ 本质上可以认为是虚光子的传播子。“虚"是指这种光子仅参与中间过程的能量传递。因此 $v(1, 3)$ 描述电子-空穴对复合产生的虚光子从位置 $3$ 传播到位置 $1$。
• 我们可以将相互作用核写作 $K^x(1, 2; 3, 4)$，它代表在 $3, 4$ 处施加扰动后，$1, 2$ 处的响应。这样的矩阵是一个四维张量，为了方便计算，我们把输入端合并、输出端合并，写作二维矩阵：

$$K_{IJ} = \begin{pmatrix} K_{11} & K_{12} & \dots \\ K_{21} & K_{22} & \dots \\ \dots & \dots & \dots \end{pmatrix} \tag{213}$$

输入复合指标 $J$ 代表初始电子-空穴对状态 $(v’, c’)$，输出指标 $I$ 代表最终电子-空穴对状态 $(v, c)$。一个激子由一个价带空穴 $v$ 和一个导带电子 $c$ 组成。

• 对其在跃迁基底下积分：

$$\langle vc | K | v'c' \rangle = \int d1 \, d2 \, d3 \, d4 \, \Psi^*_{vc}(1, 2) \, K(1, 2; 3, 4) \, \Psi_{v'c'}(3, 4) \tag{214}$$

其中双粒子波函数 $\Psi_{vc}(1, 2) = \phi_c(1) \phi^*_v(2)$（空穴在 $2$，电子在 $1$）。

• 因此对 $K^x$ 可以写成：

$$\langle vc | K^x | v'c' \rangle = 2 \int d\mathbf{r} \, d\mathbf{r}' \, \phi_c^*(\mathbf{r}) \phi_v(\mathbf{r}) \, v(\mathbf{r} - \mathbf{r}') \, \phi_{c'}(\mathbf{r}') \phi_{v'}^*(\mathbf{r}') \tag{215}$$

（前置因子 2 来源于自旋单态求和。）

• 上式的形式很像跃迁偶极矩 $\langle c | \mathbf{r} | v \rangle$。
• 我们将左边的 $\phi_c^(\mathbf{r}) \phi_v(\mathbf{r})$ 叫做左偶极、右边的 $\phi_{c’}(\mathbf{r}’) \phi_{v’}^(\mathbf{r}’)$ 叫做右偶极。它的物理本质不再是带电粒子之间的静电相互作用，而是两个振荡偶极子之间通过电磁场进行的能量交换。我们把这种效应叫做"天线效应”，因为它具有长程性质。
• 这种交换作用所产生的效果是排斥效应，因为偶极振荡产生的感应极化场会阻止电子-空穴复合过程。

式 (208) 的第二项 $\frac{\delta \Sigma}{\delta G}$ 是激子结合能的来源。根据 GW 近似，电子自能为 $\Sigma(1, 2) = i G(1, 2) W(1, 2)$，因此这一项可以写作：

根据求导法则 $(uv)’ = u’v + uv’$，可以写出两项：

• 右边第一项根据 $\delta$ 函数性质为：

  $$\frac{\delta G(1, 2)}{\delta G(3, 4)} W(1, 2) = \delta(1, 3) \delta(2, 4) W(1, 2) \tag{218}$$

• 右边第二项 $\frac{\delta W}{\delta G}$ 我们设为 0——因为 $W$ 取决于介电函数 $\epsilon^{-1}$，而 $\epsilon^{-1}$ 取决于极化率 $P$，$P \sim GG$；这一项描述了"激子的形成反过来改变了晶体的屏蔽性质"。在 BSE 中我们假设激子导致的扰动不会改变整个材料的介电常数。

因此第二项直接相互作用核为：

这一项代表电子-空穴对之间的吸引作用：因为 $\delta(1, 3)$ 和 $\delta(2, 4)$ 分别"标记"了电子和空穴的位置，$W(1, 2)$ 代表屏蔽后的库仑相互作用（形式如 $\frac{1}{\epsilon r}$）。静电势能 $U \propto -|W|$ 因电子-空穴异号而为负——表示吸引势，降低体系能量。这一项就是我们常说的激子结合能的来源。

我们写出这一项在跃迁基组下的表达：

其中 $\phi_c^(\mathbf{r}) \phi_{c’}(\mathbf{r})$ 为导带电子在 $\mathbf r$ 处的密度重叠，$\phi_v(\mathbf{r}’) \phi_{v’}^(\mathbf{r}’)$[^d6] 为价带电子（空穴）在 $\mathbf r’$ 处的密度重叠。当 $c = c’$ 且 $v = v’$ 时（对角项）最大，代表静态的电子-空穴吸引。

最后，我们写出总的相互作用核：

• 第一项交换作用（$K^x$）代表激子-激子型的远程偶极耦合，给出排斥；
• 第二项直接作用（$K^d$）代表激子内电子-空穴的屏蔽库仑吸引，给出结合能。

要写出整体的 BSE 哈密顿量，价带电子跃迁到导带还需要跨越一个准粒子带隙：

因此，完整的 BSE 哈密顿量 $H^{BSE}$ 可以写作：

BSE 激子波函数

我们现在有了 BSE 的哈密顿量，要解出激子能量还需要激子态的波函数。

激子态 $|S\rangle$ 本质是在基态 $|0\rangle$ 上激发一个电子-空穴对：

其中 $\hat a_c^\dagger$ 是导带轨道 $c$ 上的电子产生算符，$\hat a_v$ 是价带轨道 $v$ 上的电子湮灭算符（即等效地在价带中创造空穴）。

为了解出 $A_{vc}$，我们引入多体哈密顿量本征方程：

其中 $\hat{H}$ 是二次量子化的总哈密顿量算符，$E_0$ 是基态能量，$\Omega_S$ 是激子激发能量（小于带隙），即 $\Omega_S = E_{\text{gap}} - E_b$（能隙减去激子结合能）。

说明： 在二次量子化中，总能量算符的标准定义是把单粒子哈密顿量密度积分：

$$\hat{H}_0 = \int d\mathbf{r} \, \hat{\Psi}^\dagger(\mathbf{r}) \, \hat{h}(\mathbf{r}) \, \hat{\Psi}(\mathbf{r}) \tag{226}$$

• 代入场算符的定义式 (58)（59）：

$$\hat{H}_0 = \int d\mathbf{r} \left[ \sum_j \phi_j^*(\mathbf{r}) \hat{a}_j^\dagger \right] \hat{h}(\mathbf{r}) \left[ \sum_i \phi_i(\mathbf{r}) \hat{a}_i \right] \tag{227}$$

• 把求和和算符 $\hat a^\dagger, \hat a$ 提到积分外（与空间坐标 $\mathbf r$ 无关）：

$$\hat{H}_0 = \sum_{j, i} \hat{a}_j^\dagger \hat{a}_i \int d\mathbf{r} \, \phi_j^*(\mathbf{r}) \, \hat{h}(\mathbf{r}) \, \phi_i(\mathbf{r}) \tag{228}$$

• 因为 $\phi_i$ 是本征态，所以 $\hat h \phi_i = \epsilon_i \phi_i$，代入可得：

$$\hat{H}_0 = \sum_i \epsilon_i \, \hat{a}_i^\dagger \hat{a}_i \tag{229}$$

说明： 我们一般说本征值方程，其形式为：

$$\mathbf{A} \mathbf x = \lambda \mathbf x \tag{230}$$

• 举个例子[^d7]，假设算符 $\mathbf A$ 可以是：

$$\mathbf{A} = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \tag{231}$$

• 现在有一个向量 $\mathbf{x} = \begin{pmatrix} 1 \ 1 \end{pmatrix}$，验证 $\mathbf x$ 是不是 $\mathbf A$ 的本征向量——即 $\mathbf x$ 在算符 $\mathbf A$ 的作用下保持方向不变、只被整体放大或缩小。

$$\mathbf{A}\mathbf{x} = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 3 \\ 3 \end{pmatrix} \tag{232}$$

• 输出结果是原向量[^d8] $\mathbf x$ 的三倍，即：

$$\lambda \mathbf{x} = 3 \times \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 3 \\ 3 \end{pmatrix} \tag{233}$$

• 因此式 (230) 成立。它说明：在经过一个物理过程后，输出仅是输入的一个缩放，其方向（波函数）没有发生改变。

为了解出 $A_{vc}$，我们用一个探测态 $\langle 0 | \hat{a}{v’}^\dagger \hat{a}{c’}$ 左乘式 (225)[^d9]：

将 $|S\rangle$ 的展开式代进去：

其中，后面的 $\langle 0 | \hat{a}{v’}^\dagger \hat{a}{c’} , \hat{a}c^\dagger \hat{a}v | 0\rangle = \delta{vv’}\delta{cc’}$。

并且我们在前面也推导出：

即相互作用核是交换核和直接核之和。BSE 哈密顿量再加入动能项（准粒子带隙）后，可写作：

第一项是"动能"（准粒子色散）项：$\epsilon_c(\mathbf{k}) - \epsilon_v(\mathbf{k})$ 已经是 GW 修正过的准粒子能差——既包含了 $E_g$（带边能差），又包含了对 $\mathbf k$ 的色散（动能修正）。有时为了方便简写为 $(\epsilon_c - \epsilon_v)$。

因此，重新代入式 (235)，可得：

这就是 BSE 本征方程。

现在上面式子变成了与式 (230) 一致的本征方程形式 $H \cdot \vec{A} = \Omega \cdot \vec{A}$。其中，左边的 $A_{vc}$ 是输入（所有可能的跃迁对混合在一起），右边是特定的跃迁对 $A_{v’c’}$。在数值计算中，调用线性代数库（如 LAPACK）找到满足上式的 $\Omega_S$ 和 $A_{v’c’}$。

接下来我们可以求出波函数的形式，写成基组展开：

其中 $|v\mathbf{k} \to c\mathbf{k}\rangle$ 是单次跃迁，代表在某个动量点 $\mathbf k$ 把电子从价带激发到导带；$A_{vc\mathbf{k}}^S$ 代表某个特定跃迁的贡献——它的平方即为该激子在 $v\mathbf{k} \to c\mathbf{k}$ 通道的占比；然后对所有 $\mathbf k$ 点和所有能带求和。

也可以把波函数写成实空间的坐标形式：

它表明：激子在实空间的形状是由多个导带电子波函数 $\phi_{c\mathbf{k}}(\mathbf{r}e)$ 和价带空穴波函数 $\phi{v\mathbf{k}}^*(\mathbf{r}h)$ 的乘积，按照 $A{vc\mathbf{k}}^S$ 权重系数比例混合而成。

Casida 方程

场算符展开

在前面的式 (68) 中，我们用场算符表达了电子密度算符。现在考虑自旋 $\sigma$，密度算符可以写作：

现在我们使用正交归一的、定态 KS 轨道 ${\phi_{p\sigma}(\mathbf{r})}$ 作为基底，将场算符展开：

上式是湮灭场算符，取厄米共轭得到产生场算符：

其中：

• $p, q$ 是 KS 轨道指标，遍历所有 KS 轨道（占据态和非占据态均包括在内）[^e1]；
• $\hat{a}{p\sigma}$ 和 $\hat{a}{q\sigma}^\dagger$ 分别是作用在第 $p$、$q$ 个 KS 轨道上的湮灭算符和产生算符；
• $\phi_{q\sigma}^*(\mathbf{r})$ 和 $\phi_{p\sigma}(\mathbf{r})$ 是标量函数（空间坐标确定时是确定的复数）。

我们将用 KS 轨道表达的场算符代入式 (241)，可得：

将求和符号合并整理：

如果要在某个量子态 $|\Phi(t)\rangle$ 下获得可观测电子密度 $\rho(\mathbf{r}, t)$，需要取期望值：

对于左右矢而言，产生/湮灭算符会改变粒子数，因此不能提出左右矢；而两个标量函数则与态空间无关。需要注意：这里的左右矢是 Fock（粒子数占据）空间的态，并非位置表象的态：

$n_i$ 表示第 $i$ 个 KS 轨道上的电子占据数（对费米子为 0 或 1）。

我们把期望值 $\langle \Phi(t) | \hat{a}{q\sigma}^\dagger \hat{a}{p\sigma} | \Phi(t) \rangle$ 定义为单粒子约化密度矩阵在 KS 轨道基底下的矩阵元，记作 $P_{pq\sigma}(t)$：

在上式中：

• $\phi_{q\sigma}^*(\mathbf{r}) \phi_{p\sigma}(\mathbf{r})$ 表达了轨道在空间中的重叠与干涉；
• $P_{pq\sigma}(t)$ 表示了相应轨道乘积 $\phi_{q\sigma}^*(\mathbf{r}) \phi_{p\sigma}(\mathbf{r})$ 在 $t$ 时刻的权重。

我们将空间分布属性和时间演化解耦——所有的含时演化都包含在 $P_{pq\sigma}(t)$ 中。最后，电子密度可以写作：

二次量子化下的时间演化密度矩阵

首先回顾一下密度矩阵。在 FSSH 的推导中曾对密度矩阵下过定义：对于双态系统 $|\psi\rangle = c_1|1\rangle + c_2|2\rangle$，密度矩阵的矩阵元为 $\rho_{kj} = c_k c_j^*$，整体的密度矩阵为[^e2]：

现在我们需要一个更一般的定义。从现在开始我们做如下约定：

• $\rho(\mathbf{r})$ 表示空间坐标上的真实电子密度；
• $\hat{\rho}$ 表示密度算符；
• $\mathbf{P}$ 表示单粒子约化密度矩阵。

那么更一般的密度算符 $\hat{\rho}$ 可以写作：

我们可以通过这个更一般的定义式推导出式 (250)。

说明：

• 假设存在一个混合态系综：$|\Psi\rangle$ 以概率 $p_A$ 处于态 $|\Psi_A\rangle$，以概率 $p_B$ 处于态 $|\Psi_B\rangle$。
• 此时密度算符可以写作：

$$\hat{\rho} = p_A |\Psi_A\rangle\langle\Psi_A| + p_B |\Psi_B\rangle\langle\Psi_B| = \sum_k p_k |\Psi_k\rangle\langle\Psi_k| \tag{252a}$$

• 然后我们把每个态 $|\Psi_k\rangle$ 用一组正交归一的基底 ${|1\rangle, |2\rangle, \dots}$ 展开：

$$|\Psi_k\rangle = c_{k1} |1\rangle + c_{k2} |2\rangle = \sum_i c_{ki} |i\rangle \tag{252b}$$

• 在量子力学中，任何一个算符 $\hat{A}$ 在基底 ${|i\rangle}$ 下的矩阵元定义为 $A_{ij} = \langle i | \hat{A} | j \rangle$，其中 $i$ 是行指标、$j$ 是列指标。因此密度算符的矩阵元可以写作：

$$\rho_{ij} = \langle i | \hat{\rho} | j \rangle \tag{252c}$$

• 代入密度算符的定义式 (251)，可以得到：

$$\rho_{ij} = \sum_k p_k \langle i | \Psi_k\rangle \langle \Psi_k | j \rangle \tag{253}$$

• 其中 $\langle i | \Psi_k\rangle = c_{ki}$ 是展开系数；对 $\langle \Psi_k | j \rangle$ 求复共轭得到 $c_{kj}^*$。因此密度矩阵元为：

$$\rho_{ij} = \sum_k p_k \, c_{ki} \, c_{kj}^* \tag{254}$$

• 排成矩阵形式（以 2×2 情形为例）：

$$\rho = \begin{pmatrix} \rho_{11} & \rho_{12} \\ \rho_{21} & \rho_{22} \end{pmatrix} = \begin{pmatrix} \sum_k p_k \, c_{k1} c_{k1}^* & \sum_k p_k \, c_{k1} c_{k2}^* \\ \sum_k p_k \, c_{k2} c_{k1}^* & \sum_k p_k \, c_{k2} c_{k2}^* \end{pmatrix} \tag{255}$$

• 这是完整的混合态密度矩阵（同时包含量子叠加和经典混合）。

• 纯态特例： 如果系统是纯态，即只有单个 $k$（$p_1 = 1$、其余为 0），求和消失，演化为以下简单形式[^e3]：

$$\rho = \begin{pmatrix} c_1 c_1^* & c_1 c_2^* \\ c_2 c_1^* & c_2 c_2^* \end{pmatrix} = \begin{pmatrix} |c_1|^2 & c_1 c_2^* \\ c_2 c_1^* & |c_2|^2 \end{pmatrix} \tag{256}$$

• 注意： 退相干并不代表系统塌缩到了纯态，而只是指在[^e4]密度矩阵中非对角元降低到 0。

• 例子： 假设存在一个叠加态：

$$|\Psi\rangle = \frac{1}{\sqrt{2}}|0\rangle + \frac{1}{\sqrt{2}}|1\rangle \tag{257}$$

• 根据式 (256) 可以直接写出它的密度矩阵：

$$\rho = \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix} \tag{258}$$

• 非对角元 0.5 表明两个态之间存在显著的量子相干性。

• 我们使用 $\mathrm{Tr}(\rho^2)$ 来检验一个态是否是纯态——纯态 $\mathrm{Tr}(\rho^2) = 1$，混合态 $\mathrm{Tr}(\rho^2) < 1$。验证当前叠加态：

$$\rho^2 = \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix} \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix} = \begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix} \tag{259}$$$$\mathrm{Tr}(\rho^2) = 0.5 + 0.5 = 1 \tag{260}$$

• 退相干后，两个态之间的相干性消失（非对角元为 0），此时：

$$\rho_{\text{mix}} = \begin{pmatrix} 0.5 & 0 \\ 0 & 0.5 \end{pmatrix} \tag{261}$$

• 此时的纯度为：

$$\rho_{\text{mix}}^2 = \begin{pmatrix} 0.5 & 0 \\ 0 & 0.5 \end{pmatrix} \begin{pmatrix} 0.5 & 0 \\ 0 & 0.5 \end{pmatrix} = \begin{pmatrix} 0.25 & 0 \\ 0 & 0.25 \end{pmatrix} \tag{262}$$$$\mathrm{Tr}(\rho_{\text{mix}}^2) = 0.25 + 0.25 = 0.5 \tag{263}$$

• 此时已变为混合态——失去了量子相干性，只保留了"50% 概率在态 $|0\rangle$、50% 概率在态 $|1\rangle$“的经典概率。

• 如果对混合态进行一次测量并测到了态 $|0\rangle$，那么密度矩阵塌缩为：

$$\rho_{\text{post}} = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \tag{264}$$

• 再检验纯度：

$$\mathrm{Tr}(\rho_{\text{post}}^2) = 1 + 0 = 1 \tag{265}$$

• 因为观测的"投影"作用，态变成了 100% 确定在[^e5]态 $|0\rangle$ 上的纯态。

回到含时演化推导。假设一个纯态体系的含时薛定谔方程为（原子单位制 $\hbar = 1$）：

对应的左矢方程为：

前面得到了密度算符（纯态下）$\hat{\mathbf P} = |\Psi\rangle\langle\Psi|$。现在对其求时间导数（链式法则）：

然后代入式 (266) 和 (267)：

上式就是密度算符与哈密顿量的对易关系。引入对易子 $[\hat{H}, \hat{\mathbf{P}}] \equiv \hat{H}\hat{\mathbf{P}} - \hat{\mathbf{P}}\hat{H}$，因此可以写作：

这就是 Liouville–von Neumann 方程——它是 Casida 方程推导的出发点。下一步应当对 $\mathbf P(t) = \mathbf P^{(0)} + \delta\mathbf P(t)$ 做线性响应展开，并引入 KS 哈密顿量的扰动。

关于本章的整体提醒： 本章是讲义的未完结部分，缺少 Casida 方程的核心推导（线性响应展开、KS 扰动哈密顿量、Hartree + $f_{xc}$ 耦合核、最终的 $\begin{pmatrix}A&B\B^&A^\end{pmatrix}$ 矩阵形式）。如需补全可参考：Casida, M. E. Recent Advances in Density Functional Methods, Part I（1995）中第 5 章 “Time-Dependent Density Functional Response Theory”，这是 Casida 方程的标准引文。

附录

参考文章