单原子层薄膜热传导性质的晶格动力学研究(I)——声子线宽和热传导系数公式

Lattice Dynamics Study on the Thermal Conduction Properties of Single Atomic Layer Films (I)—Formulas for Phonon Linewidth and Thermal Conductivity

摘要: 本文在运用晶格动力学和Hardy能量通量公式得到单原子层薄膜的晶格振动的频率、原子位移、原子动量、晶格振动能量，非和谐势能和能量通量等公式，在此基础上运用Green函数理论和Green-Kubo公式推导单原子层薄膜的声子谱线宽度公式和热传导系数公式，结果表明单原子层薄膜的热传导系数为所有声子的热传导系数之和，而单个声子的热传导系数与其速度、声子能量、声子寿命或自由程密切相关。

Abstract: The formulas for lattice vibration frequency, atomic displacement and momentum, lattice vibration energy, an harmonic potential energy and energy flux of lattice vibration in single atomic layer film are derived in this paper on the basis of the lattice dynamics theory, and then the formulas for phonon line width and thermal conductivity are derived with the aid of Green function theory and Green-Kubo formula. The result shows that the thermal conductivity of the film is the sum of contribution from every single phonon which is closely related to phonon’s velocity, energy and lifetime or free path.

1. 引言

以石墨烯为代表的单原子层薄膜，由于独特的结构和各种优异的性质，在诸多领域中具有广阔的应用前景，已成为研究热点 [1] [2] [3] [4]。随着制造工艺与技术的进步，微米纳米电子器件的集成度进一步提高，从而导致其热效应愈加严重，散热问题更加突出。单原子层薄膜是新型纳米电子器件中的重要材料，由于它们的性质与大块材料截然不同，因此我们必须对它们的热传导性质进行深入研究，这对于纳米电子器件的散热设计有着非常重要的意义。

国内外已经有许多关于单原子层薄膜热学性质理论研究的报导，其中分子动力学模拟是一种常用方法。分子动力学模拟 [3] [4] [5] 建立在牛顿力学的基础上，通过对多原子的运动方程进行求解来计算热性质。例如李保文等人 [5] 用分子动力学模拟方法研究了石墨烯的热传导率，发现石墨烯的热导率随着其尺寸长度的增加而增加。但分子动力学存在一些局限性，如，由于纳米电子器件的材料尺寸很小，再加上低温等因素，肯定会存在量子效应，而由于分子动力学一般不考虑量子效应，必然会带来较大的计算误差；分子动力学模拟的计算量非常大，计算时间长，对计算机要求高；分子动力学只能得到数值计算结果，难以根本上从数值计算结果中分析物理机制。

本文将推导单原子层薄膜的热传导系数公式。首先求解晶格动力学问题，得到晶格振动频率、原子位移、原子动量、晶格振动能量、群速度和非和谐势能公式；然后将非和谐势能作为微扰，利用Green函数理论，推导声子线宽公式；最后根据Hardy能量通量公式推导单原子层薄膜的能量通量公式，再结合Green-Kubo公式得到其热传导系数公式，为对单原子层薄膜的热传导性质进行数值计算做准备。

2. 单原子层薄膜的晶格动力学理论

设单原子层薄膜具有正方结构，晶格常数为a。在x和y方向都有N个质量为m的原子，并采用玻恩–卡曼循环边界条件。为了方便计算，原子间的相互作用采用M-P模型 [6]，即只考虑最近邻原子间的相互作用，并且原子间径向力常数与切向力常数相同，表示为k。设原子在薄膜中所处的位置为 $la$，其中， $l=\left(l_x,l_y\right)$。在和谐近似下，求解晶格动力学问题，可得以下声子频率、晶格振动能量、原子位移和动量公式。

${\omega }_k^2=\frac{2k}{m}\left(2-\cos k_{x}a-\cos k_{y}a\right)$ (1)

$H=\underset{k\sigma }{{\displaystyle \sum }}\left(a_{k\sigma }^{+}{a}_{k\sigma }+\frac{1}{2}\right)\hslash {\omega }_k$ (2)

$u_{\sigma }\left(l\right)=\frac{1}{N}\underset{k\sigma }{{\displaystyle \sum }}\sqrt{\frac{\hslash }{2m{\omega }_k}}{A}_{k\sigma }{\text{e}}^{ia\stackrel{\to }{k}\cdot \stackrel{\to }{l}}$ (3)

$p_{\sigma }\left(l\right)=-\frac{i}{N}\underset{k\sigma }{{\displaystyle \sum }}\sqrt{\frac{\hslash m{\omega }_k}{2}}{B}_{k\sigma }{\text{e}}^{ia\stackrel{\to }{k}\cdot \stackrel{\to }{l}}$ (4)

其中， $k$ 为波矢， $\sigma$ 为晶格振动的偏振方向即声子的偏振方向， ${\omega }_k$ 为晶格振动频率，其中 $a_{k\sigma }^{+}$ 、 $a_{k\sigma }$ 表示声子的产生与湮灭算符， $a_{k\sigma }^{+}{a}_{k\sigma }$ 是声子数算符，用 $n_{k\sigma }$ 表示， $A_{k\sigma }=a_{k\sigma }+a_{-k\sigma }^{+}$， $B_{k\sigma }=a_{k\sigma }-a_{-k\sigma }^{+}$。根据晶格动力学计算可知，在M-P模型中，单原子薄膜有三个声学支，对应晶格振动的三个不同的偏振方向， $\sigma$ 可为x、y和z，但不论声子偏振方向为何，晶格振动频率即声子频率只取决于波矢 $k$ 而与声子偏振方向无关，因此根据玻色统计可知，声子布居数 ${\bar{n}}_{k\sigma }$ 也只与波矢有关而与偏振方向无关，可记为 ${\bar{n}}_k$。

热传导性质与原子间的非和谐势能有关，通常只需计算三阶非和谐势能，在M-P模型下，三阶非和谐势能为

${H^{\prime }}_{\sigma }=\frac{\delta }{6}\underset{l_x,l_y}{{\displaystyle \sum }}\left\{{\left[u_{\sigma }\left(l_x,l_y\right)-u_{\sigma }\left(l_x-1,l_y\right)\right]}^3+{\left[u_{\sigma }\left(l_x,l_y+1\right)-u_{\sigma }\left(l_x,l_y\right)\right]}^3\right\}$ (5)

其中， $\delta$ 为三阶力常数。将晶格原子位移公式代入(5)式中，可得

${H^{\prime }}_{\sigma }=\underset{k{k}^{\prime }{k}^{″}}{{\displaystyle \sum }}V\left(k,k^{\prime },k^{″}\right)A_{\sigma }\left(k\right)A_{\sigma }\left(k^{\prime }\right)A_{\sigma }\left(k^{″}\right)$ (6)

其中，当对于所有的 $\sigma =x,y,z$，都有 ${\left(k+k^{\prime }+k^{″}\right)}_{\sigma }a$ 为2π的整数倍时， $\Delta \left(k+k^{\prime }+k^{″}\right)$ 为1，其它时候为0。

$V\left(k,k^{\prime },k^{″}\right)=\frac{2{\delta }^2{\hslash }^3}{9N^2{m}^3}\frac{\Delta \left(k+k^{\prime }+k^{″}\right)}{\omega \left(k\right)\omega \left(k^{\prime }\right)\omega \left(k^{″}\right)}{\left|\underset{\eta =x,y}{{\displaystyle \sum }}\sin \frac{k_{\eta }a}{2}\sin \frac{{k^{\prime }}_{\eta }a}{2}\sin \frac{{k^{″}}_{\eta }a}{2}{\text{e}}^{-i\left(k_{\eta }+{k^{\prime }}_{\eta }+{k^{″}}_{\eta }\right)a/2}\right|}^2$ (7)

3. 声子Green函数

在考虑三阶非和谐势能的情况下声子Green函数 [7]

$G_{k\sigma ,k^{\prime }{\sigma }^{\prime }}^{AA}\left(\omega \right)=\frac{{\omega }_k{\delta }_{k,-k^{\prime }}{\delta }_{\sigma ,{\sigma }^{\prime }}}{\pi \left[{\omega }^2-{\omega }_k^2-2{\omega }_k{M}_{k\sigma }\left(\omega \right)\right]}$ (8)

其中

$M_{k\sigma }\left(\omega \right)=\frac{18\pi }{{\hslash }^2}\underset{k_1{k}_2{q}_1{q}_2}{{\displaystyle \sum }}{V}_3\left(k_1,k_2,k\right)V_3\left(q_1,q_2,-k\right)\langle A_{k_1\sigma }\left(t\right)A_{k_2\sigma }\left(t\right),A_{q_1\sigma }\left(t^{\prime }\right)A_{q_2\sigma }{\left(t^{\prime }\right)}_{\omega }\rangle$ (9)

根据解耦公式 [8]，并考虑到声子存在线宽的情况下有以下相关函数

$\langle a_{k\sigma }^{+}\left(t\right)a_{k^{\prime }{\sigma }^{\prime }}\left(0\right)\rangle ={\delta }_{k{k}^{\prime }}{\delta }_{\sigma {\sigma }^{\prime }}{n}_k{\text{e}}^{i{\omega }_{k}t-{\Gamma }_{k\sigma }\left|t\right|}$ (10)

$\langle a_{k\sigma }\left(t\right)a_{k^{\prime }{\sigma }^{\prime }}^{+}\left(0\right)\rangle ={\delta }_{k{k}^{\prime }}{\delta }_{\sigma {\sigma }^{\prime }}\left(n_k+1\right){\text{e}}^{i{\omega }_{k}t-{\Gamma }_{k\sigma }\left|t\right|}$ (11)

可得

$M_{k\sigma }\left(\omega \right)=\frac{36}{{\hslash }^2}\underset{k_1{k}_2}{{\displaystyle \sum }}\underset{\pm }{{\displaystyle \sum }}\left[\left({\bar{n}}_{k_2}+\frac{1}{2}\right)\pm \left({\bar{n}}_{k_1}+\frac{1}{2}\right)\right]\frac{{\left|V\left(k_1;k_2;-k\right)\right|}^2\left({\omega }_{k_1}\pm {\omega }_{k_2}\right)}{{\left[\omega +i\left({\Gamma }_{k_1\sigma }+{\Gamma }_{k_2\sigma }\right)\right]}^2-{\left({\omega }_{k_1}\pm {\omega }_{k_2}\right)}^2}$ (12)

可以将 $M_{k\sigma }\left(\omega \right)$ 分解为实部和虚部之和，即 $M_{k\sigma }\left(\omega \right)={\Delta }_{k\sigma }-i{\Gamma }_{k\sigma }$。声子Green函数极点的实部表示声子频谱，极点的虚部则表示声子谱线宽度 [7]

${\Gamma }_{k\sigma }=\frac{72}{{\hslash }^2}\underset{k_1{k}_2}{{\displaystyle \sum }}\underset{\pm }{{\displaystyle \sum }}\left[\left({\bar{n}}_{k_2}+\frac{1}{2}\right)\pm \left({\bar{n}}_{k_1}+\frac{1}{2}\right)\right]\frac{\omega \left({\omega }_{k_1}\pm {\omega }_{k_2}\right)\left({\Gamma }_{k_1\sigma }+{\Gamma }_{k_2\sigma }\right){\left|V\left(k_1;k_2;-k\right)\right|}^2}{{\left[{\omega }^2-\left({\omega }_{k_1}\pm {\omega }_{k_2}\right)\right]}^2+{\left[2\omega \left({\Gamma }_{k_1\sigma }+{\Gamma }_{k_2\sigma }\right)\right]}^2}$ (13)

上式是声子谱线宽度的迭代计算公式，设置声子谱线宽度的迭代初值后，经过有限次迭代计算，就可以得到声子谱线宽度的准确计算结果。由于在M-P模型中声子频率、布居数及 $\left|V\left(k_1;k_2;-k\right)\right|$ 等与偏振方向无关，所以上式迭代计算的结果 ${\Gamma }_{k\sigma }$ 即声子线宽也与偏振方向无关，而只与其波矢有关，因此可将(10)式、(11)式及上式中的 ${\Gamma }_{k\sigma }$ 、 ${\Gamma }_{k_1\sigma }$ 和 ${\Gamma }_{k_2\sigma }$ 分别记为 ${\Gamma }_k$ 、 ${\Gamma }_{k_1}$ 和 ${\Gamma }_{k_2}$。

若将(13)式中的 ${\Gamma }_{k_1}+{\Gamma }_{k_2}$ 取为0⁺，(13)式就可变为Maradudin的声子线宽公式 [6]，Maradudin在进行数值计算时，用一个常数取代0⁺，如果常数选择太小，声子之间不会进行足够的碰撞，如果选择较大，则可能让一些本来无法碰撞的声子发生碰撞，从而产生误差。Turney [9] 等人建议选择 ${\Gamma }_k+{\Gamma }_{k_1}+{\Gamma }_{k_2}$ 替代Maradudin声子线宽公式中的0⁺，相当于在(13)式中用 ${\Gamma }_k+{\Gamma }_{k_1}+{\Gamma }_{k_2}$ 代替 ${\Gamma }_{k_1}+{\Gamma }_{k_2}$，肯定会带来热传导系数计算结果的误差。(13)式是我们经过严密的推导得到的，很明显比Maradudin和Turney等人的公式更为合理。

4. 单原子层薄膜的能量通量公式

可用Green-Kubo公式 [10] 计算 方向的热传导系数 ${\kappa }_x$

${\kappa }_x=\frac{k_B{\beta }^2}{V}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}t}\langle S_x\left(t^{\prime }\right)S_x\left(t\right)\rangle$ (14)

其中， $S_x\left(t\right)$ 为x方向的能量通量。

在固体物理中广泛运用Hardy能量通量公式 [11] 计算能量通量

$S=\frac{1}{2i\hslash m}\underset{i,j}{{\displaystyle \sum }}\left[R\left(i\right)-R\left(j\right)\right]\underset{\sigma }{{\displaystyle \sum }}\left\{p_{\sigma }\left(i\right)\left[p_{\sigma }\left(i\right),V_j\right]+\left[p_{\sigma }\left(i\right),V_j\right]p_{\sigma }\left(i\right)\right\}$ (15)

其中， $V_j$ 表示与原子 $j=\left(l_x,l_y\right)$ 有关的相互作用势能，原子i的平衡位置和动量分别用 $R\left(i\right)$ 和 $p\left(i\right)$ 表示，动量在 $\sigma$ 方向分量的算符定义为 $p_{\sigma }\left(i\right)=-i\hslash \partial /\partial u_{\sigma }\left(i\right)$。

现根据Hardy公式求解x方向的能量通量 $S_x$。由于 $S_x$ 只与 $\left[R\left(i\right)-R\left(j\right)\right]$ 在x方向分量有关，所以在计算 $S_x$ 时，原子i只能是j水平方向上的相邻原子 $i_1=\left(l_x-1,l_y\right)$， $i_2=\left(l_x+1,l_y\right)$。由(15)式，得

$S_x=\frac{a}{2m}\underset{\sigma }{{\displaystyle \sum }}\left\{\underset{j{i}_1}{{\displaystyle \sum }}\left[p_{\sigma }\left(i_1\right)\frac{\partial V_j}{\partial u_{\sigma }\left(i_1\right)}+\frac{\partial V_j}{\partial u_{\sigma }\left(i_1\right)}{p}_{\sigma }\left(i_1\right)\right]-\underset{j{i}_2}{{\displaystyle \sum }}\left[p_{\sigma }\left(i_2\right)\frac{\partial V_j}{\partial u_{\sigma }\left(i_2\right)}+\frac{\partial V_j}{\partial u_{\sigma }\left(i_2\right)}{p}_{\sigma }\left(i_2\right)\right]\right\}$ (16)

其中

$V_j^{\sigma }=\frac{1}{2}\underset{\sigma }{{\displaystyle \sum }}\frac{1}{2}k\left\{{\left[u_{\sigma }\left(i_1\right)-u_{\sigma }\left(j\right)\right]}^2+{\left[u_{\sigma }\left(i_2\right)-u_{\sigma }\left(j\right)\right]}^2\right\}$ (17)

其中，第一个1/2是因为对j求和计算总的系统势能时，对所有原子间相互作用势能重复计算了一次。将(17)式代入(16)式中，并利用(3)式和声子产生和消灭算符的性质，得

$S_x={\displaystyle \underset{k\sigma }{\sum }{n}_{k\sigma }\hslash {\omega }_k{v}_k^x}$ (18)

其中， $v_k^x=\partial {\omega }_k/\partial k_x$，为声子在x方向上的群速度。将(1)式两边对 $k_x$ 求导，根据群速度的定义，可得

$v_k^x=\frac{ka}{m{\omega }_k}\sin k_{x}a$ (19)

可知，在M-P模型中声子群速只取决于其波矢，而与其偏振无关。

由(18)可知，晶格振动的能量通量等于不同模态的所有声子能量通量之和，某种模态的声子能量通量等于该模态的声子数、单个声子的能量和声子群速度的乘积。

5. 单原子层薄膜的热传导系数公式

将(18)式代入(14)式，得到

${\kappa }_x=\frac{k_B{\hslash }^2{\beta }^2}{V}\underset{k\sigma k^{\prime }{\sigma }^{\prime }}{{\displaystyle \sum }}{\omega }_k{\omega }_{k^{\prime }}{v}_k^x{v}_{k^{\prime }}^x{\displaystyle {\int }_{-\infty }^{\infty }\text{d}t}\langle n_{k\sigma }\left(t\right)n_{k^{\prime }{\sigma }^{\prime }}\left(t^{\prime }\right)\rangle$ (20)

根据解耦公式 [8]，以及声子算符性质，(20)式可以变为

${\displaystyle {\int }_{-\infty }^{\infty }\text{d}t}\langle n_{k\sigma }\left(t\right)n_{k^{\prime }{\sigma }^{\prime }}\left(t^{\prime }\right)\rangle ={\displaystyle {\int }_{-\infty }^{\infty }\text{d}t}\langle a_{k\sigma }^{+}\left(t\right)a_{k^{\prime }{\sigma }^{\prime }}\left(0\right)\rangle \langle a_{k\sigma }\left(t\right)a_{k^{\prime }{\sigma }^{\prime }}^{+}\left(0\right)\rangle$ (21)

再将(10)式和(11)式代入上式，计算可得到热传导系数公式

${\kappa }_x=\frac{k_B{\hslash }^2{\beta }^2}{V}\underset{k\sigma }{{\displaystyle \sum }}\frac{{\omega }_k^2{\left(v_k^x\right)}^2{\bar{n}}_{k\sigma }\left({\bar{n}}_{k\sigma }+1\right)}{{\Gamma }_k}$ (22)

其中， $1/{\Gamma }_k$ 为声子寿命 ${\tau }_k$， $v_k^x/{\Gamma }_k$ 为声子自由程 $L_k^x$， ${\bar{n}}_{k\sigma }$ 为声子布居数。考虑到在M-P模型中，声子频率只与波矢有关而与偏振方向无关，因此上式中声子布居数 ${\bar{n}}_{k\sigma }$ 也只与波矢有关而与偏振方向无关，可表示为 ${\bar{n}}_k$，加之考虑到前面得到的声子群速和线宽也与偏振方向无关的结论，所以，上式中各求和项只与波矢有关而与偏振方向无关，可以取消对不同的偏振方向x、y和z的求和，而将对波矢 $k$ 的求和结果直接乘以3倍，最终结果也可表示为

${\kappa }_x=\frac{3k_B{\beta }^2}{V}\underset{k}{{\displaystyle \sum }}{\hslash }^2{\omega }_k^2{L}_k^x{v}_k^x{\bar{n}}_k\left({\bar{n}}_k+1\right)$ (23)

由上式可知，声子对热传导系数的贡献取决于群速度、声子自由程、声子数和单个声子能量，其中，声子自由程与声子谱线宽度和声子寿命密切相关。

6. 结论

本文运用Hardy能量通量公式，并根据单原子层薄膜的晶格振动频率、群速度、原子位移公式和动量公式推导了单原子层薄膜的晶格振动能量通量公式，在此基础上运用Green函数理论推导了声子线宽公式，该公式相对于Maradudin和Turney等人的公式更为合理，最后用Green-Kubo公式推导了单原子层薄膜的热传导系数公式，为后续的数值计算打下了基础。由于后续的数值计算完全建立在我们已经得到的解析公式基础上，因此对计算的条件要求比较低，计算速度比较快，可以对原子或分子数较大的体系进行模拟计算，同时也适用于存在量子效应的场合，得到的结果误差较小。更重要的是借助于通过晶格动力学推导得到的有关热传导系数的公式，我们可以对各个声子对热传导系数的贡献进行分析，从而更加深入地探究单原子层薄膜中热传导性质的物理机制。