硅晶体原子间相互作用力常数的计算与负热膨胀机制的研究

Calculations of the Interatomic Force Constants and Study on the Mechanism of Negative Thermal Expansion of Silicon Crystal

摘要: 本文根据Rignanese等人基于第一性原理得到的硅晶体中原子间的力常数矩阵元，计算了两体和三体线性力常数，再将这些线性力常数对原子间距求导数得两体和三体非线性力常数，在此基础上运用硅晶体的热膨胀系数公式计算了其热膨胀系数。计算结果与实验结果很好地吻合，这表明硅晶体的各个线性和非线性力常数是正确的。计算结果还表明，硅晶体的两体非线性力常数为负，两体势引起了正热膨胀，而三体非线性力常数为正，三体势引起负热膨胀，且低温时负热膨胀效应大于正热膨胀效应，因而总体上呈现低温负热膨性质。

Abstract: The two-body and three-body linear force constants in silicon crystal were calculated based on the interatomic force constant matrix elements obtained by Rignanese with the ab initio method, and then the two-body and three-body non-linear force constants were obtained by derivate the corresponding linear force constants with respect to bond length. Finally, the thermal expansion coefficients of silicon crystal were calculated based on these force constants and formula for thermal expansion coefficients of silicon crystal, and the calculated results are in good agreement with experimental results, it means that the results of all the linear and non-linear force constants are correct. It is also found that the thermal expansion caused by two-body potential is positive because of the negative two-body non-linear force constant, the thermal expansion caused by three-body potential is negative because of the positive three-body non-linear force constant, and at low temperature the total thermal expansion is negative because the absolute value of thermal expansion caused by three-body potential is greater than thermal expansion caused by two-body potential.

1. 引言

硅单晶被广泛运用于MEMS制造，研究它的热学性质对于MEMS的正确设计和可靠使用是非常重要的。自Erfling发现低温下硅晶体低温负热膨胀性质以来 [1] ，硅晶体这种奇异的热学性质一直是研究热点 [2] ，人们一直试图搞清其物理机制 [3] [4] [5] [6] [7] 。我们曾基于微扰论推导了硅晶体的热膨胀系数公式 [8] ，发现硅晶体的热膨胀由原子间二体和三体的非线性相互作用所引起，而采用常被用于模拟硅单晶热学性质的原子间相互作用势的Stillinger-Weber模型 [8] [9] [10] ，计算原子间二体和三体的线性和非线性力常数，并计算硅晶体的热膨胀系数，发现原子两体和三体相互作用都引起正热膨胀，因而不能呈现低温负热膨胀，这表明由该模型得到的原子间的力常数并不准确。然而，准确的力常数对于计算硅晶体的热学性质来说非常重要。我们曾推测了以下的硅晶体低温负热膨胀的物理机制 [8] ：为负值的两体非线性力常数引起正热膨胀，而为正值的三体非线性力常数引起负的热膨胀，并且低温时三体势引起的负热膨胀强于两体势引起的正热膨胀，从而使硅晶体呈现低温负热膨胀性质。因为没有由第一性原理得到的正确的硅晶体原子间的线性和非线性相互作用力常数的数据，因而我们还一直没有对这一推测进行验证。本文将根据Rignanese等人 [5] 基于第一性原理得到的硅晶体中原子间的力常数矩阵元，计算两体和三体的线性和非线性力常数，再利用热膨胀系数公式计算两体及三体相互作用对热膨胀系数的贡献，计算硅晶体的热膨胀系数并与实验结果比较，从而验证硅晶体低温负热膨胀的物理机制，同时也为以后通过晶格动力学研究硅晶体的热学性质提供准确可靠的两体和三体的线性和非线性力常数。

2. 硅晶体晶格动力学与热膨胀系数公式

硅晶体中键长、静态键长及其变化量分别用 $R$ 、 $R_0$ 和 $r$ 表示，键角、静态键角及其变化量分别用 $\Theta$ 、 ${\Theta }_0$ 和 $\theta$ 表示。求解硅晶体的晶格动力学问题，可得声子频谱 ${\omega }_{kj}$ 及单位本征矢 $e\left(\chi |kj\right)$ ，其中， $k=\left(k_x,k_y,k_z\right)$ ， $j=1,2,\cdots ,6$ 。设原子 $A$ 、 $O$ 和 $B$ 的平衡位置 $R_A$ 、 $R_O$ 和 $R_B$ 分别为 $\left(l_x,l_y,l_z\right)a/2$ 、 $R_A+\left(1,1,1\right)a/4$ 和 $R_A+\left(0,1,1\right)a/2$ ，形成键角 $AOB$ ，则键长 $AO$ 和键角 $AOB$ 的变化量可分别表示为

$r_{AO}\left(l\right)={\displaystyle \underset{k}{\sum }{\eta }_{AO}\left(k\right)\left(a_{kj}+a_{-kj}^{+}\right)e^{ik\cdot R\left(l\right)}}$ (1)

${\theta }_{AOB}\left(l\right)={\displaystyle \underset{k}{\sum }{\eta }_{AOB}\left(kj\right)\left(a_{kj}+a_{-kj}^{+}\right)e^{ik\cdot R\left(l\right)}}$ (2)

其中， $a_{kj}$ 和 $a_{-kj}^{+}$ 分别是声子的消灭算符和产生算符， $l=\left(l_x,l_y,l_z\right)$ 。

${\eta }_{AO}\left(kj\right)=\sqrt{\frac{\hslash }{2Nm{\omega }_{kj}}}\left[e\left(O|kj\right)-e\left(A|kj\right)\right]\cdot I_{AO}$ (3)

${\eta }_{AOB}\left(kj\right)=\frac{2}{3a}\sqrt{\frac{\hslash }{Nm{\omega }_j\left(k\right)}}\left\{\left[\left(2,1,1\right)e^{\frac{ia}{2}\left(k_y+k_z\right)}+\left(2,-1,-,1\right)\right]\cdot e\left(A|kj\right)-4e_x\left(O|kj\right)\right\}$ (4)

其中， $I_{AO}$ 为 $AO$ 方向的单位矢量。

硅原子间相互作用势可分为两体势 $U_2\left(R\right)$ 和三体势 $U_3\left(R_1,R_2,\Theta \right)$ 两个部分。定义两体和三体线性力常数 $f_r$ 和 $f_{\theta }$ 分别为 ${\partial }^2{U}_2/\partial R^2$ 和 ${\partial }^2{U}_3/R_0^2{\partial }^2\Theta$ ，两体和三体线非线性力常数 $g_r$ 和 $g_{\theta }$ 分别为 $\partial f_r/\partial R$ 和 $\partial f_{\theta }/\partial R_{1\left(2\right)}$ 。我们运用量子力学微扰理论推导了硅晶体热膨胀系数公式 [8]

$\gamma =-\frac{8\sqrt{3}\varphi g_r}{3}{\displaystyle \underset{kj}{\sum }{\left|{\eta }_{OA}\left(kj\right)\right|}^2\frac{\partial {\bar{n}}_{kj}}{a\partial T}}-\frac{48\sqrt{3}\varphi g_{\theta }}{3}{\displaystyle \underset{kj}{\sum }{l}^2{\left|{\eta }_{AOB}\left(kj\right)\right|}^2\frac{\partial {\bar{n}}_{kj}}{a\partial T}}$ (5)

其中，右边第一和第二项分别是两体和三体势产生的热膨胀系数 ${\gamma }_2$ 和 ${\gamma }_3$ ， $\varphi$ 为正，可由下式计算

$\varphi ={\displaystyle \underset{k\text{'}j\text{'}}{\sum }\delta \left(k\text{'}\right)\frac{{\left|{\eta }_{OA}\left(k\text{'}j\text{'}\right)\right|}^2}{\hslash {\omega }_{k\text{'}j\text{'}}}}$ (6)

3. 硅晶体原子间二体和三体力常数计算公式

与OA键有关的晶格相互作用势可表示为两体势 $U_2\left(OA\right)$ 和三体势 $U_3\left(OA\right)$ 之和：

$U\left(OA\right)=U_2\left(OA\right)+U_3\left(OA\right)$ (7)

将 $U\left(OA\right)$ 在原子平衡位置处进行Taylor级数展开，其二次项为

$U^{\left(2\right)}\left(OA\right)=\frac{1}{2}{f}_r{r}_{OA}^2+\frac{1}{2}{\rho }^2{f}_{\theta }\left({\theta }_{AOB}^2+{\theta }_{AOC}^2+{\theta }_{AOD}^2+{\theta }_{OA{B}^{\prime }}^2+{\theta }_{OA{C}^{\prime }}^2+{\theta }_{OA{D}^{\prime }}^2\right)$ (8)

其中，原子 $B$ 、 $C$ 、 $D$ 、 $B^{\prime }$ 、 $C^{\prime }$ 和 $D^{\prime }$ 的平衡位置坐标分别为： $R_A+\left(0,1,1\right)a/2$ ， $R_A+\left(1,1,0\right)a/2$ ， $R_A+\left(1,0,1\right)a/2$ ， $R_A+\left(1,-1,-1\right)a/4$ ， $R_A+\left(-1,-1,1\right)a/4$ ， $R_A+\left(-1,1,-1\right)a/4$ ，键角的变化量可分别表示为

${\theta }_{AOB}=\frac{2\sqrt{2}}{3a}\left(u_B^y-u_A^y+u_B^z-u_A^z\right)+\frac{4\sqrt{2}}{3a}\left(u_A^x+u_B^x-2u_O^x\right)$ (9)

${\theta }_{OA{B}^{\prime }}=\frac{2\sqrt{2}}{3a}\left(u_O^y-u_{B^{\prime }}^y+u_O^z-u_{B^{\prime }}^z\right)-\frac{4\sqrt{2}}{3a}\left(u_O^x+u_{B^{\prime }}^x-2u_A^x\right)$ (10)

${\theta }_{AOC}=\frac{2\sqrt{2}}{3a}\left(u_C^x-u_A^x+u_C^y-u_A^y\right)+\frac{4\sqrt{2}}{3a}\left(u_A^z+u_C^z-2u_O^z\right)$ (11)

${\theta }_{OA{C}^{\prime }}=\frac{2\sqrt{2}}{3a}\left(u_O^x-u_{C^{\prime }}^x+u_O^y-u_{C^{\prime }}^y\right)-\frac{4\sqrt{2}}{3a}\left(u_O^z+u_{C^{\prime }}^z-2u_A^z\right)$ (12)

${\theta }_{AOD}=\frac{2\sqrt{2}}{3a}\left(u_D^z-u_A^z+u_D^x-u_A^x\right)+\frac{4\sqrt{2}}{3a}\left(u_A^y+u_D^y-2u_O^y\right)$ (13)

${\theta }_{OA{D}^{\prime }}=\frac{2\sqrt{2}}{3a}\left(u_O^z-u_{D^{\prime }}^z+u_O^x-u_{D^{\prime }}^x\right)-\frac{4\sqrt{2}}{3a}\left(u_O^y+u_{D^{\prime }}^y-2u_A^y\right)$ (14)

由(8)式可得力常数矩阵元 ${\partial }^{2}U\left(AO\right)/\partial u_A^x\partial u_O^x$ 和 ${\partial }^{2}U\left(AO\right)/\partial u_A^x\partial u_O^y$ ，分别记为 $\alpha$ 和 $\beta$ 。

$\alpha =-\frac{1}{3}{f}_r-\frac{8}{3}{f}_{\theta }$ (15)

$\beta =-\frac{1}{3}{f}_r+\frac{4}{3}{f}_{\theta }$ (16)

根据非线性力常数 $g_r$ 和 $g_{\theta }$ 的定义以及晶格常数与键长之间的几何关系，得

$g_r=\frac{4}{\sqrt{3}}\frac{\partial f_r}{\partial a}$ (17)

$g_{\theta }=\frac{2}{\sqrt{3}}\frac{\partial f_{\theta }}{\partial a}$ (18)

4. 硅晶体原子间力常数与热膨胀系数计算结果及分析

Rignanese等人根据第一性原理计算得到的最近邻原子 $A$ 和 $O$ 之间的力常数矩阵元 $\alpha$ 和 $\beta$ ，如表1所示，其中 $\text{B}$ 表示长度单位波尔半径， $\text{H}$ 表示能量单位哈特利。Rignanese等人也计算了第二近邻以远原子之间的力常数矩阵元，由于这些力常数矩阵元较小，因此本文不作考虑。利用上面的两体和三体线性力常数及非线性力常数的计算公式，和表1中的 $\alpha$ 和 $\beta$ 数据，得线性及非线性力常数如表1所示，它们采用国际单位。

根据表1中晶格常数为10.26波尔半径时的线性力常数计算晶格动力学矩阵，并求解其本征值问题，再在将表1中的非线性力常数代入(5)式进行计算，得硅晶体二体、三体相互作用对热膨胀系数的贡献以及硅晶体热膨胀系数随温度的变化关系如图1所示。由图1可知，硅晶体热膨胀系数随温度的变化关系的理论计算曲线与Ibach [3] 的实验数据很好地吻合，这表明各线性力常数和非线性力常数的计算结果是正确的。由表1中的数据可知，两体非线性力常数 $g_r$ 为负，由图1中据此计算得到的两体势引起的热膨胀系数随温度变化的关系曲线可知，由两体势引起的热膨胀系数总为正值，这表明为负值的两体非线性力常数 $g_r$ 引起正的热膨胀。由表1中的数据还可知，三体非线性力常数 $g_{\theta }$ 为正，由图1中据此计算得到的三体势引起的热膨胀系数随温度变化的关系曲线可知，由三体势引起的热膨胀系数总是为负值，这表明为正值的三体非线性力常数 $g_{\theta }$ 引起负的热膨胀。将两体势和三体势引起的热膨胀系数相加，即得总的硅晶体热膨胀系数，其随温度变化的关系曲线如图1所示。由该曲线可知，在低温时总的硅晶体呈现负

Table 1. The calculated results of the linear force constants and the non-linear force constants

表1. 线性力常数 $f_r$ 和 $f_{\theta }$ 与非线性力常数 $g_r$ 和 $g_{\theta }$ 的计算结果

Figure 1. The thermal expansion coefficient of silicon crystal vs temperature

图1. 硅晶体热膨胀系数与温度关系

热膨胀性质，这就表明在低温时三体势引起的负热膨胀强于两体非和谐势引起的正热膨胀，从而验证了我们曾推测的硅晶体低温负热膨胀的物理机制 [8] 。

5. 结论

本文计算了硅晶体原子间两体和三体的线性和非线性力常数，再利用热膨胀系数公式计算了硅晶体的热膨胀系数，得到了与实验数据一致的结果。考虑到硅晶体的热学性质与原子间的相互作用密切相关，因此本文得到的硅原子间的线性和非线性力常数，对于运用晶格动力学准确模拟MEMS中的硅晶体构件的热学性质具有重要意义。本文还验证了我们原来提出的硅晶体低温负热膨胀的物理机制 [8] ：为负值的两体非线性力常数引起正热膨胀，而为正值的三体非线性力常数引起负的热膨胀，并且低温时三体势引起的负热膨胀强于两体势引起的正热膨胀，从而使硅晶体呈现低温负热膨胀性质，因此彻底解决了硅晶体低温负热膨胀的机制问题。Stillinger-Weber模型目前被广泛应用于硅晶体热学性质的模拟计算。由于当初提出Stillinger-Weber这种经验模型主要是为了模拟硅晶体的融化性质 [9] ，而不是针对硅晶体的热膨胀性质模拟，因而该模型不能解释硅晶体的低温负热膨胀现象，其根本原因在于根据该模型计算得到的两体和三体非线性力常数都为负，都分别引起正的热膨胀，因此总的热膨胀系数总是正值。针对Stillinger-Weber模型存在的较大缺陷，我们可以根据本文得到的力常数进行修正，使之也能适用于模拟热膨胀和热传导等需要考虑原子间非线性相互作用的情形。