楽しく学ぶ…統計力学
調和振動子の応用(2) デバイ比熱
同じ比熱という事でアインシュタイン比熱の次に解説される。しかし導出過程は黒体輻射と酷似している。
デバイモデル
アインシュタインモデルは角振動数\(~\omega~\)が全て同じ\(~3N~\)個の調和振動子として扱った。デバイモデルでは角振動数\(~\omega~\)が全て異なる\(~3N~\)個の調和振動子として取り扱う。
\(\omega~\)が同一の系全体の分配関数は\(~Z(T)=\{Z_1(T)\}^N~\)で, その後の計算も容易であった。しかしデバイモデルでは各調和振動子の\(~\omega~\)が異なるため, 一つの粒子の分配関数の\(~N~\)乗という訳にはゆかない。相互作用は無いので, 分配関数の乗積\(~Z=\prod_{i=1}^{N}Z_i~\)として求める事も不可能ではないが, それには角振動数に関する状態密度\(~D(\omega)~\)を初めとして多くの仮定が必要となる。
ここでは, 今後頻出する
波数空間での状態数を求める手法を用いることにする。離散量(状態数)を連続量で表すという意味では, 正統派のカノニカル(正準形式)アンサンブルかもしれない。
デバイモデルの説明はどうしても長くなるが, それは弾性波の解説が理由である。周知の読者は(3)式
状態密度へ跳んで良い。
格子振動の量子化 弾性波の量子化
格子定数\(~a~\), 一辺に\(~L~\)個の原子が並ぶ立方体とする。つまり体積\(~V=(La)^3\), 原子数\(~N=L^3~\)である。
この固体を等方的な弾性体と見做し, 各点の変位を\(~f(\bm{r})~\) とすると, 通常の波動方程式
\[\frac{\partial^2f}{\partial t^2}=v^2\nabla^2f\tag{1} \]
が成り立つ。\(v~\)は音速である。一般解は難解であるが, 今の場合は仮定して解ければそれでよいので, 変数分離型の解を仮定する。
\[f(r,t)=X(x)Y(y)Z(z)T(t) \]
(1)式に当てはめて
\[XYZ\frac{d^2T}{dt^2}=v^2\left(YZT\frac{d^2X}{dx^2}+ZXT\frac{d^2Y}{dy^2}+XYT\frac{d^2Z}{dz^2}\right) \]
両辺を\(XYZT~\)で割ると
\[\frac{1}{T}\frac{d^2T}{dt^2}=v^2\left(\frac{1}{X}\frac{d^2X}{dx^2}+\frac{1}{Y}\frac{d^2Y}{dy^2}+\frac{1}{Z}\frac{d^2Z}{dz^2}\right)\tag{2} \]
各項が正だと発散してしまうので, 負と仮定し,
\[\frac{d^2T}{dt^2}=-\omega^2T,\;\frac{d^2X}{dx^2}=-k_x^2X,\;\frac{d^2Y}{dy^2}=-k_y^2Y,\;\frac{d^2Z}{dz^2}=-k_z^2Z \]
これを(2)式に代入して
\[\omega^2=v^2(k_x^2+k_y^2+k_z^2) \]
であれば解となる。分散関係は
\[\omega=v|\bm{k}|,\bm{k}=(k_x,k_y,k_z)\]
となっている(直線関係)ので, 低温付近を記述する式としては適切であろう。各座標について
\[\frac{\partial^2X}{\partial x^2}=-k_x^2X\cdots \]
であるから, 単振動の方程式を解けば良い。この一般解は量子力学ではないので実数解で
\[X(x)\propto sin(k_xx+\alpha_x) \]
としてよい。振幅は任意で良い。\(k_x~\)と\(\alpha_x~\)は境界条件で決まる。
固定端(自由端としても同じ)とすると\(~x=0,x=La~\)で\(~X(x)=0~\)であるから位相は\(~\alpha_x=0~\)である。さらに\(~X(x)\propto sink_x x~\)より
\[k_x=\frac{\pi}{La}n_x\ (n_x=1,2,\cdots) \]
であり, \(k_x~\)は等間隔\(~\displaystyle \frac{\pi}{La}~\)である。\(Y,Z~\)も同様であり, 系が大きくなる, つまり\(~L~\)が大きくなると取り得る波数は連続値に近づく。
一つ一つの弾性波は三つの波数の組\(~(k_x,k_y,k_z)~\)で指定される。この一つ一つのことをモードと呼ぶ。「～の波数を持つモード」という言い方をする。モードを一つ決めると分散関係より\(~\omega=v|\bm{k}|~\)が決まる。各モードの運動は
\[\frac{d^2T}{dx^2}=-\omega^2T \]
より求められ,\(~T=Acos(\omega t+\delta)~\)である。\(A,\delta~\)は任意である。
波動関数の各解は互いに干渉しないので, 各モードは夫々独立した調和振動子である。つまり粒子の相互作用の結果として波動方程式が成立するのであるが, 波動方程式の解である
各モードは独立した調和振動子である。
状態密度 波数\(~\bm{k}~\)空間で状態数を数える。
角振動数の分布
\[\frac{\omega^2}{v^2}=k_x^2+k_y^2+k_z^2\tag{3} \]
は波数空間で半径\(~\omega/v~\)の球面を表している。その球面上には同じ角振動数\(~\omega~\)を持つモードが存在する。
各格子点上に一つのモードが存在し, 球面の内側の格子点の数が, 角振動数\(~\omega~\)以下のモードの数である。
\(\omega~\)以下のモードの数を\(~D_0(\omega)~\)とする。\(k_x,k_y,k_z\gt 0~\), \((La)^3=V~\)を考慮すると
\[\begin{align}
D_0(\omega)&=\frac{1}{8}\x\frac{4}{3}\pi\left(\frac{\omega}{v}\right)^3/\left(\frac{\pi}{La}\right)^3 \\
&=\frac{V\omega^3}{6\pi^2v^3}
\end{align} \]
各波数ごとに横波\(~(transvers\ mode)\x 2~\), 縦波\(~(longitudinal\ mode)\x 1~\)あり, それぞれの速度は異なるが全て同じと近似し, ３倍すると,
\[D_0(\omega)=\frac{V\omega^3}{2\pi^2v^3} \]
となる。[\(\omega,\omega+\Delta\omega\)] にあるモードの数は\(~\Delta\omega~\)が小さい時,
\[D_0(\omega+\Delta\omega)-D_0(\omega)=\frac{dD_0(\omega)}{d\omega}\Delta\omega \]
と書ける。\(D(\omega)\equiv D_0(\omega)/d\omega~\)と定義すると, 角振動数に関する状態密度\(~D(\omega)~\)
\[D(\omega)=\frac{3V\omega^2}{2\pi^2v^3}\propto \omega^2\tag{4} \]
を得る。正準集団の方法(1)で求めたのは
エネルギー状態密度\(~\Omega(\varepsilon)~\)である。分配関数は\(~\Omega(\varepsilon)~\)から求める。\(D(\omega)~\)からは分配関数は求められない。
さて(4)式のままだと無限に大きい波数即ち, 無限に大きい角振動数が可能となってしまう。元々の自由度は粒子の連成振動だったので有限で\(~3N~\)であったが, これを連続体で近似した。あたかも無限大までとれるように錯覚してしまうが, 実際には上限\(~\omega_D~\)がある。
最小の\(~\omega_{min}~\)は最大の波長(結晶の大きさ程度)\(~\lambda_{max}\sim La~\)の時で,\(~\omega_{min}~\)から数えて\(~3N~\)個のモードを使うとする。全てを加えた時最後が\(~\omega_D~\)になれば良いので,
\[D_0(\omega_D)=\frac{V\omega_d^3}{2\pi^2v^3}=3N \]
\[\omega_D^3=6\pi^2v^3\left(\frac{N}{V}\right) \]
である。\(\omega_D~\)をデバイ周波数と呼ぶ。\(N/V~\)は単位体積当たりの原子数, 数密度である。
内部エネルギー 角振動数状態密度\(~D(\omega)~\)からは分配関数は算出出来ない。
格子振動による内部エネルギーは各モードを調和振動子として全て足せば良い。調和振動子(1)量子力学では内部エネルギーは
\[U(T,N)=-\dd{}{\beta}logZ=N\hbar\omega\left(\frac{1}{2}+\frac{1}{e^{\beta\hbar\omega}-1}\right) \]
零点振動を無視すれば, 各モードの内部エネルギーは
\[u_{\omega}(T)=\frac{\hbar\omega}{e^{\beta\hbar\omega}-1}\]
となる。従ってデバイモデルの格子振動の全エネルギーは上式及び(4)式を使って
\[\begin{align}
U(T,N)&=\sum_{mode}u_{\omega}(T) \\
&=\int_0^{\omega_D}D(\omega)u_{\omega}(T)d\omega \\
&=\frac{3\hbar V}{2\pi^2v^3}\int_0^{\omega_D}\frac{\hbar\omega^3}{e^{\beta\hbar\omega}-1}d\omega \\
&=\frac{9N\hbar}{\omega_D^3}\int_0^{\omega_D}\frac{\hbar\omega^3}{e^{\beta\hbar\omega}-1}d\omega
\end{align} \]
２行目は, 角周波数\(~\omega~\)の近辺のモード数\(D(\omega)d\omega~\)\(~\x~\)各モードのエネルギー\(~u_{\omega}(T)~\)を表している。
以前と同じように\(~x=\beta\hbar\omega~\)と変数変換すると,
\[U(T,N)=\frac{9N\hbar\omega_D}{(\beta\hbar\omega_D)^4}\int_0^{\beta\hbar\omega_D}\frac{x^3}{e^x-1}dx\tag{5} \]
となる。全エネルギーであるから,\(~N~\)および\(~\hbar\omega_D~\)に比例するはずであり, まとめて前に出した。
(5)式の積分は解析的には解けないので, ここでは高温領域と低温領域を見てみる。
高温領域
エネルギー尺度として\(~\hbar\omega_D~\)を用いよう。高温は\(~k_BT\gg\hbar\omega_D~\), すなわち\(~\beta\hbar\omega_D\ll 1~\)の時である。
\(e~\)をテイラー展開\(~\displaystyle e^x=1+x+\frac{1}{2}x^2+\cdots\)して(5)式に適用すると,
\[\begin{align}
U(T,N)&=\frac{9N\hbar\omega_D}{(\beta\hbar\omega_D)^4}\int_0^{\beta\hbar\omega_D}\frac{x^3}{e^x-1}dx \\
&=\frac{9N\hbar\omega_D}{(\beta\hbar\omega_D)^4}\int_0^{\beta\hbar\omega_D}\frac{x^3}{x}dx \\
&=\frac{9N\hbar\omega_D}{(\beta\hbar\omega_D)^4}\frac{1}{3}\left[x^3\right]_0^{\beta\hbar\omega_D}\\
&=3Nk_BT
\end{align} \]
である。これはエネルギー等分配を表している。比熱は\(~T~\)で微分して, \(C_V=3Nk_B~\)(デュロン・プティの式)である。
物性論では\(~k_B\Theta_D=\hbar\omega_D~\)でデバイ温度\(~\Theta_D~\)を定義し, これを温度の尺度として用いる。物質同士の比較に都合が良いからである。統計力学では\(~\hbar\omega_D~\)のままの方が良いと思う。
低温領域
\(~\beta\hbar\omega_D\gg 1~\)であるから, 積分の上限を\(~\infty~\)としてしまえば(5)式は
\[\begin{align}
U(T,N)&=\frac{9N\hbar\omega_D}{(\beta\hbar\omega_D)^4}\int_0^{\infty}\frac{x^3}{e^x-1}dx \\
&=\frac{9N\hbar\omega_D}{(\beta\hbar\omega_D)^4}\frac{\pi^4}{15} \\
&=\frac{3\pi^4}{5}N\hbar\omega_D\left(\frac{k_BT}{\hbar\omega_D}\right)^4\\
\end{align} \]
となる。\(C_V=\partial U/\partial T~\)より,
\[C_V=\frac{12\pi^4}{5}Nk_B\left(\frac{k_BT}{\hbar\omega_D}\right)^3 \]
を得る。低温領域では熱容量(比熱)は\(~T^3~\)に比例し, 実験結果との一致はアインシュタインモデルよりも良い。
波数空間 波数ベクトル\(~\bm{k}~\)の導入。周知の読者は読む必要はありません。
統計力学では\(~\bm{k}(k_x,k_y,k_z)~\)は使えれば, その意味を必ずしも知る必要は無いが, 波数ベクトルのおさらいをしておこう。
位相の等しい面が平面になっている波を平面波と呼ぶ。
速度\(~v~\)で進む進行波は
\[u(x,t)=Asin(kx-\omega t)\quad or \quad u(x,t)=f(x-\upsilon t) \]
で表される。\(~u(x,t)=Asin(kx-\omega t)~\)は, \(~x=\)一定の所で位相が等しい（同時刻\(~t~\)で同じ値を取る）ので, 平面波であり, \(x~\)方向に進む波を表す。
今, 進行方向が\(~x~\)ではなく, 任意の方向を向いているときはどうすれば良いのであろうか？
それは簡単で, 座標系をどのように設定しようと物理現象は変わらないので, 座標系を回転して\(~x~\)軸を波の進行方向に取り直せば良い。
つまり\(~x~\)軸も\(~y~\)軸も\(~z~\)軸も, 都合の良い向きに設定すれば良いのである。
原点を含む波面\(~E~\)上の点 \(\rm O(0,0,0)~\)上の媒質の変位を
\[u(0,0,0,t)=Asin(\omega\,t +\phi)\]
とする。この波面が速度\(~\bm{\upsilon}~\)で動く。時刻\(~\triangle t~\)秒後の波面を\(~E'~\)とし, \(~E'~\)上の任意の点\(~P~\)の変位\(~u(x,y,z,t)~\)を求めてみよう。
点\(~\rm O~\)の, 波面\(~\rm E'~\)上への移り先を\(~\rm O'~\)とする。点\(~P~\)の変位と点\(~\rm O'~\)の変位は等しい（平面波の定義）。
今, 時刻\(~t~\)の点\(~P~\)の変位, すなわち点\(~\rm O'~\)の変位を考える。同時刻に点\(~\rm O~\)も変位しているが, 点\(~\rm O'~\)と同じではない。
点\(~\rm O'~\)の変位は, 時刻\(~t-\triangle t~\)の点 O の変位が遅れてやって来る。
つまり\(~\triangle t~\)秒前の点\(~\rm O~\)の変位\(~u(0,0,0,t-\triangle t)~\)が点\(~\rm O'~\)の時刻\(~t~\)における変位である。
\(~\rm O~\)から\(~\rm O'~\)へ向かう単位ベクトルを\(~\bm{e}~\)とする。\(\rm O~O'~\) の距離\(~s~\)は, \(s=\bm{r}\cdot\bm{e}~\)だから\(~\triangle t~\)は
\[\triangle t=\frac{s}{\upsilon}=\frac{\bm{r}\cdot\bm{e}}{\upsilon}\]
これより
\[\begin{align}
u(x,y,z,t)
&=u(0,0,0,t-\triangle t)\\
&=u\left(0,0,0,t-\frac{\bm{r}\cdot\bm{e}}{\upsilon}\right)\\
&=Asin\left(\omega\left(t-\frac{\bm{r}\cdot\bm{e}}
{\upsilon}\right)+\phi \right)\\
&=Asin\left(\omega t-\left(\frac{\omega}
{\upsilon}\bm{e}\right)\cdot\bm{r}+\phi \right)\tag{6}\\
&=Asin(\omega t-\bm{k}\cdot\bm{r}+\phi)\tag{7}\\
\end{align}\]
(6)式から(7)式へ移る途中で
\[\frac{\omega}{\upsilon}=\frac{2\pi f}{f\lambda}=k \]
を用いた。\(k~\)は波数である。周期\(~2\pi~\)の中に, どれだけ波があるかを表している。
\[\bm{k}=k\cdot\bm{e}\]
を波数ベクトルと呼び, 波の進行方向へ向かう, 大きさ\(~k~\)のベクトルである。
激しく振動する波は波数ベクトルの大きい波, ゆっくり振動する波は波数ベクトルの小さい波である。
