固体物理中的色散关系

郭帅斐 2013301020066 物基一班

1 摘要：

总得来说，固体物理研究的是周期性势场中波的传播问题，包括晶格振动产生的格波，电子传播的布洛赫波，还有固体作为介质而在其内存在的电磁波等。量子力学告诉我们，描述一个波我们需要知道它的两个量子数：能量和动量，在布洛赫波和格波的情形下同样成立，只是动量是所谓的由于晶格周期性而产生的准动量。一个动量能量二维数就能描述一个波动的本征态，色散关系是一系列这样的二维数构成的点的集合，知道物质内某种波的色散关系，就相当于知道了这种波的全部信息，也能反映物质结构的信息。比如，知道晶体内布洛赫波的色散关系（能带结构），我们就能判断这种晶体是导体半导体还是绝缘体；同时我们根据能带结构计算群速度，有效质量等物理量，进而计算电导率等实验上可观测的非常重要的物理量。鉴于此，如果能编写程序计算晶体的格波色散关系（声子谱）和布洛赫波色散关系（能带结构）的话，将大大加深对于这些概念的理解。作者分别编写了计算一维双原子链，二维双原子链的声子谱的程序，同时编写了用平面波展开法计算能带结构的简单程序，并且得到一些初步结果。同时，我们还可能关心波函数的信息，因此一个计算一维薛定谔方程的程序也被编写了，可以展示一些特殊势的波函数信息。

2 正文：

2.1 算法介绍：

固体物理问题求解的麻烦在于体系很大，粒子很多，在阿伏伽德罗常数($10^{23}$)的量级，需要求解很大的方程组。实践中一般根据一些近似的原理，编写成熟的程序进行计算，这里我采用普通计算机可以实现的计算量来编写程序，因此格波计算只能计算到两维，因为三维的话至少需要将一个$30^3\times 30^3$ 的矩阵对角化，即使是python里现成的计算包也会把我心爱的电脑算死掉。是的，整个程序的核心就是矩阵对角化，下面我们来了解一下其中的算法。下列算法都归结为一类本征值方程问题，因而具有相似之处。

• 首先是格波色散关系的计算，我们只记及最近邻原子之间的相互作用（在编写程序中并不是必须的，因为我们可以很容易的写出次近邻原子的作用矩阵元），同时考虑简谐近似（这个近似对程序有最大的简化作用，因为考虑高阶项的话，算法就完全变了），最后应用周期性边界条件（这在程序中很容易实现，只需在矩阵顶角放入矩阵元）。至此，我们可以讲格波问题抽象为一个矩阵问题。利用晶格振动系统的哈密顿量：
  

  $$H=\frac{1}{2}\sum \dot{q}_i+\frac{1}{2}\sum \lambda_{ij}q_iq_j$$ 由正则方程得到$3N$个耦合的振动方程：
  $$\ddot{q}_i+\sum\lambda_{ij}q_j=0$$
  根据上面我们所用的近似办法，我们可以得到体系的动力学矩阵$\lambda_{ij}$的形式，将此矩阵对角话，得到简正坐标的方程：
  $$\ddot{Q}_i=\omega_i^2Q_i$$
  上述过程相当于：
  $$\frac{mw^2}{K}\vec{A}=D\cdot\vec{A}$$
  若 $m$ 也是矩阵，则只需运行相应的矩阵逆运算即可。在理论上求解时，我们一般将矩阵元对倒格矢做傅里叶变换，从而得到关于每个波矢 $\vec{q}$ 的动力学矩阵只有 $3n\times 3n$ 个矩阵元。但是在数字计算时，只要系统不太大，我们可以直接求解整个矩阵方程，得到对应于每个模式的频率 $w$ 。在python-numpy中，有把矩阵对角化的模块，直接运用可以大大节省运算量降低程序的复杂性。

• 其次是能带结构的计算，根据布洛赫定理，周期场中的波函数具有如下的形式：
  

  $$\psi_{\vec{k}}(\vec{r})=u_{\vec{k}}(\vec{r})e^{i\vec{k}\cdot\vec{r}}$$
  其中 $u_{\vec{k}}(\vec{r})$ 具有晶格的周期性，可按倒格矢展开。利用波函数具有倒格子的周期性，最终得到布洛赫波函数可表示为：
  $$\psi_{\vec{k}}(\vec{r})=\frac{1}{\sqrt{N\Omega}}\sum _{\vec{K}_h}a_{\vec{k}+\vec{K}_h}e^{i({\vec{k}+\vec{K}_h})\cdot\vec{r}}$$
  将此波函数代入定态薛定谔方程 $H\psi_{\vec{k}}(\vec{r})=E_{\vec{k}}\psi_{\vec{k}}(\vec{r})$ ，并且利用势能的晶格周期性，将其按倒格矢展开，得到势能的矩阵元为：
  $$\langle \vec{k}+\vec{K}_h|\vec{k}+\vec{K}_{h'}\rangle=\frac{1}{N\Omega}\int V(\vec{r})e^{-i({\vec{K}_h-\vec{K}_{h'}})\cdot\vec{r}}d\vec{r}$$
  这是一个傅里叶变换的形式，同样，利用 python-numpy 自带的 fft 模块可以轻松得到此矩阵元。最终得到矩阵方程，也是我们要求解的形式：
  $$D\cdot \vec{A}=E_{\vec{k}}\vec{A}$$
  其中 $D$ 的矩阵元为：当$\vec{K}_h=\vec{K}_{h'}$ 时，有 $D_{\vec{K}_h,\vec{K}_{h'}}=\frac{\hbar}{2m}(\vec{k}+\vec{K}_h)^2$ ; 当 $\vec{K}_h\ne\vec{K}_{h'}$ 时，有 $D_{\vec{K}_h,\vec{K}_{h'}}=V(\vec{K}_h - \vec{K}_{h'})$ ，这里 $\vec{K}_h$ 是一个无穷角标，实际运算中，一般只取有限个倒格矢来计算，我选择的数量是500个平面波，每次运行大约需要两分钟。有实际计算表明，计算 Cu 的能带结构得到收敛的结果，大约需要6000个平面波！我们的程序只得到了初步的结果。分别将一系列模式$\vec{k}$ 代入方程，得到一系列的矩阵，分别得到相应的能量。至此我们得到所需要对角化的矩阵，下面就是运行程序了。实践中使用的有缀加平面波方法和赝势方法，以及第一性原理的从头算法，都有相应的成熟的软件，对于我来说，肯定是写不出来的。

• 最后，我们可能也对波函数的具体形式感兴趣，下面介绍矩阵方法求解定态薛定谔方程。事实上，我们在前面求解 laplace 方程的时候用到了差分后 laplace 算子的矩阵形式，再加上势函数本身在位置表象中就是对角化的。我们得到哈密顿量的矩阵表达式：
  

  $$\begin{equation} \left( \begin{array}{ccc} 2+V(0,0) & -1 & 0 & \cdots \\ -1 & 2+V(0,1) & -1 & \cdots \\ \cdots &\cdots & \cdots & \cdots \\ \end{array} \right) \end{equation}$$
  根据特征值方程 $H\psi=E\psi$ ，解此矩阵得到特征值即为能量，相应的本征矢即为对应于该能量波函数。

2.2 程序及运行结果

• 关于格波色散关系的结果：
  首先运行一维双原子链的模块 Lattice ，运行程序得到一系列的本征值，进而得到其色散关系如下所示，图1和图2分别是双原子链的两种原子质量相等和不相等时的情形：
  
  

  $$图1,当双原子链两种原子的质量不相同时，两支频带之间出现频隙$$
  
  $$图2，当双原子链两种原子的质量相同时，两支频带之间频隙消失$$
  我们可以看到声频支格波是较低的那支，其在长波情形下近似为一条直线，相当于连续介质情形，因为其相邻原子之间距相比波长很小，近似是连续的，与声波相似。而较高的那支因其可与光波耦合而得名，光波相当于色散关系中斜率很大的直线。
  其次运行不同方向的双原子链相互耦合的情形，得到更接近于真实三维情形的色散关系（声子谱），下图3所示是两个方向的声子谱：
  
  $$图3，真实晶体中两个不同方向的色散关系$$
  最后是二维原子链的情形，我们考虑其在垂直与表面的方向的运动，在平行与表面的方向的情形与此类似（受相互作用情况相同，动力学矩阵形式一样）。我们首先考虑原子相同的情形，运行程序中的two_dlattice 模块， 得到色散关系如下图4所示：
  
  $$图4，二维原子链中所有原子相同的时的色散关系，左右分别是直接绘制的模式图与将w排序后画出的图形$$
  对于对角化后的到的数据，我们不仅能得到关于色散关系的信息，同时，特征向量能显示出各模式具体的振动情况，现画出其中的三个模式的具体形貌如下图5所示：
  
  $$图5，不同模式的具体形貌，z坐标表示振动幅度相对大小$$
  当我们改变两种原子的质量，使相邻原子的质量分别为1和2时，我们发现总的色散关系变得更加弯曲，但并未出现频隙，其结构如下图6 所示：
  
  $$图6，当相邻原子质量比为1:2时，色散曲线变得更加弯曲$$
  当相邻原子质量分别为1和8时，我们明显看到频隙出现，说明晶格振动中不可能出现相应的频率模式，其色散曲线如下图7所示：
  
  $$图7，当相邻原子质量比为1:8时，色散曲线中出现频隙$$
  至此，通过观察不同维度，不同原子种类的色散曲线，我们发现了出现带隙的原因，真正得观察到了格波的色散关系曲线。

• 关于能带结构的运行结果：
  首先需要知道，虽然用500个平面波来做线性组合，已经用去了电脑的100%的CPU，但是不敢保证结果是收敛的，只能得到关于色散关系的初步的信息。我们用500个平面波叠加，应能算出500个能量值，这里（实际上在实践中）我们只关心前两个能带，因为这两个能带是可以有电子填充的，与物质性质相关性更强的。
  当我们选取的吸引势很小时，也就是动力学矩阵的非对角元绝对值很小时，我们得到增长的色散曲线如下图8所示：
  
  

  $$图8，当吸引势较小时，得到增长的能带结构$$
  当选取的势函数大小合适时，可能会出现下面的带递减，上面的带递增的能带结构，结果如下图9所示：
  
  $$图9，当吸引势大小合适时，得到对称的能带结构$$
  当吸引势较大时，则会出现前两条能带是下降的结构，其结果如下图10所示：
  
  $$图10，当吸引势较大时，得到下降的能带结构$$
  我们只是通过改变势函数的大小，得到了能带的变化趋势，并没有得到与固体物理书本上一样的图形，这是这个程序的遗憾之处。

• 关于一维定态波函数的计算
  通过矩阵对角化的办法，我们计算了一维谐振子的不同能量本征值的波函数形式，绘制图像如下图11所示：
  
  

  $$图11，一维谐振子的不同能量本征值的波函数$$
  原则上，用这种方法可以得到所有的一维势的定态波函数，比如一维无限深势阱及氢原子的径向方程等，运行程序（更改势函数的形式）后，得到相应与氢原子径向方程的解，也就是 0阶 Bessel 函数，得到如下图12所示的定态波函数：
  
  $$图12，氢原子径向波函数（l=0）$$

结论与致谢

• 结论：
  1.首先，我们运用程序Lattice.py计算了一维双原子链的色散关系，同时计算出实际晶体在不同方向的色散关系形式，并得到色散关系中能隙出现的原因在于两种原子质量不一样的结论。计算二维原子链的色散关系，得到同样的结论。对于同类原子情形，我们画出了每种模式中的相对振幅与位置的关系，得到漂亮的曲线，说明不同模式是不同驻波叠加。
  2.其次，我们运行程序band.py计算了不同势函数情形下的色散关系，发现了色散关系趋势与吸引势强弱的关系：吸引势越强，色散曲线下降；吸引势越弱，色散曲线上升。同时验证了平面波方法的收敛性很差，需要大量的平面波才能得到正确的解。
  3.最后，通过运行程序wavefun.py分别计算了一维谐振子和氢原子径向的波函数形式，原则上我们可以用这个程序计算所有的一维定态薛定谔方程，至此，我们不仅能得到关于量子态的能量-动量色散关系，还能得到相应于每个态的波函数具体形式。

• 致谢：
  首先感谢蔡浩老师课上的教导，感谢李芳莹陈锋等同学的交流与帮助，感谢范登栋学长关于波函数计算程序的指导。计算物理结课了，感觉还是学到了很多东西，以后再不怕写程序了。

参考文献

【1】计算物理，清华大学出版社
【2】Computational Physics, J.M.Thijssen.世界图书出版社
【3】Matlab中的数值计算方法，清华大学出版社
【4】固体物理学，胡安. 高等教育出版社
【5】量子力学概论，David J.Griffiths. 机械工业出版社