固体物理学⑦ 格子比熱とアインシュタイン模型・デバイ模型の関係
格子比熱からデュロン・プティの法則の導出
固体の格子振動とフォノン
格子振動を原子同士がばねでつながっていると考えていたモデルが一般的になっています．固体物理学の格子振動の記事を参考にすると，1つのフォノンのエネルギーは，\begin{align}
E=\left(n+\dfrac{1}{2}\right)\hbar\omega
\end{align}
となります．これをたくさんのフォノンがある系に拡張すると，
\begin{align}
E=\sum_{k}\sum_{s}\left(n_{k,s}+\dfrac{1}{2}\right)\hbar\omega_{k,s}
\end{align}
となります．この\(k,s\)はそれぞれ波数と振動のモード(振動の種類)を表します．
波数として考える範囲
この波数の取り得る範囲ですが，\(-\dfrac{\pi}{a}\leq k\leq \dfrac{\pi}{a}\)だけになります．これはなぜかというと，波数が\(\frac{2\pi}{a}\)ずれると振動する原子にだけみれば同じ振動とみなせるからです．たとえば，原子が\(x=na\)(\(n\)は自然数)にあるとすると，波数の異なる二つの正弦波\(\sin{kx}\),\(\sin{\left\{\left(k+\frac{2\pi}{a}\right)x\right\}}\)は
\begin{align}
\sin{\left\{\left(k+\dfrac{2\pi}{a}\right)na\right\}}=\sin{(ka+2\pi n)}=\sin{ka}
\end{align}
となり，区別ができないことになります．という理由で波数の範囲は限定されます．
カノニカル分布によって粒子数の平均値を求める．
後で理由を補足しますが，グランドカノニカル分布ではなくカノニカル分布で(\(\mu=0\)として)考えることができます．フォノン数の期待値は\begin{align}
\braket{n_{k,s}}&=\dfrac{\sum_{n=0}^\infty n
\exp{\left\{-\frac{\left(n+\frac{1}{2}\right)\hbar\omega_{k,s}}{k_BT}\right\}}}{\sum_{n=0}^\infty
\exp{\left\{-\frac{\left(n+\frac{1}{2}\right)\hbar\omega_{k,s}}{k_BT}\right\}}}\\
&=\dfrac{\sum_{n=0}^\infty n
\exp{\left(-\frac{n\hbar\omega_{k,s}}{k_BT}\right)}}{\sum_{n=0}^\infty
\exp{\left(-\frac{n\hbar\omega_{k,s}}{k_BT}\right)}}
\end{align}
この式が見えにくいので\(x=\dfrac{\hbar\omega}{k_BT}\)とおきましょう．
\begin{align}
\braket{n_{k,s}}&=\dfrac{\sum_{n=0}^\infty ne^{-nx}}{\sum_{n=0}^\infty e^{-nx}}\\
&=-\dfrac{\sum_{n=0}^\infty \frac{d}{dx}e^{-nx}}{\sum_{n=0}^\infty e^{-nx}}\\
&=-\dfrac{\frac{d}{dx}\sum_{n=0}^\infty e^{-nx}}{\sum_{n=0}^\infty e^{-nx}}\\
&=-\dfrac{d}{dx}\ln{\left(\sum_{n=0}^\infty e^{-nx}\right)}
\end{align}
さて，この対数の引数部分は，
\begin{align}
\sum_{n=0}^\infty =\dfrac{1}{1-e^{x}}
\end{align}
となるので，
\begin{align}
\braket{n_{k,s}}=\dfrac{1}{\exp{\left(\frac{\hbar\omega_{k,s}}{k_BT}\right)}-1}
\end{align}
となります．
ボース・アインシュタイン分布関数を用いる方法
いま，導いた式は\(\mu=0\)としたボース・アインシュタイン分布関数でした．つまり，出発点をボース・アインシュタイン分布にしてもいいかもしれませんね．ボース・アインシュタイン分布関数よりエネルギーが\(E=\hbar\omega\)の準位の粒子数の期待値は，以下のように表されます．ちなみにこの場合には\(\mu=0\)とできます．\begin{align}
\braket{n}=\dfrac{1}{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}-1}
\end{align}
化学ポテンシャルはどこへ?
今回はフォノンを扱っていますが，これは粒子数が保存されません．グランドカノニカル分布は粒子数が変動するといっても系と外系のやりとりによって系の粒子数が変動するわけです．つまり，結局，系と外系を合わせた粒子数は一定です．ただし，今回の場合にはフォノンの数が一定ではない,つまり粒子浴を考えなくても良いわけです．
もし，ボース・アインシュタイン関数の導出を覚えているなら，そのラグランジュの未定乗数法を用いた解析の中に粒子数一定の束縛条件を課していたことを思い出しましょう．この部分から出てくるのが，ボース・アインシュタイン分布に現れる\(\mu\)でしたが，今回は粒子数に関する束縛条件はないのでこの項は出てきません．
エネルギーの微分から格子比熱を求める
エネルギーの平均値はフォノンの平均値\(\braket{n_{k,s}}\)を用いて，\begin{align}
\braket{E}&=\sum_{k}\sum_{s}\left(\braket{n_{k,s}}+\dfrac{1}{2}\right)\hbar\omega_{k,s}\\
&=\sum_{k}\sum_{s}\left(\dfrac{1}{\exp{\left(\frac{\hbar\omega_{k,s}}{k_BT}\right)}-1}+\dfrac{1}{2}\right)\hbar\omega_{k,s}
\end{align}
と表され，その格子比熱は，この微分で表されます．第二項は微分で消えるので，
\begin{align}
C&=\dfrac{\partial\braket{E}}{\partial T}\\
&=\dfrac{\partial}{\partial T}\left\{\sum_{k}\sum_{s}\dfrac{\hbar\omega_{k,s}}{\exp{\left(\frac{\hbar\omega_{k,s}}{k_BT}\right)}-1}\right\}
\end{align}
となります．このままだとよくわかりませんね．他の条件にも目を向けてみましょう．
周期的境界条件を課す
\(x\)軸方向に長さ\(L\)の一次元結晶を考えましょう．格子定数を\(a\)として原子数を\(N_x\)と表すと，\begin{align}
L=N_xa
\end{align}
となります．周期的境界条件，つまり，波がちょうど結晶で位相が\(2\pi\)の整数倍になればいいことになります．つまり，波数は，整数\(l\)を用いて
\begin{align}
k=\dfrac{2\pi l}{L}=\dfrac{2\pi l}{N_xa}
\end{align}
となります．ここで，辺々の微小変化を取ると，
\begin{align}
dk_x=\dfrac{2\pi}{N_xa}dl
\end{align}
つまり，波数変化に対する取り得る波数の数の変化は，
\begin{align}
\dfrac{dl}{dk_x}=\dfrac{N_xa}{2\pi}
\end{align}
これは一次元での話なので，3次元に拡張します．他の方向にも同様の格子定数で，\(y\)軸方向に原子数\(N_y\),\(z\)軸方向に原子数\(N_z\)とすれば，波数密度は
\begin{align}
\dfrac{N_xN_yN_za^3}{(2\pi)^3}
\end{align}
となります.
積分と和の記号の変換
\(k\to\infty\)で以下のように和と積分が入れ替え可能です．\begin{align}
\sum_{k}f(k)=\int dk_xdk_ydk_z \dfrac{N_xN_yN_za^3}{(2\pi)^3}f(k)
\end{align}
この式で大事なことは，
- \(N_x,N_y,N_z\)が十分大きいからとりうる\(k\)の間隔は非常に小さく連続的にみなせる
- 波数\(k\)が波数空間での密度の逆数分変化すれば積分が1変化するようにする．
\begin{align}
C=\dfrac{\partial}{\partial T}\left\{\dfrac{N_xN_yN_za^3}{(2\pi)^3}\sum_s\int dk_xdk_ydk_z\dfrac{\hbar\omega_{k,s}}{\exp{\left(\frac{\hbar\omega_{k,s}}{k_BT}\right)}-1}\right\} \label{eq:sp7.23}
\end{align}
角振動数について一定，線形の仮定をする．
さて，いま，\(\omega\)は\(k\)に依存する可能性があるのですが，格子振動の記事では\(\omega\)は\(k\)の三角関数で表されたのでした．これがネイピア数の指数部分にくることになります．これは積分できませんね...仕方ないので仮定をして乗り越えることにします．具体的には角振動数がアインシュタイン模型では一定に，デバイ模型では波数に比例すると仮定します．
アインシュタイン模型での和の処理
アインシュタインモデルでの仮定では，\(\) ここで，\(dk_xdk_ydk_z\)の積分の処理を考えましょう．いま，被積分関数にある\(\omega\)が\(k\)に依存しない定数として考えましょう．つまり，積分の外に被積分関数を出すことができます．また，振動は進行方向の振動(縦波)1成分と進行方向に垂直な振動(横波)2成分があり，この3種類のみです．つまり，\(\sum_s\)は3倍すればよく，\begin{align}
C=\dfrac{\partial}{\partial T}\left\{\dfrac{3N_xN_yN_za^3}{(2\pi)^3}\dfrac{\hbar\omega}{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}-1}\int dk_xdk_ydk_z\right\} \label{eq:sp7.24}
\end{align}
さて，一番後ろの積分項はたとえば\(x\)成分に関しては\(-\dfrac{\pi}{a}\leq k\leq\dfrac{\pi}{a}\)となります．\(y,z\)成分も同様で，積分は
\begin{align}
\int dk_zdk_ydk_z=\dfrac{(2\pi)^3}{a^3}
\end{align}
となります．原子数の総数は\(N=N_xN_yN_z\)なので，\eqref{eq:sp7.24}は以下のように変形できます．
\begin{align}
C=\dfrac{\partial}{\partial T}\left\{\dfrac{3N\hbar\omega}{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}-1}\right\}
\end{align}
この偏微分を計算すると，
\begin{align}
C=\dfrac{3N\hbar^2\omega^2}{k_BT^2}\dfrac{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}}{\left\{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}-1\right\}^2}
\end{align}
これをあえて以下のように変形してみましょう．アボガドロ定数と気体定数に関して\(N_Ak_B=R\)，\(n=\frac{N}{N_A}\)はモル数を表すので，
\begin{align}
C&=3Nk_B\left(\dfrac{\hbar\omega}{k_BT}\right)^2\dfrac{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}}{\left\{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}-1\right\}^2}\\
&=3nR\left(\dfrac{\hbar\omega}{k_BT}\right)^2\dfrac{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}}{\left\{\exp{\left(\frac{\hbar\omega}{k_BT}\right)}-1\right\}^2}
\end{align}
\(x=\dfrac{\hbar\omega}{k_BT}\)とおいて，高温極限では\(T\to\infty\)より\(x\to0\)の極限を求めれば，\(3R\)になることがわかります．高温極限で単位物質量当たりの格子比熱が\(3R\)に近づくことをデュロン・プティの法則といいます．
デバイ模型での解析
アインシュタイン模型を適用する前の式をもう一度．以下の式です．ただし原子数だけおきかえています．\begin{align}
C=\dfrac{\partial}{\partial T}\left\{\dfrac{Na^3}{(2\pi)^3}\sum_s\int dk_xdk_ydk_z\dfrac{\hbar\omega_{k,s}}{\exp{\left(\frac{\hbar\omega_{k,s}}{k_BT}\right)}-1}\right\}
\end{align}
今回は\(\omega=v|\boldsymbol{k}|\)という形，波の伝搬速度が1成分の縦波と2成分の横波で変わらないということを仮定します．つまり，\(v\)が波の速度ということですね．このとき，\(s\)についての和は3倍することで解決します．\(k=|\boldsymbol{k}|\)として，
\begin{align}
C=\dfrac{\partial}{\partial T}\left\{\dfrac{3Na^3}{(2\pi)^3}\int dk_xdk_ydk_z\dfrac{\hbar vk}{\exp{\left(\frac{\hbar vk}{k_BT}\right)}-1}\right\}
\end{align}
さて，積分を計算しましょう．ところで，被積分関数の\(k\)というのは波数ベクトルの大きさです．つまり，この体積分は方向に依存しないので\(k\)についての1変数の積分に置き換わります．\(dk_xdk_ydk_z=4\pi k^2dk\)であり，整理すると，
\begin{align}
C=\dfrac{\partial}{\partial T}\left\{\dfrac{3N\hbar va^3}{2\pi^2}\int \dfrac{k^3}{\exp{\left(\frac{\hbar vk}{k_BT}\right)}-1}dk\right\}
\end{align}
ここで，温度による偏微分をすると，
\begin{align}
C&=\dfrac{3N\hbar va^3}{2\pi^2}\dfrac{\hbar v}{k_BT^2}\int \dfrac{k^4\exp{\left(\frac{\hbar vk}{k_BT}\right)}}{\left\{\exp{\left(\frac{\hbar vk}{k_BT}\right)}-1\right\}^2}dk\\
&=\dfrac{3N\hbar^2v^2a^3}{2\pi^2k_BT^2}\int \dfrac{k^4\exp{\left(\frac{\hbar vk}{k_BT}\right)}}{\left\{\exp{\left(\frac{\hbar vk}{k_BT}\right)}-1\right\}^2}dk
\end{align}
積分を実行するために\(\dfrac{\hbar vk}{k_BT}=x\)とおいて，置換積分を行います．
\(dk=\dfrac{k_BT}{\hbar v}dx\)なので，
\begin{align}
C&=\dfrac{3N\hbar^2v^2a^3}{2\pi^2k_BT^2}\int \dfrac{\left(\frac{k_BT}{\hbar v}x\right)^4e^x}{\left(e^x-1\right)^2}\dfrac{k_BT}{\hbar v}dx \nonumber\\
&=\dfrac{3Nk_B^4a^3}{2\pi^2\hbar^3v^3}T^3\int \dfrac{x^4e^x}{\left(e^x-1\right)^2}dx \nonumber \\
&=3Nk_B\dfrac{a^3}{2\pi^2}\left(\dfrac{k_BT}{\hbar v}\right)^3\int \dfrac{x^4e^x}{(e^x-1)^2}dx \label{eq:sp7.35}
\end{align}
積分範囲を設定しなおす
ところで，積分範囲について考えてみましょう．元々の積分範囲は波数空間で各成分に対して，\(-\frac{\pi}{a}\leq k_i\leq\frac{\pi}{a},i=x,y,z\)となります．つまり，この波数空間上の体積は\begin{align}
\dfrac{(2\pi)^3}{a^3}
\end{align}
となります．ただし，さきほど，被積分関数が波数空間で球対称となるため，体積分を\(4\pi k^2dk\)としたのでした．というわけで積分範囲も球対称でないと積分できません．よって波数空間での体積が同じまま，形を球にしましょう．つまり，半径\(k_D\)の球を考え，その体積がもとの領域とおなじになるように，つまり，
\begin{align}
\dfrac{4}{3}\pi k_D^3=\dfrac{(2\pi)^3}{a^3}\ \ \therefore k_D=\left(\frac{6\pi^2}{a^3}\right)^\frac{1}{3}
\end{align}
この\(k_D\)をデバイ波数といいます．このデバイ波数を用いて，\eqref{eq:sp7.35}から\(a\)を消去しましょう．
\begin{align}
\dfrac{a^3}{2\pi^2}=\dfrac{3}{k^3_D}
\end{align}
となるので，
\begin{align}
C=9Nk_B\left(\dfrac{k_BT}{\hbar vk_D}\right)^3\int \dfrac{x^4e^x}{(e^x-1)^2}dx
\end{align}
ここで，デバイ波数に対応する角振動数\(\omega_D=vk_D\)をデバイ角振動数といいます．さて,\(k\leq k_D\)より，
\begin{align}
\dfrac{1}{xT}=\dfrac{k_B}{\hbar vk}\leq \dfrac{\hbar vk_D}{k_B}=\dfrac{\hbar\omega_D}{k_B}=\Theta_D
\end{align}
として，この\(\Theta_D\)をデバイ温度といいます．このデバイ温度を用いて，積分範囲は
\begin{align}
0\leq x \leq \dfrac{\Theta_D}{T}
\end{align}
となる．よって，
\begin{align}
C=9Nk_B\left(\frac{T}{\Theta_D}\right)^3\int_0^\frac{\Theta_D}{T} \frac{x^4e^x}{(e^x-1)^2}dx
\end{align}
デバイモデルによる低温領域の扱い
低温領域，つまり，\(T\ll \Theta_D\)の場合では，\begin{align}
\int_0^\frac{\Theta_D}{T}\frac{x^4e^x}{(e^x-1)^2}dx\to \int_0^\infty \frac{x^4e^x}{(e^x-1)^2}dx=\dfrac{4\pi^4}{15}
\end{align}
より，
\begin{align}
C=\dfrac{12\pi^4}{5}Nk_B\left(\dfrac{T}{\Theta_D}\right)^3
\end{align}
となります．つまり，低温では\(T^3\)に比例する比熱をもちます．
波数の話を本来は逆格子ベクトルで考えるべきという話
格子を今回は波数空間で...という解析をしましたが，周期的境界条件を課す際にはできれば逆格子を使う必要があります．(必ずしも格子が直交座標だとは限らないので)波数ベクトルと基本並進ベクトルの間にの関係を考えます．たとえば，基本並進ベクトル\(\boldsymbol{a_1}\)と波数ベクトルとの関係は，この方向に基本単位胞が\(N_1\)個並んで周期的になると考えて，整数\(m_1\)を用いて，
\begin{align}
N_1\boldsymbol{a_1}\cdot\boldsymbol{k}=2\pi m_1
\end{align}
となります．これは2,3成分についても同様で,
\begin{align}
N_i\boldsymbol{a_i}\cdot\boldsymbol{k}&=2\pi m_i\\
\therefore \boldsymbol{a_i}\cdot\left(\dfrac{N_i\boldsymbol{k}}{m_i}\right)&=2\pi
\end{align}
ところで，逆格子基本ベクトルと基本並進ベクトルの間には，
\begin{align}
\boldsymbol{a_i}\cdot\boldsymbol{b_j}=2\pi \delta_{ij}
\end{align}
という関係がありました．つまり，この2つの式を考えると波数ベクトルが，
\begin{align}
\boldsymbol{k}=\dfrac{m_1}{N_1}\boldsymbol{b_1}+\dfrac{m_2}{N_2}\boldsymbol{b_2}+\dfrac{m_3}{N_3}\boldsymbol{b_3}
\end{align}
だとよさそうですね．ところで，1次元の逆格子ベクトルの大きさは結晶の長さ\(L=N_1a\)に対して\(\dfrac{2\pi n}{L}=\dfrac{2\pi n}{Na}\)と表されたのでした．というわけで，たとえば，\(\boldsymbol{b_1}\)方向の波数を\(-\dfrac{\pi}{a}\leq k_1\leq \dfrac{\pi}{a}\)にするためには，
\begin{align}
-\dfrac{1}{2}N_1\leq m_1\leq \dfrac{1}{2}N_1
\end{align}
となります．他の成分も同様です．さて，というわけで波数の和を積分に入れ替えるには，
\begin{align}
\sum_k f(\boldsymbol{k})=\int_{-\frac{N_1}{2}}^{\frac{N_1}{2}} dm_1\int_{-\frac{N_2}{2}}^{\frac{N_2}{2}} dm_2\int_{-\frac{N_3}{2}}^{\frac{N_3}{2}} dm_3 f\left(\dfrac{m_1}{N_1}\boldsymbol{b_1}+\dfrac{m_2}{N_2}\boldsymbol{b_2}+\dfrac{m_3}{N_3}\boldsymbol{b_3}\right)
\end{align}
実質これは逆格子の単位胞で積分するということです．