基于晶格动力学的硅单晶热学性质研究(II)——晶格非和谐势能和格律纳森参数

Study on Thermal Properties of Silicon Single Crystal Based on Lattice Dynamics (II)—Lattice Anharmonic Potential and Gruneisen Parameter

摘要: 为用晶格动力学方法研究硅单晶的热膨胀性质，必须首先计算作为微扰哈密顿的硅单晶晶格的三阶非和谐势能，而三阶非和谐势能的计算又可以建立在已知不同晶格常数所对应的晶格和谐势能的基础上，因此本文运用他人计算得到的最近邻和次近邻原子间的相互作用力常数，推导了的硅单晶的晶格和谐势能公式，并在此基础上根据不同晶格常数下的晶格和谐势能推导了三阶非和谐势能公式，为后续硅单晶的热膨胀系数计算公式的推导做了必要准备。最后还运用本文提出的三阶非和谐势能模型，计算了不同温度下硅单晶的格律纳森参数，发现其与他人的实验结果吻合，从而验证了本文所提出的三阶非和谐势能模型的正确性。

Abstract: Abstract: In order to study the thermal expansion properties of silicon single crystal by lattice dynamics, the third order anharmonic potential of silicon single crystal lattice should be calculated as the perturbation Hamiltonian first, and the calculation of the third order anharmonic potential can be based on the lattice harmonic potential corresponding to different lattice constants. Therefore, the formula for lattice harmonic potential energy of silicon single crystal is derived first by using the interaction force constants between the nearest neighbor atoms and the next-nearest neighbor atoms calculated by others, and then the formula for the third order anharmonic potential is derived based on the lattice harmonic potential with different lattice constants. It is a necessary preparation for the subsequent derivations of formula for thermal expansion coefficient. Finally, with the third-order anharmonic potential energy model proposed in this paper, the Gruneisen parameters of silicon single crystal at different temperature are calculated, and it coincides with the experimental results of others. Thus the correctness of the third-order anharmonic potential energy model proposed in this paper is verified.

1. 引言

在文献 [1] 中，我们给出了硅单晶的晶格动力学矩阵，对该矩阵求解本征值问题，即可求得和谐晶体的晶格振动频率和晶格振动位移。用晶格动力学研究硅单晶的热膨胀性质的思路是，将和谐晶体的晶格振动状态作为未微扰态，将硅单晶晶格的三阶非和谐势能作为微扰哈密顿，利用量子力学的微扰论推导原子间距的热膨胀和热膨胀系数公式，在此基础上进行数值计算并进而分析热膨胀机制 [2] [3] [4]。根据在不同晶格常数下的和谐势能，可以得到三阶非和谐势能，即微扰哈密顿。本文将推导硅单晶的和谐势能公式和三阶非和谐势能公式，为热膨胀公式的推导做准备。

格律纳森定律给出了热膨胀系数 $\gamma$ 与单位体积的比热 $C_v$ 之间的关系 [5]

$\gamma =\frac{\rho }{3B}{C}_v$ (1)

其中， $\rho$ 为格律纳森参数，B为体积模量。如果知道格律纳森参数 $\rho$，体积模量B，即可计算得到热膨胀系数 $\gamma$，所以格律纳森参数对材料的热膨胀性质研究是十分重要的 [6]。因此，最后还将运用本文提出的三阶非和谐势能模型，计算不同温度下硅单晶的格律纳森参数，并将其与他人的实验结果相比较，从而验证本文所提出的三阶非和谐势能模型的正确性，还将基于格律纳森参数 $\rho$ 公式对不同的三阶力常数对硅单晶热膨胀性质的影响进行了分析与研究。

2. 原子间两体和谐势能

最近邻原子间存在两体相互作用力，该作用力是中心力，大小与两原子间的距离变化量成正比，其弹性系数设为k，则原子O与最近邻原子A、B、C和D之间的两体作用和谐势能可表示为：

$U_2=\frac{1}{2}k\left(r_{OA}^2+r_{OB}^2+r_{OC}^2+r_{OD}^2\right)$ (2)

其中， $r_{OA}$ 、 $r_{OB}$ 、 $r_{OC}$ 和 $r_{OD}$ 分别为原子A和O、原子B和O、原子C和O及、原子D和O间的距离变化量。

$r_{OA}=\frac{1}{\sqrt{3}}\left(x_{AO}+y_{AO}+z_{AO}\right)$ (3)

$r_{OB}=\frac{1}{\sqrt{3}}\left(-x_{BO}+y_{BO}+z_{BO}\right)$ (4)

$r_{OC}=\frac{1}{\sqrt{3}}\left(x_{CO}+y_{CO}-z_{CO}\right)$ (5)

$r_{OD}=\frac{1}{\sqrt{3}}\left(x_{DO}-y_{DO}+z_{DO}\right)$ (6)

在(3)式~(6)式中， ${\sigma }_{\xi {\xi }^{\prime }}={\sigma }_{{\xi }^{\prime }}-{\sigma }_{\xi }$ 表是原子 ${\xi }^{\prime }$ 相对于原子 $\xi$ 在 $\sigma$ 方向的位移，即 $u_{\sigma }\left({\xi }^{\prime }\right)-u_{\sigma }\left(\xi \right)$。

3. 原子间三体和谐势能

以原子A、O和B为例，它们组成了一个等腰三角形，原子O处于其顶角，原子A和O、原子B和O互为近邻原子，原子A和原子B互为次近邻原子，现只考虑这三个原子间的三体相互作用，并将原子A、O和B之间三体相互作用势按泰勒级数展开，其二阶和谐势能可表示为，

$U_3\left(AOB\right)=\frac{1}{2}{\displaystyle \underset{\sigma {\sigma }^{\prime }}{\sum }{\displaystyle \underset{\eta {\eta }^{\prime }}{\sum }{u}_{\sigma }\left(\eta \right){\Phi }_{\sigma {\sigma }^{\prime }}\left(\eta {\eta }^{\prime }\right)u_{{\sigma }^{\prime }}\left({\eta }^{\prime }\right)}}$ (7)

其中， $u_{\sigma }\left(\eta \right)$ 和 $u_{{\sigma }^{\prime }}\left({\eta }^{\prime }\right)$ 分别代表原子 $\eta$ 和 ${\eta }^{\prime }$ 的晶格振动在 $\sigma$ 和 ${\sigma }^{\prime }$ 方向上的位移， ${\Phi }_{\sigma {\sigma }^{\prime }}\left(\eta {\eta }^{\prime }\right)$ 为原子 $\eta$ 和 ${\eta }^{\prime }$ 之间的力常数， $\eta$ 和 ${\eta }^{\prime }$ 可以为A、O或B。由于O处于等腰三角形AOB的顶角，根据对称性及原子A、O或B的位置关系，可设

${\Phi }_{A0}=\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ {\delta }_1& 0& 0\\ {\delta }_1& 0& 0\end{array}\right)$ (8)

设原子A、O和B之间没有相对运动，则原子A所受合力为0，则有 [1]

${\Phi }_{AA}+{\Phi }_{A0}+{\Phi }_{AB}=0$ (9)

根据(8)式、(9)式和原子A、B之间的力常数矩阵 [1]

${\Phi }_{AB}=\left(\begin{array}{ccc}\lambda & \bar{\delta }& \bar{\delta }\\ \delta & \mu & \nu \\ \delta & \nu & \mu \end{array}\right)$ (10)

可得

${\Phi }_{AA}=-\left(\begin{array}{ccc}{\lambda }_1+\lambda & \bar{\delta }& \bar{\delta }\\ {\delta }_1+\delta & \mu & \nu \\ {\delta }_1+\delta & \nu & \mu \end{array}\right)$ (11)

由于 ${\Phi }_{AA}$ 是对称矩阵，因此上式中 ${\delta }_1=-2\delta$，所以有

${\Phi }_{AA}=-\left(\begin{array}{ccc}{\lambda }_1+\lambda & \bar{\delta }& \bar{\delta }\\ \bar{\delta }& \mu & \nu \\ \bar{\delta }& \nu & \mu \end{array}\right)$ (12)

${\Phi }_{A0}=\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 2\bar{\delta }& 0& 0\\ 2\bar{\delta }& 0& 0\end{array}\right)$ (13)

根据对称性，由(12)和(13)式可得以下两式

${\Phi }_{BB}=-\left(\begin{array}{ccc}{\lambda }_1+\lambda & \delta & \delta \\ \delta & \mu & \nu \\ \delta & \nu & \mu \end{array}\right)$ (14)

${\Phi }_{B0}=\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 2\delta & 0& 0\\ 2\delta & 0& 0\end{array}\right)$ (15)

同理根据 [1]

${\Phi }_{BO}+{\Phi }_{A0}+{\Phi }_{OO}=0$ (16)

可得

${\Phi }_{OO}=\left(\begin{array}{ccc}-2{\lambda }_1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$ (17)

将以上的力常数矩阵代入(7)式，得原子A、B和C之间的三体和谐势能

$U_3\left(AOB\right)=-\frac{1}{2}\left[{\lambda }_1\left(x_{OA}^2+x_{OB}^2\right)+\lambda x_{AB}^2\right]+\delta \left(x_{AO}+x_{BO}\right)\left(y_{AB}+z_{AB}\right)-\frac{1}{2}\mu \left(y_{AB}^2+z_{AB}^2\right)-\nu y_{AB}{z}_{AB}$ (18)

由于文献中 [7] 的数据可知， $\mu -\nu$ 非常接近于0，且可以从对称性加以证明

$y_{AO}{x}_{BO}-x_{AO}{y}_{BO}+z_{AO}{x}_{BO}-x_{AO}{z}_{BO}=0$ (19)

$U_3\left(AOB\right)=-\frac{1}{2}\left\{{\lambda }_1\left(x_{OA}^2+x_{OB}^2\right)+2\delta \left[x_{BO}\left(y_{BO}+z_{BO}\right)-x_{AO}\left(y_{AO}+z_{AO}\right)\right]+\left[\lambda x_{AB}^2+\nu {\left(y_{AB}+z_{AB}\right)}^2\right]\right\}$ (20)

其中，共有三项，分别与A原子和O原子、B原子和O原子，以及A原子和B原子之间的相对运动有关。根据硅晶体的结构，与原子O最近邻的原子还有原子B、C、D，根据对称性，由(18)式并根据对称性，得以原子O为顶点的等腰三角形AOC、AOD的势能

$U_3\left(AOC\right)=-\frac{1}{2}\left\{{\lambda }_1\left(z_{OA}^2+z_{OC}^2\right)+2\delta \left[z_{CO}\left(x_{CO}+y_{CO}\right)-z_{AO}\left(x_{AO}+y_{AO}\right)\right]+\left[\lambda z_{AC}^2+\nu {\left(x_{AC}+y_{AC}\right)}^2\right]\right\}$ (21)

$U_3\left(AOD\right)=-\frac{1}{2}\left\{{\lambda }_1\left(y_{OA}^2+y_{OD}^2\right)+2\delta \left[y_{DO}\left(x_{DO}+z_{DO}\right)-y_{AO}\left(x_{AO}+z_{AO}\right)\right]+\left[\lambda y_{AD}^2+\nu {\left(x_{AD}+z_{AD}\right)}^2\right]\right\}$ (22)

根据硅晶体的结构，与原子A最近邻的原子除O外还有 $B^{\prime }$ 、 $C^{\prime }$ 、 $D^{\prime }$，根据以上方法，同样可得以原子A为顶点的等腰三角形 $OA{B}^{\prime }$ 、 $OA{C}^{\prime }$ 和 $OA{D}^{\prime }$ 的势能

$U_3\left(OA{B}^{\prime }\right)=-\frac{1}{2}\left\{{\lambda }_1\left(x_{OA}^2+x_{A{B}^{\prime }}^2\right)+2\delta \left[x_{A{B}^{\prime }}\left(y_{A{B}^{\prime }}+z_{A{B}^{\prime }}\right)-x_{AO}\left(y_{AO}+z_{AO}\right)\right]+\left[\lambda x_{O{B}^{\prime }}^2+\nu {\left(y_{O{B}^{\prime }}+z_{O{B}^{\prime }}\right)}^2\right]\right\}$ (23)

$U_3\left(OA{C}^{\prime }\right)=-\frac{1}{2}\left\{{\lambda }_1\left(z_{OA}^2+z_{A{C}^{\prime }}^2\right)+2\delta \left[z_{A{C}^{\prime }}\left(x_{A{C}^{\prime }}+y_{A{C}^{\prime }}\right)-z_{AO}\left(x_{AO}+y_{AO}\right)\right]+\left[\lambda z_{O{C}^{\prime }}^2+\nu {\left(x_{O{C}^{\prime }}+y_{O{C}^{\prime }}\right)}^2\right]\right\}$ (24)

$U_3\left(OA{D}^{\prime }\right)=-\frac{1}{2}\left\{{\lambda }_1\left(y_{OA}^2+y_{OD}^2\right)+2\delta \left[y_{A{D}^{\prime }}\left(x_{A{D}^{\prime }}+z_{A{D}^{\prime }}\right)-y_{AO}\left(x_{AO}+z_{AO}\right)\right]+\left[\lambda y_{O{D}^{\prime }}^2+\nu {\left(x_{O{D}^{\prime }}+z_{O{D}^{\prime }}\right)}^2\right]\right\}$ (25)

以上(20)式~(25)式是共有OA键的所有三体和谐势公式。

同样可以分别计算共有BO键、共有CO键和共有DO键的六个三体势，显然，计算重复了一倍，因此在整个晶体内对它们求和之后应该将求和结果除以2。

根据对称性，把三体和谐势能的所有结果相加，得：

$U_3=\underset{{\xi }_1}{{\displaystyle \sum }}{U}_3\left({\xi }_{1}O\right)+\underset{{\xi }_1\ne {{\xi }^{\prime }}_1}{{\displaystyle \sum }}{U}_3\left({\xi }_1{{\xi }^{\prime }}_1\right)+\underset{{\xi }_2\ne {{\xi }^{\prime }}_2}{{\displaystyle \sum }}{U}_3\left({\xi }_2{{\xi }^{\prime }}_2\right)$ (26)

其中， ${\xi }_1$ 和 ${{\xi }^{\prime }}_1$ 为与原子O与最近邻的所有1类原子A、B、C和D， ${\xi }_2$ 和 ${{\xi }^{\prime }}_2$ 为与原子A与最近邻的所有2类原子O、 $B^{\prime }$ 、 $C^{\prime }$ 和 $D^{\prime }$。 $U_3\left({\xi }_{1}O\right)$ 为三体势中与原子 ${\xi }_1$ 和O间相对运动有关的相互作用势能，可表示为：

$U_3\left(\xi O\right)=-6\delta r_{O\xi }^2-\frac{1}{2}\left\{\left({\lambda }_1-4\delta \right)d_{O\xi }^2\right\}$ (27)

其中， $r_{O\xi }$ 可根据(3)~(6)式计算，为原子O与所有最近邻原子A、B、C和D间的键长变化量。 $d_{O\xi }^2$ 为原子O与所有最近邻原子A、B、C和D间位移之差的模的平方。

$d_{O\xi }^2=x_{O\xi }^2+y_{O\xi }^2+z_{O\xi }^2$ (28)

(26)式中， $U_3\left({\xi }_1{{\xi }^{\prime }}_1\right)$ 为三体势中与原子 ${\xi }_1$ 和 ${{\xi }^{\prime }}_1$ 之间相对运动有关的势能， $U_3\left({\xi }_2{{\xi }^{\prime }}_2\right)$ 为三体势中与原子 ${\xi }_2$ 、 ${{\xi }^{\prime }}_2$ 之间相对运动有关的势能，不重复计算。

$U_3\left(AB\right)$ 、 $U_3\left(CD\right)$ 、 $U_3\left(O{B}^{\prime }\right)$ 和 $U_3\left(C^{\prime }{D}^{\prime }\right)$ 可表示为

$U_3\left({\xi }_i{{\xi }^{\prime }}_i\right)=-\frac{1}{2}\left\{\lambda x_{{\xi }_i{{\xi }^{\prime }}_i}^2+\nu {\left(y_{{\xi }_i{{\xi }^{\prime }}_i}+z_{{\xi }_i{{\xi }^{\prime }}_i}\right)}^2\right\}$ (29)

$U_3\left(AC\right)$ 、 $U_3\left(BD\right)$ 、 $U_3\left(O{C}^{\prime }\right)$ 和 $U_3\left(B^{\prime }{D}^{\prime }\right)$ 可表示为

$U_3\left({\xi }_i{{\xi }^{\prime }}_i\right)=-\frac{1}{2}\left\{\lambda z_{{\xi }_i{{\xi }^{\prime }}_i}^2+\nu {\left(x_{{\xi }_i{{\xi }^{\prime }}_i}+y_{{\xi }_i{{\xi }^{\prime }}_i}\right)}^2\right\}$ (30)

$U_3\left(AD\right)$ 、 $U_3\left(BC\right)$ 、 $U_3\left(O{D}^{\prime }\right)$ 和 $U_3\left(B^{\prime }{C}^{\prime }\right)$ 可表示为

$U_3\left({\xi }_i{{\xi }^{\prime }}_i\right)=-\frac{1}{2}\left\{\lambda y_{{\xi }_i{{\xi }^{\prime }}_i}^2+\nu {\left(x_{{\xi }_i{{\xi }^{\prime }}_i}+z_{{\xi }_i{{\xi }^{\prime }}_i}\right)}^2\right\}$ (31)

以上各式中， $i=1,2$。

根据(2)式和(27)式，可以得到与原子O、 $\xi$ 间相对运动有关的总相互作用势能

$U\left(\xi O\right)=\left(\frac{1}{2}k-6\delta \right)r_{O\xi }^2-\frac{1}{2}\left({\lambda }_1-4\delta \right)d_{O\xi }^2$ (32)

根据力常数定义式 [1]，可知 ${\partial }^{2}U\left(AO\right)/\partial x_A\partial x_O=\alpha$， ${\partial }^{2}U\left(AO\right)/\partial x_A\partial y_O=\beta$，最后可将(32)式表示为

$U\left(\xi O\right)=-\frac{3}{2}\beta r_{O\xi }^2-\frac{1}{2}\left(\alpha -\beta \right)\left\{d_{O\xi }^2\right\}$ (33)

由(29)~(31)式以及(33)式可知，和谐晶体晶格振动势能与力常数存在线性关系，因此相关各项中写入相关力常数 $\alpha ,\beta ,\lambda ,\nu$。综合各式可得

$U\left(\alpha ,\beta ,\lambda ,\nu \right)=\underset{\xi }{{\displaystyle \sum }}U\left(\xi O,\alpha ,\beta \right)+\underset{{\xi }_1\ne {{\xi }^{\prime }}_1}{{\displaystyle \sum }}{U}_3\left({\xi }_1{{\xi }^{\prime }}_1,\lambda ,\nu \right)+\underset{{\xi }_2\ne {{\xi }^{\prime }}_2}{{\displaystyle \sum }}{U}_3\left({\xi }_2{{\xi }^{\prime }}_2,\lambda ,\nu \right)$ (34)

根据Rignanese等人 [7] 从第一性原理出发，计算得到的二阶力常数 $\alpha$ 、 $\beta$ 、 $\lambda$ 和 $\nu$，及在参考文献 [1] 中我们给出的晶格振动位移，利用(34)式即可计算晶格和谐势能。

4. 原子间三阶非和谐势能

在研究晶格热学性质时，一般采用泰勒级数对原子间相互作用势能进行展开，在一般无需考虑非和谐势能的情况下，只需展开至二次项，在计算热传导和热膨胀时，才需要展开至三次以上的高次项。由于四次非和谐势能相对于平衡位置是对称的，因而不会对热膨胀有贡献，因此作为很好的近似，在计算热膨胀时只需展开至三次非和谐势能项。以下根据和谐势能项公式，和力常数对晶格常数的导数，计算三阶非和谐势能。

根据Rignanese等人在参考文献 [7] 中给出的晶格常数下的二阶力常数 $\alpha$ 、 $\beta$ 、 $\lambda$ 和 $\nu$，可以计算力常数 $\alpha$ 、 $\beta$ 、 $\lambda$ 和 $\nu$ 对晶格常数的一阶导数 ${\alpha }^{\prime }$ 、 ${\beta }^{\prime }$ 、 ${\lambda }^{\prime }$ 和 ${\nu }^{\prime }$。由上式可得晶格振动势能对晶格常数的导数

$\frac{\text{d}U}{\text{d}a}=\underset{\xi }{{\displaystyle \sum }}U\left(\xi O,{\alpha }^{\prime },{\beta }^{\prime }\right)+\underset{{\xi }_1\ne {{\xi }^{\prime }}_1}{{\displaystyle \sum }}{U}_3\left({\xi }_1{{\xi }^{\prime }}_1,{\lambda }^{\prime },{\nu }^{\prime }\right)+\underset{{\xi }_2\ne {{\xi }^{\prime }}_2}{{\displaystyle \sum }}{U}_3\left({\xi }_2{{\xi }^{\prime }}_2,{\lambda }^{\prime },{\nu }^{\prime }\right)$ (35)

可见，只需将和谐势能公式中的力常数用其对晶格常数的导数替代即可得晶格振动势能对晶格常数的导数。

由各方向晶格常数的变化 $\Delta a_{\sigma }$ 引起的晶格体积的变化为

$\Delta V=a^2\left(\Delta a_x+\Delta a_y+\Delta a_z\right)$ (36)

再根据 $\Delta U=\Delta V\partial U/\partial V$ 可得

$\Delta U=\frac{1}{3}U\left({\alpha }^{\prime },{\beta }^{\prime },{\lambda }^{\prime },{\nu }^{\prime }\right)\left(\Delta a_x+\Delta a_y+\Delta a_z\right)$ (37)

上式中，

$\Delta a_x=x_{AO}+x_{BO}+x_{OC}+x_{OD}$ (38)

$\Delta a_y=y_{AO}+y_{OB}+y_{OC}+y_{DO}$ (39)

$\Delta a_z=z_{AO}+z_{OB}+z_{CO}+z_{OD}$ (40)

将(38)式~(40)式代入(37)式，得

$\Delta U=\frac{1}{\sqrt{3}}U\left({\alpha }^{\prime },{\beta }^{\prime },{\lambda }^{\prime },{\nu }^{\prime }\right)\left(r_{AO}+r_{BO}+r_{CO}+r_{DO}\right)$ (41)

上式中， $r_{\xi O}$ 为原子间距离的变化量。根据我们在参考文献 [1] 中给出的晶格振动位移公式、和根据Rignanese给出的不同晶格常数下力常数而得到的力常数对晶格常数的导数，即可根据上式计算三阶非和谐势能。

5. 格律纳森参数

在准和谐近似(QHA)下，格律纳森参数 $\rho$ 可以通过下式计算 [5]

$\rho =\frac{1}{C_v\left(T\right)}\underset{kj}{{\displaystyle \sum }}{C}_{vkj}\left(T\right){\rho }_{kj}$ (42)

其中， $C_{vkj}\left(T\right)$ 为 $kj$ 模态声子的每单位体积的比热，而 $C_v\left(T\right)$ 为其总和， ${\rho }_{kj}$ 为模态格律纳森参数， $C_{vkj}\left(T\right)$ 和 ${\rho }_{kj}$ 分别由下式定义

$C_{vkj}\left(T\right)=\frac{\hslash {\omega }_{kj}}{V}\frac{\text{d}{\bar{n}}_{kj}}{\text{d}T}$ (43)

${\rho }_{kj}=-\frac{\text{d}\ln {\omega }_{kj}}{\text{d}\ln V}$ (44)

由上式和(42)式可知，格律纳森参数ρ无量纲。将(43)式和(44)式代入(42)式，并将体积用晶格常数表示，得：

$\rho =-\frac{a}{3V{C}_v\left(T\right)}\frac{{\text{d}}^2}{\text{d}T\text{d}a}\underset{kj}{{\displaystyle \sum }}\hslash {\omega }_{kj}\left({\bar{n}}_{kj}+\frac{1}{2}\right)$ (45)

上式中对求和号内各项求和即得所有声子能量总和，即为晶格振动能量E，其包含动能 $E_k$ 和 $E_p$ 势能。与此相对应，系统哈密顿算子H也有动能哈密顿算子 $H_k$ 和势能哈密顿算子 $H_p$ (即势能U)这两部分构成，能量为对应哈密顿算子的热力学平均值，即 $E_k={\langle H_k\rangle }_T$， $E_p={\langle H_p\rangle }_T$， $E={\langle H\rangle }_T$。因此，上式变为

$\rho =-\frac{a}{3V{C}_v\left(T\right)}\frac{{\text{d}}^2}{\text{d}T\text{d}a}\left({\langle H_k\rangle }_T+{\langle H_p\rangle }_T\right)$ (46)

由于热膨胀会引起晶格常数的变化，只会导致其势能哈密顿的变化，因而由上式得，

$\rho =-\frac{a}{3V{C}_v\left(T\right)}\frac{\text{d}}{\text{d}T}{\langle \frac{\text{d}U}{\text{d}a}\rangle }_T$ (47)

将(34)式代入上式，并考虑到由于 ${\lambda }^{\prime }$ 和 ${\nu }^{\prime }$ 比 ${\alpha }^{\prime }$ 和 ${\beta }^{\prime }$ 小很多，可以忽略(34)式中第二和第三项，并考虑到由于晶体结构对称性第一项求和式中所有四项所产生的热膨胀都等于 $\langle U\left(AO,{\alpha }^{\prime },{\beta }^{\prime }\right)\rangle$，得

$\rho =-\frac{4a}{3VC\left(T\right)}\frac{\text{d}}{\text{d}T}U\left(AO,{\alpha }^{\prime },{\beta }^{\prime }\right)$ (48)

将(2)式、(28)式和(33)代入上式，得

$\rho =\frac{2a}{3C\left(T\right)}\frac{\text{d}}{\text{d}T}\langle {\beta }^{\prime }{\left(x_{AO}+y_{AO}+z_{AO}\right)}^2+\left({\alpha }^{\prime }-{\beta }^{\prime }\right){\left(x_{OA}^2+y_{OA}^2+z_{OA}^2\right)}_T\rangle$ (49)

将文献(1)中的公式(26)代入，得

$\rho =\frac{2a\hslash }{3mC\left(T\right)}\underset{kj}{{\displaystyle \sum }}\left\{\frac{{\beta }^{\prime }}{{\omega }_{kj}}{\left[\underset{\sigma }{{\displaystyle \sum }}{e}_{\alpha }\left(O|kj\right)-e_{\alpha }\left(A|kj\right)\right]}^2+\frac{{\alpha }^{\prime }-{\beta }^{\prime }}{{\omega }_{kj}}\underset{\sigma }{{\displaystyle \sum }}{\left[e_{\alpha }\left(O|kj\right)-e_{\alpha }\left(A|kj\right)\right]}^2\right\}\frac{\partial {\bar{n}}_{kj}}{\partial T}$ (50)

根据Rignanese等人在参考文献 [7] 中提供的不同晶格常数下的二阶力常数 $\alpha$ 和 $\beta$ 而计算出的它们对晶格常数的一阶导数 ${\alpha }^{\prime }=0.0200h/b^3$ 和 ${\beta }^{\prime }=0.02638h/b^3$，并代入上式可计算不同温度下的格律纳森参数 $\rho$，结果如图1所示，图中还列出了为了Ibach的实验结果 [8]，可见我们的计算结果与实验数据能较好地吻合，说明本文给出的原子间三阶非和谐势能模型及格律纳森参数 $\rho$ 公式是正确的。

由于(1)式中的体积模量B和单位体积的比热 $C_v$ 都为正值，因此，热膨胀系数的正负取决于格律纳森参数 $\rho$ 的正负，图1中硅单晶在约120 K温度下格律纳森参数为负，表明在该温度段有负热膨胀系数，而在其它温度段的格律纳森参数 $\rho$ 为正，表明有正的热膨胀系数。根据 ${\alpha }^{\prime }=0.0200h/b^3$ 和 ${\beta }^{\prime }=0.02638h/b^3$ 得知 ${\alpha }^{\prime }$ 为正、 ${\alpha }^{\prime }-{\beta }^{\prime }$ 为负，由(50)式可知，这就有可能使得某些模态的声子对格律纳森参数 $\rho$ 形成正的贡献，而又有些模态的声子对格律纳森参数 $\rho$ 形成负的贡献，并且随温度的变化，各个模态对格律纳森参数 $\rho$ 贡献的大小也会有相应的调整，从而使得格律纳森参数 $\rho$ 在不同温度段呈现不同的正负特性，例如硅单晶在约120 K温度下呈现反常的负热膨胀的性质，而在高温下呈现正常的正热膨胀性质。在后续的论文中，我们将在晶格动力学和微扰论的基础上对硅单晶的热膨胀系数公式进行严密的推导，并在数值计算的基础上更为详细探讨硅单晶低温负热膨胀性质的物理机制。

Figure 1. The gruneisen constant of silicon single crystal vs temperature

图1. 硅单晶的格律纳森常数与温度的关系

6. 结论

本文运用他人计算得到的最近邻和次近邻原子间的相互作用力常数，推导了的硅单晶的晶格和谐势能公式，并在此基础上根据不同晶格常数下的晶格和谐势能推导了三阶非和谐势能公式。本文还根据这个三阶非和谐势能模型，推导了在准和谐近似(QHA)下计算硅单晶的格律纳森参数公式，并将计算结果与他人的实验结果进行了比较，从而证明了本文提出的三阶非和谐势能模型是正确的，这为后续硅单晶的热膨胀系数的计算做了必要准备。本文最后对硅单晶的负热膨胀的物理机制进行了初步分析，认为因为 ${\alpha }^{\prime }$ 为正、 ${\alpha }^{\prime }-{\beta }^{\prime }$ 为负，这使得某些模态的声子对格律纳森参数 $\rho$ 形成正的贡献，而又有些模态的声子对格律纳森参数 $\rho$ 形成负的贡献，并且随温度的变化，各个模态对格律纳森参数 $\rho$ 贡献的大小也会有相应的调整，从而使得格律纳森参数 $\rho$ 在不同温度段呈现不同的正负特性，从而使硅单晶在约120 K温度下呈现反常的负热膨胀的性质，而在高温下呈现正常的正热膨胀性质。

NOTES

^*通讯作者。